| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
Earth Sciences Division, Lawrence Berkeley National Laboratory, University of California, Berkeley, CA 94720
* Corresponding author (K_Pruess{at}lbl.gov)
Received 8 August 2003.
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONHISTORICAL REMARKSGOVERNING EQUATIONSSETTING UP AND SOLVING...INTEGRAL FINITE DIFFERENCESCODE ARCHITECTUREAPPLICATIONSCONCLUSIONSREFERENCES |
|---|
Abbreviations: EOS, equation of state • ESTSC, Energy Science and Technology Software Center • IFD, integral finite differences • LBNL, Lawrence Berkeley National Laboratory • NAPL, nonaqueous phase liquid • NCG, noncondensible gas • PDE, partial differential equation • VOC, volatile organic compound
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONHISTORICAL REMARKSGOVERNING EQUATIONSSETTING UP AND SOLVING...INTEGRAL FINITE DIFFERENCESCODE ARCHITECTUREAPPLICATIONSCONCLUSIONSREFERENCES |
|---|
|
|
HISTORICAL REMARKS |
|---|
|
TOPABSTRACTINTRODUCTIONHISTORICAL REMARKSGOVERNING EQUATIONSSETTING UP AND SOLVING...INTEGRAL FINITE DIFFERENCESCODE ARCHITECTUREAPPLICATIONSCONCLUSIONSREFERENCES |
|---|
Since its beginnings in the early 1980s, the development of TOUGH2 was driven by a desire to model specific types of flow systems. The early development was focused on geothermal reservoir dynamics. Among the important issues for geothermal reservoir modeling are the nonisothermal nature of flow, the importance of phase change (boiling and condensation), and the highly nonlinear nature of two-phase (water-steam) flow. The first functional version of MULKOM was a single-porosity simulator that solved a mass balance for water and an energy balance—noncondensible gases (NCGs) or dissolved solids were not included. In geothermal reservoir problems the coupling between the mass and energy balance equations can be very strong, severely limiting the time step for which a sequential iteration procedure will converge. For example, for a problem of cold water injection into a vapor-dominated reservoir that would entail rapid vaporization with strong latent heat effects, it was estimated that a sequential solution of mass and energy balance equations would converge only for time steps of up to a few hours (Pruess et al., 1983). Accordingly, it was decided to implement a fully simultaneous solution of mass and energy balances and to use fully implicit time stepping to overcome impractical time step limitations. Phase change and two-phase flow introduce strong nonlinearities, requiring an iterative solution by means of Newtonian iteration. The Jacobian matrices arising in geothermal reservoir problems tend to be nondiagonally dominant and ill-conditioned. The coupling to the energy equation introduces very large off-diagonal elements. This ruled out the iterative matrix methods available at the time, and a direct solution procedure was implemented, using a sparse version of Gaussian elimination (Duff, 1977). Subsequent progress with preconditioned conjugate gradient techniques changed the situation, and the current version of TOUGH2 includes iterative solvers that can handle severely ill-conditioned problems (Moridis and Pruess, 1998).
Geofluids typically include NCGs and dissolved solids, primarily CO2 and NaCl. The needs of geothermal reservoir modeling naturally led to the development of fluid property modules for fluid mixtures, with the main focus on CO2 (O'Sullivan et al., 1985). Furthermore, most geothermal reservoirs are fracture-dominated and cannot be adequately described with single-porosity approaches. As will be discussed in more detail below, the IFD technique used for space discretization in the MULKOM and TOUGH codes offers a great deal of flexibility in the geometric description of flow systems. It turned out to be possible to implement double- and multiple-porosity techniques for fractured media simply by preprocessing geometric data, without any coding changes whatsoever (Pruess and Narasimhan, 1982, 1985).
During the 1980s, MULKOM and TOUGH were applied extensively to geothermal reservoir studies, including natural state modeling (Bodvarsson et al., 1984b), design and analysis of well tests (Bodvarsson et al., 1984a), production and injection problems in producing geothermal fields (Pruess and Bodvarsson, 1984; Bodvarsson et al., 1987), and fundamental studies of geothermal reservoir dynamics (Pruess and Narasimhan, 1982, 1985; Pruess, 1985; O'Sullivan et al., 1985; Pruess et al., 1987). Geothermal applications remain a prominent area for TOUGH2 (O'Sullivan et al., 2001) and are covered in a special issue of the journal Geothermics, to be published in 2004. Another application area that made a large impact on the direction of code development was in geologic disposal of heat-generating nuclear wastes. The fluid flow and heat transfer processes in geologic disposal of nuclear wastes have many similarities to processes in geothermal reservoirs. From a broad perspective, the main differences between the proposed repository at Yucca Mountain and a vapor-dominated geothermal reservoir, such as The Geysers, California, is in the thermodynamic regime. In a Yucca Mountain repository, pressures would remain near ambient while at The Geysers original (pre-exploitation) reservoir pressures are near 3.5 MPa. In addition, air is a major fluid component at Yucca Mountain while at The Geysers it is present only in trace quantities. The issues to be addressed for Yucca Mountain focused attention on the ambient groundwater regime, prompting development of fluid property modules that included air as a component, as well as a module that implements a Richards' equation approximation (passive gas phase; Richards, 1931) for unsaturated flow. Analysis of thermohydrologic conditions near individual waste packages in the nuclear waste disposal problem required resolution of much smaller spatial scales than would normally be addressed in geothermal reservoir studies (Pruess et al., 1990; Pruess, 1997). This in turn placed more emphasis on special effects such as formation dry-out from heating, along with strong vapor pressure lowering effects that arise on the path toward dry-out.
Our involvement in international cooperation projects sustained an interest in nuclear waste disposal in the saturated zone. Issues of corrosive gas generation prompted the development of a fluid property module for mixtures of water and hydrogen. Other issues addressed within a nuclear waste disposal context included brines of variable salinity (Oldenburg and Pruess, 1995), parent–daughter radionuclide chains (Oldenburg and Pruess, 1996), and solute tracers that may volatilize into a gas phase. A capability for Fickian hydrodynamic dispersion in two-dimensional regular grids was also developed (Oldenburg and Pruess, 1993).
Another important impetus for developing new simulation capabilities came from environmental contamination problems, especially those that involve volatile organic chemicals (VOCs), such as solvents and hydrocarbon fuels. This led to the development of capabilities for three-phase flow of water, air, and a NAPL (Falta et al., 1995). Recent advancements in this area include multicomponent NAPLs and multicomponent mixtures of noncondensible gases (Pruess and Battistelli, 2002).
Beginning in the mid 1990s, efforts were made to develop capabilities for reactive chemical transport. This was initially motivated by problems in mining engineering, such as the enrichment of copper protore during weathering processes (Xu et al., 2000), and later focused on chemical issues in geothermal systems (Xu and Pruess, 2001a), and in geologic disposal of greenhouse gases (Xu et al., 2003).
The concepts originally developed for the MULKOM code in the early 1980s have proven sufficiently flexible to accommodate all these enhancements and extensions, and development of new fluid property modules and new physical and chemical process capabilities continues.
|
GOVERNING EQUATIONS |
|---|
|
TOPABSTRACTINTRODUCTIONHISTORICAL REMARKSGOVERNING EQUATIONSSETTING UP AND SOLVING...INTEGRAL FINITE DIFFERENCESCODE ARCHITECTUREAPPLICATIONSCONCLUSIONSREFERENCES |
|---|
 | [1] |
The integration is over an arbitrary subdomain Vn of the flow system under study, which is bounded by the closed surface 
n. The quantity M appearing in the accumulation term (left-hand side) represents mass per volume, with 
= 1,..., NK labeling the components (i.e., water, air, CO2, other NCGs, solutes, volatile organic chemicals). F denotes mass flux (see below), and q denotes sinks and sources. n is a normal vector on surface element dn, pointing inward into Vn. Equation [1] expresses the fact that the rate of change of fluid mass in Vn is equal to the net inflow across the surface of Vn, plus net gain from fluid sources.
The general form of the accumulation term is
 | [2] |
In Eq. [2], the total mass of component is obtained by summing over the fluid phases ß (i.e., liquid, gas, NAPL). 
is porosity, Sß is the saturation of phase ß (i.e., the fraction of pore volume occupied by phase ß), 
ß is the density of phase ß, and Xß is the mass fraction of component present in phase ß. TOUGH2 version 2.0 allows for a more general form of the mass accumulation term that includes equilibrium sorption on the solid grains (Pruess et al., 1999). The EWASG module for water–air–salt (NaCl; Table 2) includes solid salt as an active phase (Battistelli et al., 1997).
|
 | [3] |
 | [4] |
Here uß is the Darcy velocity (volume flux) in phase ß, k is absolute permeability, krß is relative permeability to phase ß, µß is viscosity, and
 | [5] |
0). g is the vector of gravitational acceleration. TOUGH2 version 2.0 also considers diffusive fluxes in all phases and includes coupling between diffusion and phase partitioning that can be very important for volatile solutes in multiphase conditions (Pruess, 2002). Diffusive flux of component in phase ß is given by
 | [6] |
0ß is the tortuosity, which includes a porous medium dependent factor 0 and a coefficient that depends on phase saturation Sß, ß = ß(Sß), and dß is the diffusion coefficient of component in bulk fluid phase ß. Special TOUGH2 versions that include a conventional Fickian model for hydrodynamic dispersion have also been developed (Oldenburg and Pruess, 1993, 1995; Wu and Pruess, 2000). Multiphase flows often involve strong heat transfer effects, requiring consideration of an energy balance along with mass balances. Equations similar to the foregoing can be formulated for energy conservation and are given in the TOUGH2 user's guide (Pruess et al., 1999).
By applying Gauss' divergence theorem, Eq. [1] can be converted into the following partial differential equation (PDE):
 | [7] |
 | [8] |
Neglecting variations in liquid phase density and viscosity, as is appropriate for (nearly) isothermal conditions, Eq. [8] simplifies to Richards' (1931) equation:
 | [9] |
= Sl is specific volumetric moisture content, K = kkrllg/µ1 is effective hydraulic conductivity, and h = z + P1/lg is the hydraulic head. Equation [9] is far easier to solve numerically than a full multiphase treatment. Because of its simplicity and importance for modeling saturated–unsaturated flows, a special fluid property module has been developed for TOUGH2 that implements Eq. [9]. |
SETTING UP AND SOLVING DISCRETIZED EQUATIONS |
|---|
|
TOPABSTRACTINTRODUCTIONHISTORICAL REMARKSGOVERNING EQUATIONSSETTING UP AND SOLVING...INTEGRAL FINITE DIFFERENCESCODE ARCHITECTUREAPPLICATIONSCONCLUSIONSREFERENCES |
|---|
 | [10] |
 | [11] |
Here we have written the surface integral as a sum of integrals over surface segments Anm between volume elements (or grid blocks) Vn and Vm, with Fnm the average flux over a surface segment. The discretized flux corresponding to the basic Darcy flux term, Eq. [4], is expressed in terms of parameters for volume elements Vn and Vm as follows:
 | [12] |
|
 | [13] |
Time is discretized as a first-order finite difference, and the flux and sink and source terms on the right-hand side of Eq. [13] are evaluated at the new time level, tk+l = tk + 
t, to obtain the numerical stability needed for an efficient calculation of multiphase flow. This treatment of flux terms is known as "fully implicit" because the fluxes are expressed in terms of the unknown thermodynamic parameters at time level tk+l, so that these unknowns are only implicitly defined in the resulting equations (e.g., Peaceman, 1977). The time discretization results in the following set of coupled nonlinear, algebraic equations:
 | [14] |
,k+1n. For each volume element (grid block) Vn, there are NEQ equations ( = 1, 2,..., NEQ), where usually NEQ = NK + 1, with NK the number of mass components, and an additional equation set up for energy balance. For a flow system with NEL grid blocks, Eq. [14] represents a total of NEL x NEQ coupled nonlinear equations. The unknowns are the NEL x NEQ independent primary variables {xi; i = 1,..., NEL x NEQ} that completely define the thermodynamic state of the flow system at time level tk+1. To apply Eq. [14] to a specific flow system, it is necessary to express all parameters appearing in these equations as functions of the primary thermodynamic variables. This involves (i) PVT relations for fluid properties, and (ii) constitutive relations for the interaction between fluids and the permeable medium. Different equation of state (EOS) formulations have been developed for different fluid mixtures. Table 2 lists the fluid property modules that are presently included in TOUGH2 version 2.0. Additional modules have been developed and may be made available to the public in the future. The most important fluid in the subsurface is of course water, and its properties are calculated in the TOUGH codes within experimental accuracy from the steam table equations as given by the International Formulation Committee (IFC, 1967). A library of user-selectable relative permeability and capillary pressure formulations is available in the code.
The nature of the primary thermodynamic variables depends on the phase state of the fluid system. For single-phase, single-component water, natural primary variables are pressure P and temperature T. In two-phase conditions, however, P and T are no longer independent, because P is equal to the saturated vapor pressure, P = Psat(T). Two-phase water-steam mixtures may be described with (P, S) as primary variables, where S is vapor saturation. In the TOUGH codes we employ "variable switching" to treat phase transitions. For a single-component water system, if a grid block is in single-phase liquid conditions, we monitor fluid pressure P and compare it with the saturated vapor pressure. If P 
Psat single-phase liquid conditions are maintained, whereas when P < Psat we make a transition to two-phase conditions, switching primary variables to (P, S), and initializing S with a small non-zero number, S = 10–6. In a two-component water–air system as in the unsaturated zone, another primary variable is needed in addition to (P, T). In single-phase liquid conditions the additional primary variable may be chosen as air mass fraction X. To test for the possibility that a gas phase may evolve we compute the "bubble pressure" Pbub = Psat + Pair and compare it with total fluid pressure P. When Pbub > P a gas phase will evolve, and primary variable X is replaced by gas saturation S. In flow systems with spatially variable phase composition, we will in general have different kinds of primary variables in different grid blocks. Our experience has shown variable switching to be a very robust and efficient technique for treating phase appearance and disappearance. Because phase transitions are associated with highly nonlinear behavior, they often limit attainable time step sizes.
Equation [14] is solved by Newton–Raphson iteration, which is implemented as follows. We introduce an iteration index p and expand the residuals R,k+1n in Eq. [14] at iteration step p + 1 in a Taylor series in terms of those at index p.
 | [15] |
Retaining only terms up to first order, we obtain a set of NEL x NEQ linear equations for the increments (xi,p+1 – xi,p):
 | [16] |
All terms 
Rn/xi in the Jacobian matrix are evaluated by numerical differentiation to achieve maximum flexibility in the manner in which various terms in the governing equations may depend on the primary thermodynamic variables. Equation [16] is solved by sparse direct matrix methods (Duff, 1977) or iteratively by means of preconditioned conjugate gradients (Moridis and Pruess, 1998). Iteration is continued until all residuals are reduced below a preset convergence tolerance, |R,k+1n/M,k+1n| 
, typically chosen as = 10–5. We note that some numerical simulation codes monitor the maximum changes in primary variables from one Newtonian iteration to the next and diagnose convergence based on these maximum changes being "sufficiently" small. Such a procedure is hazardous in that changes in primary variables may be small not necessarily because the governing equations are satisfied within small errors, but because convergence rates are poor. The residual-based criterion used in the TOUGH codes avoids such "false convergence."
|
INTEGRAL FINITE DIFFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONHISTORICAL REMARKSGOVERNING EQUATIONSSETTING UP AND SOLVING...INTEGRAL FINITE DIFFERENCESCODE ARCHITECTUREAPPLICATIONSCONCLUSIONSREFERENCES |
|---|
Integral finite differences encourages a "physical" view of model building, analogous to assembling a laboratory experiment. It provides a very simple conceptual basis for assigning boundary conditions. A flow system can be viewed as a network of boxes that exchange mass and energy. If temperature is to be held constant at a bounding surface, all that is required is defining a grid block on the "other" side of the surface, with a "small" nodal distance, giving it a very large heat capacity (e.g., by assigning a very large rock grain density or a very large volume), and assigning it the desired boundary temperature. Such a grid block would be treated like any other block during a flow simulation. Heat and mass fluxes entering or leaving that block will not be able to effect a noticeable change in temperature because of the very large heat capacity. Boundary conditions of constant pressure (or moisture tension, or solute concentration) can be realized in analogous fashion. No-flow boundaries are modeled in IFD with simplicity itself, by simply not including any grid blocks beyond the no-flow boundary. More general and time-dependent boundary conditions can also be modeled with TOUGH2.
The IFD technique for space discretization has a number of advantages over conventional finite difference discretizations. All geometric quantities are defined locally, and the entire geometric information for the flow problem is provided through a list of volumes Vn, interface areas Anm, nodal distances Dnm, and components gnm of gravity acceleration in the direction of the line connecting the nodal points. There is no need to make reference to a global system of coordinates, and unstructured grids may be used. This provides great flexibility when dealing with spatially irregular features. One-, two-, and three-dimensional flow systems with regular or irregular grids can all be treated on the same footing. More advanced discretization approaches for fractured and highly heterogeneous media, such as double-porosity (Barenblatt et al., 1960), multiple interacting continua (Pruess and Narasimhan, 1985), and multiregion models (Gwo et al., 1996) can be implemented simply by appropriate geometric preprocessing. Higher-order discretization methods, such as nine-point differencing (Yanosik and McCracken, 1979), can also be realized by means of appropriate geometric preprocessing, without any changes in the simulation code (Pruess, 1991b). This flexibility comes at no cost to the user because for regular grids referred to a global coordinate system, IFD is mathematically equivalent to conventional finite differences (Moridis and Pruess, 1992).
Another important advantage of the IFD method is that it facilitates writing readable code. Quantities with many indices can be avoided, and variable names used in TOUGH2 correspond closely to the theoretical formulation given above.
As introduced in the previous section, the IFD technique for space discretization is very general but also entirely formal, in the sense that it does not provide specific guidance as to how discrete subdomains (grid blocks) should be defined. The fundamental IFD equations (Eq. [13] and [14]) are always valid, regardless of space discretizations; however, these equations will be practically useful only for such space discretizations for which an approximate evaluation of interface fluxes can be made from volume averages in the grid blocks involved (Eq. [12]). The IFD technique itself does not offer "intrinsic" means to estimate and minimize space discretization errors. For regular grid systems, such errors can be estimated from Taylor series expansions as used in conventional finite differences. In more general cases, physical arguments can be used to judge how well the pressure gradient driving flux across an interface can be approximated in terms of average pressures in the participating grid blocks. In isotropic media, it is the normal component of pressure gradient at an interface that drives flow between the two grid blocks. To be able to obtain this normal component from information for just two grid blocks, it is necessary that the interface be perpendicular to the nodal line connecting grid block centers. An important type of discretization that meets this constraint is (generally unstructured) Voronoi grids, where interfaces are constructed as the perpendicular bisectors of nodal lines.
As currently implemented in the TOUGH codes, a flow connection is an ordered pair of two volume elements. This is very simple and convenient from a user's viewpoint, but it provides only a single component of pressure (or other) gradients, along the nodal line. There is nothing intrinsic to IFD that would preclude providing more geometric information locally, so all components of gradients at the interface could be obtained. For example, the Fickian model for hydrodynamic dispersion is anisotropic, and requires the full concentration gradient vector at the interface, not just the normal component. For regular grid systems the full gradient can be calculated in a straightforward manner by using information from additional neighboring grid blocks (Oldenburg and Pruess, 1993, 1995). Similarly, higher-order techniques for more accurate representation of advective terms can be easily implemented in structured grids (Oldenburg and Pruess, 2000). Although ordered-pair connections allow the use of unstructured Voronoi grids, we believe that the IFD is far more powerful and flexible than has been realized in implementations to date. It should be possible to define "multiblock connections" that, in addition to geometry data of the two grid blocks between which flow is occurring, include data of additional neighboring grid blocks to allow an accurate estimation of the full pressure gradient (and other driving forces) at the interface. These more general connections should be able to deal with irregular geometric features, such as fault offsets, interfaces that are not perpendicular to nodal lines, and locally refined grids.
|
CODE ARCHITECTURE |
|---|
|
TOPABSTRACTINTRODUCTIONHISTORICAL REMARKSGOVERNING EQUATIONSSETTING UP AND SOLVING...INTEGRAL FINITE DIFFERENCESCODE ARCHITECTUREAPPLICATIONSCONCLUSIONSREFERENCES |
|---|
|
Two ancillary arrays are in use. DX holds latest (Newtonian iteration) increments to X, and DELX holds "small" increments of X, typically of order 10–8 relative to primary variables, which are used to calculate numerical derivatives of all terms in Eq. [16] with respect to the primary variables. These arrays need to be available throughout many of the calculations performed in TOUGH2. They are placed in COMMON blocks, and the action of subprograms is viewed as operating on these COMMON blocks.
The current version of TOUGH2 is written in FORTRAN 77, predating dynamic memory allocation. To achieve convenient scalability (flexibility in adjusting memory allocations to different problem sizes), the arrays X and PAR, as well as other problem-size dependent arrays, are each placed in a separate labeled COMMON block and are given actual dimension only in an easily modifiable PARAMETER statement in the main program, while they are dimensioned X(1) etc. in subroutines. The more recent TOUGH2 version 2.0 also uses INCLUDE files to facilitate flexible memory assignments.
An important design goal has been readability of the code. When reviewing a specific code section or subprogram, it should generally be transparent to the user what is being done, so the code may be easily modified to suit particular application needs.
|
APPLICATIONS |
|---|
|
TOPABSTRACTINTRODUCTIONHISTORICAL REMARKSGOVERNING EQUATIONSSETTING UP AND SOLVING...INTEGRAL FINITE DIFFERENCESCODE ARCHITECTUREAPPLICATIONSCONCLUSIONSREFERENCES |
|---|
Most applications of the TOUGH codes are currently being run on Unix workstations and on PCs. The spatial scale of flow systems simulated with TOUGH2 ranges from microscopic (Webb and Ho, 1998) to basin-scale. Time scales on which flow processes have been modeled range from a fraction of a second to geologic time. Three-dimensional problems with a few thousand grid blocks are considered modest size with current computing platforms. For the nuclear waste storage investigations at Yucca Mountain, the LBNL group is routinely running three-dimensional problems with more than 100000 grid blocks on PCs. A massively parallel version of TOUGH2 has also been developed and has been used for problems with more than 2 million grid blocks (Wu et al., 2002; Zhang et al., 2003).
Special fluid property modules that are not currently part of the publicly distributed version of TOUGH2 include capabilities for supercritical temperatures (Brikowski, 2001), freezing water and hydrate formation (Moridis, 2002), multiple tracers and colloids (Moridis et al., 1999), and four-phase mixtures of water, liquid and gaseous CO2, and solid salt (Pruess, 2003). Modules for flow of foam and other non-Newtonian fluids have also been developed.
|
CONCLUSIONS |
|---|
|
TOPABSTRACTINTRODUCTIONHISTORICAL REMARKSGOVERNING EQUATIONSSETTING UP AND SOLVING...INTEGRAL FINITE DIFFERENCESCODE ARCHITECTUREAPPLICATIONSCONCLUSIONSREFERENCES |
|---|
Most of the development of the TOUGH codes was performed at LBNL. Important contributions were made by authors from outside Berkeley. The development of the various TOUGH codes at LBNL was generally not performed as software development projects, but rather was driven by a desire to obtain an understanding of different kinds of flow systems. Whenever simulation capabilities evolved to the point where they were believed to be sufficiently accurate, efficient, and robust to be useful to others, then we would clean up and document the codes, and release them to the public.
There is no finite process by which all bugs that may be present in a large simulation program can be identified and corrected. Consequently, a user should expect that there may be bugs in his/her application. It is the user's responsibility to check the numbers. Multiphase flow systems include an infinite variety of features. It is not practically possible to exhaustively cross-check the mutual compatibility and proper performance of all simulator options. Fixing a bug may cause unanticipated problems elsewhere. Continuing vigilance and application testing are needed. From a practical viewpoint, the best approach for identifying and correcting coding errors is through extensive applications to a wide range of diverse problems, by a broad community of users and analysts. This has been one of the motivating forces behind our drive for public release and wide distribution of our simulation codes. The TOUGH2 homepage on the internet (http://www-esd.lbl.gov/TOUGH2/) includes a section with reported bugs in the version 2.0 code and suggested fixes, testifying to the value of an interactive user community for making a "better code."
The LBNL group has always offered free technical support to users of the TOUGH codes, but being a research and development organization, the time and resources we can devote to such support are limited. Public release of the TOUGH codes and their use by others are workable in practice only if users' needs for technical support are generally modest. We achieve this by releasing codes only after they are thoroughly tested, debugged, and documented. This involves a period of ß-testing by users outside the development team, so that errors or flaws in coding and documentation can be identified and corrected. Care is taken to maintain a continuous chain of code testing and verification, where newer codes are run and rerun on test problems used in older versions (Pruess et al., 1996). Adequate internal as well as external documentation is crucial for efficient code maintenance. Worked sample problems are an important aspect of this. They aid in code installation and debugging by providing benchmarks, and they serve as templates to help jump-start new applications.
|
ACKNOWLEDGMENTS |
|---|
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONHISTORICAL REMARKSGOVERNING EQUATIONSSETTING UP AND SOLVING...INTEGRAL FINITE DIFFERENCESCODE ARCHITECTUREAPPLICATIONSCONCLUSIONSREFERENCES |
|---|
This article has been cited by other articles:
|
|
 |
|
  L. Bolshov, P. Kondratenko, K. Pruess, and V. Semenov Nonclassical Transport Processes in Geologic Media: Review of Field and Laboratory Observations and Basic Physical Concepts Vadose Zone J., November 26, 2008; 7(4): 1135 - 1144. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
S. Finsterle, C. Doughty, M.B. Kowalsky, G. J. Moridis, L. Pan, T. Xu, Y. Zhang, and K. Pruess Advanced Vadose Zone Simulations Using TOUGH Vadose Zone J., May 27, 2008; 7(2): 601 - 609. [Abstract] [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]
|
|
|
|
|
|
H.-H. Liu and T. H. Illangasekare Preface: Recent Advances in Modeling Multiphase Flow and Transport with the TOUGH Family of Codes Vadose Zone J., February 25, 2008; 7(1): 284 - 286. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 N. K. C. Twarakavi, H. Saito, J. Simunek, and M. Th. van Genuchten A New Approach to Estimate Soil Hydraulic Parameters Using Only Soil Water Retention Data Soil Sci. Soc. Am. J., February 15, 2008; 72(2): 471 - 479. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
A. Boivin, J. Simunek, M. Schiavon, and M. Th. van Genuchten Comparison of Pesticide Transport Processes in Three Tile-Drained Field Soils Using HYDRUS-2D Vadose Zone J., June 21, 2006; 5(3): 838 - 849. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||
