A 3-D finite element conjugate gradient model of subsurface flow with automatic mesh generation

A 3-D finite element conjugate gradient model of subsurface flow with automatic mesh generation

Giuseppe Gambolati and Giorgio Pini Istituto di Matematica Applicata, Universita' degli Studi, Padova, Italy

Tullio Tucciarelli lstituto di ldraulica, Universita' di Palermo, Italy

The 3-D flow modelling of groundwater systems of realistic size generally requires a big effort for the preparation of the input data as well as large computational costs. A numerical finite element model (MAITHREE) is developed for the efficient analysis of the steady and unsteady behaviour of natural confined 3-D basins. Starting from an initia! triangular grid the code automatically generates a set of tetrahedral elements in each of the geologic units or subunits specified by the user in the vertical profile. The original element incidences list is then rearranged in order to provide conforming 3-D elements throughout the domain. The model is designed with a view to saving much of the labour involved in setting a 3-D grid and to providing flexibility as well as economical convenience through a high computational efficiency. The latter task is achieved by the aid of a solver based on the modified conjugate gradient (MCG) method which has proved to be an excellent technique for the solution of large linear finite element sets of sparse 3-D subsurface equations. Some examples derived from both hypothetical and real-world situations are discussed to illustrate the innovative features of MAITHREE and its computational performance.

INTRODUCTION Realistic simulations of natural 3-D hydrologic basins by numerical methods require a large number of elements and nodes (and hence much effort to prepare the input data) as well as expensive computational resources. A model which effectively deals with real practical problems has to provide an adequate response to the following points: (a) usdr facility and convenience in view of the labour involved in setting the 3-D mesh; (b) high numerical efficiency to reduce the processing time which easily becomes a critical factor for moderate to large size problems. One of the first 3-D analyses was developed by Freeze ~ with the aid of finite differences. Freeze t used a line successive over-relaxation scheme (LSOR) to solve his model in large hydrologic systems. France 2 resorted to a finite element approach with cubic isoparametric elements and to a direct solver. The strong implicit procedure (SIP) was adopted by Trescott and Larson 3 and by Winter 4 in their 3-D finite difference approximations. A 3-D model with mixed order isoparametric elements was developed by Gupta and Tanji 5 to simulate steady state problems and was recently extended to transient analyses6. In their code Gupta et al. 6 store the f.e. matrix in compressed mode with the Accepted May 1985.Discussioncloses May 1986.

34

Adv. Water Resources, 1986, Volume 9, March

elimination of all zero coefficients and use a triangular decomposition into an upper and a lower factor followed by a back substitution step. An excellent review of the major difficulties encountered with 3-D numerical analyses of groundwater flow was contributed by Frind and VergC. They pointed out that the problems are three-fold: (1) the 3-D cost is substantially higher than the 2-D cost of a comparable model; (2) the labour involved in laying out a 3-D grid may be quite formidable; (3) the computational burden of a 3-D model may be increased by the presence of the unsaturated zone or by a moving boundary. Frind and Verge T made use of mixed order isoparametric brick elements. Their solvers were modified versions of the Gaussian elimination and Cholesky band algorithms. The present paper is intended to address points (1) and (2) above. The generation of the 3-D mesh is obtained with an automatic procedure starting from a basic triangular grid defined by the user on the upper surface of the 3-D domain. The planar view of this grid is repeated in each unit or subunit into which the vertical profile of the basin has been preliminarily subdivided. To preserve a high degree of flexibility and maintain the structure of the f.e. matrix as sparse as possible tetrahedral elements are employed. The model is completely general and is particularly suited for regional systems where the 0309-1708/86/010034-08 $2.00 ~' 1986ComputationalMechanicsPublications

A 3-D finite element conjugate gradient model: G. Gambolati, G. Pini and T. Tucciarelli horizontal layer dimension is much larger than the vertical one. The numerical efficiency of the model is attained by the use of the modified conjugate gradients (MCG) which have already proven extraordinarily effective for 2-D flow analysess,9.1O, The automatic generation of the 3-D tetrahedral mesh is described first. Special consideration is given to the reordering of the element incidence list so as to obtain conforming volume elements throughout the domain. Then the 3-D transient model is briefly developed and solved with the aid of MCG. Finally the performance of the code (MAITHREE) is shown with hypothetical and real world applications of 3-D groundwater flows.

1

TetrahedronNode,s,

//// I

/

/

II

I )

4563

2 )

1243

31

2453

I

/

l

L,:'" I,' I ~ 6

5 Fig. 2. Each prismatic element is subdivided into 3 tetrahedrons. This subdivision is not unique

The equation of transient groundwater flow in 3-D basins may be written as (assuming that the coordinate directions are parallel to the principal directions of anisotropy):

~, i

1" /

4

F.E. EQUATIONS

Oh\

/

II lI

A U T O M A T I C 3-D MESH GENERATION AND

_iK Oh'~ O f

3

1

3

3 2

4

~h\

(1) where h is the hydraulic head, K=, Ky, and K= are the principal components of the hydraulic conductivity tensor, 7 is the specific weight of water, ~ and fl are the soil and water compressibility, respectively, n is the porosity and Q is point-source representing wells or distributed source representing recharge on top of the aquifer system. The term ~,(~+ nil) is often denoted by S, and referred to as specific storage coefficient. Equation (1) has to be supplemented by appropriate initial and boundary conditions. To solve equation (1) in space we adopt a mesh of tetrahedral elements which are generated automatically from a set of prismatic elements with triangular crosssection. The planar view of these sections is similar in each unit or sub-unit into which the vertical profile of the basin has been subdivided (Fig. 1). Each prism is divided into 3

,tI ~ %~ I

Jt

#1%

#

a)

#

•

•

b}

Fig. 1. Schematic portion of a 3-D groundwater system made up of prismatic elements with triangular crosssections (a); vertical set of prismatic elements with the same triangular section (b)

¢

5

6

,. ,-" "

7

7

TetrahedronNodes

e

TetrahedronNodes

!1

5763

4}

6782

2)

1523

5)

6842

3) 5 7 2 3 6) 6 4 3 2 Fi9. 3. Exploded view of two adjacent prisms. I f the decomposition of the prisms into tetrahedrons does not satisfy certain rules, non-conforming elements may be generated, as shown here: on the common .face 2 3 6 7 the continuity requirement for h is violated

tetrahedrons (Fig. 2). It is interesting to notice that the decomposition is not unique. This might lead to a violation of the continuity requirements on the common lateral side of two adjacent prisms and hence to the generation of non-conforming elements (Fig. 3). The input data are given by an initial triangular grid defined by nodal coordinates as well as by the nodal numbering and element incidences. The model automatically builds the prismatic elements in each layer (note that these need not have a uniform thickness, see Fig. 4) from which tetrahedral elements are further derived. If N, T and M denote the number of nodes and triangles of the initial 2-D grid and the number of vertical units, respectively, the 3-D mesh will be characterized by L = N ( M + 1) nodes and 3 M T tetrahedrons. Consider now a typical tetrahedron with the vertices i,j, m, n numbered in such a way that the determinant De:

l1

De =

1 1

I1

Xi

Yi

--i

xj

yj

zj

Xm

Ym

Zrn

X.

y.

z.

Adv. Water Resources, 1986, Volume 9, March

(2)

35

A 3-D finite element conjugate gradient model: G. Gambolati, G. Pini and T. Tucciarelli 2

1

1

1

ee=SseVe

1

2

1

1

20

1

1

2

1

1

1

1

2

(6)

The vectorsf andf+a, account for the forcing functions Q and the boundary conditions of Neumann type.

N O N - C O N F O R M I N G TETRAHEDRAL ELEMENTS AND PATCH-TEST

Fig. 4. The vertical layers of the model need not have a uniform thickness in the horizontal dimension is positive. Then we have V~=De~6, lie being the element volume. Taking linear basis functions inside each tetrahedraon e, the hydraulic head he may be expressed through the nodal variables hi, hi, h,, and h. as:

h e = ~ e [(ai+bix +ci5 +diz)hi + (ai + bjx + Qy + djz)h~

h=x+y+z

+ (a,. + b.,x + c,.y + d,.z)h,. +(a. +b,x +c.y+d,z)h,]

(3)

where:

ai=

Ci=

i[

yj

zj

Xm

Ym

2ra

Yn

zn

1 1

Xj Xm

Zj 2m

1

Xn

2n

h i = --

di ~ --

1

y,.

2m

1

y.

Zn

1

xj

yj

1

x,n

Ym

1

x.

y.

h,-f-f+a,

(4)

e 1 • H,~= 3~)] (Kxb, b~+Kyc,G+Kzd, d~) r,s=i,j, m, n

(5) Adv. Water Resources, 1986, Volume 9, March

x~ + y~ + zo

t8)

The analytical equation for node 6 has been derived as well but it will not be given here for the sake of brevity. The exact value of h6 in the test is 5. Table 1 provides the f.e. solution /~ for different ratio bl/b 2 (see Fig. 5). Looking at Table 1, it is seen that the patch test is not satisfied unless b l = b 2 and that the difference /~6-h, increases as the ratio b~/b2 increases. Hence nonconforming elements are not recommended when the layer thicknesses in the vertical profile are not uniform. A simple mechanism to generate conforming tetrahedraons can be derived with the aid of Fig. 2. Denote by i,j, m the nodes of a typical initial triangle. I f N is the overall number of nodes of the 2-D mesh which is taken to coincide with the upper boundary of the aquifer basin, we number the nodes underlying node i as:

i+N, i+2N, i + 3 N , etc . . . .

which provides the potential head behaviour in the 3-D basin. Matrices H and P above are the stiffness and the capacity matrices and are obtained by properly assembling the local matrices H e and Pe:

36

(7)

satisfies the steady form of the flow equation (1) for uniform values of the permeability coefficient and will be used for the patch-test. On the nodes i = 1, 2..... 12 i:~6 we prescribe h as provided by equation (7) and write the f.e. equation for the unknown variable h6 on node 6. For the elements to pass the test, the numerical solution k76 has to coincide with the theoretical hydraulic potential, i.e.:

~o =

and the remaining coefficients are determined by cyclic permutation of the indices i,j, m, n and of the determinant signs. Following the variational approach of Ritz-Galerkin t1'12 and integrating in time with the CrankNicolson t3 scheme lead to the algebraic finite element system:

= ~-H

AS was outlined in the previous section a prism can be divided into 3 tetrahedrons in more than one way. Therefore, if an appropriate subdivision strategy is not defined a priori, non-conforming elements are easily generated and the continuity of the potential head may fail to occur somewhere within the flow domain. However it is known that non-conforming grids may converge if the elements satisfy the so-called 'patchtest,12.14,~ 5. Let us consider a patch of 6 prisms as shown in Fig. 5. Following the procedure described in the previous section, each prism is partitioned into 3 nonconforming tetrahedrons (see Fig. 3). The solution:

The same scheme applies for the nodes which underly nodes j and m. Now we rearrange the triple i, j, m in ascending (descending) order. This implies that some values of the determinant (2) will turn out to be negative, so we take the absolute value of (2). It can be easily realized that, if the tetrahedrons are built according to the mechanism evidenced in Fig. 2, a partition similar to that shown in Fig. 3 cannot occur and all conforming elements will be obtained. This is illustrated with a simple example in Fig. 6.

A 3-D finite element conjugate gradient model: G. Gambolati, G. Pini and T. Tucciarelli

T H E M C G SOLVER

1i

4

Setting:

-5-2P D = H +-~-t

I

bl

I

l

8

-]-

(9)

we know that D is a sparse symmetric positive definite matrix. If m - 1 denotes the maximum number of nodes adjacent to a given node (in the sense of element connectivity) in the initial 2-D mesh, the largest number o f non-zero coefficients in a row of D is 2m+ I. On the average m = 7 and therefore D has usually 15 non-zero elements per row. The bandwidth of D is equal to N, i.e. it is rather large for realistic f.e. systems. Since the overall number L of nodes is very easily greater than 103, we conclude that D is a very sparse matrix with large bandwidth. The model has been solved with the M C G scheme. Its numerical efficiency and accuracy have already been demonstrated for 2-D f.e. subsurface problems by Table 1. Finite element solution of the patch test tbr d!(lbrent ratios between the thickness of the vertical units

t

•

i • -

I I I I I I I I I I I I

b2

bl/b 2

h~,

h~-t~,

1

5.0000

0.0000

4 10 20

5.0320 5.0675 5.0844

0.0320 0.0675 0.0844

5 4

II I

t 9

K- , , ~ ,t

V,I \'

G

12

10

i'

.-"FI.-"'.'i".. I

II

7.~'-----I----,~,.8

L.""-.l.-'""-d g

I

12

5

i 1

"

Fig. 5. Patch o f prismatic elements sectioned into nonconforming tetrahedrons for the application of the patchtest

i

i~.q

1

1

ll

L,/.- ,,,q 1

I

+

]

,

Thus the procedure to be incorporated with the finite element model can be shortly summarized as follows: (1) the input data of the initial 2-D mesh are set with the most convenient (arbitrary) ordering of the nodes (2) the nodes of each individual triangle are rearranged in ascending (descending) order (3) the nodal numeration of the various vertical projections from the initial mesh is accomplished as specified in the present section (4} from each prismatic element conforming tetrahedrons are automatically generated following the scheme depicted by Fig. 2.

4~13

+ 7

10

9

8 9

g

12

Fig. 6. I f the nodes of each triangle in the initial 2-D grid are given in ascending order and the prisms are sectioned into tetrahedrons according to the mechanism drawn in Fig. 2 conforming elements are generated throughout the 3-D flow domain

Adv. Water Resources, 1986, Volume 9, March

37

A 3-D finite element conjugate gradient model: G. Gambolati, G. Pini and T. Tucciarelli Gambolati a and Gambolati and Perdon 1°. For very sparse matrices of large size the M C G technique is superior to the other techniques commonly employed 9. Moreover, contrary to the other techniques, M C G is rather insensitive to the structure of D, in particular to its bandwidth, and this make M C G especially suited for 3-D models of the type developed in the present paper. Writing equations (4) under the standard form:

dh=b

(10)

the M C G recursive equations read1°: hi+ l = h i + c t i P i Pi+l =K-lri+ ri + l =

O~i =

piirri p-fDpi

b-

1+fliPi

D h i + I : ri - o~iDpl

fli =

(11)

rr÷IK- tDp i plOp i

Po = K - tr o = K - l(b - D h o ) In the above equations rl is the residual at the ith iteration and K- ~= (LL r)- i is the preconditioning matrix of the M C G scheme, L being the lower incomplete factor of D 1°A6. The initial guessed solution h o is taken equal to: ho=K-lb i.e. equal to the result from the first iteration of the residual correction scheme 8. The rate of convergence of relationships (11) is controlled by the spectral behaviour of the M C G iterative matrix: E = D K -~

(12)

It may be shown 9 that the closer to 1 the ratio 2~/2 L between the maximal and the minimal eigenvahte of E, the faster the M C G convergence.

Fig. 8. Detail of the upper portion of the 3-D prismatic mesh around the pumping well Table 2. I f the stiffness matrix H is stored in a double precision vector no ill-conditioning occurs and the material balance is satisfied rery well in stead), simulations

Stiffness matrix H

Relative error on the material balance

single precision double precision

0.703 0.213.10-7

ILLUSTRATIVE RESULTS The numerical performance of the model MAITHREE is first analysed with a theoretical example representing a circular 6 layer aquifer system as shown in Fig. 7. The

Figl 7. Schematic diagram of the 3-D 6 layer aquifer system. The initial 2-D triangular grid is also drawn

38

Adv. Water Resources, I986, Volume 9, March

outer radius and the overall thickness are assumed to be 1000 m and 10 m, respectively (regional formation). The aquifer has a different vertical as well as horizontal hydraulic conductivity in each subunit and is pumped at a constant rate from a partially penetrating well with radius 0.1 m. The initial triangular grid is made up of 416 triangles and 221 nodes and is schematically drawn in Fig. 7. Note that the size of the triangles is not uniform and varies quite significantly as one moves from the inner to the outer boundary. The 3-D mesh automatically generated from the original 2-D grid has 7488 tetrahedrons and 1547 nodes. As was described earlier prismatic elements are built (see in Fig. 8 a small representative detail of the prismatic mesh around the abstraction well) and from each one of them 3 conforming tetrahedrons are derived. The model has been tested in both steady and unsteady simulations on a medium size IBM computer. When the time integration step At is very large some ill-conditioning problems may arise. This is accounted for by the behaviour of the stiffness matrix H which dominates over the capacity matrix and tends to become ill-conditioned if its coefficients, although computed in double precision accumulators, are stored in a single precision vector. Generally matrix H is sensitive to the geometry of the mesh and ill-conditioning may occur in 2-D steady analyses as well. Switching to double precision leads to quite satisfactory results as may be seen from Table 2

A 3-D finite element conjugate gradient model: G. Gambolati, G. Pini and T. Tucciarelli

rm

j 0

5

10 NUMBER

15 20 OF ITERATIONS

25

Fig. 9. Convergence rate of the MCG solver in steady (At = ~ , solid line) and unsteady (At = 10 -3 day, dotted line) simulations

which shows the errors on the overall material balance for the two situations. In transient problems the illconditioning may disappear depending on the size of At. However it is recommended that H is computed and stored in any case using double precision. Figure 9 shows the behaviour of the M C G solver in steady and unsteady simulations. The mean residual

(provided in Table 3 as well) and becomes closer to 1 for smaller values of the time integration step. Thus the number of multiple eigenvalues of E grows and the desired solution is obtained in fewer iterations 1° A second application of the model M A I T H R E E is concerned with the determination of the hydraulic conductivity of a water table aquifer located at Piana dei Colli, north of Palermo, Sicily, Italy (Fig. 11) by trial and error. A schematic vertical cross-section of the formation may be seen in Fig. 11. The aquifer rests on an impervious bedrock at 15 m depth and is made up of an upper layer A which is much more permeable than a lower layer B. The subunits A and B are 12 m and 3 m thick, respectively (Fig. 11). Initially the water table is 3 m beneath the ground surface. The aquifer is pumped at a constant rate from a fully penetrating well with radius equal to 0.1 m. Three pumping tests have been performed with different withdrawal rates until a steady state profile was practically obtained. The observed drawdowns have been matched using the results from the steady version of MAITHREE. Layers A and B have been vertically divided into 7 and 1 prismatic blocks, respectively. The planar view of the mesh is the same as that shown in Figs 7 and 8 while a portion of the vertical cross-section close to the pumping well is given in Fig. 12. The final steady state was achieved through an iterative procedure where the elements of Table 3. The eigenvalues of the MCG iteration matrix E are clustered around 1 and .fall into a much smaller interval than those 01 the .finite element matrix D. I f At decreases the length of the spectral interval of E and D also decreases Matrix E

Matrix D

At

rM=X/

(days)

is plotted against the number of iterations. Note that the M C G convergence is extraordinarily fast and an accurate solution is achieved after less than x / L iterations. This also holds true for At---* oc where the f.e. matrix (H in this case) is not diagonally dominant as is evidenced by Fig. 10 providing the histogram of the ratios:

zo 0.1 0.01 0.001

2t 1.276 1.276 1.276 1.276

2L

2t/)4.

0.0162 0.0248 0.0337 0.0498

78.8 51.5 37.9 25.6

2~

2L

2j/2L

387.5 0.542"10-5 7.1 ' 107 387.5 0.808'10-s 4.8"107 387.5 1.05"10-5 3.7'10 7 387.5 1.33"10-5 2.9'107

100 % 80

Hi

hi i

between the sum of the absolute values of the extradiagonal terms and the diagonal element of H. Although the diagonal dominance was proven to be a desirable condition vis-a-vis the performance of the M C G scheme 1°, the previous example emphasizes the robustness of the M C G solver which turns out to converge very fast in 3-D subsurface flow models irrespective of the diagonal dominance of the f.e. matrix. Generally the convergence is improved in transient problems (Fig. 9) as is also shown by the spectral behaviour of the iteration matrix E (equation (12)) whose extreme eigenvalues 21 (E) and 2L(E) are given in Table 3. The minimal eigenvalue has been assessed by an efficient reverse power-conjugate scheme recently developed by Gambolati and Perdon 1~. It is worth noting that the ratio 2x(E)/2L(E) is much smaller than the corresponding ratio between the maximal and minimal eigenvalues of D

60

40

20

0 0

!

0.4

i

0.8

1.2

I

1.6

,

11,

2.0

Fig. 10. In 3-0 finite element models of regional aquifer systems the stiffness matrix is not diagonally dominant

Adv. Water Resources, 1986, Volume 9, March

39

A 3-D finite element conjugate gradient model: G. Gambolati, G. Pini and T. Tucciarelli

TYRRHENIAN SEA

~

~

~

~ ? ~

SEA

model to reproduce the observed water decline in the well and the withdrawal rates. The values of K4 have been determined by trial and error. The outer radius R e of the aquifer system is assumed to be equal to 410m. F r o m Table 4 it appears that the response for KA is not unique but we succeeded in restricting significantly the range of variation of K4 between 2 0 . 1 0 - 6 m / s and 40-10-6m/s. The results from a sensitivity analysis on the influence o f R e shows (Table 5) t.hat R e does not play a significant role on the determination of KA and hence the differences between t he maximal and minimal estimates of K4 are to be accounted for by the errors affecting the field measurements and by the well losses. CONCLUSION A new 3-D finite model of subsurface flow has been developed and implemented in a computer program called M A I T H R E E . The most innovative elements of the model consist of:

El

(a) automatic generation of a 3-D tetrahedral mesh in a

¥ INITIALWATERTABLE FINALWATERTABLE

! 1

T E

QUATERNARY CALCARENITE

"1 |

• I ! :In,

B

QUATERNARY GREY CLAY NUMIDIAN FLYSCH ( OLIGOCE.NE) Fig. 11. Schematic cross-section and geologic profile of the water table aquifer at Piana dei Colli, Palermo, Sicily, Italy

±

,,*/'/~','7//~,'///f4 ~/'/)~'~'//,'/'/'//, ;"/S

Fig. 12. Cross-section of the 3-D finite element model in the steady configuration for the water table aquifer at Piana dei Colli, Palermo, Sicily (test 1) Table 4. Well drawdowns and pumping rates as measured in the field and as reproduced by the 3-D model. The lowering o f the seepage face and the isotropic permeability o f layer A are also given

layer A were allowed to deform to fit the moving position of the free surface (Fig. 12). The mechanism by which the water table and the height of the seepage face are adjusted at each iteration is similar to that described by Isaacs ~8. Possible variants of this approach may be found in the works by N e u m a n and Witherspoon ~9, Verruijt 2° and France et al. 21 Under the assumption that the permeability Ka of layer B is significantly less than that of layer A, the model shows that the potential profiles and the discharge rates are not very sensitive to KB. The latter has therefore been set equal to 1 0 - 6 m / s as suggested by geological considerations. Table 4 gives the drawdowns in the well and the corresponding pumping rates as observed in the field. Table 4 also provides the computed lowering of the intersection between the water table and the well wall and the hydraulic conductivities KA of layer A which allow the

40

Adv. Water Resources, 1986, Volume 9, March

Well drawdown (m) Pumping rate (l/s) Lowering of the seepage face (m) Permeability of layer A as predicted by the model (m/s)

test 1

test 2

test 3

7.51 0.67

2.88 0.45

0.93 0.26

3.91

2.06

0.67

20.10 -6

25.10 -6

4 0 . 1 0 -6

Table 5. Sensitivity of the calibrated hydraulic conductivity KA to the outer radius R e

KA (m/s)

Re [m)

167

410

test 1 test 2 test 3

18' 10-6 22- 10-6 37.10 -6

25.10 -6

20.10 40.10

954 -6 -6

22.10 - 6 27.10 -6 43" 10-6

A 3-D finite element conjugate gradient model: G. Gambolati, G. Pint and T. Tucciarelli multilayer aquifer system starting from an initial 2-D t r i a n g u l a r grid on one surface of the b o u n d a r y of the domain: (b) solution to the 3-D f.e. e q u a t i o n s by the modified c o n j u g a t e gradient ( M C G ) technique. F e a t u r e (a) implies a s u b s t a n t i a l saving of the l a b o u r involved in p r e p a r i n g the i n p u t data while feature (b) ensures a very high c o m p u t a t i o n a l performance. A suitable strategy of n o d a l r e n u m e r a t i o n provides c o n f o r m i n g v o l u m e elements t h r o u g h o u t the flow domain. Applications Of M A I T H R E E have shown that the M C G solver achieves the desired solution after a n u m b e r •of iterations which is less than the square root of the p r o b l e m size, b o t h in steady a n d u n s t e a d y simulations, irrespective of the d i a g o n a l d o m i n a n c e of the f.e. matrix. T h e model has also been successfully used to interpret u n c o n f i n e d field data. T o this purpose a n a p p r o p r i a t e iterative procedure to adjust the position of the m o v i n g free b o u n d a r y has been i n c o r p o r a t e d with the code. It m a y therefore be c o n c l u d e d that M A I T H R E E represents a very efficient a n d practical tool for the analysis of regional 3-D g r o u n d w a t e r flows whenever 2-D schemes are unrealistically simple.

6 7 8 9

10

11 12 13 14 15 16

REFERENCES 1 2 3 4 5

Freeze, R. A. Three-dimensional, transient, saturatedunsaturated flow in a groundwater basin, Water Resour. Res., 1971, 7, 347-366 France, P. W. Finite element analysis of three-dimensional groundwater flow problems, J. Hydrology, 1974, 21,381-398 Trescon, P. C. and Larson, S. P. Solution of the threedimensional groundwater flow equations using the strongly implicit procedure, Journal tgf Hydrolooy, 1977. 35, 49-60 Winter, T. C. Numerical simulation of steady state threedimensional groundwater flow near lakes, Water Resour. Res., 1978, 14, 245-254 Gupta, S. K. and Tanji. K. K. A three-dimensionalGalerkin finite

17 18 19 20 21

element solution of flow through multiaquifers in Sutter Basin, California. Water Resour. Res., 1976, 12, 155-162 Gupta, S. K., Cole, C. R, and Pinder, G. F. A finite-elementthreedimensional groundwater (FE3DGW) model for a multiaquifer system, Water Resour. Res., 1984. 20, 553 563 Frind, E. O. and Verge, M. J. Three-dimensiona] modelling of groundwater flow systems, Water Resour. Res,, 1978, 14. 844-856 Gambolati, G. Fast solution to finite element flow equations by Newton iteration and modified conjugate gradient method, Int. J, Num. Methods Engng., 1980, 15, 6614575 Gambolati. G. and Volpi, G. Analysisof the performance of the modified conjugate gradient method for solution of sparse linear sets of finite element equations, Third Int. Conf. on Finite Elements in Flow Problems, Banff (Canada), 10 pp., !980 Gambolati, G. and Perdon, A. M. The conjugate gradients in flow and land subsidencemodelling,Fundamentals of Transport Phenomena in Porous Media, J. Bear and Y. Corapcioglu Eds, Martinus Nijoff B.V., The Hague, 953-984, [984 Pinder,G. F. and Gray, W. G. Finite element simulation in surface and subsurface hydrology, Academic Press, New York, 1977 Zienkiewicz,O. C. The finite element method, McGraw-Hill, New York, 1977 Crank, J. and Nicolson, P. A practical method for numerical evaluation of solutions of partial differentialequations of the heat conduction type, Proc. Cambridge Philos. Soc., 40, 50-67, 1947 Mitchell A. R. and Wait, R. The finite element method in partial differential equations, John Wiley & Sons, New York, 1977 Dhatt, G. and Thouzot. G. The finite element method displayed, John Wiley & Sons, New York, 1984 Kershaw,D. S. The incomplete Cholesky conjugate gradient method for the iterative solution of systemsof linear equations, J. Comp. Phys., 1978, 26, 43~5 Gambolati, G. and Perdon, A. M. Minimal eigenvalue of large sparse matrices by an efficientreverse power-conjugate gradient scheme, Comp. Methods Applied Mech. Engng., 1983, 41, 1-10 Isaacs,L. T. Location of seepage free surface in finite element analyses, Civil Engineering Transaction, Australia, CE 22, 1980, No 1, 9-17 Neuman, S. P. and Witherspoon, P. A. Finite element method for analyzing steady seepage with a free surface, Water Resour. Res,, 1970, 6, 889-897 Verruijt,A. Theory of groundwater flow, Macmillan,New York, 1970 France, P. W., Parekh, C. J., Peters, J. C. and Taylor, C. Numerical analysis of free surface seepaage problems, d. lrri9. Drain. Div., ASCE, 97, No IRI, 1971. 165-179

Adv. Water Resources, 1986, Volume 9, March

41