VZJ Journal of Natural Resources and Life Sciences Education
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (9)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Walvoord, M. A.
Right arrow Articles by Striegl, R. G.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Walvoord, M. A.
Right arrow Articles by Striegl, R. G.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Walvoord, M. A.
Right arrow Articles by Striegl, R. G.
Related Collections
Right arrow Vadose Zone Processes and Chemical Transport
Right arrow Coupled Flow/Transport Models
Published in Vadose Zone Journal 3:502-512 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

SPECIAL SECTION: UNCERTAINTY IN VADOSE ZONE FLOW AND TRANSPORT PROPERTIES

Constraining the Inferred Paleohydrologic Evolution of a Deep Unsaturated Zone in the Amargosa Desert

Michelle A. Walvoord*,a, David A. Stonestromb, Brian J. Andraskic and Robert G. Striegla

a U.S. Geological Survey, Lakewood, CO 80225
b U.S. Geological Survey, Menlo Park, CA 94025
c U.S. Geological Survey, Carson City, NV 89706
* Corresponding author (walvoord{at}usgs.gov).

Received 3 April 2003.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 Hydraulic, Chemical, and...
 Site Description
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Natural flow regimes in deep unsaturated zones of arid interfluvial environments are rarely in hydraulic equilibrium with near-surface boundary conditions imposed by present-day plant–soil–atmosphere dynamics. Nevertheless, assessments of water resources and contaminant transport require realistic estimates of gas, water, and solute fluxes under past, present, and projected conditions. Multimillennial transients that are captured in current hydraulic, chemical, and isotopic profiles can be interpreted to constrain alternative scenarios of paleohydrologic evolution following climatic and vegetational shifts from pluvial to arid conditions. However, interpreting profile data with numerical models presents formidable challenges in that boundary conditions must be prescribed throughout the entire Holocene, when we have at most a few decades of actual records. Models of profile development at the Amargosa Desert Research Site include substantial uncertainties from imperfectly known initial and boundary conditions when simulating flow and solute transport over millennial timescales. We show how multiple types of profile data, including matric potentials and porewater concentrations of Cl–, {delta}D, {delta}18O, can be used in multiphase heat, flow, and transport models to expose and reduce uncertainty in paleohydrologic reconstructions. Results indicate that a dramatic shift in the near-surface water balance occurred approximately 16000 yr ago, but that transitions in precipitation, temperature, and vegetation were not necessarily synchronous. The timing of the hydraulic transition imparts the largest uncertainty to model-predicted contemporary fluxes. In contrast, the uncertainties associated with initial (late Pleistocene) conditions and boundary conditions during the Holocene impart only small uncertainties to model-predicted contemporaneous fluxes.

Abbreviations: ADRS, Amargosa Desert Research Site • CMB, chloride mass balance [age] • EVT, enhanced vapor transport • FEHM, Finite Element Heat and Mass Transfer Code • SMOW, Standard Mean Ocean Water • UZ, unsaturated zone


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
Hydraulic, Chemical, and...
 Site Description
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
ARID ECOSYSTEMS ARE particularly sensitive to climate. High sensitivity to near-surface water balance coupled with slow hydraulic response at depth renders deep unsaturated zones potentially useful as archives of past climatic and vegetational conditions in arid regions (Edmunds and Walton, 1980; Edmunds and Tyler, 2002). Deep unsaturated-zone (UZ) vertical profiles record major shifts in both climate (Cook et al., 1992; Phillips, 1994) and vegetation (Walvoord et al., 2002a, 2002b) that extend back >100 kyr (kyr = thousand years) (Tyler et al., 1996). Such paleorecords include the physical, chemical, and isotopic composition of porewater. These parameters vary systematically with depth, and thus time. Extracting paleohydrologic information from UZ profiles requires making assumptions with large associated uncertainties. Techniques of paleohydrologic reconstruction have advanced significantly in the last several years (Edmunds and Tyler, 2002). Even so, obtaining a unique paleohydrologic reconstruction from UZ data still remains challenging (Walvoord et al., 2002b; Scanlon et al., 2003).

The present study introduces an approach for integrating deep UZ stable isotope data into the paleohydrologic reconstruction. Incorporating multiple lines of UZ data, including independent paleoclimate proxy records, appears to be an effective approach to reducing and exposing nonuniqueness in reconstructions. A major goal of the integrated modeling is to constrain the timing of the last major paleohydrologic transition (mesic to arid) at the study site, since this variable strongly influences modeled UZ fluxes at all times (past, present, and projected). We examine the uncertainty associated with (i) specifying initial (Pleistocene) conditions derived from current chemical and isotopic conditions at depth and (ii) specifying boundary conditions throughout the Holocene on the basis of current values.


   Hydraulic, Chemical, and Isotopic Tracers
TOP
 ABSTRACT
 INTRODUCTION
 Hydraulic, Chemical, and...
Site Description
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Matric potential describes the freely exchangeable energy of porewater, due to capillary and adsorptive forces, relative to bulk water of the same composition, temperature, and elevation. Matric potential is negative for unsaturated materials, lower values indicating drier conditions. Instruments typically measure the sum of matric and osmotic (or solute) potentials. The gradient of the total potential, which is the sum of matric, osmotic, thermal, and gravitational potentials, indicates the direction of water movement. The slow hydraulic response of deep unsaturated zones following transitions from pluvial periods to xeric conditions enables matric-potential profiles to be used as a paleohydrologic tool, supplementing environmental tracers (Walvoord et al., 2002a).

Chloride, a commonly used conservative hydrologic tracer, is continually deposited on the land surface and moved into the subsurface with infiltrating rain (Allison, 1988; Phillips, 1994). Evapotranspiration concentrates porewater Cl–; higher concentrations denote drier conditions. Profiles of Cl– concentrations in arid interfluvial regions typically display a "bulge" in the upper 10 to 20 m, with peaks observed between 2 and 5 m (Phillips, 1994). The integrated amount of Cl– in the profile between the land surface and the base of the bulge divided by the average Cl– input rate yields the chloride mass balance (CMB) age of the accumulation (e.g., Phillips, 1994; Tyler et al., 1996). These calculations acquire the uncertainty of the average Cl– deposition assumed for millennial time scales (Scanlon, 2000). Uncertainties are linear; an assumed input that is twice the true average underestimates the CMB age by a factor of two. The CMB calculations also provide paleoflux estimates, using Cl– concentrations in the deepest (oldest) part of the profile and assumed paleoinputs of Cl– and precipitation.

The isotopic composition of porewater is a useful tracer of unsaturated flow in arid environments (e.g., Stonestrom et al., 1999). Heavier isotopic species of water such as HDO and H218O concentrate in liquid relative to gas during changes in phase (Dinçer and Davis, 1984). Here unsuperscripted H is 1H and unsuperscripted O is 16O. Stable isotopic compositions are reported as per mil deviations from Standard Mean Ocean Water (SMOW):

[1]
where Rsample is the isotopic ratio between heavy (HDO or H218O) and light (ordinary H2O) species in the sample, and RSMOW is the corresponding ratio for SMOW. Profiles of {delta}D and {delta}18O have been used to infer evaporation rates (Barnes and Allison, 1988) as well as shifts in climate over decadal time scales (Gaye and Edmunds, 1996). Higher temperatures and increased evaporation at the land surface lead to enriched (heavier) isotopic compositions. Within arid UZ profiles, heavier isotopic species diffuse downward against an upward evaporative flux. Considerable research has focused on interpreting shallow stable isotope profiles. In contrast, comparatively little attention has been given to deep profiles. Tyler et al. (1996) speculated that isotopically depleted water deep in the UZ at Frenchman Flat, NV represents Pleistocene-age infiltration during a cooler, wetter climate. Our modeling invokes and tests this concept.


   Site Description
TOP
 ABSTRACT
 INTRODUCTION
 Hydraulic, Chemical, and...
 Site Description
MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
The Amargosa Desert Research Site (ADRS) is located in an alluvial basin immediately east of Death Valley National Park, 17 km south of Beatty, NV (Fig. 1 ; Andraski and Stonestrom, 1999). The present climate is arid-thermic, moderately seasonal with respect to precipitation and strongly seasonal with respect to temperature (Bull, 1991). Precipitation averages 108 mm yr–1 (based on 1981–2001 records), of which about 70% arrives in frontal systems October through April. Daily temperatures average 33°C for the hottest month (August) and 3°C for the coolest month (December). The Amargosa Desert is part of the Mojave Desert ecosystem. The dominant plant is creosote bush [Larrea tridentata (Sessé & Moc. ex DC.) Coville] (Smith et al., 1997). Depth to groundwater at the ADRS ranges from 85 to 115 m (Fischer, 1992). Sediments consist almost entirely of unconsolidated and poorly sorted debris-flow, fluvial, and alluvial-fan deposits (Nichols, 1987; Stonestrom et al., 2003). The unsaturated zone beneath the surface soil is almost entirely sand and gravel, with minor amounts of silt and clay. Cores from boreholes used in this study that were analyzed for texture all classified as gravelly to very gravelly sand to sandy loam (Fig. 2) .



View larger version (53K):
[in this window]
[in a new window]
  		Fig. 1. Location of Amargosa Desert Research Site in the Mojave Desert of southern Nevada.

 


View larger version (40K):
[in this window]
[in a new window]
  		Fig. 2. Textural analyses of the <2-mm fraction of samples from Borehole UZB2 at the Amargosa Desert Research Site (Prudic et al., 1997) at depths of 2, 3, 5, 6, 9, 15, 18, 21, 27, 36, 48, 60, and 85 m.

 
High-stands of now-dry pluvial lakes together with floral and faunal evidence indicate that the current hot, arid climate has been widely prevalent in the Great Basin–Mojave Desert for at least the last 9 to 11 kyr (Grayson, 1993). Paleogroundwater levels at Devils Hole, 50 km from the ADRS (Fig. 1), suggest that regional drying began as early as 18 to 20 kyr ago (Szabo et al., 1994). Vegetation and glacial histories indicate that temperatures during the late Pleistocene were cooler throughout the Mojave Desert by 5 to 7°C (Grayson, 1993). Precipitation was substantially higher, but by uncertain amounts.


   MATERIALS AND METHODS
TOP
 ABSTRACT
 INTRODUCTION
 Hydraulic, Chemical, and...
 Site Description
 MATERIALS AND METHODS
RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Field and Laboratory
Water content profiles were determined from core samples and neutron logging (Andraski, 1997; Prudic et al., 1997). Water potential profiles were determined by a water-activity meter (Gee et al., 1992) in core-sample material and by in situ thermocouple psychrometers (Stonestrom et al., 1999; Tambusch and Prudic, 2000; Andraski and Scanlon, 2002). The stable isotopic composition of pore water was determined by isotope-ratio mass spectroscopy using cryodistilled pore water from core samples and condensed water vapor from unsaturated-zone gas samples (Prudic et al., 1997). Pore water Cl– profiles were determined by ion chromatography (Stonestrom et al., 2003). Hydraulic properties were determined as described in Andraski (1996) and Smith et al. (1999).

Modeling
To interpret the UZ profiles, one-dimensional numerical simulations were run using the Finite Element Heat and Mass Transfer Code (FEHM) (Zyvoloski et al., 1997) following the general approach of Walvoord et al. (2002b). FEHM simulates nonisothermal, multiphase, multicomponent flow and transport in porous media. The version of FEHM used for ADRS simulations included {delta}D and {delta}18O transport, as described in Wolfsberg and Stauffer (2003). Individual isotopic species of water are treated as separate chemical components that partition between vapor and liquid phases and move by advection and diffusion. An assumption inherent to the model and supported by data at the ADRS (Prudic et al., 1997) is that liquid and vapor in the unsaturated zone are in isotopic equilibrium with respect to D and 18O. We conducted a series of simulations for an 80- to 110-m unsaturated vertical column consisting of the dominant sediment type in the profile, corresponding most closely to the coarse sandy loam of Andraski (1996). A ladder-type finite-element mesh with variable nodal spacing ranging from 0.1 to 1.0 m was used to achieve finer resolution near the boundaries. Soil properties were based on values provided in Andraski (1996) except that the saturated hydraulic conductivity for our base case is one order of magnitude greater than the value reported in Andraski (1996). The latter was determined using repacked <2-mm material from fairly shallow depth, and values were empirically corrected for rock-fragment content. An analysis of in situ air-pressure responses to barometric fluctuations at depths ranging from 5 to 108 m resulted in an average air permeability that corresponded to a saturated hydraulic conductivity 13 times greater than this laboratory value (Smith et al., 1999).

The assumption of a uniform soil column enables straightforward evaluation of the uncertainty related to initial and boundary conditions, which was the intent of our study. Sediment heterogeneity introduces dominant uncertainties in certain situations (e.g., Fogg et al., 1998). Yet, a modeling study of unsaturated flow at the nearby Nevada Test Site found that uncertainty in fluxes related to sediment heterogeneity was generally smaller than that related to boundary conditions (Wolfsberg and Stauffer, 2003).

Base-case simulations were run using the soil parameters given in Table 1 and applying the initial and boundary conditions described below and given in Table 2. The change from initial (Pleistocene) to current (Holocene) boundary conditions represents a major transition in near-surface hydraulic conditions that is linked to both paleoclimate and vegetational changes (Phillips, 1994; Walvoord et al., 2002a). Uniform percolation supported by a wetter, cooler climate in the past (Grayson, 1993) describes the initial flow regime. The change in boundary conditions at t = 0, subsequently referred to as the hydraulic transition, reflects the onset of drying induced by the shift to xeric climate and vegetational succession near the end of the Pleistocene (Spaulding, 1990; Forester et al., 1999). Multiple types of profile data are used to evaluate the effects of uncertainties in (i) the hydraulic-transition time and (ii) initial and boundary conditions. To evaluate the effects of these uncertainties on profile development and model-predicted water fluxes, initial and boundary conditions are varied within reasonable ranges.


View this table:
[in this window]
[in a new window]
  		Table 1. Site characteristics and average soil properties measured and used in model simulations.

 

View this table:
[in this window]
[in a new window]
  		Table 2. Initial and boundary conditions for base simulation.

 
Initial conditions are estimated using deep UZ data (Table 2; data in Fig. 3) . The nearly uniform Cl– concentrations below the 14-m depth and the uniform {delta}D and {delta}18O values below the 40-m depth reflect late Pleistocene values for a relatively small period of time ({approx}3 kyr based on Cl– mass) before the major Pleistocene–Holocene climate shift. These values are specified as initial compositions for the entire unsaturated zone. The greater depth to uniform values of {delta}D and {delta}18O reflects the significant vapor-phase transport of water, which has no counterpart for Cl–. Initial percolation rates of 2.4 and 4.8 mm yr–1 are evaluated. These rates are estimated by dividing 82 and 164 mg Cl– m–2 yr–1, the low and high Cl– deposition estimates of Prudic (1994), by 34 mg L–1, the average Cl– concentration deep in the profile (Fig. 3A).



View larger version (50K):
[in this window]
[in a new window]
  		Fig. 3. Model-predicted (A) Cl– (B) matric potential (C) {delta}18O, and (D) {delta}D profiles at various simulation times since the initiation of drying conditions. Amargosa Desert Research Site data are shown for comparison. The matric potential data in (B) are calculated from water and osmotic potential measurements. Isotopic values at 72.5 m (C and D) are suspect due to a crack in the paraffin covering of the sample and anomalously low water content noted at the time of analysis, both suggesting evaporative loss of water during storage (Prudic et al., 1997).

 
For the base-case simulation, upper and lower hydraulic, thermal, chemical, and isotopic concentration boundary conditions for the Holocene are assigned by applying contemporary values (Table 2). The upper boundary for fixed values of matric potential, {delta}D, and {delta}18O is positioned at 2.3 m below ground surface to avoid temporal fluctuations at shallower depths. Because FEHM does not explicitly represent the osmotic component of total water potential, matric-potential boundary conditions are prescribed. For comparison between model results and observations, water potential measurements are converted to matric potentials by subtracting the osmotic potential. Osmotic components of water potential at the ADRS comprise up to 11% in the upper 10 m for current conditions but <1% deeper in the profile. Continuous monitoring from 1987 to 1992 at the ADRS indicated that wetting fronts from infiltrating rain penetrated <1.2 m under native vegetation (Andraski, 1997; Tambusch and Prudic, 2000). The data also show that plants maintain extremely low matric potentials at the base of the root zone (less than or equal to –400 m), preventing deeper percolation. The cessation of percolation deeper than 2.3 m is a key feature of the simulations; however, the location of the peak Cl– concentration at 2.3 m requires occasional penetration of liquid water to that depth. Bimonthly measurements of {delta}18O and {delta}D in shallow soil water, water vapor, and creosote stem water indicate a considerable range in {delta}18O and {delta}D (E. McConnaughey, 1995, unpublished data; Fig. 4) . The temporal variability expressed in the bimonthly data suggests that {delta}18O and {delta}D fluctuations propagate below 1 m, but the exact extinction depth of seasonal variations is uncertain given the paucity of deeper data. We assume a propagation depth of 2.3 m to coincide with the sub-root zone boundary and assign interpolated values of 0{per thousand} {delta}18O and –60{per thousand} {delta}D at this boundary (Table 2). A simple step change in stable isotope concentration is assumed in the base-case simulation, consistent with increased temperatures and evaporation associated with the hydraulic transition from Pleistocene to Holocene conditions.



View larger version (34K):
[in this window]
[in a new window]
  		Fig. 4. Shallow liquid-phase (A) {delta}18O values, and (B) {delta}D values measured at, or adjacent to, the ADRS at specified times between 1994 and 1999. Creosote stem water reflects a mixture of isotopic water compositions in the root zone.

 

   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 Hydraulic, Chemical, and...
 Site Description
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
CONCLUSIONS
 REFERENCES
 
Base Simulations
Results from applying initial and boundary conditions in Table 2 for different lengths of time illustrate the simultaneous development of Cl–, matric potential, and stable isotope profiles (Fig. 3). Simulated profiles mimic the general shape and magnitude of measured profiles, lending support to the overall conceptual model. The strong dependence of profile evolution on time since the hydraulic transition imposes multiple constraints on its occurrence, exposing uncertainty that otherwise is hidden. Considering time as the only variable, the best fit between measured and modeled data differs for isotopic, hydraulic, and chemical data. For example, simulation times between 8 and 16 kyr yield reasonable matches to the {delta}18O and {delta}D data, but modeled matric potentials profiles better match observed data given 32 kyr (and longer) of continuous drying (Fig. 3B–3D). Given the range of Cl– deposition estimates (Prudic, 1994), modeled Cl– profiles best match observations between 16 and 32 kyr (Fig. 3A). Within these broad ranges there is no single transition time that wholly satisfies all lines of data, arguing against a simple, synchronous transition.

Independent paleoclimate proxies supply supplemental information about the hydraulic transition in the region. These records also differ on the exact timing of changes in effective moisture (Szabo et al., 1994) because of a general reduction in precipitation and increase in temperature that was not necessarily synchronous (Van Devender and Spaulding, 1979). Despite differences in timing, most paleoclimate proxies from southern Nevada consistently indicate late Pleistocene temperature that was cooler by 3 to 7°C, precipitation that was higher by 110 to >200%, and vegetation that was dominated by pinyon-juniper woodlands at comparable elevations in the northern Mojave Desert (Van Devender and Spaulding, 1979, Van Devender et al., 1987; Spaulding, 1990; Grayson, 1993; Forester et al., 1999). Ancient packrat middens reveal a shift from juniper woodlands to xeric desert scrub between 17 and 14 kyr ago at elevations relevant to the ADRS (Spaulding, 1990). The presence of woodlands until at least 17 kyr ago reduces the likelihood of sub-root zone drying conditions before this time, as suggested by the matric potential modeling results. Alternatively, Scanlon et al. (2003) suggested that enhanced vapor transport (EVT) should be included in simulations of unsaturated zone flow. Experimental vapor fluxes reportedly exceed those predicted by Fick's Law by factors of 3 to 10 at water contents relevant to the ADRS (Philip and de Vries, 1957; Cass et al., 1984). To test whether EVT could produce better model fits to the matric potential data given drying times consistent with paleovegetation records, EVT was incorporated into FEHM (Walvoord et al., 2002a). A fixed enhancement factor of four produces a matric potential profile after 16 kyr that is very similar to the profile produced after 32 kyr without EVT (Fig. 5A) . Matches between modeled and measured matric potentials are also sensitive to representation of hydraulic properties (e.g., Fig. 5B). Due to the extreme nonlinearity between matric potential and water content, it is not surprising that matric potential profiles are difficult to match. Nevertheless, it is important to capture the general trend of the data, which the model achieves.



View larger version (32K):
[in this window]
[in a new window]
  		Fig. 5. Comparison of model-predicted matric potentials (A) with and without enhanced vapor transport (EVT), and (B) with variable saturated hydraulic conductivities (Ksat). The 32-kyr profile without EVT and the 16-kyr profile with 4x EVT completely overlap (A). t = 16 kyr for all profiles shown in (B).

 
The shorter development time for the {delta}18O and {delta}D profiles relative to the Cl– and matric potential profiles (Fig. 3) reflects the abrupt transition assumed in the model. A recent paleotemperature study in the Bonneville Basin of Utah indicates that mean annual Holocene temperatures reached maximal values only in the last 5.8 kyr (Kaufman, 2003). If this result applies at the ADRS, the shallow isotopic composition is more enriched at present than it has been during most of the time since the hydraulic transition. The stable-isotope upper-boundary condition is highly dependent on surface temperatures and may not coincide with the hydraulic transition expressed in Cl– and matric potential profiles. A delayed shift in boundary temperature would be a dominant factor in determining present-day isotopic profiles. Chloride and matric potentials mainly reflect a dramatic shift in near-surface water balance that could be initiated by decreased precipitation and botanical succession to xeric vegetation, with only slight increases in temperature. The later arrival of notably warmer temperatures apparently initiated the enrichment in shallow soil {delta}18O and {delta}D compositions subsequent to the near-surface water-balance change.

Estimation of elapsed time since the hydraulic transition influences the prediction of contemporary UZ water fluxes (Fig. 6) . Unsaturated zone fluxes decrease in absolute magnitude through time to approach an upward net (liquid plus vapor) water flux of approximately 0.01 mm yr–1 at the ADRS, which is controlled mainly by the geothermal vapor flux (Walvoord et al., 2002a). Figure 7 illustrates the transition in direction of net flux just above the water table from downward at 8 kyr to approximately 0 at 13 kyr to upward at 16 kyr. These results have implications for water-resource availability. Given that the area is almost entirely dependent on groundwater, constraining the timing and amount of recharge is important.



View larger version (31K):
[in this window]
[in a new window]
  		Fig. 6. Model-predicted sub-root zone water (liquid and vapor) fluxes for simulation times of (A) 8 kyr, (B) 13 kyr, and (C) 16 kyr since the hydraulic transition to drying conditions. Negative values indicate downward fluxes. Net = liquid + vapor.

 


View larger version (38K):
[in this window]
[in a new window]
  		Fig. 7. Model-predicted recharge (downward flux across the water table) with time since the hydraulic transition. (A) Graph with a logarithmic ordinate shows the rapid decline in water table recharge. (B) Graph with a linear ordinate shows the transition from a downward to an upward net flux.

 
The base-case simulations indicate that the simple strategy of applying modern boundary conditions to the past 16 kyr yields approximate matches to observed hydraulic and tracer profiles. Assumed initial and time-invariant boundary conditions are among the largest sources of uncertainty in evaluating current arid-region UZ water fluxes (Scanlon, 2000). To examine this issue further, we conducted sensitivity analyses on initial and boundary conditions.

Effects of Initial-Condition Uncertainties
Initial (late Pleistocene) conditions represented in the model include the depth to the water table and the percolation rate, which in turn determine the steady matric and water content initial profiles. They also include the initial temperature profile as well as the initial distributions of Cl–, {delta}D, and {delta}18O. Assuming that UZ profiles preserve paleowater compositions at depth, the initial Cl–, {delta}D, and {delta}18O distributions are well constrained by measurements. Thus, our analysis of sensitivity to initial conditions focuses on the initial depth of the water table, the initial temperature profile, and the initial percolation rate.

Water-table elevations near the ADRS in the late Pleistocene almost certainly exceeded current levels. Paces and Whelan (2001) concluded that late Pleistocene water tables were <30 m above modern levels 12 km southeast of the ADRS. Szabo et al. (1994) inferred that water-table elevations at Devils Hole, about 50 km from the ADRS (Fig. 1), were 9 m higher approximately 20 kyr ago and 6 m higher about 17 kyr ago, followed by a rapid decline to near modern levels. We therefore assume a maximum water table elevation of 30 m above modern during the late Pleistocene. Simulations using this level until 16 kyr ago (t = 0 in the model) show nearly the same matric potentials and water fluxes as in the base-case simulation. Sedimentologic, palynologic, and packrat-midden data in the Las Vegas valley, 130 km southeast of the ADRS, suggest that even though the decline in effective moisture began about 14 kyr ago, approximately modern conditions were not attained until about 8 kyr ago (Van Devender and Spaulding, 1979; Quade, 1986). An extreme model scenario, with a 30-m water table elevation persisting until 8 kyr ago, causes a more rapid reduction in recharge (downward flux across the water table) during the high-water-table period (0 ≤ t < 8 kyr in the model) owing to the faster response of the thinner UZ (Fig. 8A) . Following the simulated drop of the water table 8 kyr ago, recharge increases due to the drainage of UZ water from the 80- to 110-m interval. Contemporary UZ fluxes (t = 16 kyr, Fig. 8B) and matric potential profiles (not shown) for the extended +30-m water table simulation resemble the t = 13 kyr of profile drying in the base-case simulation (Fig. 6B). Other less-extreme scenarios, with smaller water table increases and increases of shorter durations yield contemporary fluxes that fall between these end members.



View larger version (37K):
[in this window]
[in a new window]
  		Fig. 8. (A) Model-predicted recharge (downward flux across the water table) with time since the hydraulic transition for the base simulation (dotted line) and for the simulation with a +30-m water table between t = 0 and 8 kyr (dash-dot line). (B) Contemporary modeled sub-root zone water (liquid and vapor) fluxes for the +30-m water table simulation are nearly equivalent to t = 13 kyr fluxes for the base simulation. Negative values in (B) indicate downward fluxes. Net = liquid + vapor.

 
An initial temperature profile with a surface temperature ≤7°C cooler than the present average produces virtually no effect in model-predicted matric potentials, porewater composition, or net water fluxes as long as the temperature change precedes or coincides with the hydraulic transition. Perturbations in thermal-vapor fluxes due to the initially different thermal gradients persist <100 yr, the time required for temperatures to re-equilibrate to nearly linear profiles corresponding to the posttransition boundary values (Fig. 9) . Changes in thermal vapor fluxes are small relative to changes in liquid fluxes during the first 100 yr of the simulation.



View larger version (27K):
[in this window]
[in a new window]
  		Fig. 9. Temperature profiles following a step change in specified boundary values.

 
Varying initial percolation rates from 0.5 to 50 mm yr–1 affects model-predicted fluxes and matric potentials only 1 to 2 kyr after the hydraulic transition. At t > 2 kyr, modeled liquid fluxes, vapor fluxes, and matric potentials no longer depend on initial percolation rate within the tested range (one order of magnitude above and below the estimated paleopercolation rate). Figure 10 illustrates the effect of varying the initial percolation rate. Initial pore water composition is independently specified and hence insensitive to initial percolation rates.



View larger version (52K):
[in this window]
[in a new window]
  		Fig. 10. Model-predicted recharge (downward flux across the water table) with time for variable initial percolation rates.

 
Effects of Boundary Condition Uncertainties
In the model, a fixed sub-root zone matric potential ({psi}rz) at a depth of 2.3 m controls the near-surface water balance (Walvoord et al., 2002a, 2002b; Scanlon et al., 2003). Although matric potentials near the base of the root zone have remained in the range of–400 to –650 m during more than 5 yr of monitoring at the ADRS (Andraski, 1997; Tambusch and Prudic, 2000), it is likely that greater variations have occurred over the course of the Holocene. Figures 11 and 12 show the model sensitivity to {psi}rz. Discrepancies with measured matric potentials intensify for {psi}rz > –400 m compared with {psi}rz ≤ –400 m (Fig. 11). Even so, contemporary UZ fluxes, including recharge rates, exhibit minimal sensitivity to {psi}rz values in the range –600 to –100 m (Fig. 12). In contrast, simulated sub-root zone liquid fluxes switch from upward to downward when {psi}rz exceeds a threshold of about –100 m (Fig. 12A). Likewise, the direction of fluxes across the water table switches at the threshold {psi}rz and downward fluxes (recharge rates) increase rapidly with increasing {psi}rz values (Fig. 12B). Chloride profiles do not vary significantly from the base-case simulation until {psi}rz exceeds the threshold. Once {psi}rz exceeds –100 m, simulations yield increasingly dispersive Cl– profiles that deviate markedly from the measured data (Fig. 13) . On the basis of these results, we conclude that soil water has not percolated below the 2.3-m depth at the ADRS for many thousands of years. The combination of a highly developed drying profile (expressed, e.g., in the matric potential data of Fig. 3B), together with the slow response of sub-root zone matric potentials to episodic percolation reinforces this conclusion.



View larger version (35K):
[in this window]
[in a new window]
  		Fig. 11. Contemporary model-predicted matric potential profiles for variable specified sub-root zone matric potentials ({psi}rz).

 


View larger version (42K):
[in this window]
[in a new window]
  		Fig. 12. (A) Contemporary model-predicted shallow (just below the root zone) fluxes, and (B) water table fluxes for variable specified sub-root zone matric potentials ({psi}rz). Negative values in (A) and (B) indicate downward fluxes.

 


View larger version (28K):
[in this window]
[in a new window]
  		Fig. 13. Contemporary model-predicted Cl– profiles in the upper 20 m for variable specified sub-root zone matric potentials.

 
The time-invariant value assigned to average Cl– deposition following the hydraulic transition directly influences the paleohydrologic reconstruction, since we rely on the accumulation time of the Cl– bulge to determine its timing. As demonstrated above, it is this timing that establishes drying conditions beneath the root zone and largely controls modeled fluxes (Fig. 6 and 7). The Cl– deposition rate before the transition influences paleopercolation values used as initial conditions in our modeling approach. Paleopercolation values exert no influence on the timing of the hydraulic transition and little influence on modeled fluxes beyond about 2 kyr after the transition. While there are important uncertainties associated with Cl– deposition since the transition (Scanlon, 2000), studies of 36Cl/Cl– ratios in 14C-dated packrat middens suggest that Cl– deposition rates have remained fairly constant for at least the past 9.3 kyr (Fabryka-Martin et al., 2000). Furthermore, fluctuations in Cl– deposition since the hydraulic transition are effectively damped out within the zone of Cl– accumulation. Below this zone, diffusion (rather than advection) controls Cl– transport.

The assumption of time invariance for {delta}D and {delta}18O input at the surface subsequent to the hydraulic transition appears to be a less realistic simplification. The available data suggest that annual fluctuations in near-surface porewater {delta}D and {delta}18O values are not damped out entirely at 2.3 m (Fig. 4), as assumed. Fluctuations in precipitation patterns and surface temperature complicate the isotopic boundary condition, with effects propagating downward in the vapor phase as well as the liquid phase. Nine- to 10-kyr-old plant fossils from Pintwater Cave in southern Nevada are enriched 5 to 6{per thousand} in {delta}18O above present, suggesting an increased frequency of summer monsoons in the early Holocene (Jahren et al., 2001). We expect such low-frequency fluctuations to propagate much deeper than 2.3 m. However, observed deep {delta}D and {delta}18O profiles change monotonically with depth, within the errors of determination. Regardless of near-surface details, the deep {delta}D and {delta}18O profiles record a definite shift from a climate that was cool and wet to a climate that was warmer and had increased evaporation sometime between late-Pleistocene and early-mid Holocene times. Whether the shift was abrupt or gradual cannot be resolved with our approach. A series of incremental steps in the isotopic upper-boundary condition spread across the Holocene produces {delta}D and {delta}18O profiles that all reasonably match the data, falling between the 8- and 16-kyr profiles for the one-step simulations shown in Fig. 3. Given the limited constraints on stable-isotopic boundary conditions during the Holocene, various scenarios are consistent with the data. Future analyses will benefit from the development of near-surface stable-isotope histories extending throughout the Holocene (e.g., Jahren et al., 2001; Kaufman, 2003).


   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 Hydraulic, Chemical, and...
 Site Description
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
REFERENCES
 
Reconstructing the paleohydrologic evolution of deep unsaturated zones involves a high degree of subjectivity. Simultaneous use of hydraulic, chemical, and isotopic profiles helps reduce and expose uncertainties in conceptual models and their parameterization. We explored the ability of simple climatic conceptualizations, represented by stepwise changing boundary conditions, to match currently observed UZ profiles using a multiphase model of fluid, solute, and energy transport that accounts for isotopic species and exchanges between vapor and liquid phases. Our study focused primarily on uncertainties related to the timing of the climatically induced transition between the cooler, wetter climate of the late Pleistocene and the hotter, drier climate of today. It also assessed uncertainties due to assumed initial and boundary conditions. Modeling did not explore uncertainties related to sediment heterogeneity and multi-dimensional flow, which may be important in some settings. While the latter are felt to be smaller than the uncertainties considered for the ADRS, better quantification of relative uncertainties will require three-dimensional characterization and continuing modeling efforts.

Unknowns in paleoconditions include water-table elevation, subsurface temperatures, and percolation rates, yet only a small fraction of these uncertainties are imparted to simulated present-day fluxes. The long duration of UZ drying in interfluvial areas substantially diminishes the influence of Pleistocene conditions. Uncertainties in post-transitional boundary conditions, particularly in the upper-boundary matric potential and in the Cl– deposition rate, strongly affect the paleohydrologic reconstruction. Uncertainties in other Holocene boundary conditions produce relatively small uncertainties in modeled fluxes for equivalent elapsed times since the initiation of UZ drying. Because the timing of the hydraulic transition imparts the largest uncertainty to model-predicted contemporary fluxes, the most effective approach in constraining past and current fluxes is to combine multiple lines of UZ data with paleoenvironmental data that independently constrain the timing of the transition.

Application of this approach to the Amargosa Desert Research Site indicates that development of observed Cl– and matric potential profiles requires 16 to 32 kyr of drying under current conditions, whereas stable isotope profiles require higher temperatures and increased evaporation for only the last 8 to 16 kyr. Paleoclimate and paleovegetation data argue against UZ drying before about 16 kyr ago. Climatic proxies suggest that the major shift from cool, mesic conditions to warm, arid conditions occurred 10 to 16 kyr ago. Vegetational evidence suggests that a shift from juniper woodlands to desert scrub occurred 13 to 16 kyr ago. We conclude that the transition from downward flow throughout the UZ to upward flow below root zone occurred about 16 kyr ago due to the combination of drier conditions and resulting vegetational changes. This transition initiated the accumulation of Cl– and the propagation of a drying front that is observed today in the matric potential profile. The hydraulic transition was followed somewhat later by markedly warmer temperatures resulting in the observed isotopically enriched shallow porewater. Estimates of contemporary UZ fluxes depend sensitively on the timing of the hydraulic shift. The direction of net water (liquid + vapor) flux across the water table reversed from downward to upward 8 to 16 kyr after the hydraulic transition. At present, corresponding to 16 kyr after the hydraulic transition, the model predicts virtually no downward liquid flux across the water table and an upward net water flux of <0.01 mm yr–1 just above the water table. The absence of groundwater recharge under current conditions has important implications for groundwater resources as well as contaminant transport under natural-gradient conditions.


   ACKNOWLEDGMENTS
 
We gratefully acknowledge Chris Milly, Bridget Scanlon, John Tauxe, Isaac Winograd, and Andrew Wolfsberg for their insightful comments on earlier versions of this manuscript. This research was performed while M.A.W. held a National Research Council Research Associateship Award at the U.S. Geological Survey in Lakewood, CO. Financial support was provided by the U.S. Geological Survey through the National Research Program and the Toxic Substances Hydrology Program.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 Hydraulic, Chemical, and...
 Site Description
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 




This article has been cited by other articles:


Home page
 Vadose Zone JHome page
J. J. Gurdak, R. T. Hanson, P. B. McMahon, B. W. Bruce, J. E. McCray, G. D. Thyne, and R. C. Reedy
Climate Variability Controls on Unsaturated Water and Chemical Movement, High Plains Aquifer, USA
Vadose Zone J., August 23, 2007; 6(3): 533 - 547.
[Abstract] [Full Text] [PDF]

Home page
 Vadose Zone JHome page
E. M. Kwicklis, A. V. Wolfsberg, P. H. Stauffer, M. A. Walvoord, and M. J. Sully
Multiphase, Multicomponent Parameter Estimation for Liquid and Vapor Fluxes in Deep Arid Systems Using Hydrologic Data and Natural Environmental Tracers
Vadose Zone J., August 24, 2006; 5(3): 934 - 950.
[Abstract] [Full Text] [PDF]

Home page
 Vadose Zone JHome page
C. J. Mayers, B. J. Andraski, C. A. Cooper, S. W. Wheatcraft, D. A. Stonestrom, and R. L. Michel
Modeling Tritium Transport Through a Deep Unsaturated Zone in an Arid Environment
Vadose Zone J., October 10, 2005; 4(4): 967 - 976.
[Abstract] [Full Text] [PDF]

Home page
 Vadose Zone JHome page
B. J. Andraski, D. A. Stonestrom, R. L. Michel, K. J. Halford, and J. C. Radyk
Plant-Based Plume-Scale Mapping of Tritium Contamination in Desert Soils
Vadose Zone J., August 16, 2005; 4(3): 819 - 827.
[Abstract] [Full Text] [PDF]

Home page
 Vadose Zone JHome page
R. M. Holt and M. J. Nicholl
Uncertainty in Vadose Zone Flow and Transport Prediction
Vadose Zone J., May 1, 2004; 3(2): 480 - 484.
[Full Text] [PDF]

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
<