Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

A mechanistic ecohydrological model to investigate complex interactions in cold and warm water‐controlled environments: 2. Spatiotemporal analyses

A mechanistic ecohydrological model to investigate complex interactions in cold and warm... 1. Introduction The studies of dynamics of hydrological and vegetation processes have been historically more focused on the understanding of mechanisms and controls at the point/plot spatial scale rather than at the scale of topographically complex landscapes. At the plot scale, numerous theoretical [e.g., , 1999 ; , 2001 ; , 2005 ], field experimental work [e.g., , 2004 ; , 2004 , 2010 ; , 2010 ], and novel modeling approaches [ , 2008 , 2009 ; , 2010 ] have been presented. Recently, there was an emergence of interest in extending ecohydrological analysis to larger spatial scales in complex landscapes [ , 2008b ; , 2008 , 2010 ], where the two‐way coupling between vegetation dynamics and water cycle is far less studied. For instance, this is the case for mountainous environments, where topographic, vegetation, and climatic gradients play a crucial role in controlling hydrological, ecological, and biogeochemical behaviors [ , 2008 ; , 2008 ; , 2009 ; , 2009 ; , 2011 ]. The effects of a spatially varying structure of partition and magnitude of the incoming radiative energy, the role of water redistribution, and other terrain‐related hydrological processes are much less understood at the watershed scale than at the local scale [ , 2010 ]. This can be explained by the spatially and temporally heterogeneous conditions of watersheds and the difficulty of observing or accounting for such differences in a parsimonious manner in models. A distributed representation of ecohydrological dynamics and a subsequent aggregation to the watershed scale are thus warranted. Both analytical [ , 2004 , 2005 ] and numerical methods [ , 2001 ; , 2008b ; , 2009 ; , 2010 ] have been proposed to extend ecohydrological analysis to larger spatial scales. In this study, a numerical approach is pursued using Tethys‐Chloris, a novel ecohydrological model presented in the companion paper [ , 2012 ]. Tethys‐Chloris (T&C) is an ecohydrological model that reproduces all essential components of hydrological cycle, resolving the mass and energy budgets at the hourly scale. It includes a description of energy and mass exchanges in the atmospheric surface layer accounting for up to two layers of vegetation, a module of saturated and unsaturated soil water dynamics, and a snowpack evolution module. Vegetation dynamics are also reproduced and include all essential plant life‐cycle processes, i.e., photosynthesis, phenology, carbon allocation, and tissue turnover [ , 2012 ]. The first confirmation of the model was presented in the companion paper for two plot‐scale case studies [ , 2012 ]. In this work, the analysis is extended to the watershed scale. The two cases are the Lucky Hills experimental watershed in Arizona (U.S.A.) and the Tollgate/Reynolds Creek Mountain East watersheds in Idaho (U.S.A.). Both represent semiarid environments but while the former is characteristic of a partially vegetated desert shrub system, the latter is vegetated with sagebrush, firs, and aspens and characterized by a seasonal occurrence of snow. Semiarid ecosystems represent a significant fraction of global terrestrial surface area [ , 2005 ] and they have an important role for the Earth's climate [ , 2008 ]. The model application for semiarid environments is also consistent with the assumption of water‐limited vegetation, since nutrient dynamics are currently neglected in Tethys‐Chloris. The model is applied at the distributed spatial scale with the overall objective to assess the capability of reproducing physical and biophysical metrics of dynamics. These are energy fluxes, soil moisture, snowpack, vegetation production, and seasonality of Leaf Area Index (LAI). A distributed analysis at the spatial scale of a watershed presents numerous issues for a confirmation of any mechanistic ecohydrological model. They arise due to the common lack of spatially distributed observations at the appropriate scales for most of the simulated variables, and because of the high uncertainty in specifying the boundary conditions for a system, e.g., meteorological inputs, bedrock depth distribution, and vegetation spatial variations. In the case of Lucky Hills basin and the Tollgate/Reynolds Creek Mountain East watersheds, we identified two suitable catchment systems that have long‐term, quality‐controlled, diverse (energy flux, snowpack, streamflow, etc.) records. Nonetheless, the spatially distributed comparison for the Lucky Hills basin remains severely limited due to its small size that does not allow a sound comparison with remote sensing products. The present paucity of spatially‐distributed data permits only a limited evaluation using remote sensing observations, with more emphasis on an overall qualitative consistency of spatial patterns rather than their explicit quantitative confirmation. The plausibility of simulated variables and insights on the mechanisms of simulated dynamics are also discussed, as indirect metrics of model confirmation. Despite limited calibration (or virtually none as the total number of test runs was less than a dozen) of Tethys‐Chloris, the partial confirmation of the model yielded satisfactorily results. This further encourages the development and improvement of fully‐coupled ecohydrological models for applications in a spatially distributed fashion. The possibility offered by such models for studies of internal functioning of a river basin will be useful for designing virtual experiments, testing hypotheses, and focusing questions of scientific inquiry. Promising directions of research that current applications of Tethys‐Chloris open for further investigation are described in the analysis of results and the concluding discussion. 2. Issues of Model Confirmation Spatially distributed data are typically available for a limited range of metrics over small areas. Specifically, while streamflow data are available for numerous watersheds, comprehensive data on energy flux, soil moisture, snow, and vegetation productivity, exist only for small‐scale experimental watersheds (10 −3 –10 1 km 2 ) and are at a coarse resolution or entirely unavailable at larger scales, as also demonstrated in the presented examples. In this regard, this work further advocates the urgent need for novel distributed data to better validate complex mechanistic models [ , 2011 ; , 2011 ; , 2012 ]. An assessment of performance skill of ecohydrological models of “new generation” [ , 2008b ; , 2009 ] cannot rely only on aggregated metrics, such as streamflow discharge. In addition to a spatial extent, an evaluation of a mechanistic ecohydrological model requires a diverse suit of observational data. Unfortunately, they are rarely captured by a single experimental field design. For example, the interdisciplinary nature of T&C applications requires meteorological, hydrological, vegetation productivity, energy and carbon fluxes measurements to be co‐located and collected at the same period(s) of time. Furthermore, vegetation physiological and structural attributes, as well as soil texture profiles have to be known. The conventional scarcity of interdisciplinary data makes it difficult, or sometimes even impossible, testing all of the desired behavioral aspects of process‐based models [e.g., , 2008a ]. A significant contribution to the problem of meaningful confirmation of mechanistic ecohydrological models has been offered by the “FLUXNET” monitoring program ( www.fluxnet.ornl.govfluxnetindex.cfm ) [ , 2001 ; , 2007 ]. Further progress in unification of measurement protocols and comprehensiveness of collected data will be addressed by the “NEON” network [ , 2008 ] Another promising source of information is offered by remote sensing data, such retrievals from the sensor Moderate Resolution Imaging Spectroradiometer (MODIS) ( http:modis.gsfc.nasa.gov ) [ , 1998 ; , 1998 ; , 2002 ]. Remote sensing observations can be used to test the capability of reproducing structural properties of vegetation canopy, such as Leaf Area Index (LAI) and vegetation phenology. In general, combining observational data from flux towers with measurements of subsurface states (soil moisture) and vegetation attributes allows one to effectively evaluate individual model components as well as the model overall performance. Although the data contributions provided by FLUXNET and MODIS databases are significant, the demand for multidisciplinary, spatially distributed data for testing and improvement of Earth‐science models exceeds data availability (see also discussion in [2011] ). There is a clear need to quantitatively check whether model structures are suitable for representing distributed ecohydrological dynamics across a range of spatial and temporal scales. The generation of appropriate data would require an effort of the entire community and experimental watersheds with multi‐scale and long‐term monitoring plans. In the presented study, the objective of confirmation is to underline that T&C behaves consistently across different scales, as can be inferred from the available data. The term “confirmation” used here is generally preferred to “verification”, “validation”, or “corroboration”. The first two processes are inherently impossible for numerical models of open natural systems [ , 1994 ]. Models can only be evaluated in relative terms and confirmed by a demonstration of agreement between observations and predictions. This keeps the issue of development of more rigorous methods of assessment open [ , 1994 ]. Despite a careful analysis of the model performance and consistency of simulations carried out in this study, some of the model results will remain critical to confirm further. For instance, this is the case for the coupled modeling of vegetation, radiation transfer, and snowpack dynamics. Detailed field experiments to confirm these dynamics are rarely available [ , 2008 ; , 2008 ; , 2009 ], although the scope and aims of intercomparison projects, such as the recent SNOWMIP2 [ , 2009 ; , 2009 ; , 2010 ], provide a promising framework for testing these process interactions at the plot scale. A number of uncertainties can be also identified in the modeling of vegetation components, such as the carbon allocation scheme, the process of tissue turnover, and drought effects on stomatal closure and plant mortality. Such crucial aspects are only partially testable with data [ , 2006 ; , 2009 ; , 2009 ; , 2009 ] and a substantial progress is required to further their ecological and physiological understanding [e.g., , 2001 ; , 2002 ; , 2003 ; , 2007 ; , 2008 ; , 2010 ]. 3. Case Study Two case studies characterized by different vegetation types and climatic conditions were simulated using T&C. The first case study is the Lucky Hills experimental watershed in Arizona, U.S.A. (Section 3.1). The second case study is represented by two nested basins located in the mountainous area of the Reynolds Creek experimental watershed (Idaho, U.S.A.), specifically, the Tollgate and Reynolds Creek Mountain East basins (Section 3.2.2). 3.1. Lucky Hills Watershed The Lucky Hills monitoring site (110.30 W , 31.44 N ; elevation 1372 [ m a.s.l. ]) is located within the Walnut Gulch Experimental Watershed, near Tombstone in the south‐eastern Arizona, U.S.A. Walnut Gulch is a long‐term experimental watershed managed by the Southwest United States Department of Agriculture‐Agriculture Research Service (USDA‐ARS), where research on hydrology and soil erosion started in 1953 [ , 2008 ]. The studied ecosystem is a plant community of desert shrubs. The climate at the Lucky Hills site can be classified as semiarid, hot, with dry winter. The mean annual temperature from meteorological observations (July 1996 through December 2009) is 17.2 [° C ] and the mean annual precipitation is 333 [ mm ] [ , 2008 ]. The Lucky Hills experimental watershed has an area of 0.037 [ km 2 ]. The elevation range is about 10 [ m ] ( Figure 1a ). This small watershed has shallow to moderate slopes of <0.1 [−], except for the central part, where steep fragments of terrain can have slopes of ≈0.4 [−] ( Figure 1b ). 1 A representation of topographic attributes of the Lucky Hills experimental watershed. (a) Digital Elevation Model. (b) Slope fraction [−] calculated with the steepest descent method [ , 2008 ]. Long‐term measurements of runoff and sediment transport have been collected at Lucky Hills over the period of 1963 through 2008 [ , 2007 ; , 2008 ; , 2008 ]. The mean long‐term runoff is estimated to be about 18.2 [ mm yr −1 ] and the mean long‐term sediment transport is 0.42 [ kg m −2 yr −1 ] [ , 2007 ; , 2008 ]. In order to study carbon dioxide and water fluxes in the Walnut Gulch area, two flux towers have been in operation since 1997 [ , 2008 ]. One of the flux towers is located in the Lucky Hills catchment and is a part of the FLUXNET network. The flux tower footprint includes a shrub plant community, mainly composed of evergreen shrubs, such as creosotebush ( Larrea tridentata ), tarbush ( Flourensia cernua ), and deciduous shrubs, such as whitethorn acacia ( Acacia constricta ) [ , 2008 ; , 2008 ]. Average vegetation height is around 0.3–0.5 [ m ] and vegetation cover fraction is approximately 0.25–0.4 [ , 1994 ; , 2001 ]. The soil at this site is coarse loam [ , 2008 ]. The surface horizon (0–6 [ cm ]) contains 650 [ g kg −1 ] of sand, 290 [ g kg −1 ] of silt, and 60 [ g kg −1 ] of clay, with 290 [ g kg −1 ] of coarse fragments larger than 2 [ mm ], 8 [ g kg −1 ] of organic carbon, and 21 [ g kg −1 ] of inorganic carbon [ , 2008 ]. The saturated conductivity has been observed to consistently diminish exponentially with the soil depth [ , 2000 ]. Tethys‐Chloris was used to simulate ecohydrological dynamics that occurred in the Lucky Hills basin over the period of July 1996 through December 2009, or 13.5 years in total. The grid resolution was of 10×10 m 2 with the total number of basic computational elements equal to 371. A two meter soil column with a free drainage boundary condition at the bottom was assumed for all computational elements composing the watershed. The vertical mesh consisted of 18 soil layers with increasing thickness with depth, starting with 10 [ mm ] at the soil surface, and ending with 400 [ mm ] at 2m depth. Spatially homogenous soil hydraulic properties were assumed; they were derived from the pedotransfer functions of [2006] using a 0.75 fraction of sand and 0.10 fraction of clay. The saturated hydraulic conductivity, K v,s , was parameterized to decay with soil depth, K v,s ( z ) = K v,s (0) e (−0.0011 z ) , where z [ mm ] is the soil depth [ , 2000 ], and the soil anisotropy factor was assumed to be a r = 1 [ , 2012 ]. Spatially uniform crown area fractions, C crown , were assigned to be 0.25 for deciduous shrubs, and 0.10 for evergreen shrubs (see [2012] for the definition of C crown ). Physiological and structural characteristics of whitethorn acacia and creosote bush derived from a literature survey were used in the choice of the parameters and properties of deciduous and evergreen shrubs [ , 1986 ; , 1990 ; , 1994 ; , 2006 ; , 2008 ; , 2009 ]. 3.2. Reynolds Creek, Idaho The Reynolds Creek Experimental Watershed (RCEW) managed by the USDA‐ARS Northwest Research center is located in the Owyhee Mountains 80 km southwest of Boise, ID, U.S.A. [ , 2001 ]. The RCEW area was described in a series of papers by [2001] , [2009] , [2010] , and [2011] . The collected dataset is particularly appealing for hydrologic modeling because of high quality, high frequency, and long‐term measurements of climatic and hydrological variables, including snow water equivalent, snow depth, and soil moisture [ , 2001 ; , 2001 ; , 2001 ; , 2001b ; , 2001 ; , 2002 ; , 2002 ; , 2010 ; , 2011 ]. Energy and mass fluxes have recently started being collected [ , 2010 ]. The RCEW domain encompasses an area of 239 [ km 2 ], with an elevation range of about 1200 [ m ] (1101–2241 [ m a.s.l. ]), where a strong topographic precipitation gradient (236–1123 [ mm ]) exists. A strong precipitation seasonality is characteristic of the area: 70% of precipitation occurs during November through May period [ , 2001 ; , 2002 ]. Snow is the dominant form of precipitation, especially in the headwater part of the catchment. All details on data collection can be found in the above references; only a succinct description is presented in the following. The T&C model was used to simulate two nested sub‐watersheds of RCEW ( Figure 2 ): the Reynolds Creek, Mountain East (RME) and the Tollgate watershed (TOL) [ , 2000 , 2001a ; , 2001 ]. 2 A representation of topographic and vegetation attributes of the upper part of the Reynolds Creek experimental watershed. (a) The Digital Elevation Model of the Tollgate watershed. (b) The vegetation map of the Tollgate watershed. (c) The Digital Elevation Model of the Reynolds Creek, Mountain East watershed. The symbols in Figures 2a and 2c indicate the locations of meteorological ( clim. ), precipitation ( prec. ), and snow ( snow ) measurement stations. 3.2.1. Reynolds Creek Mountain East Watershed The Reynolds Mountain East (RME) is a small headwater catchment (0.389 [ km 2 ]) in the southwestern corner of RCEW that ranges in elevation from 2024 to 2139 [ m a.s.l. ] ( Figure 2c ). Hourly meteorological variables, snow depth, snow water equivalent, soil moisture, and discharge for a 25‐year period (October 1983 through October 2008) were collected at two locations ( Figure 2c ). The stations represent major landscape units within the RME catchment: one is a sheltered site, within a clearing of an aspen/fir grove near the center of the catchment; the other is an exposed site dominated by mixed sagebrush near the western catchment divide [ , 2000 ; , 2001 ; , 2001 ; , 2002 ; , 2002 ; , 2009 ; , 2010 ; , 2011 ]. Hourly observations of air temperature, humidity, wind, solar radiation, and precipitation from the wind‐exposed ridge site, along with the hourly data of air temperature, wind, and precipitation, from the grove site were used in the simulation. Snow density and depth observations were available from snow pillow and manual measurements; soil moisture measurements at different depths and locations were also available [ , 2000 , 2001b ; , 2000 , 2001; , 2011 ]. At the sheltered site, during an average water year (October through September), shielded precipitation is 964 [ mm ]. At the exposed site, it is 550 [ mm ] [ , 2011 ]. Prevailing winds blow from the southwest direction with the mean winter‐time direction of 230°. Avalanching does not play a role in redistributing snow in this catchment [ , 2002 ]. The soil in Reynolds Creek, Mountain East area is sandyloam and loam, with elevated contents of rocks and organic fraction [ , 2001b ]. Vegetation of the RME watershed ( Figure 8a ) is characterized by patches of Douglas fir ( Psuedotsuga menziesii ) (4.6%), aspen ( Populus tremuloides ) (26.1%), and sagebrush plants (69.3%), such as low sagebrush ( Artemisia arbuscula ), vaseyana sagebrush, and Wyoming big sagebrush ( Artemisia tridentata ) [ , 2001b ; , 2002 ; , 2002 ; , 2010 ; , 2011 ]. [2010] provide estimates of the maximum leaf area index and vegetation heights: aspens have an average height of 9.5 [ m ] and LAI = 1.35; LAI of fir trees is about 2.0; sagebrush are approximately 60 [ cm ] in height with LAI of 0.77. 8 (a) The vegetation map of the Reynolds Creek Mountain East watershed. (b) The volumetric soil water content averaged over the simulation period of October 1983 through October 2008. Spatial variability of meteorological forcing is a dominant feature impacting hydrology within the Reynolds Creek Experimental watershed [ , 2001 ; , 2010 ]. It justifies the need to accurately interpolate meteorological inputs over the area of the RME watershed. Air temperature and relative humidity are linearly interpolated between the two available stations, accounting for the elevation effect and using the observed gradients of these variables at each hour. Precipitation and wind speed are interpolated accounting for terrain structure and vegetation, following the approach outlined by [2002] and [2002] . Specifically, an accumulation factor for precipitation, and a wind terrain factor are used to spatially interpolate these variables using four classes of terrain and vegetation shelter computed with the parameters S x (maximum upwind slope) and S b (slope break). Lastly, incident solar radiation in the watershed is interpolated accounting for the local and remote shading effects of topography [ , 2011 ]. The ecohydrological dynamics at RME were simulated for the period of October 1983 through October 2008, for a total duration of 25 years. The imposed grid resolution was 30×30 m 2 that led to 433 computational elements. A uniform one meter soil column and nearly impermeable bedrock (conductivity of bedrock equal to 0.01 [ mm h −1 ]) were assumed for all computational elements of the grid. The chosen soil column depth was assumed to be a representative average for the RME watershed, where typical soil depths are between 15 and 200 [ cm ] [ , 2001b , 2009 ]. A spatially uniform value was assumed because of the lack of distributed information. The value of bedrock conductivity was derived from in situ substrate geology [ , 2001a ; , 2011 ]. The vertical mesh was composed of 10 soil layers of increasing thickness, dz , from the surface ( dz = 10 [ mm ]) to a 1m depth ( dz = 200 [ mm ]). Soil hydraulic properties were derived from the pedotransfer functions of [2006] using the spatially variable fractions of sand and clay derived from the soil map [ , 2000 , 2001a ]. Vertical homogeneity of soil hydraulic properties was assumed because of lack of information. An horizontal anisotropy coefficient, a R = 140, was used for all of the grid elements (see Section 4). The vegetation map was used to assign a vegetation type to each element, i.e., sagebrush, fir, or aspen [ , 2002 ; , 2010 ]. A single plant functional type and a vegetated fraction of 0.9 were assumed for each cell, considering that in such an environment small patches of bare soil are very likely. Physiological and structural characteristics of low sagebrush, fir, and aspen inferred from literature were used to parameterize vegetation physiological parameters [ , 2000 ; , 2007 ]. 3.2.2. Tollgate Watershed The Tollgate watershed represents the mountainous part of the RCEW and has an area of 54.8 [ km 2 ]. Its elevation ranges from 1400 to 2241 [ m a.s.l. ] ( Figure 2a ). The Tollgate watershed representing the mountainous part of the RCEW and containing the RME basin, is characterized by a vegetation cover similar to that of RME, i.e., firs occupy 11.9% of the area, aspens occupy 8.7%, and sagebrushes 79.4% ( Figure 2b ). For this watershed, two additional meteorological stations and 10 precipitation gauges located within and in the proximity to the catchment ( Figure 2a ) were used in the spatial interpolation of meteorological inputs [ , 2001 ; , 2001 ]. However, the topographic effects on wind speed and snow accumulation were neglected, given the lack of appropriate information to extend the [2002] approach to such large spatial scales. Temperature, relative humidity, and precipitation were interpolated, summing the residual and gradient fields. The gradient field was obtained from a linear regression of meteorological observations with altitude. Consequently, the residuals, i.e., the deviations of the measurements from the gradient field, were spatially interpolated using the Inverse Distance Weight (IDW) method [ , 2005 ]. This resulted in a residual field representing the local deviations from the gradient field [ , 2008 ]. The available dataset includes a period of 8 years: August 1985 through July 1993. The simulation period for the TOL catchment corresponds to the entire period of meteorological observations. The grid resolution is 50×50 m 2 , resulting in 21,924 computational elements. The description of soil properties and the representation of vegetation for the TOL catchment are analogous to that of the previously described RME catchment. 4. Model Parameter Adjustment Given the relatively high computational cost of distributed simulations and the high dimension of the parameter space of Tethys‐Chloris, no formal, traditional, or advanced calibration of the model can be carried out [ , 2002 ; , 2002 ; , 2008 ]. The choice of model parameter values has been made subjectively, based only on available data or literature information (see the companion paper of [2012] ). The model was first used for the plot scale simulations [ , 2012 ]. This permitted to control the overall consistency and reliability of simulations, using various metrics of model performance, such as the coefficient of determination, the root mean square error, the Nash‐Sutcliffe efficiency, etc. The estimated performance metrics together with qualitative considerations of plausibility of results at different temporal scales (e.g., annual, monthly and daily cycles) in comparison to other studies [e.g., , 2004 ; , 2004 ; , 2010 ] were used to adjust the model parameters, while keeping them within a realistic range. The parameters that resulted in a satisfactory performance at the plot scale were subsequently used for spatially distributed simulations. The transition from the plot scale to a distributed domain raises the question of transferability and significance of other parameters that directly affect the spatial dynamics (e.g., the conductivity anisotropy ratio, bedrock leakage). More than one attempt is usually necessary to improve the model results at the watershed scale. In the analyzed cases, the parameters adjusted at the plot scale led to a good performance when applied to spatial domains. The results presented in the following were obtained by running the model two times for the Lucky Hills watershed, where the parameter adjustments only concerned the parameterization of the soil sealing process. Five model simulations were carried out for the Reynolds Creek Mountain East watershed: major modifications concerned the hydraulic conductivity and the anisotropy ratio, as well as physiological parameters of vegetation, specifically, the root depth, the empirical parameter linking photosynthesis to stomatal conductance, a , and the maximum Rubisco capacity, V max L [ , 2012 ]. The results obtained for the Tollgate watershed were not calibrated and obtained with the same parameter values used for the Reynolds Creek Mountain East watershed. Only manual parameter adjustment was used. It implies an expert evaluation of the model performance with a subsequent modification of parameters keeping them within a physically realistic range. The overall objective was to obtain results consistent with observations or independent estimates of the simulated variables. Therefore, when compared to traditional hydrological modeling, the model performance skill relies more on the model structure rather than on the skill of a calibration procedure. Undoubtedly, uncertainties in meteorological inputs and boundary conditions limit a thorough testing of the T&C model structure. Nonetheless, the background rationale is that the structure of a physically‐based/mechanistic model should lead to a satisfactory performance with a narrow range of possible outcomes. This is possible because the physically‐based (i.e., theoretically measurable) nature of the parameters implies narrow bounds of values that they can assume, when used in the context of representing an actual physical media/ecosystem. The direct consequence is a “tighter” range of possible model results. This further implies that a mismatch between model simulations and observations can be attributed to the model structure or uncertainty of boundary conditions, rather than to a poor calibration procedure. Furthermore, the narrow feasible region of model parameters diminishes the importance of subjectivity in choosing parameter values, i.e., the expected outcome of a simulation cannot be too different as a consequence of a different parameter value. Calibration in our study, thus represents a final adjustment to refine the simulation skill, which is mainly dictated by the model structure and boundary conditions. 5. Distributed Results The results obtained with T&C for the Lucky Hills watershed and Reynolds Creek Mountain East/Tollgate watersheds are discussed in the following. 5.1. The Lucky Hills Watershed In a spatially distributed application, hydrometeorological inputs can vary among the computational elements due to topographic or local meteorological conditions. Given the small watershed area of Lucky Hills, the only input that exhibits spatial distribution is the incoming shortwave radiation (Section 3.1). Both local and remote terrain effects are accounted for in the simulation. Figure 3a shows that the inclusion of topographic effects modifies the distribution of shortwave incoming radiation. As expected, the steep slope exposed to the south receives shortwave radiation that over the long‐term is about 30% higher than that for the steepest element exposed to the north. This uneven distribution of radiation is directly reflected in the long‐term distribution of simulated transpiration fluxes ( Figure 3b ) and bare soil evaporation ( Figure 3c ). The relatively smaller available energy in the north‐facing elements leads to a lower photosynthetic activity and reduced transpiration fluxes. Related to that effect is a slight increase in bare soil evaporation (although it is fairly constant across the domain). Water availability also plays a significant role in shaping the spatial pattern of transpiration. For example, transpiration and, to some extent, bare soil evaporation are higher in the convergent part of the catchment, where higher soil water content is a result of lateral water redistribution. 3 The results of hourly spatially‐distributed ecohydrological simulations averaged over the simulation period (July 1996 through December 2009) for the Lucky Hills experimental watershed. (a) Incoming shortwave radiation. (b) Transpiration flux. (c) Bare soil evaporation flux. Water redistribution in such an environment is expected to be controlled by runon‐re‐infiltration effects [ , 2002 ; , 2003 ; , 2003 ], rather than by subsurface flows. This is confirmed by the model results. The computed lateral subsurface redistribution of water between neighboring cells is a small fraction of the annual hydrological budget ( Figure 4a ). It is generally less than 2.0 [ mm yr −1 ] and peaks in steepest slopes. The spatial differences among the mean annual infiltration rates are significant and they are a consequence of localized runon‐re‐infiltration process ( Figure 4b ). Most of water redistribution is produced by short intense events during monsoon seasons that can lead to infiltration excess runoff. The process representation is accommodated by using disaggregated rainfall at 5 [ min ] time step [ , 2005 , 2008 ; , 2009 ; , 2010 ] and accounting for the soil sealing formation [ , 2012 ]. Runoff produced after intense precipitation events is routed through the watershed drainage paths towards the outlet. Although a large fraction of generated runoff is “lost” as streamflow/overland flow, there are topographic locations with shallow slopes that facilitate reinfiltration of runon. This process is most efficient at several locations along the hollow part of the catchment, where the mean annul infiltration rates are higher than the annual precipitation ( Figure 4b ). Overall, runon subsurface flows, and evapotranspiration fluxes all contribute to creating a heterogenous distribution of mean soil water content ( Figure 4c ). The wettest parts are located along the trough, where the process of runon occurs [ , 2003 ] and the north slope, where transpiration fluxes are smaller due to lower surface irradiance. 4 The results of hourly spatially‐distributed ecohydrological simulations averaged over the simulation period (July 1996 through December 2009) for the Lucky Hills experimental watershed. (a) Total lateral subsurface flow. (b) Infiltration flux. (c) Soil moisture content. The mean canopy characteristics during the 13.5‐year simulation period are the result of vegetation adaptation to local topographic and environmental conditions: a larger LAI is simulated in the trough and a lower LAI in steep hillslopes ( Figure 5a ). Note that both the north and south facing slopes have smaller LAI values, as compared to the other parts of the catchment. This is due to the relatively drier conditions in the hillslope exposed to the south and lower levels of shortwave radiation in the north‐facing slope. Because of the fairly gentle topography and limited water redistribution, the relative spatial variability of LAI is small (≈10%). Above ground Net Primary Productivity of vegetation (ANPP) has a spatial distribution similar to that of LAI. However, it has a larger relative variability, with the maximum difference of ≈20% ( Figure 5b ), highlighting a non‐linear relationship between LAI and productivity. Simulated ANPP values of 70–85 [ gC m −2 yr −1 ] are regarded as very plausible for creosotebush in southern Arizona and match values reported in literature very well [ , 1965 ]. As can be observed in Figure 5c , the distribution of surface radiative temperature reflects the forcing of incoming shortwave radiation, with only secondary effects due to vegetation cover and soil moisture distribution (not appreciable from the figure). 5 The results of hourly spatially‐distributed ecohydrological simulations averaged over the simulation period (July 1996 through December 2009) for the Lucky Hills experimental watershed. (a) Leaf Area Index. (b) Above Ground Net Primary Productivity. (c) Surface radiative temperature, T S . The simulated dynamics of LAI and Gross Primary Production (GPP) have been compared with estimates based on remote sensing MODIS data ( Figure 6 ). A single MODIS cell contains the entire Lucky Hills catchment. The T&C model is able to reproduce the LAI seasonality and its overall magnitude quite well. The GPP is slightly overestimated with the seasonal maximum following the monsoon rainfall, especially in the relatively wetter years. The small size of the Lucky Hills watershed as compared to the resolution of MODIS vegetation products (1 km 2 ) precludes any spatial comparison between the simulated and observed patterns. However, it is important to note that the distributed application produces spatially integrated fluxes and states that are nearly identical to those simulated in the plot scale application [ , 2012 ]. The determination coefficients of the relationships between the time series obtained at the plot scale and spatially integrated LAI and GPP are 0.96 and 0.94, respectively. This is due to the small dimension of the watershed, its gentle topography, and the limited impact of surface‐subsurface flows. It is very likely that extending the analysis to a larger watershed area in the order of tens of square kilometers in this climatic zone would still provide very similar results. Consequently, for the examined case study, the plot scale simulations can be representative of a much larger area (with the apparent exception of streamflow magnitude). 6 Spatially integrated (a) Leaf Area Index, (b) Gross Primary Production, (GPP), based on the simulated (SIM.) and remote sensing estimates for the Lucky Hills watershed. The simulated annual partition of water budget among hydrological fluxes is detailed in Table 1 . Evaporation and transpiration represent the largest components of the hydrological budget. The average annual evapotranspiration is ≈310 [ mm yr −1 ], which represents 93% of the annual precipitation. Recharge to deeper soil layers is estimated to be ≈9 [ mm yr −1 ]; it is highly discontinuous during the simulation period and concentrated in a few wet periods. Infiltration excess runoff can be observed during summer months of monsoon periods due to heavy precipitation and effects of soil sealing. The streamflow at the catchment outlet is ≈23 [ mm yr −1 ]. This value is very close to the observed values for small nested basins of Walnut Gulch [ , 2008 ]. 1 The Partition Among the Hydrological Budget Components for the Lucky‐Hills Watershed for a 13.5‐Year Simulation Period (July 1996–December 2009). a Component Flux [ mm yr −1 ] P r 333.8 E G 109.3 T 191.4 E In 9.4 E snow 0.0 L K 9.0 Q 23.2 a The hydrological budget terms are: precipitation, P r , soil evaporation, E G , transpiration, T , evaporation from interception storage, E In , evaporation/sublimation of snow, E snow , deep leakage, L K , and outlet discharge, Q . The imbalance is due to a slight storage change between the beginning and the end of the simulation. The cumulative runoff simulated by the model for the last nine years of the simulation is compared with the observations in Figure 7a . There are differences between the simulated and observed streamflows, especially during the years of 2006–2007. The overall result, however, is surprisingly good. For the period of 2000 through 2009, the simulated annual runoff is 21.1 [ mm yr −1 ], while the observed runoff is 23 [ mm yr −1 ]. The error is essentially negligible, given the minimal calibration of the model, the uncertainties associated with the soil hydraulic properties, bedrock conditions, soil seal formation, and the adopted rainfall disaggregation technique. Similarly, small differences can be also observed in the flow‐duration curve ( Figure 7b ), where larger runoff values are somewhat overestimated. A very good agreement is exhibited for the durations of ≈0.5–0.8 [ day ], above which the two curves show essentially zero runoff. The simulation of runoff that is related to sporadic and intense events (typical of semiarid systems) is a challenging problem for most hydrological models. This is especially true when the objective is to truthfully simulate all of the involved physical processes. Therefore the presented simulation results can be considered to be “effectively” corroborated by the data. 7 A comparison between the observed (“OBS.”) and simulated (“SIM.”) (a) hourly cumulative runoff and (b) flow duration curves for the Lucky Hills watershed. 5.2. Reynolds Creek Mountain East Watershed The map of vegetation of the RME watershed is shown in Figure 8a . As seen, vegetation cover is rather diverse: aspen groves are located in the convergent part of the catchment and near the southeastern boundary. Patches of fir border the upstream side of the aspen grove in the center. The remaining part of the catchment is occupied by sagebrush. 5.2.1. Hydrological Partition The distribution of the top 1‐m soil moisture averaged over the 25‐year simulation period (1983 through 2008) is shown in Figure 8b . Two principal controls can be discerned, related either to soil texture or topography. Specifically, the sandy southwestern part of the catchment has consistently lower soil water contents, as compared to the rest of the domain characterized by a finer soil texture. The convergent part of the watershed has the highest soil moisture, especially in the downslope areas affected by snow drifts, where accumulation of snowpack allows the persistence of wet conditions over longer periods (see later Figure 9b ). 9 The results of spatially‐distributed ecohydrological simulations averaged over the simulation period (October 1983 through October 2008) for the Reynolds Creek, Mountain East watershed. (a) Annual transpiration flux. (b) Annual fraction of time with snow cover. The simulated lateral fluxes are relatively high due to the pronounced slopes and the high anisotropy ratio chosen for the RME watershed to mimic preferential flows associated with topography [ , 2005 ; , 2010 ]. These fluxes permit an efficacious redistribution of water that tends to converge in the central part of the watershed, ultimately leading to saturated areas downstream of the central aspen grove during melting seasons. The redistribution effect is clearly appreciable in the spatial distribution of transpiration fluxes ( Figure 9a ) that are higher in the convergent part of the catchment. The simulated spatial distribution of transpiration is due to a combined effect of moisture availability during growing season and uptake characteristics imposed for a given vegetation type. For example, taller trees, such as aspens and firs with deeper roots have access to deeper soil water and tend to transpire over longer periods with higher rates. A drastic reduction of transpiration can be observed for areas with significant snow accumulation. This abrupt decrease is partly due to the model assumption that neglects transpiration fluxes, when a cell is covered by snow, and because of a delayed start of the growing season in the snow‐covered elements. A patch of relatively small transpiration (around 300–500 m north and 400–600 m east) can be observed in Figure 9b . This is an interesting result because that area has been previously identified by field studies as a drier location [ , 2009 ; , 2011 ]. This has been now consistently confirmed by the model results. The signature of vegetation is also evident in the duration of snow cover on the ground, shown in Figure 9b . The simulated patterns are also strongly influenced by meteorological inputs. For example, snow tends to accumulate in drifts in sheltered areas (i.e., these areas correspond to the values higher than 0.6 in Figure 9b ), while precipitation tends to be less and sublimation rates are higher at wind‐exposed locations (south‐eastern part of the watershed). For the same topographic aspect and slope, aspens reduce the duration of snow cover due to snow interception effects and because of a larger aerodynamic roughness that facilitates snow sublimation. The integrated effect of firs is expressed by both enhanced turbulent exchange and a decreased incoming energy to the ground snowpack due to canopy shadow effects [see also 2006 ]. The spatial distribution of the fraction of time with snow cover agrees well with observations and model simulations previously carried out for the RME watershed by [2002] and [2002] for specific years. The mean annual partition into hydrological components is presented in Table 2 . The average simulated runoff is 543 [ mm yr −1 ], which represents about 59% of annual precipitation and compares favorably to the mean observed streamflow of 518 [ mm yr −1 ]. The latent heat flux is partitioned into bare soil evaporation, 105 [ mm yr −1 ], transpiration, 203 [ mm yr −1 ], and snow sublimation, 59 [ mm yr −1 ]. Given the small imposed bedrock conductivity, the leakage losses are negligible in the overall balance and constitute less than 1 [ mm yr −1 ]. As expected, the transpiration flux is higher than bare soil evaporation and snow sublimation. However, the magnitudes of these three fluxes are comparable in this mountain environment, implying that none can be neglected without introducing significant errors. 2 The Partition Among the Hydrological Budget Components for the TOL Watershed for a 8‐Year Simulation Period (August 1985–July 1993) and for the RME Watershed for a 25‐Year Period (October 1983–October 2008). a Component Flux (TOL) [ mm yr −1 ] Flux (RME) [ mm yr −1 ] P r 653.0 912.4 E G 100.7 105.0 T 232.6 203.8 E In 7.4 0.3 E snow 112.7 59.7 L K 0.2 0.7 Q 202.0 543.5 a The hydrological budget terms are: precipitation, P r , soil evaporation, E G , transpiration, T , evaporation from interception storage, E In , evaporation/sublimation of snow, E snow , deep leakage, L K , and outlet discharge, Q . The imbalance is due to a storage change between the beginning and the end of the simulation. 5.2.2. Streamflow The overall consistent and representative performance of T&C is further corroborated by the simulated discharge at the outlet of the RME watershed, as shown in Figure 10 . The cumulative discharge is reproduced with a sufficiently high accuracy over the entire period ( Figure 10a ). A marginal deviation occurs in the last 7 years. This is also testified by an agreement between the simulated and observed annual runoff for each water‐year (October through September) of the simulation ( Figure 10b ). The global Nash‐Sutcliffe efficiency, (NS), for the hourly discharge values is 0.59, and the coefficient of determination between the observed and simulated series is R 2 = 0.65. The result is somewhat poorer in terms of the flow‐duration curve ( Figure 10c ). The model is able to reproduce the higher discharge rates with the proper duration but tends to underestimate the flow minima. For instance, for durations larger than 100 days, the simulated discharge is typically lower than the observed one. In other words, the simulated flow rates exhibit dampening of the recession curve, while observations suggest persisting flows. 10 A comparison between the observed and simulated (a) cumulative runoff, (b) water‐year values of runoff and precipitation, and (c) flow duration curves for the Reynolds Creek, Mountain East watershed. Q obs and Q sim denote the observed and simulated annual discharge, respectively; P r denotes the annual precipitation averaged over the watershed. The simulated discharge series are illustrated in Figure 11 for each of the 25 water‐years. In general, the discharge is simulated in a realistic fashion in most water‐years but the overall quality of simulations changes from year to year. This is not surprising, given the possible uncertainty of the input data such as snow, or the capability of the model to better reproduce hydrological conditions of a given year. In quantitative terms, the model skill for the water‐year 1996 is extremely good, ( NS = 0.83), while for the water‐year 2004 is very poor ( NS = −0.25). There is no a clear time‐evolution of the model performance and a negative efficiency is obtained only in 3 out of 25 years; for most years, NS >0.5. It can be concluded that the simulation skill reflects the physically‐based nature of represented processes responding in a consistent fashion to hydrometeorological conditions of a given year. 11 A comparison between the hourly observed (“OBS.”) and simulated (“SIM.”) discharge at the RME watershed outlet for each of the 25 water‐years (October to September). 5.3. Tollgate Watershed Spatially distributed results for the TOL watershed were averaged over the 8‐year simulation period of August 1985 through July 1993, that corresponds to the period of data availability for the entire watershed (Section 3.2.2). Over this period, MODIS data were not yet available. Remote sensing estimates of snow cover, LAI, and GPP (MODIS products: MOD10A2, MOD15A2, and MOD17A2), for the period of January 2001 through December 2010 were used to confirm the model performance. Given the two different periods (i.e., the period of watershed data availability and the period of remote sensing observations), results can only be compared in terms of average spatial patterns and seasonal dynamics under the assumptions of climate and vegetation stationarity. While the assumption of vegetation stationarity can be considered rather acceptable for a fifteen year temporal lag (1985–2011), certain climatic trends have been observed for the Reynolds Creek watershed [ , 2010 ]. Specifically, when climate metrics for the periods of October 1985 through October 1993 and October 2000 through October 2008 are compared for the RME stations (the only stations for which data were available for the period of 1983–2008, see Section 3.2.2), an increase of air temperature by 0.16 [° C ] and precipitation by 122 [ mm yr −1 ] is observed over the 2000–2008 period. Therefore, in a comparison of T&C results and MODIS observations, we need to qualitatively account for the fact that the period of MODIS operation corresponds to a somewhat warmer and wetter climate (see also Figure 10 ). In order to present a “mechanistic” insight on an internal structure of the watershed response, model performance is also illustrated in terms of spatial patterns for bare soil evaporation, transpiration, snow sublimation, and effective saturation. While these comparisons are only qualitative, they further corroborate the effective consistency of T&C in reproducing realistic dynamics. They further argue for the need of detailed spatially‐distributed observations of states and fluxes of a land‐surface for a rigorous model confirmation. 5.3.1. Hydrological Partition Figure 12a shows the simulated distribution of bare soil evaporation across the TOL watershed. Bare soil evaporation typically ranges between 60 and 150 [ mm yr −1 ], with higher values in the convergent topographic areas or the areas exposed to stronger winds; these two corresponding controls are the emerging features of the simulated pattern. In the northwestern part of the catchment, bare soil evaporation tends to be lower. This is the result of enhanced activity of vegetation ( Figure 12b ) that decreases the amount of water available for soil evaporation. The distribution of transpiration flux in Figure 12b shows that larger rates are typical for areas covered with fir but also that moisture distribution exerts an important control on this flux. Specifically, higher transpiration rates in the western and northern parts of the catchment are the result of a larger amount of precipitation [ , 2000 ]. Being wetter, the western part of the catchment favors vegetation covers and an adaptation of taller vegetation types ( Figure 2b ). The control of soil moisture is clearly appreciable in the downstream part of the watershed: the simulated transpiration fluxes are different in the western and eastern sides of the catchment, despite a uniform sagebrush fractional cover. Furthermore, a north‐south elevation gradient of transpiration can be appreciated for areas vegetated with sagebrush. Low elevations that have shorter snow cover periods and higher temperatures favor photosynthetic activity and, consequently, transpiration. 12 The results of spatially‐distributed ecohydrological simulations averaged over the simulation period (August 1985 through July 1993) for the Tollgate watershed. (a) Mean annual bare soil evaporation flux. (b) Mean annual transpiration flux. Figure 13a illustrates the simulated average fraction of time that each computational element is covered by snow. As already pointed out for the RME watershed, the signature of vegetation is the predominant feature of the simulated spatial pattern. The canopy shadow effect exerted by vegetation on snow is not as significant as the impact of interception and higher roughness resulting in larger snow evaporation/sublimation in aspen and fir areas (not explicitly presented). The total snowpack that accumulates on the ground below aspens and firs is smaller, when compared to open areas with sagebrush. Although melting of undercanopy snow is slower due to the less available energy, the duration of snow season is shorter because of the significantly reduced snow water equivalent. When compared with MODIS observations the simulated patterns of temporal fraction of snow cover are very consistent ( Figures 13a and 13b ). The major appreciable differences are at low watershed elevations, where T&C simulates a longer snow cover, and in the fir‐aspen vegetated areas at high elevations. The differences at the lower elevations can be the result of aforementioned warmer climate conditions corresponding to the MODIS observation period. The lack of a vegetation pattern in the data might be attributed to the MODIS resolution (about 500 m×500 m), which is clearly insufficient to capture the smaller‐scale variability due to topographic features and vegetation effects. Despite these differences the fraction of time with snow cover averaged over the watershed is 0.35 in T&C simulations, and 0.37 in MODIS observations, a surprisingly good match. A similarly successful comparison is also obtained in reproducing the time dynamics of snow cover fraction: Figure 13c shows the space‐time average seasonal dynamics over the Tollgate watershed. 13 The fraction of time with snow cover for the Tollgate watershed. (a) The simulated pattern over the period of August 1985 through July 1993; (b) remote sensing estimations over the period of January 2001 through December 2010; (c) average seasonal cycles from remote sensing estimates and simulations. The signature of vegetation distribution is also evident in the patterns of snow evaporation/sublimation ( Figure 14a ). Taller plants significantly increase the rates of snow sublimation because of snow interception by canopies and turbulence enhancement (due to the larger roughness of vegetated surfaces). This is expressed especially well for fir‐vegetated areas. Therefore, vegetation has a very important role in controlling the amount and dynamics of snowpack, exerting both positive (radiation interception) and negative (snow interception, turbulence enhancement) feedbacks. The net outcome, however, is dependent on both vegetation structure and local meteorological conditions, especially wind speed. 14 The results of spatially‐distributed ecohydrological simulations averaged over the simulation period (August 1985 through July 1993). (a) Mean annual snow evaporation/sublimation. (b) Effective saturation. The two or three orders of magnitude of difference of the characteristic roughness of snow‐covered sagebrush vs. aspen/fir grove areas lead to drastic differences in snow evaporation rates. They are typically 2–4 times larger for aspen/fir areas: 70–100 [ mm yr −1 ] for sagebrush areas and 200–300 [ mm yr −1 ] for aspen/firs. In conditions of lower wind speeds (for instance in the northern valleys of the basin), the difference in sublimation between the two land cover types is smaller, e.g., 40 vs. 60 [ mm yr −1 ]. Obviously, these results are affected by large uncertainties, such as the representation of the within‐plant turbulence effects or the spatial distribution of wind speed fields because of a limited number of meteorological stations. Nonetheless, the simulated snow evaporation/sublimation rates are comparable to those obtained in other modeling studies for mountainous areas [ , 2008 ; , 2010 ]. The spatial distribution of effective saturation clearly highlights that the northeastern part of the watershed is relatively drier than the rest of the basin. It can also be inferred that vegetation, especially firs, create patches of lower water content ( Figure 14b ). The presence of firs near the channel network leads to a somewhat smaller soil moisture in the elements containing channels. The areas of lowest effective saturation are due to a combination of soil texture and vegetation related effects. Specifically, these areas are typical for coarse sandy soils with higher hydraulic conductivities and covered by firs. The coarse soil texture facilitates moisture uptake by deep roots of tall vegetation. Simulations results indicate that subsurface lateral flows are the major source of water redistribution in this watershed and they typically peak in steepest slopes that are in the southern part of the watershed. The annual partition among the simulated hydrological components is presented in Table 2 . The average simulated runoff is 202 [ mm yr −1 ], this represents ≈31% of annual precipitation and compares well with the observed discharge that is 188 [ mm yr −1 ]. The magnitude of bare soil evaporation is 100 [ mm yr −1 ], and of transpiration is 232 [ mm yr −1 ]. The rate of snow evaporation is quite high, 112 [ mm yr −1 ], due to a high average interpolated wind speed in the TOL watershed, i.e., 3.47 [ m s −1 ] (this is relative to the RME watershed, where the average value is 2.94 [ m s −1 ]), and a larger areal fraction covered by evergreens (Section 3.2.2). 5.3.2. Vegetation Metrics The spatial distribution of vegetation Leaf Area Index (LAI) is presented in Figure 15 . The simulated mean LAI reflects adaptation of vegetation to long‐term environmental conditions, i.e., although at the scale of each computational element the vegetative fraction is assumed to be constant (equal to 0.9), soil moisture availability, radiation, air temperature, etc., all contribute to the canopy dynamics and therefore to the spatial distribution of average LAI ( Figure 15a ). A comparison between the simulated and the MODIS derived LAI ( Figure 15b ) is significantly hampered by the MODIS resolution for LAI (1 km×1 km). Nonetheless, qualitatively similar spatial distribution and magnitude of LAI values can be observed in Figure 15 . Specifically, lower LAI in the northeastern part and higher LAI in the central‐southwest part of the watershed (≈1–2 kilometers from the divide) can be inferred. The catchment simulated average LAI is 0.61 and the MODIS‐derived LAI is 0.66. 15 The Leaf Area Index for the Tollgate watershed. (a) The mean annual simulated pattern over the period of August 1985 through July 1993; (b) the mean annual LAI based on remote sensing estimations over the period of January 2001 through December 2010; (c) average seasonal cycles from remote sensing estimates and simulations. The highest LAI values are simulated for evergreen trees since they do not shed leaves. Firs have values of LAI ≈2.5 in the wettest, southeastern part of the catchment, and ≈1.5 in the drier northern part. Conversely, aspens exhibit a pattern of LAI that is inversely related to the elevation, e.g., the multi‐annual average values of LAI are around 0.5–0.6 at higher elevations of the southern part and 0.8–0.9 at some locations in the northern part of the catchment. Since aspens are deciduous trees, their summer LAI is typically 2–2.5 higher than the annual average. Lastly, the LAI of sagebrush is influenced by many factors, such as elevation (air temperature), moisture availability, and soil texture. The sandy northwestern part of the catchment exhibits LAI magnitudes around 0.5–0.9, conversely, the northeastern part and some of the locations closest to the mountainous divides or in topographically shaded areas exhibit the lowest values of 0.25–0.4. Similarly to aspens, the summer peaks of sagebrush LAI are typically 1.5–2 higher than the annual average. The simulated LAI peaks compares surprisingly well with field observations of maximum LAI carried out by [2010] in the RME watershed for the three represented vegetation types (see Section 3.2.2). The average seasonal dynamics of LAI over the Tollgate watershed are shown in Figure 15c , based on simulations and MODIS data. The simulated LAI dynamics seemingly overestimate winter LAI, which is simply an artifact due to the presence of evergreen firs and shrubs. This “discrepancy”, however, can be due to poor reliability of MODIS estimates, during periods when vegetation is totally or partially buried by snow (the case for the Tollgate watershed). The difference during the summer peak can be related to the generally wetter and warmer conditions of the MODIS observational period. Further data would be necessary to corroborate such a statement. Note also that the differences are of the order of 0.1–0.3 LAI, which is below typical uncertainties associated with LAI estimation, even when it is made from the ground [ , 2003 ]. The spatial pattern of the long‐term average Gross Primary Production (GPP) is shown in Figure 16 . As seen, features present in the distribution of LAI can be also appreciated in the distribution of GPP. The spatial organization of GPP is even more pronounced than that of LAI. Evergreen firs have GPP mainly in the range of ≈350–650 [ gC m −2 yr −1 ], with a spatially average value of 485 [ gC m −2 yr −1 ]; deciduous aspens have a somewhat lower productivity with GPP of ≈250–450 [ gC m −2 yr −1 ], and a spatially average of 329 [ gC m −2 yr −1 ]. Sagebrush plants exhibit a wider range of environmental conditions with the corresponding variability of GPP in the range of 250–550 [ gC m −2 yr −1 ] and a spatially average value of 236 [ gC m −2 yr −1 ]. Remote sensing estimates mainly confirm these spatial results, except for the northwestern part of the watershed where a higher productivity is simulated ( Figure 16b ). The available resolution of MODIS product (1 km×1 km) prevents further spatial confirmation of the modeling results. However, MODIS inferred average annual GPP (for the period 2001–2010) is 270 [ gC m −2 yr −1 ]; the corresponding simulated value averaged over the TOL watershed is 274 [ gC m −2 yr −1 ]. The average seasonal dynamics of GPP are also well captured as demonstrated in Figure 16c . The two average yearly cycles are quite similar during most of the season with the difference during late summer and early fall, where MODIS estimates are double that of the simulated results. This difference can probably be attributed to the warmer and wetter climate over the period of MODIS deployment. A possible outcome of such conditions is a delayed onset of summer drought and lengthening of the growing season. Note how soil moisture limitation effects are pronounced in GPP during the second part of the year both in remote sensing estimates and simulations. MODIS inferences are not sufficiently detailed to confirm the productivity for the three individual PFTs. However, the simulated values are in the expected range for such an ecosystem [ , 2003 ; , 2005 ; , 2011 ]. 16 The Gross Primary Production for the Tollgate watershed. (a) The mean annual simulated pattern over the period of August 1985 through July 1993; (b) the mean annual GPP based on remote sensing estimations over the period of January 2001 through December 2010; (c) average seasonal cycles from remote sensing estimates and simulations. 5.3.3. Streamflow and Snow Depth The time series of streamflow at the outlet of the TOL watershed are shown in Figure 17 , where the runoff rates ( Figure 17a ) and the cumulative discharge ( Figure 17b ) are illustrated. While the total runoff rate is predicted satisfactorily, the model tends to overestimate the discharge in the early part of the melting season and, conversely, to underestimate it later in the summer ( Figure 17b ). Two main reasons can create such a mismatch: (i) the absence of wind induced snow re‐distribution in the model simulation for the TOL watershed; and (ii) the uniform soil depth of 1 m assumed for the entire domain. Snow redistribution can delay snow melt, creating areas with large amounts of accumulated snow. A variable soil depth can produce deep storages of water at the beginning of the melting season that can feed the stream later during the summer. Such processes are not accounted for in the simulation. However, the overall coefficient of determination computed using hourly data of the simulated period is R 2 = 0.59, which can be considered as an acceptable result, given that no effort was applied to calibrate the model. 17 A comparison between the observed and simulated series for the Tollgate watershed (a) hourly outlet discharge and (b) cumulative runoff. Finally, Figure 18 presents an explicit confirmation of the distributed simulation using observed data on snow water equivalent collected at different sites of the catchment ( Figure 2a ). The comparison was carried out for four locations. Note that confirmation is entirely “blind”, implying no effort was applied to tune the simulation results. The model is able to reproduce satisfactorily the distributed dynamics of snow water equivalent, although it appears to underestimate the ground snowpack in years characterized by smaller snow precipitation. The model is spatially consistent across different precipitation conditions (note the differences in scale of the ordinate axes in Figure 18 ). The observed underestimation of T&C can be related to the model structure but it can be also affected by the spatial interpolation of precipitation and wind speed, or by the fact that snow water equivalent measurements are typically collected at sheltered sites, where snow drift accumulation might occur. 18 A comparison between the observed (“OBS.”) and simulated (“SIM.”) snow water equivalent at four monitoring locations (UTM zone 11, coordinates: first location: x = 515864, y = 4771970, second location x = 514041, y = 4769438, third location x = 521613, y = 4769718, fourth location x = 520055, y = 4768117) within the Tollgate watershed. 6. Discussion and Conclusions The modeling approach presented here relies on essential principles of land‐surface physics and life‐cycle vegetation processes. It allows one to reasonably approximate the entire range of system dynamics without summarizing information in conceptual, spatially lumped variables or using assumptions that distort the first principles of physics. Despite certain simplifications, the model provides results that can be confirmed with a large spectrum of metrics that are theoretically observable at different scales [see also , 2012 ]. The result of simulations conducted so far suggest that the model satisfactorily reproduces principal mechanisms that control the system's response. Since data required for a detailed validation/verification of a model such T&C are typically unavailable at the landscape scale of application, a community effort is warranted. This study further advocates [ , 2004 ; , 2006 ; , 2008a ; , 2011 ] the need to develop methods for collecting data on multiple variables at different spatial and temporal scales. This will facilitate the characterization of spatial‐temporal processes of hydrology and vegetation dynamics and provide better constraints of model performance. This should lead to a more robust confirmation of a new generation of Earth‐system models [ , 2006 ; , 2010 ]. At present, the scarcity of interdisciplinary data makes it difficult or sometimes even impossible to test all of the desired behavioral aspects of the model. For example, see the discussion on snow dynamics below vegetation canopies presented in this study. In this regard, the model can be also used as an exploratory tool to design interdisciplinary field campaigns. The capability of Tethys‐Chloris to produce consistent results in terms of many hydrological and ecological metrics has been demonstrated to a degree allowed by available data for two very different environments at the watershed scale. The performance obtained for the Lucky Hills site, a desert shrub ecosystem, and for the Reynolds Creek watershed exhibiting a complex ecohydrological system with different soil textures, vegetation types, and climatic gradients is considered to be satisfactory. Although, we cannot solely rely on discharge to confirm a complex spatially distributed ecohydrological model, simulation efficiencies obtained in this study with no, or limited calibration were comparable to those of calibrated conceptual models. This can be regarded as a significant result. Further, the satisfactory performance for the RME watershed for a 25‐year simulation with static parameters and a lack of any calibration is considered to be a noteworthy result. Generally, all of the simulations are obtained without significant (or, as presently popular, “automated”) calibration efforts, despite the large phase‐space dimensionality of the model. The satisfactory performance is thus a consequence of a significant investment into a physically‐realistic structure of the model, and information content used by the model for the analyzed case studies. The latter includes both a priori information on the feasible range of model parameters and detailed meteorological and catchment hydrological data. The aforementioned shortcomings, such as the limited skill in reproducing the recession curve in the Reynolds Creek subwatersheds, is likely due to the assumptions regarding the soil depth and leakage into the fractured bedrock aquifer. Presently, these are nearly impossible to verify. The shallow soil mantle used in the simulations cannot capture water table dynamics at deeper locations and their effect on streamflow [ , 2011 ]. Nevertheless, we should note that very small flow rates, around 1–10 [ l s −1 ], are typically difficult to capture with any deterministic model. A synthesis with stochastic/deterministic models of flow in fractured media [ , 2002 ] is therefore warranted. Undoubtedly, the model performance for the larger watershed (TOL), especially in terms of discharge, can be ameliorated to better match the observations. The small amount of information available in terms of inputs and boundary conditions is likely to be one of the reasons of poorer model results in terms of streamflow. Nonetheless, we speculate that while the revealed patterns do not necessarily reproduce the actual ones with a high precision, they are highly plausible, given the imposed boundary conditions and the available data. For example, the corroboration between the basin‐averaged LAI, GPP, fractional snow cover, and inferences based on MODIS data for the Tollgate watershed are a surprisingly good confirmation of the simulated vegetation and snow dynamics. Note that LAI dynamics emerge only from the imposed vegetation physiological parameters and environmental conditions. Furthermore, some of the discussed findings represent interesting insights into the quantitative analysis of the dominant hydrological processes at the landscape/watershed scale. For instance, the role of lateral subsurface flow has been found to be unimportant for water redistribution, as compared to the surface runoff‐runon mechanism, in the semiarid Lucky‐Hills watershed. However, it appears to be a very important mechanism for redistributing water in the Reynolds Creek watershed. A peculiar insight is that vegetation imposes an important signature on the snowpack dynamics by modifying precipitation, radiative energy fluxes, and turbulent exchange. Micrometeorological conditions and plant structure can have a significant role in the overall water budget at the watershed scale. This implies that spatial vegetation processes, such as tree encroachment in mountainous areas [e.g., , 2007 , 2008 ], could significantly affect water fluxes into the soil and, consequently, streamflow and deep recharge, in the long‐term influencing aquifers and water budget of these areas. Acknowledgments Simone Fatichi wishes to thank the support of the International joint Ph.D. program on “Mitigation of Risk due to Natural Hazards on Structures and Infrastructures” between the University of Firenze (Italy) and the T.U. Braunschweig (Germany). Authors are extremely grateful to all the persons involved in the collection of data and information in the two analyzed locations, Lucky‐Hills (AZ) and Reynolds Creek (ID). The datasets for Lucky Hills were provided by the USDA‐ARS Southwest Watershed Research Center. Funding for these datasets was provided by the United States Department of Agriculture, Agricultural Research Service. The datasets for Reynolds Creek were provided by USDA‐ARS Northwest Watershed Research Center. The technical support of the Center for Advanced Computing at the University of Michigan is also acknowledged. Valeriy Yu. Ivanov was partially supported through the NSF grant 0911444. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Journal of Advances in Modeling Earth Systems Wiley

A mechanistic ecohydrological model to investigate complex interactions in cold and warm water‐controlled environments: 2. Spatiotemporal analyses

Loading next page...
 
/lp/wiley/a-mechanistic-ecohydrological-model-to-investigate-complex-5DjQZ0ftCy

References (128)

Publisher
Wiley
Copyright
© 2012 by the Chinese Geophysical Society
ISSN
1942-2466
eISSN
1942-2466
DOI
10.1029/2011MS000087
Publisher site
See Article on Publisher Site

Abstract

1. Introduction The studies of dynamics of hydrological and vegetation processes have been historically more focused on the understanding of mechanisms and controls at the point/plot spatial scale rather than at the scale of topographically complex landscapes. At the plot scale, numerous theoretical [e.g., , 1999 ; , 2001 ; , 2005 ], field experimental work [e.g., , 2004 ; , 2004 , 2010 ; , 2010 ], and novel modeling approaches [ , 2008 , 2009 ; , 2010 ] have been presented. Recently, there was an emergence of interest in extending ecohydrological analysis to larger spatial scales in complex landscapes [ , 2008b ; , 2008 , 2010 ], where the two‐way coupling between vegetation dynamics and water cycle is far less studied. For instance, this is the case for mountainous environments, where topographic, vegetation, and climatic gradients play a crucial role in controlling hydrological, ecological, and biogeochemical behaviors [ , 2008 ; , 2008 ; , 2009 ; , 2009 ; , 2011 ]. The effects of a spatially varying structure of partition and magnitude of the incoming radiative energy, the role of water redistribution, and other terrain‐related hydrological processes are much less understood at the watershed scale than at the local scale [ , 2010 ]. This can be explained by the spatially and temporally heterogeneous conditions of watersheds and the difficulty of observing or accounting for such differences in a parsimonious manner in models. A distributed representation of ecohydrological dynamics and a subsequent aggregation to the watershed scale are thus warranted. Both analytical [ , 2004 , 2005 ] and numerical methods [ , 2001 ; , 2008b ; , 2009 ; , 2010 ] have been proposed to extend ecohydrological analysis to larger spatial scales. In this study, a numerical approach is pursued using Tethys‐Chloris, a novel ecohydrological model presented in the companion paper [ , 2012 ]. Tethys‐Chloris (T&C) is an ecohydrological model that reproduces all essential components of hydrological cycle, resolving the mass and energy budgets at the hourly scale. It includes a description of energy and mass exchanges in the atmospheric surface layer accounting for up to two layers of vegetation, a module of saturated and unsaturated soil water dynamics, and a snowpack evolution module. Vegetation dynamics are also reproduced and include all essential plant life‐cycle processes, i.e., photosynthesis, phenology, carbon allocation, and tissue turnover [ , 2012 ]. The first confirmation of the model was presented in the companion paper for two plot‐scale case studies [ , 2012 ]. In this work, the analysis is extended to the watershed scale. The two cases are the Lucky Hills experimental watershed in Arizona (U.S.A.) and the Tollgate/Reynolds Creek Mountain East watersheds in Idaho (U.S.A.). Both represent semiarid environments but while the former is characteristic of a partially vegetated desert shrub system, the latter is vegetated with sagebrush, firs, and aspens and characterized by a seasonal occurrence of snow. Semiarid ecosystems represent a significant fraction of global terrestrial surface area [ , 2005 ] and they have an important role for the Earth's climate [ , 2008 ]. The model application for semiarid environments is also consistent with the assumption of water‐limited vegetation, since nutrient dynamics are currently neglected in Tethys‐Chloris. The model is applied at the distributed spatial scale with the overall objective to assess the capability of reproducing physical and biophysical metrics of dynamics. These are energy fluxes, soil moisture, snowpack, vegetation production, and seasonality of Leaf Area Index (LAI). A distributed analysis at the spatial scale of a watershed presents numerous issues for a confirmation of any mechanistic ecohydrological model. They arise due to the common lack of spatially distributed observations at the appropriate scales for most of the simulated variables, and because of the high uncertainty in specifying the boundary conditions for a system, e.g., meteorological inputs, bedrock depth distribution, and vegetation spatial variations. In the case of Lucky Hills basin and the Tollgate/Reynolds Creek Mountain East watersheds, we identified two suitable catchment systems that have long‐term, quality‐controlled, diverse (energy flux, snowpack, streamflow, etc.) records. Nonetheless, the spatially distributed comparison for the Lucky Hills basin remains severely limited due to its small size that does not allow a sound comparison with remote sensing products. The present paucity of spatially‐distributed data permits only a limited evaluation using remote sensing observations, with more emphasis on an overall qualitative consistency of spatial patterns rather than their explicit quantitative confirmation. The plausibility of simulated variables and insights on the mechanisms of simulated dynamics are also discussed, as indirect metrics of model confirmation. Despite limited calibration (or virtually none as the total number of test runs was less than a dozen) of Tethys‐Chloris, the partial confirmation of the model yielded satisfactorily results. This further encourages the development and improvement of fully‐coupled ecohydrological models for applications in a spatially distributed fashion. The possibility offered by such models for studies of internal functioning of a river basin will be useful for designing virtual experiments, testing hypotheses, and focusing questions of scientific inquiry. Promising directions of research that current applications of Tethys‐Chloris open for further investigation are described in the analysis of results and the concluding discussion. 2. Issues of Model Confirmation Spatially distributed data are typically available for a limited range of metrics over small areas. Specifically, while streamflow data are available for numerous watersheds, comprehensive data on energy flux, soil moisture, snow, and vegetation productivity, exist only for small‐scale experimental watersheds (10 −3 –10 1 km 2 ) and are at a coarse resolution or entirely unavailable at larger scales, as also demonstrated in the presented examples. In this regard, this work further advocates the urgent need for novel distributed data to better validate complex mechanistic models [ , 2011 ; , 2011 ; , 2012 ]. An assessment of performance skill of ecohydrological models of “new generation” [ , 2008b ; , 2009 ] cannot rely only on aggregated metrics, such as streamflow discharge. In addition to a spatial extent, an evaluation of a mechanistic ecohydrological model requires a diverse suit of observational data. Unfortunately, they are rarely captured by a single experimental field design. For example, the interdisciplinary nature of T&C applications requires meteorological, hydrological, vegetation productivity, energy and carbon fluxes measurements to be co‐located and collected at the same period(s) of time. Furthermore, vegetation physiological and structural attributes, as well as soil texture profiles have to be known. The conventional scarcity of interdisciplinary data makes it difficult, or sometimes even impossible, testing all of the desired behavioral aspects of process‐based models [e.g., , 2008a ]. A significant contribution to the problem of meaningful confirmation of mechanistic ecohydrological models has been offered by the “FLUXNET” monitoring program ( www.fluxnet.ornl.govfluxnetindex.cfm ) [ , 2001 ; , 2007 ]. Further progress in unification of measurement protocols and comprehensiveness of collected data will be addressed by the “NEON” network [ , 2008 ] Another promising source of information is offered by remote sensing data, such retrievals from the sensor Moderate Resolution Imaging Spectroradiometer (MODIS) ( http:modis.gsfc.nasa.gov ) [ , 1998 ; , 1998 ; , 2002 ]. Remote sensing observations can be used to test the capability of reproducing structural properties of vegetation canopy, such as Leaf Area Index (LAI) and vegetation phenology. In general, combining observational data from flux towers with measurements of subsurface states (soil moisture) and vegetation attributes allows one to effectively evaluate individual model components as well as the model overall performance. Although the data contributions provided by FLUXNET and MODIS databases are significant, the demand for multidisciplinary, spatially distributed data for testing and improvement of Earth‐science models exceeds data availability (see also discussion in [2011] ). There is a clear need to quantitatively check whether model structures are suitable for representing distributed ecohydrological dynamics across a range of spatial and temporal scales. The generation of appropriate data would require an effort of the entire community and experimental watersheds with multi‐scale and long‐term monitoring plans. In the presented study, the objective of confirmation is to underline that T&C behaves consistently across different scales, as can be inferred from the available data. The term “confirmation” used here is generally preferred to “verification”, “validation”, or “corroboration”. The first two processes are inherently impossible for numerical models of open natural systems [ , 1994 ]. Models can only be evaluated in relative terms and confirmed by a demonstration of agreement between observations and predictions. This keeps the issue of development of more rigorous methods of assessment open [ , 1994 ]. Despite a careful analysis of the model performance and consistency of simulations carried out in this study, some of the model results will remain critical to confirm further. For instance, this is the case for the coupled modeling of vegetation, radiation transfer, and snowpack dynamics. Detailed field experiments to confirm these dynamics are rarely available [ , 2008 ; , 2008 ; , 2009 ], although the scope and aims of intercomparison projects, such as the recent SNOWMIP2 [ , 2009 ; , 2009 ; , 2010 ], provide a promising framework for testing these process interactions at the plot scale. A number of uncertainties can be also identified in the modeling of vegetation components, such as the carbon allocation scheme, the process of tissue turnover, and drought effects on stomatal closure and plant mortality. Such crucial aspects are only partially testable with data [ , 2006 ; , 2009 ; , 2009 ; , 2009 ] and a substantial progress is required to further their ecological and physiological understanding [e.g., , 2001 ; , 2002 ; , 2003 ; , 2007 ; , 2008 ; , 2010 ]. 3. Case Study Two case studies characterized by different vegetation types and climatic conditions were simulated using T&C. The first case study is the Lucky Hills experimental watershed in Arizona, U.S.A. (Section 3.1). The second case study is represented by two nested basins located in the mountainous area of the Reynolds Creek experimental watershed (Idaho, U.S.A.), specifically, the Tollgate and Reynolds Creek Mountain East basins (Section 3.2.2). 3.1. Lucky Hills Watershed The Lucky Hills monitoring site (110.30 W , 31.44 N ; elevation 1372 [ m a.s.l. ]) is located within the Walnut Gulch Experimental Watershed, near Tombstone in the south‐eastern Arizona, U.S.A. Walnut Gulch is a long‐term experimental watershed managed by the Southwest United States Department of Agriculture‐Agriculture Research Service (USDA‐ARS), where research on hydrology and soil erosion started in 1953 [ , 2008 ]. The studied ecosystem is a plant community of desert shrubs. The climate at the Lucky Hills site can be classified as semiarid, hot, with dry winter. The mean annual temperature from meteorological observations (July 1996 through December 2009) is 17.2 [° C ] and the mean annual precipitation is 333 [ mm ] [ , 2008 ]. The Lucky Hills experimental watershed has an area of 0.037 [ km 2 ]. The elevation range is about 10 [ m ] ( Figure 1a ). This small watershed has shallow to moderate slopes of <0.1 [−], except for the central part, where steep fragments of terrain can have slopes of ≈0.4 [−] ( Figure 1b ). 1 A representation of topographic attributes of the Lucky Hills experimental watershed. (a) Digital Elevation Model. (b) Slope fraction [−] calculated with the steepest descent method [ , 2008 ]. Long‐term measurements of runoff and sediment transport have been collected at Lucky Hills over the period of 1963 through 2008 [ , 2007 ; , 2008 ; , 2008 ]. The mean long‐term runoff is estimated to be about 18.2 [ mm yr −1 ] and the mean long‐term sediment transport is 0.42 [ kg m −2 yr −1 ] [ , 2007 ; , 2008 ]. In order to study carbon dioxide and water fluxes in the Walnut Gulch area, two flux towers have been in operation since 1997 [ , 2008 ]. One of the flux towers is located in the Lucky Hills catchment and is a part of the FLUXNET network. The flux tower footprint includes a shrub plant community, mainly composed of evergreen shrubs, such as creosotebush ( Larrea tridentata ), tarbush ( Flourensia cernua ), and deciduous shrubs, such as whitethorn acacia ( Acacia constricta ) [ , 2008 ; , 2008 ]. Average vegetation height is around 0.3–0.5 [ m ] and vegetation cover fraction is approximately 0.25–0.4 [ , 1994 ; , 2001 ]. The soil at this site is coarse loam [ , 2008 ]. The surface horizon (0–6 [ cm ]) contains 650 [ g kg −1 ] of sand, 290 [ g kg −1 ] of silt, and 60 [ g kg −1 ] of clay, with 290 [ g kg −1 ] of coarse fragments larger than 2 [ mm ], 8 [ g kg −1 ] of organic carbon, and 21 [ g kg −1 ] of inorganic carbon [ , 2008 ]. The saturated conductivity has been observed to consistently diminish exponentially with the soil depth [ , 2000 ]. Tethys‐Chloris was used to simulate ecohydrological dynamics that occurred in the Lucky Hills basin over the period of July 1996 through December 2009, or 13.5 years in total. The grid resolution was of 10×10 m 2 with the total number of basic computational elements equal to 371. A two meter soil column with a free drainage boundary condition at the bottom was assumed for all computational elements composing the watershed. The vertical mesh consisted of 18 soil layers with increasing thickness with depth, starting with 10 [ mm ] at the soil surface, and ending with 400 [ mm ] at 2m depth. Spatially homogenous soil hydraulic properties were assumed; they were derived from the pedotransfer functions of [2006] using a 0.75 fraction of sand and 0.10 fraction of clay. The saturated hydraulic conductivity, K v,s , was parameterized to decay with soil depth, K v,s ( z ) = K v,s (0) e (−0.0011 z ) , where z [ mm ] is the soil depth [ , 2000 ], and the soil anisotropy factor was assumed to be a r = 1 [ , 2012 ]. Spatially uniform crown area fractions, C crown , were assigned to be 0.25 for deciduous shrubs, and 0.10 for evergreen shrubs (see [2012] for the definition of C crown ). Physiological and structural characteristics of whitethorn acacia and creosote bush derived from a literature survey were used in the choice of the parameters and properties of deciduous and evergreen shrubs [ , 1986 ; , 1990 ; , 1994 ; , 2006 ; , 2008 ; , 2009 ]. 3.2. Reynolds Creek, Idaho The Reynolds Creek Experimental Watershed (RCEW) managed by the USDA‐ARS Northwest Research center is located in the Owyhee Mountains 80 km southwest of Boise, ID, U.S.A. [ , 2001 ]. The RCEW area was described in a series of papers by [2001] , [2009] , [2010] , and [2011] . The collected dataset is particularly appealing for hydrologic modeling because of high quality, high frequency, and long‐term measurements of climatic and hydrological variables, including snow water equivalent, snow depth, and soil moisture [ , 2001 ; , 2001 ; , 2001 ; , 2001b ; , 2001 ; , 2002 ; , 2002 ; , 2010 ; , 2011 ]. Energy and mass fluxes have recently started being collected [ , 2010 ]. The RCEW domain encompasses an area of 239 [ km 2 ], with an elevation range of about 1200 [ m ] (1101–2241 [ m a.s.l. ]), where a strong topographic precipitation gradient (236–1123 [ mm ]) exists. A strong precipitation seasonality is characteristic of the area: 70% of precipitation occurs during November through May period [ , 2001 ; , 2002 ]. Snow is the dominant form of precipitation, especially in the headwater part of the catchment. All details on data collection can be found in the above references; only a succinct description is presented in the following. The T&C model was used to simulate two nested sub‐watersheds of RCEW ( Figure 2 ): the Reynolds Creek, Mountain East (RME) and the Tollgate watershed (TOL) [ , 2000 , 2001a ; , 2001 ]. 2 A representation of topographic and vegetation attributes of the upper part of the Reynolds Creek experimental watershed. (a) The Digital Elevation Model of the Tollgate watershed. (b) The vegetation map of the Tollgate watershed. (c) The Digital Elevation Model of the Reynolds Creek, Mountain East watershed. The symbols in Figures 2a and 2c indicate the locations of meteorological ( clim. ), precipitation ( prec. ), and snow ( snow ) measurement stations. 3.2.1. Reynolds Creek Mountain East Watershed The Reynolds Mountain East (RME) is a small headwater catchment (0.389 [ km 2 ]) in the southwestern corner of RCEW that ranges in elevation from 2024 to 2139 [ m a.s.l. ] ( Figure 2c ). Hourly meteorological variables, snow depth, snow water equivalent, soil moisture, and discharge for a 25‐year period (October 1983 through October 2008) were collected at two locations ( Figure 2c ). The stations represent major landscape units within the RME catchment: one is a sheltered site, within a clearing of an aspen/fir grove near the center of the catchment; the other is an exposed site dominated by mixed sagebrush near the western catchment divide [ , 2000 ; , 2001 ; , 2001 ; , 2002 ; , 2002 ; , 2009 ; , 2010 ; , 2011 ]. Hourly observations of air temperature, humidity, wind, solar radiation, and precipitation from the wind‐exposed ridge site, along with the hourly data of air temperature, wind, and precipitation, from the grove site were used in the simulation. Snow density and depth observations were available from snow pillow and manual measurements; soil moisture measurements at different depths and locations were also available [ , 2000 , 2001b ; , 2000 , 2001; , 2011 ]. At the sheltered site, during an average water year (October through September), shielded precipitation is 964 [ mm ]. At the exposed site, it is 550 [ mm ] [ , 2011 ]. Prevailing winds blow from the southwest direction with the mean winter‐time direction of 230°. Avalanching does not play a role in redistributing snow in this catchment [ , 2002 ]. The soil in Reynolds Creek, Mountain East area is sandyloam and loam, with elevated contents of rocks and organic fraction [ , 2001b ]. Vegetation of the RME watershed ( Figure 8a ) is characterized by patches of Douglas fir ( Psuedotsuga menziesii ) (4.6%), aspen ( Populus tremuloides ) (26.1%), and sagebrush plants (69.3%), such as low sagebrush ( Artemisia arbuscula ), vaseyana sagebrush, and Wyoming big sagebrush ( Artemisia tridentata ) [ , 2001b ; , 2002 ; , 2002 ; , 2010 ; , 2011 ]. [2010] provide estimates of the maximum leaf area index and vegetation heights: aspens have an average height of 9.5 [ m ] and LAI = 1.35; LAI of fir trees is about 2.0; sagebrush are approximately 60 [ cm ] in height with LAI of 0.77. 8 (a) The vegetation map of the Reynolds Creek Mountain East watershed. (b) The volumetric soil water content averaged over the simulation period of October 1983 through October 2008. Spatial variability of meteorological forcing is a dominant feature impacting hydrology within the Reynolds Creek Experimental watershed [ , 2001 ; , 2010 ]. It justifies the need to accurately interpolate meteorological inputs over the area of the RME watershed. Air temperature and relative humidity are linearly interpolated between the two available stations, accounting for the elevation effect and using the observed gradients of these variables at each hour. Precipitation and wind speed are interpolated accounting for terrain structure and vegetation, following the approach outlined by [2002] and [2002] . Specifically, an accumulation factor for precipitation, and a wind terrain factor are used to spatially interpolate these variables using four classes of terrain and vegetation shelter computed with the parameters S x (maximum upwind slope) and S b (slope break). Lastly, incident solar radiation in the watershed is interpolated accounting for the local and remote shading effects of topography [ , 2011 ]. The ecohydrological dynamics at RME were simulated for the period of October 1983 through October 2008, for a total duration of 25 years. The imposed grid resolution was 30×30 m 2 that led to 433 computational elements. A uniform one meter soil column and nearly impermeable bedrock (conductivity of bedrock equal to 0.01 [ mm h −1 ]) were assumed for all computational elements of the grid. The chosen soil column depth was assumed to be a representative average for the RME watershed, where typical soil depths are between 15 and 200 [ cm ] [ , 2001b , 2009 ]. A spatially uniform value was assumed because of the lack of distributed information. The value of bedrock conductivity was derived from in situ substrate geology [ , 2001a ; , 2011 ]. The vertical mesh was composed of 10 soil layers of increasing thickness, dz , from the surface ( dz = 10 [ mm ]) to a 1m depth ( dz = 200 [ mm ]). Soil hydraulic properties were derived from the pedotransfer functions of [2006] using the spatially variable fractions of sand and clay derived from the soil map [ , 2000 , 2001a ]. Vertical homogeneity of soil hydraulic properties was assumed because of lack of information. An horizontal anisotropy coefficient, a R = 140, was used for all of the grid elements (see Section 4). The vegetation map was used to assign a vegetation type to each element, i.e., sagebrush, fir, or aspen [ , 2002 ; , 2010 ]. A single plant functional type and a vegetated fraction of 0.9 were assumed for each cell, considering that in such an environment small patches of bare soil are very likely. Physiological and structural characteristics of low sagebrush, fir, and aspen inferred from literature were used to parameterize vegetation physiological parameters [ , 2000 ; , 2007 ]. 3.2.2. Tollgate Watershed The Tollgate watershed represents the mountainous part of the RCEW and has an area of 54.8 [ km 2 ]. Its elevation ranges from 1400 to 2241 [ m a.s.l. ] ( Figure 2a ). The Tollgate watershed representing the mountainous part of the RCEW and containing the RME basin, is characterized by a vegetation cover similar to that of RME, i.e., firs occupy 11.9% of the area, aspens occupy 8.7%, and sagebrushes 79.4% ( Figure 2b ). For this watershed, two additional meteorological stations and 10 precipitation gauges located within and in the proximity to the catchment ( Figure 2a ) were used in the spatial interpolation of meteorological inputs [ , 2001 ; , 2001 ]. However, the topographic effects on wind speed and snow accumulation were neglected, given the lack of appropriate information to extend the [2002] approach to such large spatial scales. Temperature, relative humidity, and precipitation were interpolated, summing the residual and gradient fields. The gradient field was obtained from a linear regression of meteorological observations with altitude. Consequently, the residuals, i.e., the deviations of the measurements from the gradient field, were spatially interpolated using the Inverse Distance Weight (IDW) method [ , 2005 ]. This resulted in a residual field representing the local deviations from the gradient field [ , 2008 ]. The available dataset includes a period of 8 years: August 1985 through July 1993. The simulation period for the TOL catchment corresponds to the entire period of meteorological observations. The grid resolution is 50×50 m 2 , resulting in 21,924 computational elements. The description of soil properties and the representation of vegetation for the TOL catchment are analogous to that of the previously described RME catchment. 4. Model Parameter Adjustment Given the relatively high computational cost of distributed simulations and the high dimension of the parameter space of Tethys‐Chloris, no formal, traditional, or advanced calibration of the model can be carried out [ , 2002 ; , 2002 ; , 2008 ]. The choice of model parameter values has been made subjectively, based only on available data or literature information (see the companion paper of [2012] ). The model was first used for the plot scale simulations [ , 2012 ]. This permitted to control the overall consistency and reliability of simulations, using various metrics of model performance, such as the coefficient of determination, the root mean square error, the Nash‐Sutcliffe efficiency, etc. The estimated performance metrics together with qualitative considerations of plausibility of results at different temporal scales (e.g., annual, monthly and daily cycles) in comparison to other studies [e.g., , 2004 ; , 2004 ; , 2010 ] were used to adjust the model parameters, while keeping them within a realistic range. The parameters that resulted in a satisfactory performance at the plot scale were subsequently used for spatially distributed simulations. The transition from the plot scale to a distributed domain raises the question of transferability and significance of other parameters that directly affect the spatial dynamics (e.g., the conductivity anisotropy ratio, bedrock leakage). More than one attempt is usually necessary to improve the model results at the watershed scale. In the analyzed cases, the parameters adjusted at the plot scale led to a good performance when applied to spatial domains. The results presented in the following were obtained by running the model two times for the Lucky Hills watershed, where the parameter adjustments only concerned the parameterization of the soil sealing process. Five model simulations were carried out for the Reynolds Creek Mountain East watershed: major modifications concerned the hydraulic conductivity and the anisotropy ratio, as well as physiological parameters of vegetation, specifically, the root depth, the empirical parameter linking photosynthesis to stomatal conductance, a , and the maximum Rubisco capacity, V max L [ , 2012 ]. The results obtained for the Tollgate watershed were not calibrated and obtained with the same parameter values used for the Reynolds Creek Mountain East watershed. Only manual parameter adjustment was used. It implies an expert evaluation of the model performance with a subsequent modification of parameters keeping them within a physically realistic range. The overall objective was to obtain results consistent with observations or independent estimates of the simulated variables. Therefore, when compared to traditional hydrological modeling, the model performance skill relies more on the model structure rather than on the skill of a calibration procedure. Undoubtedly, uncertainties in meteorological inputs and boundary conditions limit a thorough testing of the T&C model structure. Nonetheless, the background rationale is that the structure of a physically‐based/mechanistic model should lead to a satisfactory performance with a narrow range of possible outcomes. This is possible because the physically‐based (i.e., theoretically measurable) nature of the parameters implies narrow bounds of values that they can assume, when used in the context of representing an actual physical media/ecosystem. The direct consequence is a “tighter” range of possible model results. This further implies that a mismatch between model simulations and observations can be attributed to the model structure or uncertainty of boundary conditions, rather than to a poor calibration procedure. Furthermore, the narrow feasible region of model parameters diminishes the importance of subjectivity in choosing parameter values, i.e., the expected outcome of a simulation cannot be too different as a consequence of a different parameter value. Calibration in our study, thus represents a final adjustment to refine the simulation skill, which is mainly dictated by the model structure and boundary conditions. 5. Distributed Results The results obtained with T&C for the Lucky Hills watershed and Reynolds Creek Mountain East/Tollgate watersheds are discussed in the following. 5.1. The Lucky Hills Watershed In a spatially distributed application, hydrometeorological inputs can vary among the computational elements due to topographic or local meteorological conditions. Given the small watershed area of Lucky Hills, the only input that exhibits spatial distribution is the incoming shortwave radiation (Section 3.1). Both local and remote terrain effects are accounted for in the simulation. Figure 3a shows that the inclusion of topographic effects modifies the distribution of shortwave incoming radiation. As expected, the steep slope exposed to the south receives shortwave radiation that over the long‐term is about 30% higher than that for the steepest element exposed to the north. This uneven distribution of radiation is directly reflected in the long‐term distribution of simulated transpiration fluxes ( Figure 3b ) and bare soil evaporation ( Figure 3c ). The relatively smaller available energy in the north‐facing elements leads to a lower photosynthetic activity and reduced transpiration fluxes. Related to that effect is a slight increase in bare soil evaporation (although it is fairly constant across the domain). Water availability also plays a significant role in shaping the spatial pattern of transpiration. For example, transpiration and, to some extent, bare soil evaporation are higher in the convergent part of the catchment, where higher soil water content is a result of lateral water redistribution. 3 The results of hourly spatially‐distributed ecohydrological simulations averaged over the simulation period (July 1996 through December 2009) for the Lucky Hills experimental watershed. (a) Incoming shortwave radiation. (b) Transpiration flux. (c) Bare soil evaporation flux. Water redistribution in such an environment is expected to be controlled by runon‐re‐infiltration effects [ , 2002 ; , 2003 ; , 2003 ], rather than by subsurface flows. This is confirmed by the model results. The computed lateral subsurface redistribution of water between neighboring cells is a small fraction of the annual hydrological budget ( Figure 4a ). It is generally less than 2.0 [ mm yr −1 ] and peaks in steepest slopes. The spatial differences among the mean annual infiltration rates are significant and they are a consequence of localized runon‐re‐infiltration process ( Figure 4b ). Most of water redistribution is produced by short intense events during monsoon seasons that can lead to infiltration excess runoff. The process representation is accommodated by using disaggregated rainfall at 5 [ min ] time step [ , 2005 , 2008 ; , 2009 ; , 2010 ] and accounting for the soil sealing formation [ , 2012 ]. Runoff produced after intense precipitation events is routed through the watershed drainage paths towards the outlet. Although a large fraction of generated runoff is “lost” as streamflow/overland flow, there are topographic locations with shallow slopes that facilitate reinfiltration of runon. This process is most efficient at several locations along the hollow part of the catchment, where the mean annul infiltration rates are higher than the annual precipitation ( Figure 4b ). Overall, runon subsurface flows, and evapotranspiration fluxes all contribute to creating a heterogenous distribution of mean soil water content ( Figure 4c ). The wettest parts are located along the trough, where the process of runon occurs [ , 2003 ] and the north slope, where transpiration fluxes are smaller due to lower surface irradiance. 4 The results of hourly spatially‐distributed ecohydrological simulations averaged over the simulation period (July 1996 through December 2009) for the Lucky Hills experimental watershed. (a) Total lateral subsurface flow. (b) Infiltration flux. (c) Soil moisture content. The mean canopy characteristics during the 13.5‐year simulation period are the result of vegetation adaptation to local topographic and environmental conditions: a larger LAI is simulated in the trough and a lower LAI in steep hillslopes ( Figure 5a ). Note that both the north and south facing slopes have smaller LAI values, as compared to the other parts of the catchment. This is due to the relatively drier conditions in the hillslope exposed to the south and lower levels of shortwave radiation in the north‐facing slope. Because of the fairly gentle topography and limited water redistribution, the relative spatial variability of LAI is small (≈10%). Above ground Net Primary Productivity of vegetation (ANPP) has a spatial distribution similar to that of LAI. However, it has a larger relative variability, with the maximum difference of ≈20% ( Figure 5b ), highlighting a non‐linear relationship between LAI and productivity. Simulated ANPP values of 70–85 [ gC m −2 yr −1 ] are regarded as very plausible for creosotebush in southern Arizona and match values reported in literature very well [ , 1965 ]. As can be observed in Figure 5c , the distribution of surface radiative temperature reflects the forcing of incoming shortwave radiation, with only secondary effects due to vegetation cover and soil moisture distribution (not appreciable from the figure). 5 The results of hourly spatially‐distributed ecohydrological simulations averaged over the simulation period (July 1996 through December 2009) for the Lucky Hills experimental watershed. (a) Leaf Area Index. (b) Above Ground Net Primary Productivity. (c) Surface radiative temperature, T S . The simulated dynamics of LAI and Gross Primary Production (GPP) have been compared with estimates based on remote sensing MODIS data ( Figure 6 ). A single MODIS cell contains the entire Lucky Hills catchment. The T&C model is able to reproduce the LAI seasonality and its overall magnitude quite well. The GPP is slightly overestimated with the seasonal maximum following the monsoon rainfall, especially in the relatively wetter years. The small size of the Lucky Hills watershed as compared to the resolution of MODIS vegetation products (1 km 2 ) precludes any spatial comparison between the simulated and observed patterns. However, it is important to note that the distributed application produces spatially integrated fluxes and states that are nearly identical to those simulated in the plot scale application [ , 2012 ]. The determination coefficients of the relationships between the time series obtained at the plot scale and spatially integrated LAI and GPP are 0.96 and 0.94, respectively. This is due to the small dimension of the watershed, its gentle topography, and the limited impact of surface‐subsurface flows. It is very likely that extending the analysis to a larger watershed area in the order of tens of square kilometers in this climatic zone would still provide very similar results. Consequently, for the examined case study, the plot scale simulations can be representative of a much larger area (with the apparent exception of streamflow magnitude). 6 Spatially integrated (a) Leaf Area Index, (b) Gross Primary Production, (GPP), based on the simulated (SIM.) and remote sensing estimates for the Lucky Hills watershed. The simulated annual partition of water budget among hydrological fluxes is detailed in Table 1 . Evaporation and transpiration represent the largest components of the hydrological budget. The average annual evapotranspiration is ≈310 [ mm yr −1 ], which represents 93% of the annual precipitation. Recharge to deeper soil layers is estimated to be ≈9 [ mm yr −1 ]; it is highly discontinuous during the simulation period and concentrated in a few wet periods. Infiltration excess runoff can be observed during summer months of monsoon periods due to heavy precipitation and effects of soil sealing. The streamflow at the catchment outlet is ≈23 [ mm yr −1 ]. This value is very close to the observed values for small nested basins of Walnut Gulch [ , 2008 ]. 1 The Partition Among the Hydrological Budget Components for the Lucky‐Hills Watershed for a 13.5‐Year Simulation Period (July 1996–December 2009). a Component Flux [ mm yr −1 ] P r 333.8 E G 109.3 T 191.4 E In 9.4 E snow 0.0 L K 9.0 Q 23.2 a The hydrological budget terms are: precipitation, P r , soil evaporation, E G , transpiration, T , evaporation from interception storage, E In , evaporation/sublimation of snow, E snow , deep leakage, L K , and outlet discharge, Q . The imbalance is due to a slight storage change between the beginning and the end of the simulation. The cumulative runoff simulated by the model for the last nine years of the simulation is compared with the observations in Figure 7a . There are differences between the simulated and observed streamflows, especially during the years of 2006–2007. The overall result, however, is surprisingly good. For the period of 2000 through 2009, the simulated annual runoff is 21.1 [ mm yr −1 ], while the observed runoff is 23 [ mm yr −1 ]. The error is essentially negligible, given the minimal calibration of the model, the uncertainties associated with the soil hydraulic properties, bedrock conditions, soil seal formation, and the adopted rainfall disaggregation technique. Similarly, small differences can be also observed in the flow‐duration curve ( Figure 7b ), where larger runoff values are somewhat overestimated. A very good agreement is exhibited for the durations of ≈0.5–0.8 [ day ], above which the two curves show essentially zero runoff. The simulation of runoff that is related to sporadic and intense events (typical of semiarid systems) is a challenging problem for most hydrological models. This is especially true when the objective is to truthfully simulate all of the involved physical processes. Therefore the presented simulation results can be considered to be “effectively” corroborated by the data. 7 A comparison between the observed (“OBS.”) and simulated (“SIM.”) (a) hourly cumulative runoff and (b) flow duration curves for the Lucky Hills watershed. 5.2. Reynolds Creek Mountain East Watershed The map of vegetation of the RME watershed is shown in Figure 8a . As seen, vegetation cover is rather diverse: aspen groves are located in the convergent part of the catchment and near the southeastern boundary. Patches of fir border the upstream side of the aspen grove in the center. The remaining part of the catchment is occupied by sagebrush. 5.2.1. Hydrological Partition The distribution of the top 1‐m soil moisture averaged over the 25‐year simulation period (1983 through 2008) is shown in Figure 8b . Two principal controls can be discerned, related either to soil texture or topography. Specifically, the sandy southwestern part of the catchment has consistently lower soil water contents, as compared to the rest of the domain characterized by a finer soil texture. The convergent part of the watershed has the highest soil moisture, especially in the downslope areas affected by snow drifts, where accumulation of snowpack allows the persistence of wet conditions over longer periods (see later Figure 9b ). 9 The results of spatially‐distributed ecohydrological simulations averaged over the simulation period (October 1983 through October 2008) for the Reynolds Creek, Mountain East watershed. (a) Annual transpiration flux. (b) Annual fraction of time with snow cover. The simulated lateral fluxes are relatively high due to the pronounced slopes and the high anisotropy ratio chosen for the RME watershed to mimic preferential flows associated with topography [ , 2005 ; , 2010 ]. These fluxes permit an efficacious redistribution of water that tends to converge in the central part of the watershed, ultimately leading to saturated areas downstream of the central aspen grove during melting seasons. The redistribution effect is clearly appreciable in the spatial distribution of transpiration fluxes ( Figure 9a ) that are higher in the convergent part of the catchment. The simulated spatial distribution of transpiration is due to a combined effect of moisture availability during growing season and uptake characteristics imposed for a given vegetation type. For example, taller trees, such as aspens and firs with deeper roots have access to deeper soil water and tend to transpire over longer periods with higher rates. A drastic reduction of transpiration can be observed for areas with significant snow accumulation. This abrupt decrease is partly due to the model assumption that neglects transpiration fluxes, when a cell is covered by snow, and because of a delayed start of the growing season in the snow‐covered elements. A patch of relatively small transpiration (around 300–500 m north and 400–600 m east) can be observed in Figure 9b . This is an interesting result because that area has been previously identified by field studies as a drier location [ , 2009 ; , 2011 ]. This has been now consistently confirmed by the model results. The signature of vegetation is also evident in the duration of snow cover on the ground, shown in Figure 9b . The simulated patterns are also strongly influenced by meteorological inputs. For example, snow tends to accumulate in drifts in sheltered areas (i.e., these areas correspond to the values higher than 0.6 in Figure 9b ), while precipitation tends to be less and sublimation rates are higher at wind‐exposed locations (south‐eastern part of the watershed). For the same topographic aspect and slope, aspens reduce the duration of snow cover due to snow interception effects and because of a larger aerodynamic roughness that facilitates snow sublimation. The integrated effect of firs is expressed by both enhanced turbulent exchange and a decreased incoming energy to the ground snowpack due to canopy shadow effects [see also 2006 ]. The spatial distribution of the fraction of time with snow cover agrees well with observations and model simulations previously carried out for the RME watershed by [2002] and [2002] for specific years. The mean annual partition into hydrological components is presented in Table 2 . The average simulated runoff is 543 [ mm yr −1 ], which represents about 59% of annual precipitation and compares favorably to the mean observed streamflow of 518 [ mm yr −1 ]. The latent heat flux is partitioned into bare soil evaporation, 105 [ mm yr −1 ], transpiration, 203 [ mm yr −1 ], and snow sublimation, 59 [ mm yr −1 ]. Given the small imposed bedrock conductivity, the leakage losses are negligible in the overall balance and constitute less than 1 [ mm yr −1 ]. As expected, the transpiration flux is higher than bare soil evaporation and snow sublimation. However, the magnitudes of these three fluxes are comparable in this mountain environment, implying that none can be neglected without introducing significant errors. 2 The Partition Among the Hydrological Budget Components for the TOL Watershed for a 8‐Year Simulation Period (August 1985–July 1993) and for the RME Watershed for a 25‐Year Period (October 1983–October 2008). a Component Flux (TOL) [ mm yr −1 ] Flux (RME) [ mm yr −1 ] P r 653.0 912.4 E G 100.7 105.0 T 232.6 203.8 E In 7.4 0.3 E snow 112.7 59.7 L K 0.2 0.7 Q 202.0 543.5 a The hydrological budget terms are: precipitation, P r , soil evaporation, E G , transpiration, T , evaporation from interception storage, E In , evaporation/sublimation of snow, E snow , deep leakage, L K , and outlet discharge, Q . The imbalance is due to a storage change between the beginning and the end of the simulation. 5.2.2. Streamflow The overall consistent and representative performance of T&C is further corroborated by the simulated discharge at the outlet of the RME watershed, as shown in Figure 10 . The cumulative discharge is reproduced with a sufficiently high accuracy over the entire period ( Figure 10a ). A marginal deviation occurs in the last 7 years. This is also testified by an agreement between the simulated and observed annual runoff for each water‐year (October through September) of the simulation ( Figure 10b ). The global Nash‐Sutcliffe efficiency, (NS), for the hourly discharge values is 0.59, and the coefficient of determination between the observed and simulated series is R 2 = 0.65. The result is somewhat poorer in terms of the flow‐duration curve ( Figure 10c ). The model is able to reproduce the higher discharge rates with the proper duration but tends to underestimate the flow minima. For instance, for durations larger than 100 days, the simulated discharge is typically lower than the observed one. In other words, the simulated flow rates exhibit dampening of the recession curve, while observations suggest persisting flows. 10 A comparison between the observed and simulated (a) cumulative runoff, (b) water‐year values of runoff and precipitation, and (c) flow duration curves for the Reynolds Creek, Mountain East watershed. Q obs and Q sim denote the observed and simulated annual discharge, respectively; P r denotes the annual precipitation averaged over the watershed. The simulated discharge series are illustrated in Figure 11 for each of the 25 water‐years. In general, the discharge is simulated in a realistic fashion in most water‐years but the overall quality of simulations changes from year to year. This is not surprising, given the possible uncertainty of the input data such as snow, or the capability of the model to better reproduce hydrological conditions of a given year. In quantitative terms, the model skill for the water‐year 1996 is extremely good, ( NS = 0.83), while for the water‐year 2004 is very poor ( NS = −0.25). There is no a clear time‐evolution of the model performance and a negative efficiency is obtained only in 3 out of 25 years; for most years, NS >0.5. It can be concluded that the simulation skill reflects the physically‐based nature of represented processes responding in a consistent fashion to hydrometeorological conditions of a given year. 11 A comparison between the hourly observed (“OBS.”) and simulated (“SIM.”) discharge at the RME watershed outlet for each of the 25 water‐years (October to September). 5.3. Tollgate Watershed Spatially distributed results for the TOL watershed were averaged over the 8‐year simulation period of August 1985 through July 1993, that corresponds to the period of data availability for the entire watershed (Section 3.2.2). Over this period, MODIS data were not yet available. Remote sensing estimates of snow cover, LAI, and GPP (MODIS products: MOD10A2, MOD15A2, and MOD17A2), for the period of January 2001 through December 2010 were used to confirm the model performance. Given the two different periods (i.e., the period of watershed data availability and the period of remote sensing observations), results can only be compared in terms of average spatial patterns and seasonal dynamics under the assumptions of climate and vegetation stationarity. While the assumption of vegetation stationarity can be considered rather acceptable for a fifteen year temporal lag (1985–2011), certain climatic trends have been observed for the Reynolds Creek watershed [ , 2010 ]. Specifically, when climate metrics for the periods of October 1985 through October 1993 and October 2000 through October 2008 are compared for the RME stations (the only stations for which data were available for the period of 1983–2008, see Section 3.2.2), an increase of air temperature by 0.16 [° C ] and precipitation by 122 [ mm yr −1 ] is observed over the 2000–2008 period. Therefore, in a comparison of T&C results and MODIS observations, we need to qualitatively account for the fact that the period of MODIS operation corresponds to a somewhat warmer and wetter climate (see also Figure 10 ). In order to present a “mechanistic” insight on an internal structure of the watershed response, model performance is also illustrated in terms of spatial patterns for bare soil evaporation, transpiration, snow sublimation, and effective saturation. While these comparisons are only qualitative, they further corroborate the effective consistency of T&C in reproducing realistic dynamics. They further argue for the need of detailed spatially‐distributed observations of states and fluxes of a land‐surface for a rigorous model confirmation. 5.3.1. Hydrological Partition Figure 12a shows the simulated distribution of bare soil evaporation across the TOL watershed. Bare soil evaporation typically ranges between 60 and 150 [ mm yr −1 ], with higher values in the convergent topographic areas or the areas exposed to stronger winds; these two corresponding controls are the emerging features of the simulated pattern. In the northwestern part of the catchment, bare soil evaporation tends to be lower. This is the result of enhanced activity of vegetation ( Figure 12b ) that decreases the amount of water available for soil evaporation. The distribution of transpiration flux in Figure 12b shows that larger rates are typical for areas covered with fir but also that moisture distribution exerts an important control on this flux. Specifically, higher transpiration rates in the western and northern parts of the catchment are the result of a larger amount of precipitation [ , 2000 ]. Being wetter, the western part of the catchment favors vegetation covers and an adaptation of taller vegetation types ( Figure 2b ). The control of soil moisture is clearly appreciable in the downstream part of the watershed: the simulated transpiration fluxes are different in the western and eastern sides of the catchment, despite a uniform sagebrush fractional cover. Furthermore, a north‐south elevation gradient of transpiration can be appreciated for areas vegetated with sagebrush. Low elevations that have shorter snow cover periods and higher temperatures favor photosynthetic activity and, consequently, transpiration. 12 The results of spatially‐distributed ecohydrological simulations averaged over the simulation period (August 1985 through July 1993) for the Tollgate watershed. (a) Mean annual bare soil evaporation flux. (b) Mean annual transpiration flux. Figure 13a illustrates the simulated average fraction of time that each computational element is covered by snow. As already pointed out for the RME watershed, the signature of vegetation is the predominant feature of the simulated spatial pattern. The canopy shadow effect exerted by vegetation on snow is not as significant as the impact of interception and higher roughness resulting in larger snow evaporation/sublimation in aspen and fir areas (not explicitly presented). The total snowpack that accumulates on the ground below aspens and firs is smaller, when compared to open areas with sagebrush. Although melting of undercanopy snow is slower due to the less available energy, the duration of snow season is shorter because of the significantly reduced snow water equivalent. When compared with MODIS observations the simulated patterns of temporal fraction of snow cover are very consistent ( Figures 13a and 13b ). The major appreciable differences are at low watershed elevations, where T&C simulates a longer snow cover, and in the fir‐aspen vegetated areas at high elevations. The differences at the lower elevations can be the result of aforementioned warmer climate conditions corresponding to the MODIS observation period. The lack of a vegetation pattern in the data might be attributed to the MODIS resolution (about 500 m×500 m), which is clearly insufficient to capture the smaller‐scale variability due to topographic features and vegetation effects. Despite these differences the fraction of time with snow cover averaged over the watershed is 0.35 in T&C simulations, and 0.37 in MODIS observations, a surprisingly good match. A similarly successful comparison is also obtained in reproducing the time dynamics of snow cover fraction: Figure 13c shows the space‐time average seasonal dynamics over the Tollgate watershed. 13 The fraction of time with snow cover for the Tollgate watershed. (a) The simulated pattern over the period of August 1985 through July 1993; (b) remote sensing estimations over the period of January 2001 through December 2010; (c) average seasonal cycles from remote sensing estimates and simulations. The signature of vegetation distribution is also evident in the patterns of snow evaporation/sublimation ( Figure 14a ). Taller plants significantly increase the rates of snow sublimation because of snow interception by canopies and turbulence enhancement (due to the larger roughness of vegetated surfaces). This is expressed especially well for fir‐vegetated areas. Therefore, vegetation has a very important role in controlling the amount and dynamics of snowpack, exerting both positive (radiation interception) and negative (snow interception, turbulence enhancement) feedbacks. The net outcome, however, is dependent on both vegetation structure and local meteorological conditions, especially wind speed. 14 The results of spatially‐distributed ecohydrological simulations averaged over the simulation period (August 1985 through July 1993). (a) Mean annual snow evaporation/sublimation. (b) Effective saturation. The two or three orders of magnitude of difference of the characteristic roughness of snow‐covered sagebrush vs. aspen/fir grove areas lead to drastic differences in snow evaporation rates. They are typically 2–4 times larger for aspen/fir areas: 70–100 [ mm yr −1 ] for sagebrush areas and 200–300 [ mm yr −1 ] for aspen/firs. In conditions of lower wind speeds (for instance in the northern valleys of the basin), the difference in sublimation between the two land cover types is smaller, e.g., 40 vs. 60 [ mm yr −1 ]. Obviously, these results are affected by large uncertainties, such as the representation of the within‐plant turbulence effects or the spatial distribution of wind speed fields because of a limited number of meteorological stations. Nonetheless, the simulated snow evaporation/sublimation rates are comparable to those obtained in other modeling studies for mountainous areas [ , 2008 ; , 2010 ]. The spatial distribution of effective saturation clearly highlights that the northeastern part of the watershed is relatively drier than the rest of the basin. It can also be inferred that vegetation, especially firs, create patches of lower water content ( Figure 14b ). The presence of firs near the channel network leads to a somewhat smaller soil moisture in the elements containing channels. The areas of lowest effective saturation are due to a combination of soil texture and vegetation related effects. Specifically, these areas are typical for coarse sandy soils with higher hydraulic conductivities and covered by firs. The coarse soil texture facilitates moisture uptake by deep roots of tall vegetation. Simulations results indicate that subsurface lateral flows are the major source of water redistribution in this watershed and they typically peak in steepest slopes that are in the southern part of the watershed. The annual partition among the simulated hydrological components is presented in Table 2 . The average simulated runoff is 202 [ mm yr −1 ], this represents ≈31% of annual precipitation and compares well with the observed discharge that is 188 [ mm yr −1 ]. The magnitude of bare soil evaporation is 100 [ mm yr −1 ], and of transpiration is 232 [ mm yr −1 ]. The rate of snow evaporation is quite high, 112 [ mm yr −1 ], due to a high average interpolated wind speed in the TOL watershed, i.e., 3.47 [ m s −1 ] (this is relative to the RME watershed, where the average value is 2.94 [ m s −1 ]), and a larger areal fraction covered by evergreens (Section 3.2.2). 5.3.2. Vegetation Metrics The spatial distribution of vegetation Leaf Area Index (LAI) is presented in Figure 15 . The simulated mean LAI reflects adaptation of vegetation to long‐term environmental conditions, i.e., although at the scale of each computational element the vegetative fraction is assumed to be constant (equal to 0.9), soil moisture availability, radiation, air temperature, etc., all contribute to the canopy dynamics and therefore to the spatial distribution of average LAI ( Figure 15a ). A comparison between the simulated and the MODIS derived LAI ( Figure 15b ) is significantly hampered by the MODIS resolution for LAI (1 km×1 km). Nonetheless, qualitatively similar spatial distribution and magnitude of LAI values can be observed in Figure 15 . Specifically, lower LAI in the northeastern part and higher LAI in the central‐southwest part of the watershed (≈1–2 kilometers from the divide) can be inferred. The catchment simulated average LAI is 0.61 and the MODIS‐derived LAI is 0.66. 15 The Leaf Area Index for the Tollgate watershed. (a) The mean annual simulated pattern over the period of August 1985 through July 1993; (b) the mean annual LAI based on remote sensing estimations over the period of January 2001 through December 2010; (c) average seasonal cycles from remote sensing estimates and simulations. The highest LAI values are simulated for evergreen trees since they do not shed leaves. Firs have values of LAI ≈2.5 in the wettest, southeastern part of the catchment, and ≈1.5 in the drier northern part. Conversely, aspens exhibit a pattern of LAI that is inversely related to the elevation, e.g., the multi‐annual average values of LAI are around 0.5–0.6 at higher elevations of the southern part and 0.8–0.9 at some locations in the northern part of the catchment. Since aspens are deciduous trees, their summer LAI is typically 2–2.5 higher than the annual average. Lastly, the LAI of sagebrush is influenced by many factors, such as elevation (air temperature), moisture availability, and soil texture. The sandy northwestern part of the catchment exhibits LAI magnitudes around 0.5–0.9, conversely, the northeastern part and some of the locations closest to the mountainous divides or in topographically shaded areas exhibit the lowest values of 0.25–0.4. Similarly to aspens, the summer peaks of sagebrush LAI are typically 1.5–2 higher than the annual average. The simulated LAI peaks compares surprisingly well with field observations of maximum LAI carried out by [2010] in the RME watershed for the three represented vegetation types (see Section 3.2.2). The average seasonal dynamics of LAI over the Tollgate watershed are shown in Figure 15c , based on simulations and MODIS data. The simulated LAI dynamics seemingly overestimate winter LAI, which is simply an artifact due to the presence of evergreen firs and shrubs. This “discrepancy”, however, can be due to poor reliability of MODIS estimates, during periods when vegetation is totally or partially buried by snow (the case for the Tollgate watershed). The difference during the summer peak can be related to the generally wetter and warmer conditions of the MODIS observational period. Further data would be necessary to corroborate such a statement. Note also that the differences are of the order of 0.1–0.3 LAI, which is below typical uncertainties associated with LAI estimation, even when it is made from the ground [ , 2003 ]. The spatial pattern of the long‐term average Gross Primary Production (GPP) is shown in Figure 16 . As seen, features present in the distribution of LAI can be also appreciated in the distribution of GPP. The spatial organization of GPP is even more pronounced than that of LAI. Evergreen firs have GPP mainly in the range of ≈350–650 [ gC m −2 yr −1 ], with a spatially average value of 485 [ gC m −2 yr −1 ]; deciduous aspens have a somewhat lower productivity with GPP of ≈250–450 [ gC m −2 yr −1 ], and a spatially average of 329 [ gC m −2 yr −1 ]. Sagebrush plants exhibit a wider range of environmental conditions with the corresponding variability of GPP in the range of 250–550 [ gC m −2 yr −1 ] and a spatially average value of 236 [ gC m −2 yr −1 ]. Remote sensing estimates mainly confirm these spatial results, except for the northwestern part of the watershed where a higher productivity is simulated ( Figure 16b ). The available resolution of MODIS product (1 km×1 km) prevents further spatial confirmation of the modeling results. However, MODIS inferred average annual GPP (for the period 2001–2010) is 270 [ gC m −2 yr −1 ]; the corresponding simulated value averaged over the TOL watershed is 274 [ gC m −2 yr −1 ]. The average seasonal dynamics of GPP are also well captured as demonstrated in Figure 16c . The two average yearly cycles are quite similar during most of the season with the difference during late summer and early fall, where MODIS estimates are double that of the simulated results. This difference can probably be attributed to the warmer and wetter climate over the period of MODIS deployment. A possible outcome of such conditions is a delayed onset of summer drought and lengthening of the growing season. Note how soil moisture limitation effects are pronounced in GPP during the second part of the year both in remote sensing estimates and simulations. MODIS inferences are not sufficiently detailed to confirm the productivity for the three individual PFTs. However, the simulated values are in the expected range for such an ecosystem [ , 2003 ; , 2005 ; , 2011 ]. 16 The Gross Primary Production for the Tollgate watershed. (a) The mean annual simulated pattern over the period of August 1985 through July 1993; (b) the mean annual GPP based on remote sensing estimations over the period of January 2001 through December 2010; (c) average seasonal cycles from remote sensing estimates and simulations. 5.3.3. Streamflow and Snow Depth The time series of streamflow at the outlet of the TOL watershed are shown in Figure 17 , where the runoff rates ( Figure 17a ) and the cumulative discharge ( Figure 17b ) are illustrated. While the total runoff rate is predicted satisfactorily, the model tends to overestimate the discharge in the early part of the melting season and, conversely, to underestimate it later in the summer ( Figure 17b ). Two main reasons can create such a mismatch: (i) the absence of wind induced snow re‐distribution in the model simulation for the TOL watershed; and (ii) the uniform soil depth of 1 m assumed for the entire domain. Snow redistribution can delay snow melt, creating areas with large amounts of accumulated snow. A variable soil depth can produce deep storages of water at the beginning of the melting season that can feed the stream later during the summer. Such processes are not accounted for in the simulation. However, the overall coefficient of determination computed using hourly data of the simulated period is R 2 = 0.59, which can be considered as an acceptable result, given that no effort was applied to calibrate the model. 17 A comparison between the observed and simulated series for the Tollgate watershed (a) hourly outlet discharge and (b) cumulative runoff. Finally, Figure 18 presents an explicit confirmation of the distributed simulation using observed data on snow water equivalent collected at different sites of the catchment ( Figure 2a ). The comparison was carried out for four locations. Note that confirmation is entirely “blind”, implying no effort was applied to tune the simulation results. The model is able to reproduce satisfactorily the distributed dynamics of snow water equivalent, although it appears to underestimate the ground snowpack in years characterized by smaller snow precipitation. The model is spatially consistent across different precipitation conditions (note the differences in scale of the ordinate axes in Figure 18 ). The observed underestimation of T&C can be related to the model structure but it can be also affected by the spatial interpolation of precipitation and wind speed, or by the fact that snow water equivalent measurements are typically collected at sheltered sites, where snow drift accumulation might occur. 18 A comparison between the observed (“OBS.”) and simulated (“SIM.”) snow water equivalent at four monitoring locations (UTM zone 11, coordinates: first location: x = 515864, y = 4771970, second location x = 514041, y = 4769438, third location x = 521613, y = 4769718, fourth location x = 520055, y = 4768117) within the Tollgate watershed. 6. Discussion and Conclusions The modeling approach presented here relies on essential principles of land‐surface physics and life‐cycle vegetation processes. It allows one to reasonably approximate the entire range of system dynamics without summarizing information in conceptual, spatially lumped variables or using assumptions that distort the first principles of physics. Despite certain simplifications, the model provides results that can be confirmed with a large spectrum of metrics that are theoretically observable at different scales [see also , 2012 ]. The result of simulations conducted so far suggest that the model satisfactorily reproduces principal mechanisms that control the system's response. Since data required for a detailed validation/verification of a model such T&C are typically unavailable at the landscape scale of application, a community effort is warranted. This study further advocates [ , 2004 ; , 2006 ; , 2008a ; , 2011 ] the need to develop methods for collecting data on multiple variables at different spatial and temporal scales. This will facilitate the characterization of spatial‐temporal processes of hydrology and vegetation dynamics and provide better constraints of model performance. This should lead to a more robust confirmation of a new generation of Earth‐system models [ , 2006 ; , 2010 ]. At present, the scarcity of interdisciplinary data makes it difficult or sometimes even impossible to test all of the desired behavioral aspects of the model. For example, see the discussion on snow dynamics below vegetation canopies presented in this study. In this regard, the model can be also used as an exploratory tool to design interdisciplinary field campaigns. The capability of Tethys‐Chloris to produce consistent results in terms of many hydrological and ecological metrics has been demonstrated to a degree allowed by available data for two very different environments at the watershed scale. The performance obtained for the Lucky Hills site, a desert shrub ecosystem, and for the Reynolds Creek watershed exhibiting a complex ecohydrological system with different soil textures, vegetation types, and climatic gradients is considered to be satisfactory. Although, we cannot solely rely on discharge to confirm a complex spatially distributed ecohydrological model, simulation efficiencies obtained in this study with no, or limited calibration were comparable to those of calibrated conceptual models. This can be regarded as a significant result. Further, the satisfactory performance for the RME watershed for a 25‐year simulation with static parameters and a lack of any calibration is considered to be a noteworthy result. Generally, all of the simulations are obtained without significant (or, as presently popular, “automated”) calibration efforts, despite the large phase‐space dimensionality of the model. The satisfactory performance is thus a consequence of a significant investment into a physically‐realistic structure of the model, and information content used by the model for the analyzed case studies. The latter includes both a priori information on the feasible range of model parameters and detailed meteorological and catchment hydrological data. The aforementioned shortcomings, such as the limited skill in reproducing the recession curve in the Reynolds Creek subwatersheds, is likely due to the assumptions regarding the soil depth and leakage into the fractured bedrock aquifer. Presently, these are nearly impossible to verify. The shallow soil mantle used in the simulations cannot capture water table dynamics at deeper locations and their effect on streamflow [ , 2011 ]. Nevertheless, we should note that very small flow rates, around 1–10 [ l s −1 ], are typically difficult to capture with any deterministic model. A synthesis with stochastic/deterministic models of flow in fractured media [ , 2002 ] is therefore warranted. Undoubtedly, the model performance for the larger watershed (TOL), especially in terms of discharge, can be ameliorated to better match the observations. The small amount of information available in terms of inputs and boundary conditions is likely to be one of the reasons of poorer model results in terms of streamflow. Nonetheless, we speculate that while the revealed patterns do not necessarily reproduce the actual ones with a high precision, they are highly plausible, given the imposed boundary conditions and the available data. For example, the corroboration between the basin‐averaged LAI, GPP, fractional snow cover, and inferences based on MODIS data for the Tollgate watershed are a surprisingly good confirmation of the simulated vegetation and snow dynamics. Note that LAI dynamics emerge only from the imposed vegetation physiological parameters and environmental conditions. Furthermore, some of the discussed findings represent interesting insights into the quantitative analysis of the dominant hydrological processes at the landscape/watershed scale. For instance, the role of lateral subsurface flow has been found to be unimportant for water redistribution, as compared to the surface runoff‐runon mechanism, in the semiarid Lucky‐Hills watershed. However, it appears to be a very important mechanism for redistributing water in the Reynolds Creek watershed. A peculiar insight is that vegetation imposes an important signature on the snowpack dynamics by modifying precipitation, radiative energy fluxes, and turbulent exchange. Micrometeorological conditions and plant structure can have a significant role in the overall water budget at the watershed scale. This implies that spatial vegetation processes, such as tree encroachment in mountainous areas [e.g., , 2007 , 2008 ], could significantly affect water fluxes into the soil and, consequently, streamflow and deep recharge, in the long‐term influencing aquifers and water budget of these areas. Acknowledgments Simone Fatichi wishes to thank the support of the International joint Ph.D. program on “Mitigation of Risk due to Natural Hazards on Structures and Infrastructures” between the University of Firenze (Italy) and the T.U. Braunschweig (Germany). Authors are extremely grateful to all the persons involved in the collection of data and information in the two analyzed locations, Lucky‐Hills (AZ) and Reynolds Creek (ID). The datasets for Lucky Hills were provided by the USDA‐ARS Southwest Watershed Research Center. Funding for these datasets was provided by the United States Department of Agriculture, Agricultural Research Service. The datasets for Reynolds Creek were provided by USDA‐ARS Northwest Watershed Research Center. The technical support of the Center for Advanced Computing at the University of Michigan is also acknowledged. Valeriy Yu. Ivanov was partially supported through the NSF grant 0911444.

Journal

Journal of Advances in Modeling Earth SystemsWiley

Published: Feb 1, 2012

There are no references for this article.