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 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 Google Scholar
Google Scholar
Right arrow Articles by Glassley, W. E.
Right arrow Articles by Grant, C. W.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Glassley, W. E.
Right arrow Articles by Grant, C. W.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Glassley, W. E.
Right arrow Articles by Grant, C. W.
Related Collections
Right arrow Isotopes
Right arrow Software
Right arrow Vadose Zone Processes and Chemical Transport
Vadose Zone Journal 1:3-13 (2002)
© 2002 Soil Science Society of America

REVIEWS & ANALYSES

The Impact of Climate Change on the Chemical Composition of Deep Vadose Zone Waters

William E. Glassley*, John J. Nitao and Charles W. Grant

Lawrence Livermore National Laboratory, Livermore CA 94550
* Corresponding author (glassley1{at}llnl.gov)

Received 25 January 2002.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 Background
 SIMULATIONS
 RESULTS
 DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Chloride mass balance, and stable (deuterium and 18O) and radiogenic (3H, 36Cl) isotope studies of deep vadose zone pore waters have generally concluded that variations in moisture flux can account for the observed variations in abundance of these approximately conservative tracers. It can be inferred, on the basis of these observations and interpretations, that a climate change record is preserved in these vadose zone waters. In arid regions where thick (>100 m) vadose zones persist, it has been concluded that this record may extend back more than 100 000 yr. Consideration of the mechanisms that control reactive transport led to the conclusion that such climate-driven effects will also be evident as chemical reactions involving dissolution and/or precipitation of mineral phases along the flow pathway. As a result, there should also be variations in the concentrations of nonconservative chemical species that correspond to changes in the concentrations of the conservative tracers. Simulations of this reactive transport, in a regime typical of the arid U.S. Southwest, demonstrate that these changes can modify pore water chemistry by factors of up to 200%, but the changes take place slowly, requiring thousands of years to achieve steady-state conditions. This suggests that a very rich archive of climate change history is preserved in this type of setting. However, extracting that history is currently hampered by limitations in data and models (e.g., effective mineral reactive surface areas, fluid flow pathways, and quantified models of wetted fracture surface in unsaturated, fractured systems). This challenge may be overcome if coordinated efforts are undertaken that exploit the power of detailed studies of isotope systematics, microscale rock characterization, and high performance computing.


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
Background
 SIMULATIONS
 RESULTS
 DISCUSSION
 CONCLUSIONS
 REFERENCES
 
DURING THE LAST QUARTER CENTURY, research that has examined the distribution of "conservative" tracers in vadose zone pore waters has concluded that this geological setting preserves a record of local moisture flux (Allison et al., 1985; Barnes and Allison, 1988; Cook et al., 1989, 1994; Dettinger, 1989; Scanlon, 1991; Walker et al., 1991). In arid environments, where the vadose zone can be hundreds of meters thick, and infiltration fluxes are low, it has been inferred that this record may extend as much as 100 000 yr into the past (Tyler et al., 1996). If this conclusion is correct, the potential exists to reconstruct moisture flux records on continental land masses in many parts of the world. If moisture flux can be treated as a proxy for changes in precipitation, then the possibility exists that long-term climate records may exist in some vadose zone settings.

The chemical composition of vadose zone pore water is a product of those physical and chemical interactions that take place among water, mineral surfaces, and pore gas in the unsaturated zone (Faybishenko, 2000). If moisture flux varies in the vadose zone as a result of variation in climate, it is likely that the dissolved load in vadose zone pore waters will also vary in sympathetic ways. It is the purpose of this paper to evaluate the magnitude and extent of this potential sympathetic variation using published data and simulations of reactive transport in the vadose zone. The conclusion is reached that significant changes in the composition of vadose zone pore waters should occur in response to climate-forced changes in infiltration flux and surface temperature. As a result, there is a potentially rich archive of regional-scale climate change data preserved within large regions of most continental land masses. However, interpreting empirical pore water chemical compositions in terms of past climate change is currently problematic because of uncertainties in key parameters that control the chemical kinetics and the hydrologic history.


   Background
TOP
 ABSTRACT
 INTRODUCTION
 Background
SIMULATIONS
 RESULTS
 DISCUSSION
 CONCLUSIONS
 REFERENCES
 
That vadose zone pore waters may preserve a nonconservative chemical signature of climate change is implicit in two lines of reasoning. One line of reasoning, which considers time variation of water flux across the land surface boundary layer, suggests that changes in climate-driven moisture flux are preserved in the vadose zone as slowly migrating variations in concentration of the solute load. The second line of reasoning, which considers the kinetics of dissolution and precipitation reactions in the vadose zone, implies that the reactions that are responsible for the evolution of pore water chemical characteristics may be responsive to changes in the land surface temperature and infiltration flux boundary conditions. Both of these considerations are discussed below.

Variability of Moisture Flux Through the Vadose Zone
The analysis of diverse chemical tracers in pore waters has allowed development of conceptual models useful for quantifying the temporal record of water movement rates and residence times in the vadose zone. Chloride was among the first tracers to support such an effort. As the chemical characteristics of Cl-bearing aerosols became better understood (e.g., Winchester and Duce, 1967), and deposition of chloride on land surfaces was conceptualized, use of chloride as a tracer for water movement in vadose zones became relatively common (e.g., Allison et al., 1985; Barnes and Allison, 1988; Bresler, 1973; Cook et al., 1989, 1994; Dettinger, 1989; Ginn and Murphy, 1997; Jolly et al., 1989; Nativ et al., 1995; Peck et al., 1981; Sharma and Hughes, 1985; Tyler and Walker, 1994). As an outgrowth of this work, the first quantitative estimates of the ages of deep vadose zone waters were made. In some studies, the ages of the waters were computed to be tens of thousands to hundreds of thousands of years old (e.g., Allison et al., 1985; Scanlon, 1991; Tyler et al., 1996). Variation in chloride concentration with depth was interpreted to record changes in moisture flux in response to changes in regional climate patterns as well as changes in land use and root zone effects (e.g., Scanlon, 1991; Tyler and Walker, 1994; Walker et al., 1991).

Simultaneously, stable isotope systematics of hydrogen and oxygen were utilized for developing a more rigorous understanding of water movement through the soil zone and into the underlying vadose zone environment (e.g., Zimmerman et al., 1967; Edmunds and Walton, 1980; Barnes and Allison, 1983; Allison et al., 1985; Christmann and Sonntag, 1987). These studies provided additional data regarding water and solute movement, elucidating transport mechanisms that affect both liquid and gas phase chemistries. Complicating interpretation of water ages was clear evidence that chloride deposition rates can be highly variable, even within relatively small regions, and that the mechanisms responsible for chloride transport through the root zone can be complex.

Use of radiogenic isotopes to unravel some of this complexity has provided more detailed information regarding transport processes. Both tritium (3H) and 36Cl were produced in atmospheric nuclear weapons tests. Although analysis of 3H and 36Cl in vadose zone pore waters has not been extensive, insight has been gained from the limited data sets available (e.g., Allison and Hughes, 1978; Cook et al., 1994; Gee and Hillel, 1988; Phillips et al., 1998; Scanlon, 1992; Scanlon and Milly, 1994; Smith et al., 1970; Walker et al., 1991). Striking is the lack of correspondence between ages derived from isotope systems and chloride mass balance calculations (Cook et al., 1994; Gvirtzman et al., 1986; Scanlon, 1991; Van De Pol et al., 1977). In some instances, anion exclusion (Biggar and Nielsen, 1962; Krupp et al., 1972; Van De Pol et al., 1977; Wierenga and van Genuchten, 1989) is inferred to play a role in the computed flux estimates from chloride mass balance (in those cases where chloride appears to migrate faster than 3H), while in those instances where 3H travels faster than chloride, vapor phase transport of tritiated water is inferred. The balance between these processes is sensitive to hydraulic properties and conditions that determine hydraulic conductivity (Scanlon and Milly, 1994). Also, as pointed out by Cook et al. (1994) and Tyler and Walker (1994), processes within the root zone affect transport differently, when considering either the center of mass of the respective isotope systems, or the chloride mass balance. Furthermore, evidence of preferential flow pathways playing a significant role in vadose zone transport adds to the difficulty of quantifying mass fluxes from these sparse isotope data sets (Dixon, 2001; Fabryka-Martin et al., 1998).

Nevertheless, the interpretation that the empirical evidence suggests climate change effects extend to the deepest levels of arid region vadose zones is consistent with analysis of the depth to which surface temperature perturbations would extend in soil systems. (Hillel, 1998) noted that the damping depth (the depth at which a temperature change signal is reduced by 1/e) can be described by

[1]
where Dh is the thermal diffusivity in units of cal/(cm s °C), and {omega} is the radial frequency (2{pi} s-1). The damping depth increases as the square root of the time duration of the cycle. For diurnal changes in temperature, the damping depth is approximately 0.1 m, while annual temperature cycles have a damping depth of slightly more than 2 m. Although not strictly cyclic, climate change patterns which possess characteristics typical of a time series extend over hundreds to tens of thousands of years. Damping depths for these changes would extend to depths in excess of hundreds of meters. Hillel (1998) noted, on the basis of this analysis, that temperature signals associated with long-term climate changes may be discernible. Cook et al., (1992) provide a similar analysis for the transport of 36Cl and 3H tracers in the unsaturated zone. They conclude that, on the basis of an analysis of dispersion- and diffusion-controlled transport, isotopic tracers may preserve a record of "antecedent climatic conditions", provided that certain constraints regarding infiltration flux, thickness of the vadose zone, and transport rate are met. Their analysis suggests that such signals of earlier climate conditions may persist for more than 1000 yr.

Despite the complexities of the chemical and isotopic systematics, these data sets are important for several reasons. Most importantly for the topic considered here, it is clear that the tracer and isotope data support the view that many deep vadose zone environments contain water that is very old, that predates the present climate regime (Phillips, 1994), and that records variation in groundwater recharge rates. Second, changes in groundwater recharge must reflect an interplay between temporal variation in water flux through the vadose zone and the hydrologic saturation state. Variation in either of these system variables will affect the local pore water chemical composition by either changing the residence time of water in contact with local mineral surfaces, or changing the water/rock ratio, which effectively leads to concentration or dilution effects.

Chemical Kinetics in the Vadose Zone
Norton (1984) was among the first to rigorously describe the nonlinear coupling between fluid movement and chemical evolution of an advecting hydrothermal system. Most of the changes he described were expressed as changes in solution composition, and dissolution and precipitation of mineral phases, in response to externally driven perturbations such as the thermal pulses associated with emplacement of magma bodies. He noted that the interactions expressed in such a system form a feedback loop in which the changes in the physical framework due to chemical dissolution or precipitation of minerals affect fluid movement and heat transport via modification of the permeability. This, in turn, affects where fluids migrate and thus modifies the locus of sites where chemical change takes place. The result is a strongly coupled, highly nonlinear dynamic system.

The same principles described by Norton apply equally well to a vadose zone environment. In this case, however, the perturbation would be changes in surface temperature and infiltration flux driven by climate change, rather than thermal perturbation driven by emplacement of a magma body. As noted above, only long-term variations in average properties would be capable of being propagated any significant distance through a thick (>100 m) vadose zone. For those changes in climate that occur over time scales of hundreds to thousands of years, however, propagation of chemical change through the vadose zone environment would be expected to occur in response to the propagating changes in temperature and moisture flux.

The magnitude of any chemical signal, and the rate at which it changes, would be small, however. The magnitude and rate of change of a chemical signal reflect the thermodynamic properties and reaction kinetics of fluid–rock interaction dominated by mineral dissolution and/or precipitation. Lasaga et al. (1994) examined in detail the role of reaction kinetics on weathering rates and geochemical cycles. Mineral dissolution rates at near-surface temperatures (i.e., <40°C) range from about 10-7 to 10-14 mol m-2 of specific surface area per second (Blum and Lasaga, 1988; Burch et al., 1994; Carroll et al., 1998; Carroll and Walther, 1990; Carroll-Webb and Walther, 1988; Chou and Wollast, 1985; Knauss and Wolery, 1986, 1988; Lin and Clemency, 1981; Nagy et al., 1991; Nagy and Lasaga, 1992; Rimstidt and Barnes, 1980; Rimstidt and Dove, 1986; Rose, 1991; Schott et al., 1989; Schweda, 1989; Tole et al., 1986). These rates are controlled by temperature, solution composition and ionic strength, exposed reactive mineral surface areas, and the net free energy change of the relevant dissolution and precipitation reactions. The relevant rate equation proposed by Lasaga et al. (1994) is

[2]
where Ko is the net rate constant, S the effective surface area for the i dissolving or precipitating phases, E the net activation energy for the overall reaction, ai the activities of the relevant aqueous species, and f({Delta}Gr) the functional relationship that describes the dependence of the rate on the deviation from equilibrium. By definition, equilibrium is the {Delta}Gr = 0 state.

Activation energies for reactions involving common mineral phases range between about 4.4 and 125.6 kJ mol-1 (2 and 30 Kcal mol-1), and average about 63 kJ mol-1 (15 Kcal mol-1). Lasaga et al. (1994) pointed out that such values for the average activation energy will result in approximately an order of magnitude change in reaction rate for every ~30°C temperature change at near-surface conditions. However, they also noted that the broad range in reaction rate constants as well as activation energies can lead to complex mineral paragenetic relationships in which metastable persistence of mineral phases will occur.

The impact of changing water flux on the concentration (c) of solutes, can be approximated, for systems far from chemical equilibrium (Lasaga et al., 1994) by the expression

[3]
where k is the relevant dissolution rate constant, and A/V is the kinetic water/rock ratio, which is a measure of the time-integrated exposed surface area to the rock volume. Changes in flux change A/V, thus, affect the evolution of the solute load.

Taken together, these system properties require that a change in either the infiltration flux or the surface temperature will propagate through the vadose zone as chemical change in the migrating pore water, and as precipitation or dissolution of minerals. In thick vadose zones composed of rocks with low hydraulic conductivities, the low infiltration fluxes and slow reaction kinetics imply that chemical changes initiated by climate changes hundreds or thousands of years ago are still propagating through the system. However, the nature of that response will be a reflection of the local, unique rock mineralogy, exposed surface areas, fluid flow pathways (e.g., matrix vs. fracture), hydraulic properties, and history of prior fluid changes.


   SIMULATIONS
TOP
 ABSTRACT
 INTRODUCTION
 Background
 SIMULATIONS
RESULTS
 DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Vadose Zone Water Chemistry and Mineralogical Change in Response to Climate Change
The considerations outlined above imply that climate-driven chemical effects may be present in deep vadose zone settings. However, the nature and magnitude of these chemical effects, and the response times of these geological systems remain unclear.

In order to quantitatively evaluate the characteristic response of a typical vadose zone to climate change, we conducted a series of simulations using the reactive transport simulator NUFT-C (Non-isothermal Unsaturated-saturated Flow and Transport—Chemistry; Glassley et al., 2001, 2003; Nitao, 1998). NUFT-C simulates multicomponent, multiphase, saturated and unsaturated reactive transport for nonequilibrium, nonisothermal systems. Account is taken of heat transfer via conduction and convection, pore water vaporization, evaporation and condensation as temperature changes, gas phase migration, humidity, saturation, fracture and matrix pore water and gas chemistry, and modification of the porosity and permeability due to dissolution and precipitation of primary and secondary minerals. Mineral dissolution and precipitation are modeled using temperature-dependent thermodynamic data (Johnson and Lundeen, 1994; Johnson et al., 1992) for the relevant hydrolysis and speciation reactions. Reaction rates are modeled using a transition state theory approach. These capabilities provide a quantitative description of the evolution of the permeability field as transport and chemical reactions occur, thus providing a direct feedback that modifies fluid flux. A more extensive description of the code is available in Glassley et al. (2003). Although designed to run on a massively parallel IBM SP-2 computer with 1200 processors (Mirin, 1998), allowing us to conduct high-resolution simulations of complex geochemical and hydrological systems, versions of the code are available that run on single processor work stations and personal computers.

The geological framework used in these simulations was that of Yucca Mountain, Nevada. This volcanic stratigraphic sequence was selected because detailed thermal–hydrological and mineralogical properties of the rock units are readily available in Yucca Mountain Project reports (Sonnenthal and Spycher, 2001 and references therein) and in Flint et al. (2001). The specific units considered in this study were the thermal–hydrological units (from top to bottom) Tcw13, Ptn21, Ptn22, Ptn23, Ptn24, Ptn25, Ptn26, Tsw31, Tsw32, Tsw33, Tsw34, Tsw35, Tsw36, Tsw37, Tsw38, and Ch2z, which compose nearly all of the about 550-m volcanic stratigraphy. In addition, a 2-m-thick overburden of alluvial sand was assumed to act as the interface between the volcanic stratigraphy and the atmosphere.

A dual porosity–dual permeability continuum model was employed in which each of the 17 lithologic units in the model is represented by unique fracture and matrix continua. Interactions between matrix and fracture continua are evaluated through functional relationships that consider bulk density, porosity, water retention (capillary pressure–saturation) behavior, the unsaturated soil hydraulic conductivity (permeability), and fracture abundance, volume, and effective exposed surface area. This approach allows independent tracking of fracture and matrix responses throughout the more than 500-m thickness of this vadose zone. The chemical system considered consisted of the components CaO–Na2O–K2O–Al2O3–MgO–SiO2–Cl–CO2–SO4–H2O, with and without the additional components Fe2+, Fe3+ and O2, and up to 34 minerals and 58 aqueous species. All of the water compositions were restricted to waters that would be classified as chloride-bearing calcium or sodium bicarbonate waters. In this report we discuss the key results from a subset of these simulations.

The protocol followed for these simulations was to select one of the reported pore water chemistries for Yucca Mountain rocks (Browning et al., 2000; Yang et al., 1996, 1998). An initialization run was then conducted, allowing the system to achieve steady-state conditions for saturation, temperature, pressure, and liquid and gas chemistries in all lithologic units. For our purposes, we defined "steady state" to be that condition for which the rate of change of the variables being considered was <1%, relative, over 10000 yr. After steady-state conditions persisted for more than 20000 yr, a new run was undertaken that used these steady-state conditions, but which then, after running for a few thousand years longer, was perturbed. The system was perturbed either by changing the infiltration flux (from 87.5 to 8.75 mm yr-1 over a period of 500 yr at a constant surface temperature of 17.7°C) or the surface temperature (from 12.7 to 17.7°C over a period of 200 yr at a constant infiltration flux of 8.75 mm yr-1). The final conditions in each simulation (8.75 mm yr-1 infiltration flux and 17.7°C surface temperature) were selected based on earlier simulations of Yucca Mountain thermal–hydrology behavior (e.g., Buscheck et al., 1999). The starting conditions (an elevated infiltration flux of 87.5 mm yr-1 that is 10 times the current value and/or lower surface temperature of 12.7°C, which is 5° cooler than present) were chosen so as to illustrate potential vadose zone behavior when the surface conditions evolve from those of a cooler, wetter climate such as that present during the Pleistocene pluvial period. We emphasize, however, that these initial temperature and flux conditions are selected solely for illustrative purposes and are not presented with the implication that the actual conditions at this location were, in fact, those used here.

In all simulations we assumed a constant temperature water table, which introduces a small error into the derived geothermal gradients described below. The infiltrating water was assumed to have the chemical composition of local rainwater (McKinley and Oliver, 1993). These simulations were run until steady state persisted for more than 20000 yr after perturbation. An additional simulation was run in which the change from pluvial conditions to the present climate was approximated by simultaneously changing the infiltration flux from 87.5 to 8.75 mm yr-1, and temperature from 12.7 to 17.7°C over a 10000-yr period. We selected these conditions on the assumption that they would be illustrative of the character of change that may have occurred during that time period (Spaulding, 1985). However, we recognize that the precise conditions 10000 yr ago are unknown, and changes in conditions would not have followed a smooth evolution.


   RESULTS
TOP
 ABSTRACT
 INTRODUCTION
 Background
 SIMULATIONS
 RESULTS
DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Although all properties in all cells of the computational mesh were monitored, for illustrative purposes, only the results at depths of 200 and 400 m are presented here.

Effects of Decreasing the Infiltration Flux Rate
A decrease in the infiltration flux results in an increase in overall temperature. A new steady-state condition is achieved over a period of approximately 10000 yr (Fig. 1). However, more than 90% of the total temperature change occurs within the first 5000 yr. The initial temperature rate of change is between 0.0012 and 0.0022°C yr-1, and decreases to about 0.0001 °C yr-1 by 5 000 yr. The new steady-state temperature is approximately 3°C warmer at 200 m, and approximately 3.5°C warmer at 400 m, resulting in a steeper vertical temperature gradient (2.0°C per 200 m) than at the higher flux rate (1.5°C per 200 m). These effects reflect the fact that infiltrating water tends to suppress the natural geothermal gradient by absorbing and advecting heat downwards. Our assumption of a constant water table temperature has virtually no effect on the geothermal gradient we compute from the change in temperature over the 200- to 400-m depth interval, but would introduce a few percent error in an average geothermal gradient computed from the ground surface to the water table.



View larger version (19K):
[in this window]
[in a new window]
  		Fig. 1. Temperature as a function of time at depths of about 200 and 400 m. Dashed red lines labeled {Delta}Temp indicate temperature change due to an increase in surface temperature from 12.7 to 17.7°C, at a constant infiltration flux of 8.75 mm yr-1. Solid blue lines labeled {Delta}Flux indicate temperature change due to a decrease in infiltration flux from 87.5 to 8.75 mm yr-1, at constant surface temperature of 17.7°C. The rate of temperature change for the portion of the curve indicated by the arrow is 0.002°C yr-1. The curves have been plotted such that the time at which the perturbations were initiated exactly correspond. Approximately 4000 yr of the initial steady-state conditions are shown.

 
At the new infiltration flux of 8.75 mm yr-1, complete replacement of in situ fracture water is accomplished within ~2000 yr at both depth intervals (~200 and ~400 m). The largest fluid reservoirs in the system, however, are the matrix blocks, accounting for more than 90% of the total volume. Exchange of water between the matrix and fracture continua is primarily due to capillary imbibition. As a result, complete pore water replacement takes much longer in the matrix, requiring more than 13000 yr at ~200 m and 21000 yr at ~400 m.

Infiltration flux affects solution chemistry by changing the ambient temperature and the residence times and volumes of solutions in fracture and matrix pores. As temperatures evolve, aqueous speciation reactions and dissolution or precipitation of mineral phases occurs as a function of the relevant, temperature-dependent mass action expressions. The rate at which chemical change occurs is dependent upon the exposed surface areas of the solid phases in contact with the solution, as well as the degree of mineral supersaturation or undersaturation. Given that the dissolution and precipitation rate constants for the minerals in this study range between 10-7 and 10-14 mol m-2 s-1, and that initial mineral abundances and surface areas in each lithologic unit are different, solution evolution will be highly complex, and differ among and within lithologic units. The time required for solution composition to recover to a steady state will, therefore, be no shorter than the time required to achieve steady-state temperature and liquid saturation conditions.

Typical examples of this complex behavior can be observed for the evolution of aqueous SiO2 and bicarbonate. Aqueous SiO2 concentrations reach steady-state values somewhat after temperature steady-state conditions are achieved (compare Fig. 1 and 2). At the 200-m depth, the change in conditions leads to fracture and matrix pore waters attaining the same silica concentration, while the deeper rock unit at 400 m maintains a large difference in fracture and matrix pore water silica concentration. However, for both depth intervals, the absolute changes are small, with the molality increasing by <10-5.



View larger version (23K):
[in this window]
[in a new window]
  		Fig. 2. Negative log of the aqueous SiO2 concentration in fractures and matrix as a function of time at about 200 and 400 m, for the conditions represented in Fig. 1. The negative log of the aqueous SiO2 concentrations associated with surface temperature change are plotted as red lines, and those associated with surface flux change are plotted as blue lines. Fracture waters are represented as dashed lines, and matrix waters by solid lines. The curves have been plotted such that the time at which the perturbations were initiated exactly correspond, and ~4000 yr of the initial steady-state conditions are shown.

 
Bicarbonate concentrations follow a different evolutionary path as infiltration flux changes. The HCO-3 concentrations require about twice as much time to achieve new steady-state values (Fig. 3), compared with the time required for changes in aqueous silica, and the increases are up to an order of magnitude larger (2 x 10-4 or greater).



View larger version (23K):
[in this window]
[in a new window]
  		Fig. 3. Negative log of the aqueous HCO-3 concentration in fractures and matrix as a function of time at about 200 and 400 m, for the conditions represented in Fig. 1. The color scheme, line patterns and layout are the same as for Fig. 2.

 
Controls on solution chemistry are complex. This can be seen by examining changes in the saturation index, Q/K, where Q = / and K is the mass action constant for the relevant hydrolysis reaction at the temperature of interest, for a suite of primary and secondary minerals. Before and after the change in flux, pore waters maintain approximately constant degrees of undersaturation with respect to analcime, cristobalite, quartz, and dolomite at 200 and 400 m (Fig. 4), and remain strongly supersaturated with respect to kaolinite. Calcite, on the other hand, goes from undersaturated to supersaturated in fractures, but remains undersaturated in the matrix at 200 m, while, at 400 m, fractures remain supersaturated and the matrix becomes undersaturated after initially being saturated (Fig. 4). Potassium feldspar at 200 m remains supersaturated in the fractures and undersaturated in the matrix, but is undersaturated in both at 400 m. Complex responses can also be seen for illite, mordenite, and smectite, all of which are strongly supersaturated in fractures and matrix at 200 m, but have contrasting behaviors at 400 m.



View larger version (14K):
[in this window]
[in a new window]
  		Fig. 4. Changes in mineral saturation indices (Q/K) at 200- and 400-m depths for selected minerals, for the change in surface infiltration flux scenario. Saturation indices for fracture waters are indicated by square symbols, and matrix water by round symbols. Open black symbols indicate the initial saturation indices, and the red filled symbols indicate the new steady-state saturation indices. The minerals represented in the figure are labeled along the horizontal axis. Values of Q/K <1.0 indicate solutions that are undersaturated with respect to the indicated minerals, and which would dissolve the mineral if present, while values >1.0 indicate supersaturated waters from which the indicated minerals would precipitate. Illite, mordenite, and smectite are strongly supersaturated (Q/K {approx} >2.0) in matrix and fractures at 200 m, but exhibit more complex relationships at 400 m. Note that Q/K for illite and smectite initially exceeds 2.0 in the matrix at 400 m. Q/K for kaolinite is consistently >>1.0 under all conditions.

 
Temperature Effect
Temperature changes at the surface affect deeper level thermal conditions via thermal conduction through the matrix continuum and via the heat capacity of the advecting fracture and matrix water. The time scale for achieving steady-state conditions is the same as for the changing flux case, since the temperature change mechanisms are the same. The magnitude of temperature change is depth-dependent. For this simulation in which the surface temperature was increased by 5°C, {Delta}T at ~200 m is about +4.4°C, and at ~400 m is about +3.2°C (Fig. 1). This results in an overall decrease in the thermal gradient over this depth interval for a period of time exceeding 20 000 yr.

Changes in pore water compositions are greater for this scenario than for the change in flux condition. This reflects the consequences of several simultaneous processes, including differences in initial temperature states, the impact of changing temperature on reaction rates, and in initial saturations. The net result is that changing the surface temperature by 5°C, from 12.7 to 17.7°C at an infiltration flux of 8.5 mm yr-1, results in changes in aqueous SiO2 and HCO-3 molalities of about 300% (Fig. 2 and 3). Although the overall chemical behavior of the different rock units at the ~200- and ~400-m depths are similar, the magnitude of the responses is different. For the shallower level rock, the contrasts between fracture and matrix aqueous SiO2 ({Delta}logf-m SiO2) decreases from ~0.16 to ~0.03, while {Delta}logf-m HCO-3 decreases from ~0.025 to ~0.015. For the deeper rock, {Delta}logf-m SiO2 increases from ~0.05 to ~0.12 and {Delta}logf-m HCO-3 decreases from ~0.05 to ~0.03. This results in greater differences between fracture and matrix pore waters in aqueous SiO2 for the deeper rock, but smaller differences for the shallower rock once the new steady-state conditions are achieved. The chemical changes exhibited by the pore waters result from a combination of dissolution and precipitation reactions involving the local mineral assemblages. Although these changes in mineral abundance would be difficult to detect through standard analytical means, since the net abundance change for an individual participating mineral phase is <10-4 volume fraction over about 10000 yr, microscale surface textures may show features consistent with dissolution or precipitation.

The rate of change of the carbonate system contrasts significantly between fracture and matrix, compared with that for aqueous SiO2. For the latter, within each rock unit the response times for fracture and matrix are similar, whereas in the carbonate system, the response times for the fracture water chemistries are similar to that of the aqueous SiO2 system, while that of the matrix is dramatically slower, taking ~25000 yr to reach new steady-state conditions (Fig. 3).

Multiple Long-Term Changes
The results of the numerical experiment in which simultaneous changes in the infiltration flux (85 to 8.5 mm yr-1) and surface temperature (12.7 to 17.7 °C) were assumed to occur over a 10000-yr period are shown in Fig. 5Go through 7. The end point surface conditions are the same in this case as for the numerical experiments described above.



View larger version (15K):
[in this window]
[in a new window]
  		Fig. 5. The rate of temperature change (°C) as a function of time (years) vs. depth, for the case in which surface temperature and infiltration flux change simultaneously over 10000 yr, after 20000 yr of steady-state conditions (line labeled 19900 yr). The line labeled 30160 yr indicates the rate of temperature change 160 yr after climate change has ceased, and approximately defines the maximum values achieved for dT/dt. Note that even after about 10000 yr after perturbation at the surface, the rate of temperature change is still >0 (line labeled 39940 yr).

 


View larger version (23K):
[in this window]
[in a new window]
  		Fig. 6. Negative log of the aqueous SiO2 concentration in fractures (dashed lines) and matrix (solid lines) as a function of time at about 200 (red) and 400 (blue) m, for the conditions and time span represented in Fig. 5. Approximately 10000 yr of the initial steady-state conditions are shown. The period of time over which the infiltration flux and surface temperature changes occur is indicated by the shaded region labeled Perturbation Period.

 


View larger version (24K):
[in this window]
[in a new window]
  		Fig. 7. Negative log of the aqueous HCO-3 concentration in fractures (dashed lines) and matrix (solid lines) as a function of time at about 200 (red) and 400 (blue) m, for the conditions and time span represented in Fig. 5. Approximately 10000 yr of the initial steady-state conditions are shown. The period of time over which the infiltration flux and surface temperature changes occur is indicated by the shaded region labeled Perturbation Period.

 
The temperature evolution over this period is slow and progressive, requiring ~10000 yr after surface temperature change has ended to reach a new steady state (Fig. 5). The maximum rate of temperature change is ~0.0009°C yr-1, which is the same maximum rate achieved for the case in which only temperature was perturbed.

Changes in fracture saturations occur within a few hundred years after perturbation is initiated, and achieve a new steady state within a few hundred years after the perturbation ceases.

These combined changes lead to complex and divergent evolutionary pathways for the fracture and matrix chemistries. There is a small decrease in the log of the aqueous SiO2 concentration in the fractures ({Delta}logf SiO2 ~0.02–0.06), but a relatively large increase ({Delta}logm SiO2 ~0.25) in the matrix (Fig. 6). The net result, at steady state, is that the aqueous SiO2 concentration ratio between fracture and matrix pore waters changes from ~1.5 to 0.8. The primary change affecting the matrix water is an increase in temperature, since matrix saturations are always very high, thus leaving temperature as the principal physical change affecting matrix waters. Fractures, on the other hand, experience a significant decrease in saturation (nearly 50%), but an increase in temperature, resulting in a more complex interaction in controlling parameters.

The log of the bicarbonate concentration ({Delta}logf and m HCO-3) increases in both fractures and matrix by ~0.15 (Fig. 7). At the point steady-state conditions are realized, the ratio between fracture and matrix bicarbonate concentrations recovers to approximately its original value. However, for more than half of the transition period, the fracture and matrix concentrations reverse their relative values. Note that the end point steady-state concentrations for this system are different from those in the previous simulations, even though the end point surface conditions are all identical.


   DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 Background
 SIMULATIONS
 RESULTS
 DISCUSSION
CONCLUSIONS
 REFERENCES
 
The results of these simulations demonstrate that climate change will affect the long-term thermal, hydrological, and geochemical characteristics of deep vadose zone systems, such as those commonly found in the American Southwest. The characteristic response times for these systems to achieve new steady-state conditions after a brief perturbation are on the order of thousands of years, but differ between fracture and matrix. Recent studies of climate variability in the U.S. Southwest documents significant climate fluctuations on many different time scales, some as short as about 100 yr (Allen and Anderson, 2000; Dean et al., 1996; Hostetler and Benson, 1990; Krider, 1998; Oviatt, 1997; Spaulding, 1985; Waters, 1989). This implies that the entire vadose zone in the Southwest may well be in the process of chemically and thermal–hydrologically adjusting to relatively recent, postglacial climate changes, and is not at steady state. This result supports previous suggestions that deep vadose zone settings may be in a nonequilibrium, dynamic chemical state (Phillips, 1994). Inverse modeling using data derived from thick vadose zone temperature, saturation, and pore water chemistry measurements to reconstruct past climate may well yield a rich and detailed history.

Reconstructing past climate conditions from vadose zone date, however, faces significant obstacles. These obstacles derive from limitations in available data and numerical models.

Data Limitations
Vadose zone temperature conditions, saturation states, and pore water chemistries are strongly influenced by surface temperature and infiltration flux. Quantifying these parameters requires detailed knowledge of a host of thermal, hydrological, and geochemical and mineralogical rock properties (e.g., lithology-specific heat capacity, thermal conductivity, porosity, permeability, saturation-dependent moisture retention characteristics, and pore-lining vs. bulk rock mineralogy) that are usually not well characterized on the small scale needed to accomplish such work. Some of the chemical and mineralogical changes that are evident in the simulations we have conducted are too small to be analytically detected using currently available analytical methods. There is a paucity of rate constant data for most mineral phases of interest that would be thermodynamically stable under these vadose zone conditions. Mineral reactive surface areas are currently impossible to quantify for in situ geological materials, and will change during reaction in ways that are currently unpredictable. The actual flow pathways will be controlled by complex, and thus far indeterminate, three-dimensional permeability structures, and these will ultimately control saturation states, reaction progress, and the distribution of secondary mineral assemblages. Finally, developing accurate and precise initial and boundary conditions for simulations will be difficult to accomplish. Each stratigraphic unit possesses a unique history of mineralogical and geochemical evolution. This evolutionary pathway includes an inventory of mineral phases and their respective reactive surface areas that have developed in response to a multi-million year history of fluid–rock interaction.

Model Limitations
Representing complex, heterogeneous rock systems by homogeneous continua introduces errors in the rate at which chemical steady-state conditions are achieved and the spatial distribution of mineral precipitation and dissolution (Glassley et al., 2002), and in the hydrological interactions between fracture and matrix (P. Lichtner, personal communication, 2000; Pruess et al., 1999). Developing the ability to more realistically represent properties of complex geological systems would significantly improve the realism of computer simulations. Key enhancements include: implementation of discrete fracture models that take into account the mineralogical and physical properties of fracture-lining and fracture-bounding "skins", discrete representation of mineral and hydrological heterogeneity in matrix blocks, and soil permeability changes in response to climate change that impact percolation and flow pathways and water–mineral reactions and reaction rates. Soil–atmosphere interactions are particularly important under those conditions in which carbonate-silica precipitation and/or dissolution occur in response to climate change. These changes in soil properties can lead to relative rapid and large changes in soil permeability that may strongly impact infiltration flux. Also needed is a better understanding of the physics controlling fluid movement at low saturations. Although water fingering along fracture surfaces has been observed experimentally, it is not currently possible to represent this process in simulations. The same thing can be said regarding the development of preferential flow pathways in matrix blocks. Such processes have important implications for how best to treat effective reactive surface areas and fluid fluxes.

Even so, although there is much work that needs to be done to achieve tightly constrained reconstructions of past climates from vadose zone environments, current attempts to extract climate records from these settings can still be fruitful. Detailed stable and radiogenic isotope studies can be used to elucidate rates of surface reactions, and effective exposed surface areas provide some insight into rock–water interaction histories and fluid flow pathways (Johnson and DePaolo, 1997; Marshall and Mahan, 1994; Paces et al., 1996, 1997; Whelan et al., 1994, 1996; Whelan and Stuckliss, 1992). This information, coordinated with high resolution reactive transport simulations, detailed characterization of mineral surface textures to elucidate dissolution and precipitation relative histories, and detailed analytical studies of microsamples of pore fluids could lead to the development of detailed chronologies of moisture flux and fluid chemistry and may provide the means to constrain regional patterns of surface temperatures up to about 100000 yr in the past.


   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 Background
 SIMULATIONS
 RESULTS
 DISCUSSION
 CONCLUSIONS
REFERENCES
 
The results of the simulations we present here provide quantitative support for the view that the vadose zone may possess a rich paleoclimate record that may encompass tens to hundreds of thousands of years. But the results also underscore the substantial data needs and improvements in conceptual models that must be satisfied to decipher these records. Multiphase reactive transport simulators, such as the one employed here that rigorously accounts for coupling between thermal, hydrological, and geochemical processes under unsaturated conditions, are essential in order to interpret pore fluid chemistries, water saturation, and temperature conditions in terms of surface temperature and infiltration flux histories. However, such codes rely on exhaustive empirical data that are often not routinely obtained in the course of vadose zone studies (e.g., reactive surface areas, fracture and matrix hydraulic properties, and detailed mineralogies). In addition, these results suggest that the water characteristics are path dependent, thus posing a problem for defining initial and boundary conditions for inverse modeling. The solution to this problem lies in collaborative efforts in which isotopic studies and microcharacterization of mineral surface textures along three-dimensional flow pathways are combined with detailed simulations to guide construction of conceptual models and define constraints for application of inverse modeling techniques. Particularly important will be development of higher resolution studies of permeability changes associated with climate effects in the shallow soil zone where mineral dissolution and precipitation can significantly impact infiltration flux.


   ACKNOWLEDGMENTS
 
Comments from Nina Rosenberg and Charles Carrigan on an early version of the manuscript, and thorough reviews by Scott Tyler and an anonymous reviewer greatly improved this paper and are gratefully acknowledged. Editorial handling by Bridget Scanlon refined the manuscript significantly and is gratefully acknowledged. This work was performed under the auspices of the U.S. Department of Energy by University of California Lawrence Livermore National Laboratory under contract number W-7405-Eng-48.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 Background
 SIMULATIONS
 RESULTS
 DISCUSSION
 CONCLUSIONS
 REFERENCES
 





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