A stable-and-convergent numerical groundwater model is developed herein. Stability being a necessary condition for model operation, the emphasis is on convergence, i.e., whether the numerical model can reproduce the differential equation with sufficient accuracy.
A significant feature of the model is its simplicity, which makes it attractive and
useful for a wide variety of applications, including water quality. The numerical model is
an explicit formulation of the two-dimensional (
A specific focus of the model development is the assessment of stability and convergence. Explicit models are subject to a strict stability criterion, which must be satisfied
if the model is to simulate natural conditions in a realistic way (O'Brien
After more than 20 years of implicit model development of groundwater flow, the
choice of an explicit model at this juncture requires further explanation. It is well known
that implicit models are unconditionally stable. Lesser known is the fact that implicit
models are subject to a convergence criterion which effectively places an upper limit on
the time step. Examples of lack of convergence in surface-water and sediment-transport
modeling have been documented by Ponce
The three-dimensional diffusion equation for flow through porous media is (Rushton and Redshaw 1979; Kresic 1997):
in which K, and _{y}K are the hydraulic conductivities along the _{z}x , y , and z coordinate axes, respectively [ L T ^{-1}]; W = volumetric flux per unit volume, representing sources (+) or sinks (-) [ T ^{-1}]; S = specific storage of the porous material [ L_{s}^{-1}]; and t = time.
In two dimensions in plan, Eq. 1 simplifies to:
Assuming a homogeneous isotropic aquifer of hydraulic conductivity
Defining the hydraulic diffusivity of the aquifer (Freeze and Cherry 1979) as
leads to
Equation 5 represents two-dimensional flow in a homogeneous isotropic aquifer of properties
The spatial interval is set as
Eq. 6 reduces to
In which
The time-averaged source term in velocity units is
in which ^{-1}] (due to rainfall and/or irrigation) at spatial node (j, k), and b = aquifer depth [L].
The time-averaged sink term in discharge units is
in which ^{3} T^{-1}] at spatial node (j,k).
Since transmissivity
and diffusivity
in which
With Eqs. 7 and 13, Eq. 14 reduces to
For the special case
Equation 15 calculates piezometric head at the unknown node (
The boundary conditions can be of Dirichlet or Neumann type (Kinzelbach 1986).
Dirichlet conditions specify the head Neumann type A boundaries may be specified by linear extrapolation from the two immediately neighboring interior nodes. Linear extrapolation results in a finite gradient, which amounts to a permeable boundary. Neumann type B boundaries may be specified by relocation of the neighboring interior node. Relocation results in a zero gradient, which amounts to an impermeable boundary.
Mixed Neumann type A/type B conditions, appropriate for semipermeable boundaries or temporally varying inflows, can be specified as follows:
in which
The model was subjected to the following four tests to examine its numerical properties:
This condition tests the model's ability to return to the steady-equilibrium baseline condition following the specification of (a) a depleted water table as initial condition, (b) permeable boundaries, and (c) no pumping.
This condition tests the model's ability to achieve a steady-equilibrium cone of depression following the specification of (a) the baseline steady-equilibrium initial condition, (b) permeable boundaries, and (c) pumping.
This condition tests the model's ability to return to a steady-equilibrium nonbaseline condition following the specification of (a) a depleted water table as initial condition, (b) impermeable boundaries, and (c) no pumping.
This condition tests the model's ability to simulate the linear depletion of the aquifer
following the specification of (a) the baseline steady-equilibrium initial condition,
(b) impermeable boundaries, and (c) pumping.
First, the groundwater model was tested for convergence (O'Brien et al. 1949). For this purpose, a permeable hot-start test was used, wherein the head recovery to steady-equilibrium baseline condition is a proof of the model's convergence. The test prototype specification was a 100-km Stability requires that
A series of sixteen (16) runs were performed to test the convergence of the groundwater
model under varying cell Reynolds number h = (_{0}h - 100) in a square area of 5 km × 5 km located in the middle of the computational field (51 × 51 = 2601 grid nodes).
_{ref}
Four cell Reynolds numbers ( h = 0 m. From Eq. 18, the time interval for _{bottom}D = 1.0 is Δt = 6.944 hr.
Therefore, the time intervals to be used in the test series are: Δt = 6.944, 3.472, 1.736, and 0.868 hr, corresponding to D = 1.0, 0.5, 0.25, and 0.125, respectively.
Figures 2 (a) to (d) show head recovery at centerfield node (50, 50) for reference heads D = 1.0, 0.5, 0.25, and 0.125. For values other than D = 1, a head recovery deficit is present, signaling the model's inability to return to steady-equilibrium baseline condition, i.e., lack of convergence. Table 1 shows head recovery deficit h as a function of reference head _{def}h and cell Reynolds number _{ref}D.
The following conclusions are drawn from Fig. 2 and Table 1: The head recovery deficit *h*= 0, i.e., the head recovery is complete, only for_{ref}*D*= 1.Head recovery deficit increases nonlinearly with decreases in cell Reynolds number from 1.0 to 0.125. Head recovery deficit increases linearly with increases in reference head from 100 to 800 m. The condition *D*= 1 is not only the most accurate but also the most economical. For smaller values of*D*, the model is less convergent and the run takes more processor and total time to complete.
The permeable/impermeable hot/cold tests were carried out by specifying the same
geometry and aquifer properties as for the convergence test runs. The cell Reynolds
number was set at
For the h = (_{0}h - 100) in a square area of 5 km × 5 km located in the middle of the computational field. The aquifer bottom is set at _{ref}h = 0 m. Then, the steady-equilibrium aquifer volume is 50,000 hm_{bottom}^{3}.
For the ^{3}. A battery of 17 pumps is symmetrically positioned across the field, at a distance of 1,414.2 m apart, making the shape of an X (Fig. 3), and each pumping p = 250 L s^{-1}.
The initial aquifer volume is 47,399 hm
The initial aquifer volume is 50,000 hm
The initial aquifer volume is 47,399 hm
The initial aquifer volume is 50,000 hm
A stable-and-convergent two-dimensional groundwater model of a homogeneous isotropic aquifer has been developed and tested under a wide range of flow conditions. The model is explicit and can account for sources/sinks and head/flux or permeable/impermeable boundaries. A significant feature of the model is that it sets the cell Reynolds number Testing of the model under hypothetical conditions showed that the model is stable, convergent, and mass conservative. The permeable hot-start test converged to the steady-equilibrium baseline condition after 15.85 yr. The permeable cold-start test converged to a steady-equilibrium cone of depression after 12.43 yr. The impermeable hot-start test converged to a steady-equilibrium non-baseline condition after 8.34 yr. The impermeable cold-start test properly simulated the linear depletion of the aquifer.
Douglas, Jr., G. 1956.
On the relation between stability and convergence in the numerical solution of linear parabolic and hyperbolic equations.
Freeze, and Cherry. 1979.
Kilzenbach, W. 1986.
Kresic, N. 1979.
O'Brien, G. G. Hyman, and M. A. Kaplan. 1950. A study of the numerical solution of partial differential equations.
Ponce, V. M., H. Indlekofer, and D. B. Simons. 1978. Convergence of four-point implicit water wave models.
Ponce, V. M., H. Indlekofer, and D. B. Simons. 1979. The convergence of implicit bed transient models.
Roache, P. 1972.
Rushton, K. R., and S. C. Redshaw. 1979. |

190130 14:20 |