PMFCT-2D: A transport simulator for various grid peclet numbers

PMFCT-2D: A transport simulator for various grid peclet numbers

Aimo, N J


In most natural aquifers, the heterogeneity of the permeability is primarily responsible for the scale dependency of the dispersivity. Gueven et al. (1986) among others noted that one reason why several calibrated models, based on the advection-dispersion equation, yielded large apparent longitudinal dispersivities was that the seepage velocity was assumed to be uniform, at least throughout vertical sections of the flow field. Based on experimental and modeling results, researchers have recently come to realize that the better the flow field is represented, the smaller the dispersivities become that are needed to match observed and predicted data (Domenico and Schwartz, 1990). Moltyaner and Killey (1988) found that the values of dispersivities, determined from time-concentration profiles, were virtually equal to those obtained from laboratory columns. They used detailed knowledge of the flow field at the site of a tracer experiment.

As more detail on the heterogeneity of the subsurface is implemented in numerical models, transport algorithms have to be able to handle high Peclet numbers (Pe). The Peclet number is defined as Pe = v Delta x/D, where v is the seepage velocity, Delta x is the grid spacing, and D is the hydrodynamic dispersion. Most currently available codes, however, exhibit significant numerical dispersion when the Peclet number is larger than 2. When such codes are used to model advection-dominated transport, a large number of nodes may be necessary to obtain accurate results. In practice, unrealistically high dispersivities are often used to decrease the demand on the simulator being used (Gelhar et al., 1985).

In addition, Dougherty and Bagtzoglou (1993) have indicated that the use of classical simulation algorithms could create numerical dispersion. In general, numerical errors tend to spread (disperse) the solute more than physical processes alone. Dougherty and Bagtzoglou (1993) point out that for many contaminants, regulatory and design decisions are based on small dimensionless concentrations; thus, decisions based on classical numerical model solutions at low concentrations must be taken cautiously.

Several papers have been published in the last few years describing transport algorithms for high Peclet numbers (e.g., Illangasekare and Doll, 1989; Hills et al., 1994; Putti et al., 1990; Sudicky, 1989). However, most of the published algorithms have not been implemented in widely available software. In this paper, we present a brief overview of PMFCT-2D, a computer code that is based on the Flux Corrected Transport (FCT) technique. FCT was introduced by Boris and Book (1973; 1976) and Book et al. (1975), and was later enhanced by Zalesak (1979). Even though the FCT technique has been in use within the fluid dynamics community, it has not yet been used in regular geohydrologic practice (Dougherty and Bagtzoglou, 1993). The PMFCT-2D code, an extension of the algorithm presented by Hills et al. (1994), is able to simulate one- and two-dimensional transport in heterogeneous variably saturated porous media for any Peclet number, under transient or steady-state flow conditions. Two-dimensional plan-view simulations with nonuniform aquifer thickness are also possible with PMFCT-2D. The general FCT technique is discussed in the next section, followed by an overview of the PMFCT-2D code. The algorithm is demonstrated by two illustrative problems.

The FCT Technique

Originally, the FCT technique was developed as a new approach to numerically solve the continuity equation (advective transport) in the presence of steep density gradients and hydrodynamic shocks. In these cases, standard algorithms usually have failed to provide accurate solutions (Boris and Book, 1973). Classical finite-difference schemes are inherently inefficient for the simulation of problems with steep gradients (Book et al., 1975). However, FCT algorithms can be applied to solute transport problems where steep gradients need to be advected, or to advection-dominated problems with superior results. The Eulerian FCT finite-difference algorithms developed by Boris and Book (1973; 1976) are able to represent steep gradients and shock fronts quite accurately, without the smearing caused by numerical dispersion, which is usually present in the classical algorithms. For problems where dispersion dominates advection (low Peclet numbers), the results are essentially the same whether or not an FCT algorithm is used (Book et al., 1975).

It should be noted that FCT is a technique, not a finite-difference nor a finite-element scheme, or a specific algorithm. It is possible to derive many different numerical algorithms, involving different finite-difference, finite-element, or other schemes, all incorporating the FCT technique.

An FCT algorithm consists of three main stages: an advective or transport stage (Stage I), a diffusion stage (Stage II), and an antidiffusive or correcting stage (Stage III). Here the term diffusion is used, following usage by Boris and Book (1973), for numerical smearing of the component being transported. Stages I and II can be applied separately or as a single operation, depending on the numerical scheme. All stages are mass conservative. The particular choices of schemes used for each stage allow FCT algorithms to accurately simulate shocks and strong gradients. The desired result of applying Stages I and II is characterized by significant diffusion of the transported field, and the lack of overshoot or undershoot ripples (typical of high-order schemes). This result often can be achieved by using low-order schemes such as the upwind/donor-cell scheme. Stage III corrects the excess diffusion by calculating and limiting the so-called antidiffusion fluxes. The correction is done by first finding a high-order solution for the advective fluxes and then calculating the difference between these and the low-order fluxes of Stage II. The result is the antidiffusive fluxes, which are then applied to the low-order solution with the aid of the flux-limiter.

The flux-limiter is a critical element of the FCT technique, and different flux-limiters are possible. The implementation of the FCT technique in PMFCT-2D uses the flux-limiter introduced by Zalesak (1979), which is multidimensional and superior to earlier flux-limiters.

In essence, the technique involves applying sufficient numerical diffusion everywhere to avoid artificial ripples (Stages I and II) and then canceling it with an equal amount of antidiffusion where it is not needed (Stage III). The criteria by which Stage III is carried out is built into the flux-limiter. The flux-limiting process is a nonlinear operation.

The PMFCT-2D Code

The PMFCT-2D code uses the algorithm presented by Hills et al. (1994), with modifications that allow the inclusion of nodal sources and variable cell thickness in the third dimension. The first modification allows solute to enter the simulated domain at points other than the boundaries, while the latter enables the simulation of two-dimensional areal transport in aquifers with variable thickness (depth-averaged concentration).

PMFCT-2D solves the advection-dispersion equation in a variably saturated porous medium with changing thickness in the third dimension:

(equation omitted)–(1)

where C is the time (t) dependent and space (x) dependent solute concentration, R is the retardation factor, theta is the water content, q is the Darcy velocity, b is the thickness perpendicular to the two-dimensional plane being simulated, D is the dispersion tensor, and Q is a sink/source term.

When the time derivative is replaced by a finite-difference approximation, an expression is obtained that relates the solutions of two consecutive time-steps, n and n + 1, separated by a period of time, Delta t:

(equation omitted)–(2)

Purely advective transport of the solute can be represented by introducing the quantity (R theta bC) sup adv , such that

(equation omitted)–(3)

Substitution of equation (3) into equation (2) yields an expression for the dispersive transport

(equation omitted)–(4)

(equation omitted)–(5)

where S sup adv is a source term that incorporates the effects of the purely advective transport.

As pointed out by Hills et al. (1994), an important advantage of this FCT algorithm is that different numerical schemes may be used to solve the advective and dispersive phases. In PMFCT-2D, the advective transport is solved using the FCT technique with an explicit predictor-corrector scheme. The dispersive transport is solved using a fully implicit, conventional, finite-difference scheme. These particular choices result in a very fast, stable, and efficient code.

Stages I and II of the FCT technique are achieved in PMFCT-2D using the explicit upwind finite-difference scheme:

(equation omitted)–(6)

(equation omitted)–(7)

where Delta f is the difference in f between the right and left faces of the cell, and Delta g is the difference in g between the top and bottom faces of the cell. This scheme is first-order accurate in space.

The predictor step begins with C sup in containing the solute concentrations of the previous time-step C sup n , which are used to compute Delta f and Delta g, and then C sup out . For the corrector step, C sup in and C sup out are averaged, and the procedure is repeated, with the final values in C sup out representing the low-order solution C sup L .

Stage III involves the use of a fourth-order approximation to the cell-averaged mass-fluxes, (character omitted) and (character omitted). The high-order fluxes, denoted as F and G, are:

(equation omitted)–(8)

(equation omitted)–(9)

(equation omitted)–(10)

where the subscripts i and j are nodal indices, and the + 1/2’s imply the quantity is defined at the cell face instead of the cell center (node).

A first approximation of F and G (predictor) is done with C sup in = C sup n , and the resulting values of F and G are used to find C sup out :

(equation omitted)-(11)

where Delta F and Delta G are found in the same way as Delta f and Delta g above. Averaging C sup in and C sup out allows new values for F and G to be calculated (corrector), which then are used to find the antidiffusive fluxes A sup x and A sup y at each cell face

A sup x = F – f and A sup y = G – g–(12)

The antidiffusive fluxes need to be limited, to avoid introducing the ripples characteristic of the high-order approximation. In PMFCT-2D, the flux-limiter of Zalesak (1979) is used; the reader is referred to his work for the details.

The final results from the advective transport C sup adv are found from

(equation omitted)–(13)

where C sup L is the low-order solution (Stages I and II).

Dispersive transport is performed on the concentration field from the previous time-step C sup n using C sup adv , resulting from the advective transport, as a cell source. A conventional, fully implicit, finite-difference scheme, which is unconditionally stable, is used to solve for the final concentration field.

The only boundary condition implemented in PMFCT-2D is third type (also known as Cauchy, Robin, or Danckwerts). This boundary condition is consistent with the volume-averaged calculation of the solute concentrations (Parker and van Genuchten, 1984). The third-type boundary condition is written as

(equation omitted)–(14)

and is enforced during the advective transport. In other words, mass enters the simulated domain only during the advective stage. For the dispersive transport stage,

(equation omitted)–(15)

is enforced.

One restriction needs to be observed. To maintain positivity in the resulting concentration fields, the maximum Courant number, Cu, should be less than 0.5 (Boris and Book, 1973). The Courant number is defined as Cu = v delta t/delta x. The Courant number relates the size of the time-step Delta t to the grid spacing. In practice, keeping Cu

The FCT algorithm implemented in PMFCT-2D may seem to require excessive calculations. However, the explicit formulation for the advective stage is very efficient and fast because no matrix inversions need to be performed. In addition, for calculations involving steep gradients, a given error tolerance can be achieved by an FCT algorithm two to eight times faster than standard implicit methods (Boris and Book, 1973). PMFCT-2D has been written with efficiency in mind. One-dimensional arrays are used exclusively to avoid overhead array pointer calculations, and nodal calculations are coded to avoid memory page faults.

A generalized flow chart of the algorithm used in PMFCT-2D is shown in Figure 1. (All figures omitted) Before starting a PMFCT-2D simulation, the user needs to solve the flow problem and save the appropriate velocity, saturation, and cell thickness fields for input to PMFCT-2D. The input file describing the simulation domain and transport problem is read first, followed by a file containing the velocity, saturation, and cell thickness information. The low-order solution is computed next, followed by the high-order solution. Subsequently, the antidiffusive fluxes are computed using the low- and high-order solutions, and then applied to the low-order solution using the flux-limiter. At this point, the advective solution for the time-step is complete. This solution then is used in combination with the final solution of the previous time-step to calculate source terms for inclusion in the dispersive stage. The dispersive stage is calculated next, completing the final solution for the time-step. At this point, concentrations for the next time-step may be calculated, using the same velocity, saturation, and thickness fields, in the case of steady-state flow, or new fields may be read, in the case of transient flow.

Sample Problems

In this section, two illustrative examples are discussed. In the first example, the accuracy of the PMFCT-2D simulator is verified by comparison to an analytical solution for the one-dimensional propagation of a rectangular solute wave. In addition, the performance of the PMFCT-2D code for this problem is compared to results obtained with the STOMP (White and Oostrom, 1996) and CFEST (version SC-01, Feb. 1988, Gupta et al., 1987) codes. The two latter simulators use the classical integrated finite-difference and Galerkin finite-element techniques, respectively. In the second example, a two-dimensional contaminant plume emanating from a leaking tank is simulated by STOMP and PMFCT-2D. A comparison of the results obtained by both simulators is shown.

One-Dimensional Transport Comparison

The analytical solution for the one-dimensional, transient, rectangular wave, transport problem was modified from a solution given by van Genuchten and Alves (1982), to account for translation of the initial conditions along the x-axis. Assuming the solute to be conservative, and given the initial and boundary conditions,

(equation omitted)–(16)

where x sub 1 and x sub 2 define the edges of the rectangular wave, the solution of the one-dimensional advection-dispersion equation is

(equation omitted)–(17)

where v is the pore-water velocity (LT sup -1 ), and D is the dispersion coefficient (L sup 2 T sup -1 ).

The numerical parameter values used in the simulations of this problem are the same as those used by Dougherty and Bagtzoglou (1993), as recommended by the Convection-Diffusion Forum during the VII International Conference on Computational Methods in Water Resources (Baptista et al., 1988). This choice was made with the aim of preserving a common comparison. The computational domain ranges from 0

Results for the analytical and the three numerical solutions (PMFCT-2D, STOMP, and CFEST) are shown in Figure 2 for Pe = 2, and in Figure 3 for Pe = 50. The comparison on Figure 2 shows a close match between the analytical and the PMFCT-2D solutions. The STOMP results already show significant smearing as a result of numerical dispersion. When advection is the dominant transport mechanism (Figure 3), the PMFCT-2D solution still produces a reasonable match with the analytical solution. As expected, the two other simulators perform poorly when the Peclet number is large. For both cases, Pe = 2 and Pe = 50, the PMFCT-2D simulator performs better than the STOMP and CFEST codes.

Two-Dimensional Transport Comparison

The second sample problem deals with a hypothetical contaminant plume developing from a leaking storage tank in a steady-state flow field. The conceptual model is depicted in Figure 4. In this example, the composition of the subsurface is assumed to be well-known, allowing the use of small dispersivity values in transport simulations. The subsurface consists mainly of a medium sand. Two heterogeneities are distinguished: a 20-m long by 1-m thick clay layer with a low permeability and a 15-m long by 1-m thick gravel layer with a high permeability. All porous materials are considered to be homogeneous and isotropic. Flow in the unconfined aquifer is from left to right, driven by a hydraulic gradient of 0.005. For the vadose zone, the following retention relation was used to describe the saturation as a function of the air-water capillary pressure head (van Genuchten, 1980):

(equation omitted)–(18)

where S sub e , S, and S sub m are the effective, actual, and residual saturations; alpha, n, and m are saturation curve shape parameters; and h sub aw is the air-water capillary pressure head. The shape parameters m and n are related as follows: m = 1 – 1/n. The values of alpha, n, and S sub m used in the simulations are listed in Table 1. (Table 1 omitted) The relative permeabilities, k sub r , were approximated by (Mualem, 1976):

(equation omitted)–(19)

The tank sits on level ground and has been leaking for an extended period of time with a constant leakage rate of 20 l/day, resulting in a steady-state flow field. After the steady-state flow field was established, a certain contaminant started to leak from the tank. The steady-state velocities and saturations were computed with the STOMP code using a uniform cell size of 1 m by 0.5 m (Figures 4 and 5). From Figure 5, it is clear that the leaking water is diverted by both the clay and gravel layers. The water diversion on top of the clay layer is a direct result of the low permeability of the clay, while the diversion on top of the gravel layer can be mainly attributed to the relatively high alpha value of the gravel. The water pressures within and adjacent to the gravel layer are not large enough to allow significant amounts of water to flow into the gravel.

The computed velocities at the faces of the cells, and the nodal saturations, are subsequently used as input data for transport simulations with both the PMFCT-2D and STOMP codes. For the simulation with the PMFCT-2D code, the maximum Peclet number was 40, corresponding to a longitudinal dispersivity of 0.025 m. For the STOMP simulation, a maximum Peclet number of 2 was allowed, resulting in a longitudinal dispersivity of 0.5 m. For both simulations, the transverse dispersivity was 10% of the longitudinal dispersivity. Solute plumes at various times, computed with the PMFCT-2D and STOMP codes, are shown in Figures 6 and 7, respectively. As expected, the STOMP results show more mixing (dilution) of the contaminant than the PMFCT-2D results in both the vadose and saturated parts of the computational domain. The STOMP simulation indicates that, after two years, the saturated region downstream of the tank is almost completely contaminated (Figure 7a). At the same point in time, the PMFCT-2D simulation predicts that only a small part of the saturated zone is polluted (Figure 6a). At five years, the PMFCT-2D code predicts the contaminant is transported in relatively high concentrations in the upper part of the saturated region (Figure 6b), while STOMP predicts mixing over the entire thickness of the saturated region (Figure 7b). This comparison demonstrates that the choice of the transport simulator can have an important impact on the resulting predictions and subsequent regulatory decisions.

Conclusions As more detail about subsurface heterogeneity is incorporated in models, the dispersivities needed to match observed and predicted data become smaller. With smaller dispersivities, transport codes are required to handle high Peclet numbers. For many codes based on classical numerical techniques, it is not practical to use small dispersivities because of Peclet number limitations. In addition, studies have shown that the use of classical simulation algorithms can create numerical errors that are largest where concentrations are lowest because of excessive numerical dispersion.

Although several papers have been published in the last few years on algorithms to simulate high Peclet number transport, most of these algorithms have not yet been incorporated into readily available codes. A general-purpose transport code, PMFCT-2D, has been introduced, that is able to simulate advection-dominated transport characterized by high Peclet numbers. PMFCT-2D splits the advective-dispersive transport into two stages: advection and dispersion. The advective transport is solved using the FCT technique with an explicit predictor-corrector finite-difference scheme, while the dispersive transport is solved using a conventional, fully implicit, finite-difference scheme. The code is designed to accurately model solute transport in variably saturated, heterogeneous porous media for steady-state and transient flow conditions. Steep concentration gradients can be simulated without the numerical dispersion found in classical numerical schemes.

The performance of the PMFCT-2D code was compared against the CFEST and STOMP codes for a one-dimensional transport problem. The results show that PMFCT-2D performs significantly better than CFEST and STOMP, which are based on classical numerical schemes. In another example, a contaminant plume in a well-defined, variably saturated, heterogeneous subsurface was simulated using the PMFCT-2D code with a high Pe and using the STOMP code with Pe = 2. The results show remarkable differences in the predicted contaminant plumes, especially in the saturated region, which could lead to different regulatory and clean-up decisions.


Development and testing of PMFCT-2D was done under the Integrated Environmental Monitoring Initiative, a Pacific Northwest National Laboratory (PNNL) directed research and development initiative. PNNL is managed and operated for the U.S. Department of Energy by Battelle Memorial Institute under Contract DE-AC06-76RL0 1830.


Baptista, A., P. Gresho, and E. Adams. 1988. Reference Problems for the Convection-Diffusion Forum. VII International Conference on Computational Methods in Water Resources, Cambridge, MA.

Book, D. L., J. P. Boris, and K. Hain. 1975. Flux corrected transport. II: Generalizations of the method. J. of Computational Physics. v. 18, pp. 248-283.

Boris, J. P. and D. L. Book. 1973. Flux corrected transport. I: SHASTA, A fluid transport algorithm that works. J. of Computational Physics. v. 11, pp. 38-69.

Boris, J. P. and D. L. Book. 1976. Flux corrected transport. III: Minimal-error FCT algorithms. J. of Computational Physics. v. 20, pp. 397-431.

Domenico, P. A. and F. W. Schwartz. 1990. Physical and Chemical Hydrogeology. John Wiley and Sons, New York.

Dougherty, D. E. and C. Bagtzoglou. 1993. A caution on the regulatory use of numerical solute transport models. Ground Water. v. 31, pp. 1007-1010.

Gelhar, L. W., A. Mantoglou, C. Welty, and K. R. Rehfeldt. 1985. A review of field-scale physical solute transport processes in saturated and unsaturated porous media. Electric Power Research Report EA-4190, Project 2485-5, Palo Alto, CA.

Gupta, S. K., C. R. Cole, C. T. Kincaid, and A. M. Monti. 1987. Coupled Fluid, Energy and Solute Transport Model: Formulation and Users Manual. Battelle Memorial Institute, Columbus, OH. Report BMI/ONWI-660.

Gueven, O., F. J. Molz, and J. G. Melville. 1986. Comment on an advection-diffusion concept for solute transport in heterogeneous unconsolidated geological deposits by Gillham et al. Water Resour. Res. v. 22, pp. 89-91.

Hills, R. G., K. A. Fisher, M. R. Kirkland, and P. J. Wierenga. 1994. Application of flux-corrected transport to the Las Cruces Trench site. Water Resour. Res. v. 30, pp. 2377-2386.

Illangasekare, T. H. and P. Doll. 1989. A discrete kernel method of characteristics model of solute transport in water table aquifers. Water Resour. Res. v. 25, pp. 857-867.

Moltyaner, G. L. and R.W.D. Killey. 1988. Twin Lake tracer tests: Longitudinal dispersion. Water Resour. Res. v. 24, pp. 1613-1627.

Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. v. 12, pp. 413-522.

Parker, J. C. and M. Th. van Genuchten. 1984. Flux-averaged concentrations in continuum approaches to solute transport. Water Resour. Res. v. 20, pp. 866-872.

Putti, M., W.W.-G. Yeh, and W. A. Mulder. 1990. A triangular finite volume approach with high-resolution upwind terms for the solution of groundwater transport equations. Water Resour. Res. v. 26, pp. 2865-2880.

Sudicky, E. A. 1989. The Laplace transform Galerkin technique:

time-continuous finite element theory and application to mass transport to groundwater. Water Resour. Res. v. 25, pp. 1833-1846.

van Genuchten, M. Th. 1980. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Science Soc. of America J. v. 44, pp. 892-898.

van Genuchten, M. Th. and W. J. Alves. 1982. Analytical solutions of the one-dimensional convective-dispersive solute transport equation. USDA/ARS Technical Bulletin 1661.

White, M. D. and M. Oostrom. 1996. Subsurface Transport Over Multiple Phases (STOMP). User’s Guide. PNL Report (in press).

Zalesak, S. T. 1979. Fully multidimensional flux-corrected transport alporithms for fluids. J. of Computational Physics. v. 31, pp. 335-362.

N. J. Aimo and M. Oostrom, Environmental Technology Division, K9-36, Pacific Northwest National Laboratory, P.O. Box 999, Richard, Washington 99352 (. J. Aimo is the corresponding author).

Copyright Ground Water Publishing Company Jan/Feb 1997

Provided by ProQuest Information and Learning Company. All rights Reserved