POLF: Two-dimensional finite-element model for predicting the areal flow of pollutant in confined and unconfined aquifers

PII:

Computers & Geosciences Vol. 24, No. 6, pp. 509±522, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain S0098-3004(98)00074-0 0098-3004/98/$ - see front matter

POLF: TWO-DIMENSIONAL FINITE-ELEMENT MODEL FOR PREDICTING THE AREAL FLOW OF POLLUTANT IN CONFINED AND UNCONFINED AQUIFERS G. GOTTARDI* Dipartimento di Ingegneria Chimica Mineraria e delle Tecnologie Ambientali, UniversitaÁ di Bologna, Viale del Risorgimento 2, 40136 Bologna, Italy (Received 13 February 1998; accepted 14 May 1998) AbstractÐA ®nite-element program for modeling two-dimensional (2D) pollutant transport in con®ned and uncon®ned aquifers is presented. The code accounts for the absorption and biological decay of the pollutant, which is supposed to behave as a tracer. The ¯ow and pollutant transport equations are solved by using either a Galerkin ®nite-element (GFE) or a control-volume ®nite-element (CVFE) formulation. The discretization of the ¯ow domain is obtained by using triangular ®nite elements and linear shape functions, whereas integration in time is performed by a ®nite-di€erence (FD) Picard scheme. For two problems of which the analytical solution is known, the results of the simulations carried out by the CVFE and by GFE formulations are compared, taking as evaluation parameters the global mass-balance errors and the average deviation between the numerical and analytical solutions. The listing of the FORTRAN 77 source code, along with the user manual and two input ®les for solving the above mentioned problems are available by anonymous FTP from IAMG.ORG. # 1998 Elsevier Science Ltd. All rights reserved Code available at http://www.iamg.org/CGEditor/index.htm Key Words: Numerical code, Saturated ¯ow, Pollutant transport, Galerkin formulation, Control-volume ®nite-element model.

INTRODUCTION

The physical phenomena governing pollutant transport in con®ned and uncon®ned aquifers are, as well known, advection, dispersion, absorption and biological decay (Friend, 1975; Kinzelbach, 1986; Bear and Verruijt, 1987; Gottardi and Dall'olio, 1994). In the hypothesis that the pollutant is miscible with water and that the viscosity and density of the water±pollutant mixture are not altered in a relevant way by the pollutant concentration, the ¯ow and transport equations are disjointed and may be solved separately. The numerical solution of transport problems with sharp transition or high advection are surprisingly dicult to solve using traditional ®nite element (FE), ®nite di€erence (FD) or integrated ®nite di€erence (IFD) techniques (Gray and Pinder, 1976; Narasimham and Witherspoon, 1976; Dick, 1983; Yu and Heinrich, 1987; Ferraresi and Gottardi, 1990). When the advective term becomes dominant, the transport equation tends to become hyperbolic and the solutions based on traditional FE and FD methods are always plagued with either spurious oscillations (overshoot) or excessive slope of the concentration front (smearing), due to nu*E-mail: [email protected]. 509

merical dispersion. The two phenomena of overshoot and of smearing are closely related. When a scheme is used to minimize the numerical dispersion, overshoot is encountered and when a scheme is used to control overshoot, this is obtained at the expenses of numerical dispersion. In general, an upward scheme is numerically dispersive but limits overshoot and, on the other hand, a centered scheme is free of dispersion error but induces overshooting. Acceptable numerical solutions by FE and FD methods can be obtained by using a ®ne discretization in space and time, to satisfy the stability requirements expressed by the criterion of the Courant and Peclet numbers (Gray and Pinder, 1976; Kinzelbach, 1986). Other methods used to obtain acceptable solutions of the transport equation are based on preconditioning techniques of ``upwind'' methods (Heinrich and others, 1977; Hughes, 1978; Thomas, 1978; Huyakorn and Nilkuha, 1979; Sun and Yeh, 1983; Chavent and others, 1984; Yeh, 1986) or moving ®nite element methods (Gottardi and Venutelli, 1992, 1994). The presentation of this work follows the steps: Ð Description of the governing equations for ¯ow and transport mathematical models in a 2D domain. Ð GFE and CVFE formulation of the elemental matrices used in the 2D code POLF.

510

G. Gottardi

Ð Presentation of the global equations obtained by assembling the elemental matrices for the GFE and CVFE formulations. Ð Presentation of the results of the numerical solutions obtained by the code POLF for two problems whose the analytical solution is known. Ð Concluding remarks. A brief description of the subroutines of the code POLF is given in Appendix A.

THE MATHEMATICAL MODEL

The equations to be fed into a pollutant transport model, in which the pollutant behaves as a tracer, are the mixture ¯ow equation, the pollutant transport equation and the equations expressing the boundary and initial conditions. The mixture ¯ow equation, arising from the material balance equation of the mixture of water and pollutant, allows the distribution to be obtained in the ¯ow domain and in time of the Darcy velocity of the mixture, whereas the integration of the transport equation, arising from the material balance of the pollutant, allows the distribution of the pollutant concentration to be obtained in space and time. If O indicates the ¯ow domain in the two-dimensional Cartesian space x1 and x2 and G is the boundary of the ¯ow domain, the equations of the mathematical 2D model are: (1) In the hypothesis that the pollutant behaves as a tracer (i.e. in the hypothesis that density and the viscosity of the mixture are independent from pressure and pollutant concentration), the mixture ¯ow equation is: ÿ

@ @c …hvi † ÿ hqw ˆ a @xi @t

on O

…1†

where c, xi, vi, h, qw and a are, respectively, the velocity potential of the mixture (L), the ith Cartesian coordinate (L), the Darcy velocity of the mixture in the xi direction (LTÿ1), the thickness of the aquifer (L), a source/sink term (positive for source) (Tÿ1) and the accumulation term coecient, which is di€erent for con®ned and uncon®ned situations. The Darcy velocity of the mixture is: vi ˆ ÿkij

@c , @xj

…1a†

kij being the permeability tensor (LTÿ1) of the porous medium. The velocity potential of the mixture is: p cˆ ‡z …1b† g where p, g and z are, respectively, the pressure of the mixture (MLÿ1Tÿ2), the speci®c weight of the mixture (MLÿ3) and z is the depth (L) measured in the gravity direction and taken positive upward. The accumulation term coecient and the thickness

of the aquifer for the con®ned and uncon®ned cases are: a ˆ Ss h;

h ˆ zt ÿ zb for a confined aquifer

…1c†

a ˆ Sy ;

h ˆ c ÿ zb for an unconfined aquifer …1d†

where Ss an Sy are, respectively, the speci®c storage coecient (Lÿ1) and the speci®c yield (ÿ) and zt and zb are, respectively, the depth (L) of the top and bottom of the aquifer. (2) The transport equation, obtained by considering absorption and biological decay (Kinzelbach, 1986; Gottardi and Dall'olio, 1994), is: ÿ

@ @ …hJi †ÿ …hvi c† @xi @xi ÿ hlne cfr ÿ hqc ˆ ne fr

@ …hc† @t

on O

…2†

where Ji, c, ne, fr, l and qc are, respectively, the speci®c mass ¯ow-rate of pollutant (MTÿ1Lÿ2), the pollutant concentration (mass of pollutant/volume of mixture) (MLÿ3), the porosity of the medium (ÿ), the retardation factor (ÿ), the biological decay velocity of pollutant (Tÿ1) and a source/sink term for the pollutant (MTÿ1Lÿ3). The speci®c pollutant mass ¯ow rate is given by the Fick's law: Ji ˆ ÿne Dij

@c @xj

…2a†

Dij being the hydrodynamic dispersion tensor (L2Tÿ1). The hydrodynamic dispersion tensor is a function of the porous medium and of the pore velocity (Bear, 1979): ui uj Dij ˆ dm dij ‡ aT udij ‡ …aL ÿ aT † …2b† u where dm, aL, aT and dij are, respectively, the coecient of molecular di€usion of pollutant (L2Tÿ1), the longitudinal and transversal dispersivity coecients (L) and the Kronecker operator. u and ui are, respectively, the component of the pore velocity vector (ui=vi/ne) in the xi direction (LTÿ1) and the pore velocity vector module (LTÿ1). The retardation factor is: fr ˆ 1 ‡ wrm

…1 ÿ ne † ne

…2c†

where rm and w are, respectively, the mass density of the porous medium (MLÿ3) and the absorption coecient (L3Mÿ1). (3) The boundary conditions of pre®xed ¯ow for the mixture and pollutant are: hvi
ni ˆ hq~ w …x 1 , x 2 , t†

on G1

…3a†

hJi
ni ˆ hq~d …x 1 , x 2 , t† on G2

…3b†

POLF: 2D ®nite element model

hvi c
ni ˆ hq~c …x 1 , x 2 , t†

on G3

…3c†

where ni is the ith component of the unitary outward normal to a boundary line in the xi direction. qÄw, qÄ d and qÄc are, respectively, the pre®xed speci®c mixture ¯ow rate (LTÿ1) on boundary G1, the pre®xed speci®c di€usion mass rate (MTÿ1Lÿ2) on boundary G2, and the speci®c pre®xed advection mass rate of pollutant (MTÿ1Lÿ2) on boundary G3. (4) The boundary conditions of pre®xed pollutant concentration and velocity potential are: ~ 1 , x 2 , t† c ˆ c…x

on G4

…3d†

~ 1 , x 2 , t† c ˆ c…x

on G5

…3e†

where the functions cÄ and c~ are, respectively, the pre®xed pollutant concentration (MLÿ3) on boundary G4 and the pre®xed velocity potential (L) on boundary G5. For con®ned aquifers, when the thickness h may be considered as constant, Equation (1) is linear, whereas for uncon®ned systems the thickness h depends on c(h = c ÿ zb) and Equation (1) becomes non-linear. Equations (1) and (2) are only coupled by the mixture velocity. GFE FORMULATION OF THE MODEL

In the following a brief sketch of the procedure used to formulate the model by the classical GFE method, along with the structure of the elemental matrices obtained by this method, are presented. To integrate the partial di€erential Equations (1) and (2) the ¯ow domain O is partitioned by using triangular elements Oe with nodes n, p and q. Now we assume that on Oe (see Fig. 1), the variables c(x1, x2, t), c(x1, x2, t) and h (x1, x2, t) can be approximated by a linear combination of their nodal values (Connor and Brebbia, 1976): ^ 1 , x 2 , t† ˆ Cm …t†N e …x 1 , x 2 †, c…x 1 , x 2 , t† ' c…x m ^ 1 , x 2 , t† ˆ Cm …t†N em …x 1 , x 2 †, c…x 1 , x 2 , t† ' c…x ^ 1 , x 2 , t† ˆ hm …t†N e …x 1 , x 2 †, h…x 1 , x 2 , t† ' h…x m m ˆ n, p, q:

511

Oe of the mesh, i.e. N el …x 1 , x 2 † ˆ

1 …al ‡ bl x 1 ‡ cl x 2 † 2De

l ˆ n, p, q:

e

…5† e

Where D is the area of the eth triangle O and the coecients al, bl and cl are obtained from the corner node coordinates xl1, xl2 (l = n, p, q) by anticyclic rotation of the indices n, p and q using the following formulae: an ˆ x p1 x q2 ÿ x q1 x p2 ; cn ˆ x q1 ÿ x p1 :

bn ˆ x p2 ÿ x q2 ; …6†

The geometrical meaning of the coecients al and bl is shown in Figure 2. For the interpolation functions Nen, Nep and Neq associated with nodes n, p and q of Oe the following conditions hold: m N el …x m 1 , x 2 † ˆ dlm

N el …x 1 , x 2 † ˆ 0,

l, m ˆ n, p, q

l ˆ n, p, q 8…x 1 , x 2 †= 2O e :

…7a† …7b†

^ and L2(hÃ, cÃ) be the residuals obtained Let L1(hÃ, c) from substituting the linear approximations of Equation (4) into Equations (1) and (2), respectively. According to Galerkin's method, the orthogonality condition between residuals and interpolation functions is imposed on Oe, yielding: … ^ c†N ^ e …x 1 , x 2 † do ˆ 0 l ˆ n, p, q, …8a† L1 …h, l Oe

… Oe

^ c†N ^ el …x 1 , x 2 † do ˆ 0 l ˆ n, p, q: L2 …h,

…8b†

Each time, Equations (8a)±(b) are equivalent to a system of six ordinary di€erential equations in the nodal unknowns cl(t) and cl(t) of Oe. By substituting in the residual expressions L1 and L2 in Equations (8a)±(b) the linear approximations,

…4†

Nel (x1, x2) are the usual linear interpolation functions, associated with the nodes of a generic element

Figure 1. Basic triangular element Oe along with subcontrol volumes One, Ope, Oqe associated with nodes n, p and q

Figure 2. Graphical representation of coecients a, b and c used in de®nition of interpolations functions Nne, Npe, Nqe of linear triangular element Oe

512

G. Gottardi

Equation (4), and applying the Gauss±Green lemma, thus eliminating second-order derivatives and introducing the pre®xed ¯ux boundary conditions, the following set of equations is obtained: …  e e ^ ij @ N l @ N m do Cm hk @xi @xj Oe …  e e ‡ aN l N m do Cm ˆ Oe

… ÿ qw …

Oe

^ e do ÿ q~ hN w l

… Ge

^ e ds hN l

…9a†

 @ N el @ N em do Cm @xi @xj Oe …  @ N em e ‡ vi h^ N l do Cm @xi Oe …  @ h^ e ‡ vi N em N l do Cm @xi Oe …  ^ e N e do Cm ‡ lne fr hN m l

GFE matrices of ¯ow equation For l, m = n, p, q the structures of the elemental matrices for the ¯ow equation, Equation (9a), for a con®ned aquifer are: … e e ^ ij @ N l @ N m do hk Aelm ˆ e @xi @xj O ˆ

h …k11 bm bl ‡ k12 bl cm ‡ k21 bm cl ‡ k22 cm cl † 4De …11†

^ ij ne hD

B elm ˆ Ss

@ h^ e e N N do Cm @t m l O …  ^ e @ vi N e do Cm ‡ hN m @xi l Oe …  ^ e N e do C_ m ˆ ‡ ne fr hN m l e

…

ne fr

Oe

Oe

^ c N e do ‡ hq l

… Ge

h^q~ d N el ds:

Oe

^ l Nm do hN lˆm

…12a†

l 6ˆ m

applying the ``lumping approximation'', the capacity matrix becomes … e ^ l do ˆ Ss D …2hl ‡ hp ‡ hq †dlm …12b† B elm ˆ Ss hN 12 Oe



‡

…

8 De > > < Ss …3hl ‡ hp ‡ hq † 30 ˆ e > D > : Ss …2hl ‡ 2hm ‡ hk6ˆl6ˆm † 60

Oe

…

ÿ

contributions of all elements connected with node m are summed.

…9b†

This set of ordinary di€erential equations are the Galerkin weak approximations of Equations (1) and (2) on the generical ®nite element Oe. Equations (9a)±(b) for Oe can be written in the following form: _ m ˆ Qe Aelm Cm ‡ B elm C l

l, m ˆ n, p, q

…10a†

T elm Cm ‡ Relm C_ m ˆ G el

l, m ˆ n, p, q

…10b†

where Cm and Cm are (3  1) vectors of the velocity potential and tracer concentration at element nodes, _ m and CÇ m are the (3  1) vectors of their correC sponding time derivatives. Alme and Tlme are, respectively, the (3  3) sti€ness matrices in ¯ow and in transport equations; Blme and Rlme are the corresponding (3  3) capacity matrices; Qle and Gle are the (3  1) vectors that account for source/sink terms and boundary conditions. It must be noted that in Equations (10a)±(b) the _ m , Cm and CÇ m are only the contriquantities Cm, C butions of Oe to the values of these quantities in the generic node m of the global ¯ow domain O. The true nodal values of these quantities at node m is obtained by the assembling procedure, by which the

It can be seen that from Equations (12a)±(12b), the lumping approximation is obtained by assuming that the time derivative of the velocity potential may be considered as a constant on Oe. So doing, we note that the capacity matrix becomes a diagonal matrix. In general, this procedure reduces the smearing but increases the overshoot. The vector accounting for source/sink terms and boundary ¯uxes is: … … ^ e do ÿ q~ ^ e ds Qel ˆ ÿ qw hN hN w l l Oe

ˆ ÿ

Ge

e

Lpq D qw …2hl ‡ hp ‡ hq † ÿ q~ …2hl ‡ hp † 12 6 w

…13†

Lpq being the length of the side pq of Oe. For an uncon®ned aquifer the capacity matrix is 8 De > > … lˆm < Sy 6 B elm ˆ Sy Nl Nm do ˆ …14a† e > Oe > : Sy D l 6ˆ m 12 By applying the lumping approximation, the capacity equation for an uncon®ned aquifer becomes: … De B elm ˆ Sy Nl do ˆ Sy dlm : …14b† 3 Oe Matrices Aelm and Belm are symmetric. GFE matrices of transport equation The elemental matrices equation, Equation (9b), are:

for

the

transport

POLF: 2D ®nite element model

T elm ˆ T e1lm ‡ T e2lm ‡ T e3lm ‡ T e4lm ‡ T e5lm ‡ T e6lm …15† where T e1lm

ing the time derivative of h must be modi®ed, i.e. T e5lm ˆ ne fr

e e  ^ ij @ N l @ N m do ˆ ne h …D11 bm bl ˆ ne hD e e 4D @ x @ x i j O

T e2nm

by

@ N em e ˆ vi h^ N do @xi l Oe

_^ e N e do hN l m

applying

the

lumping

…17†

T e5lm ˆ ne fr h_ m

… O

e

8 1 > > < ‰v1 …hl bl ‡ hp bp ‡ hq bq † ‡ v2 …hl cl ‡ hp cp ‡ hq cq †Š @ h^ e e 12 ˆ vi N m N l do ˆ > Oe @ x i > : 1 ‰v1 …hl bl ‡ hp bp ‡ hq bq † ‡ v2 …hl cl ‡ hp cp ‡ hq cq †Š 24 8 De > > < lne fr …3hl ‡ hp ‡ hq † 30 ^ e N e do ˆ ˆ lne fr hN l m e > D Oe > : lne fr …2hl ‡ 2hm ‡ hk6ˆl6ˆm † 60 …

_ in the con®ned case h=0, as a consequence: … _^ e N e do ˆ 0 T e5lm ˆ ne fr hN l m Oe

…20†

By using triangular linear elements the Darcy velocity is constant inside each element Oe, giving: … ^ e @ vi N e do ˆ 0 hN T6lm ˆ …21† m @xi l Oe For what concerns the capacity matrix: … ^ e N e do Relm ˆ ne fr hN l m Oe

8 De > > < ne fr …3hl ‡ hp ‡ hq † 30 ˆ e > D > : ne fr …2hl ‡ 2hm ‡ hk6ˆl6ˆm † 60

lˆm

…22a†

l 6ˆ m

by applying the lumping approximation we obtain: … e ^ e do ˆ ne fr D …2hl ‡ hp ‡ hq †dlm : Relm ˆ ne fr hN l 12 Oe …22b† The vector accounting for source/sink terms and boundary ¯uxes of pollutant is: … … ^ e do ÿ q~ ^ e ds G el ˆ ÿ qc hN hN d l l Oe

lˆm

…24a†

l 6ˆ m

approximation

N el do ˆ ne fr

…

T e4lm

ˆ ÿ

Oe

Te5lm

becomes:

1 ˆ …v1 bm ‡ v2 cm †…2hl ‡ hp ‡ hq † 24

T e3lm

…16†

…

…

8 De > > < ne fr …3h_ l ‡ h_ p ‡ h_ q † 30 ˆ e > D > : ne fr …2h_ l ‡ 2h_ m ‡ h_ k6ˆn6ˆm † 60

…

‡ D12 bl cm ‡ D21 bm cn ‡ D22 cm cl †

513

Ge

Lpq De qc …2hl ‡ hp ‡ hq † ÿ q~ …2hl ‡ hp † 12 6 d

…23†

For an uncon®ned aquifer only the terms contain-

De _ hm dlm 3

lˆm

…24b†

…18†

l 6ˆ m

lˆm

…19†

l 6ˆ m

All the others matrices remain unchanged with respect to the con®ned case. In all the previous equations, the terms with indices n, p and q represent quantities computed on the corresponding node of Oe; moreover, these terms must be taken in cyclical order n, p and q (e.g. for l = p, then p = q and q = n). The overlined quantities represent values averaged on the three node of Oe [e.g. h=(hn+hp+hq)/3]. The Kronecker operator dlm is used to indicate the diagonal structure of the capacity matrices Belm and Relm, when the lumping approximation is used. The parameters Ss, Sy, l, fr, ne, qw and qÄ d are considered to be constant inside each element Oe and along each side Ge. CVFE FORMULATION OF THE MODEL

As in GFE formulation, in CVFE formulation (Gottardi and Venutelli, 1992), the ¯ow domain O is discretized by triangular elements Oe. On each element Oe with nodes n, p and q, we assume that the variables c, c and h can be approximated by a linear combination of their nodal values, as shown in Equations (4). To each node of Oe, as shown in Figure 1, is associated a subdomain (subcontrol volume), on which Equations (1) and (2) are integrated. The subcontrol volumes Oen, Oep and Oeq associated, respectively, with the nodes n, p and q of Oe are obtained by drawing the three segments (en, ep, eq) connecting the barycentre G of Oe with the midpoint of each side (points N, P, Q) of the triangle. By integrating Equations (1) and (2) on the

514

G. Gottardi

subcontrol volume Oel (l = n, p, q), we obtain: … … … @ @c do …25a† ÿ …hvi † do ÿ hqw do ˆ a e @xi e e Ol Ol Ol @ t …

… @ @ …hJi † do ÿ …hvi c† do e @xi e @xi Ol Ol … … … @ …hc† ÿ hne lfr c do ÿ do : hqc do ˆ ne fr @t Oel Oel Oel

ÿ

…25b† By applying the divergence theorem to the ®rst terms of Equations (25a)±(b), and using Equations (3a)±(b), we obtain: "… # "… # @ N em _mˆ ÿ hkij ni ds Cm ‡ aNm do C @xj Gel Oel … … ÿ hqw do ÿ hq~ w ds …26a† Oel

S el

" … # "… # @ N em @ N em ^ ^ ÿ ne hDij ni ds Cm ‡ vi h do Cm @xj @xi Gel Oel "…

‡

vi N em e

Ol

# "… # @ h^ @ v i e ^ hN do Cm ‡ do Cm m @xi @xi Oel

# # " … … e e _^ ^ ‡ ne fr l hN m do Cm ‡ ne fr N m h do Cm "

" ‡ ne fr

Oel

… Oel

The structure of the elemental matrices for the ¯ow Equation (26a) are: … e ^ ij @ N m
ni ds Aelm ˆ ÿ hk e @xj Gl ˆ ÿ

1 ‰…k11 bm ‡ k12 cm †…hq2 ÿ hp2 † 2De

ÿ …k21 bm ‡ k22 cm †…hq1 ÿ hp1 †Š where

…

hp1 ˆ … hq1 ˆ

1 h^ dx 1 ˆ …5hl ‡ 2hp ‡ 5hq †…cq ÿ cl † …27a† 72 ep

1 h^ dx 1 ˆ …5hl ‡ 5hp ‡ 2hq †…cl ÿ cp † …27b† 72 eq

… hp2 ˆ … hq2 ˆ

ep

1 h_ dx 2 ˆ …5hl ‡ 2hp ‡ 5hq †…bl ÿ bq † …27c† 72

1 h^ dx 2 ˆ …5hl ‡ 5hp ‡ 2hq †…bp ÿ bl †: …27d† 72 eq

For a con®ned aquifer, the capacity matrix Belm is: … ^ e do B elm ˆ noSs hN m Oel

Oel

# … … e ^ ^ c do ÿ hN m do C_ m ˆ ÿ hq h^q~ d ds Oel

CVFE matrices of ¯ow equation

S el

…26b† where Gel and Sel are, as shown in Figure 1, respectively, the internal and external boundary of the subcontrol volume Oel (Gen=epÿeq, Gep=eqÿen, and Sen=enq+enp, Sep=epn+epq, Geq=enÿep e Sq=eqp+eqn), en, ep and eq being the segments connecting the barycentre G with the midpoint of each side (points N, P and Q). These segments are considered positive when oriented from the middle point of the side to the barycentre G of Oe. The ¯ow across the segments enp and enq may be di€erent from zero only for nodes on the boundary of O, whereas for nodes internal to O they contribute to the mass ¯ow across the boundary Sel (expressed by the second integral on the right side of Eq. (26b)) is zero. In Figure 3 is shown the control volume associated with each internal node of the mesh used to discretize the ¯ow domain O. From Equations (26a)±(26b), using for c, c and h the approximations shown in Equations (4), we obtain for each subelement Oen a set of di€erential equations in the same form as Equations (10a)± (10b).

8 > < …170hn ‡ 47hp ‡ 47hq † ˆ Bnn ˆ Bpn ˆ Bqn D ˆ Ss …47hn ‡ 23hp ‡ 14hq † ˆ Bnp ˆ Bpp ˆ Bqp 1296 > : …47hn ‡ 14hp ‡ 23hq † ˆ Bnq ˆ Bpq ˆ Bqq …28a† whereas for the capacity matrix, calculated by the lumping approximation, we obtain: … De e B lm ˆ Ss h^ do ˆ dlm Ss …22hl ‡ 7hp ‡ 7hq † 108 Oel …28b†

Figure 3. Control volume associated with node in triangular ®nite-element mesh O

POLF: 2D ®nite element model

…

Qel ˆ ÿ

Oel

qw h^ do ÿ

… S el

_ In the con®ned case h=0, as a consequence: … _^ e do ˆ 0 T e5lm ˆ ne fr hN m

q~ w h^ ds

Oel

De ˆ ÿ qw …2hl ‡ hp ‡ hq † 108 q~w ‰…3hl ‡ hp †Llnp ‡ …3hl ‡ hq †Llnq Š: 8

ÿ

…29†

For an uncon®ned aquifer, all the previous equations continue to hold valid, provided we take h = c ÿ zb. In this case the capacity matrix Belm is: 8 22 > > … Sy De l ˆ m < 108 e e B lm ˆ Sy N m do ˆ > Oel > : 7 Sy De l 6ˆ m 108 whereas, for the lumping case, we obtain: … De B elm ˆ Sy do ˆ Sy : 3 Oel

…30b†

CVFE matrices of transport equation With reference to Equation (26b), the elemental matrices for the transport equation are: T

e lm

ˆT

e 1lm

‡T

where T e1lm ˆ ÿ ne

…

e 2lm

‡T

e 3lm

‡T

e 4lm

‡T

e 5lm

…31†

ÿ …D21 bm ‡ D22 cm †…hq1 ÿ hp1 †Š

ˆ

…

…31a†

e

@Nm vi h^ do e @xi Ol

1 …v1 bm ‡ v2 cm †…22hl ‡ 7hp ‡ 7hq † 216

T e3lm ˆ

…

vi N em

Oel

where

…31b†

@ h^ do ˆ b‰v1 …hl bl ‡ hp bp @xi

‡ hq bq † ‡ v2 …hl cl ‡ hp cp ‡ hq cq †Š

… Oel

…

…31e†

e ^ m do ˆ D ne fr hN 1296 Oel

8 > < …170hn ‡ 47hp ‡ 47hq † ˆ R4nn ˆ R4np ˆ R4qn …47hn ‡ 23hp ‡ 14hq † ˆ R4np ˆ R4pp ˆ R4qp > : …47hn ‡ 14hp ‡ 23hq † ˆ R4nq ˆ R4pq ˆ R4qq …32a† by applying the lumping approximation: … De Relm ˆ ne fr h^ do ˆ ne fr …22hl ‡ 7hp ‡ 7hq †dlm : e 108 Ol …32b† The vector accounting for source/sink and boundary ¯uxes is: … … G el ˆ ÿ q~ d h^ ds qc h^ do ÿ Oel

ˆ ÿ

S el

e

q~ D qc …2hl ‡ hp ‡ hq † ÿ d ‰…3hl ‡ hp †Llp 108 8 …33†

In an uncon®ned aquifer all the previous equations continue to hold valid, provided we take h = c ÿ zb. For the matrix Te5lm we obtain: … De _^ e T e5lm ˆ ne fr hN do ˆ ne fr m 1296 Oel 8 …170h_ n ‡ 47h_ p ‡ 47h_ q † ˆ T5nn ˆ T5pn ˆ T5qn > > < …47h_ n ‡ 23h_ p ‡ 14h_ q † ˆ T5np ˆ T5pp ˆ T5qp > > : …47h_ n ‡ 14h_ m ‡ 23h_ q † ˆ T5nq ˆ T5pq ˆ T5qq …34a†

…31c†

8 22 > > lˆm < 216 bˆ > > : 7 lˆ 6 m 216

T e4lm ˆ ne lfr

Relm ˆ ne fr

‡ …3hl ‡ hq †Llnq Š:

e

^ ij @ N m
ni ds hD e @xj Gl

ne ˆ ÿ e ‰…D11 bm ‡ D12 cm †…hq2 ÿ hp2 † 2D

T e2lm ˆ

515

^ e do ˆ D ne lfr hN m 1296

8 > < …170hn ‡ 47hp ‡ 47hq † ˆ T4nn ˆ T4pn ˆ T4qn …47hn ‡ 23hp ‡ 14hq † ˆ T4np ˆ T4pp ˆ T4qp > : …47hn ‡ 14hp ‡ 23hq † ˆ T4nq ˆ T4pq ˆ T4qq …31d†

whereas, for the case of lumping approximation, we obtain: … De T e5lm ˆ ne fr h_^ do ˆ ne fr h_ m dlm : …34b† 3 Oel In all the previous equations the terms with indices n, p and q are quantities computed on the corresponding node of Oe; these terms must be taken in cyclical order n, p and q (e.g. for l = p, then p = q and q = n). Symbols with an over dot indicate the time derivative of the those quantities. The Kronecker operator dlm is used to indicate the diagonal structure of the capacity matrices Belm and Relm when the lumping approximation is used. The parameters Ss, Sy, l, fr and ne are considered to be constant inside each element.

516

G. Gottardi GLOBAL EQUATIONS e

Going from element region O to the whole ¯ow region O, by summing the contribution of all the elements (assemblage), in which the ¯ow region was discretized, the following set of ordinary di€erential equations is obtained: _ m ˆ Ql Alm Cm ‡ Blm C

l, m ˆ 1, . . . , nd

…35a†

Tlm Cm ‡ Rlm C_ m ˆ Gl

l, m ˆ 1, . . . , nd

…35b†

This set of equations is formally identical to that of Equations (10a)±(b), but now, Cm and Cm are the (nd1) vectors of the velocity potential and tracer _m concentration at the nd nodes of O, whereas, C and CÇ m are the (nd1) vectors of their corresponding time derivatives. Alm and Tlm are, respectively, the (ndnd) sti€ness matrices in ¯ow and in transport equations; Blm and Rlm are the corresponding (ndnd) capacity matrices; Ql and Gl are the (nd1) vectors that account for source/sink terms and boundary conditions. By discretizing in time the set of ordinary di€erential equations (Eqs (35a)±(b)) we obtain the following set of algebraic di€erence equations:   Bn B nlm n n Anlm ‡ lmn Cn‡1 C l, m ˆ 1, nd m ˆ Ql ‡ Dt Dtn m …36a†   Rn Rnlm n n C T nlm ‡ lmk C n‡1 m ˆGl ‡ Dt Dtn m n+1

Test case 1 This test case models the ¯ow of a pollutant into a semi-in®nite strip of saturated soil. The one dimensional form of the advection±di€usion equation governing pollutant ¯ow (Kinzelbach, 1986; Bear and Verruijt, 1987), written using the pore velocity (u = v/ne), is: @ c D @ 2c u @ c ˆ ÿ lc ÿ @ t fr @ x 2 fr @ x

n

where Dt =t ÿt and superscripts n and n + 1 indicate time levels. By solving, for each time step, Equation (36a), after having imposed in the boundary nodes either the Dirichlet or the Neuman boundary condition, we obtain the spatial distribution of the C, by which the Darcy velocities are computed using Equation (1a). By using the distribution of the Darcy velocities, so obtained, it is possible to assemble and solve Equation (36b), obtaining the space distribution of the pollutant concentration. This sequential procedure, that allows these two equations to be solved separately, is permissible due to the fact that, as previously pointed out, when the pollutant is supposed to behave as a tracer, the ¯ow and transport equations are disjointed.

TEST PROBLEMS

In order to show the use and the performances of the code, two test cases, for which the analytical solution is known, are presented here. In the ®rst test case, the results of the simulations for a con®ned system, with di€erent Peclet and Courant numbers, obtained by the GFE and CVFE

…36†

The solution of Equation (36) with the following initial and boundary conditions: c…x, 0† ˆ 0;

c…1, t† ˆ 0;

c…0, t† ˆ c0

is:

      c0 x ÿxg x ÿ utg=fr c…x, t† ˆ exp exp erfc p 2aL 2aL 2 2 aL ut=fr     xg x ‡ utg=fr ‡ exp …37† erfc p 2aL 2 aL ut=fr where

l, m ˆ 1, nd …36b†

n

methods, are compared on the basis of the massbalance errors and of the shifting from the analytical solutions. For the second test case only the shifting from the analytical solution is considered. The two test cases are the ones presented by Sun and Yeh (1983).

gˆ

p 1 ‡ 4laL fr =u

The Peclet and Courant numbers are: Pe ˆ

uDx Dx ˆ ; D aL

Cr ˆ

uDt Dx

…37a†

where Dx is the spacing of nodes along the x-axis and, by neglecting the molecular di€usion dm, the dispersion coecient becomes: D = aLu. The FE mesh used for this simulation is given in Figure 4. The following data were used: c0=10 kg/m3, u = 1, 2, 4 m/d, h = 10 m, Dx = 5 m, l = 0, fr=1, dm=0. The problem was solved by using the GFE and CVFE methods for di€erent values of Cr and Pe numbers. The di€erent values of Pe (5, 10, 15, 20, 25, 50 and 100) for the di€erent values of Cr (0.20, 0.4 and 0.8), were obtained by keeping constant the pore velocity u (1, 2 and 4 m/d) and assuming di€erent values for aL (1, 0.5, 0.33, 0.25, 0.20, 0.1 and 0.05 m). In order to use the same mesh for the di€erent pore velocities, the simulation times corresponding to the three values of Cr (0.2, 0.4 and 0.8) were, respectively, 50, 25 and 12.5 days. In this way, the results given for the di€erent Cr make reference to the same position of the pollutant concentration front. The results of the simulations are shown in Table 1. The parameters used to evaluate the reliability of the di€erent runs were, respectively, the mass-balance error

POLF: 2D ®nite element model

517

Figure 4. FE element mesh used for test case 1

eb%, and the two errors ef% and et% linked to the analytical solution. The mass balance error (eb%) is calculated as: eb % ˆ 100 

cpin ‡ cpi ÿ cpo ÿ cpp …cpi ÿ cpo †

ting''), takes into account the smearing of the numerical solution in comparison with that of the analytical one, and is calculated as

…38a†

ef % ˆ 100 

cpin being the initial cumulative pollutant, cpi and cpo, respectively, the cumulative pollutant in input and output from the system and cpp is the cumulative mass of pollutant present in the system. The mass balance error et% (we named it ``theoretical mass balance error'') is calculated as: …L …ca …x† ÿ c…x†† dx et % ˆ 100  0 … L …38b† ca …x† dx

nd 1X jca …i † ÿ c…i †j nd iˆ1 c0

…38c†

nd being the number of grid points used to discretize the ¯ow domain, ca(i) and c(i) are, respectively, the analytical and numerical values of pollutant concentration at the ith node of the mesh and c0 is the imposed pollutant concentration at node i = 1. This error takes into account the smearing of the numerical solution in comparisons with that of the analytical one and increases by the increase of the Pe and Cr numbers. The results of simulations at times t = 50, 25 and 12.5 d, for di€erent values of Cr (0.20, 0.40 and 0.80) and of Pe (5, 10, 15, 20, 25, 50, 100), are given in Table 1. In Table 1 we have indicated, respectively, with CV and FE the mass-balance errors

0

ca(x) and c(x) being the pollutant concentration distributions at time t, calculated respectively, by Equation (37) and by the numerical procedure. The ef% error (we named it ``reliability of ®t-

Table 1. Global errors in simulation of test case 1 for di€erent Peclet and Courant numbers Cr=0.20 aL

Pe

1

5

0.5

10

0.33

15

0.25

20

0.20

25

0.10

50

0.05

100

(%) eb ef et eb ef et eb ef et eb ef et eb ef et eb ef et eb ef et

CV ÿ3

5.73e 1.42 2.10 1.05eÿ2 1.93 1.05 1.22eÿ2 2.29 6.94eÿ1 1.20eÿ2 2.53 0.51 1.12eÿ2 2.69 4.09eÿ1 7.27eÿ3 3.25 1.95eÿ1 4.25eÿ3 3.56 9.31eÿ2

FE ÿ3

2.78e 1.33 2.11 6.18eÿ3 1.73 1.05 9.04eÿ3 1.96 6.97eÿ1 1.03eÿ2 2.13 0.52 1.07eÿ2 2.26 4.11eÿ1 8.74eÿ3 2.61 1.94eÿ1 5.36eÿ3 2.79 9.02eÿ2

Cr=0.40 CVL FEL ÿ2

1.68e 2.02 2.09 2.81eÿ2 3.05 1.03 3.13eÿ2 3.68 6.77eÿ1 3.21eÿ2 4.07 0.49 3.21eÿ2 4.34 3.90eÿ1 3.25eÿ2 4.99 1.75eÿ1 3.35eÿ2 5.35 6.84eÿ2

CV ÿ3

6.54e 2.37 2.10 1.13eÿ2 3.11 1.05 1.29eÿ2 3.52 6.94eÿ1 1.25eÿ2 3.77 5.16eÿ1 1.15eÿ2 3.94 4.08eÿ1 6.52eÿ2 4.38 1.95eÿ1 2.53eÿ3 4.64 9.37eÿ2

FE ÿ3

3.38e 2.28 2.10 6.67eÿ3 2.92 1.05 9.28eÿ3 3.28 6.98eÿ1 1.04eÿ2 3.51 5.19eÿ1 1.06eÿ2 3.68 4.11eÿ1 8.53eÿ3 4.10 1.94eÿ1 5.28eÿ3 4.31 9.11eÿ2

Cr=0.80 CVL FEL ÿ2

1.75e 2.85 2.10 2.31eÿ2 3.85 1.04 2.39eÿ2 4.35 6.83eÿ1 2.35eÿ2 4.72 5.05eÿ1 2.29eÿ2 4.98 3.98eÿ1 2.17eÿ2 5.61 1.85eÿ1 2.19eÿ2 5.94 8.00eÿ2

CV ÿ2

1.29e 3.89 2.09 1.47eÿ2 4.83 1.04 1.55eÿ2 5.30 6.91eÿ1 1.49eÿ2 5.59 5.13eÿ1 1.38eÿ2 5.79 4.10eÿ1 8.87eÿ3 6.28 1.90eÿ1 4.90eÿ3 6.53 9.20eÿ2

FE ÿ3

8.35e 3.78 2.10 9.27eÿ3 4.68 1.05 1.11eÿ2 5.13 6.97eÿ1 1.18eÿ2 5.42 5.19eÿ1 1.17eÿ2 5.62 4.10eÿ1 9.04eÿ3 6.10 1.90eÿ1 5.55eÿ3 6.35 9.21eÿ2

CVL FEL 2.88eÿ2 4.27 2.08 2.37eÿ2 5.29 1.03 2.02eÿ2 5.81 6.84eÿ1 1.74eÿ2 6.13 5.10eÿ1 1.53eÿ2 6.37 4.00eÿ1 1.12eÿ2 6.93 1.90eÿ1 9.89eÿ3 7.23 9.21eÿ2

518

G. Gottardi

Figure 5. Pollutant concentration fronts after 50 days with Pe=10

associated, respectively, with the CVFE and GFE methods. By CVL and FEL we indicate the massbalance errors obtained, respectively, by CVFE and GFE methods when the lumping approximation is imposed. It must be noted from the elemental matrices that, in the situation of thickness constant on Oe (hn=hp=hq), the CVL and FEL solutions are identical. From the results of simulations shown in Table 1 we see that: (1) The mass-balance errors (eb%) are small for all the three methods and depend little on the Pe and Cr numbers. Moreover, the GFE method has shown to give better results. (2) The et% error for constant Pe number does not depend on Cr and on the numerical method used.

(3) The ef% error, which is linked to the smearing of the numerical solution, increases with the increase of both the Pe and Cr numbers. In Figures 5 and 6 the numerical solutions obtained by the CVFE and GFE methods, with two di€erent values of Peclet numbers (Pe=10, 100), are compared with the corresponding analytical solutions. Figures 5 and 6 give a visual con®rmation of the above remarks about the role of the Peclet number and the application of the lumping procedure.

TEST CASE 2

This test case simulates the pollutant transport across a semi-in®nite fracture of a fractured porous medium. The analytical solution of this problem

Figure 6. Pollutant concentration fronts after 50 days with Pe=100

POLF: 2D ®nite element model

519

has been presented by (Tang and others, 1981). With reference to the FE mesh shown in Figure 7, we used for the simulation of this problem the following initial and boundary conditions: c…x 1 , x 2 , 0† ˆ 0; 



Figure 7. FE mesh used for test case 2

@c @x2

@c @x1

c…0, x 2 , t† ˆ c0 

x 2 ˆAC

ˆ

ˆ0

x 1 ˆCD





…42a†

@c @x2

…42b†

 x 2 ˆBD

ˆ0

…42c†

The following data were used for the simulation of this problem: fracture width 2b = 120 mm; in the matrix we used: ne=0.35, u1=u2=0, and a molecular di€usion coecient dm variable from 10ÿ10 to 10ÿ6 cm2/s. The parameters used for the fracture were: ne=1, aL=0.76 m, u1=0.75 m/day and u2=0. The dimensions (396 nodes and 700 elements) of the FE mesh shown in Figure 7, and used for this problem, are the same used by Sun and Yeh (1983), that is: AB = 2 cm (subdivided into 10 intervals: Dx1=0.006, 0.009, 0.015, 0.025, 0.045, 0.07, 0.13, 0.25, 0.45, 1 cm), and AC = 341 cm (subdivided into 35 intervals: Dx2=4(1), 3(1.5), 3(2.5), 15(4.0), 6.0, 8, 10, 16, 25, 30, 35, 40, 45, 50 cm).

Figure 8. Comparison among analytical (continuous lines) and numerical solutions (dots) of pollutant concentration along fracture after 4 days and for di€erent values of molecular di€usion coecient dm

520

G. Gottardi

In Figure 8 is shown the analytical solution (continuous line) of the pollutant distribution along the ®rst 100 cm of the fracture and after 4 days, calculated for di€erent values of the di€usion coecient dm of the porous medium. The corresponding solutions obtained by CVFE, GFE methods (dots) are very similar (not distinguishable in Fig. 8) and in good agreement with the analytical solution until dm>10ÿ10 cm2/s (dmmin). For dmRdmmin the pollutant concentration at the exit face of the fracture is no more negligible, and then the spatial discretization used does not respect the boundary condition expressed by Equation (42b).

CONCLUDING REMARKS

The CVFE and GFE formulations implemented in the numerical code POLF for the numerical integration of a two-dimensional ¯ow and pollutant transport problems have been presented. From the results of the numerical experiments the following main conclusions can be drawn: (1) The results of simulations obtained by GFE and CVFE methods show mass-balance errors of the same order of magnitude. (2) The GFE seems to be a little more mass conservative than the CVFE method. (3) As can be seen in Table 1, the application of the lumping procedure, for both the CVFE and GFE methods, does not improve any of the three mass-balance errors used to evaluate the goodness of the simulation.

REFERENCES Bear, J. (1979) Hydraulics of Groundwater. MacGraw-Hill Book Co., New York, 568 pp. Bear, J. and Verruijt, A. (1987) Modeling Groundwater Flow and Pollution. D. Reidel Publishing Co., Dordrecht, 414 pp. Chavent, G., Cohen, G. and Ja€re, J. (1984) Discontinuous upwinding and mixed ®nite elements for two-phase ¯ows in reservoir simulation. Computer Methods in Applied Mechanics and Engineering 47, 93± 118. Connor, J. J. and Brebbia, C. A. (1976) Finite Element Techniques for Fluid Flow. Newnes-Butterworths, London, 310 pp. Dick, E. (1983) Accurate Petrov±Galerkin methods for transient convective di€usion problems. International Journal for Numerical Methods in Engineering 19 , 1425±1433. Ferraresi, M. and Gottardi, G. (1990) Galerkin's ®nite element matrices for transport Phenomena in porous media. Proceedings of the 8th International Conference on Computational Methods in Water Resources, Venice, Italy, June 11±15. Computational Mechanics Publications, Southampton Boston, pp. 271±276. Friend, J. J. (1975) Groundwater Pollution. Elsevier, New York, 330 pp.

Gottardi, G. and Dall'olio, D. (1994) Modello per il trasporto di inquinanti in falde acquifere. Geoengineering Environment and Mining 2, 39±45. Gottardi, G. and Venutelli, M. (1992) Moving ®nite element model for one-dimensional in®ltration in unsaturated soil. Water Resources Research 28(12) , 3259± 3267. Gottardi, G. and Venutelli, M. (1994) One-dimensional moving ®nite-element model for solute transport. Ground Water 32(4), 645±649. Gray, W. G. and Pinder, G. F. (1976) An analysis of the numerical solution of the transport equation. Water Resources Research 12(3), 547±555. Heinrich, J. C., Huyakorn, P. S., Zienkiewicz, O. C. and Mitchell, A. R. (1977) An ``upwind'' ®nite element scheme for two dimensional convective transport equation. International Journal for Numerical Methods in Engineering 11, 131±143. Hughes, T. J. R. (1978) A simple scheme for developing upwind ®nite element. International Journal for Numerical Methods in Engineering 12, 1359±1365. Huyakorn, P. S. and Nilkuha, K. (1979) Solution of transient transport equation using upstream ®nite element schemes. Applied Mathematical Modelling 3, 7±17. Kinzelbach, W. (1986) Groundwater Modelling: An Introduction with Sample Programs in BASIC. Elsevier, Amsterdam, 333 pp. Narasimham, T. N. and Witherspoon, P.-A. (1976) An integrated ®nite di€erence method for analyzing ¯uid ¯ow in porous media. Water Resources Research 17(6), 1665±1675. Sun, N.-Z. and Yeh, W. W.-G. (1983) A proposed upstream weight numerical method for simulating pollutant transport in groundwater. Water Resources Research 19(6), 1489±1500. Thomas, J. R. H. (1978) A simple scheme for developing ``upwind'' ®nite elements. International Journal for Numerical Methods in Engineering 12, 1359±1365. Tang, D. H., Frind, E. O. and Sudicky, E. A. (1981) Contaminant transport in fractured porous media: analytical solution for a single fracture. Water Resources Research 17(3), 555±564. Yeh, G. T. (1986) An orthogonal-upstream ®nite element approach to modeling aquifer contaminant transport. Water Resources Research 22(6), 952±964. Yu, C. C. and Heinrich, J. C. (1987) Petrov±Galerkin method for multidimensional, time dependent, convective±di€usion equations. International Journal for Numerical Methods in Engineering 24, 2201±2215.

APPENDIX A Description of the Program The computer program POLF is written in FORTRAN 77 and consists of a main and 34 subroutines. Two test cases are provided complete of input and output data. The listing of the program along with the user manual and the data ®les for running the two test cases are available by anonymous FTP from IAMG.ORG. In the following a brief description of the subroutines of the program is presented. main: The main program reads some control data, opens the ®les, and drives the ¯ow of the model by calls to the subroutines. analit: This subroutine computes, when kana = 1, and for the node speci®ed by the user, the analytical solution of the pollutant transport equation by using Equation (37). ass1: For each time step, this subroutine sets to zero the global matrix and the right-hand side vector and assembles the global ¯ow equation by calling, following the directive

POLF: 2D ®nite element model of the parameter kmeth1, one of the following subroutines: elm1-cv, elm1-fe. ass2: For each time step, this subroutine sets to zero the global matrix and the right-hand side vector and assembles the global pollutant transport equation by calling, following the directive of the parameter kmeth2, one of the following subroutines: elm2-cv, elm2-fe, dlm2-cv and dlm2-fe. balance: When called, this subroutine computes the mass-balance errors of both water and pollutant on the basis of the initial and present mass of water and pollutant and on the basis of the cumulative mass of water and of pollutant that have entered and left the ¯ow domain. For both water and pollutant three types of mass-balance errors (err1, err2 and err3) are computed: err1 ˆ 100  Dm=mf err2 ˆ 100  Dm=mp err3 ˆ 100  Dm=…mi ÿ mo † where Dm ˆ mf ‡ mi ÿ mo ÿ ma ÿ md ÿ mp and where mf is the initial mass, mi and mo are, respectively, the cumulative-masses which entered and left the system in the considered time, ma and md the cumulative mass absorbed and decayed and mp the mass present in the system at the considered time. bnd1-cv: Imposes the Dirichlet condition of pre®xed total head for solving the ¯ow equation by the CVFE formulation. bnd1-fe: Imposes the Dirichlet condition of pre®xed total head for solving the ¯ow equation by the GFE formulation. bnd2: This subroutine imposes the Dirichlet condition of pre®xed pollutant concentration for solving the pollutant transport equation by both the CVFE and GFE formulations. cmnt: Reads from logical unit 14 the commented user input data ®le and copies, on logical unit 15, all the lines of this data ®le which are not empty and not preceded by a semicolon. This allows the user to put as many comment lines he wishes to facilitate the preparation of the input data ®le. convert: Converts input data from ®eld to SI system of measure. darcy: Computes, by an implicit way, the Darcy velocities for each node of the mesh by using the total heads. delm2-cv: Is used in the integration of the pollutant transport equation for computing the elemental conductivity matrix, Tlm, the capacity matrix, Rlm, and the righthand side vector, Gl, using the CVFE formulation. This is done by using for each node the Darcy velocities computed implicitly by the subroutine darcy. Moreover, it assembles the global matrix and right-hand side vector for the computation of the pollutant concentration at each node of the mesh. delm2-fe: Is used in the integration of the pollutant transport equation for computing the elemental conductivity matrix, Tlm, the capacity matrix, Rlm, and the righthand side vector, Gl, using the GFE formulation. This is done by using for each node the Darcy velocities computed implicitly by the subroutine darcy. Moreover, it assembles the global matrix and right-hand side vector for the computation of the pollutant concentration at each node of the mesh. elm1-cv: Is used in the integration of the ¯ow equation for computing the elemental conductivity matrix, Alm, the capacity matrix, Blm, and the right-hand side vector, Ql, using the CVFE formulation. Moreover, it assembles the

521

global matrix and right-hand side vector for the computation of the total heads at each node of the mesh. elm1-fe: Is used in the integration of the ¯ow equation for computing the elemental conductivity matrix, Alm, the capacity matrix, Blm, and the right-hand side vector, Ql, using the GFE formulation. Moreover, it assembles the global matrix and right-hand side vector for the computation of the total heads at each node of the mesh. elm2-cv: Is used in the integration of the pollutant transport equation for computing the elemental conductivity matrix, Tlm, the capacity matrix, Rlm, and the right-hand side vector, Gl, using the CVFE formulation. Moreover, it assembles the global matrix and right-hand side vector for the computation of the pollutant concentration at each node of the mesh. elm2-fe: Is used in the integration of the pollutant transport equation for computing the elemental conductivity matrix, Tlm, the capacity matrix, Rlm, and the right-hand side vector, Gl, using the GFE formulation. Moreover, it assembles the global matrix and right-hand side vector for the computation of the pollutant concentration at each node of the mesh. exitt: Is called by the subroutine gdata for printing the error messages during the reading of the input data ®le. ¯x1: Computes the total volumes of water entering and leaving the system through the external boundaries, through the wells and the distributed sources. ¯x2: Computes the total mass of pollutant entering and leaving the system through the external boundaries, through the wells and the distributed sources. gdata: Reads in free format the user input data ®le. If the automatic grid generation is required (krete = 1, 2, 3, 4), this subroutine calls the subroutine grete. grete: Allows the automatic generation of the four types of meshes shown in Figure 9. The input parameter (krete = 1, 2, 3, 4) is used for choosing the type of mesh. gsolve: Solves, by Gauss elimination technique without pivoting, a set of linear equations whose matrix is nonsymmetric and band structured. In order to save memory, the band matrix is stored in a rectangular array. gsolvs: Solves, by Gauss elimination technique without pivoting, a set of linear equations whose matrix is symmetric and band structured. In order to save memory,

Figure 9. Four types of FE grids that can be generated automatically

522

G. Gottardi

only the main diagonal and the upper diagonals of the band matrix are stored in a rectangular array. openf: Opens the ®les according the user speci®cations. pdata: Prints the input data, according to the print options speci®ed by the user. report: Prints, at pre®xed times, a report showing the water heads and pollutant concentrations at all nodes of the mesh, input and output volume of water and mass of pollutant and mass balance errors of water and pollutant. rread: Reads the data in the case of a restart run. rewrite: Saves at intervals of time speci®ed by the user the data necessary for a restart run. test: Checks the convergence of the iterative procedure used in solving the set of algebraic equations.

timer: Is a system-dependent intrinsic subroutine to compute execution CPU time. update: This subroutine reads the data that may be changed during the run time, that is: the step control parameters, the data specifying wells behaviour and the sources of pollutant. velox: Computes, by an explicit way, the components of the Darcy velocity for each element using the heads at each node of the element. vol: Computes the total volume of water and mass of pollutant present in the system on the basis of the water depth and of the geometry of each ®nite element used to discretize the ¯ow domain.