A numerical analysis of quasi-3D models for linear multiaquifer systems
P. Teatini, G. Gambolati
Dept. Mathematical Methods and Models for Scientific
Applications, University of Padova, Padova, Italy
INTRODUCTION
When water movement takes place in a multiaquifer system due to
groundwater withdrawal or changes in hydrogeological conditions,
the flow direction is predominantly horizontal in the aquifers
and vertical in the intervening aquitards.
Excluding a fully three-dimensional analysis, which solves the problem
without taking advantage of the specific features of the flow field, numerical
models, known as ''quasi three-dimensional'' models, can be built
based on the assumption that flow conforms with the pattern recalled above.
The error introduced by this simplified
approach has been numerically investigated by
Neuman and Witherspoon (1969), and found to be less than 5%
if the permeability contrast between
aquifer and aquitard exceeds two orders of magnitudes.
Recently Follou et al. (1992),
using a perturbation technique to calculate the subsidence
due to pumping from layered soils, have shown that the error introduced with
the quasi 3-D scheme is of the same order of the ratio
Kz/Kx
between the vertical aquitard and horizontal aquifer permeability.
''Quasi-3D'' linear models have been developed using two theoretical formulations.
The first formulation, followed by
Herrera and Figueroa (1969), Herrera (1970) and
Herrera and Rodarte (1973), is based on Dhuamel's theorem, and leads to a system
of integro-differential equations in which convolution integrals account for
leakage between sandy and clayey layers and exactly link the
hydraulic potential head of each semi-permeable unit to that of the
underlying and overlying aquifers. The convolution integrals
involve two kernels, called ''memory'' and ''influence function'', which
are expressed in the form of an infinite series
(Herrera and Rodarte, 1973). The resulting set of integro-differential
equations is solved by de Marsily et al. (1978) and
Premchitt (1981) with the aid of finite differences,
and the overrelaxation (SOR) method.
Herrera (1976), Herrera and Yates (1977) and
Hennart et al. (1981) used instead the finite element method
and a direct approach with the triangular matrix factorization (LDU).
Following Herrera's approach, Gambolati and Perdon (1984) and
Gambolati et al. (1986) have combined
in an original way the integro-differential
equations with the modified conjugate gradient (MCG) method
(Gambolati, 1980). In Gambolati's numerical implementation,
the convolution integrals are not analytically calculated with the complex
technique of Herrera and Yates (1977). The ''memory'' and ''influence'' functions
are assessed a priori
with a high accuracy for a number of discrete values of the dimensionless
time factor DTv
and used directly in the execution of the integrals with the aid
of simple and inexpensive interpolation schemes.
As has been demostrated by Gambolati and Volpi (1980), the MCG method is more
efficient than the most traditional solution methods,
including triangular factorization with optimal
preliminary reordering and SOR, especially for large finite element systems.
The second approach to solve the quasi 3-D model relies on the full
numerically integration of the original partial differential equations
of the flow. This method has been
used by Fujinawa (1977) and Chorley and Frind (1978).
In Fujinawa's work, aquifers are discretized with a
two-dimensional grid of triangular finite elements. The overlying and
underlying aquitards are discretized by a single
vertical element, along which the head is assumed to vary linearly. The
resulting set of equations is then solved by
Gauss elimination. By distinction, in the model
of Chorley and Frind the aquitards are
discretized by several one-dimensional linear finite
elements. In the neighborhood of a sandy layer, where vertical hydraulic
gradients are steeper, more complex elements and
a higher order interpolation
scheme are used. The discrete overall solution is obtained by iterating between the
aquifer and aquitard solution locally approximated
using the Cholesky factorization.
Also Neuman et al. (1982) developed a numerical model with
1-D finite elements to simulate flow in aquitards. Their most
significant contribution consists of the implementation of
an original adaptive explicit-implicit solver.
The scheme takes advantage of the relatively small hydraulic
diffusivity of the aquitards as compared to that of the aquifers, so that for
a small time step value Dt, most of the equations associated with the
aquitard nodes can be solved explicity with a very low computational cost. The
remaining equations are solved implicitly with a direct scheme. For small
Dt's, the system of linear equations turns out to be virtually decoupled
into a number of smaller subsystems equal to the number
of aquifers. This adaptive strategy brings the efficiency of
Neuman's numerical approach to a level which is approximately equal to that of
the integro-differential method of Gambolati, even if,
as time progresses and Dt increases, a growing number of nodes,
including aquitard nodes, shift from explicit to implicit condition
and the model performance tends to reduce significantly.
Recently, Rivera et al. (1990, 1991) have used the
fully numerical formulation solved with finite differences.
The flow domain is discretized into 2D grids for the aquifers
and 1-D grids for the aquitards. Nested square meshes of
variable size are connetted through the intervening aquitard to properly
simulate steep gradient regimes.
The resulting system of algebraic equations is solved iteratively using SOR.
A few comparisons between these two different formulations have been made.
Herrera et al. (1980) have compared the numerical model of
Chorley and Frind (1978) with their own semi-analytical method
and have concluded that the latter is much more efficient requiring up to
7-10 times less computer time. Neuman et al. (1982) raised
the concern that the result of this comparison could not be entirely meaningful
pointing out that their adaptive explicit-implicit scheme may significantly
reduce the distance between the performances of the semi-analytical and
the fully numerical model.
Gambolati et al. (1986) finally noted that ''the theory underlying
a model and the solver employed for its solution should be considered separatel''
as ''the conceptual superiority of a theoretical formulation can
be partially or totally offset by an inefficient solver''.
Taking the last remark as a starting point, the work described in this paper
is intended to provide a systematic comparison between the
fully numerical approach and the integro-differential approach of
Gambolati et al. (1986). The main goal is to evaluate
the relative efficiency of the two schemes in relation to the general
hydrogeological setting and parameter contrast, and
indipendently of the solver used to solve the final
system of algebraic equations.
Two numerical codes are implemented, one, entitled ''MAINFLW'', based on the
integro-differential approach of Gambolati et al. (1986) and the other,
''MULTSTR'', based on the fully numerical model of
Chorley and Frind (1978) and Neuman et al. (1982).
The two codes are applied to the same test problems and compared in terms of
both CPU time and storage requirement, trying to detect
the breakeven point of the two approaches. Runs are performed on the
IBM RISC6000/560.
