Transient simulations of an arid mountain system were conducted using a fully integrated model of subsurface and overland water flow and the land surface energy balance (ParFlow). An hourly atmospheric time series was constructed and used to force a two-dimensional model containing detailed, stochastic descriptions of subsurface heterogeneity for alluvial and fractured units. Two cases were simulated for multiple years of simulation: a dry case based on a year with below-average precipitation, and a wet case with above-average precipitation. For each case, four realizations of hydraulic conductivity were simulated, along with a homogeneous domain and a domain with lower fracture density. Detailed analysis of water balances from each simulation indicated that recharge (change in subsurface water storage) over multiple years was significant even for the dry case. The fractured model of subsurface heterogeneity, along with episodic infiltration, combined to create localized regions of fully saturated (perched) water. This was shown to be even more substantial in the wet case. By contrast, homogeneous simulations did not exhibit this behavior and estimated lower recharge than the heterogeneous simulations. Geospatial statistics were subsequently used to evaluate correlations in land-energy fluxes. The land-energy fluxes exhibited clear spatial correlation and were influenced by both the shallow and deeper soil structures. This indicates that land–atmosphere fluxes may provide an integrated measure of subsurface heterogeneity. Finally, the vertical spatial structure of soil saturation in the fractured system was shown to exhibit multiscale behavior.

Detailed geostatistical models of alluvial and fracture permeability and an integrated hydrologic model were used to demonstrate focused recharge and spatial correlation of land-energy fluxes in an arid mountain system

The terrestrial water balance is a complicated set of interrelated processes. Water and energy fluxes from the land surface, notably evapotranspiration and latent heat flux, are the primary mechanisms of communication between the land and atmosphere and have been shown to depend strongly on soil saturation (e.g., Chen and Avissar, 1994; Famiglietti and Wood, 1994; Baldocchi and Xu, 2007; Maxwell et al., 2007; Chen et al., 2008; Kollet and Maxwell, 2008; Kollet, 2009; Kollet et al., 2009; Siqueira et al., 2009). In arid environments, understanding land–atmosphere fluxes is complicated by the episodic nature of precipitation, extremely dry soil moisture conditions, very small potential recharge (i.e., the difference between precipitation and evapotranspiration), vegetative controls, and spatial complexity and scaling (e.g., Scanlon et al., 2002, 2003, 2005; Walvoord et al., 2002; Newman et al., 2006). Subsurface heterogeneity, particularly in fractured systems, further complicates the distribution and movement of moisture in the subsurface, which subsequently affects the spatial and temporal distribution of land–atmosphere fluxes.

Individual components of these arid, fractured systems have been studied previously. For example, Duke et al. (2007) used field-scale tracer tests to study flow and transport in a fractured system consisting of a surficial alluvium underlain by fractured basalt. They found flow in the fractured basalt to be focused in the fractures, with very low increased matrix saturation. Zhou et al. (2006) investigated water seepage through a fractured system at the mesoscale by artificially wetting the soil surface and observing seepage fluxes in a tunnel below. They also found significant variability in fluxes due to fracture patterns. Similar results were obtained in a study of fractured tuff (Salve, 2005) and in a shallow bedrock aquifer in a humid climate (Gleeson and Novakowski, 2009).

Other studies have attempted to upscale infiltration in fractured and arid systems. Heilweil et al. (2007) used a simple regression model calibrated to environmental tracer observations to study net infiltration rates into permeable bedrock under arid conditions. They identified many localized zones of focused infiltration, then attempted to upscale their results to larger areas based on their regression model. Although very different in approach but similar in goal, Ye et al. (2008) used expert elicitation to derive recharge probabilities across a large, arid region in Nevada that encompasses both Yucca Mountain and the Nevada Test Site for five simple recharge models. Faybishenko (2007) used a simplified water balance model (modified Budyko with Penman potential evaporation) to provide bounds for infiltration estimates for Yucca Mountain under current and future climate conditions. In addition to providing a range of recharge estimates, a key conclusion of his work was that uncertainties due to runoff, run-on, heterogeneity, and plant processes need to be included in future studies. Struthers et al. (2007a,b) developed a simple model of fractured infiltration using a two-domain approach. They applied this model to single and multiple storms. Finally, net bedrock infiltration for a watershed on Yucca Mountain was estimated using a watershed model (Woolhiser et al., 2006). The saturated hydraulic conductivity of the subsurface was found to be a highly sensitive parameter.

Recent modeling studies have demonstrated spatial and temporal variability in land–atmosphere water and energy fluxes for different levels of model integration (fully or loosely coupled), scale (local to regional), and subsurface complexity (heterogeneous or homogeneous). These studies include fine-scale simulations of uncorrelated heterogeneous systems with shallow water tables (Zhu and Mohanty, 2003; Mohanty and Zhu, 2007; Kollet et al., 2009), as well as fully coupled homogeneous simulations (Maxwell et al., 2007; Kollet and Maxwell, 2008; Maxwell and Kollet, 2008; Kollet et al., 2009). Those studies that included explicit descriptions of heterogeneity, including both simplified, steady-state simulations (Zhu and Mohanty, 2003; Mohanty and Zhu, 2007) and integrated, fully coupled simulations (Kollet, 2009), demonstrated interdependence between heterogeneity and shallow soil moisture and thus evapotranspiration. Studies that included effective subsurface descriptions also demonstrated correlations between land–atmosphere fluxes and both near-surface saturation and saturated groundwater storage, with topographic influences contributing to the majority of the spatial heterogeneity (Maxwell et al., 2007; Kollet and Maxwell, 2008; Maxwell and Kollet, 2008; Kollet et al., 2009). Of these studies, one in particular used geostatistical techniques to demonstrate spatial autocorrelation in land energy fluxes at the hillslope scale (∼7 km) driven primarily by the topography in the system (Kollet and Maxwell, 2008).

The work presented here both integrates and builds on the approaches used previously to directly investigate water and energy fluxes at fine spatial scales in an arid system characterized by a deep vadose zone and complicated subsurface heterogeneity. Multiyear, transient simulations were performed with a fully integrated model of variably saturated flow, overland flow, and land surface processes (including snow and vegetation). Explicit, structured models of heterogeneity—correlated, Gaussian random fields for a surficial alluvial layer and a stochastic continuum fracture model for a deeper tuff unit—provided a realistic, complex distribution of subsurface properties. Two central questions were addressed by this study:

  1. How does complicated subsurface heterogeneity influence infiltration rates and subsurface moisture distributions for arid systems?

  2. What are the influences of correlated subsurface heterogeneity on spatial patterns of soil saturation and land–atmosphere water and energy fluxes?

A number of techniques were used to understand these relationships, including integrated water balances and spatial and nonspatial statistical measures.

The complete water balance for a hillslope element can be developed from a control volume that encompasses the land cover (vegetation) and a deep vadose zone (Fig. 1 ). Vertical fluxes into and out of the control volume include precipitation (Qprecip) and evapotranspiration (ET) across the upper surface and deep percolation or recharge (Qrechage) across the lower surface; note that deep percolation is conceptualized at a depth much deeper than the root zone. Lateral fluxes include surface fluxes (Qrun-on and Qrunoff) and subsurface fluxes (QS-in and QS-out). Fluxes are balanced by changes in surface storage (Ssurface), composed of ponded surface water and snow, and subsurface storage (Ssubsurface).

The complete conceptual water balance is then given by (Fig. 1a)
This formulation is similar to that proposed in, e.g., Scanlon et al. (2002) but with fewer explicit subdivisions of snowmelt and other surface fluxes (which would be lumped into Qrecharge under the current formulation) but with more explicit subdivision of subsurface lateral flows. At some point, if an element is large enough that it crosses a sufficient number of topographic and watershed divides, lateral subsurface flow may be neglected and run-on along the boundaries may also be neglected and Eq. [1] simplifies to (Fig. 1b)

In this study, I used the flow code ParFlow to simulate the coupled surface and subsurface water balance for a hillslope element. The ParFlow code is a variably saturated groundwater flow code with overland flow represented by a fully coupled overland flow boundary condition and land surface processes represented by an integrated land surface model. A brief summary of the flow equations is presented below; complete details were given by Kollet and Maxwell (2006).

The ParFlow code solves the Richards equation in three spatial dimensions, which may be written as
where ψp is the subsurface pressure head [L], z is the depth below the surface [L], Ks(x) is the saturated hydraulic conductivity [L T−1], kr is the relative permeability (dimensionless) (a function of pressure head ψp), Ss is the specific storage coefficient [L−1], ϕ is the porosity (dimensionless), Sw is the degree of saturation (dimensionless), qs is a general source–sink term [T−1], and t is time [T]. The van Genuchten (1980) relationships were used to relate saturation and relative permeability to the pressure head.
Shallow overland flow is represented in ParFlow by the two-dimensional kinematic wave equation, which appears in the overland flow boundary condition after applying continuity conditions of pressure and flux:
where v is the depth-averaged velocity vector [L T−1], hs is the surface ponding depth [L], and qr(x) is a general source–sink (e.g., rainfall) rate [L T−1]. Note that ‖hs,0‖ indicates the greater of the two quantities, and that the overland flow condition assumes that hs = ψp at the ground surface under saturated conditions (Kollet and Maxwell, 2006). If diffusion terms are neglected, the momentum equation can be written as
which is commonly referred to as the kinematic wave approximation. In Eq. [5], So,i is the bed slope (gravity forcing term) (dimensionless), which is equal to the friction slope Sf,i [L], and i denotes for the x and y direction.
Manning's equation is used to establish a flow depth–discharge relationship:
where n [T L−1/3] is Manning's coefficient.

ParFlow is also coupled to the land surface model CLM (Dai et al., 2003; Maxwell and Miller 2005; Kollet and Maxwell 2008), which solves the full energy budget, including snow and plant processes, at the land surface. These models are coupled such that CLM depends on the surface and subsurface water availability provide by ParFlow. The CLM model, in turn, provides ParFlow with precipitation and evapotranspiration fluxes (via qs). The CLM model is forced by an atmospheric time series of precipitation, solar radiation, wind, temperature, humidity, and pressure. The coupled model, PF.CLM directly simulates the terrestrial components of the hydrologic cycle under transient conditions, including the diurnal cycle, seasonal trends, and a rapid response to meteorological conditions. Complete details of the equations and coupling between ParFlow and CLM have been discussed previously (Kollet et al., 2009) and are not repeated here for brevity; however, some details of the vegetation processes in CLM that are pertinent to the current study are provided here. The CLM model calculates a complete canopy water and energy balance and represents plant transpiration using a photosynthesis-based approach. Stomatal resistance depends on soil suction (in addition to other factors), which is characterized by a linear increase in resistance up to ψp = −150 m, at which point transpiration is assumed to be completely moisture limited and set to zero. The leaf area index is dynamically calculated between minimum and maximum values by plant type and is parameterized as a function of deep soil temperature, creating a seasonal variation in vegetation cover. The root zone is assumed to cover the top 10 model cells below the land surface and a root zone distribution is used to withdraw transpiration fluxes nonuniformly throughout that region. Bare soil evaporation is withdrawn from the model cell at the land surface.

Precipitation, temperature, pressure, wind, and relative humidity data used to force ParFlow in this study were obtained from the Meteorological Data Acquisition (MEDA) Station 12, located on Rainier Mesa at the Nevada Test Site courtesy of NOAA. The data span the years 1983 to 2007 at 15-min intervals. Observed downward solar radiation (short- and longwave) was obtained from the NOAA desert rock SurfRad site located in Mercury, NV. Radiation data were available from 1989 to 2008 at 3-min intervals. Before use with ParFlow, hourly climatologies (values for each hour of the year computed from the observed 1989–2008 record) were calculated for each variable for the available record. Data were processed for missing values and outliers. These data were then aggregated from 15 and 3 min to hourly series. Hours for which no data were available were filled by interpolating values from the previous and following hours. If values from the previous and following hours were also unavailable (i.e., there were large gaps in the observations), then the data was interpolated from the same time the previous and following days to preserve the diurnal cycle; if data gaps were longer than 1 d, then the data were interpolated from the calculated climatological values discussed above.

During this data processing procedure, the precipitation data were aggregated monthly and compared with climatological monthly values from a nearby precipitation station (<100 m, essentially co-located) as reported in Soulé (2006). During this check, it was discovered that precipitation values measured in the winter months were much lower for MEDA Station 12 than that reported in Soulé (2006) due to freezing of the tipping-bucket precipitation gauge in the MEDA station. A companion station uses a heated weighing design and provides more accurate measurements of precipitation that falls in the form of snow. As the more accurate precipitation was not available at finer than monthly temporal resolution, it was not appropriate for the hourly forcing time series needed for this study. The year where monthly precipitation from Station 12 and the nearby station most closely matched, 2006, was used as forcing. With annual precipitation of 253 mm, this year was drier than the climatological average for the site, approximately 326 mm, and will thus be referred to as the “DRY” simulation forcing. Interannual variability in precipitation is significant at this site and more than 25 yr of observations ranges from a minimum of 85 mm to a maximum of 682 mm. To account for the variation between wet and dry years, this time series was perturbed, doubling the precipitation rate while leaving timing unchanged, to 506 mm. This case will be hereafter referred to as the “WET” simulation. While the timing of the wetter years on record does not vary significantly on a monthly scale, it would probably vary on the hourly scale used to force the coupled model. This is a limitation of doubling the precipitation intensity and leaving the timing unchanged. Because much of the precipitation at this site falls in the form of snow, not rain, however, the water flux at the land surface depends on snowmelt more than snowfall. The coupled model used in these simulations (PF.CLM) includes a multilayer snow model that has been shown to accurately represent observed snow water equivalent (Maxwell and Miller, 2005), making the approach more applicable to this system. Time series of precipitation, air temperature, and solar radiation are shown in Fig. 2 .

Because specific parameter values were unavailable for the specific site, the simulated domain was hypothetical, two dimensional, 1000 m in length (x dimension), and 50 m in height (z dimension). The domain was discretized using Δx = 5 m and Δz = 0.25 m for 200 by 200 grid cells. This hypothetical domain was constructed with a slanted land surface with a 21-m drop across the 1000-m domain for a slope, Sf,x = 0.021 m/m. The land surface was characterized as open shrubland in CLM and vegetation parameters were assigned using the International Geosphere–Biosphere Program standard. The subsurface was characterized as two heterogeneous units: a 2-m-thick alluvium that follows the ground surface and a fractured tuff unit below. The saturated hydraulic conductivity of the alluvium was represented as a correlated, Gaussian random field, which was generated using a parallel turning bands approach (Tompson et al., 1989, 1998; Ashby and Falgout, 1996). In this approach, which is common in the stochastic hydrogeology community (e.g., Rubin, 2003), the spatial variation of hydraulic conductivity is represented as a statistically stationary, random field where
and 〈lnKs(x)〉 = F, 〈f(x)〉 = 0 with 〈f(x)2〉 = σF2 = σlnK2, and Kg = exp(F). The mean or expected value is signified by angle brackets, Kg is the geometric mean of the hydraulic conductivity, and σ2 is the variance. The correlation structure (Rffe)for any two hydraulic conductivity values separated by a distance ξ, with a correlation scale λ, is represented by an exponential of the form (e.g., Rubin, 2003)

The heterogeneous parameters used for the alluvium were λx = 50 m, λz = 1.75 m, Kg = 0.0417 m/h, and σF2 = 1.0. For all cases, the numerical turning bands parameters were specified: 150 lines were used with the minimum grid spacing used in the line process rζ = 5.0 and the normalized frequency increment Δk = 0.2 (Tompson et al., 1989). Additionally, the principal direction heterogeneity was dipped to follow the 2% slope of the land surface.

The fractured tuff unit was represented using a stochastic continuum approach, where discrete fractures are generated and then mapped to create a heterogeneous saturated hydraulic conductivity field. Seven hundred fractures were distributed throughout the domain for an average fracture density of 0.7 fractures/m. Fracture lengths were probabilistic and were calculated using the approach of Reeves et al. (2008):
where the probability of a fracture of length l is dependent on the minimum fracture length w and the exponent a. A minimum fracture length of w = 24 m and exponent of a = 1 were specified. Because the true fracture orientation at the site was not known, all fractures were oriented at a ±45° angle with an equally likely chance of falling in either direction. A 45° fracture angle was used to provide realistic fracture networks; other fracture orientations are possible in this approach and might change local flow paths (see, e.g., Başağaoğlu et al., 2009) although other studies have found fracture orientation and density to be less significant variables (e.g., Gleeson et al., 2009). An example fracture realization is shown in Fig. 3 . Using the techniques of Reeves et al. (2008) and Botros et al. (2008), fractures were mapped to a single continuum of saturated hydraulic conductivity. The matrix was assigned a saturated hydraulic conductivity of Kmatrix = 0.000001 m/h. The fractures were assigned Kfracture = 0.001 m/h. Hydraulic conductivity was assigned as a simple linear combination of matrix and fracture permeability, depending on the number of fractures present in each 5- by 0.25-m cell. The fracture locations were chosen randomly and multiple “realizations” of fracture networks, generated by a change of random seed, had the same global statistics. Each fracture was traced sequentially to ensure connection of the high-conductivity region (as in Reeves et al., 2008, Fig. 2–3; Botros et al., 2008, Fig. 2). This is the simplest approach used to upscale the discrete fractures to the continuum grid that does not modify subsurface variables according to fracture dimension. Other unsaturated upscaling approaches have been explored in the recent literature (e.g., Finsterle, 2000; McKenna et al., 2003; Ohman and Niemi, 2003) including those that modify pressure–saturation relationships according to fracture geometry (e.g., Liu and Bodvarsson, 2001) and those that provide an analysis of error in the representation of the fracture network (e.g., Wellman et al., 2009). A realization composed of the combined alluvial, tuff, and explicit hillslope land surface is shown in Fig. 4 .

Realistic but hypothetical porosity and van Genuchten pressure–saturation parameters (van Genuchten, 1980) were assigned separately for each of three zones, one for the alluvium and two for the tuff unit, corresponding to fracture and matrix regions. Effective porosity was assigned as ϕ = 0.25 for the alluvial layer, ϕ = 0.1 for the matrix (i.e., fracture-free) portion of the tuff, and ϕ = 0.001 for any tuff units containing fractures. The van Genuchten parameters were also assigned to these three zones with αalluvium = 2 (1/m), nalluvium = 2, and Sres_alluvium = 0.0001; αmatrix = 0.1 (1/m), nmatrix = 3, and Sres_matrix = 0.01; and αfracture = 2 (1/m), nfracture = 3, and Sres_fracture = 0.0001. This approach was intended to capture the major variability in these parameters but did not vary the van Genuchten properties as a function of K or ln(K) within the alluvial layer. From this base-case configuration, four equally likely realizations of subsurface properties were generated by changing the random seed used to determine fracture locations. Each of these realizations was then used in a fully coupled simulation that was forced by the 2006 forcing data set repeatedly for five contiguous years. All simulations were initialized to very dry conditions using a hydrostatic pressure gradient with a pressure of −10 m specified at the bottom of the domain.

Two sensitivity test cases were also performed with perturbed parameter values. A homogeneous two-layered case was generated by assigning each layer the geometric mean of its hydraulic conductivity field from the base case described above (Kalluvium = 0.0417 m/h and Ktuff = 2.3242 × 10−5 m/h) and van Genuchten properties averaged between the two units (αtuff = 1.1 [1/m], ntuff = 3), and Sres_tuff = 0.0005, and ϕ = 0.05. All properties except tuff porosity and the van Genuchten functions were specified the same as in the base case. An additional fractured case with lower fracture density (0.35 fractures/m) was also generated, with all other properties kept the same as the base case. As detailed above, these three cases were also simulated with the wet-year forcing data with double the rainfall (but the same timing), with all other variables the same. This wet year forcing was used for 2 yr of simulation time.

The WET and DRY cases were used to understand the patterns in the system for drier and wetter than average years. The real system is dynamic and subject to interannual variability in meteorology and is not forced uniformly by a single, transient year of forcing and the real system is also not initialized completely dry at the start of a year. While spun-up or equilibrium simulations might have been used here (e.g., Maxwell and Miller, 2005; Maxwell et al., 2007; Kollet and Maxwell, 2008; Maxwell and Kollet, 2008; Kollet, 2009; Ferguson and Maxwell, 2010), work has shown that arid systems with very deep vadose zones are not in this type of equilibrium (e.g., Walvoord et al., 2002, 2004). Therefore, the system here was not brought into a complete equilibrium, given the range of interannual variability at the site, because the true system is transient and probably always out of equilibrium and provides a range of initial states (very dry to relatively wet for this system as produced by the multiple repeated years) for both WET and DRY forcing.

Water Balance Components

Figure 5 shows the saturation fields after 4.5 yr of simulation for two realizations of the heterogeneous base case along with the homogenous case for the dry simulation forcing. Although the spatial patterns of saturation are different between the two base-case realizations, the overall behavior is quite similar. Both cases demonstrate highly focused recharge in the fractured regions and very little recharge in the matrix, something noted by previous studies (e.g., Salve, 2005; Zhou et al., 2006). Because this figure is from June in the simulation year, the effects of summer ET are evident, as the top of this layer (which is in connection with the atmosphere) is close to residual saturation. Regions of approximately 50% saturation in the deeper alluvium, however, are evidence of the effects of heterogeneity in the alluvial layer. Despite the dry conditions at the land surface, very wet conditions are present deeper in the fractured system. This is shown by the isolated blue and green regions, which indicate saturations close to 100% (perched water) in many locations. By contrast, the unfractured regions exhibit no appreciable change in saturation (recall that the system was initialized as hydrostatic and very dry) and some fractured tuff regions remain at their initial saturation.

This figure also clearly shows differences in infiltration between the heterogeneous (Fig. 5A and 5B) and homogeneous (Fig. 5C) cases. The homogeneous case exhibits very dry conditions in the alluvial layer with very little infiltration into the tuff unit. There is a region of increased soil saturation (approximately 50%) at the alluvial–tuff interface, as shown by the yellow region in Fig. 5C. In contrast to the heterogeneous cases, infiltration in the homogeneous case is limited to the alluvial system.

Plots of soil saturation are shown in Fig. 6 for base-case Realization 1 under WET forcing conditions on 1 February (Fig. 6A), 1 March (Fig. 6B), and 1 April (Fig. 6C) of the second simulation year. Note that this is the same subsurface realization shown in Fig. 5A. Figure 6 clearly shows the pattern of recharge in the heterogeneous fractured system. On 1 February (Fig. 6A), the system is predominately dry, with isolated regions of higher saturation in the fractured tuff. By 1 March (Fig. 6B), the system is being recharged primarily by spring snowmelt, as shown by the change in atmospheric temperature and precipitation during this time frame in Fig. 2. The effects of the Gaussian heterogeneity are evident in Fig. 6B, as shown by the variations in soil moisture in the alluvial unit, and significant recharge is occurring in the tuff unit, as shown by the fully saturated (blue) regions. By 1 April (Fig. 6C), spring precipitation occurs as rain (as shown by the precipitation and warmer temperatures in Fig. 2) and contributes to significant recharge in both the alluvial and tuff units. While the upper portions of the alluvial unit are dry, isolated, saturated conditions exist at the tuff interface. Significant, focused recharge is also evident in the fractured tuff regions, with infiltration fronts moving up to 40 m vertically. No perceptible change in recharge is seen in the matrix regions of the tuff unit.

The change in water storage (surface and subsurface storage terms in Eq. [1–2]) is shown in Fig. 7 for all DRY simulations. It is important to note that due to the different saturation–pressure relationships and porosities, each simulation had a different initial storage, despite having the same pressure initialization. Because of this, it was important to calculate and plot only the change in storage. Figure 7 shows an increase in total subsurface storage during the 5-yr simulation in all cases. Storage also exhibits a clear seasonal cycle, with considerable recharge during the late winter and spring months and strong drying in the summer, with (temporally) isolated summer recharge associated with summer precipitation events.

For the first half of the first year, Fig. 7 shows very little difference in changes in storage between realizations or between the base and sensitivity cases. During the summer and fall of the first year, however, we see significant reduction in storage overall and a divergence between cases. This divergence increases during the remaining four simulation years, resulting in large differences in subsurface storage between cases after the third year. The homogeneous case demonstrates the least change in storage (and thus the least recharge), while three realizations of the heterogeneous, fractured base case show the greatest recharge, with >100 mm during the 5 yr, compared with <80 mm for the homogeneous case. This is due to the role of subsurface heterogeneity, as is clearly shown in Fig. 5. The fractured system allowed deeper percolation of infiltrated water below the root zone (2.5 m in this case). Infiltration in the homogeneous case was limited to only the shallow alluvial layer and remained within the root zone. It should be noted that the entire subsurface is connected hydraulically: even though water is below the root zone, fluxes across the land surface may pull water from regions deeper than the root zone itself. This is moderated, however, by the relative permeability function (as the fractures get drier they become less permeable) and by limitations on ET by water stress in the shallow subsurface.

Figure 8 plots the change in storage (as in Fig. 7) for the WET simulations. This plot shows a significant change in storage across all cases, with a maximum of >300 m in 2 yr. Figure 8 suggests that while wet years are infrequent, they contribute significantly to recharge. This figure shows similar differences between cases as seen in Fig. 7, starting after 0.5 yr of simulation time. Similar to the previous figure, the heterogeneous realizations show the greatest recharge and the homogeneous case shows the least. Differences between cases are smaller under WET conditions, however, with <20 mm of difference between simulations.

It is important to emphasize that all cases show clear infiltration, even for the DRY forcing, an unexpected result: even under DRY conditions, the average yearly storage increases by approximately 20 mm, indicating that approximately 10% of the 253 mm of annual precipitation is recharging the subsurface. Note that all 5 yr of the simulation do not show equal recharge amounts, indicating that the initial system state influences the infiltration and land-energy fluxes. The WET cases recharged approximately 170 mm/yr, or approximately 33% of the 506 mm of precipitation. This demonstrates a strong nonlinearity, further emphasizing the role of infrequent, but above-average precipitation and recharge, important because the precipitation record for the MEDA Station 12 shows significant annual variability. No overland, lateral outflow was generated by any simulation case, indicating that all precipitation contributes to either recharge or evaporation.

Time series of the monthly mean and variance of latent heat flux (spatial mean and variance across the simulated domain) are shown in Fig. 9 for the four base-case heterogeneous realizations under DRY forcing conditions. There is almost no difference between realizations except during the late summer months when there are slight differences, when the shallow subsurface is moisture limited. Differences are similar to those shown in Kollet (2009), where the largest differences in latent heat flux were also seen in late summer. The largest differences between the homogeneous and heterogeneous cases also occurred during this time (not shown). Spatial variability in latent heat flux is quite large. This is the variance calculated spatially, not temporally, because the mean monthly value for each location was calculated before spatial statistics. The maximum variance occurred in July, while the maximum latent heat flux occurred in June. Figure 9 also shows larger differences between realizations in the latent heat flux variance than for means.

Time series of the monthly mean and variance of latent heat flux for the four base-case realizations forced under WET conditions are shown in Fig. 10 . Contrasting Fig. 9 and 10 demonstrates that the WET cases were much less water limited than the DRY cases. In Fig. 9, there are many points throughout the year when latent heat flux is reduced due to moisture constraints, and monthly mean latent heat flux is closely tied to precipitation. Figure 10 shows two seasonal peaks in latent heat flux, one after the spring snowmelt and one during midsummer. In general, however, the peak latent heat flux is similar between the WET and DRY cases, as shown by a maximum of 50 W/m2 in Fig. 10 and just over 40 W/m2 in Fig. 9. This indicates that during June, both cases had sufficient soil moisture such that ET was close to potential values. The variance in latent heat flux for the WET case is also quite large, but appears to peak when the mean reaches a maximum value (Fig. 10), not later as is seen in Fig. 9. This shift in peak ET mean and variance contributes to the nonlinear differences in recharge estimates between DRY and WET cases discussed above. Peak ET (as shown by the latent heat flux) was very similar between cases, while precipitation was doubled.

The variance in shallow soil moisture should explain the variance in latent heat flux, with the greatest variances in latent heat flux occurring during moisture-limited times. Monthly averaged shallow soil saturation for the DRY cases (not shown) was generally very low, never >30%, and demonstrates a general pattern that is similar to the mean latent heat flux shown in Fig. 9. The saturation variance, however, was quite small, with a maximum under 0.3%. The saturation variance also peaked when the saturation peaked in April, not later in the year as seen for the latent heat flux. The monthly mean and variance of shallow soil moisture for base-case heterogeneous realizations under WET conditions (not shown) showed a peak in April, with large differences between realizations during this time, and greater saturation than the DRY cases in the second half of the year.

Spatial Correlations

Mean latent heat flux (Fig. 9 and 10) and mean soil saturation exhibited similar seasonal trends; however, the magnitude and timing of their variances differed significantly. If it is assumed that shallow soil saturation closely follows heterogeneity, then spatial statistics, namely the calculated experimental semivariogram, should follow the theoretical semivariogram, Eq. [8], used to generate the random alluvial saturated hydraulic conductivity fields. To explore this further, experimental semivariograms, γ, and unit semivariograms, γ/σ2, were calculated.

Figure 11 shows the experimental semivariogram calculated for the yearly average soil saturation for the last year of simulation for the uppermost model cells. This figure displays curves for each of the four base-case realizations under DRY forcing conditions, the average across the four realizations, and the theoretical, exponential semivariogram with horizontal correlation length specified for the alluvial layer, 50 m. Figure 11 shows close agreement between the theoretical semivariogram and those calculated based on the soil saturations, indicating that indeed the shallow soil moisture is strongly influenced by the near-surface heterogeneity in K.

Figure 12 shows the experimental semivariograms calculated for ground surface temperature (Fig. 12A) and latent heat flux (Fig. 12B), again averaged across the final DRY simulation year as in Fig. 11. This figure also plots the ensemble mean and a theoretical semivariogram, but in this case the theoretical semivariogram has a much shorter correlation length of 20 m. Note that the theoretical variogram shown in this figure has a correlation length much shorter than the theoretical variogram used to generate the alluvial correlation structure. This indicates that, while the imposed spatial correlation in the alluvial heterogeneity strongly influences soil saturation, the ground surface temperature and latent heat flux are influenced by other factors than just the near-surface geologic heterogeneity. Figure 13 plots the variogram for monthly averaged latent heat flux in August of the second and fifth simulation years under DRY forcing conditions. The second simulation year demonstrates a correlation structure similar to the soil saturation and subsurface heterogeneity, with a correlation length very close to the theoretical value of 50 m. This correlation structure changed with time, however, as indicated by Fig. 13B, which shows a shift to the shorter correlation length of 20 m as seen in Fig. 12.

Differences between the experimental and theoretical correlation lengths in Fig. 11 and 12 and between Fig. 13A and 13B is an interesting result and suggests that shallow soil moisture exerts some, but not all, influence over energy fluxes. To explore these results further, unit semivariograms (normalized by the mean for comparison) of yearly averaged soil saturation for DRY forcing are plotted for the 10 uppermost model layers in Fig. 14 . These cells represent the root zone region where CLM evapotranspiration fluxes were extracted from the ParFlow soil column. At a vertical model discretization of 0.25 m, this represents eight cells in the alluvial layer (gray lines) and two cells in the fractured tuff unit (black lines). This figure shows two base-case realizations under DRY forcing conditions along with two experimental unit semivariograms, one with a correlation length of 20 m (dashed) and one with a correlation length of 50 m (solid). Correlation structure in the alluvial layers is very different than in the fractured tuff unit. Notably, the alluvial layers all exhibit longer correlation lengths, while the tuff units much shorter correlation lengths. By comparing Fig. 12 and 14, it is clear that the ground surface temperature and the latent heat flux are influenced by a combination of shallow and deeper soil saturations. As mentioned above, the root zone in CLM as used here is 10 cells deep (2.5 m) and clearly the influence of the fractured tuff unit is evident here. It is also apparent from Fig. 13 that latent heat flux evolves with time as soil moisture infiltrates into the subsurface.

These results are interesting given the temporal and episodic nature of recharge that we see in the model. As discussed above, the soil saturation distributions for all cases (WET and DRY) in the fractured system generated “pockets” of fully saturated water, even in the dry case, which is unexpected. Figure 15 plots monthly, vertical, unit semivariograms for the tuff unit. This figure plots four base-case realizations under WET (gray) and DRY (black) forcing conditions at a log–log scale. The variograms in this figure, unlike the previous ones, are plotted at log–log scale to emphasize correlations across a range of spatial scales. Log–log variograms are commonly used (Gelhar, 1993; Murato and Saito, 1999; Mela and Louie, 2001) along with other methods (e.g., power spectra) to diagnose spatial correlations that are indicative of nonstationary or fractal behavior. Inspection of this figure indicates that the spatial correlation (over the vertical) of soil moisture distributions exhibits multiple scales. The signature of these patterns also appears to vary seasonally, with the most linear (in log–log space) semivariograms seen in the wetter months (March–May). This figure shows that there is similarity between the WET and DRY cases in most months. The plots for the months of March and April, however, when a majority of recharge occurs, show the largest differences between the WET and DRY cases. During these months, the WET cases indicate the strongest multiscale behavior, primarily due to the strongly local infiltration down the fractured regions demonstrated in Fig. 6B and 6C.

This work used coupled model simulations, forced by real meteorological data at fine time scales across multiple years and for different realizations of a detailed, stochastic fracture model, to investigate infiltration and water balance in an arid setting. These simulations were interrogated using temporal and geospatial statistics to determine correlations between subsurface heterogeneity, soil saturation, and land-energy fluxes. The simulations presented here are for a two-dimensional, hypothetical system, with realistic heterogeneity and observed forcing at high temporal resolution for a fractured-rock mesa system at moderately high altitude. The coupled land surface (radiation, energy balance, and vegetation fluxes) and hydrologic processes represented here, however, provide a test bed for land-energy and hydrologic interactions in an arid system.

There are a number of conclusions that may be derived from these results regarding infiltration that confirm the findings of previous studies using more traditional modeling and observational approaches. This study also shows that infiltration into the subsurface is localized and highly transient, with the majority of recharge occurring during short-duration recharge events at a range of temporal scales (e.g., Seyfried and Wilcox, 1995; Seyfried, 1998; Salve, 2005; Sandvig and Phillips, 2006; Scanlon et al., 2006; Zhou et al., 2006; Duke et al., 2007). The relatively large recharge rates seen in these results (∼20 mm/yr), however, are greater than those presented previously (∼0.1 mm/yr) for the same, general, geographic region (Walvoord et al., 2002, 2004), which were arrived at through modeling and observations. It is quite important to note, however, that the previous, lower estimates of recharge were obtained in a very different setting (Yucca Flat) than the current study and the sites exhibit differences in geology, elevation, and form of precipitation. Thus further underscores the variability in recharge across elevation and geologic features within a region, as this range of recharge values falls within estimates for similar settings (Heilweil et al., 2007).

Additionally, soil moisture, ground temperature, and energy fluxes at the land surface were shown to exhibit coherent spatial patterns. Correlation lengths for soil moisture were shown to agree very closely with those imposed for hydraulic conductivity. Land-surface temperature and latent heat flux, however, demonstrated different, shorter correlation lengths that were more complicated. Moreover, correlation lengths for latent heat flux were shown to evolve with time, decreasing with subsequent simulation years. Ultimately, it appears that these correlation structures depend on both shallow and deeper soil moisture and patterns of heterogeneity.

Finally, the model of this fractured system appears to create complex vertical subsurface soil moisture patterns that exhibit spatial correlations at multiple scales. This is due to both the underlying fracture model and the episodic nature of the recharge. For these correlations, differences between realizations were shown to be small; however, there were some differences between the simulations forced with WET and DRY climate conditions, particularly during months of great recharge. This is also in line with other models of fracture flow (e.g., Pruess, 1998; Liu et al., 2005) and recent observations of spatial patterns of recharge (Green et al., 2009).

This work implies that fully integrated watershed models that explicitly represent complicated subsurface heterogeneity may be used to develop insights into spatial and temporal patterns of infiltration, recharge, and land–atmosphere interactions in arid environments. In particular, variations in land surface fluxes are quite complicated. These complications need to be addressed in observational studies in addition to developing relationships to more accurately represent these fluxes at larger scales or in simpler models. Finally, the fact that the land surface fluxes contain information about deeper soil structure is very intriguing and might mean that these fluxes may be inverted to provide estimates of heterogeneity below the shallow land surface.

This research was supported in part by the Golden Energy Computing Organization at the Colorado School of Mines using resources acquired with financial assistance from the National Science Foundation and the National Renewable Energy Laboratory. I also wish to acknowledge the review and discussions of this manuscript by I. Ferguson.

Freely available online through the author-supported open access option.