Published online 24 August 2006
Published in Vadose Zone J 5:934-950 (2006)
DOI: 10.2136/vzj2006.0021
© 2006 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: PARAMETER IDENTIFICATION AND UNCERTAINTY ASSESSMENT IN THE UNSATURATED ZONE
Multiphase, Multicomponent Parameter Estimation for Liquid and Vapor Fluxes in Deep Arid Systems Using Hydrologic Data and Natural Environmental Tracers
Edward M. Kwicklisa,
Andrew V. Wolfsberga,*,
Philip H. Stauffera,
Michelle A. Walvoordb and
Michael J. Sullyc
a Earth and Environmental Sciences Division, Los Alamos National Laboratory, Los Alamos, NM 87545
b USGS, Box 25046, Denver Federal Center, Mail Stop 413, Lakewood, CO 80225-0046
c Neptune and Company, Inc., 1505 15th St., Suite B, Los Alamos, NM 87544
* Corresponding author (awolf{at}lanl.gov)
Received 9 February 2006.
 |
ABSTRACT
|
|---|
Multiphase, multicomponent numerical models of long-term unsaturated-zone liquid and vapor movement were created for a thick alluvial basin at the Nevada Test Site to predict present-day liquid and vapor fluxes. The numerical models are based on recently developed conceptual models of unsaturated-zone moisture movement in thick alluvium that explain present-day water potential and tracer profiles in terms of major climate and vegetation transitions that have occurred during the past 10 000 yr or more. The numerical models were calibrated using borehole hydrologic and environmental tracer data available from a low-level radioactive waste management site located in a former nuclear weapons testing area. The environmental tracer data used in the model calibration includes tracers that migrate in both the liquid and vapor phases (
D, 18O) and tracers that migrate solely as dissolved solutes (Cl), thus enabling the estimation of some gas-phase as well as liquid-phase transport parameters. Parameter uncertainties and correlations identified during model calibration were used to generate parameter combinations for a set of Monte Carlo simulations to more fully characterize the uncertainty in liquid and vapor fluxes. The calculated background liquid and vapor fluxes decrease as the estimated time since the transition to the present-day arid climate increases. However, on the whole, the estimated fluxes display relatively little variability because correlations among parameters tend to create parameter sets for which changes in some parameters offset the effects of others in the set. Independent estimates on the timing since the climate transition established from packrat midden data were essential for constraining the model calibration results. The study demonstrates the utility of environmental tracer data in developing numerical models of liquid- and gas-phase moisture movement and the importance of considering parameter correlations when using Monte Carlo analysis to characterize the uncertainty in moisture fluxes.
Abbreviations: ET, evapotranspiration • GMWL, Global Meteoric Water Line • LLW, low-level radioactive waste • NTS, Nevada Test Site • RWMS, Radioactive Waste Management Site
|
INTRODUCTION
|
|---|
SUBSIDENCE CRATERS created from underground nuclear weapons tests in Area 3 of the Nevada Test Site (NTS) are being used to dispose of low-level radioactive waste (LLW) generated from defense and environmental activities at the NTS. To help evaluate the long-term potential for migration of these wastes from the disposal areas, estimates of liquid and vapor fluxes are necessary. However, direct measurement of these quantities is essentially impossible. Therefore, we have constructed numerical models to predict these fluxes by matching observable concentrations of natural tracers that reflect long-term liquid- and gas-phase transport processes near the disposal sites. The ability of these models to predict reliably the isolation capability of these disposal areas depends, in part, on the inclusion of the relevant transport processes in the models and on the accurate estimation of the model parameters and their uncertainties.
Recent studies of thick unsaturated zones in arid regions have demonstrated the utility of using natural environmental tracers in conjunction with hydrologic data to characterize long-term hydrologic processes in these environments (e.g., Phillips, 1994; Tyler et al., 1996; Walvoord et al., 2002a, 2002b; Scanlon, 2000; Scanlon et al., 2003; Wolfsberg and Stauffer, 2003). The concentrations and subsurface distribution of these tracers reflect the influences of past and present climates and related vegetation communities, as well as subsurface transport processes in both the liquid and gas phases. By calibrating a transient flow and transport model to field observations, we identify parameters necessary to calculate present day subsurface liquid and vapor fluxes in undisturbed locations of Area 3 and subsequent liquid- and gas-phase contaminant migration.
Current strategies for the management of the waste-disposal sites in Area 3 call for the restoration of long-term natural surface conditions (Bechtel Nevada, 1998). Filled disposal craters have been covered and revegetated, and open sites now receiving waste will be similarly closed. Therefore, it is meaningful to evaluate the character and magnitude of processes that have affected liquid water and water vapor movement under natural conditions as an indicator of future long-term conditions in the waste-disposal areas. Under the present-day arid climate (approximately the last 10000 yr), liquid water movement through the thick unsaturated zones typical of alluvial basins in the American southwest is estimated to be very small (Tyler et al., 1996; Walvoord et al., 2002a, 2002b; Wolfsberg and Stauffer, 2003). In these environments, vapor flux may be of a magnitude comparable to, or exceeding, liquid flux. Consequently, vapor transport processes are examined in this study as indicators of the gas-phase contaminant migration associated with waste-disposal sites.
To calibrate numerical models of liquid- and gas-phase transport under natural conditions, this study relies on the subsurface measurements of hydrologic variables, such as water potential, saturation, and water content, as well as naturally occurring groundwater tracers, such as chloride (Cl) and the isotopes of hydrogen (2H or D) or oxygen (18O) that become part of the water molecule (e.g., as HD16O and H218O). The models also rely on laboratory measurements of permeability, bulk- and grain-density, porosity, and parameters of the soil moisture functions relating water potential, liquid saturation, and relative permeability (Bechtel Nevada, 1998, Table 1–1). Throughout this paper, pluvial and postpluvial conditions are discussed. The time scales for these conditions correspond to the wetter climate at the end of the Pleistocene epoch and the dryer climate during the Holocene epoch. The transition between these epochs occurred approximately 10 000 yr ago. The uncertainty in the timing of the climate transition and the corresponding change in vegetation are addressed in this paper and have been the subject of numerous related publications (e.g., Tyler et al., 1996; Walvoord et al., 2002a, 2002b, 2004).
|
SITE DESCRIPTION
|
|---|
The Area 3 Radioactive Waste Management Site (RWMS) encompasses about 0.52 km2 (
128 acres) and originally included eight subsidence craters created from underground nuclear tests. Two craters were combined into one disposal cell (U-3ax/bl) that is now completely filled, capped, and closed (Fig. 1
). Two other craters were also combined into one disposal cell (U-3ah/at) that is still being filled as of today (2006). Of the remaining craters, U-3bh is receiving waste and Craters U-3az and U-3bg are reserved for future use (Bechtel Nevada, 1998).
View larger version (52K):
[in this window]
[in a new window]
|
Fig. 1. Map of the NTS showing the location of the Area 3 RWMS with waste-disposal areas and selected boreholes.
|
|
Meteorological data collected during a 36-yr period from a station approximately 5 km to the northwest of the Area 3 RWMS at a 1241-m elevation indicate an average annual precipitation of 161 mm (Bechtel Nevada, 1998), with considerable year-to-year variability. Because of its similar elevation (1225 m) and close proximity to the weather station, the Area 3 RWMS is assumed to have the same average precipitation. Lysimeter studies conducted in Frenchman Flat about 22.4 km south of the Area 3 RWMS indicate that most or all of the soil moisture that results from the infiltration of winter precipitation and occasional summer thunderstorms is removed during the summer months when evapotranspiration (ET) rates are large because of high air temperatures and plant activities (Bechtel Nevada, 1998).
The unsaturated zone at the Area 3 RWMS is 492 m thick and consists of approximately 305 to 460 m of alluvium underlain by tuff aquifers and confining units (Bechtel Nevada, 1998, p. 5). Seven boreholes ranging in depths from 50 to 147 m were drilled beneath the Area 3 RWMS using air as the drilling fluid. Core samples taken at approximately 3-m intervals in these boreholes provide the subsurface hydrologic and geochemical data used here. To avoid disturbing the waste in the open or retired disposal cells, boreholes, except those in the U-3bh cell and U-3bl-U1, were drilled at an angle beneath the waste from undisturbed areas adjacent to the disposal cells (Fig. 1). Borehole U-3bl-U1 was drilled at an angle away from the disposal cell to sample undisturbed sediments. All of the angled boreholes originate in undisturbed material covered by sparse xeric vegetation, and those intersecting disturbed zones do so only at depth. Therefore, data from the upper parts of these angled boreholes can be used in models of the natural hydrologic system.
To supplement the Area 3 data, this study also relies on the selective use of unsaturated-zone data from Borehole UE-6s in Area 6 of the NTS (Fig. 1). Data from this 45-m-deep borehole were analyzed by Walvoord et al. (2002b) in conjunction with data from the nearby 1200-m-deep Borehole UE-6e. The UE-6e borehole indicates that the unsaturated zone is approximately 464 m thick, with the alluvium–tuff contact located at a depth of 409 m (Walvoord et al., 2002b).
|
CONCEPTUAL MODEL
|
|---|
The basic conceptual model governing the development and implementation of the numerical models discussed here is described in Walvoord et al. (2002a, 2002b) and Wolfsberg and Stauffer (2003). In brief, this conceptual model maintains that in thick vadose zones found in alluvial basins of the arid southwest, the present-day distributions of moisture and of dissolved solutes such as Cl reflect the influences of past wetter climates, as well as the slow, ongoing approach to equilibrium with drier conditions established in the soil zone during the transition to the present-day climate. The change from pluvial to postpluvial time periods occurred approximately 10 000 yr ago, roughly coincident with the transition between the Pleistocene and Holocene epochs.
The conceptual model presented in Walvoord et al. (2002a) proposes that, although increasing aridity at the time of the pluvial to postpluvial climate transition affected soil moisture, the most important factor governing soil zone moisture is vegetation type. As climate became increasingly arid, vegetation changed from a juniper woodlands (mesic) type to a desert shrub (xeric) type. Although climate change was probably synchronous over large regions, packrat midden data indicate the vegetation response was elevation dependent, occurring earlier at low elevation sites (13 ka at 900-m elevation) and later at higher elevations sites (9.5 ka at 1200-m elevation) (Spaulding, 1990). The new desert shrub type vegetation was more efficient at scavenging moisture and created low water potentials (large capillary suctions) in the root zone. These low water potentials prevent the downward movement of infiltration past the root zone and create upward hydraulic and vapor-pressure gradients that cause upward liquid and water-vapor fluxes for some distance below the root zone. (Vapor-pressure gradients toward the root-zone are created because vapor pressures are less where water is bound to grain surfaces by large capillary forces.) At greater depths, below the so-called "zero-flux" plane, the drainage of water that infiltrated during earlier pluvial periods continues, eventually recharging the saturated zone. Even below the zero-flux plane, however, vapor fluxes are upward because of vapor-pressure gradients created by the geothermal gradient.
The effects of past and present climates and plant communities on moisture boundary conditions are manifested in depth profiles of pore-water Cl concentrations and water potentials. The Cl dissolved in vadose zone pore waters originates from salts dissolved in precipitation and from windblown dust. The concentration of Cl in pore water is a function of the total deposition rate and the subsequent increases in dissolved Cl concentrations as water (but not Cl) is removed by ET processes in the soil zone. Pore-water data from desert environments, including those from the NTS, often show very high Cl concentrations near the root zone and low Cl concentrations at greater depths. Under the assumption that total Cl deposition rates have been relatively constant with time, these observations have been interpreted to reflect the timing and duration of major transitions in climate and plant communities. Water potential profiles characterized by extremely low values in the root zone (–6 to –8 MPa) and asymptotically increasing water potentials at greater depths have likewise been interpreted to reflect climate transitions, with secondary effects superimposed by vapor-transport processes affecting the curvature of these profiles (Walvoord et al., 2002a). Specifically, a multiphase, nonisothermal model is needed to simulate this system because gravity drainage dominated moisture redistribution immediately after the climate change, but vapor-phase transport became increasingly important as the system continued to dry (Walvoord et al., 2002a).
|
DATA SUMMARY
|
|---|
Data used to evaluate long-term moisture and chemical transport processes at the Area 3 RWMS originate from Boreholes U-3at-D1, U-3bl-U1, and U-3bl-D1 (Fig. 1). Borehole U-3bl-U1 was drilled at an angle beneath an undisturbed part of Area 3 RWMS; the other two boreholes originate in undisturbed material and eventually enter the chimney beneath active and closed waste-disposal trenches at depths of about 30 m. (Note that all borehole depths have been converted to depths below land surface in the figures and discussions presented here.) Boreholes shown in Fig. 1 but not included in this study include (i) those drilled into the bottom of the U-3bh crater, because they clearly indicate the presence of recent infiltration as a result of focused surface runoff in the crater, and (ii) U-3at-D2 and U-3bl-D2, because no shallow data from the undisturbed zone were collected in these boreholes.
Water-Potential Data
Water-potential data were collected from the upper 50 m of the three boreholes (Fig. 2
). The water-potential profiles show the low root-zone values, concave-downward curvature, and asymptotic approach to wetter conditions with depth that are characteristic of profiles measured elsewhere at the NTS in vertical boreholes (Walvoord et al., 2002a, 2002b). Although data from Borehole U-3at-D1 show considerably more scatter than data from the other boreholes, the data as a group are remarkably consistent. Water potentials from the U-3at-D1 and U-3bl-D1 boreholes do not show substantial differences from U-3bl-U1 data, where the U-3at-D1 and U-3bl-D1 boreholes penetrate the subsidence craters at a depth of about 30 m. The consistency of the water-potential data from the three boreholes suggests that data from within the U-3ax/bl and U-3ah/at craters have not been influenced by runoff generated within the craters, a process that has increased recent infiltration rates above long-term values in the U-3bh crater (Kwicklis et al., 2006). Unlike the vertical boreholes in the U-3bh crater, Boreholes U-3at-D1 and U-3bl-D1 are located along the periphery of the craters, rather than at low areas in the center of the crater where surface runoff would tend to pond.
Pore-Water Isotope Data
Pore-water isotope data from Area 3 include deuterium (2H or D) and oxygen-18 (18O). The relative abundances of these isotopes compared with the more abundant isotopes 1H and 16O are typically expressed in delta () notation as
 | [1] |
where "isotope" is either D or 18O, Rsample is the mass ratio of the heavier and lighter isotope (D/1H or 18O/16O) in the sample, and Rstandard is the ratio of these isotopes in the standard (Vienna Standard Mean Ocean Water or VSMOW) (Clark and Fritz, 1997). The units of D and 18O are parts per thousand or permil. The exact values of D and 18O in precipitation depend on many factors, including the temperature and isotopic composition of the source water (typically the ocean), the distance and topographic features between the moisture source and the point of precipitation, and the condensation temperature of the precipitation, which varies with season, elevation, and as a result of climate change (Clark and Fritz, 1997; Ingraham et al., 1991; Friedman et al., 2002a, 2002b). The D and 18O of precipitation worldwide generally follows the Global Meteoric Water Line (GMWL) (Clark and Fritz, 1997):
 | [2] |
The intercept (10) in this equation is known as the deuterium excess and is a function of the relative humidity over the moisture source area (Merlivat and Jouzel, 1979). Groundwater data suggest that the deuterium excess may have been as small as 5 during the late Pleistocene, when relative humidity over the oceans may have been higher because of globally cooler conditions (Benson and Klieforth, 1989).
The isotopic composition of precipitation can be modified by partial evaporation at the ground surface and in the shallow subsurface, so that soil water and recharge can have isotopic values that plot below the GMWL along a line with slope that can vary between 3 and 5, depending on relative humidity (Clark and Fritz, 1997). More evaporated water plots farther to the right of the GMWL along the evaporation trend than does less evaporated water. The original isotopic composition of the recharge before evaporation can be estimated by projecting the measured isotope compositions back along the evaporation trend to the GMWL or a similar line that describes the isotopic composition of local precipitation at the time of recharge (if this line is known or can be estimated).
Studies of the isotopic composition of modern precipitation and groundwater in southern Nevada and elsewhere in the Great Basin have noted that groundwater often has much lighter (more negative) values of D and 18O than modern precipitation, an observation that has led to the understanding that much of the groundwater in southern Nevada is fossil groundwater recharged during a period when the climate was cooler and wetter than today (Benson and Klieforth, 1989; Smith et al., 2002). Records of paleoclimatic change in southern Nevada preserved in the 18O of calcite precipitated at a regional discharge area south of the NTS suggest that recharge during glacial maxima has a 18O that is about 2
lighter than recharge during interglacial periods (Winograd et al., 1992), confirming the basic premise that fossil groundwater that was recharged during glacial (or pluvial) periods is isotopically distinct from modern precipitation.
Pore-water isotopic data from the inclined boreholes in Area 3 show features that can be interpreted in the context of the discussion presented in the preceding paragraphs. Systematic decreases in isotopic values are evident with increasing depth, with pore-water 18O values decreasing from near-surface values of about –9 to approximately –15 at a 50-m depth (Fig. 3
) and pore-water D values decreasing from near-surface values of approximately –92 to about –113 at a 50-m depth (Fig. 4
). A scatter plot of D versus 18O values from the three boreholes, along with the present-day GMWL and a possible paleometeoric water line for southern Nevada (White and Chuma, 1987), suggest that pore waters with D values greater than –108 or 18O values greater than –14 have been influenced by evaporation (Fig. 5
). Although pore-water data near the paleometeoric water line show some scatter that may reflect earlier, short-duration climate fluctuations, the trend of the evaporation line suggests that pre-evaporated waters may have been derived from precipitation with initial D and 18O values of about –110 and –14.5, respectively. As mentioned above, the D and 18O of precipitation in southern Nevada during the late Pleistocene is believed to have been lighter than at present because of cooler temperatures that existed at that time and differences in the paths of storms bringing moisture to southern Nevada (Benson and Klieforth, 1989).
Although D and 18O variations with depth are not specifically addressed in the conceptual model described, the observed trends in the pore-water isotopic values have interpretations that are consistent with that model. Namely, precipitation with D and 18O values considerably lighter than present-day values infiltrated with relatively little enrichment by evaporation during the pluvial period. As climate became increasingly arid and warmer, evaporation of the pore water in the shallow subsurface became more important. Downward liquid and vapor diffusion of isotopically enriched water, combined with deeper evaporation and in situ enrichment of pore-water below the root zone once the xeric vegetation was in place, caused the shallower part of the D and 18O profiles to undergo isotopic shifts, even in the absence of deep infiltration beneath the root zone.
Water Content and Saturation Data
Paired gravimetric moisture content and dry bulk-density data were available from Boreholes U-3bl-U1 and U-3bl-D1, allowing the conversion from gravimetric to volumetric water content at these boreholes. The resulting volumetric water-content values show considerable scatter and do not indicate obvious trends with depth (Fig. 6a
). The average moisture content for the combined data from Boreholes U-3bl-U1 and U-3bl-D1 is 0.108 ± 0.022 (1 
). Water saturation data display similar scatter within the upper 50 m (Fig. 6b).
View larger version (12K):
[in this window]
[in a new window]
|
Fig. 6. Comparison of (a) water content and (b) saturation data from the Area 3 RWMS with model results.
|
|
Pore-Water Chloride Concentration Data
Pore-water Cl concentration data were obtained from Boreholes U-3bl-U1 and U-3bl-D1 (Fig. 7
). These data suggest a shallow Cl bulge typical of desert alluvial basins exists in the root zone above about an 8-m depth in Area 3. However, the sampling frequency was insufficient to fully define the bulge or the mass of Cl contained within it. The mass of Cl contained within the bulge is the product of the Cl deposition rate and the elapsed time since the shift to xeric vegetation. If one of these factors can be determined through independent means, the other can be calculated, provided the total mass of Cl in the root zone is well characterized.
Hydrologic and Tracer Data from Area 6 Borehole UE-6s
To compensate for the incomplete Cl inventory provided by the Area 3 data, this study relied on additional data from Area 6 (Fig. 8
) reported in Hartsough et al. (2001). Most of the Cl profiles from Area 6 show peak Cl concentrations of about 600 to 3700 mg L–1 at a depth of about 2 m. The Cl concentration profile from Borehole UE-6s is fairly representative of the Area 6 data and was included in the present analysis to help constrain estimates of Cl deposition rates and the timing of the climate transition at Area 3. Hydrologic and natural tracer data from the UE-6s borehole in Area 6 of the NTS were previously analyzed by Walvoord et al. (2002b). The water potential data from Borehole UE-6s show many of the same features as the Area 3 water potential profile (Fig. 2), but have a slightly less negative water potential (–4.4 MPa) in the root zone (Fig. 9
). The pore-water D and 18O profiles from Borehole UE-6s (Fig. 10
and 11
) have a more complex shape than the isotope data from Area 3 (Fig. 4 and 3), but nonetheless appear to be derived by evaporation from a source water with a D of about –116 and a 18O value of about –15.2 (Fig. 12
). These observations suggest that conditions in the vicinity of Borehole UE-6s may have undergone periodic changes that caused evaporation rates to vary. One possible explanation, although not the only one, is that periodic flooding of a nearby playa lake locally raised the relative humidity near Borehole UE-6s, thereby temporarily suppressing evaporation rates.
|
MODELING APPROACH
|
|---|
The flow and transport models described in this paper were created with the computer code FEHM (Zyvoloski et al., 1997; Tseng and Zyvoloski, 2000), which is able to simulate the simultaneous movement of liquid, vapor, and heat, and consider multiphase, multicomponent transport. Numerical models of water, heat, and solute transport were created for Areas 3 and 6 for both the pluvial and postpluvial periods associated with the end of the Pleistocene and the Holocene, respectively. For each area, the model results at the end of the pluvial period formed the initial conditions to the model for the postpluvial period to represent the Pleistocene–Holocene transition. The consistent values and trends in the water potential and tracer concentration data from three boreholes in Area 3 suggest that flow and transport processes are dominantly one-dimensional and relatively uniform throughout Area 3. Hence, a single, one-dimensional model can adequately simulate the conditions at the three inclined boreholes in Area 3. Similar uniformity and dimensionality were assumed for the unsaturated zone in Area 6 based on the observations of Hartsough et al. (2001). Each area is modeled with a 226-node, one-dimensional finite-element grid with fine resolution (0.1 m) near the upper boundary, where large water-potential, water-vapor pressure, and chemical gradients are expected to develop, and coarser resolution elsewhere (maximum 4 m). The entire unsaturated-zone thickness of approximately 500 m is represented by hydraulic properties typical of undisturbed alluvium at these sites (Table 1).
View this table:
[in this window]
[in a new window]
|
Table 1. Base-case hydrologic properties used in flow and transport simulations at the Area 3 RWMS and Borehole UE-6s.
|
|
The computer code PEST (Doherty, 2002) was coupled to FEHM to estimate selected parameter values and their uncertainties and to derive statistical relations among the parameters. PEST obtains these quantities using the Levenberg–Marquardt algorithm for adjusting uncertain parameters in a systematic way until an optimal match between the FEHM model results and the observations are obtained. Optimal parameters are considered to be those that result in the smallest value for the objective function, which is calculated as the sum of the weighted-residuals squared, 
(
r)2.
A complete description of the models and methods contained in FEHM is beyond the scope of this paper and is described elsewhere (Zyvoloski et al., 1997; Tseng and Zyvoloski, 2000). However, some aspects of the modeling presented in this paper are relatively innovative, and these are described in some detail in the following subsections.
Phase Partitioning and Fractionation of Stable Isotopes
To simulate D and 18O transport, the isotopic water species, H216O, H218O, and HD16O, are treated as separate volatile species with different mass-based diffusion coefficients and different phase partitioning according to Henry's Law. By considering the differences among these species, both kinetic and equilibrium fractionation are accommodated by the model. H218O and HD16O constitute only trace fractions of the total number of water molecules in any sample. In natural waters, the proportions of H216O to H218O and H216O to HD16O are 1000:2 and 1000:0.32, respectively (Merlivat and Coantic, 1975). Thus, as in the reactive geochemical transport model PHREEQC (Thorstenson and Parkhurst, 2002; Parkhurst and Appelo, 1999), the equation structure in FEHM assumes the presence of only one solvent, liquid H216O. To simulate D and 18O transport, the isotopic water species, H218O and HD16O, are treated as electrically neutral, volatile solutes.
The concentration of water molecules in the gas phase is approximated using Henry's Law equilibrium partitioning coefficients for the dominant water molecule H216O (Jury et al., 1991). The molecules with the heavier trace isotopes, H218O and HD16O, do not partition between the liquid and gas phases with the same equilibrium ratios as the dominant H216O. Their concentrations in the vapor phase are estimated by converting equilibrium fractionation factors into Henry's Law coefficients, KH218OgH and <KHD16OgH, similar to Thorstenson and Parkhurst (2002), as follows:
 | [3] |
 | [4] |
where KH216OgH is the Henry's Law constant for water and 
i are the equilibrium fractionation factors for 18O and D. The equilibrium fractionation factors can be defined as ratios of saturation vapor pressures (Friedman and O'Neil, 1977; Merlivat and Coantic, 1975) as follows:
 | [5] |
 | [6] |
Henry's Law constants in FEHM are specified in units of MPa, to solve the equation
 | [7] |
where Pv is the vapor pressure of either H218O or HD16O, X is the mole fraction in water, and KH is the Henrys Law constant. The Henry's Law constants used in this study are listed in Table 2.
Kinetic fractionation is captured in the model as a result of different vapor diffusivities for H216O, H218O, and HD16O. Here, we use the measured ratios of the vapor diffusion coefficients provided in Merlivat (1978) as listed in Table 2. Although small kinetic fractionation processes due to evaporation and diffusion in the vadose zone are captured with this model, the dominant effect associated with the stable isotope profiles is due to the enriched ratios resulting from substantial evaporation in the shallow vadose zone. The effects of shallow evaporation on pore-water isotopic values are represented in the model by specifying isotopic compositions for the infiltrating water at the upper boundary that are similar to those measured in pore water from the root zone.
Liquid and Vapor Diffusion Models
The diffusive fluxes of liquid- and vapor-phase components are governed by concentration gradients as described by Fick's Law. In porous media, the grains cause solutes to take tortuous pathways and reduce the cross-sectional area across which diffusion can take place. Therefore, effective diffusion coefficients are calculated which reduce the mass flux by scaling the free-air or free-liquid diffusion coefficients for individual solutes. Although FEHM includes several empirical models that calculate effective liquid and vapor diffusion coefficients from moisture content and porosity (e.g., the diffusion model of Millington, 1959), none of the liquid diffusion models and only one of the vapor-diffusion models contain adjustable parameters that permit them to be optimized to fit the Area 3 data. To permit flexibility in adjusting effective liquid diffusion coefficients and yet limit the number of adjustable parameters, all liquid chemical (Cl) and isotopic species (H216O, H218O, and HD16O) were assumed to have the same constant (but adjustable) effective liquid diffusion coefficients. The absence of pronounced trends in the Area 3 moisture content data (Fig. 6) and the relatively uniform moisture contents generated by the models for Areas 3 and 6 (see Model Results) suggest that the use of variable moisture-content diffusion models is not necessary for liquid-phase solutes in this system. Likewise, because free-air diffusion coefficients are about four orders of magnitude larger than free-water diffusion coefficients for the isotopic species that are partitioned by Henry's Law, approximations of the liquid diffusion coefficients can be made with negligible effect.
For the vapor-phase species H216O, H218O, and HD16O, the effective vapor diffusion coefficients are calculated from the relation (Vargaftik, 1975)
 | [8] |
where Dva is the effective diffusion coefficient of vapor in air (m2 s–1), 
is the tortuosity (the linear displacement divided by the actual flow path length) (L L–1), 
is the porosity (L3 L–3), Sg is the gas saturation, Dva0 is the free-air vapor diffusion coefficient at standard pressure and temperature, P is the ambient total gas pressure (MPa), P0 is standard pressure (0.1 MPa), T is temperature (°C), and 
is an empirically determined dimensionless factor (1.81) that accounts for the increase in the kinetic energy of water-vapor molecules at higher temperatures. The values of Dva0 differ depending on the isotopic species undergoing transport (Table 2).
Boundary Conditions
For both the Area 3 and Area 6 models, the pressures on bottom boundary nodes were fixed at atmospheric pressure to represent water-table conditions. Total gas pressures throughout the unsaturated zone were chosen so that the gas column remains approximately static; however, volume changes associated with the evaporation and condensation of vapor induce a very small advective gas flux. During the postpluvial period, changes in the air-filled pore volume associated with the drainage of paleoinfiltration also induce a small amount of gas flow. Nonetheless, the movement of vapor is dominated by diffusion associated with the vapor pressure differences caused by the geothermal gradient and, during the postpluvial period, by vapor-pressure lowering that results from large capillary forces in the root zone. Fixed temperatures of 15 and 24°C are specified at the top and bottom boundaries of the model to establish an average temperature gradient of 0.018°C m–1 in the unsaturated zone for both the pluvial and postpluvial periods. This gradient, combined with assumed thermal conductivity values of 1.5 J s–1 m–1 °C–1, resulted in an average heat flux of 27.0 x 10–3 J s–1 m–2.
The moisture and tracer boundary conditions at the top of the model are specified differently for pluvial and postpluvial periods (Table 3). During the pluvial period, the infiltration rate and tracer concentrations represent the deep infiltration that penetrates below the root zones of mesic vegetation. Paleoinfiltration rates were estimated from deep pore-water Cl concentrations, assuming that Cl deposition rates during both the pluvial and postpluvial periods were the same as those estimated by Tyler et al. (1996) for the present-day in Frenchman Flat (105 mg m–2 yr–1). These assumptions result in estimated paleoinfiltration rates of 7.6 mm yr–1 in Area 3 and 2.3 mm yr–1 in Area 6 (Table 3). The Cl concentrations and isotopic values of the deep infiltration were determined from the measured Cl concentrations and D and 18O values of the pore water far beneath the Cl bulge. Pluvial conditions are assumed to have been in effect before an abrupt transition to the present-day arid climate. The duration for the pluvial period was set sufficiently long to flush the system of prepluvial water, none of which is indicated in the Area 3 or Area 6 data sets, unlike Area 5, where boreholes show a secondary deep Cl bulge indicative of incomplete pluvial flushing (Tyler et al., 1996; Walvoord et al., 2002b; Wolfsberg and Stauffer, 2003). Because it is assumed that flushing of prepluvial pore water was complete at the end of the pluvial period in Areas 3 and 6, pore water throughout the unsaturated zone has the chemical and isotope characteristics listed for the pluvial period in Table 3 at the beginning of the postpluvial period of the simulation.
View this table:
[in this window]
[in a new window]
|
Table 3. Summary of boundary conditions used for the Area 3 and Area 6 simulations during pluvial and postpluvial periods.
|
|
Initial conditions for the postpluvial period were determined from the simulation results at the end of the pluvial periods in Areas 3 and 6. As discussed in the Conceptual Model section, the xeric vegetation that replaced the mesic vegetation at the end of the last pluvial was very efficient at removing moisture from the soil, and extremely low water potentials developed in the root zone. The low water potentials in the root zone prevent the downward movement of water past the roots and create upward hydraulic and vapor-pressure gradients that draw moisture back up toward the roots in both the liquid and gas phases. Transpiration processes are represented in the models through the specification of low water potentials in the root zone at the beginning of the postpluvial period. These low water potentials cause water to flow downward from land surface and upward from depth toward the root-zone nodes, where it is removed from the model. Based on measured water-potential profiles (Bechtel Nevada, 1998; Walvoord et al., 2002b), a root-zone water potential of –7.6 MPa was imposed at a 2.0-m depth in the Area 3 model and a root-zone water potential of –4.7 MPa was specified at a 2.75-m depth in the Area 6 model. These continuous-in-time water potential boundary conditions are supported by Andraski's (1997) 5-yr record for a similar system in which temporal water potential variations are damped out at the bottom of the root zone where near-constant water potentials are preserved. Water-vapor diffusing upward toward the root-zone node condenses because of the low equilibrium vapor pressures associated with large capillary forces, and the condensate exits the model through the root-zone node as liquid. In accordance with the view that transpiration is a nonfractionating process (Clark and Fritz, 1997), the various isotopes of water (H218O and HD16O) exit the model through the root-zone node along with the bulk of the water (H216O). The fractionating effects of shallow evaporation evident in the scatter plots of D versus 18O (Fig. 5 and 12) were considered by specifying values of D and 18O in infiltration that already account for the effects of this process. The infiltration applied to the top of the model during the postpluvial period was set equal to the estimated precipitation rate of 161 mm yr–1 at Area 3 (Bechtel Nevada, 1998, Fig. 1
–6) and 167 mm yr–1 at Area 6 (Walvoord et al. (2002b, p. 27-2). Because all of the infiltration leaves the model at the root-zone node, the main purpose of applying infiltration is to introduce Cl into the soil zone and set the upper boundary conditions for D and 18O compositions. Downward liquid diffusion provides for the spread of the Cl bulge to depths of 6 to 10 m in Fig. 7 and 8. Gas diffusion of H218O and HD16O molecules permits the much deeper penetration of the 18O and D fronts following climate change (Fig. 3 and 4). The Cl arriving at the ground surface as wet and dry fallout is carried downward toward the root-zone nodes with the infiltration, but unlike H218O and HD16O, the dissolved Cl is assumed to be left behind in the root zone, where with time it accumulates to form the Cl bulge characteristic of the desert environments.
Precalibration Scoping Simulations with Base-Case Parameters
A numerical model for the Area 3 RWMS was run with the base-case hydrologic parameters given in Table 1, an assumed Cl deposition rate of 105 mg m–2 yr–1 for both the pluvial and postpluvial period, and the empirical Millington diffusion model (Millington, 1959) for liquids and gases encoded in FEHM (Wolfsberg and Stauffer, 2003). Results of these simulations for a postpluvial period of up to 20000 yr duration indicate that this combination of parameters and diffusion models produces results that are internally inconsistent in terms of the probable duration of the postpluvial period (Kwicklis et al., 2006). The observation that numerical models of hydrologic data and various natural tracer data from the unsaturated zone appear to point to different durations for the postpluvial period was also noted by Walvoord et al. (2004), who suggested that drying and vegetation shifts may have preceded the major warming trend. It also is possible that the apparent discrepancies in the timing of the climate transition indicated by the models result from uncertainties in the embedded submodels and their parameters (e.g., the van Genuchten effective hydraulic conductivity model or the Millington diffusion model), and in other input parameters like Cl deposition rates. In this study, it was assumed that although the climate transition was gradual, the change from mesic to xeric vegetation was relatively abrupt at any particular location, and that this transition exerts the dominant effect on soil moisture conditions. The problem of reconciling the apparent discrepancies in the postpluvial duration is therefore considered as a problem of parameter identification in this paper.
Calibration Targets and Optimization of Selected Model Parameters with PEST
Data from the Area 3 boreholes used as calibration targets with PEST include the U-3bl-U1 water potential and saturation profiles (Fig. 2 and 6), the pore-water Cl concentration profiles from Boreholes U-3bl-U1 and U-3bl-D1 (Fig. 7), and the D and 18O values estimated at approximately 2-m intervals to a 50-m depth by fitting a sixth-order polynomial to the combined 18O and D from Boreholes U-3bl-U1, U-3bl-D1, and U-3at-D1 (Fig. 3 and 4; Kwicklis et al., 2006). The different treatment of the water-potential and isotopic data from the Area 3 boreholes arose because, whereas the water potential data from U-3bl-U1 can be considered representative of the water-potential profiles from the other two Area 3 boreholes (Fig. 2), the isotope data show considerable scatter at a given depth, and it is difficult to identify a simple set of target observations in the absence of curve fitting.
Preliminary attempts to calibrate the model using only data from Area 3 indicated that the lack of detail in the upper parts of the Cl profile prevent the reliable estimation of parameters related to the time of climate change (duration of the postpluvial period, thol) and Cl deposition rates (Cl concentration of postpluvial precipitation, Clin). To compensate for the lack of detail in Cl measurements at Area 3, more detailed Cl profiles from Borehole UE-6s in southern Yucca Flat were used in the optimization processes. The use of Cl data from borehole UE-6s to estimate parameters at Area 3 is justified on the basis of the similar hydrologic settings and elevations of the UE-6s and Area 3 boreholes. Although water potential, saturation, and isotope data from Borehole UE-6s are not used as calibration targets in the PEST optimization, a comparison between these data and simulation results was made to ensure that reasonable predictions of these quantities are obtained at Borehole UE-6s.
If no water leaves the root zone except through evapotranspiration, the mass of Cl (mg m–2) left behind in the root-zone Cl bulge can be expressed as
 | [9] |
where P is the precipitation rate (L/yr·m2), Clin is the Cl mass arriving at land surface as wet and dry fallout per liter of precipitation (mg L–1), and thol is the duration of the postpluvial period or, more precisely, the time since modern vegetation was established (yr). Thus, if either Clin or thol is known, the other term can be estimated can be estimated from the mass of Cl in the root zone, but they cannot both be estimated independently from the mass of root-zone Cl alone. In this paper, it is assumed that at the 1200-m elevation of the study area, the Cl bulge in the root zone began accumulating only after xeric vegetation was established approximately 9500 yr ago (Walvoord et al., 2002a, 2002b). However, rather than treat thol as a completely certain parameter, it is implemented as an adjustable parameter in PEST with a small range of allowable values (thol = 9500 ± 500 yr) so that its uncertainty and its effects on other adjustable parameters can be calculated. The potential impacts of constraining thol in this manner are evaluated in the Discussion section.
In addition to Clin and thol, other adjustable parameters considered during parameter estimation include the effective liquid diffusion coefficient (Deff), the saturated permeability of the alluvium (ksat), and the value of the tortuosity () used in Eq. [8] for calculating the effective vapor diffusion coefficient (tort). Although ksat might seem to be one of the more certain parameters, it is used in the van Genuchten (1980) model to estimate effective hydraulic conductivity at conditions far from saturation. Therefore, it was adjusted to compensate for the uncertainty in the van Genuchten (1980) effective hydraulic conductivity function at low water potentials and saturations. Initial scoping calculations indicated that increases in effective hydraulic conductivity at low moisture contents and water potentials were necessary to match the water-potential data. These increases were most readily accomplished through increases in ksat. A good match to the water-potential data was considered essential for estimating both liquid and vapor fluxes and reproducing the depth of the zero-flux plane indicated by the water-potential profile.
In summary, five parameters (ksat, thol, Clin, Deff, and tort) were estimated during the calibration process by simultaneously fitting the water potential, moisture content, Cl, and stable isotope data. In recognition that some pairs of parameters such as thol and Clin were likely to be highly correlated, some constraints were imposed on the possible duration of the postpluvial period (thol) using information from packrat midden studies. Other simplying assumptions were made to reduce the number of parameters to be estimated during calibration. These include the assumption that a single diffusion coefficient (Deff) is adequate to model the diffusion of all liquid-phase solutes and isotopic tracers. Additionally, it was assumed that gas-phase tortuosity (tort) is the most uncertain parameter in the equation governing gas-phase diffusion and that uncertainty in the effective permeability at low saturations could be best addressed by adjusting the saturated permeability (ksat) rather than by making adjustments to parameters in the moisture-retention function.
|
RESULTS OF MODEL CALIBRATION
|
|---|
Parameter Values and Uncertainties
The estimated parameter values and their lower and upper 95% confidence limits are given in Table 4. The value of ksat estimated by PEST is about two orders of magnitude greater than the measured permeability values at Area 3 (Table 1). Again, it is unlikely that ksat is actually this large; the parameter ksat is used in this model to increase the effective hydraulic conductivities predicted by the van Genuchten function at low water potentials (and saturations), thereby allowing the simulated water-potential profile to respond more rapidly to the low root-zone water potentials imposed at the beginning of the postpluvial period. The estimated Cl concentration of precipitation in the postpluvial period is 1.03 x 10–5 mol kg–1 (0.364 mg L–1), which combined with an estimated precipitation rate of 161 mm yr–1, yields Cl deposition rate of 59 mg m–2 yr–1 (with lower and upper limits of 36 and 95 mg m–2 yr–1). These rates are lower than the value of 105 mg m–2 yr–1 assumed by Tyler et al. (1996) for Area 5 in Frenchman Flat. However, the estimated Cl concentration is nearly identical to the average Cl concentration of precipitation (0.35 mg L–1) measured over a 5-yr period in the Kawich Range about 120 km northwest of our study area (Meijer, 2002). The duration of the postpluvial period (thol) was estimated to be at its permissible upper bound of 10 000 yr, a reasonable value in light of the packrat middens studies cited previously. However, to account for the effects of uncertainty in thol on other estimated parameters, parameter sensitivities and uncertainties were calculated with thol unconstrained by these parameter bounds. In this case, the upper and lower 95% confidence limits on thol were estimated to be between 6443 and 15522 yr when the best estimate of thol was fixed at 10 000 yr. The value for Deff of 4.54 x 10–12 m2 s–1 is somewhat smaller than the value for Deff of 7.2 x 10–12 m2 s–1 estimated from the Millington diffusion function (Millington, 1959) using a moisture content of 0.10, a porosity of 0.382, and a free-water diffusion coefficient (D0) for H216O of 2.26 x 10–9 m2 s–1. However, unlike the other estimated parameters, the lower and upper uncertainty limits for Deff are extremely large (Table 4), most likely because gas-phase diffusion dominates the migration of HD16O and H218O, and the lower bound on the Cl bulge results from the competing effects of downward Cl diffusion and upward advection. The lack of sensitivity of Deff could also be due to (i) the limited amount and resolution of pore-water Cl data in Area 3 (Fig. 7) and (ii) the high frequency variations in the data in UE-6s (Fig. 13
).
Parameter Sensitivities
To varying degrees, the observations used as calibration targets provide information about the values of the estimated parameters. One measure of the importance of the observations or subset of observations for estimating a parameter is provided by the parameter's sensitivity with respect to those observations. The sensitivity of each parameter to all observations used in the calibration is provided in Table 5, and the sensitivities of each parameter to individual groups of observations (e.g., water potential, pore-water Cl concentration) are given in Table 6. (Note that for log-transformed parameters, relative sensitivities are obtained by multiplying the sensitivity by the absolute value of the log-transformed parameter. Relative sensitivities are a measure of the composite changes in model outputs that result from a fractional change in the value of the parameter.) Overall, the duration of the postpluvial period (thol) has the highest sensitivity and relative sensitivity, whereas the effective liquid diffusion coefficient (Deff) is least sensitive to the observations (Table 5). The low sensitivity of the model to Deff helps to explain the absence of meaningful confidence limits for Deff (Table 4).
Individual parameters have variable sensitivities to different types of observations (Table 6). The value of ksat is influenced most strongly by observations of water potential and pore-water isotope values. The value of Clin depends almost entirely on the pore-water Cl concentrations, whereas the value of thol is influenced strongly by pore-water isotope and water-potential values, with smaller influences from Cl concentration data. Not surprisingly, the value of Deff is influenced most strongly by pore-water tracer concentrations and the value of tort depends on pore-water isotopic and water potential data. The pore-water saturation data do not strongly influence the estimates of any of the parameters. This reflects, in part, the low weights assigned to these observations, as well as the scatter in the saturation data (Fig. 6).
Parameter Correlations
The parameter correlation coefficient matrix (Table 7) identifies the interdependencies among the different model parameters. As Table 7 indicates, strong correlations exist between ksat, Clin, thol and tort. The strong negative correlations (–0.90 or lower) between thol and ksat and tort indicate that, as the time available for liquid and vapor phase moisture transport increases, the rates at which transport occurs become lower. The model can match the water potential and stable isotope profiles only by making opposing changes in this pair of parameters. Likewise, as the time available for Cl accumulation in
TOP
ABSTRACT
INTRODUCTION
