| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
a Soil and Water Science Department, University of Florida, 2169 McCarty Hall, Gainesville, FL 32611-0290
b USDA-ARS-NA-GPSR, 301 South Howes Street, Fort Collins, CO 80521
* Corresponding author (rsm{at}mail.ifas.ufl.edu)
Received 26 December 2001.
 |
ABSTRACT |
|---|
 TOPNOTES ABSTRACT INTRODUCTIONPertinent Features of...Adaptive Grid Refinement MethodsWater Flow in Porous...Solute Transport in Porous...Examples of Local Adaptive...CONCLUSIONSREFERENCES |
|---|
Abbreviations: CPU, central processing unit • CVFEM, control volume finite element numerical method • FCT, flux-corrected transport • FDM, finite difference method • FEM, finite element method • LAGR, local adaptive grid refinement • LEM, Lagrangian-Eulerian method • LUGR, local-uniform-grid refinement • MMOC, modified method of characteristics • MSRPT, modified single-step reverse particle tracking • PDE, partial differential equation • SRPT, single-step reverse particle tracking
|
INTRODUCTION |
|---|
|
TOPNOTESABSTRACTINTRODUCTIONPertinent Features of...Adaptive Grid Refinement MethodsWater Flow in Porous...Solute Transport in Porous...Examples of Local Adaptive...CONCLUSIONSREFERENCES |
|---|
A number of specific problems of water flow and solute transport in soils provide cases where pregenerated, fixed grids may provide inaccurate numerical results. Critical problems such as penetration of sharp wetting fronts during infiltration in initially dry soil or sharp concentration fronts during convection-dominant solute transport impose local spatial resolution requirements that typically change with time and space. Although a fixed grid with a fine mesh may provide accurate numerical solutions for such problems, high computational cost may be prohibitive for complex multidimensional simulations of large-scale problems. One such problem involves the challenge of modeling variably saturated water flow in highly heterogeneous soils, particularly at relatively large hydrologic scales.
Three types of errors are associated with any numerical solution of partial differential equations: (i) round off, (ii) truncation (or discretization error), and (iii) inherited errors (Adey and Brebbia, 1983). Round off errors occur as a consequence of finite precision arithmetic in numerical calculations. These errors are normally random and generally negligible in high precision computers. Truncation errors occur due to truncation of the specific expansion used. Assessment of the order of the error is extremely difficult for expansions other than the finite difference type. Inherited errors represent accumulation of total errors from previous steps. Determination of errors is critical since numerical simulations rather than exact solutions are involved (Yeh, 2000).
"Scientific and engineering computation has become so complex that traditional numerical computation on uniform meshes is generally not possible or too expensive" (Bern et al., 1999). Although faster computers have emerged in recent years, faster computers tend to be used to solve even more difficult problems (Finlayson, 1992). Therefore, developing efficient numerical methods that accurately solve a problem with a minimum of computer time continues to provide a worthy challenge. This is particularly important for three-dimensional simulations of large geographical areas of the vadose zone and associated groundwater zone.
Bern et al. (1999) emphasized that numerical solutions should utilize computable estimates of discretization errors to dynamically assign mesh grids. Local adaptive grid refinement algorithms provide a powerful means to accomplish that purpose through the use of automated grid assignment based upon error estimates. Incorporation of LAGR algorithms into numerical models offers potential for accuracy enhancement of numerical approximations for such cases. Automatic adaptive grid refinement procedures in LAGR algorithms provide enrichment of an initial mesh with the goal of providing solutions with prescribed accuracy specifications in an optimal manner. Local adaptive grid refinement algorithms provide a gain in efficiency by the use of fewer nodes (Finlayson, 1992).
Adaptive grid refinement is one of three broad categories of approaches reported in the literature to ensure accurate numerical solutions of PDEs for critical cases of water flow and solute transport in porous media. These categories are (i) mathematical alteration of governing PDEs, (ii) incorporation of adaptive grid refinement or LAGR algorithms, and (iii) combinations of the first two categories. Each category provides approaches designed to improve inadequate matches between the mathematics given by the governing PDEs and the numerical approximations presented in model computer codes. Representative examples of approaches in the first category for water flow include incorporation of mathematical transformations (Kirkland et al., 1992; Johnsen et al., 1995; Pan and Wierenga, 1995) and interpolation (Pan et al., 1996). Eulerian (Chang and Slattery, 1990) and combined Eulerian-Lagrangian (Zhang et al., 1993) approaches in this category have been reported for solute transport. Mesh refinement (Clausnitzer et al., 1998) and moving mesh (Dane and Mathis, 1981) LAGR approaches have been reported for water flow. For solute transport mesh refinement (Trompert, 1993; Wolfsberg and Freyberg, 1994), Lagrangian (Zegeling et al., 1992), and combined Eulerian-Lagrangian (Yeh, 2000) LAGR approaches have been reported. Cao and Kitanidis (1999) present an approach for multidimensional water flow that falls into the third category. Examples of all three categories are discussed here. Although the effectiveness of many numerical approaches based on the mathematical alteration of governing PDEs (first category) have been clearly demonstrated for many problems, more complex problems often benefit from the other two approaches. Thus, the primary emphasis in this review is to critically examine benefits and limitations for using LAGR algorithms (numerical approaches including incorporation of adaptive grid refinement or LAGR algorithms and combinations of the first two categories) in both FDM and FEM numerical models for critical cases of flow or transport in soils.
|
Pertinent Features of Traditional Numerical Methods |
|---|
|
TOPNOTESABSTRACTINTRODUCTIONPertinent Features of...Adaptive Grid Refinement MethodsWater Flow in Porous...Solute Transport in Porous...Examples of Local Adaptive...CONCLUSIONSREFERENCES |
|---|
Finite difference methods and finite element methods are two common approaches for transforming a set of PDEs that constitute a mathematical model into a set of algebraic equations that constitute a discrete model (Wang and Anderson, 1982). The flow domain is discretized into cells to provide FDM spatial grids and elements to provide FEM spatial grids. An approximate solution for the mathematical problem is then obtained by use of iterative techniques or direct matrix methods to solve the set of algebraic equations.
Accurate and reasonable approximations to the solution of the governing PDEs require three conditions: (i) convergence of the approximating solutions to that of the PDEs, (ii) stable decay of errors in arithmetic operations, and (iii) compatibility between the differential and approximating equations (Yeh, 1999). Differences between the approximating and differential equations is designated as the truncation error. Although compatibility and stability often imply convergence, this relationship has only been established for a limited number of specific differential equations (Yeh, 1999).
The FEM and FDM numerical approaches differ in a philosophical sense. In FDM, a value for the unknown variable (e.g., pressure head h in water flow) is computed at a node, which also is the average for the cell that surrounds the node; whereas, in FEM, variation in the variable within an element is precisely defined using linear or nonlinear interpolation functions (Anderson and Woessner, 1992). Positive features of FDM relative to FEM include being simpler to understand and program and generally requiring fewer input data. Mesh-centered FDM grids are convenient for problems where unknown parameters such as water pressure head h and solute concentration C are specified; whereas, block-centered grids are advantageous for problems where water or solute flux is specified across a boundary (Yeh, 1999). The FEM provides a distinct advantage over FDM for approximating irregular geometrical regions. Other positive features for FEM include ready incorporation of heterogeneous and anisotropic porous media properties, better capacity to handle internal boundaries, and providing better simulation of point sources and sinks, seepage faces, and moving water tables (Anderson and Woessner, 1992; Friedel, 2001).
Taylor series expansion in FDM is commonly used to provide difference approximation for derivatives in differential equations. Implicit (i.e., backward-difference scheme) and explicit (i.e., forward-difference scheme) finite difference approximations are common for spatial derivatives. For example, a weighted average of approximations between current tn = n
t and advanced times tn+1 = [n + 1]t (where t is the time step) is given for the second partial derivative of water pressure head, h, that occurs in the one-dimensional water flow equation
 | [1] |
is a weighting parameter such that 0 
1 (Wang and Anderson, 1982). Subscripts and superscripts designate the spatial location of nodes and time levels, respectively. For = 1, the spatial derivative is approximated solely at the advanced time level tn+1 to give a fully implicit scheme which is unconditionally stable; that is, the time step t can be chosen independently of node spacing x. For = 0, the spatial derivative is approximated solely at the old or current time level tn to give a fully explicit scheme that tends not to be computationally efficient due to required very small values for both x and t. For = 1/2, the spatial derivative occurs midway between times tn and tn+1 to give the Crank–Nicolson method. The Crank–Nicolson method is both computationally efficient and unconditionally stable. Generally for the Crank–Nicolson method, the smaller the truncation error, the faster is the convergence of difference equations to the differential equation (Remson et al., 1971). |
Adaptive Grid Refinement Methods |
|---|
|
TOPNOTESABSTRACTINTRODUCTIONPertinent Features of...Adaptive Grid Refinement MethodsWater Flow in Porous...Solute Transport in Porous...Examples of Local Adaptive...CONCLUSIONSREFERENCES |
|---|
In an adaptive grid, physics of a problem at hand must ultimately direct the grid points to congregate so that a functional relationship on these points can represent the physical solution with sufficient accuracy (Thompson, 1985). The mathematics controls the points by sensing gradients in the evolving physical solution, evaluating accuracy of discrete representation of the solution, communicating needs of the physics to the points, and providing mutual communication among points as they respond to the physics.
Two important target capabilities of adaptive grid refinement strategies are to optimize overall computational effort (i.e., providing the best possible results for a fixed computational effort) and to provide indicators of residual error as a measure of simulation reliability (Oden, 1989). Adaptivity involves automated "restructuring" a numerical scheme to improve the quality of the resolution (Oden and Demkowicz, 1987). Such restructuring may include changing the number of cells or elements, increasing the local order of the approximation, moving nodal points, changing algorithm structure during the solution of governing PDEs on a changing numerical model. Effective LAGR methods should require minimal mesh structure (i.e., essentially coordinate-independent), demonstrate stability during changes and distortions of the mesh, and utilize parallelism or vectorization to deal with large data management requirements in an efficient way (Oden and Demkowicz, 1987). The LAGR finite element methods meet these requirements as effective adaptive schemes for very complex three-dimensional problems (Oden and Demkowicz, 1987).
Adaptive grid methods can be simply defined as "numerical schemes which automatically adjust themselves to improve solutions" (Oden and Demkowicz, 1987). Successful numerical models with adaptive mesh schemes contain three components: a flow solver, a strategy for identifying regions for refinement and coarsening, and a mechanism for dynamically altering the mesh. A posteriori local error estimates are commonly used to assess the quality of a numerical solution after an initial calculation on a trial mesh is made. Calculation of very precise and sharp error estimates provides the only general method for assessing the actual quality of a numerical solution (Oden and Demkowicz, 1987). First and second derivatives of domain variables are often used as simple means to identify high-error regions for mesh adaption. The trial solution is then used to compute an indication of the distribution of error and then to change the approximation to reduce the error. Numerical schemes are restructured so as to optimize resolution.
Local adaptive grid resolution approaches typically include three major tasks: grid maintenance, integration, and communication (Diaz et al., 1989). Grid maintenance determines placement and/or removal of high resolution grids. Placement or removal of high resolution grids is very important for modeling problems having both global and localized phenomena. Ultimately, the solution obtained with such numerical methods is as good as the method's ability to track localized phenomena over the domain.
An adaptive grid approach is most effective when it is dynamically coupled with the solution, so that the solution and the grid are solved together in a single continuous problem (Thompson, 1985). Points must not move independently, but rather, each point much somehow be coupled at least to its neighbors.
A number of major computational challenges relate to the use of LAGR algorithms, however. These challenges include:
Parallel computation procedures become important when models include such complexities as nonlinearity, multidimensionality, and transcience.
|
Water Flow in Porous Media: Special Resolution Requirements |
|---|
|
TOPNOTESABSTRACTINTRODUCTIONPertinent Features of...Adaptive Grid Refinement MethodsWater Flow in Porous...Solute Transport in Porous...Examples of Local Adaptive...CONCLUSIONSREFERENCES |
|---|
The continuity equation for three-dimensional transient water flow in unsaturated subsurface media is typically expressed in three forms:
 | [2] |
 | [3] |
(iii) and the mixed-form (i.e., the Richards equation)
 | [4] |
is volumetric water content, h is water pressure head, K is soil hydraulic conductivity, D() = K()/C() is soil water diffusivity, and 
(h) = d/dh is the specific water-holding capacity (Celia et al., 1990). Numerical solutions for h-based forms can be used for both saturated and unsaturated conditions (Kirkland et al., 1992). The h-based form is sometimes associated with poor accuracy (e.g., mass balance error or numerical oscillations) due primarily to implementation of a capacity term in a numerical scheme that may be inconsistent with that given in the PDE (Eq. [4]) for highly nonlinear cases. Such accuracy problems may occur even with the use of iteration for numerical solutions of the h-based form. The -based form degenerates under fully saturated conditions and has restricted applicability in heterogeneous subsurface porous media because material discontinuities produce discontinuous profiles. The mixed form of the continuity equation in Eq. [4] is generally considered superior to the other two forms because of robustness with respect to mass balance. Note that taking water content as a secondary dependent variable in the mixed-form (Eq. [4]) such that (h) effectively converts this PDE to an h-based form. Although mass balance is often used as an indicator for the success of numerical solutions of Richards' equation, Warrick (1991) showed that mass balance alone does not always assure an acceptable solution. He demonstrated that finite difference approximations to Darcy flow may lead to erroneous water velocities for unsaturated conditions, particularly for flow across large pressure gradients. The choice of appropriate intrablock approximations for unsaturated hydraulic conductivity was shown to be a critical step in estimating Darcy velocity during one-dimensional simulations of constant flux infiltration into a very dry clay loam soil. Simulations with effective conductivities based upon traditional arithmetic and geometric averages provided relatively diffuse and sharp wetting fronts, respectively, in comparison with that for a quasi-analytical solution. A weighted average provided intermediate results. In spite of these differences in simulated water content profiles, mass was conserved for all cases.
The three differential forms of the water continuity equation (Eq. [2], [3], and [4]) for water flow provide mixed hyperbolic–elliptic behavior under conditions of variable water saturation of the soil. Hyperbolic behavior occurs during unsaturation, and elliptic behavior occurs during saturation.
An example of mathematical alteration of governing PDEs [category (i) numerical approaches] was presented by Kirkland et al. (1992) where two numerical algorithms were implemented which retain advantages of -based finite difference methods for fully unsaturated heterogeneous problems but still permit fully saturated conditions. One algorithm was based upon transforming Richards' equation with a variable parameter 
that has characteristics of water content when the soil is unsaturated and of pressure head h when the soil is at or near saturation. The other algorithm was an h-based method (Eq. [3]), which uses flux updating to achieve mass balance by updating for unsaturated nodes using
 | [5] |
 | [6] |
-based algorithm, has the advantage of significantly improving performance relative to conventional h-based algorithms. It was also reported to be very easy to incorporate into existing h-based codes. Water fluxes were reported to be more sensitive to spatial step size in solving Richards' equation than were water contents or pressure heads. They suggest that the use of water flux rather than pressure head may be preferable in judging convergence if solute transport is to be included with the water flow model.
A further example of category (i) (mathematical alteration of governing PDEs) approaches where mathematical transformations have been reported for enhanced numerical solution of water continuity equations can be found in Pan and Wierenga (1995). A fast, numerically robust numerical scheme utilizing PDE transformation was introduced for solving the h-based equation for variably saturated, heterogeneous media. (Pan and Wierenga, 1995). The one-dimensional form of Eq. [3] was transformed using a simple nonlinear transform of h (e.g., pressure Pt)
 | [7] |
-0.04 cm-1) independent of both K(h) and (h) relationships. The resulting three-dimensional form of Richards' equation was
 | [8] |
* is a transformed specific water capacity given as
 | [9] |
 | [10] |
A critical feature of Pt-based Eq. [8] is that near saturation ßh << 1 such that Pt = h for h 
0. For the specific case when h = 0, * = 
= 0 and K* = K = Ksat. Thus, the continuity of both Pt(h) and 
Pt/t is guaranteed at the joint point h = 0 linking conditions of positive and negative soil water pressure head. For large negative values of h, Pt/z < h/z, resulting in faster simulations for vertical flow with less mass-balance error for conditions involving large gradients of h. For zero and positive values of h, Pt/z = h/z and Pt/z > h/z, respectively. The Pt transform is free from difficulties commonly experienced for heterogeneous soils and with hysteresis for the case of the -based continuity equation.
The Pt transform greatly reduces the steepness of wetting fronts in initially dry soil, thus mitigating the need for LAGR. Consequently, the transformed pressure head gradient needed to move water through dry soil was reported to be less, the numerical scheme more robust, and decreased central processing unit (CPU) time. Excellent agreement was observed for simulations of constant-flux infiltration into soil with very dry initial conditions using Pt-based and -based (Kirkland et al., 1992) models, although the Pt-based model used much less CPU time. Smeared wetting fronts resulted for similar simulations using an h-based model strongly indicated the need for local adaptive grid refinement. Incorporation of the Pt-based numerical model into existing h-based codes was reported to be easy.
A general Eulerian FDM with modified Picard iteration for the mixed form of Richards' equation for one-dimensional water flow has been expressed in tridiagonal matrix form as
 | [11] |
t/(z)2 for a fixed spatial grid; = 1 specifies backward FDM, which is fully implicit; = 1/2 specifies the Crank–Nicolson FDM method, which is centrally differenced in space; and = 0 specifies fully explicit or forward FDM (Johnsen et al., 1995). Superscripts n and m denote time and iteration levels, respectively. Celia et al. (1990) used a fully implicit form of Eq. [7] without the source–sink term.
The fully implicit, mixed form of Richards' equation is generally considered to be an accurate model for both saturated and unsaturated flow (Johnsen et al., 1995). The next to last term on the right-hand side of Eq. [11] provides a mass-balance adjustment since the term on the left-hand side of Eq. [4] was approximated using a Taylor series expansion to give
 | [12] |
The discrete analog of (h)h/t in the h-based Eq. [3] is not equivalent to that of /t in Richards' Eq. [4] due to the highly nonlinear nature of the specific capacity term (h), even though they are mathematically equivalent (Pan and Wierenga, 1995).
Conservation of mass alone, however, does not guarantee good numerical solutions for water flow using either FDM or FEM. For the special case of sharp wetting front development during water infiltration into initially dry soil, conventional mass-distributed FEM approximations have been reported to produce oscillatory solutions even as mass is conserved (Celia et al., 1990). Low initial water contents (Celia et al., 1990) and large grid spacing (van Genuchten, 1982) enhance such undesirable oscillation. Pan et al. (1996) provided a physical interpretation of numerical oscillation when FEM is used to model water flow into dry soil. They emphasized that mass conservation in FEM is applied to each element, whereas mass conservation in FDM is applied directly to each cell. Element size represents the FEM spatial resolution because the nodal values represent the final solutions. Detailed distributions of h, , and K within an element can be assumed arbitrarily without violating the mass conservation law at the element level. Linear or nonlinear interpolation functions are commonly used to obtain distributions for these variables. Water storage and flux within each element in FEM are split into several components in the function space, each of which corresponds to one component of the boundary flux of the element. Therefore even with correct application of physical laws at the element level, traditional mass-distribution schemes may generate an incorrect response for a neighboring node due to the highly nonlinear properties for water flow in unsaturated soil causing numerical oscillation (Pan et al., 1996). Although FEM offers the flexibility to arbitrarily choose both element shape and interpolation functions within an element, the cost of that flexibility is that relationships between nodes (at the point level) may be incorrect physically, even though their integral (at the element level) is correct. Mass-lumped schemes have been shown to improve numerical convergence and eliminate oscillations when simulating infiltration into dry soils (Neuman, 1972). However, smearing of the wetting front may occur when a mass-lumped FEM scheme is used for infiltration into a dry soil (Pan et al., 1996).
Pan et al. (1996) introduced two FEM mass-distribution schemes that are free of oscillation and decrease smearing in the vicinity of a sharp wetting front. Freedom from oscillation occurs with those schemes because the distribution of /t correctly expresses the physical relationship between nodes. Interpolation functions utilized are dependent upon the nodal values of pressure head at the previous time. The introduced mass-distributed schemes provided less smearing and smaller global mass balance errors than a traditional mass-lumped scheme. Accuracy improvements of the new schemes were obtained at slightly higher CPU time requirements.
An example of category (ii) (incorporation of adaptive grid refinement or LAGR algorithms) numerical approaches for water flow was presented by Dane and Mathis (1981) where a self- adaptive LAGR algorithm was combined with FDM for the one-dimensional h-based form of the Richards equation (Eq. [3]) to describe transient water flow in unsaturated porous media. A fixed number (N) of spatial grid nodes were dynamically arranged such that fine grids are locally imposed in those flow regions where large changes in total head gradients (
H = h - 1) occur. Water loss as evaporation from the surface of moist soil and water infiltration into initially dry soil profiles tend to create regions with large changes in H that tend to cause problems with numerical convergence in model simulations.
Numerical tests were performed by Dane and Mathis (1981) for constant flux infiltration of water into a vertical sand column with initial, uniform water content and for simultaneous evaporative and drainage loss in uniform coarse porous media with a water table located at 3 m depth. Space nodes were redistributed by the code as a relatively sharp wetting front moved downward. The smallest z values were dynamically placed in the vicinity of the front with the largest values placed both above and below the front. An important advantage of this adaptive method is that no preliminary simulations are necessary to obtain the optimum combination of t and z for a given problem. For evaporation, an initial constant flux boundary condition at the soil surface was changed to a constant head when h reached -15 000 cm of water. This change in the surface boundary condition occurred between 4.97 and 6.17 h, resulting in a decrease in evaporation rate as the water content of the surface soil decreased. For higher evaporation rates occurring at 0.27 and 4.97 h, the smallest z values were dynamically placed in the soil below the soil surface, as well as in the soil in the vicinity of the water table. For transient water flow with locally large changes in H due to infiltration, evaporation, and water uptake by roots, the LAGR approach for water flow offers more accurate and efficient computation of water flux q and volumetric content than nonadaptive schemes. Accurate values for q and provide critical input for accurate simulation of associated solute transport during transient water flow.
A more recent form of the self-adaptive scheme by Dane and Mathis (1981) was incorporated into the WAFLOWM model, where hydraulic head H rather than pressure head h was the dependent variable in the mixed-form of Richards' equation (Johnsen et al., 1995), such that
 | [13] |
 | [14] |
 | [15] |
 | [16] |
Dane et al. (1982) demonstrated the utility for the use of second-order spatial derivative of H as an effective a posteriori error estimator during LAGR-FDM simulation of water infiltration into initially dry soil. Figure 1 was adapted from Dane et al. (1982) and provides a graphic plot of the one-dimensional spatial distribution of simulated soil water content (z) and the corresponding second-order spatial derivative 2H/z2. This error estimator justifies the fine grid automatically selected by the LAGR algorithm for the immediate vicinity of the sharp wetting front.
|
The RZWFLO model, extracted from the Root Zone Water Quality Model (RZWQM) reported by Ahuja et al. (2000) was modified to allow simulation of saturated flow. Within RZWFLO, the Green and Ampt approach was used to describe the water infiltration phase and the FDM given by Eq. [11] was used to solve the water redistribution phase. A 1-cm node spacing was used in RZWFLO during infiltration and variable node spacing up to 5 cm with variable time steps up to 1 h were used during water redistribution. Following rainfall, extra grid points were manually added within regions of steep H for the RZWFLO to prevent numerical convergence problems. Most of the time 34 grid nodes were used in RZWFLO for 15-m tile spacing, although up to 21 extra nodes were temporarily added during times with steep H gradients. The adaptive, optimizing grid scheme in WAFLOWM maintained 50 total nodes throughout simulations, but dense grids were dynamically placed in regions with steep gradients.
Transpiration, T, was calculated as potential evapotranspiration, PET, minus potential evaporation, E. A simple Thornewaite approach was used to approximate PET. T was distributed as a sink–source term over a constant rooting depth of 15 cm. Although results from the two models predicted fluctuations in groundwater table fairly well, estimated depths to water table were generally overpredicted for all 3 yr. The overprediction was attributed to inaccuracies in approximations of both evapotranspiration (ET) and flow rates from tile drains. Simple source–sink terms were used to describe drainage. In most cases, results from WAFLOM more closely approximated observed data than RZWFLO. Both WAFLOM and RZWFLO models maintained perfect soil water mass balance throughout the simulation, thus paralleling observations by Pan et al. (1996).
An adaptive-grid finite FEM solver RES-1D was developed by Clausnitzer et al. (1998) for solving a one-dimensional form of the h-based Richards equation. The adaptive-grid algorithm recreates the finite-element grid, depending on the position of the wetting front, to maintain the highest node density in those domain regions where |2h/z2| is maximal. Consequently, for this solver, the total number of nodes changes with time during a simulation. The solver incorporates a finite-element, Picard time-iterative model using linear elements and a Galerkin formulation at each time step. Simulations were reported for ponded (5-cm ponding depth) infiltration into clay and sandy loam soils initially in hydrostatic equilibrium with groundwater tables at 15 m depth for clay and 5 m depth for sandy loam soil (Fig. 1). The RES-1D solver effectively prevented numerical instabilities for the extreme head gradients associated with advancing wetting fronts in these two soil profiles.
Cao and Kitanidis (1999) presented an example of category (iii) approaches (combinations of mathematical alteration of governing PDEs and incorporation of adaptive grid refinement or LAGR algorithms) for numerical solutions of multidimensional water flow. Modeling single-phase fully saturated flow in highly heterogeneous formations represents a challenge due to the tendency for flow to channel into preferential flow paths. In the case of very large conductivity variation, the flow regime resembles a network. Relatively fine grids are required to maintain low truncation errors and to prevent numerical dispersion associated with simulation of such flows. However, the use of uniform fine element meshes to capture the high-flow zones tends to be computationally inefficient since a coarser mesh may adequately describe typically slow flow rates that exist in most of the domain. Cao and Kitanidis (1999) applied a mesh-adaptive approach to a dual-flow FEM formulation to simulate steady groundwater flow in heterogeneous isotropic porous media. Two uncoupled elliptic equations for steady flow were utilized, one describing hydraulic potential head H
 | [17] |
 | [18] |
was used to adapt the mesh through bisection refinements, unrefinements, edge swaps, and node moving. The error was defined as = - l where l is a piecewise linear solution for Eq. [18] using linear elements. The mesh adaptation procedure produces a triangulation on which granularity varies gradually, improving its geometric quality and reducing the error related to the discretization. The procedure by Cao and Kitanidis (1999) starts by obtaining a solution utilizing a fairly coarse initial mesh. A higher-order approximation is used to estimate the error distribution, which serves as a guide to locate regions of high flow rates and consequently adjusts the grid granularity. A hierarchical-style local refinement scheme bisects the longest edge of the triangle of interest. Areas with higher flow rates are discovered by a series of mesh adaptations. An auxiliary unrefinement tool is used to delete unnecessary previously set points. An optional node-moving process can further improve the grid. Every refinement or unrefinement process is followed by a node-moving adjustment that improves the mesh quality. Also, error estimates are computed on the basis of updated FEM solution. Refinement of the grid is imposed for high-flow areas and coarsened for low-flow areas depending on accuracy levels specified. Thus, adaptive refinement of the grid ensures both quality of the flow solution and significant reduction of overall processing time.
Cao and Kitanidis (1999) utilized a numerical experiment to demonstrate that successive mesh adaptation efficiently enhanced accuracy for two-dimensional groundwater flow in heterogeneous aquifers. The experiment involved simulation of groundwater flow in a 1 m by 1 m square domain where the logarithm of hydraulic conductivity is assumed to be randomly distributed. On the adapted meshes, the grid is highly nonuniform in node density, with higher density in higher conductivity areas. Large differences between the first and fifth adaptive meshes are clearly obvious in Fig. 2. The first adaptive mesh was coarse, having regular triangulation with unit length of 0.05 m in the x and y directions (1916 finite elements). Error estimates indicated that the first adaptive mesh was inadequate. The fifth adaptive mesh had 9845 finite elements. Nodes were distributed unevenly, so that density is higher where needed as determined with a posteriori error estimate analysis. The CPU percentage utilized for grid optimization in the simulation was <20%. Streamlines computed for the first and fifth adaptive meshes are also shown in Fig. 2. Streamlines from the initial coarse grid tend to miss the channeling of flow in preferential flow paths. A numerical solution based upon a very fine uniform mesh with 0.01 m in both x and y directions (20 000 finite elements) was observed to be less than sufficient to resolve high conductivity areas. The adaptive approach described the groundwater flow more effectively and efficiently than did the very fine uniform mesh.
|
|
Solute Transport in Porous Media: Special Resolution Requirements |
|---|
|
TOPNOTESABSTRACTINTRODUCTIONPertinent Features of...Adaptive Grid Refinement MethodsWater Flow in Porous...Solute Transport in Porous...Examples of Local Adaptive...CONCLUSIONSREFERENCES |
|---|
The three-dimensional continuity equation for convective-dispersive transport of multiple solutes during water flow in fully saturated subsurface porous media can be expressed as
 | [19] |
is porosity, D is the dispersion tensor, 
i is a sink–source term, 
q/ is the seepage velocity, q is the Darcy flow velocity, and Nc is the number of species (Wolfsberg and Freyberg, 1994). For transient water flow in variably saturated subsurface porous media, the continuity equation for multi-species transport becomes
 | [20] |
is volumetric water content and 
q/ is the pore velocity. A simplified form of Eq. [2] was utilized in obtaining Eq. [20]. For two-dimensional systems, the dispersion tensor D may be defined by the components
 | [21] |
j and k are the velocity components in two Cartesians directions, and L and T are the longitudinal and transverse dispersivities. The first and third terms on the right-hand side of Eq. [20] represent hydrodynamic dispersion, and the second term represents convection.
The i term in Eq. [20] represents chemical processes such as sorption on solid surfaces, stochastic exchange on solid surfaces, solution-phase reactions, and precipitation–dissolution of solids (Wolfsberg and Freyberg, 1994). For a single reactive solute species that undergoes simple instantaneous nonlinear Freundlich sorption, the corresponding sink term becomes
 | [22] |
b is assumed constant for a given strata of porous media and the exponent 
< 1.0 for many solute species. For steady flow, the corresponding form of Eq. [20] simplifies to
 | [23] |
 | [24] |
When multiple interacting species are present, subsurface Eq. [20] represents a system of coupled equations where coupling occurs through the sink terms i. For competitive instantaneous ion exchange between multiple cations, the i terms are typically given as functions of concentrations C of both the ith species and competing species (Mansell et al., 1993). Typical units for C in the solution phase and S in the exchange phase are moles of charge per cubic meter and per megagram, respectively, where associated units of b are megagrams per cubic meter. Exchange isotherms relating S and C with fixed exchange coefficients are assumed for binary combinations of species. The cation exchange capacity defined as the summation
 | [25] |
Equation [20] provides mixed parabolic–hyperbolic behavior for solute transport during transient water flow in variably saturated soil. For conditions of dispersion-dominant transport parabolic behavior prevails, but during convection-dominant transport hyperbolic behavior prevails. Sharp fronts of concentration may appear during convection-dominant transport, creating a need for increased resolution localized in numerical models. For the case of convection only, transport occurs along characteristic lines that follow water flow; whereas for dispersion only, transport occurs along and between characteristic lines (Yeh, 2000). Three categories of numerical methods are utilized to solve for solute transport: Eulerian (computations at nodes of a fixed grid), Lagrangian (computations at nodes moving with the fluid [particles]), and Eulerian-Lagrangian methods. When convection dominates transport, Eulerian methods tend to generate numerical oscillation (wiggles) and spreading. Numerical schemes that resolve numerical oscillation tend to yield excessive numerical spreading, whereas schemes that resolve spreading tend to create severe oscillation (Yeh, 2000). Oscillation is often the beginning of instability, whereas excessive spreading gives inaccurate solutions and may cause nonconvergence for nonlinear problems. The most severe limitation of Eulerian methods is that none of them can be safely applied for large Courant numbers (Cr = t/z).
An example of category (i) numerical approaches for solute transport is the flux-corrected transport (FCT) algorithm, which provides an Eulerian method for overcoming both numerical oscillation and dispersion (Boris and Book, 1973). The FCT algorithm, which utilizes a weighting average of the flux computed by a low- and a high-order scheme, has been successfully used to simulate unstable solute transport (Chang and Slattery, 1990). The FCT employs a clipping mechanism by forcing calculated concentration values to stay between the given maximum and minimum values for each node.
Wolfsberg and Freyberg (1994) presented an example of category (ii) numerical approaches for solute transport. A LAGR algorithm was utilized in conjunction with an Eulerian solver that utilized an explicit second-order central finite difference spatial discretization and a second-order Runge–Kutta time stepping scheme to demonstrate efficient simulation of solute transport in groundwater. The LAGR tracks error-prone regions and supplies high-resolution subgrids locally as needed and maintains relatively few nodes elsewhere on a coarse base grid. The algorithm proceeds through five steps during each integration on the base grid: (i) truncation error estimation, (ii) node flagging and subgrid generation, (iii) interpolation of initial and boundary conditions on subgrids, (iv) integration through multiple steps on subgrids, and (v) updating coarse grid solution. Important features include a unique method for detecting a priori where the numerical error is unacceptable, variable time step control, which allows smaller time steps on subgrids than on the base grid, and a modular framework that allows easy exchange of PDE solvers to accommodate different problem formulations. The process of computing local truncation error at every node is inexpensive relative to the simulation cost because only coarse grid solutions are required for estimates. Subgrids with higher resolution than the original base grid are created and placed in regions of unacceptable numerical error. Speed gained through enablement of parallel computations outweighed the time step limitation required of the explicit numerical method used by Wolfsberg and Freyberg (1994). Any FDM using a uniform discretization and of the same order in space and time can be coupled with the LAGR routines used by Wolfsberg and Freyberg (1994).
An example of multiple species transport of three cations in a uniform flow field demonstrated that the LAGR by Wolfsberg and Freyberg (1994) achieved numerical solutions with accuracies comparable to those achieved with a uniform fine grid but with only 28% of the computational cost. The controlling PDEs for transport of the three cations are coupled through reaction terms. Due to iterations among multiple PDEs, solutions to this reactive multispecies problem cost significantly more CPU time per node than nonreactive single-species cases. Transport of three cation species that undergo instantaneous competitive exchange was simulated during steady flow into a two-dimensional unidirectional flow field with homogeneous medium (x = 0.18 m d-1, y = 0, and dispersivities L = T = 0.45 m). Initially, the solution and exchange phases were assumed saturated with species 1 with C1 = CT = 0.5 mg L-1 and S1 = ST = 0.5 mg L-1, where CT is the total solution normality and ST is the concentration of exchange sites. Species 2 and 3 were applied as a strip source centered on the upstream boundary (source length is 0.6 of the 400-m domain width) with solution-phase concentrations C2 = C3 = 1.0 mg L-1. Equilibrium partitioning coefficients (Kij = [Sj/Si]/[Cj/Ci]), K12 = 2 and K13 = 4, were specified such that invader species 3 exchanges more readily than invader species 2, which in turn exchanges more readily than the resident species 1. As the two invader species move into the flow domain, they preferentially displace the exchange phase of species 1. A fixed coarse spatial grid was imposed with 4-m spacing and a 5-d time step. A fixed fine grid was imposed with 1-m spacing and 1.25-d time step. Simulated spatial distributions of species concentrations for the fine grid, coarse grid, and LAGR were compared for an elapsed time of 1500 d. The numerical error in the coarse grid solution is greatest in the region where competitive exchange occurs between all three cation species. The LAGR subgrids and the fixed fine grid provided greater resolution, which reduced numerical error. The leading concentration front for species 1 was not covered with a subgrid because it was beyond the influence of the reactive invader species. Imposing subgrids for only two species rather than three resulted in savings in CPU time. Early in the simulation, the front for species 1 was sharp, but with time dispersion moderated the front sharpness so that the coarse grid resolution was adequate at 1500 d.
Trompert (1993) presented another example of category (ii) numerical approaches for solute transport. A LAGR solver was utilized to investigate nonlinear, brine transport in heterogeneous porous media located in groundwater. Disposal of hazardous waste in salt formations results in special problems related to brine transport. In the vicinity of rock salt formations (e.g., salt domes), salt concentrations may become very large, so as to influence groundwater flow (Zegeling et al., 1992). Adaptive grid methods are especially valuable for such problems where locally steep concentration fronts commonly occur. Trompert (1993) demonstrated that adaptive grid methods can compute a solution to such problems with locally the same resolution as on a very fine uniform grid, but with less computational cost. Equations for transient groundwater flow and salt transport were assumed to be coupled through fluid density and dynamic viscosity µ, which were assumed to be related to salt concentration C through the functions
 |
 | [26] |
(log[2.0] used here) is an experimental coefficient, and o (1000 kg m-3) and µo (10-3 kg m-1 s-1) are reference density and dynamic viscosity, respectively, for solute-free water. Thus, the coupled flow and transport equations were
 |
 | [27] |
The parameter J is the dispersive salt flux, q is Darcy water flux, p is water pressure, D is the 2 by 2 dispersion tensor, is porosity, g (9.81 m s-2) is acceleration of gravity, and k is permeability of the porous media. For time integration of these equations, an implicit Euler method was used for the first time step and a second-order implicit Gear-type BDF stiff ODE solver with variable coefficients for following time steps where variable step sizes are taken (Trompert, 1993). Standard second-order central finite differences were used for space discretization. Linear interpolation was used for obtaining initial and boundary conditions.
Local-uniform-grid refinement (LUGR) methods generally start from a coarse-base grid covering the entire domain. Finer and finer uniform subgrids are created locally in a nested manner in regions of large spatial variations. The LUGR method used by Trompert (1993) performs integration on a series of nested, local-uniform finer and finer subgrids. Subgrids are created up to a level of refinement where sufficient spatial accuracy is reached, and their location and shape are adjusted after each time step. The space domain is considered to be a rectangle, and all grids in use are uniform and Cartesian. Interfaces demarking inhomogeneities are assumed to coincide with cell edges in the numerical approximation. Heuristic error monitors are used to control mesh refinement and the variable-time step sizes.
For each time step, the computation starts at the coarse-base grid by using the most accurate solution available since coarser-grid-solution values are always updated by the finer grid solution. Eight operations are performed during each time step:
Trompert (1993) demonstrated his LUGR method with brine displacement of fresh water (i.e., salt-free water) for a vertical two-dimensional reservoir. Four distinct regions (I, II, III, and IV) in the reservoir were identified with different intrinsic permeabilities and dispersivities. Permeabilities for Regions I, II, III, and IV were 10-13, 10-15, 10-10, and 10-13 m, respectively. Lateral and transversal dispersivities for these four regions were 0.008, 0.005, 0.010, and 0.008 m, respectively, and 0.0016, 0.0010, 0.0020, and 0.0016 m, respectively. The top of the reservoir was open. Saturated brine was injected through the opening on the left-hand side at the bottom. During injection, brine moved slowly into Region III toward the top of the reservoir, initially bypassing Region IV with 1000-fold smaller permeability than Region III such that a very sharp transition in salt concentration developed at the interface between these two regions. With time, dispersion tended to smooth the moving front. Later, the front passed the interface between Regions III and I and then moved into Region I with 1000-fold smaller permeability than Region III, bypassing Region II with 100-fold smaller permeability than Region I. Much later, salt penetrated Regions II and IV from III at approximately the same time. At first, salt penetrated Region IV at the top-left corner and later at the entire interface. Steady state occurred when eventually the saturated brine has spread out over the entire domain. Salt concentrations in the reservoir were simulated using LUGR at t = 60 000 s when two grid levels were used in the model. The coarsest of the two grids was 20 x 20 m. User specified time and space tolerances were 0.10 and 0.25, respectively. The two-grid LUGR computations were 1.8 times faster than those with a 40 x 40 m fine uniform grid.
The Lagrangian approach or moving-grid method is generally less favored in operational transport models than the Eulerian approach because of complexity of application. Computation of partial derivatives becomes very complicated in distorted Lagrangian-point networks (Yeh, 2000).
Zegeling et al. (1992) applied a Lagrangian method [category (ii) numerical approach] to simulate nonlinear, brine transport for groundwater flows in the vicinity of salt domes. Coupled equations for groundwater flow and salt transport are those of Eq. [27]. An implicit FDM approach based on the method of lines (MOL) method for time-dependent PFDs was utilized. With evolving time the spatial grid is automatically refined in regions with large spatial transitions. One grid-smoothing procedure was employed to generate a spatially smooth grid and another for avoiding temporal grid oscillations. Two types of automatic grid-adaptation were used. For one type, integration occurs on grids that spatially equidistribute a relevant measure of the error. The equidistribution is realized in a dynamic Lagrangian approach where the grid is adapted continuously in time. This feature makes it possible to accurately and efficiently follow steep traveling fronts of concentration. A second type of adaptation serves to cope with rapid temporal transitions and involves the use of variable step sizes in the numerical integration. An advantage of this Lagrangian approach is that the grid movement often softens the solution behavior in time, so that larger time steps can be taken than on a fixed spatial grid. The moving-grid finite difference method was demonstrated for two one-dimensional numerical problems with sharp traveling concentration fronts. The first derivative C/x was used as a monitor to place a number of points just within the front where large values occurred.
Simulation of groundwater flow in heterogeneous aquifers with high contrast in equation coefficients can be a difficult task (Cao and Kitanidis, 1999). An accurate description of the velocity field has been shown to be a prerequisite for a realistic analysis of convective-dispersive transport of reactive solutes (Cao and Kitanidis, 1997). In mixing-limited reactions, use of inaccurate velocity causes miscalculation of the rate of mixing and thus the rate of chemical reactions (Cirpka et al., 1999). Neglecting spatial variability may be inappropriate when heterogeneity has a strong influence on mixing. Cao and Kitanidis (1999) reported a procedure for high-accuracy computation of flow in heterogeneous aquifers that utilized a dual-flow formulation and adaptive gridding [an example of category (iii) numerical approaches].
Lagrangian-Eulerian methods (LEMs) combine the best aspects of Eulerian and Lagrangian approaches. Lagrangian-Eulerian methods provide further examples of category (iii) numerical methods. All LEMs solve the Lagrangian form of the transport equation at the nodes of a fixed grid (Yeh, 2000). Particle tracking methods form a variant of LEMs in which particles are introduced into the domain. Each particle is associated with a spatial coordinate and a discrete quantity of mass. The number of particles and location of the introduced particles depend upon the initial concentration field, boundary conditions, and artificial sources or sinks. During each time step, the particles are moved forward with convection and dispersion, the number and location of particles is converted back to concentrations for the fixed grid nodes, concentration changes resulting from biochemical reactions are computed, and mass associated with each particle is recomputed. Particle tracking methods can accurately handle sharp gradients and small sources of mass. However, a very large number of particles and a fine support grid are required to conserve mass in particle tracking approaches.
Backward methods of characteristics (i.e., hybrid LEMs) form another variant of LEMs. Unfortunately, peak clipping is not eliminated by the hybrid LEM approach. Yeh (2000) utilized a Lagrangian-Eulerian decoupling with an adaptive zooming hidden fine mesh approach to solve transport equations to minimize peak clipping. The algorithm included "adaptive local grid refinement based on exact peak capture and oscillation free behavior" (i.e., ALGR-EPCOF). The ALGR-EPCOF algorithm was designed to capture peaks and valleys and high curvature areas to within a specified error tolerance to eliminate spurious oscillation, numerical diffusion, and peak clipping. This algorithm can be used with both FDM and FEM numerical approaches. It can also be used for multidimensional cases.
Spatial distributions of solute concentration obtained wi
