- Email: [email protected]
A sensitivity analysis of ‘2D-FED’, a model for seawater encroachment in leaky coastal aquifers
Journal of Hydrology, 118 (1990) 343-356
343
Elsevier S c i e n c e Publishers B.V., Amsterdam - - Printed in The N e t h e r l a n d s
[2] A S E N S I T I V I T Y A N A L Y S I S OF '2D-FED', A M O D E L FOR S E A W A T E R E N C R O A C H M E N T IN LEAKY C O A S T A L A Q U I F E R S
MOHSEN M. SHERIF ', VIJAY P. SINGH ' and ABDELWAHAB M. AMER 2
1Department of Civil Engineering, Louisiana State University, Baton Rouge, LA 70803-6405 (U.S.A.) 2Kuwait Institute for Scientific Research, 13109 Safat (Kuwait) (Received February 21, 1989; accepted after revision September 11, 1989)
ABSTRACT Sherif, M.M., Singh, V.P. and Amer, A.M., 1990. A sensitivity analysis of '2D-FED', a model for seawater encroachment in leaky coastal aquifers. J. Hydrol., 118: 343-356. In a previous study, a two-dimensional finite element model for dispersion (2D-FED) was developed and verified for some previously published problems. The model was applied to predict the shape of equiconcentration and equipotential lines in coastal aquifers. In this paper, the 2D-FED model is used to study the effect of several parameters on the simulated dispersion zone and flow pattern in leaky coastal aquifers. The aquifer thickness greatly affects the dispersion zone as well as the equipotential lines. Except for certain ranges, the effects of isotropic and anisotropic permeability coefficients and anisotropic dispersivities are relatively small. Cyclic flow is found near the shore boundary in many cases. Selection of a proper mesh size especially near the seaside is found to be very important to avoid instability problems.
INTRODUCTION
The well-knownGhyben-Herzberg(1889-1891) static equilibriumprinciple (Herzberg, 1901) was the first quantitative formulation for the extent of saltwater intrusion into coastal aquifers. Since that time, many studies have been made in this area. Some investigators (Senio, 1951; Henry, 1960; Shamir and Dagan, 1971; Mualem and Bear, 1974; Sa da Costa, 1981) assumed an abrupt interface between freshwater and saltwater. In reality, however, the freshwater and saltwater merge in a dispersion zone in which concentration and hence the fluid density vary. Where this zone is narrow, the sharp interface approach may provide an appropriate simulation, but an extensive zone of dispersion will alter the flow pattern. Therefore, it is prudent to use a simulation mode] that accounts for the hydrodynamic dispersion effects. One of the first experimental studies made on dispersion was by Kitagawa (1934). Beran (1955) distinguished between three cases of flow: (i) when the effects of molecular diffusion and randomness of the flow pattern are negligible; (2) when the effects of molecular diffusion are dominant; (3) when
0022-1694/90/$03.50
© 1990 Elsevier Science Publishers B.V.
344
M.M.SHERIFETAL.
the randomness of the flow pattern is of equal importance to molecular diffusion in the mixing process. Ogata and Banks (1961) presented a complete solution for the unsteady one-dimensional convective~iispersion equation. It was shown that the simplified form of the equation could be used provided that the dispersion coefficient was small and the region near the tracer source was not considered. Pinder and Cooper (1970) investigated the movement of saltwater in homogeneous confined aquifers. It was assumed that a constant freshwater flux enters the aquifer along the vertical boundary at the land side and a constant saltwater head at the sea side was maintained. The method of characteristics was used to solve the transport equation and the alternating direction iterative finite difference procedure was used to solve the groundwater flow equation for two-dimensional problems. Segol et al. (1975) used the Galerkin finite element method to solve the set of non-linear partial differential equations describing the movement of a saltwater front in a coastal confined aquifer. A layered aquifer was modeled either with a functional representation of permeability or with a constant value of permeability over each element. Umari et al. (1979) addressed the problem of identifying unknown aquifer dispersivities in two-dimensional transient groundwater contaminant transport from given observations on the concentration field. This inverse algorithm became the solution of the contaminant transport problem. The finite element method in conjunction with finite differencing was used to discretize the governing differential equations. Anand et al. (1980) indicated that the selection of a proper mesh size is extremely important to enhance convergence. It was found that solution instability could occur owing to high velocities near the sea side or because of the presence of sharp concentration gradients in a region away from it. For certain values of Peclet number and discharge parameters, saltwater intruded far enough inland that an aquifer length of twice its depth could not be considered adequate for boundary condition of freshwater concentration to be imposed on the landside. Several lengths of aquifers using appropriate finite element grids were used to illustrate this effect. Anand and Pandit (1982) studied the transient conditions of saltwater intrusion in coastal aquifers for a wide range of field parameters. The influence of selected time sizes on computational effort was invstigated. The movement of the saltwater front was found to be far more rapid in the beginning than near the steady state. Volker and Rushton (1982) obtained the location of the abrupt interface in homogeneous isotropic aquifers using the boundary integral method. They compared the location of the interface with the distribution of salt in the dispersion zone where the two fluids meet as determined by a finite-difference solution of the flow and convective~ispersion equations. They found that, for a constant dispersion coefficient, the interface becomes more diffused and migrates inland as the fresh water discharge to the sea decreases. Pandit and Anand (1984) did a parametric study on an idealized confined coastal aquifer. Their results indicated that if the selected region size was not long enough, the flow boundary condition at the freshwater face was not met. The depth of the
SENSITIVITY ANLYSIS OF 2D-FED MODEL
345
aquifer influenced the extent of saltwater intrusion into the aquifer. When the value of the vertical permeability coefficient, K=, was changed, the equiconcentration lines did not change significantly. Huyakorn et al. (1987) employed the picard sequential solution algorithm to simulate the saltwater intrusion problem in single and multiple coastal aquifers. Simulation results were compared with several previously published solutions. Convergence difficulties were encountered when dealing with more complex three-dimensional simulations. Sherif et al. (1988) developed a two-dimensional finite element model for dispersion (2D-FED) to simulate the same problem in confined and leaky aquifers under steady state conditions. The Galerkin technique was employed with computational scheme based on the relaxation method. The model was verified on two previously published problems. Cyclic flow was predicted at the shore boundary. This paper applies the 2D-FED model (Sherif et al., 1988) to study the effect of some hydraulic and geometric parameters on the extent of saltwater encroachment into coastal leaky aquifers and on the shape of equipotential lines. The effect of the element size on improving the convergence behavior is also investigated. AQUIFER CHARACTERIZATION
A domain of 400.0 x 40.0m, represented by a uniform mesh with 1280 triangular elements and 697 nodes, is considered. The rationale for simulating such a small region is that a bigger domain may require a non-uniform grid with smaller elements near the shore boundary to avoid numerical problems. For a wide range of the same parameter, it may be necessary to reconstruct a finer grid to enhance numerical convergence. Parameter values are based on the Nile Delta aquifer in Egypt (Farid, 1985). The aquifer porosity, n, is considered to be 0.3. Hydraulic conductivity in the x and z directions (Kxx and Kzz) are equal to K0 = 150.0m day -1. The upper semi-pervious layer (aquitard thickness is 10.0m), vertical hydraulic conductivity for the upper semi-pervious layer, K'z, is equal to 0.05 m day 1. There are no data available about the dispersivity for the Nile Delta aquifer, but Volker and Rushton (1982) indicated that values of longitudinal dispersivity aL that have been reported in field investigations range from 6 to 150m. Lateral dispersivity, aT, has usually been taken to be 10-30% of the longitudinal value (Pickens and Grisak, 1980). In laboratory experiments, longitudinal dispersivities of 0.1 10cm were found for different types of materials (Klotz, 1973). Longitudinal dispersivities in the field are usually much larger than those in laboratory experiments (Lenda and Zuber, 1970). The main reason for this is the influence of small-scale inhomogeneities in the permeability of the aquifer, which means the onset of macrodispersion. Owing to this scale behavior of dispersion, no generally valid values for longitudinal dispersivities can be given (Kinzelbach, 1986). Values between 0.1 and 500m can be found in the literature for the porous aquifer; Kinzelbach took a L 100 m and aT = 10 m in =
346
M.M. SHERIF ET AL.
the case of one permanent pollutant source of constant strength in a flow field parallel to the x-axis. In light of this discussion, longitudinal and transverse dispersivities % and a T are set equal to 100.0 and 10.0 m, respectively. The water table is set to be 6.0 m above sea level at the land side and 4.0 m at the sea side; between these two nodal points, the water table is assumed to be a straight line having a gradient of 0.05. Therefore, the water table is defined for nodal points on the upper boundary. The piezometric head is set to be 5.80 m at the land side and 3.9 m at the sea side; between these two nodal points, the piezometric head, like any node inside the domain, is not known a priori. Therefore, there will be a leakage flux into the aquifer through the upper semi-pervious layer (aquitard) all over this boundary. The aquitard storage is negligible and the leakage flux is varied from one iteration to another, depending on the difference in levels between the free water table which is constant, and aquifer piezometric head which is varied. On the seaward boundary, where the aquifer is exposed to the open sea, a mixed boundary is used, i.e. on the segment where the fluid leaves the aquifer (window) the concentration gradient is zero and on the segment where the seawater enters into the aquifer, the concentration is equal to seawater concentration C,. The exact length through which saltwater and mixed water enter and leave the domain is calculated by checking velocity directions, as shown in Fig. 1. This boundary is subject to hydrostatic pressure distribution, 7h. Over the window, the specific weight 7 is a function of the mixed water concentration that leaves the system. Below the window, the specific weight is equal to that of seawater 7,. The numerical convergence criterion is set equal to 0.001. Flow and boundary conditions are given in Fig. 1. SEA SIDE
LAND SIDE ~/// ~.,_ ~ t/ ,///, , i// ~,./ - - / ~ / / / / ////////,.//,//
"4./O/I~//,/'.PIZOMETRIC " '.."7/ , I / f l .
4..V ~ T ~c = 0 ~--~
xi.'v
"If
,q~..,~..
.~
WINDOW
j_ (4)
c = cs e =~h
.~
--
~l
I
I I I I I
.l~.'i'..,li//.i
wl Wl
SUBDOMAIN
"/'///////1"/
J
J~
lO0-O m
~1--]-
,~
[
I
(Z)
AQUIFER
I
".1"/
I00.0
m
V¢ l ~ II
~c = 0 ~n I~
•
i.i'ili.'ll''i~J~
[I l I l I I
(3)
l~//'/'//'.,',//., ,a~-/--//,
AQU|TARtD / / Q WATER TABLE :/d 6.0m
)
--'~
-"~
~
~///
4 .A.~iiy
1
~k..,.
SEA
iv I ~1
HEAD
.A.'I
~iS~l/////rf////~.l ,. / / / / / / . / /
"11/i///////////////
q.n= 0 ~:
I00.0 m
I I I I I I
I
~
"¥
4~ ql-- FRESHWATER ~
: ~
"///////
IMPERMEABLE _-It_-
400.0 m
Fig. 1. Boundary conditionsand subdomains in the studied region.
c = cF
~oo.o
m
_I, m]
P = ~f '
SENSITIVITY ANLYSIS OF 2D-FED MODEL
347
GOVERNING EQUATIONS AND SOLUTION TECHNIQUE
The governing equations describing the saltwater intrusion phenomena and the flow pattern in coastal aquifers in steady state conditions are as follows. (1) Darcy equation for groundwater flow: q =
- K(V~ + prVz)
(1)
where K is the coefficient of permeability (LT-1), ~h is the freshwater head (L), Pr is the relative density and is equal to p/pf - 1 (dimensionless), with p being the fluid density (ML 3), and pf being freshwater density (ML 3) and equal to 1000 kg m- 3. (2) The mass balance equation for the fluid: V.pq
= 0
(2)
(3) The mass balance equation for the dissolved salt or the hydrodynamic dispersion equation: nV.(Dh.VC)
-
V'qC
=
0
(3)
where n is the porosity (dimensionless), C is the concentration (ML-3), and D h is the hydrodynamic dispersion coefficient (L2T 1) which can be expressed as follows (Bear, 1979): Dxx =
V:
aL ~
V:
Dxz =
Dzx
=
V~2
D*
V:
D*
+ aT-~-~ +
(aL--aT)--
(3.1)
vxE
where aL and aT are the longitudinal and transverse dispersivities (L), Vx and Vz are the pore velocities in the x and z directions (LT-1), 1~ is the resultant pore velocity, and D* is the molecular diffusion coefficient (L2T - ' ) which is considered to be very small and can be neglected. (4) The constitutive equation relating fluid density to salt concentration is: p
=
pf + a ( C -
Cf)
(4)
where a is a known constant (dimensionless) and equal to 0.7246 for seawater having a density, Ps, of 1025.0 kg m- 3 and a concentration, C,, of 35.0 kg m 3. Cf is the freshwater concentration (ML 3) and is equal to 0.5kgm 3. The above four equations are combined together to form two non-linear partial differential equations in two unknowns, namely, the concentration and the equivalent freshwater piezometric head. The Galerkin finite element technique is adopted to solve these two non-linear equations with the appropriate boundary conditions. Triangular elements are used with linear inter-
348
M.M. SHERIF ET AL.
polation functions. The seaward boundary is updated in each iteration by checking velocity directions at this boundary to locate the exact length through which saltwater and mixed water enter and leave the domain. The computational scheme is based on local linearization of nonlinear terms. A complete discussion of the governing equations, boundary conditions and solution technique is given by Sherif et al. (1988). SENSITIVITY ANALYSIS The 2D-FED model is used to study the effect of several parameters on the width of the dispersion zone, and on the shape of the equipotential lines in a leaky coastal aquifer. The width of the dispersion zone (the extent to which the influence of hydrodynamic dispersion is expected) is measured by the length Lo defining the distance along the bottom of the aquifer from the sea boundary to isochlor line 1.0. The shape of the equipotential lines is measured by the length L' defined as the distance along the bottom of the aquifer from the sea boundary to equipotential line 4.7, which is found to be the first line that intersects the bottom boundary. Equipotential lines 4.6 and less did not meet this boundary in most of the runs. These lengths, L D and L', are normalized with respect to the total length of the aquifer LT which is 400.0m. The variation of the dimensionless ratios (LD/LT) and L'/LT) are studied with various parameter values. For each group of runs, one parameter is varied, and all others remain at the standard values mentioned above. Figure 2 presents the basic run to which other runs are compared.
Effect of isotropic permeability coefficient There are four different groups of runs (Table 1). In Group A, the coefficient of permeability is defined as: Kxx = Kzz = 50.0, 100.0, 150.0, 200.0, 250.0m day -1. It is found that as the isotropic permeability coefficient increases, the length of the dispersion zone also increases while the appearance of cyclic flow at the sea side decreases. Cyclic flow occurs when the velocity field at the sea boundary is rotational in character. Saltwater traverses away from the sea at the bottom of the aquifer and then circulates back into the sea through the upper segment of the boundary. The length of the dispersion zone and the equipotential lines is most sensitive to isotropic permeability when the permeability coefficient lies between 150.0 and 200.0m day -1, otherwise its effect is relatively small, as shown in Fig. 3(a). Increasing the isotropic permeability from 150.0 to 200.0 m day 1 increased the width of the dispersion zone from 164.0-195.0m (about 19%). Kawatani (1980) concluded that the existence of cyclic flow at the sea side prevents saltwater from intruding deeply into the aquifer. This result is confirmed in this study. In run 1 where Kxx = Kzz = 50.0m day -1, there is a strong cyclic flow at the sea side, as shown
349
SENSITIVITY ANLYSIS OF 2D-FED MODEL
BASIC RUN L :/.O0.0m
b = 40.0 m
n=0.3
K =K = K = 150.Orn day xx
zz
~0 SEA SIDE
LAND SIDE 2O
Q 1,
L 100
{,, 200
a) EQUI-CONCENTRATION
SEA SiDE
300
400
LINES.
t30 tttL 20
LAND 5~DE
t
~
100t
2OO
3OO
zOO
b) E~JI-POT~NTIAL LINES. Fig. 2. Basic run-run 3. G r o u p A.
in Fig. 3(b), the length of the dispersion zone, LD, is 152.0m. On the other hand, in run 5 where K=~ = Kz~ = 250 m day-~, the cyclic flow is not as strong as in run 1 and the saltwater intruded the bottom of the aquifer to a distance of 200.0 m as shown in Fig. 3(c).
Anisotropy of permeability coefficient In Group B, the effect of anisotropy of the permeability coefficient is studied. To that end, K~x is kept constant at 150.0m day-1 and Kz~ is varied as 25.0, 50.0, 100.0, 200.0, 250.0 m day- ~(in nature, K~z is always smaller than K=). The length of the dispersion zone, LD, decreases rapidly as (KJKx=) increases. After a certain ratio of (Kz~/Kx=),the anisotropic permeability coefficient has no effect on the dispersion zone and on the equipotential lines, as shown in Fig. 4(a). Decreasing the vertical hydraulic conductivity, K~, from 150 to 25.0 m day increased the width of the dispersion zone from 164.0 to 276.0 m (about 68%). It should also be noticed that the gradient in run 1, as shown in Fig. 4(b), is almost
350
M.M. SHERIF ET AL
TABLE 1 C o m p u t e r r u n s for s e n s i t i v i t y a n a l y s i s
Group A (Permeability effect)
Group B (Non-isotropic permeability effect)
R u n No.
Kxx (m d 1)
Kz~ (m d 1)
1 2
50.0 100.0
50.0 100.0
3 4 5
Ko 200.0 250.0
Ko 200.0 250.0
K0
25.0 50.0 100.0
1 2 3 4 5
Group C (Non-isotropic dispersivity effect)
1 2 3
~L (m)
aT (m)
b (m)
aLo
~L0
bo
aLo
aT0
b0
200.0 250.0
Ko
Ko
aLo
4 5 Group D (Aquifer t h i c k n e s s effect)
5.0 20.0 50.0 75.0 100.0
1 2
3 4
bo
20.0 48.0
K0
Ko
aLo
aT0
60.0 80.0
K 0 = 1 5 0 . 0 m d -1, aLo = 100.0m, aT0 = 10.0m, a q u i f e r t h i c k n e s s , b 0 = 40.0m.
constant except near the shore boundary. In run 5, the gradient is varied gradually from mild at the land side to steep near the sea side, as shown in Fig.
4(c). Transverse dispersivity In Group C, the longitudinal dispersivity aL is kept constant at 100.0m, while the transverse dispersivity is assigned values 5.0, 20.0, 50.0, 75.0 and 100.0 m. It is found that the shape of isochlor lines and equipotential lines are sensitive to anisotropic dispersivities if the ratio (aW/~L) is less than 0.2. Otherwise, the effect of this ratio is negligible. Decreasing the transverse dispersivity, aT, from 10.0 to 5.0m increased the width of the dispersion zone from 164.0 to 224.0m (about 36%). Figure 5 represents the effect of anisotropic dispersivity.
351
SENSITIVITY ANLYSIS OF 2D-FED MODEL
t x • 4oo.om
,,,4o.om
• .o.8
x e -150.0 m / d l 7
LT
j,
~ L T
t
o .ec
o 4e
~ o - - ~ - - - - - ~ - - - o
(
)
o. am
o .zo
~
0.15
CT
0 )VARIATION
j,
L.
OF
( L O / L T) AND
.;"
( L/ / L T ) W i T H
~
b ) EQUI - C O N C E N T R A T I O N RUN I. GROUP A.
AND
L"
(Kxx
/ K O)
*
EQUI- POTENTIAL
L I N E S , Kxx =50.Ore/DAY
LD
¢)EQUI-CONCENTRATION RUN 5 . GROUP A.
AND
EQUI-
POTENTIAL
LINES,Kxx
:ESaOm/OAY
Fig. 3. Permeability effect. Aquifer thickness in hydraulic contact with saltwater Finally, in Group D, the aquifer thickness, b, is varied between 20.0, 48.0, 60.0, 80.0 m, while the length of the studied region Lw is kept at 400.0 m. The same number of elements and nodes are used for all runs. It is found that the aquifer thickness, b, greatly affects the dispersion zone as well as the equipotential lines. The length of the dispersion zone, LD, rapidly increases with the increase of the aquifer thickness as shown in Fig. 6(a). When the aquifer depth is only 20.0 m, cyclic flow does not exist at the shore boundary and the mixed water finds its way out from the aquifer to the sea through the entire length of the boundary. Although there is no flow from the saltwater into the aquifer, salt ions still encroach the aquifer up to some limited distance ( ~ 85.0 m). Salt ions moved from higher concentrations (seawater) to lower concentration inside the aquifers under the effect of dispersion. As the depth of the aquifer
352
M.M. SHER1F E T AL.
LT
lrL
-400.Om
I% +'lOOm
,IOO.O
n = 03
KellBOO~/dly
.__T+eL • --,01
m x x • m,
~P
o.~o
+----
•
~"
o .os o.so o~4SoAo
~
o.al o.lc
old. T I
e.+
O.18
AP ~'~
---
~'e--------
e"~
0,2C
- - e t L" )
0111 OK i
o)
•
,
'
VARIATION
'
~
'
OF
'
t
{L D /L
,
.
+
,
T ) AND
i
.
,
,
(L ¢ /L T )
, --
(~Zz)
WITH
( Kzz/K
0 )
Lo b ) EQUI - CONCENTRATION RUN
,- D
t.
GROUP
AND
EQUI-POTENTIAL
L IN E $ , K Z Z = 2 5 . 0
m/
DAY
E.
..+
-+--~"' +'+
c ) EQUIRUN
CONCENTRATION 5-
GROUP
ANO
EQUI-
POTENTIAL
LINESp
KZZ
=250+0
m/DAY
B.
Fig. 4. A n i s o t r o p i c p e r m e a b i l i t y
effect.
increases, cyclic flow becomes stronger at the shore boundary. When the depth of the aquifer is 80.0 m, the length of the dispersion zone is ~ 270.0 m and a very strong cyclic flow exists near the shore, as shown in Fig. 6(c). Equipotential line 4.7 does not intersect with the aquifer bottom. Volker and Rushton (1982) concluded that the relative density, Pr = P/Pf - 1, of seawater affects the pressures at lower levels in the aquifer near the seaward boundary, and hence influences the length of intrusion into the aquifer. This result is confirmed in this study. As the aquifer depth increases, the effect of the relative density becomes more significant. For example, an aquifer depth of 100.0m and a relative density of 0.025 causes an additional 2.5 m pressure head at the bottom of the aquifer compared with an aquifer filled with freshwater. However, it should be noted that the relative density at the shore boundary will be 0.025 only over the segment where seawater enters into the aquifer (below the
353
S E N S I T I V I T Y A N L Y S I S O F 2D-FED M O D E L
LT
• 4oo.om
k,,40.O
- ~ . ,ooo . . . .
m
A .O
5
Kez I5OOm
/clay
"1
• Kzz - . .
i O.OB o iso 0.4B
C~
~30
O,2~
l't------
o ( L.~L ) LT
C~
- - . - O - -- - - - - I P - - - - - - e l c'~ )
OZC
OIc
a
O)
d
~.
d
VARIATION
o
O
d
=
-
-"
(L D / L T) AND { ~ / L T ) WITH ( ~ T
OF
/~L
/:l:l= ~.
Lo
4
b)
~.
LD
~,
EQUICONCENTRATION R U N 4, GROUP C.
AND
4" RUH
5. OROUP
AND
,l
EQUI-POTENTIAL
, c)EQUI- CONCENTRATION
L,
[OUl.
L-
LIHESI
S.Om
tf
4,
POTENTIAL
~T:
LINES ,
~'T = IOO.0 m
c.
Fig. 5. Anisotropic dispersivity effect. window), as shown in Fig. 1. Over the window, where the fluid leaves the aquifer, the relative density will be varied but less than 0.025 depending on the c on cen tr atio n of the fluid. When the 2D-FED model was applied to predict saltwater intrusion in the Nile Delta aquifer in Egypt (Sherif et al., 1988), it was found th at isochlor 35.0, which represents the seawater, intruded inland to 65.0 km from the sea boundary, measured at the bottom of the aquifer. The main reason for this severe intrusion is that the aquifer depth at the Mediterranean Sea boundary is ~ 680.0 m. A trial was made to increase the depth of the aquifer to 100.0 m, but convergence was not attained in this case because saltwater intruded all the way to the landward boundary. The length of the studied region would have to be extended to simulate these conditions.
354
M.M. S H E R I F E T AL.
L .r- 4O O O m
bl,4oo
"{~ ,io(~o •
NXX-KZZ
, ,o.s ~
K , . J~.o m / e a y ,~
•
"CT-T,
/
OTO o.aa
~
~
o.ee o.5~ o.e~ oAe o.4( o.la o.ea
LT
"~
k
o.|o N
4. 8
a) V A R I A T I O N
b)
q
q
OF
.
(L D
~. q
~
d
/ L T ) WITH
EGUI-
CONCENTRATION
RUN I
GROUP
AND
(b /bo)
EQ~.POTENT(AL
LINESj b=20.Om.
D.
c ) I=OUI . C O N C E N T R A T I O N N U N 4 G R O U P O.
AND E Q U I - P O T E N T I A L
LiNb'~#D • 80.Ore
Fig. 6. Aquifer thickness effect.
Mesh size To study the effect of mesh size on the c o n v e r g e n c e criterion, the region is divided into four subdomains of equal area (100.0 × 40.0 m) as s h o w n in Fig. 1. Four runs are carried out. In each run, the number of elements in o n l y one subdomain is reduced from 320 elements to 200 elements, i.e. bigger elements are used. In the first run, the number of elements in subdomain (1) is reduced to 200 elements. Standard values for the hydraulic parameters are considered, and c o n v e r g e n c e criterion is set equal to 0.001. Convergence was achieved after 5 iterations. It s h o u l d be noted that w h e n the number of elements in all subdomains was kept at 320, c o n v e r g e n c e was achieved after 4 iterations. In the second run, the number of elements in subdomain (2) was reduced. Convergence was achieved after 8 iterations. In the third run, the same procedure was
SENSITIVITY ANLYSISOF 2D-FEDMODEl,
355
applied to subdomain (3). Convergence was achieved after 14 iterations. The same solution was obtained for the above three runs but with different numbers of iterations. When the number of elements was reduced in subdomain (4), some numerical diffusion problems occurred and no convergence was achieved. This numerical problem usually occurs near the shore boundary where cyclic flow exists and the concentration gradient is relatively high requiring fine discretization (a large number of elements to properly represent the steep gradient). CONCLUSIONS
The 2D-FED model is used to study the effect of several hydraulic parameters on the equiconcentration and equipotential lines in coastal leaky aquifers under steady state conditions. It is found that aquifer thickness greatly affects the length of the dispersion zone as well as the flow pattern. The effect of the other parameters is relatively small except for a certain range. Cyclic flow is confirmed at the shore boundary except when the aquifer depth is relatively small and, yet, salt ions find their way into the aquifer to some limited distance under the effect of dispersion. Selection of proper mesh is very important. Generally, smaller elements accelerate the convergence. Near the shore boundary where cyclic flow exists, an intensive grid is required to avoid numerical problems. REFERENCES Anand, S.G. and Pandit, A., 1982. Finite element solution of coupled groundwater flow and transport equations under transient conditions including the effect of the selected time step size. Proc. 4th Int. Conf. Finite Elements in Water Resources, University of Hannover, Germany. Anand, S.G., Pandit, A. and Sill, B.L., 1980. Some considerations in finite element solutions to coupled groundwater flow and convective~tispersion equations. Proc. 3rd Int. Conf. Finite Elements in Water Resources, 19~23 May, University cf Mississippi, Oxford, MS. Bear, J., 1979. Hydraulics of Groundwater, McGraw-Hill, New York, pp. 22~299. Beran, M.J., 1955. Dispersion of soluble matter in slowly moving fluids. Ph.D. Thesis, Harvard University, Cambridge, MA. Farid, M.S., 1985. Management of groundwater system in the Nile Delta. Ph.D. Thesis, Fac. Eng., Cairo University. Henry, H.R., 1960. Saltwater intrusion in coastal aquifers. Int. Assoc. Sci. Hydrol., 52: 478-487. Herzberg, B., 1901. Die Wasserversorgung einiger Nordseebader. J. Gasbeleuchtung und Wasserversorgung, 44: 815-819, 842-844. Huyakorn, P.S., Anderson, P.F., Mercer, J.W. and White, H.O., 1987. Saltwater intrusion in aquifers: development and testing of a three-dimensional finite element model. Water Resour. Res., 23(2): 293-312. Kawatani, T., 1980. Behavior of seawater intrusion in layered coastal aquifers. Proc. 3rd Int. Conf. on Finite Elements in Water Resources, 19-23 May, University of Mississippi, Oxford, MS. Kinzelbach, W., 1986. Groundwater Modeling. Elsevier, New York, 333 pp. Kitagawa, K., 1934. Sur le dispersement et l'~cart moyen de l'~coulement des eaux souterraines: I. Experiences avec un module de laboratoire. Mem. Coll. Sci., Kyoto Imperial University, Series A, 17: 3742. Klotz, D., 1973. U n t e r s u c h u n g e n zur dispersion in porosen Medien. Z. Deutsch. Geol. Ges., 124: 523 533.
356
M.M. SHERIFET AL.
Lenda, A. and Zuber, A., 1970. Tracer dispersion in groundwater experiments. Isotope Hydrol. IAEA-SM-129/37, pp. 129-134. Mualem, Y. and Bear, J., 1974. The shape of the interface in steady flow in a stratified aquifer. Water Resour. Res., 10(7): 1207-1215. Ogata, A. and Banks, R.B., 1961. A solution of the differential equation of longitudinal dispersion in porous media. Prof. Pap. No. 411-A, U.S. Geol. Surv., Washington, DC. Pandit, A. and Anand, S.C., 1984. Groundwater flow and mass transport by finite elements - - A parametric study. Proc. 5th Int. Conf. on Finite Elements in Water Resources, University of Vermont. Pickens, F.J. and Grisak, E.G., 1980. Scale-dependent dispersion in stratified granular aquifers. Water Resour. Res., 17(4): 1191-1211. Pinder, G.F. and Cooper, H.H., 1970. A numerical technique for calculating the transient position of the seawater front. Water Resour. Res., 6(3): 875-882. Sa da Costa, A., 1981. Seawater intrusion in aquifers: the modeling of interface discontinuities. Proc. of Euromech 143, Delft, The Netherlands. Segol, G., Pinder, G.F. and Gray, W.G., 1975. A Galerkin finite element technique for calculating the transient position of saltwater front. Water Resour. Res., 11(2): 343-347. Senio, K., 1951. On the groundwater near sea shore. General Assembly at Brussels, Int. Assoc. Sci. Hydrol., 11: 175-177. Shamir, U. and Dagan, G., 1971. Motion of the seawater interface in coastal aquifers: a numerical solution. Water Resour. Res., 7(3): 644-657. Sherif, M.M., Singh, V.P. and Amer, A.M., 1988. A two-dimensional finite element model for dispersion (2D°FED) in coastal aquifers. J. Hydrol., 103: 11-36. Umari, A., Willis, R. and Liu, P.L.F., 1979. Identification of aquifer dispersivities in two-dimensional transient groundwater contaminant transport: an optimization approach. Water Resour. Res., 15(4): 815-831. Volker, R.E. and Rushton, K.R., 1982. An Assessment of the importance of some parameters for seawater intrusion in aquifers and a comparison of dispersive and sharp-interface modelling approaches. J. Hydrol., 56: 239-250.
343
Elsevier S c i e n c e Publishers B.V., Amsterdam - - Printed in The N e t h e r l a n d s
[2] A S E N S I T I V I T Y A N A L Y S I S OF '2D-FED', A M O D E L FOR S E A W A T E R E N C R O A C H M E N T IN LEAKY C O A S T A L A Q U I F E R S
MOHSEN M. SHERIF ', VIJAY P. SINGH ' and ABDELWAHAB M. AMER 2
1Department of Civil Engineering, Louisiana State University, Baton Rouge, LA 70803-6405 (U.S.A.) 2Kuwait Institute for Scientific Research, 13109 Safat (Kuwait) (Received February 21, 1989; accepted after revision September 11, 1989)
ABSTRACT Sherif, M.M., Singh, V.P. and Amer, A.M., 1990. A sensitivity analysis of '2D-FED', a model for seawater encroachment in leaky coastal aquifers. J. Hydrol., 118: 343-356. In a previous study, a two-dimensional finite element model for dispersion (2D-FED) was developed and verified for some previously published problems. The model was applied to predict the shape of equiconcentration and equipotential lines in coastal aquifers. In this paper, the 2D-FED model is used to study the effect of several parameters on the simulated dispersion zone and flow pattern in leaky coastal aquifers. The aquifer thickness greatly affects the dispersion zone as well as the equipotential lines. Except for certain ranges, the effects of isotropic and anisotropic permeability coefficients and anisotropic dispersivities are relatively small. Cyclic flow is found near the shore boundary in many cases. Selection of a proper mesh size especially near the seaside is found to be very important to avoid instability problems.
INTRODUCTION
The well-knownGhyben-Herzberg(1889-1891) static equilibriumprinciple (Herzberg, 1901) was the first quantitative formulation for the extent of saltwater intrusion into coastal aquifers. Since that time, many studies have been made in this area. Some investigators (Senio, 1951; Henry, 1960; Shamir and Dagan, 1971; Mualem and Bear, 1974; Sa da Costa, 1981) assumed an abrupt interface between freshwater and saltwater. In reality, however, the freshwater and saltwater merge in a dispersion zone in which concentration and hence the fluid density vary. Where this zone is narrow, the sharp interface approach may provide an appropriate simulation, but an extensive zone of dispersion will alter the flow pattern. Therefore, it is prudent to use a simulation mode] that accounts for the hydrodynamic dispersion effects. One of the first experimental studies made on dispersion was by Kitagawa (1934). Beran (1955) distinguished between three cases of flow: (i) when the effects of molecular diffusion and randomness of the flow pattern are negligible; (2) when the effects of molecular diffusion are dominant; (3) when
0022-1694/90/$03.50
© 1990 Elsevier Science Publishers B.V.
344
M.M.SHERIFETAL.
the randomness of the flow pattern is of equal importance to molecular diffusion in the mixing process. Ogata and Banks (1961) presented a complete solution for the unsteady one-dimensional convective~iispersion equation. It was shown that the simplified form of the equation could be used provided that the dispersion coefficient was small and the region near the tracer source was not considered. Pinder and Cooper (1970) investigated the movement of saltwater in homogeneous confined aquifers. It was assumed that a constant freshwater flux enters the aquifer along the vertical boundary at the land side and a constant saltwater head at the sea side was maintained. The method of characteristics was used to solve the transport equation and the alternating direction iterative finite difference procedure was used to solve the groundwater flow equation for two-dimensional problems. Segol et al. (1975) used the Galerkin finite element method to solve the set of non-linear partial differential equations describing the movement of a saltwater front in a coastal confined aquifer. A layered aquifer was modeled either with a functional representation of permeability or with a constant value of permeability over each element. Umari et al. (1979) addressed the problem of identifying unknown aquifer dispersivities in two-dimensional transient groundwater contaminant transport from given observations on the concentration field. This inverse algorithm became the solution of the contaminant transport problem. The finite element method in conjunction with finite differencing was used to discretize the governing differential equations. Anand et al. (1980) indicated that the selection of a proper mesh size is extremely important to enhance convergence. It was found that solution instability could occur owing to high velocities near the sea side or because of the presence of sharp concentration gradients in a region away from it. For certain values of Peclet number and discharge parameters, saltwater intruded far enough inland that an aquifer length of twice its depth could not be considered adequate for boundary condition of freshwater concentration to be imposed on the landside. Several lengths of aquifers using appropriate finite element grids were used to illustrate this effect. Anand and Pandit (1982) studied the transient conditions of saltwater intrusion in coastal aquifers for a wide range of field parameters. The influence of selected time sizes on computational effort was invstigated. The movement of the saltwater front was found to be far more rapid in the beginning than near the steady state. Volker and Rushton (1982) obtained the location of the abrupt interface in homogeneous isotropic aquifers using the boundary integral method. They compared the location of the interface with the distribution of salt in the dispersion zone where the two fluids meet as determined by a finite-difference solution of the flow and convective~ispersion equations. They found that, for a constant dispersion coefficient, the interface becomes more diffused and migrates inland as the fresh water discharge to the sea decreases. Pandit and Anand (1984) did a parametric study on an idealized confined coastal aquifer. Their results indicated that if the selected region size was not long enough, the flow boundary condition at the freshwater face was not met. The depth of the
SENSITIVITY ANLYSIS OF 2D-FED MODEL
345
aquifer influenced the extent of saltwater intrusion into the aquifer. When the value of the vertical permeability coefficient, K=, was changed, the equiconcentration lines did not change significantly. Huyakorn et al. (1987) employed the picard sequential solution algorithm to simulate the saltwater intrusion problem in single and multiple coastal aquifers. Simulation results were compared with several previously published solutions. Convergence difficulties were encountered when dealing with more complex three-dimensional simulations. Sherif et al. (1988) developed a two-dimensional finite element model for dispersion (2D-FED) to simulate the same problem in confined and leaky aquifers under steady state conditions. The Galerkin technique was employed with computational scheme based on the relaxation method. The model was verified on two previously published problems. Cyclic flow was predicted at the shore boundary. This paper applies the 2D-FED model (Sherif et al., 1988) to study the effect of some hydraulic and geometric parameters on the extent of saltwater encroachment into coastal leaky aquifers and on the shape of equipotential lines. The effect of the element size on improving the convergence behavior is also investigated. AQUIFER CHARACTERIZATION
A domain of 400.0 x 40.0m, represented by a uniform mesh with 1280 triangular elements and 697 nodes, is considered. The rationale for simulating such a small region is that a bigger domain may require a non-uniform grid with smaller elements near the shore boundary to avoid numerical problems. For a wide range of the same parameter, it may be necessary to reconstruct a finer grid to enhance numerical convergence. Parameter values are based on the Nile Delta aquifer in Egypt (Farid, 1985). The aquifer porosity, n, is considered to be 0.3. Hydraulic conductivity in the x and z directions (Kxx and Kzz) are equal to K0 = 150.0m day -1. The upper semi-pervious layer (aquitard thickness is 10.0m), vertical hydraulic conductivity for the upper semi-pervious layer, K'z, is equal to 0.05 m day 1. There are no data available about the dispersivity for the Nile Delta aquifer, but Volker and Rushton (1982) indicated that values of longitudinal dispersivity aL that have been reported in field investigations range from 6 to 150m. Lateral dispersivity, aT, has usually been taken to be 10-30% of the longitudinal value (Pickens and Grisak, 1980). In laboratory experiments, longitudinal dispersivities of 0.1 10cm were found for different types of materials (Klotz, 1973). Longitudinal dispersivities in the field are usually much larger than those in laboratory experiments (Lenda and Zuber, 1970). The main reason for this is the influence of small-scale inhomogeneities in the permeability of the aquifer, which means the onset of macrodispersion. Owing to this scale behavior of dispersion, no generally valid values for longitudinal dispersivities can be given (Kinzelbach, 1986). Values between 0.1 and 500m can be found in the literature for the porous aquifer; Kinzelbach took a L 100 m and aT = 10 m in =
346
M.M. SHERIF ET AL.
the case of one permanent pollutant source of constant strength in a flow field parallel to the x-axis. In light of this discussion, longitudinal and transverse dispersivities % and a T are set equal to 100.0 and 10.0 m, respectively. The water table is set to be 6.0 m above sea level at the land side and 4.0 m at the sea side; between these two nodal points, the water table is assumed to be a straight line having a gradient of 0.05. Therefore, the water table is defined for nodal points on the upper boundary. The piezometric head is set to be 5.80 m at the land side and 3.9 m at the sea side; between these two nodal points, the piezometric head, like any node inside the domain, is not known a priori. Therefore, there will be a leakage flux into the aquifer through the upper semi-pervious layer (aquitard) all over this boundary. The aquitard storage is negligible and the leakage flux is varied from one iteration to another, depending on the difference in levels between the free water table which is constant, and aquifer piezometric head which is varied. On the seaward boundary, where the aquifer is exposed to the open sea, a mixed boundary is used, i.e. on the segment where the fluid leaves the aquifer (window) the concentration gradient is zero and on the segment where the seawater enters into the aquifer, the concentration is equal to seawater concentration C,. The exact length through which saltwater and mixed water enter and leave the domain is calculated by checking velocity directions, as shown in Fig. 1. This boundary is subject to hydrostatic pressure distribution, 7h. Over the window, the specific weight 7 is a function of the mixed water concentration that leaves the system. Below the window, the specific weight is equal to that of seawater 7,. The numerical convergence criterion is set equal to 0.001. Flow and boundary conditions are given in Fig. 1. SEA SIDE
LAND SIDE ~/// ~.,_ ~ t/ ,///, , i// ~,./ - - / ~ / / / / ////////,.//,//
"4./O/I~//,/'.PIZOMETRIC " '.."7/ , I / f l .
4..V ~ T ~c = 0 ~--~
xi.'v
"If
,q~..,~..
.~
WINDOW
j_ (4)
c = cs e =~h
.~
--
~l
I
I I I I I
.l~.'i'..,li//.i
wl Wl
SUBDOMAIN
"/'///////1"/
J
J~
lO0-O m
~1--]-
,~
[
I
(Z)
AQUIFER
I
".1"/
I00.0
m
V¢ l ~ II
~c = 0 ~n I~
•
i.i'ili.'ll''i~J~
[I l I l I I
(3)
l~//'/'//'.,',//., ,a~-/--//,
AQU|TARtD / / Q WATER TABLE :/d 6.0m
)
--'~
-"~
~
~///
4 .A.~iiy
1
~k..,.
SEA
iv I ~1
HEAD
.A.'I
~iS~l/////rf////~.l ,. / / / / / / . / /
"11/i///////////////
q.n= 0 ~:
I00.0 m
I I I I I I
I
~
"¥
4~ ql-- FRESHWATER ~
: ~
"///////
IMPERMEABLE _-It_-
400.0 m
Fig. 1. Boundary conditionsand subdomains in the studied region.
c = cF
~oo.o
m
_I, m]
P = ~f '
SENSITIVITY ANLYSIS OF 2D-FED MODEL
347
GOVERNING EQUATIONS AND SOLUTION TECHNIQUE
The governing equations describing the saltwater intrusion phenomena and the flow pattern in coastal aquifers in steady state conditions are as follows. (1) Darcy equation for groundwater flow: q =
- K(V~ + prVz)
(1)
where K is the coefficient of permeability (LT-1), ~h is the freshwater head (L), Pr is the relative density and is equal to p/pf - 1 (dimensionless), with p being the fluid density (ML 3), and pf being freshwater density (ML 3) and equal to 1000 kg m- 3. (2) The mass balance equation for the fluid: V.pq
= 0
(2)
(3) The mass balance equation for the dissolved salt or the hydrodynamic dispersion equation: nV.(Dh.VC)
-
V'qC
=
0
(3)
where n is the porosity (dimensionless), C is the concentration (ML-3), and D h is the hydrodynamic dispersion coefficient (L2T 1) which can be expressed as follows (Bear, 1979): Dxx =
V:
aL ~
V:
Dxz =
Dzx
=
V~2
D*
V:
D*
+ aT-~-~ +
(aL--aT)--
(3.1)
vxE
where aL and aT are the longitudinal and transverse dispersivities (L), Vx and Vz are the pore velocities in the x and z directions (LT-1), 1~ is the resultant pore velocity, and D* is the molecular diffusion coefficient (L2T - ' ) which is considered to be very small and can be neglected. (4) The constitutive equation relating fluid density to salt concentration is: p
=
pf + a ( C -
Cf)
(4)
where a is a known constant (dimensionless) and equal to 0.7246 for seawater having a density, Ps, of 1025.0 kg m- 3 and a concentration, C,, of 35.0 kg m 3. Cf is the freshwater concentration (ML 3) and is equal to 0.5kgm 3. The above four equations are combined together to form two non-linear partial differential equations in two unknowns, namely, the concentration and the equivalent freshwater piezometric head. The Galerkin finite element technique is adopted to solve these two non-linear equations with the appropriate boundary conditions. Triangular elements are used with linear inter-
348
M.M. SHERIF ET AL.
polation functions. The seaward boundary is updated in each iteration by checking velocity directions at this boundary to locate the exact length through which saltwater and mixed water enter and leave the domain. The computational scheme is based on local linearization of nonlinear terms. A complete discussion of the governing equations, boundary conditions and solution technique is given by Sherif et al. (1988). SENSITIVITY ANALYSIS The 2D-FED model is used to study the effect of several parameters on the width of the dispersion zone, and on the shape of the equipotential lines in a leaky coastal aquifer. The width of the dispersion zone (the extent to which the influence of hydrodynamic dispersion is expected) is measured by the length Lo defining the distance along the bottom of the aquifer from the sea boundary to isochlor line 1.0. The shape of the equipotential lines is measured by the length L' defined as the distance along the bottom of the aquifer from the sea boundary to equipotential line 4.7, which is found to be the first line that intersects the bottom boundary. Equipotential lines 4.6 and less did not meet this boundary in most of the runs. These lengths, L D and L', are normalized with respect to the total length of the aquifer LT which is 400.0m. The variation of the dimensionless ratios (LD/LT) and L'/LT) are studied with various parameter values. For each group of runs, one parameter is varied, and all others remain at the standard values mentioned above. Figure 2 presents the basic run to which other runs are compared.
Effect of isotropic permeability coefficient There are four different groups of runs (Table 1). In Group A, the coefficient of permeability is defined as: Kxx = Kzz = 50.0, 100.0, 150.0, 200.0, 250.0m day -1. It is found that as the isotropic permeability coefficient increases, the length of the dispersion zone also increases while the appearance of cyclic flow at the sea side decreases. Cyclic flow occurs when the velocity field at the sea boundary is rotational in character. Saltwater traverses away from the sea at the bottom of the aquifer and then circulates back into the sea through the upper segment of the boundary. The length of the dispersion zone and the equipotential lines is most sensitive to isotropic permeability when the permeability coefficient lies between 150.0 and 200.0m day -1, otherwise its effect is relatively small, as shown in Fig. 3(a). Increasing the isotropic permeability from 150.0 to 200.0 m day 1 increased the width of the dispersion zone from 164.0-195.0m (about 19%). Kawatani (1980) concluded that the existence of cyclic flow at the sea side prevents saltwater from intruding deeply into the aquifer. This result is confirmed in this study. In run 1 where Kxx = Kzz = 50.0m day -1, there is a strong cyclic flow at the sea side, as shown
349
SENSITIVITY ANLYSIS OF 2D-FED MODEL
BASIC RUN L :/.O0.0m
b = 40.0 m
n=0.3
K =K = K = 150.Orn day xx
zz
~0 SEA SIDE
LAND SIDE 2O
Q 1,
L 100
{,, 200
a) EQUI-CONCENTRATION
SEA SiDE
300
400
LINES.
t30 tttL 20
LAND 5~DE
t
~
100t
2OO
3OO
zOO
b) E~JI-POT~NTIAL LINES. Fig. 2. Basic run-run 3. G r o u p A.
in Fig. 3(b), the length of the dispersion zone, LD, is 152.0m. On the other hand, in run 5 where K=~ = Kz~ = 250 m day-~, the cyclic flow is not as strong as in run 1 and the saltwater intruded the bottom of the aquifer to a distance of 200.0 m as shown in Fig. 3(c).
Anisotropy of permeability coefficient In Group B, the effect of anisotropy of the permeability coefficient is studied. To that end, K~x is kept constant at 150.0m day-1 and Kz~ is varied as 25.0, 50.0, 100.0, 200.0, 250.0 m day- ~(in nature, K~z is always smaller than K=). The length of the dispersion zone, LD, decreases rapidly as (KJKx=) increases. After a certain ratio of (Kz~/Kx=),the anisotropic permeability coefficient has no effect on the dispersion zone and on the equipotential lines, as shown in Fig. 4(a). Decreasing the vertical hydraulic conductivity, K~, from 150 to 25.0 m day increased the width of the dispersion zone from 164.0 to 276.0 m (about 68%). It should also be noticed that the gradient in run 1, as shown in Fig. 4(b), is almost
350
M.M. SHERIF ET AL
TABLE 1 C o m p u t e r r u n s for s e n s i t i v i t y a n a l y s i s
Group A (Permeability effect)
Group B (Non-isotropic permeability effect)
R u n No.
Kxx (m d 1)
Kz~ (m d 1)
1 2
50.0 100.0
50.0 100.0
3 4 5
Ko 200.0 250.0
Ko 200.0 250.0
K0
25.0 50.0 100.0
1 2 3 4 5
Group C (Non-isotropic dispersivity effect)
1 2 3
~L (m)
aT (m)
b (m)
aLo
~L0
bo
aLo
aT0
b0
200.0 250.0
Ko
Ko
aLo
4 5 Group D (Aquifer t h i c k n e s s effect)
5.0 20.0 50.0 75.0 100.0
1 2
3 4
bo
20.0 48.0
K0
Ko
aLo
aT0
60.0 80.0
K 0 = 1 5 0 . 0 m d -1, aLo = 100.0m, aT0 = 10.0m, a q u i f e r t h i c k n e s s , b 0 = 40.0m.
constant except near the shore boundary. In run 5, the gradient is varied gradually from mild at the land side to steep near the sea side, as shown in Fig.
4(c). Transverse dispersivity In Group C, the longitudinal dispersivity aL is kept constant at 100.0m, while the transverse dispersivity is assigned values 5.0, 20.0, 50.0, 75.0 and 100.0 m. It is found that the shape of isochlor lines and equipotential lines are sensitive to anisotropic dispersivities if the ratio (aW/~L) is less than 0.2. Otherwise, the effect of this ratio is negligible. Decreasing the transverse dispersivity, aT, from 10.0 to 5.0m increased the width of the dispersion zone from 164.0 to 224.0m (about 36%). Figure 5 represents the effect of anisotropic dispersivity.
351
SENSITIVITY ANLYSIS OF 2D-FED MODEL
t x • 4oo.om
,,,4o.om
• .o.8
x e -150.0 m / d l 7
LT
j,
~ L T
t
o .ec
o 4e
~ o - - ~ - - - - - ~ - - - o
(
)
o. am
o .zo
~
0.15
CT
0 )VARIATION
j,
L.
OF
( L O / L T) AND
.;"
( L/ / L T ) W i T H
~
b ) EQUI - C O N C E N T R A T I O N RUN I. GROUP A.
AND
L"
(Kxx
/ K O)
*
EQUI- POTENTIAL
L I N E S , Kxx =50.Ore/DAY
LD
¢)EQUI-CONCENTRATION RUN 5 . GROUP A.
AND
EQUI-
POTENTIAL
LINES,Kxx
:ESaOm/OAY
Fig. 3. Permeability effect. Aquifer thickness in hydraulic contact with saltwater Finally, in Group D, the aquifer thickness, b, is varied between 20.0, 48.0, 60.0, 80.0 m, while the length of the studied region Lw is kept at 400.0 m. The same number of elements and nodes are used for all runs. It is found that the aquifer thickness, b, greatly affects the dispersion zone as well as the equipotential lines. The length of the dispersion zone, LD, rapidly increases with the increase of the aquifer thickness as shown in Fig. 6(a). When the aquifer depth is only 20.0 m, cyclic flow does not exist at the shore boundary and the mixed water finds its way out from the aquifer to the sea through the entire length of the boundary. Although there is no flow from the saltwater into the aquifer, salt ions still encroach the aquifer up to some limited distance ( ~ 85.0 m). Salt ions moved from higher concentrations (seawater) to lower concentration inside the aquifers under the effect of dispersion. As the depth of the aquifer
352
M.M. SHER1F E T AL.
LT
lrL
-400.Om
I% +'lOOm
,IOO.O
n = 03
KellBOO~/dly
.__T+eL • --,01
m x x • m,
~P
o.~o
+----
•
~"
o .os o.so o~4SoAo
~
o.al o.lc
old. T I
e.+
O.18
AP ~'~
---
~'e--------
e"~
0,2C
- - e t L" )
0111 OK i
o)
•
,
'
VARIATION
'
~
'
OF
'
t
{L D /L
,
.
+
,
T ) AND
i
.
,
,
(L ¢ /L T )
, --
(~Zz)
WITH
( Kzz/K
0 )
Lo b ) EQUI - CONCENTRATION RUN
,- D
t.
GROUP
AND
EQUI-POTENTIAL
L IN E $ , K Z Z = 2 5 . 0
m/
DAY
E.
..+
-+--~"' +'+
c ) EQUIRUN
CONCENTRATION 5-
GROUP
ANO
EQUI-
POTENTIAL
LINESp
KZZ
=250+0
m/DAY
B.
Fig. 4. A n i s o t r o p i c p e r m e a b i l i t y
effect.
increases, cyclic flow becomes stronger at the shore boundary. When the depth of the aquifer is 80.0 m, the length of the dispersion zone is ~ 270.0 m and a very strong cyclic flow exists near the shore, as shown in Fig. 6(c). Equipotential line 4.7 does not intersect with the aquifer bottom. Volker and Rushton (1982) concluded that the relative density, Pr = P/Pf - 1, of seawater affects the pressures at lower levels in the aquifer near the seaward boundary, and hence influences the length of intrusion into the aquifer. This result is confirmed in this study. As the aquifer depth increases, the effect of the relative density becomes more significant. For example, an aquifer depth of 100.0m and a relative density of 0.025 causes an additional 2.5 m pressure head at the bottom of the aquifer compared with an aquifer filled with freshwater. However, it should be noted that the relative density at the shore boundary will be 0.025 only over the segment where seawater enters into the aquifer (below the
353
S E N S I T I V I T Y A N L Y S I S O F 2D-FED M O D E L
LT
• 4oo.om
k,,40.O
- ~ . ,ooo . . . .
m
A .O
5
Kez I5OOm
/clay
"1
• Kzz - . .
i O.OB o iso 0.4B
C~
~30
O,2~
l't------
o ( L.~L ) LT
C~
- - . - O - -- - - - - I P - - - - - - e l c'~ )
OZC
OIc
a
O)
d
~.
d
VARIATION
o
O
d
=
-
-"
(L D / L T) AND { ~ / L T ) WITH ( ~ T
OF
/~L
/:l:l= ~.
Lo
4
b)
~.
LD
~,
EQUICONCENTRATION R U N 4, GROUP C.
AND
4" RUH
5. OROUP
AND
,l
EQUI-POTENTIAL
, c)EQUI- CONCENTRATION
L,
[OUl.
L-
LIHESI
S.Om
tf
4,
POTENTIAL
~T:
LINES ,
~'T = IOO.0 m
c.
Fig. 5. Anisotropic dispersivity effect. window), as shown in Fig. 1. Over the window, where the fluid leaves the aquifer, the relative density will be varied but less than 0.025 depending on the c on cen tr atio n of the fluid. When the 2D-FED model was applied to predict saltwater intrusion in the Nile Delta aquifer in Egypt (Sherif et al., 1988), it was found th at isochlor 35.0, which represents the seawater, intruded inland to 65.0 km from the sea boundary, measured at the bottom of the aquifer. The main reason for this severe intrusion is that the aquifer depth at the Mediterranean Sea boundary is ~ 680.0 m. A trial was made to increase the depth of the aquifer to 100.0 m, but convergence was not attained in this case because saltwater intruded all the way to the landward boundary. The length of the studied region would have to be extended to simulate these conditions.
354
M.M. S H E R I F E T AL.
L .r- 4O O O m
bl,4oo
"{~ ,io(~o •
NXX-KZZ
, ,o.s ~
K , . J~.o m / e a y ,~
•
"CT-T,
/
OTO o.aa
~
~
o.ee o.5~ o.e~ oAe o.4( o.la o.ea
LT
"~
k
o.|o N
4. 8
a) V A R I A T I O N
b)
q
q
OF
.
(L D
~. q
~
d
/ L T ) WITH
EGUI-
CONCENTRATION
RUN I
GROUP
AND
(b /bo)
EQ~.POTENT(AL
LINESj b=20.Om.
D.
c ) I=OUI . C O N C E N T R A T I O N N U N 4 G R O U P O.
AND E Q U I - P O T E N T I A L
LiNb'~#D • 80.Ore
Fig. 6. Aquifer thickness effect.
Mesh size To study the effect of mesh size on the c o n v e r g e n c e criterion, the region is divided into four subdomains of equal area (100.0 × 40.0 m) as s h o w n in Fig. 1. Four runs are carried out. In each run, the number of elements in o n l y one subdomain is reduced from 320 elements to 200 elements, i.e. bigger elements are used. In the first run, the number of elements in subdomain (1) is reduced to 200 elements. Standard values for the hydraulic parameters are considered, and c o n v e r g e n c e criterion is set equal to 0.001. Convergence was achieved after 5 iterations. It s h o u l d be noted that w h e n the number of elements in all subdomains was kept at 320, c o n v e r g e n c e was achieved after 4 iterations. In the second run, the number of elements in subdomain (2) was reduced. Convergence was achieved after 8 iterations. In the third run, the same procedure was
SENSITIVITY ANLYSISOF 2D-FEDMODEl,
355
applied to subdomain (3). Convergence was achieved after 14 iterations. The same solution was obtained for the above three runs but with different numbers of iterations. When the number of elements was reduced in subdomain (4), some numerical diffusion problems occurred and no convergence was achieved. This numerical problem usually occurs near the shore boundary where cyclic flow exists and the concentration gradient is relatively high requiring fine discretization (a large number of elements to properly represent the steep gradient). CONCLUSIONS
The 2D-FED model is used to study the effect of several hydraulic parameters on the equiconcentration and equipotential lines in coastal leaky aquifers under steady state conditions. It is found that aquifer thickness greatly affects the length of the dispersion zone as well as the flow pattern. The effect of the other parameters is relatively small except for a certain range. Cyclic flow is confirmed at the shore boundary except when the aquifer depth is relatively small and, yet, salt ions find their way into the aquifer to some limited distance under the effect of dispersion. Selection of proper mesh is very important. Generally, smaller elements accelerate the convergence. Near the shore boundary where cyclic flow exists, an intensive grid is required to avoid numerical problems. REFERENCES Anand, S.G. and Pandit, A., 1982. Finite element solution of coupled groundwater flow and transport equations under transient conditions including the effect of the selected time step size. Proc. 4th Int. Conf. Finite Elements in Water Resources, University of Hannover, Germany. Anand, S.G., Pandit, A. and Sill, B.L., 1980. Some considerations in finite element solutions to coupled groundwater flow and convective~tispersion equations. Proc. 3rd Int. Conf. Finite Elements in Water Resources, 19~23 May, University cf Mississippi, Oxford, MS. Bear, J., 1979. Hydraulics of Groundwater, McGraw-Hill, New York, pp. 22~299. Beran, M.J., 1955. Dispersion of soluble matter in slowly moving fluids. Ph.D. Thesis, Harvard University, Cambridge, MA. Farid, M.S., 1985. Management of groundwater system in the Nile Delta. Ph.D. Thesis, Fac. Eng., Cairo University. Henry, H.R., 1960. Saltwater intrusion in coastal aquifers. Int. Assoc. Sci. Hydrol., 52: 478-487. Herzberg, B., 1901. Die Wasserversorgung einiger Nordseebader. J. Gasbeleuchtung und Wasserversorgung, 44: 815-819, 842-844. Huyakorn, P.S., Anderson, P.F., Mercer, J.W. and White, H.O., 1987. Saltwater intrusion in aquifers: development and testing of a three-dimensional finite element model. Water Resour. Res., 23(2): 293-312. Kawatani, T., 1980. Behavior of seawater intrusion in layered coastal aquifers. Proc. 3rd Int. Conf. on Finite Elements in Water Resources, 19-23 May, University of Mississippi, Oxford, MS. Kinzelbach, W., 1986. Groundwater Modeling. Elsevier, New York, 333 pp. Kitagawa, K., 1934. Sur le dispersement et l'~cart moyen de l'~coulement des eaux souterraines: I. Experiences avec un module de laboratoire. Mem. Coll. Sci., Kyoto Imperial University, Series A, 17: 3742. Klotz, D., 1973. U n t e r s u c h u n g e n zur dispersion in porosen Medien. Z. Deutsch. Geol. Ges., 124: 523 533.
356
M.M. SHERIFET AL.
Lenda, A. and Zuber, A., 1970. Tracer dispersion in groundwater experiments. Isotope Hydrol. IAEA-SM-129/37, pp. 129-134. Mualem, Y. and Bear, J., 1974. The shape of the interface in steady flow in a stratified aquifer. Water Resour. Res., 10(7): 1207-1215. Ogata, A. and Banks, R.B., 1961. A solution of the differential equation of longitudinal dispersion in porous media. Prof. Pap. No. 411-A, U.S. Geol. Surv., Washington, DC. Pandit, A. and Anand, S.C., 1984. Groundwater flow and mass transport by finite elements - - A parametric study. Proc. 5th Int. Conf. on Finite Elements in Water Resources, University of Vermont. Pickens, F.J. and Grisak, E.G., 1980. Scale-dependent dispersion in stratified granular aquifers. Water Resour. Res., 17(4): 1191-1211. Pinder, G.F. and Cooper, H.H., 1970. A numerical technique for calculating the transient position of the seawater front. Water Resour. Res., 6(3): 875-882. Sa da Costa, A., 1981. Seawater intrusion in aquifers: the modeling of interface discontinuities. Proc. of Euromech 143, Delft, The Netherlands. Segol, G., Pinder, G.F. and Gray, W.G., 1975. A Galerkin finite element technique for calculating the transient position of saltwater front. Water Resour. Res., 11(2): 343-347. Senio, K., 1951. On the groundwater near sea shore. General Assembly at Brussels, Int. Assoc. Sci. Hydrol., 11: 175-177. Shamir, U. and Dagan, G., 1971. Motion of the seawater interface in coastal aquifers: a numerical solution. Water Resour. Res., 7(3): 644-657. Sherif, M.M., Singh, V.P. and Amer, A.M., 1988. A two-dimensional finite element model for dispersion (2D°FED) in coastal aquifers. J. Hydrol., 103: 11-36. Umari, A., Willis, R. and Liu, P.L.F., 1979. Identification of aquifer dispersivities in two-dimensional transient groundwater contaminant transport: an optimization approach. Water Resour. Res., 15(4): 815-831. Volker, R.E. and Rushton, K.R., 1982. An Assessment of the importance of some parameters for seawater intrusion in aquifers and a comparison of dispersive and sharp-interface modelling approaches. J. Hydrol., 56: 239-250.