| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
Lawrence Berkeley National Lab., Earth Sciences Division, 1 Cyclotron Rd., MS 90-1116, Berkeley, CA 94720
* Corresponding author (SAFinsterle{at}lbl.gov).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Received 28 March 2007.
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONVadose Zone Modeling ChallengesSelect Capabilities and...ConclusionsREFERENCES |
|---|
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONVadose Zone Modeling ChallengesSelect Capabilities and...ConclusionsREFERENCES |
|---|
|
|
|
Vadose Zone Modeling Challenges |
|---|
|
TOPABSTRACTINTRODUCTIONVadose Zone Modeling ChallengesSelect Capabilities and...ConclusionsREFERENCES |
|---|
In addition to considering unsaturated flow of water and transport of tracers or nonreactive components, vadose zone hydrology has evolved to study complex processes of nonisothermal, multiphase, multicomponent fluid flow coupled to biogeochemical reactions. The multiscale, multiphysics nature of vadose zone processes poses significant conceptual and numerical challenges. In addition, advanced measurement techniques and analysis methods are needed to obtain the characterization data for site-specific models. Finally, the vadose zone has to be integrated into a larger system, which may include not only other natural subsystems, but also engineering components as well as economic and regulatory aspects.
The modeling challenges mentioned above can be grouped into three categories: (i) challenges resulting from the heterogeneity and complexity of the physical and biogeochemical processes occurring in the vadose zone; (ii) challenges resulting from the position of the vadose zone at the interface between different spheres (atmosphere, biosphere, and hydrosphere); and (iii) computational challenges resulting from the coupling of subsystems and processes involving many spatial and temporal scales.
Challenge 1: Complexity
The first group of challenges includes the characterization and incorporation of heterogeneity, which affects flow behavior (fingering, flow focusing, lateral diversion, preferential and fast flow, and water perching) and transport behavior (channelization or plume spreading). Moreover, not only do biogeochemical reactions depend on heterogeneity in chemical and mineralogical composition, but they are also affected by details of the flow patterns and transport processes to which they are strongly coupled. Also note that under unsaturated or multiphase flow conditions, the impact of heterogeneity on these processes may be self-enhancing or self-controlling, which affects the conceptualization of such processes in a numerical model.
The characterization of heterogeneity is further complicated because the soil structure may exhibit both systematic and random features on multiple scales. Additionally, sensors used for monitoring the vadose zone and for measuring its properties have different support scales, and averaging these properties often leads to effective parameters that are process specific, scale dependent, and anisotropic.
The TOUGH simulators have been used extensively to analyze the impact of heterogeneity on vadose zone processes (Pruess, 1998, 1999, 2004a; Pruess et al., 2002). Such studies include the use of high-resolution simulations with a massively parallelized version of the code (Wu et al., 2002; Zhang et al., 2003) and the development of a joint hydrologic–geophysical inversion approach to characterize heterogeneity and to reduce systematic modeling errors (Kowalsky et al., 2005; Finsterle and Kowalsky, 2008; Lehikoinen et al., 2007).
The classic approach for describing unsaturated flow based on the Richards equation has been successful but limited. One limitation stems from the difficulty in obtaining characterization data for specifying the relative permeability and capillary pressure curves. In addition, the functional form of these curves is difficult to assess and may greatly impact modeling results. For example, if a standard hysteretic capillary pressure model is used, key phenomena (such as hysteresis in the relative permeabilities and history-dependent gas entrapment) may be neglected. The hysteresis model implemented in TOUGH was presented in Doughty (2007) and will be summarized below. Similarly, fast flow through macropores or fractures may be described by an effective continuum model with multimodal characteristic curves. As demonstrated by Doughty (1999), however, the applicability of an effective continuum model is limited to cases where the assumption of thermodynamic equilibrium between the continua is valid. In addition to the active fracture model by Liu et al. (1998), the TOUGH simulators include dual- and triple-continuum (Wu et al., 2004) formulations with multiple interacting continua to accurately capture exchanges of fluids and heat between regions of the porous medium with vastly different flow and transport characteristics.
Challenge 2: Coupling
The vadose zone, being bounded by the land surface and the groundwater table, is in direct contact with the biosphere and strongly affected by human activities. Interaction among these various interfaces leads to a second group of modeling challenges. The processes at the vadose zone boundaries are often highly dynamic (infiltration events, evapotranspiration, heat transfer, water table fluctuations); more importantly, they may need to be fully coupled to account for feedback mechanisms (e.g., evaporation and root water uptake). The TOUGH codes accurately account for phase transitions (evaporation and condensation, including the related latent heat effects) as well as multiphase diffusion (e.g., of water vapor), allowing the simulation of evaporation and drying processes. Simplified approaches to simulate evaporation using the Richards equation have recently been implemented (Ghezzehei et al., 2004). TOUGH has also been coupled to the National Center for Atmospheric Research Community Land Model CLM3 (www.cgd.ucar.edu/tss/clm/index.html; verified 29 Mar. 2008) for the simulation of energy and moisture dynamics between the atmosphere and the subsurface (Pan et al., 2008). Finally, integration of TOUGH2 and iTOUGH2 into the GoldSim model (GoldSim Technology Group, 2006) provides the means to study the interaction between the subsurface environment and related natural and engineered components on a system level (Zhang et al., 2007).
Challenge 3: Numerics
The third group of challenges is related to the computational difficulties in dealing with highly nonlinear, coupled processes that occur on multiple scales and with potentially different characteristic time constants. Spatial discretization in the TOUGH codes is based on the integral finite difference approach (Narasimhan and Witherspoon, 1976), in which the basic mass conservation equations are directly discretized in their integral form (Pruess et al., 1999; Pruess, 2004b). This scheme provides great flexibility in handling complex geometries, multiregion approaches, and a variety of boundary conditions. Higher order schemes have also been implemented (Pruess, 1991a; Oldenburg and Pruess, 2000). Nonlinearities are treated using Newton–Raphson iterations with a residual-based convergence criterion. Parallelization of the forward (Wu et al., 2002) and inverse (Finsterle, 1998) runs help make tractable the solution of larger simulation and optimization problems. Nevertheless, phase changes, strong nonlinearities in characteristic curves near residual saturations, sharp saturation and reaction fronts, and the coupling of counteracting processes on disparate scales are among the features causing numerical difficulties that require careful attention by the code developer and sometimes case-by-case intervention by the user.
The following examples address diverse issues from all three groups of challenges related to vadose zone modeling mentioned above.
|
Select Capabilities and Applications |
|---|
|
TOPABSTRACTINTRODUCTIONVadose Zone Modeling ChallengesSelect Capabilities and...ConclusionsREFERENCES |
|---|
During periods of CO2 injection into brine-saturated formations, the resulting CO2 plume grows continuously; that is, all locations follow the primary drainage branch of the capillary pressure curve at all times, and this branch can be replicated using a nonhysteretic formulation. For post-injection periods, however, when the CO2 plume moves upward due to buoyancy forces, different locations experience drainage and imbibition at different times, necessitating the use of a hysteretic formulation.
In a hysteretic model, some parameters depend only on the process (drainage or imbibition) that is occurring, so it is convenient to subdivide the characteristic curves into drainage branches and imbibition branches. The parameters describing gas entrapment, however, depend on the value of the saturation when the grid block makes a transition from drainage to imbibition or vice versa, the so-called turning-point saturations. Because turning-point saturations are different in each grid block, these parameters are spatially variable and history dependent. The most critical parameter is the residual gas saturation, denoted Sgr
, which is the saturation below which gas is immobile (i.e., the saturation below which immiscible CO2 is trapped). For the primary drainage curve, Sgr = 0, but for imbibition, Sgr increases as the liquid saturation at the drainage-to-imbibition turning point, denoted Sl, decreases. Thus, grid blocks that once contained the most CO2 are those that trap the most CO2. The maximum possible value of Sgr is Sgr max, which is obtained for the minimum possible value of Sl; this minimum turning-point saturation is generally equal to the irreducible liquid saturation Slr. The value of Sgr determines the history- and location-dependent amount of CO2 that can be trapped. Note that gas entrapment cannot be adequately simulated using a model that is based solely on the Richards equation, in which gas is considered a passive bystander; a two-phase formulation is needed (Faybishenko, 1995).
We consider an application in which CO2 is injected into a porous formation 100 m thick located at a depth of 1000 m. This system is represented by an axisymmetric model with the injection well in the center. The porosity of the formation is 25%, horizontal permeability is 2 x 10–13 m2, and vertical permeability is 1 x 10–13 m2. The overburden above the storage formation extends to the surface and is assumed to have the same properties as the storage formation itself; thus, there is no low-permeability caprock. Initially, the brine saturation is 100% everywhere in the model, pore pressure is hydrostatic with a pressure of 100 kPa at the surface, and temperature follows the geothermal gradient of 30°C km–1, with the temperature at the surface and base of the model held constant at 15 and 45°C, respectively. The salinity of the pore water is assumed to be 100,000 mg L–1. Fluid and heat flow are fully coupled in the simulations.
The numerical simulations begin with injection of 900,000 t of CO2 into the porous formation at a constant rate of 30,000 t d–1 for 30 d. (This rate of CO2 corresponds roughly to emissions of a 1000 MW coal-fired power plant.) After injection stops, the only driving force in the model tending to cause movement of the CO2 is buoyancy. Simulations continue for 1000 yr.
Figure 1 shows the CO2 plume at a series of times. During the 1-mo injection period, the CO2 plume is expanding in all directions, and the hysteretic model results do not differ from what would be obtained with a nonhysteretic model. Thereafter, the leading edge of the plume, where drainage occurs and Sgr = 0, continues to advance (reaching the surface at about 700 yr), while the trailing edge of the plume, where imbibition occurs and Sgr is large, remains largely trapped. This combination of processes cannot be replicated with a nonhysteretic model or a hysteretic model that does not include phase trapping. If a nonhysteretic model with a small value of Sgr is used, most of the plume escapes through the ground surface within 10 yr; whereas a nonhysteretic model with a large value of Sgr produces a plume that never reaches the surface and remains entirely trapped indefinitely.
|
|
(shown by arrows); as |Pc| 
0 and krg 0 on the imbibition curves, Sl (1–Sgr) (shown by black-outlined dots). Thus, grid blocks near the plume center, which get much drier during the injection period and therefore have a small Sl, have a much larger Sgr, and consequently trap more CO2 than do grid blocks barely reached by the plume. More details about this modeling capability and the related application can be found in Doughty (2007).
Hydrological Processes in Permafrost
Subsurface flow is governed by the characteristics of the porous medium, the fluid properties, and constitutive relations describing the interaction between the porous medium and the fluids. The equation-of-state modules of the TOUGH codes provide an accurate description of the thermophysical properties of the considered fluid mixtures. Specifically, the properties of water are calculated within experimental accuracy from steam table equations as given by the International Formulation Committee (1967). As part of the ongoing reengineering of TOUGH, the temperature and pressure range available for the calculation of water properties has been expanded to include the phase transition from liquid or gas to ice and vice versa. The enthalpy, sublimation pressure, as well as fusion and melting pressures of ice (on the ice–vapor and ice–liquid equilibrium lines of the water phase diagram) are computed using regression equations from data obtained by the National Institute for Science and Technology. Within the ice phase (to temperature T = 50 K and pressure P 
200 MPa), ice densities are determined using the ice compressibility model of Marion and Jakubowski (2004) and the thermal expansivity data from Tanaka (1999). The ice enthalpy is computed using the heat capacity polynomial equation with the coefficients reported in Yaws (1999). Changes in fluid capillary pressure and relative permeability relations as a function of ice saturation are also taken into account.
Simulating freezing and melting processes is essential to improve our understanding of permafrost regions (e.g., to study the hydrologic response to climate change). These processes also occur during the formation and dissociation of gas hydrates deposited in permafrost and suboceanic formations (Moridis, 2003; Moridis et al., 2005). The extended equation of state for water is implemented in the reengineered version of TOUGH2, named TOUGH+. To demonstrate the capability of the TOUGH+GasH2O code to handle melting ice and steam fronts in a single model, we simulated the extreme thermal conditions encountered when a heat-generating power source (e.g., from a landing vehicle) is buried into permafrost on Mars. A heat source of 250 W was assumed to be embedded in the shallow Martian subsurface. Below a depth of 0.2 m, the pore space contained 50% ice and 50% CO2 at an atmospheric pressure of only 713 Pa and a temperature of –80°C; above the permafrost is a layer of completely dry sediment. Figure 3 shows the distribution of temperature, water and ice saturations after 10 sols (Martian days, 1 sol = 88,620 s). In the immediate vicinity of the heat source, the ice has been melted and the water vaporized. Liquid water only exists in a relatively thin region with temperatures between about 0 and 60°C. At the outer edge of this region, the water freezes, creating an ice barrier that effectively prevents convective heat transport in the radial direction. Water vapor (and associated thermal energy) is redirected upward, where it escapes to the Martian atmosphere. The simulations predict that the zone disturbed by a buried heat source is very limited. Specifically, liquid water is confined to a narrow region, reducing the risk that terrestrial microorganisms that may have been carried by the landing vehicle could survive and proliferate. More details about this simulation study can be found in Moridis and Pruess (2006).
|
Since most bacteria grow within a relatively immobile biofilm on solid surfaces, a three-region model is proposed that consists of (i) a mobile region, (ii) an immobile region, which includes stagnant water and biomass, and (iii) a solid particle region where mineral dissolution or precipitation and surface reactions may occur (Fig. 4 ). Instead of explicitly considering every pore and particle, the three regions are lumped into three overlapping continua, which are discretized separately and connected to each other. If necessary, each region can be further discretized to better resolve steep gradients and sharp reaction fronts. This concept is similar to the MINC model (multiple interacting continua; Pruess and Narasimhan, 1985) and the triple-continuum approach of Wu et al. (2004), both implemented in the TOUGH simulators. These models resolve global fluid and heat flow through a fractured (or macropore) system, while accounting for the local exchange between the fractures and the matrix. The extension of the MINC method to reactive geochemical transport was described in Xu and Pruess (2001).
|
|
Parameter Estimation and Model Structure Identification
An accurate description of physical processes, biogeochemical reactions, and fluid properties is a prerequisite for a reliable prediction of vadose zone systems. A major source of prediction uncertainty, however, lies in our incomplete knowledge of the subsurface structure and related soil properties. The complexity of soil formation as well as depositional and erosional processes result in an equally complex, usually highly heterogeneous, but not entirely random subsurface structure. While several stochastic methods exist to describe and simulate subsurface structures (using, for example, geostatistics [Deutsch and Journel, 1992] or lithofacies modeling [Carle and Fogg, 1996]), the determination of the actual structure at a given site remains challenging. In addition to identifying the geometry of the architecture of the soil strata, hydrogeologic and geochemical parameters need to be assigned. These parameters are often scale dependent, specific to the process being simulated, and related to other aspects of the conceptual model. Finally, initial and boundary conditions are either unknown or uncertain, so they need to be parameterized and estimated along with the structure and soil properties. Inverse modeling is a means to obtain parameters that in some sense can be considered optimal for a given model. It is important to realize, however, that an error in the conceptual model can yield parameter estimates that are biased or even meaningless; it is therefore essential to address the impact of systematic errors (in the data and the conceptual model) when calibrating a vadose zone model.
The iTOUGH2 code (Finsterle, 2004; www-esd.lbl.gov/iTOUGH2 [verified 29 Mar. 2008]) provides inverse modeling capabilities for most modules of the TOUGH suite of simulators. The iTOUGH2 solves a nonlinear optimization problem to determine TOUGH input parameters based on measured data for which corresponding TOUGH output is calculated. This means that any aspect of the model that can be parameterized—including boundary conditions and the soil structure—can be subjected to estimation by inverse modeling. Moreover, standard TOUGH output variables can be used to calculate new quantities that can then define the objective for optimization—including geophysical data or even costs.
This flexibility has been exploited to develop a joint inversion approach, in which hydrologic and geophysical data are inverted to estimate hydrogeologic, geostatistical, and petrophysical parameters (Kowalsky et al., 2004, 2005; Finsterle and Kowalsky, 2008). In this approach, the potential estimation bias introduced by a systematic error in the conceptual model (specifically in the representation of subsurface heterogeneity) is reduced by adjusting the model structure itself during the optimization process. The increased need for data that contain information about the subsurface structure is met by adding geophysical data. The inclusion of hydrogeologic data in the joint inversion approach ensures that a relation can be established between the geophysical signals and soil properties that determine fluid flow (note that this may require the estimation of additional parameters, such as those of a petrophysical relationship). An example is shown in Fig. 6 , where a highly heterogeneous, anisotropic soil structure and its related properties are determined by the joint inversion of hydrologic data (infiltration rates and neutron probe measurements) and geophysical data (arrival times from a cross-hole ground penetrating radar survey). Details about this specific application can be found in Finsterle and Kowalsky (2008).
|
Coupling with Other Models
Since the vadose zone is in direct contact with the land surface, it is essential to consider hydrologic and chemical interactions between the atmosphere, land surface, and subsurface. In addition, anthropogenic activities (such as land cultivation, irrigation, hydraulic engineering projects, and waste disposal) may profoundly affect the vadose zone and other related natural systems. These coupled impacts on the water and chemical cycles also influence (and are influenced by) the economic and regulatory environment.
Considering the vadose zone as a disconnected subsystem raises difficult questions regarding the boundary conditions to be used at the land surface. While net infiltration can be estimated using experimentally based parameterization models (for a review, see Faybishenko, 2007), evapotranspiration (a potentially significant contributor to the moisture balance at the surface, specifically in arid regions) is a coupled process that requires a physical description of both surface and subsurface conditions. Conversely, land hydrologic responses to meteorological forcing involve complicated exchanges of moisture and energy between soil, vegetation, snowpack, groundwater, and the overlying atmospheric boundary layer. The National Center for Atmospheric Research Community Land Model CLM3 is a model primarily developed to meet the needs of regional climate modeling. In CLM3, radiation, sensible and latent heat transfer, zonal and meridional surface stresses, and ecological and hydrologic processes are simulated as interrelated subprocesses; however, the subsurface moisture flow in CLM3 is only considered in a simplified manner. The CLM3 model has been coupled to TOUGH2 through an internal interface that includes flux and state variables shared by the two submodels. Specifically, TOUGH2 uses infiltration, evaporation, and root uptake rates, calculated by CLM3, as sink or source terms, while CLM3 uses saturation and capillary pressure profiles, calculated by TOUGH2, as state variables.
Simulation results show (see Fig. 7 ) that the coupled model greatly improves the prediction of water table elevation, evapotranspiration, surface temperature, and soil moisture, as evaluated using 18 yr of data observed at a real site (Pan et al., 2008). The new model can be extended to include an atmospheric simulation model to simulate hydraulic processes from the top of the atmosphere to deep in the ground.
|
|
|
Conclusions |
|---|
|
TOPABSTRACTINTRODUCTIONVadose Zone Modeling ChallengesSelect Capabilities and...ConclusionsREFERENCES |
|---|
|
ACKNOWLEDGMENTS |
|---|
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONVadose Zone Modeling ChallengesSelect Capabilities and...ConclusionsREFERENCES |
|---|
This article has been cited by other articles:
|
|
 |
|
  J. Simunek and S. A. Bradford Vadose Zone Modeling: Introduction and Importance Vadose Zone J., May 27, 2008; 7(2): 581 - 586. [Full Text] [PDF]
|
|
|
|
|
|
S. Panday and P. S. Huyakorn MODFLOW SURFACT: A State-of-the-Art Use of Vadose Zone Flow and Transport Equations and Numerical Techniques for Environmental Evaluations Vadose Zone J., May 27, 2008; 7(2): 610 - 631. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
R. W. Healy Simulating Water, Solute, and Heat Transport in the Subsurface with the VS2DI Software Package Vadose Zone J., May 27, 2008; 7(2): 632 - 639. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| The SCI Journals | Agronomy Journal | Crop Science | |||
| Journal of Natural Resources and Life Sciences Education |
Soil Science Society of America Journal | ||||
| Journal of Plant Registrations | Journal of Environmental Quality |
The Plant Genome | |||
