Journal of Hydrology, 78 (1985) 279--289 Elsevier Science Publishers B.V., Amsterdam -- Printed in The Netherlands
279
[2] E S T I M A T I N G I N T E R F A C E A D V A N C E D U E T O A B R U P T C H A N G E S IN R E P L E N I S H M E N T IN U N C O N F I N E D C O A S T A L A Q U I F E R S
L.T. ISAACS Department of Civil Engineering, University of Queensland, St. Lucia, Qld. 4067 (Australia) (Received September 7, 1984; revised and accepted December 4, 1984)
ABSTRACT Isaacs, L.T., 1985. Estimating interface advance due to abrupt changes in replenishment in unconfined coastal aquifers. J. Hydrol., 78: 279--289. Dimensional analysis is used to establish a functional relationship between the interface toe position and the other parameters for a one-dimensional, unconfined, coastal aquifer subjected to a rapid decrease in replenishment. Non-dimensional curves for the motion of the toe of the salt-water wedge for some cases are presented. The curves were developed from accurate numerical solutions of the equations for motion of the interface and can be used for reliable estimates of the response of interfaces under the conditions assumed for the analyses. An approximate formula for the rate of movement of the toe is derived for the special case when replenishment reduces to zero. INTRODUCTION A l t h o u g h a n u m e r i c a l s o l u t i o n o f the governing e q u a t i o n s for m o v e m e n t o f the salt--fresh interface in u n c o n f i n e d coastal aquifers is possible (Shamir and Dagan, 1 9 7 1 ) , m a n y scientists and engineers w o u l d n o t have i m m e d i a t e access t o the necessary c o m p u t e r p r o g r a m w h e n p r e d i c t i o n s o f interface m o t i o n are required. There is a need for a simple, reliable m e t h o d for estim a t i n g t h e response o f the interface t o a change in c o n d i t i o n s , particularly w h e n t h e change causes the salt-water wedge to m o v e inland. V a p p i c h a and Nagaraja ( 1 9 7 6 ) o b t a i n e d an analytical solution f o r o n e special case b y using simplifying a s s u m p t i o n s , including a quasi-steady a s s u m p t i o n relating t h e interface profile t o c o n d i t i o n s at the origin where the freshwater discharges t o the sea. C o m p a r i s o n s b e t w e e n results o b t a i n e d f r o m quasi-steady analyses a n d results o b t a i n e d f r o m the m e t h o d o f S h a m i r and Dagan indicate t h a t quasi-stead.y models m a y p r o d u c e significant errors in their results (Isaacs, 1983). Because quasi-steady m o d e l s have serious disadvantages, an alternative a p p r o a c h has been developed. I t is based o n generalized, nond i m e n s i o n a l curves p l o t t e d f r o m values c o m p u t e d b y the m e t h o d o f Shamir a n d Dagan. 0022-1694/85/$03.30
© 1985 Elsevier Science Publishers B.V.
280 MATHEMATICAL PROBLEM
Description of problem The mathematical problem for which a solution is sought is that of onedirectional groundwater flow as shown in Fig. 1. It is assumed that the Dupuit approximation applies, and that there is a well-defined interface between the salt- and freshwater zones. The aquifer geometry and properties, the boundary flow qf(L, t) and the recharge R ( x , t), are assumed to be known. The unknowns for which a solution is sought are the water table elevation s (x, t) and the depth-to-interface ~"(x, t).
Governing equations The governing equations (Shamir and Dagan, 1971; Bear, 1972) are for
x <~xt, freshwater zone: f~s ~" n ~ t + ns 3t
~ [ ~s] 3x Kf(~ +s)~xx
= R
(1)
x ~ xt, salt-water zone: =
for
0
(2)
x ~ xt (only one zone, freshwater): O [Kf(D+s) ~s ]
nf~s Ot
Ox
~xx = R
(3)
Boundary conditions The boundary conditions (assumed known) are s(0, t), ~'(0, t) and either qf (L, t) or s(L, t). The simplified conditions adopted at x = 0 are: s(0, t) = 0
and
~(0, t) = 0
(4) i
Y~ sea Level fresh
satf
"~ce~
wafer
i
. ~
~ue
wafer
,
impermeableIs) L
Fig. 1. Vertical s e c t i o n t h r o u g h an u n c o n f i n e d coastal aquifer.
281 and this simplification will n o t i n t r o d u c e a n y significant e r r o r in the results p r o v i d e d D / L ~ 1. A n o t h e r required b o u n d a r y c o n d i t i o n is t h a t s(x, 0) and ~(x, 0) be k n o w n . In this s t u d y , the a s s u m p t i o n has been m a d e t h a t steady-state c o n d i t i o n s prevail for t < 0.
SPECIAL ASSUMPTIONS F u r t h e r assumptions are i n t r o d u c e d t o r e d u c e the wide range o f possible c o n d i t i o n s to a range for which a single, general solution is feasible. In essence, t h e y imply t h a t aquifer p r o p e r t i e s are i n d e p e n d e n t o f the p r o p e r t i e s o f the w a t e r (salt or fresh) in the aquifer, t h a t spatial variations o f these p r o p e r t i e s are negligible and t h a t a single, i n s t a n t a n e o u s e v e n t is responsible for a n y changes. T h e assumptions are: (1) K ( x ) = K ( a l s o implies K ~ = K s) (2) n ~ = n s = n (3) D ( x ) = D (4) R ( x ; t) -- R ( t ) (5)/3 = (~/s -- 3'f)/7 f = 40 (6) S t e a d y - s t a t e c o n d i t i o n s for t < 0 with an i n s t a n t a n e o u s change (e.g., r e d u c t i o n ) in R (t) a n d / o r qL (t) at t = 0. R (t) and qL (t) remain c o n s t a n t for t>0. The instantaneous reduction implied by assumption (6) is not only convenient for the mathematical analyses, but also, and more importantly, represents a limiting case for a rapid reduction in replenishment and generally gives conservative estimates for the rate of response of the salt--fresh interface.
NON-DIMENSIONAL RELATIONSHIPS
Non-dimensional equations If the governing e q u a t i o n s (eqns. 1--3) are r e w r i t t e n in t e r m s of the n o n - d i m e n s i o n a l values s h o w n in Table 1, t h e y have e x a c t l y the same f o r m as the d i m e n s i o n e d equations. T h e r e f o r e , any p r o g r a m w r i t t e n t o solve the d i m e n s i o n e d e q u a t i o n s can be used to c o m p u t e n o n - d i m e n s i o n a l results, p r o v i d e d the d a t a are i n p u t as n o n < i i m e n s i o n a l values. In Table 1, D O and K 0 are characteristic values f o r aquifer d e p t h and h y d r a u l i c c o n d u c t i v i t y . With t h e assumptions a d o p t e d , D ' -- 1 and K ' = 1.
Dimensional analysis T h e p o s i t i o n o f t h e t o e o f the salt-water wedge at a n y t i m e t, xt (t), will
282 TABLE 1 Non-dimensional variables D' (x! = x = ~"(x, t) = R'(x, t) = q'(x, t) =
D(x)/Do x/D o ~(x, t)/Do R(x, t)/Ko q(x, t)/DoKo
L' s'(x, t) t' g ' (x)
= = = =
L/Do s(x, t)/Do t Ko/Do K(x)/Ko
depend on the initial conditions [represented by xt(0)] , on the aquifer parameters, and on the histories of the recharge and the boundary flow, i.e.: x t (t) = fl [t, x t (0), R ( t ) ,
qL (t), L, D, K, n, ~]
(5)
Dimensional analysis yields:
x~(t)
=
f2 [t', X~(0), R '(t), q~ (t), L', n, ~]
(6)
Equations (5) and (6) do not imply that the value o f x t (t) depends only on the instantaneous values of R (t) and qL (t) at time t. The relationships are true only if the histories of R (t) and qL (t) a r e similar between 0 and t. This restriction is satisfied when R (t) and qL (t) a r e constant for all t > 0. The use of x t ( O ) to represent initial conditions in the aquifer is not strictly correct. The initial conditions (assuming steady state for t < 0) depend on R ( t < ~ 0) and qL (t ~ 0), because different combinations of R and qL can produce the same value for x t (0). Equation (6) is in a form consistent with the non-dimensional variables of Table 1. Other forms are possible because any non-dimensional variable can be replaced by a combination of that variable with other non-dimensional variables. Therefore, the functional relationship for the relative extent of salt-water intrusion can be written: xt(t)/L
t
t
= fa [t', x't(0), R ' ( t ) , qL (t), L , n, fl]
(7)
When R (t) = 0 and qL (t) = 0, the corresponding terms are removed from the functional relationship. When R (t) and qL (t) are constant, but not zero, the interface will move towards a new equilibrium position which is reached as t ~ o~. Provided xt(oo) < L , the value for x t (oo) can be used to describe conditions in the aquifer after the change. As with the use of x t ( O ) to represent initial conditions, the use of x t (oo) is not a strictly correct representation of conditions following an abrupt change. However, it allows a reduction in the number of independent variables. The effects of different combinations of R and qL producing the same steady-state values for x t are considered later. A further reduction in the number of independent variables is obtained if the inspectional analysis m e t h o d (Bear, 1972) is applied to the governing equations. The relationships implied by the governing equations require
283 that tKD/nL 2 must be the same for all systems. This means that the three parameters, t', L' and n, in eqn. (7), can be replaced by the single parameter
tKD/nL z . With these simplifications and combinations, the functional relationship reduces to:
xt(t) L
- f4
(xt(O) Xt(°°) tKD ) L ' L , nL 2,13
(8)
and, if/3 is constant (fi = 40 is adopted):
x tO( t))' -xfts((°x°L) L~
L
, eL 2tKD)
(9)
in which xt (0) represents the steady-state conditions for t < 0 and xt (oo) represents the replenishment conditions (constant for t > 0) following the change at t = 0.
NUMERICAL COMPUTATIONS
Steady-sta te values of x t/L Under steady conditions, the toe position depends on R and qL. The relevant equations (written in terms of non-dimensional values defined in Table 1 b u t with the primes omitted) are: (1) for R 4= 0, and qL -¢ 0:
Xt _ L
qL + 1 - RL
(2) f o r R =~0, qL x_jt 1--
L
(1
[(
)2
_ qL + 1
RL
=
~ +1
~2RLZ
]1J2
(10)
0:
f i + 1~
f2RI?]
(11)
(3) for R = 0, qL :/: 0: Xt
L
~+1
--
2~2qL L
(12)
Unsteady values for xt (t)/L Numerical values of x t (t)/L are produced under the following conditions: (1) A program written to solve the dimensioned equations (eqns. 1--3) is used. (2) Values for R and qL prior to and following the change are chosen so
284 TABLE 2 Non-dimensional values used in analyses for Fig. 2 Curve
(1)
L/D (2)
n (3)
xt(O)/L (4)
R'(t < O) (5)
q't(t < O) (6)
la
50
0.10
0.05
Eqn. (11)
0
lb
50 20
0.10 0.15
0.05 0.05
0 0
Eqn. (12) Eqn. (12)
2a
50 50 100 200
0.10 0.20 0.05 0.30
0.10 0.10 0.10 0.10
Eqn. Eqn. Eqn. Eqn.
2b
50 100
0.10 0.15
0.10 0.10
0 0
Eqn. (12) Eqn. (12)
3a
50
0.10
0.20
Eqn. (11)
0
3b
50
0.10
0.20
0
Eqn. (12)
4a
50 500
0.10 0.05
0.30 0.30
Eqn. (11) Eqn. (11)
0 0
4b
50
0.10
0.30
0
Eqn. (12)
5a
50
0.10
0.40
Eqn. (11)
0
5b
50
0.10
0.40
0
Eqn. (12)
(11) (11) (11) (11)
0 0 0 0
that xt (0) and x t (oo) have predetermined values. (In some cases the specification ofxt(oo ) is replaced by R = qL = 0 for t > 0.) (3) Data are input to and o u t p u t from the program as non-dimensional values (see Table 1). (4) Curves of x t / L versus t K D / n L 2 for specified initial and final states are automatically drawn by a plotter which receives scaled values of x t ( t ) / L a n d t K D / n L 2 directly from the program used to compute solutions to the nondimensional equations. TABLE 3 Non-dimensional values used in analyses for Fig. 3 Curve (1)
L/D
n
xt(O)/L
(2)
(3)
(4)
R'(t < 0) and R ' ( t > O) (5)
q'L(t < 0) and q'L(t> O) (6)
la
50 200
0.10 0.20
0.10 0.10
Eqn. (11) Eqn. (11)
0 0
2a
50
0.10
0.20
Eqn, (11)
0
3a
50 200
0.10 0.05
0.30 0.30
Eqn. (11) Eqn. (11)
0 0
285 RESULTS
Computer-generated curves T h e results f r o m a n u m b e r o f analyses are p r e s e n t e d in t h r e e figures (Figs. 2--4) w i t h details o f t h e n o n - d i m e n s i o n a l values i n p u t to t h e p r o g r a m given in T a b l e s 2--4. In t h e tables c o l u m n s 5 a n d 6 are used to s h o w w h e t h e r R ' or q~ is zero, and, if n o n - z e r o , to indicate t h e e q u a t i o n s b y w h i c h r e l e v a n t values are d e t e r m i n e d . r ! In these analyses, e i t h e r qL ---- 0 f o r all t (curves labeled " a " ) or R = 0 for all t (curves labeled " b " ) . It is a s s u m e d t h a t analyses w i t h a c o m b i n a t i o n o f n o n - z e r o qL a n d R ' w o u l d give results b e t w e e n t h e curves p r o d u c e d . This a s s u m p t i o n was t e s t e d b y an analysis w i t h xt(O)/L = 0.10 a n d a f u n c t i o n o f r r ! b o t h qL a n d R ' f o r t ~ 0 w i t h qL 0 a n d R = 0 f o r t ~ 0. T h e resulting curve did lie b e t w e e n curves 2a a n d 2b o f Fig. 2. T h e t h r e e figures a p p l y to t h r e e d i f f e r e n t c o n d i t i o n s f o r t ~ 0. F o r Fig. 2, R ' = 0 a n d qL = 0. F o r Fig. 3, xt(oo)/L = x t ( O ) / L + 0.3, and for Fig. 4, xt(oo)/L = xt(O)/L + 0.2. Within each figure, curves are p r e s e n t e d f o r a range o f xt(O)/L values. ----
Evaluation T h e results s h o w t h a t the f u n c t i o n a l r e l a t i o n s h i p given b y eqn. (9) is essentially correct. T h e r e is a d e p e n d e n c e on t h e s o u r c e o f r e p l e n i s h m e n t , i.e. R or qL or s o m e c o m b i n a t i o n . In aquifers w i t h a l o w e r value of xt(O)/L, the i n t e r f a c e r e s p o n d s slightly m o r e q u i c k l y w h e n t h e source o f replenishm e n t is recharge. In aquifers w i t h higher values o f xt(O)/L, t h e i n t e r f a c e r e s p o n d s m o r e q u i c k l y if t h e s o u r c e o f r e p l e n i s h m e n t is t h e b o u n d a r y flow. H o w e v e r , the d e p e n d e n c e on t h e s o u r c e o f r e p l e n i s h m e n t is relatively insignificant. TABLE 4 Non-dimensional values used in analyses for Fig. 4 Curve (1)
L/D (2)
n (3)
xt(O)/L (4)
R ' ( t ( O ) and R'(t~O) (5)
q ~ ( t ~ 0 ) and q~{t~0) (6)
la
50
0.10
0.10
Eqn. (11)
0
lb
50 200
0.10 0.10
0.10 0.10
0 0
Eqn. (12) Eqn. (12)
2a
50
0.10
0.20
Eqn. (11)
0
2b
50
0.10
0.20
0
Eqn. (12)
3a
50
0.10
0.30
Eqn. (11)
0
3b
50
0.10
0.30
0
Eqn. (12)
286
I/'/J'"' J~ 0 ...... 1__
2
I
4
I
6 fKD/nL2
8
I
10
A
12
Fig. 2. Non-dimensional response curves for toe movement after replenishment reduced to zero. Given x t ( O ) / L , the rate at w h i c h the interface r e s p o n d s increases as r e p l e n i s h m e n t rates for t ~ 0 reduce. In some real cases, r e p l e n i s h m e n t rates could b e c o m e negative as a result o f losses and extractions. However, in this s t u d y , R ' = 0 and q~ = 0, f o r t :> 0, is used to represent the limiting c o n d i t i o n s f o r a s u d d e n r e d u c t i o n in replenishment. An approximate
formula
If R ' = 0 and q [ = 0, f o r t :> 0, is t a k e n as the limiting c o n d i t i o n , and if a straight line is d r a w n t h r o u g h x t ( O ) / L t a n g e n t to the higher of the a or b curve in Fig. 2, the m o v e m e n t p r e d i c t e d by points on the line w o u l d be 0.5-
0,4 --
3a 2a
0,3
0,2 i 0,1 i
i 2
I 4.
I 6 I-KDInL 2
I 8
l 10
I" 12
Fig. 3. Non-dimensional response curves for total toe movement of 0.3L.
287
0'51
3b
0,~.
0
,
O--
2
I
I
2
L,
I
I
6 8 tKD/nL 2
__
I
10
__i
12
Fig. 4. Non-dimensional response curves for total toe movement of 0.2L.
greater than or equal to the movement predicted by an analysis based on eqns. (1--3). If V* = d ( x t / L ) / d ( t K D / n L 2) is obtained from the slope of the line, a conservative estimate of the velocity at which the toe moves, Vt, is given by:
Vt = V * K D / n L
(13)
Values of V*, derived from Fig. 2, are plotted as a function of xt(O)/L in Fig. 5. An empirical formula which gives a reasonable match to the points is: V* =
0.015 xt (O)/L + 0.12
(14)
0,1-+Values from figure2
Equation ~> 0,5--
I
0,1
I
0,2 xt(o)/L
I
0,3
I
O,Z,
Fig. 5. Non-dimensional toe velocity as a function of initial toe position.
288
Therefore, in an aquifer which satisfies the stated assumptions and for which 0.05 ~ x t ( O ) / L ~ 0.40, a generally conservative estimate of the rate at which the toe will move following a sudden reduction in replenishment is given by:
Vt
0.015
=
KD
xt(O)/L + 0.12 nL
(14)
and the toe position at time t following the change by: x t ( t ) -- x t ( O ) + V t t
(15)
CONCLUSION
A functional relationship: xt(t)L
- F [ xt(O)L ' xt(~)L ' nL tgD1] 2
has been proposed for estimating the toe position of the salt--fresh interface in unconfined coastal aquifers subjected to a rapid decrease in replenishment rates. Numerical results have been used to demonstrate the validity of the functional relationship and to define its form for some cases. For the special case when replenishment reduces to zero results indicate that the toe velocity may be estimated using the simple formula: Vt =
0.015
xt(O)/L + 0.12
KD nL
The use of a functional relationship for which the form can be defined by accurate values from the numerical solution has the advantages that general results are produced and the results can be readily applied to a wide range of situations. The approach used here is not restricted to the cases described. With inclusion of appropriate additional non-dimensional parameters, the same methodology could be used to establish the effects, for example, of spatial variations in aquifer properties, of a sloping impermeable layer or of different patterns of recharge and boundary flow. NOTATION
Symbol
Meaning
D = D(x) f g
thickness of aquifer from sea level down to impervious layer (as superscript) refers to fresh zone acceleration due to gravity
h =P----+ y Pg K L
piezometric head hydraulic conductivity length of aquifer
289 tl f ns P q = q (x, t) R = R ( x , t) s = s(x, t) s
t T = T(x, Vt u
t)
x x t = xt(t ) Y
P 7 A7
f = f ( x , t)
effective porosity for movement of the phreatic surface effective porosity for movement of the interface pore water pressure two-dimensional seepage discharge nett discharge into aquifer from above elevation of phreatic surface above sea level (as superscript) refers to salt zone time transmissivity velocity at which toe moves Darcy velocity horizontal distance measured from origin at intersection of sea level and land surface distance to toe of interface vertical distance measured from origin 1+~ (7 s - - 7 f ) / 7 ~ density of water specific weight of water 7 s -- ,},f distance from sea level to the interface
REFERENCES Bear, J., 1972. Dynamics of Fluids in Porous Media. Environmental Science Series, Elsevier, New York, N.Y. Isaacs, L.T., 1983. Quasi,steady models for dynamic salt--fresh interface analyses. Rep. No. CE 47, Dep. Civ. Eng., Univ. of Queensland, St. Lucia, Qld, 20 pp. Shamir, V. and Dagan, G., 1971. Motion of the seawater interface in coastal aquifers: A numerical solution. Water Resour. Res., 7(3): 644--657. Vappicha, V.N. and Nagaraja, S.N., 1976. An approximate solution for the transient interface in a coastal aquifer. J. Hydrol., 31: 161--173.