Water Resources Management 7: 3-37, 1993 © 1993 Kluwer Academic Publishers. Printed in the Netherlands.

3

Review Article:

Mathematical Models for Saltwater Intrusion in Coastal Aquifers A. G. B O B B A

Rivers Research Branch, National Water Research Institute, Burlington, Ontario, Canada, L7R 4,46 (Received: 27 October 1992) Abstract. Flow of freshwater and saltwater intrusion in coastal aquifers has drawn the attention of many investigators. Several laboratory, as well as mathematical models have been developed to study the pattern of flow of groundwater in coastal aquifers. Mathematical models have wider range of application and are the concern of this paper. Due to the complex nature of the problem, each of these mathematical models are based on certain simplifying assumptions and approximations. This paper presents a critical review of various methods of solution which have been proposed. The validity of the results abtained and the limitations of these models are also discussed. Key words. Coastal aquifers, saltwater intrusion, groundwater, mathematical models, analytical and numerical models.

1. Introduction Coastal aquifers that have their end boundaries in contact with saltwater bodies are invaded by saltwater. A saltwater wedge exists in the aquifer because the freshwater will float above saltwater. However, such a wedge will change in size and shape due to a change in the amounts of natural freshwater flow. The extent of this invasion also depends upon climatic conditions, the characteristics of the natural flow within these aquifers, and the manner of groundwater usage. The invasion of saltwater restricts the agricultural, industrial, and domestic use of groundwater from the invaded aquifer. If the groundwater withdrawal is moderate, no problems should arise, although the saltwater may be expanding its territory. Once the groundwater is extensively withdrawn, the quality of the water may deteriorate, dictating expensive remedies unless proper management is considered. Unfortunately, groundwater management has been neglected in coastal areas. Groundwater management is a highly complex field which involves various social, economic, and technical problems. Although much recent work has been done in this respect, no successful results can be expected until and unless water resources workers or managers recognize the urgent need for basic research in this field. Analysis of some of the unsolved technical problems in the field of saltwater intrusion indicates that basic research is required. It is the main purpose of this

4

A . G . BOBBA

report to review mathematical models for simulating saltwater intrusion in coastal aquifers, and identify areas in which research is required. This basic research must be undertaken to aid in the development of both sound management practices and regulatory measures. Serious problems of saltwater invasion exist in many coastal areas of the world, including the Maritimes. There is very little work done on saltwater intrusion in Canada, although potential areas for concern include the Maritimes and Western coast of Canada (Carr, 1969).

2. Mechanisms of Saltwater Intrusion Saline intrusion of aquifers may originate from various sources, including (1) encroachment of sea water in coastal areas, (2) sea water that entered aquifers during deposition or during a high stand of the sea in past geologic time (connate water), (3) salt in salt domes, thin beds, or disseminated in the geologic formations, (4) slightly saline water concentrated by evaporation in tidal lagoon, playas, or other enclosed areas, (5) return flows to streams from irrigated lands, and (6) man's saline wastes. A major factor in determining the movement of freshwater and saltwater is the relative density between the two. The density of water in natural aquifer systems ranges from 0.9982 kg L -1 at 20 °C for pure freshwater to greater than the 1.345 kg L -1 reported for the Salado brine of New Mexico. The density of 'average' surface sea water ranges between 1.022 and 1.028 kg L-l; this value

v lff

Saltwater (Density 1.025 gm/cm ~) - -

Freshwater (Density 1.000g/m 3) ~

40 ff

/

__,L .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

1 a =1.025

- 1.00

- 40.00

Fig. 1. U-Tube concept illustrating the hydrostatic balance between fresh water and salt water (modified from Todd, 1980).

SALTWATER INTRUSION IN COASTAL AQUIFERS

5

is partly dependent on the temperature as well as the solute concentration. Because freshwater is less dense than saltwater, the former will tend to float above the latter. F o r the hydrostatic case shown in Figure 1, the fresh water stands 1 ft higher than the salt water. Writing an equation for the equality of pressure on each side of the U-tube yields

hs-

Pf

Ps - Pf

hf,

(1)

where pf and Ps are the freshwater and saltwater densities, respectively, hf is the freshwater head above the saltwater, and hs is the distance from the saltwater surface to the interface between the fresh and saline waters. For sea water with a density of 1.025 g / c m 3, Equation (1) reduces to hs = 4 0 × h f

(2)

® GROUNDSURFACE .~_._WATERTABLE.. ............ SEALEVEL

.7. SEA. L ~

///~WA~T~ FRESHWATER •

(~

.7.

EXTRACTIONFIELD IN BASIN

/

...................

/

/

7

..............

.......

SEAWATER

2;; "~ WATERTABLE

' ~ H

rI~ER

Fig. 2. Hydrologic conditions in an unconfined coastal aquifer. (A) Not subject to sea water intrusion. (B) Subject to sea water intrusion.

6

A.G. BOBBA

@ Ground surface

~ V

~

Piezometdclevel ..................

~Y_.___.___

_s_eaJ_ev~. . . . . . . . . . . . . . . . . . . . . .

fresh water ~

~

confined aquifeL___~-

® ExtractionField in Basin

V

~ ! e _ v _ e ! _

/ /

~

J___ 1_ __ L _ _ J . . . . . . . . . . . . . . . . . .

~

I

le~,~,~~

se~-wate r - - - - - - - - - - - ~ ~

Fig. 3. Hydraulic conditions in a confined aquifer in continuity with the ocean. (A) Not subject to

sea water intrusion. (B) Subject to sea water intrusion. Two early European investigators (Ghyben (1989) and Herzberg (1901)) noted that the depth of saltwater in coastal aquifers could be expressed by Equation (2), where hf represented the elevation of the water table above mean sea level. Although not a hydrostatic situation, subsequent investigations have shown that Equation (2) is a good first approximation to the depth for nearly horizontal flow conditions. The interface, consisting of the surface between the fresh and saline waters, has a parabolic form, with the saltwater tending to underride the less dense freshwater. Figure 2, shows cross-sections of an unconfined aquifer for two conditions; Figure 2A, represents an equilibrium between the seaward-flowing fresh water and sea water, whereas Figure 2B indicates the intrusion of sea water into the aquifer when extractions reduce the freshwater flow. Similar situations for a confined coastal aquifer are shown in Figure 3. Whatever the flow situation, a dynamic equilibrium tends to develop. It can be shown that the length of the intruded wedge of saltwater varies inversely with the magnitude of the freshwater flow to the sea. Thus, even a reduction of freshwater flow is sufficient to cause intrusion. A diagnostic approach should precede the solution of any problem; it is, therefore,

SALTWATER INTRUSION IN COASTAL AQUIFERS

7

V J ~

. \

7

vn

ho

/

SEA LEVEL

..FRESH WATER ~(" DISCHARGE ~ FACE

/

t

salt water

~'~

Li

a, b, c, & d, e, g: impervious boundaries a 1, b 1, fl: piezometric surface along d e f a 1, b 1, c1: piezometric surface along a b c a b c g e d: Artesian Aquifer

Fig. 4A. Groundwater flow in artesian coastal aquifers (modified from Kashef, 1986).

Fig. 4B. Groundwater flow in artesian aquifer (modified from Kashef, 1986). of great value to define the modes of saltwater invasion. There are always intricate cases that defy any simple analyses, and for this reason it is usually necessary in applied sciences to idealize the natural conditions in order to have a better comprehension of the problem. Experience and judgement are needed for application of the results of the theoretical ideal cases. In the following discussion, various idealized cases and the availability of their solutions are explained; the practical

8

A.G. BOBBA

V

"'

Dxu~ ~ i x

................. Fresh water

>

i___~_.___~____"~2~_____~ ~--I. . . . . . . . . . . . . - - - - I - ~ ° .~ ~ j °°xo

Impervious Boundary

V Sea Level

Li-

Fig. 5. Groundwater flow in unconfined coastal aquifers (modified from Kashef, 1986).

unknown features are emphasized in the hope of stimulating research workers. It should be noted in this respect that some of these solutions cannot satisfy practitioners of average training. On the other hand, some of these solutions may include mathematical approximations or certain assumptions that may result in errors of unknown magnitudes. In order to understand the modes of saltwater intrusion, we should differentiate between two main coastal aquifers, artesian and phreatic. The artesian aquifer (also know as a confined aquifer) is a pervious soil or rock formation fully saturated with water under pressure. When a bore hole is made, the water rises above the highest saturation level in such a formation (Figures 4A and 4B). On the other hand, a phreatic aquifer (also known as water table aquifer or unconfined aquifer) is a formation in which a groundwater table is established. The water level in a bore hole drilled in such a formation will be approximately the same as that of the water table (Figure 5). In most of the common analyses, the capillary fringe zone and its flow are neglected; also, the formation is assumed to consist of homogeneous and isotropic material, in nature are usually stratified, even if they are more or less homogeneous. The permeability in the direction of stratification is several times greater than the permeability in a normal direction. Some investigators solved many groundwater flow problems where the formation is homogeneous and anisotropic. However, there is not as yet a satisfactory means of determining the ratios of the various permeabilities in nature in such formations. Rock aquifers, on the other hand, may consist of uneven patterns of fissures and cracks; yet they may be considered as homogeneous in the analytical solutions. An idealized case of an artesian aquifer in which the water is at rest (equilibrium conditions) is shown diagrammatically in Figures 4A and 4B. From this idealized figure, it can be shown that within a certain region L, the aquifer may be entirely fresh, entirely salty, or partially fresh, depending on the relationship between dp,

SALTWATERINTRUSION IN COASTALAQUIFERS

9

hf, e~, ho, and e, where e = n(l+a-1). The distances dp, ho, and hfare self-explanatory from Figure 4B. a is constant: -

PY

(3)

Ps - Pf " If dp is equal to or less than e, the region L will be invaded by saltwater. If dp is greater than e but less than (e + a 1 No), the region L will be partially invaded by saltwater and a zone of freshwater of height hf will develop. Finally, if dp >_ e + a-1 h0, the entire region will be free from saltwater. For example, if h0 = 200 ft, e = 100 ft, Ps = 1.025, and pf = 1.00, then o~-1 = 0.025, e = 100 × 1.025 = 102.5 ft, a -~ h 0 = 200 × 0.025 = 5ft. Then if dp < 102.5 ft, there will be no fresh water in region L. If dp > 102.50 ft but < 107.50, a fresh water zone will develop. If d p > 107.50 ft, no salt water would exist and region L will store only fresh water. In this idealized example it is clear that the piezometric surface (or pressure surface) should always be above the sea level by at least 7.5 ft (2.5 m) to prevent salt water intrusion. This transition from freshwater to entirely saltwater takes place when the piezometric surface fluctuates within a range of 5 ft (1.6 m). It is also common in an unconfined aquifer for groundwater to flow seaward at a very low velocity. The natural flow qn varies according to the characteristics of the aquifer and the net magnitudes of water recharge. Thus, due to this flow of the groundwater in an unconfined aquifer, the water table curves toward the sea, and an interface between the freshwater and saltwater develops. Actually, the interface is not a sharp line but rather it is a band of decreasing salinity from its bottom to its top. The case indicated in Figure 5 is mainly an undeveloped area, where the groundwater has not been extensively used and is allowed to be wasted to the sea. However, it as an important basic condition which should be studied before determining the effect of various imbalances such as well pumping. The solution of this problem should include the length L i of the intruded wedge, the location and the shape of the interface, the water table curve, and the size of the discharge faces Dxu and D x above and below the sea level, repectively. The magnitude of Dxu is generally very small compared to D x (Dxu ~ 2.5 % of Dx) (Kashef, 1968). For this reason a large portion of the fresh water would flow below the sea water level. In the available solutions, the discharge face D o is considered either vertical or horizontal in order to simplify the mathematical analysis. The problem may be considered practically solved except for the sloping seaward boundaries. Also, the salt diffusion has yet to be determined by simple means. This simple problem indicated in Figure 5, can be highly complex if the following factors are included: (1) nature of salt diffusion at the interface band, (2) natural or artificial change in qn, (3) effect of precipitation, (4) effect of unsaturated zones above the water table, (5) tidal effect, (6) saltwater movement within the saltwater wedge due to the changes in qn or other factors, (7) evaporation effects, (8) soil anisotropy, (9) actual aquifer boundary, (10) barometric pressure changes, and (11)

10

A . G . BOBBA

artificial water withdrawals by wells and artificial recharge. Most of these mentioned factors and others were studied as individual problems by various investigators. However, more work is needed to digest the results for direct practical use by average practitioners. In artesian aquifers (Figure 4A), a saltwater wedge is formed due to the effect of the natural flow qn. The extent o f intrusion depends upon the permeability of the aquifer, the magnitude of qn, the thickness of the aquifer h0, and the constant a. It should be noted that the geometry of the intruded wedge is independent of the height of the sea water level. The shape of the interface, the length Li, and the discharge face can be determined by various available theories. However, in most of the analytic solutions, the discharge face was assumed either vertical or horizontal (Figure 6). Moreover, the shape of the interface can be found by analogy to the free surface of earth dam problems. Basic research is needed to determine the interface and the discharge face for the case when the actual sloping seaward boundary of the aquifer is considered (Figures 4A and 6). It is already established that the interface in this case is analogous to the free surface on an earth dam with a rock-fill toe. However, there are no direct simple practical solutions for this case although it has been studied extensively by various investigators. Various investigators studied effects produced by the fluctuating magnitudes of

V -~ . ,,~r[eslan A ° euifr

~ , '

| .I J I

Vertical r~' ~h~r~ ~IS~aoe ~ !=1~;

I

Aquifer

~ H o r

izo ntal 2g rge

i

®

®

Vertical Earth Dam (or Core)

/ (

Artesian A~esi_an Aqulrer

, i I

/

/,,

;

, I.L. ..~i \

t

'

i

I

Earth

~

V !nclined Discharge Face

Horizontal Filter (Blanket)

%;~.

~ ' ~ ' ~ ....

~

~

Toe

Fig. 6. Interface analogy between coastal aquifer discharge faces and earth dams.

SALTWATERINTRUSION IN COASTALAQUIFERS

11

qn (nonsteady o r transient flow). Also, extensive research has been done in studying the salt diffusion band along the interface. From a practical point of view, design approaches are still inadequate and more efforts are needed in these areas. The cases just discussed are the simplest possible examples in the field of saltwater invasion. They represent cases of undeveloped coastal regions where the groundwater flow is naturally replenished or depleted. Artificial water withdrawal or recharge by water wells or other means would further complicate these problems. To the writer's best knowledge, the effect of a simple condition of one water well pumping freshwater from a coastal aquifer where the natural flow is known has yet to be determined. In order to understand the complexities of such problems, a water well (Figure 7) is considered to be pumping fresh water from an aquifer of still freshwater with a horizontal interface (as in the case of the aquifer previously discussed and illustrated in Figure 4B). In this case an upward saltwater cone, a phenomenon known scientifically as 'upconing', would develop. Due to the nature of the pressure distribution during moderate pumping, the saltwater would be in a practical state of stagnation, and withdrawal would be prevented by a column of freshwater within the well casing. A small amount of saltwater may be withdrawn from the well when the saltwater seeps through the discharge face. Most likely, it would be mixed with a relatively large amount of freshwater, and it may be unnoticed in the pumped water. However, if the rate of pumping increases, the column of freshwater 'within the well casing would decrease until ultimately the freshwater would be barred from entering the well, and pure saltwater would be withdrawn. In the case shown in Figure 7, the well is assumed to be screened freshwater ' ~

/

limiting case

salt water

I

/

o~~Dw

2i IO~"1hw

r

/

-1

\N

sea level PPS ~ f " Pf

I~i.~I

E =n(l+Of."1)

~\\\\\\\\\\\ hf

~

--©

q

.~ \ Dw

~

,\\\\

/,

Iho \

\\\\\\\\\\\\

\

salt water \\\\\\\\\\\\

Fig. 7. Salt water-fresh water relationships caused by pumping.

\\

12

A.G. BOBBA /

VariablePumping

\.\~\\\\\\\\~\~'N\\\\\\\"

\

).

A,te~la.

V

4/"

Aquifer

J k --2

-- ~ 0

/ 3 % / 1 | ~ -- 1~~'~ "" ~

~-~x,

X

Sea Level

\

Fig. 8. Interface movement caused by pumping a coastal aquifer.

throughout the entire depth of the aquifer. If an artesian aquifer partially penetrating well is in operation in conjunction with the natural flow (Figure 8), more complicated conditions would develop. Both cases given in Figures 7 and 8 have yet to be solved. Solution of these problems must precede the more complicated yet practical problems of multiple well systems. With increased rate of pumping various stages of the interface are likely to occur in the case shown in Figure 8. It should be expected that these interfaces would develop in nonsymmetrical mound forms (Bear and Dagan, 1968). In unconfined coastal aquifers, two approximate cones are likely

GroundSurface

Fig. 9. Draw down curves and upconing caused by pumping.

13

SALTWATER INTRUSION IN COASTAL AQUIFERS

V SEA LEVEL

F R E S H WATER

~- / t

_

~

~

"SALT

W!

Fig. 10. Moderate pumping of a well located beyond the salt water wedge.

to develop around the wells because of the effects of natural flow and well withdrawal (Figure 9). Many workers in the field of saltwater intrusion think that pumping water of the aquifer and beyond the domain of the natural intruded wedge would yield freshwater without any danger of salt invasion. It can easily be shown that, if a well is located beyond the natural intruded saltwater wedge and pumped moderately, there will be no change in the intruded salt wedge (Figure 10). However, if the pumping rate is increased (Figure 11), saltwater upconing is likely to occur (Schmorak and Mercado, 1969). Thus, it is not uncommon to withdraw brackish

........

_<._<.--;

....

;-...-...

Fig. 11. Heavy pumping of a well located beyond the salt water wedge.

14

A.G. BOBBA

or saltwater when the well is located inland beyond the natural saltwater wedge. On the other hand, if the pumping rate is controlled, a well installed within the natural saltwater wedge would yield freshwater only. There are no available satisfactory solutions for the above simple cases. The problem is more complex in a multiple well system. The mechanics of water flow toward artesian wells in coastal aquifers is analogous to gravity wells and can not be the same as that in freshwater aquifers where the superposiition of effects of individual wells and natural flow are exercised. Saltwater invasion has attracted the attention of research workers to analyze new problems in the field of hydraulics of freshwater wells (Kashef, 1986). Whether the practitioners like it or not, there is no hope for satisfactorily solving the practical problems of salt invasion without tremendous efforts in basic research.

3. The Transition Z o n e

The phenomenon of freshwater underlain by saltwater is quite a complex one. In reality the two fluids are miscible and are separated by a transition zone (Figure 12) with a continuous upward gradient of salt concentration from saltwater below to uncontaminated above. This results in a slight upward density gradient in the transition zone. The two fluids tend to mix by molecular diffusion and macroscopie dispersion (due to bulk flow). Because saltwater intrusion represents immiscible

j

GROUND SURFACE

S WATER TABLE V

OCEAN ~ _ ~

~

.~

.q

Fig. 12. Schematicvertical cross-sectionshowing fresh water and salt water.

..

SALTWATER INTRUSION 1N COASTAL AQUIFERS

15

displacement of liquids in porous media, the processes of hydrodynamic dispersion and molecular diffusion tend to mix the two fluids. As a result, the idealized interfacial surface becomes a transition zone (Figure 12). The thickness of the zone depends upon the physical structure of the aquifer, the history of extractions from the aquifer, the variability of recharge, and the tides. Steady flows minimize the transition zone thickness, but variable influences such as pumping, recharge, and tides increase the thickness. Measured thickness of transition zones range from a few meters (feet) in undeveloped sandy aquifers to hundreds of meters (feet) in over-pumped basalt aquifers. Flow within the transition varies from that of the freshwater body at the upper surface to near zero at the lower surface. The velocity in the transition zone transports salt to the sea. Continuity considerations suggest that the salt discharge must come from the underlying saltwater dispersing upward into the zone. It follows that there must be landward saltwater flow, as shown in Figure 12. The transition zone becomes wider and more diffused as the interface moves in response to flow in one or both of the fluids. The two fluids tend to mix by molecular diffusion and macroscopic dispersion (due to bulk flow). A rigorous study of the flow through porous media of freshwater underlain by saltwater, therefore, calls for consideration of hydrodynamic dispersion (Bear and Dagan, 1968, Kashef, 1968 and Reddel and Sunada, 1970) and simultaneous solution of equations of hydrodynamic dispersion and bulk flow with appropriate boundary conditions. However, these equations constitute a complex system of nonlinear equations. Fortunately, in many cases the intermediate transition zone between the uncontaminated freshwater and undiluted saltwater below is small compared to the total freshwater thickness. In such cases, from an engineering point of view, the transition zone can be regarded as an interface separating two fluids of different densities. Most of the work to date concerning flow of freshwater in aquifers wherein it is underlain by saltwater has followed this approach. Dagan and Bear (1968) have taken the average position of the transition zone as an abrupt interface. However, it is important to note that this interface is not an oil-water interface and that unlike in the case immiscrible multiphase flow there is no pressure discontinuity across the freshwater and saltwater interface. To understand and solve the complicated practical saltwater intrusion problems, mathematical modelling is very important before there can be exploitation of groundwater along the coastal areas and small islands.

4. Mathematical Modelling The mathematical models so far developed have been broadly classified in this report as: (a) Piston flow models, (b) Miscible flow models. All the models discussed in this paper refer to the conditions that the water table or piezometric surface lies above sea level and that it slopes towards the ocean. Without these conditions, sea water will advance directly inland.

16

A.G. BOBBA

PISTON FLOWMODELS The piston flow models refer to all those theoretical solutions which employ the concept of sharp interface separating freshwater and saltwater with different densities. The displacement of one fluid by the other can thus be visualised as a piston flow. The mathematical formulation of the problem becomes simpler with the approach and it is possible to obtain approximate solutions for several flow sitiations and boundary conditions. In all the models in this section, unless otherwise stated, the following assumptions have generally been made in order to further simplify the problem mathematically: The aquifer is isotropic and homogeneous, the flow is two dimensional, incompressible and is confined to the freshwater zone, and the surface tension effects can be neglected.

STEADY-STATESOLUTIONS One of the first theoretical analyses of the freshwater and saltwater interface was made by Ghyben (1889) and Herzberg (1901). Each of them assuming that a hydrostatic equilibrium exists between freshwater and saltwater, and that the interface, sea level, and the water table meet at one point, obtained a simple relationship that along any vertical line the depth to the interface hs below the sea level was a constant times the elevation of water table above the sea and is given by Equation (1). The hydrostatic equilibrium assumed by Ghyben and Herzberg implies no flow, yet groundwater flow invariably takes place near the coast. If static conditions existed, the water table would have zero slope and the interface would become horizontal with freshwater overlying saltwater due to the density difference. The Ghyben-Herzberg analysis particularly fails at points close to where the aquifer discharges into the sea, because the former implies zero depth to freshwater at the junction of the aquifer with the sea which is mathematically wrong and physically impossible. The junction, in fact, must be a line rather than a point which can be regarded as an equipotential surface (in three dimensions). The freshwater streamlines must approach the sea perpendicular to this junction. Furthermore, because of the converting streamlines, the flux must increase toward the sea. The equipotential surfaces in the freshwater zone will be curved, their curvature increasing toward the discharge end. Hubbert (1940) considered the dynamic rather than hydrostatic equilibrium of the freshwater and saltwater interface and obtained a relation for the location of the interface which mathematically is identical to Equation (1). However, Hubbert's equation gives the ratio between the elevation of water table above sea level in a particular freshwater equipotential surface and the depth of the interface in the same equipotential surface, whereas in Ghyben-Herzberg relation hf and hs are measured in the same vertical line. For low potential gradients the difference may

17

SALTWATER INTRUSION IN COASTAL AQUIFERS

/ .....

i

tO

f

j

j

I1)

h

-I" i-ll) .O e-

rQ~ I1) 121

J s

~A

",,, ~ , , ~

J

..........................

Salt water

/

................................

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

"~ ~ ~ ~ " X Interlace Fig. 13. Interface between salt water and fresh water•

be negligible but for large gradients the static relationship of Equation (1) may not even approximately be correct. A more realistic picture based on hydrodynamic consideration is given in Figure 13. Using the basic concepts described in this section, some approximate but useful relations have been derived for this shape and location of the interface (Bear, 1979).

EXACT SHAPE OF THE INTERFACE

Considering a general case with simultaneous flow in freshwater as well saltwater equation of the interface can be expressed as (Sahni, 1973) 1

7 = psO~[psChs Pf@]

(4)

where T is the elevation of a point on the interface above some arbitrary datum and 4)s and qSf are the velocity potentials at that point in saltwater and freshwater zones, respectively. Also, the slope of the interface can be given by Sin 6

z

1

k~p-p~[P/qf- Psqs],

(5)

where qs and qf are the volume flux along the interface at the point in question in saltwater and freshwater zones. Thus, in order to determine the exact shape and location of the interface, one must know the actual potential distribution along the interface in both freshwater and saltwater. If, however, the saltwater be treated as static, the problem is simplified. For obtaining the potential distribution, one

18

A.G. BOBBA

must solve the flow equations with suitable boundary conditions pertaining to a particular problem. However, one of the boundary conditions is the interfaces the location and shape of which is not known a priori. Different schemes for solution of such a problem have, therefore, been proposed which are described below. Glovers Modek Glover (1959) developed an equation, based on dynamic considerations, for determining the shape and position of the interface and pattern of freshwater flow near a beach with a known rate of discharge of freshwater under one set of boundary conditions. He adapted, in his analysis, a solution obtained by Kozeny (1953), and by Muskat (1946) with the aid of complex potential for an analogous problem of free surface gravity flow through a semi-infinite underdrained earth dam resting on a horizontal bedrock. The analogy, however, is not perfect as there is no exact correspondence between the boundary conditions between the two problems. Glover in his analysis also assumed the flow through the seepage surface above the sea level to be negligible. He took the water table to be perfectly horizontal standing at the sea level. It was then possible to match the boundary conditions for the freshwater and saltwater interface and the free surface in the dam. Glover's expression for the interface can be written as y2

2Qx

Q2 - -

aK

c~2K

= 0

(6)

and the width of the gap through which the freshwater escapes to the sea as Q x° -

(7a)

2aK '

where x represents a distance measured horizontally landward from the shore line, y is a distance measured vertically downward from sea level, Q represents the freshwater flow per unit length of shore line, and K is the permeability of strata carrying the freshwater flow. It is interesting to note that Glover's model implies assumption of a horizontal water table and yet be obtains the following expression for the height of the water table at a distance x back from the shore line. t

~

hf = ~/(~Q-~] V\

l~

!

(7b)

Charmonman's Solution: An exact analytic solution to the problem of freshwater flow to the sea through a horizontal seepage face in an unconfined coastal aquifer was obtained by Charmonman (1965) by analyzing the governing equation and boundary conditions on the complex potential plane. He obtains the following equations for the free surface and interface respectively, each representing a parabola. The geometry of the upper free surface:

SALTWATERINTRUSION IN COASTALAQUIFERS x' - (a + 1)y'2 2a 2

19

(8)

and, the geometry of the lower free surface: x'

(a + 1)Y '2 2

( a - 1) 2

(9)

where Y'-

Ka Y

Q

and

X'-

KaX

Q

'

where X' and Y' are dimensionless X and Y, K is the coefficient of permeability of the aquifer, X and Y are the coordinates of the physical plane, and Q is the net freshwater discharge per unit width of aquifer. The above analysis gives an exact flow net for the problem and does not assume a horizontal water table unlike Glover's model (1959). However, comparisons of the results of exact solutions of the flow in an unconfined coastal aquifer and the approximate solution obtained by assuming it to be a confined aquifer show (Charmonman, 1965) that the approximate solution is satisfactory for practical purposes. Models using Hodograph Techniques." The basic principle underlying these models

can be outlined as follows. The governing equation for two-dimensional steady state flow through porous media is the Laplace equation in terms of velocity potential 4, V24, = 0 ,

(10)

where 4, is related to specific discharge q by q = -V4,.

(11)

In terms of complex variable z = x+iy and complex velocity potential f(z) = 4, + i ~ the complex velocity is w -

df dz

u-iv,

(12)

where u and v are the specific discharges in x and y direction given by u -

o4, Ox '

o4,

v - --

Oy

(13)

The transformation of the region of flow from the physical plane (z) onto the W or W (conjugate of W) plane gives the velocity hodograph. The utility of the velocity hodograph stems from the fact that, although the shape of the boundaries such as free surface or interface and the limit of the surface of seepage are not known initially in the plane, in the W or W plane their hodographs are completely defined.

20

A.G. BOBBA

It can be shown that the boundary conditions in the flow region on a free surface and an interface are represented in the hodograph W plane by U 2 q- 1,'2 q- K v

= 0

(14)

u2 +

= 0,

(15)

and v2 -

K'v

respectively, where, K'

=

(16)

aK

Both Equation (14) and (15) represent circles in hodograph plane. The step wise procedure for solving the flow problem by hodograph technique can be outlined as (a) the velocities on the boundaries of the physical flow problem (z plane) are mapped on the hodograph plane (W or W), (b) the complex potential for the boundaries of the flow region is determined; (c) a function is found where by the hodograph diagram can be transformed into thefdiagram; and (d) the above function is substituted in Equation (12) and the latter integrated to obtain the solution• Henry (1959) employed the above hodograph technique and conformal mapping to derive equations for the shape and location of the interface in confined aquifers for several sets of boundary conditions. Configurations both with a vertical seepage face and a horizontal seepage face were considered. One of the results obtained by Henry (1959) is of special interest in the present discussion. Consider the flow situation shown in Figure 14 and the corresponding hodograph in Figure 15. For the case of a horizontal strip of aquifer, semi-infinite in length, the points E and F (Figure 14) are at an infinite distance to the left, the two corresponding points in the hodograph are coincident• If in addition to L becoming infinite, the vertical thickness S becomes infinite also, the points D,

V

Salt water

~-i!iiiiiiiiiii!ii!iiiiiiiiiilV/~

A

/ B

/ •

Fresh water

/

/ / •

S

•

°

•

•

/

,

E

7

.~,/ •

>..

/

/

~

~y

D Z Plane, Z = X + iY

Fig. 14. Horizontal confined aquifer with horizontal outflow face (after Henry, 1959)•

SALTWATER INTRUSION IN COASTAL AQUIFERS

F

E

21

D 1"0 tAB B

- Plane W=u+iv

W1 " Plane A

Fig. 15. Hodograph for a horizontal confined aquifer with horizontal outflow face (after Henry, 1959).

E, and F will coincide in the hodograph. Then upon division by Q, the f plane becomes identical to the WI plane (W 1 = -o~K/W), yielding f_ Q

K c~ W

(17)

For this case, Equation (32) of Henry (1959) Z --

- - ~ + constant

(18)

can be integrated to give Z-

f2 2o~KQ "

(19)

This is mathematically identical with the seepage problem solved by Kozeny (Muskat, 1946) for a free surface flow in a semi-infinite aquifer draining downward through a horizontal outflow face with substituted for K. This correspondence between the Kozeny solution and the saltwater intrusion problem for the semiinfinite confined aquifer has been noted and utilized by Glover (1959) for studying freshwater flow in a water table coastal aquifer. Henry (1959) did not, however, carry his solution to the end to present the actual shapes of the interface in a confined aquifer extending to infinite length from the coast. An extension of Henry's analysis for this flow situation was provided by Bear and Dagan (1964). They also showed that under certain conditions, the approximate solution based on the Dupuit assumption is sufficiently accurate. In another study, Henry (1964) obtained an approximate solution for steady flow in a water table aquifer in a long oceanic island using the Dupuit assumption and also a semi-analytic solution with the aid of the hodograph technique and numerical integrations. Like Glover (1959), he neglected all flow above the ocean level and assumed a horizontal discharge face. His results indicate for many purposes the Dupuit assumption is sufficiently accurate in predicting the interface and that when the recharge parameter R / K (R is the rate of uniform vertical recharge per unit area) was as small as 0.0062, the interface position calculated by the two methods

22

A.G. BOBBA

coincide. The Dupuit assumption causes appreciable error only in the situation where the island is wide and the subject of interest is the position of the interface near the shore.

Stratified Aquifer Models: Many coastal aquifers are subdivided by thin impervious or semi-impervious layers (SPL) into a sequence of sub-aquifers. Field observations and laboratory experiments have shown that under these conditions a discontinuity in the shape of the interface occurs such that a fresh water region exists under the semi-impervious layer while immediately above it, saline or mixed water is present in the aquifer. In order to dertermine theoretically the shape of the interface in such coastal aquifers one needs to solve the flow equations simultaneously for all subregions. Recently, some approximate solutions to this problem have been proposed. Ruiner and Shiau (1968) study the shape of the interface in a layered coastal aquifer and showed how the interface, similar to the streamline, is refracted when passing from one layer to the next. Mualem (1973) suggests the special refraction laws that should be applied in such cases. Collins and Gelhar (1971) have probably offered the first analytic description of saltwater intrusion into an aquifer with the SPL, which continuous for a large distance under the sea. They assumed in their analysis that the phreatic surface has a constant slope. The latest solution presented is that by Mualem and Bear (1974). Because of the complexity of the problem, all these models have employed the Dupuit assumption in order to derive an approximate analytic solution. Mualem and Bear (1974) have shown by comparison of their theoretical results with experimental data that the use of the Dupuit assumption does not introduce any significant discrepancy.

Moving Interface Models:

The unsteady movement of the interface is a very complicated problem. The basic equation governing the movement of an interface may be developed (Bear and Dagan, 1964) by assuming the interface to be material surface. These equations are complicated nonlinear partial differential equations for which exact analytic solutions are not possible. Consequently, resort to simplified procedures has been made to obtain some approximate solutions for the moving interface, these are described below. Based on potential theory, Hantush (1968) obtained approximate time dependent analytic solutions for the growth and decay of freshwater lenses and the location of the interface for several flow situations. To make the problem mathematically tractable, the actual three-dimensional flow problem was reduced to an approximate two-dimensional one by assuming that the depth to the interface H is proportional to the head averaged over the vertical line between the water table and the interface. Actually, this depth is proportional to the head at the fresh and saltwater interface. The resulting equation is VZH2 + B =

( SY ) OH2 K-~ Ot

(20) '

SALTWATERINTRUSION IN COASTALAQUIFERS

23

with /3 -

2R K, a(l+(~)

(21)

where Sy is the specific yield of the aquifer and H is the depth to the interface below sea level. This is a nonlinear differential equation and was linearized by replacing H by H, a weighted mean of the areal distribution of the depth of the interface thus making the coefficient in the right-hand side of the above equations a constant. The resulting linearized approximate differential equation is

~72H2 + /3 = ( ~) OH2ot '

(22)

where

_ KaH Sy The mean value H was to be obtained during a time period prior to the time that the location of the interface was desired. Equation (22) was then solved for several different boundary conditions, all of which have both finite and semi-finite areas of recharge. As a consequence of Huntush's first approximation mentioned above, the actual depth to the interface would be less than that estimated by this theory in zones beneath the recharging and the immediately surrounding areas, and is greater than that estimated by the theory in zones beneath the remaining areas. However, the actual depth to the interface during pumping from a fresh water lens, whether by wells or drains, would be generally greater than the estimated value in the regions influenced by pumping, since in these regions the average head along any vertical line in the lens is less than the corresponding head at the interface. Bear and Dagan (1964) studied the case of unsteady flow during the transition period following a sudden change in seaward water of freshwater in a confined aquifer before it reaches its new equilibrium. Such changes may result for example from changes in rates of pumping or recharge in inland areas. The equation for the moving interface even after using the Dupuit assumption was nonlinear. An approximate solution, therefore, was arrived at by employing a static method (Polubarinova-Kochina, 1962). The method proposed by Bear and Dagan requires that the flow rate distribution along the interface be known a priori or asssumed. The accuracy of the solution depends entirely on the accuracy of the assumed distribution. Once such a distribution is assumed, it will give a relationship between flow above the interface, and its shape and position. This relationship corresponds at each instant to a steady interface and flow pattern. Therefore, the period of unsteady flow may be regarded as a series of successive steady states. Each stage may be used as an initial stage for the following interval of time. This in essence is the approach used by Bear

24

A.G. BOBBA

and Dagan to arrive at an approximate solution for the moving interface problem in confined coastal aquifers. Their results could also be applied to water table aquifers, as a further approximation, when the water table has only slight elevations above the sea level. Shima (1969) used more or less the same procedure as outlined above for analyzing the problem of saltwater intrusion when the freshwater aquifer is suddenly exposed to sea water. Approximate analytical solutions such as those above which have been developed to date, apply only to specific and simplified cases. In practical cases, the nature and thickness of the aquifer may vary from point to point. Also, the boundary conditions may not be simple. The natural recharge and the freshwater discharge to the sea and the artifical discharge through pumping, etc. may also vary both in space and time. In such cases even approximate analytical solutions may not be possible and one needs to resort to numerical solutions of the problem. Shamir and Dagan (1971) have presented one such scheme for motion of interface in confined aquifers based on the Dupuit assumption and linearization of equations and treating the motion as two-dimensional. Another approximate numerical solution has been presented by Vappicha and Nagaraja (1976) based on quasi-steady-state analysis. The outflow face in their analysis has been taken as vertical and the aquifer assumed to be standing on a horizontal impermeable bed. The flow boundaries thus are much simpler than those studied by Shamir and Dagan (1971). It has, therefore, been possible to do away with the Dupuit assumption which was necessary in the Shamir and Dagan (1971) model as it would have required a much more complex numercial scheme and more computer time.

MISCIBLE FLOW MODELS This class refers to models which also include the effects of dispersion at the contact of freshwater and saltwater. Problems of such a kind are mathematically more complicated and require special techniques to solve. Owing to the complexity of the problem, only a few approximate solutions for rather simple boundary conditions are available. Where the zone of diffusion is narrow the solutions given in the previous section may provide satisfactory results for the freswater flow pattern, but when the zone of diffusion is rather extensive such as in the Biscayne aquifer of south eastern Florida (Kohout, 1960); 'interface' or immiscible flow models do not yield satisfactory results. In such cases the dispersion effects must be taken into account in the theory employed to describe the nature of saltwater encroachment. Cooper (1959) was the first to advance a hypothesis that under dynamic conditions the saltwater is not static but circulates perpetually from the floor of the sea into the zone of dispersion and back to the sea. This process of returning circulation lessens the extent of sea-water encroachment in the aquifer in comparison with the extent predicted by an interface model according to a mechanism described

SALTWATERINTRUSIONIN COASTALAQUIFERS

25

by Cooper (1959). The diffusion process is enhanced by tidal action and both convection and molecular diffusion are important parts of this mixing process, convection in producing large transfers and molecular diffusion in completing the mixing. Any mathematical model describing the sea-water encroachment problem including dispersion requires the simultaneous solution of the equation of groundwater flow and the equation governing the transport of the dissolved salt. There are two basic distributed-parameter approaches to numerical modelling of aquifers that contain a freshwater lens and a transition zone to saltwater. One approach assumes a sharp interface between fresh and saltwater, so that the system is modelled as two coupled, immiscible fluids. The various areal and cross-sectional versions of sharp-interface approaches are reviewed by Reilley and Goodman (1985). Sharp interface models have either assumed static equilibrium in the saltwater (Fetter, 1972; Anderson, 1976; Ayers and Vacher, 1983; Voss, 1984), or have additionally solved for head or pressure changes in the saltwater (Shamir and Dagan, 1971; Mercer et al., 1980; Wilson and Sa da Costa, 1982; Contractor, 1983). The second modelling approach uses coupled groundwater flow and solute transport to simulate a dispersed transition zone (Pinder and Cooper, 1970; Lee and Cheng, 1974; Segol and Pinder, 1976; Frind, 1982). It is not possible to obtain an axact analytic solution of the above problem, all these mathematical formulations, therefore, make the following simplifying assumptions: An aquifer is homogeneous and isotropic, the fluid motion is two-dimensional, viscosity of the fluids is constant in space and time and although in reality the dispersion coefficient is a function of the average pore velocity of the fluid and is a tensor, it is assumed to be a scalar, constant in space and time. Shamir and Harleman (1967) presented a finite difference method for solution of dispersion problems in a steady, three-dimensional, potential flow field in porous media, in which the miscible fluids had the same density and viscosity. Pinder and Cooper (1970) determined the movement of a saltwater front in confined coastal aquifers, including the effect of dispersion. The method of characteristics was used to solve the solute transport equation. Pinder and Frind (1972) used Galerkin's procedure in conjunction with the finite element technique to simulate seawater intrusion into confined coastal aquifers. Lee and Cheng (1974)formulated a finite element model using a stream function to obtain a steady-state solution for the convective-dispersive transport equation. Segol et al. (1975) and Segol and Pinder (1976) applied the finite element technique to the solution of the transport equation. They used the fluid pressure and concentration as dependent variables. Taylor and Huyakorn (1978) extended their work for a three-dimensional case. The application of this model was limited because of the computational difficulties associated with solution of matrix equations. Sa da Costa and Wilson (1989) also applied the finite element method. Frind (1980) studied the case of a fully confined aquifer overlain by an aquitard that extends out under the sea. He applied the Galerkin finite element technique with linear rectangular elements. The solution was based on a liner interpolation for the concentration between nodal points.

26

A.G. BOBBA

Anand et al. (1980) used Galerkin's finite element technique to solve the flow and transport equations for confined aquifers with appropriate boundary conditions. Several grid sizes were tried using linear triangular elements. It was found that solution instability could occur due to high velocities near the sea water face or due to the presence of sharp concentration gradients in a region away from it at high P6clet numbers. Katwatni (1980) presented a two-dimensional model for the behaviour of sea water intrusion in layered coastal aquifers. He addressed the instability problems which occurred when the convective terms exceeded a certain criterion related to the dispersion coefficients and the size of the finite elements. Katwatni took the longitudinal dispersion coefficient as a constant to avoid the instability and other numerical problems. Wikramaratna and Wood (1982) solved the nonlinear coupled steady-state differential equations governing groundwater flow and movement of pollutant in confined aquifers. The numerical scheme combined the perturbation method and standard Galerkin finite element technique. Pandit and Anand (1984) did a parametric study on an idealized confined aquifer with the same boundary conditions as ultilized by Henry (1964). Their results indicated that the depth of the aquifer influenced the extent of saline water intrusion into the aquifer. When the value of the vertical permeability coefficient changed, the equi-concentration 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 to simulate the same problem in confined and leaky aquifers under steady-state conditions. The Galerkin technique was employed with a computational scheme based on the relaxation method. The model was verified on two-dimensional two previously published problems. Cyclic flow was predicted at the shore boundary. It was also concluded that cyclic flow existed when the characteristic velocity at the seaward boundary was bigger than longitudinal dispersivity. Voss (1984) employed a twodimensional hybrid finite element method and integrated finite difference method to simulate the fluid-density dependent saturated or unsaturated groudwater flow and transport of solute in groundwater. The simulation may be employed for areal and cross-sectional modelling of saturated groundwater flow systems, and for crosssectional modelling of unsaturated zone flow. Huyakorn et al. (1987) developed a three-dimensional finite element model for simulation of saltwater intrusion in single and multiple coastal aquifers with either a confined or phreatic top aquifer. The model employed the Picard sequential solution algorithm with special provisions to enhance convergence of the iterative solution. Reilly and Goodman (1986) analyzed saltwater upconing beneath a pumping well. Although the upconing phenomena are not exactly similar to salwater intrusion, their analysis is instructive. The upconing is in response to the pressure reduction due to drawdown of the water table around the well. If the bottom of the well is close to the saline water level or the well

SALTWATERINTRUSION IN COASTALAQUIFERS

27

discharge is relatively high, the saline water may reach the well. In some cases, the saline water cone reaches the bottom of the well with a sudden jump, indicating conditions of instability. The dispersion phenomenon is not the reason for the upconing phenomenon. On the other hand, the saltwater intrusion may occur without any water extraction from the aquifer, the dispersion phenomenon will take place, and a transition zone will exist. Some of the important numerical models are discussed in the following sections.

Numerical Solutions: The governing equations are: 7.q=0,

(23)

where the seepage velocity q is given by Darcy's law: K q = - - -#( V p - o g )

(24)

D72C = - V.(qC).

(25)

and

Here the symbols carry the same meaning as before. The density p can be expressed in terms of salt concentrations as

p = pf l + o e

,

(26)

where Of and Cs denotes the density of freshwater and the salt concentration of seawater. It can be shown (Lee and Cheng, 1974) that the components of seepage and 7 along Y and ~ can be written as

0'~ -

Kfd @ -

,

(27)

Qs V -

Ogr _ 07c

K f d _0~ + ( l + o ~ C ) ,

(28)

QS

where the bars denote the dimensionless variables and

g

a' _c

p

p

Cs '

pf

pfgd

with Qfas the average discharge of freshwater per unit width and das the characteristic thickness of the aquifer. Eliminating the pressure terms from the above two equations by cross-differentiation, we get the following systems of equations for fluid motion and distribution of salt:

28

A. G. BOBBA

V2~ -

1 06 a 0~

(30)

'

[ Ogr OC Ogr OC ] V2--C= Pe Oy O~C O.~ @

(31)

where the dimensionless parameters 'a' and Pe are the discharge parameter Qf/Kfdx and the seepage P6clet number Qf/D, respectively. This seepage P6clet number plays a role here similar to that of the Reynolds number in viscous flows. The governing Equations (30) and (31) of ~0 and C are of elliptic type. In order to solve them, appropriate boundary conditions for the entire boundary flow region in a particular problem of interest must be prescribed. Solutions presented by Henry (1960) and Lee and Cheng (1974) will now be discussed.

Henry's Model." In

addition to the general assumptions listed earlier for miscible flow formulation, Henry (1960) further simplified the problem by considering a confined rectangular aquifer. Schematic representation of Henry's problem for this idealized aquifer together with boundary conditions in terms of non-dimensional variables above is given in Figure 16. A semi-analytical solution to this problem was presented in the form of a Fourier double series, the coefficients of which were obtained numerically. With a -- 0.263 and Pe =10.0, Henry's results confirmed the existence of saltwater circulation in coastal aquifers induced by dispersion which was first hypothesised

V

CONSTANT SALT WATER HEAD m

X

a~l'

=0

Qf

~=1

~=-1,

:.................................................................................................. Y

Fig. 16. Schematicrepresentationof Henry'sproblem and its boundaryconditions.

SALTWATERINTRUSIONIN COASTALAQUIFERS

29

by Cooper (1959). But this numerical scheme was unable to predict the encroachment for large P6clet numbers such as describing the actual field conditions of Biscayne aquifer (Kohout, 1960) due to numerical difficulties.

Lee and Cheng Model." Lee and Cheng (1974) solved the system of coupled nonlinear conservation equations of mass and of salt, namely Equations (30) and (31), for the steady-state two-dimensional seawater encroachment problem by the finite element method with the aid of iteration and under-relaxation. An inspection of Equations (30) and (31) reveals that they are very similar to the Navier-Stokes equations when the latter are written in terms of stream function and vorticity. The same numerical scheme was used in this problem as was employed by Cheng (1972) for obtaining numerical solutions of the Navier-Stokes equation. The method of solution can be summarized as follows: The nonlinearity of the equations was circumvented by using an iteration procedure so that Equation (30) is considered to be Equation (30) for O and Equation (3l) for C, the right-hand sides of each equation being the forcing terms. These equations are then solved alternatively by the finite element method using triangular elements and with the aid of an under-relaxation imposed on until the convergence is reached. At each iteration the right sides of Equations (30) and (31) are replaced by the solutions of O and C obtained from the preceding iteration. The advantage of the finite element method is that it can be satisfactorily applied to flow regions with boundaries of any shape. The numerical model presented by Lee and Cheng (1974), therefore, has the special advantage over Henry (1960) in that the former can easily be adapted to actual field problems. To verify the validity of their model, Lee and Cheng solved Henry's idealized problem under exactly the same boundary conditions and for the same values of parameters 'a' and Pe by their numerical method using linear triangular elements. The 0.5 isochor (C = 0.5) obtained by this method was compared with that reported by Henry (1960). They were found to be in good agreement. To test the applicability of their model to actual field problems of sea water encroachment, they studied the flow in unconfined Biscayne aquifer. A comparison of the numerical results with the actual field observations showed a good qualitative agreement, but some quantitative differences were observed. This Lee and Cheng suggest, could be due to the uncertainties in the actual values of °a' and Pe and also to the simplifying assumptions used in their model. In particular, the aquifer modeled is probably highly anisotropic nature of the dispersion coefficient may not be negligible. Actually, the dispersion coefficient, as explained earlier, is a secondorder tensor and is not only a property of the porous medium but is also a function of the dynamics of the fluid through this medium. Pinder and Cooper Model." A numerical solution for determining the transient position of the saltwater front and the pattern of flow under the effects of dispersion has been developed by Pinder and Cooper (1970). Under the same general assumptions

30

A.G. BOBBA

as listed earlier in this section on miscible flow models, the equation governing transport of the dissolved salt for the transient case can be written as DVc 2

-

q"

(32)

VC -- OC

0

Ot"

The other governing equation for groundwater flow will be same as Equation (23). The density P is also related to C as expressed by Equation (26). Pinder and Cooper (1970) used the method of characteristics to solve the systems of equations for the transient development of the saltwater front in a rectangular confined aquifer. The accuracy of the numerical technique proposed was verified by comparing the numerical and analytical solutions for two problems treated by Henry. To verify the technique for the transient problem they used Henry's solution for the motion of a sharp interface in a U-tube, a special case of Henry's (1962) solution for the movement of saltwater front in a sloping artesian aquifer. To verify the technique for the two-dimensional flow, since analytical solution for the transient problem with the dispersion effects was not available, they used the steady-state solution given by Henry (1960) for an idealized confined aquifer. They obtained numerical solutions for these two problems with exactly the same boundary conditions. In the first problem their results show excellent agreement with analytic solution. For the second problem, Pinder and Cooper used two different initial conditions, namely, the aquifer was either filled with fresh water or assumed to have a steady interface initially, and calculated the movement of the salt front. They have shown that the unsteady salt front as calculated by their method tends to approach, from both sides, Henry's steady-state solution as time 't' increases (Figure 17). Pinder and Cooper thus verified the validity of their model by demonstrating good agreement reached between their solutions and the analytic solutions for the above two simple 10_

i

-

~S ~

-

5_

.<.;.',.: ..<.-,;.'. ....

--

_.-

-

s~r,~,

_

.-

.....

.......... ....... .

......... ,,,........

t 7

0

,0

I

I

I

I

I

I

I

l

~/~,

PINDER & C O O P E R

// LEE & C H A N G

I

I 10

I

I

I

I

I

I

I

I -1~ 20

Fig. 17. The 0.5 isochors o f Henry's problem c o m p u t e d by various m e t h o d for d = 1,0; a = 0.263, P e - 10.0.

SALTWATERINTRUSION IN COASTALAQUIFERS

31

cases. However, they did not show if their scheme works satisfactorily at higher values of seepage P6clet number. It would also have been interesting to see, when applied to a real field problem, how far the results predicted by their numerical scheme agree with actual field observations.

SUTRA Model." (S_aturated, Unsaturated, TRAnsport) is finite-element groundwater flow and energy (heat) and solute (contaminant) transport code which was developed by U.S. Geological survey (Voss, 1984). The SUTRA model employs a twodimensional, finite-element method to simulate variable density groundwater flow, and transport of the total dissolved solids (TDS) in the groundwater. This model has many of the advantages associated with finite element models. It produces stable solutions with a relatively coarse mesh and introduces small numerical dispersion. The two primary variables upon which the model is based are fluid pressure, and TDS concentration expressed as mass fraction. The solute mass balance per unit aquifer volume at a point (x,z) in an aquifer with variable-density fluid is given by (Voss, 1984): 0C

Op ff~- + Opv'VC-V_'[Op(DmI+D)'VC ] = Qp(C*-C),

(33)

where C(x,z,t) is solute concentration as a mass fraction (mass solute/mass fluid) [Ms/M] where [Ms] are units of solute mass and [M] are units of fluid mass; O(x,z) is aquifer volumetric porosity; v(x,z,t) is fluid velocity [L/T]; D m is molecular diffusivity of solute in pure fluid including aquifer material tortuosity effects [L2/T]; I is the identity tensor; D(x,z,t) is the dispersion tensor [LZ/T]; Qp (x,z,t) is a fluid mass source (mass fluid/aquifer volume/time) [M/L3.T]; C*(x,z,t) is concentration of solute as a mass fraction in the source fluid [Ms/M]; 0 (x,z,t) is fluid density [M/L3], where [L 3] is fluid volume and density is given as a linear of concentration: 0o P = P0 + ~-~(C- Co),

(34)

where P0 is fluid density when C = Co; Co is a base solute concentration; and Op/OC is a constant coefficient of density variability. Darcy's law gives the massaverage fluid velocity as =

-

.(_Vp

-

m),

(35)

where k(x,z) is permeability tensor [L2];/1 is fluid viscosity [M/L'T]; g is the gravity vector [L/T2]; p(x,z,t) is the fluid pressure [M/L-T2]. The mass balance of fluid per unit aquifer volume at a point in the aquifer, assuming a rigid solid matrix may be given by

Op ~O0 OC pSop ff[ + ~ - ~ [ - + V" (Opv) : Qp,

(36)

32

A . G . BOBBA

where Darcy's law, Equation (34), for mass-average velocity may be substituted in the fluid mass balance Equation (36). The specific pressure storivity, Sop, for a rigid solid aquifer matrix [M/(L.T2)] -1 is given by Sop = (1-0)~ + 0/3,

(37)

where o~is porous matrix compressibility [M/(L.T2)] 1, and/3 is fluid compressibility [M/(L-T2)] -1. A general boundary condition for the fluid mass balance that applies at stationary boundaries as well as approximately at a water table (Bear, 1979) is S y 0p

Igl Ot + Opv.n +

Q.

= 0,

(37)

where Sy(x,z) is the water table specific yield (volume fluid released/aquifer volume) for unit drop in hydraulic head; Igl is the magnitude of gravitational acceleration, [L/T2]; and Q* is a fluid mass source due to flow across boundaries (mass fluid recharged per unit area of boundary/time), [M/L2T]. This condition is accurate when at a water table ]Pl = P]gt, that is, when pressure gradients are nearly hydrostatic at the water table. This condition is usually satisfied because commonly vertical pressure gradients are often nearly hydrostatic and horizontal gradients are much smaller than vertical gradients. The mechanical dispersion tensor in two spatial dimensions is given by

D = Dxx Dxz Dzx Dzz '

(39)

Dzz = ( ~2) (dTVx2 + dLVz2) ,

(41)

Dxz = Dzz = ~T (dL - dT) (VxVz) ,

(42)

dL = O~LlVl,

(43)

d T = aT~Vl,

(44)

where Iv[ is the magnitude of velocity and dL and dT are longitudinal and transverse coefficients of mechanical dispersion, and, O~L(X,Z) is longitudinal dispersivity [L] and ~T(X,Z) is transverse dispersivity [L]. The intrinsic structure of finite element models provides substantial flexibility in designing the mesh and conforming to the boundary conditions as well as making the placement of material property boundaries such as hydraulic conductivity boundaries more realistic. Although the code does not allow triangular elements,

SALTWATER INTRUSION IN COASTAL AQUIFERS

33

it compensates for this by providing pinch nodes. Pinch nodes are nodes that end on an element side. The input also does not allow direct variation of the storage coefficient but allows variability of this coefficient through changes in porosity, saturated thickness, and compressibility. This model has been applied to real field data and observed the favourable results (Bush, 1988; Souza and Voss, 1987; Voss and Souza, 1987; Bobba, 1991, 1992).

5. Summary From the foregoing review of the serious present problems in the field of groundwater pollution by saltwater intrusion, it is obvious that many problems remain unsolved and more extensive basic research is required. The following is an attempt to summarize the main problems. Some of these may have been investigated partially or entirely following the usual scientific pattern; some others may need further evaluation to render practical useful applications. (a) Proper management of groundwater in coastal aquifers is very vital as it forms the main source of water supply in these regions. A systematic study of the behaviour of these aquifers under different conditions is essential in scientific design and planning of water resources systems of coastal regions. This requires adequate tools for the prediction of the motion of a saltwater body in the aquifer under different flow conditions. (b) The laboratory models are usually expensive to construct, the tests made with them are time-consuming and the results may be applicable to only certain specific field conditions. Mathematical models both analytical as well as numerical, are therefore, generally preferred as these are less expensive and faster and can be suitably modified to encompass wider range of flow situations. (c) The mathematical models developed so far to solve the problems of freshwater flow in a coastal aquifer have, broadly speaking, used two approaches, namely, the immiscible flow or interface concept and the miscible flow concept. The former have been referred to in this report as Piston flow models. The formulation of the problem and its solution in the Piston flow models is mathemetically much simpler than is the case with miscible flow models which also include the effects of dispersion. (d) Even when the sharp interface concept is employed together with some simplifying assumptions, the system of equations to be solved turns out to be nonlinear one. Inclusion of the dispersion effects makes the problem even more complicated. With the available analytic techniques and computers, it has been possible to solve only a few rather simple and specific flow problems when dispersion is included. The interface models, though less vigorous in treatment, on the other hand, have been able to analyze the problems with several different boundary conditions. (e) Study of saltwater problems in unconfined aquifers considering the actual sloping seaward boundaries is still in the preliminary stages and should receive further attention.

34

A.G. BOBBA

(f) In all the theoretical analyses, it is always safer to use the hydrodynamic equilibrium approach, although in a few flow situations the simpler GhybenHerzberg approach may yield nearly the same results. (g) The hodograph technique is one of the most powerful analytical methods for solving steady state two dimensional interface problems especially in unconfined aquifers, where two boundaries of flow along phreatic surface and the interface are not known a priori. (h) The transient interface problems can probably be best studied by use of quasisteady state approximation and numerical analysis. (i) Often, the actual field problems are such that the flow boundaries can not be simulated by a simple geometrical configuration without introducing appreciable error in theoretical analysis. In such cases, it would be better to use finite element techniques rather than difference techniques. (j) The effects of tides, barometric pressures, rainfall, evaporation, and leakage on the development of saltwater mounds need more attention. (k) No satisfactory model is presently available for simulating and analysing real field problems of freshwater flow including dispersion effects in coastal aquifers. However, the numerical scheme developed by Voss (1984) for the saltwater intrusion problem is reliable to apply but this model needs more computer memory and computer time. More work is needed in developing methods for solving flow problems with dispersion effects, both for steady state as well as transient flow conditions. This model has to be coupled with a stochastic model to estimate hydrogeological parameters. (1) In Canada, there is very little work done on the saltwater intrusion problem (Carr, 1969). It is necessary and important to do basic and applied research due to climate change effects and urban development on coastal areas.

Acknowledgements The author would like to thank Allan Crow and William Botty for comments on the manuscript and John Carey and Suzane Lesage, and Kent Novakowski who encouraged me to prepare this manuscript.

References Anand, S. G., Pandit, A., and Sill, B. L., 1980, Some considerations in finite elements solutions to coupled groundwater flow and convective-dispersion equations, Proc. 3rd Int. Conf. Finite Elements Water Resour., Univ. of Mississippi, Oxford, Miss. Anderson, M. R, 1976, Unsteady groundwater flow beneath strip oceanic islands, Water Resour. Res. 12, 640-644. Ayers,. J. F., and Vacher, H. L., 1983, A numerical model describing unsteady flow in a fresh water lens, Water Resour. Bull. 19, 785-792. Bear, J., 1979, Hydraulics of Groundwater, McGraw-Hill, New York. Bear, J. and Dagan, G., 1964, Some exact solutions of interface problems by means of the hodograph method, J. Geophys. Res. 69, 1563-1572.

SALTWATER INTRUSION IN COASTAL AQUIFERS

35

Bear, J. and Dagan, G., 1968, Solving the problem of local interface upcoming in a coastal aquifer by the method of small perturbations, J. Hydraul. Res. 1. Bear, J., Zaslavsky, D., and Irmay, S., 1989, Physical Principles of Water Percolation and Seepage, UNESCO, Paris. Bobba, A. G., 1991, Application of digital simulation model to freshwater aquifer of Lambton County, Ontario, Canada. Paper presented at International Conference on Hydrology and Hydrogeology in the '90s; Orlando, Florida, U.S.A., 3-7 Nov., 1991. Bobba, A. G., 1992, Field validation of 'SUTRA' groundwater flow model to Lambton County, Ontario, Canada. NWRI contribution 92-141. Bush, W. R, 1988, Simulation of saltwater movement in the Floridan aquifer system, Hilton Head Island, South Caroline., U.S.G.S Water Supply paper 2331. Carr, E A., 1969, Salt-water intrusion in Prince Edward Island, Canadian J. Earth Sci. 6, 63-74. Charmonman, S., 1965, A solution of pattern of fresh water flow in an unconfined coastal aquifer, J. Geophys. Res. 70, 2813-2819. Cheng, R. T., 1972, Numerical solution of the Navier-Stokes equations by the finite element method, Phys. Fluids. 15, 2098-2105. Collins, M. A. and Gelhar, L. W., 1971, Seawater intrusion in layered aquifers, Water Resour. Res. 7, 971-979. Contractor, D. N, 1983, Numerical modeling of saltwater intrusion in the Northern Guam lens, Water Resour. Bull. 19, 745-751. 'Cooper, H. H. Jr., 1959, A hypothesis concerning the dynamic balance of freshwater and saltwater in a coastal aquifer, J. Geophys. Res. 64, 461-467. Dagan, G. and Bear, J. 1968, Solving the problem of local interface upcoming in a coastal aquifer by the method of small perturbations, J. Hydraul. Res. 6, 15-44. Drabbe, J. and Baden Ghyben, W., 1888-1889, Nota in verband met de voorgenomen putboring nabij Amsterdam, Tijdschrift van het Koninklijk Instituut van Ingenieurs, The Hague, Netherlands, 8-22. Fetter Jr. C. W., 1972, Position of the saline water interface beneath oceanic islands, Water Resour. Res. 8, 1307-1314. Erind, E. O., 1980. Seawater intrusion in continuous coastal aquifer-aquitard system., Proc. 3rd Int. Conf. Finite Elements Water Resour., Univ. of Mississippi, Oxford, Miss. Vol. 2, pp. 177-198. Frind, E. O., 1982, Simulaton of long term transient density dependent transport in groundwater, Adv. Water Resour. 5, 73-78. Frind, E. O., 198i, Seawater intrusion in continuous coastal aquifer-aquitard systems, Adv. Water Resour. 5, 89-97. Ghyben, W. B., 1889, Notes on the probable results of well drilling near Amsterdam (in Dutch). Tijdschrift van het KoninkHjk Inst. van Ingenieur , The Hague 21. Glover, R. E., 1959, The pattern of fresh-water flow in a coastal aquifer, J. Geophys. Res. 64, 457459. Hantush, M. S., 1968, Unsteady movement of fresh water in thick unconfined saline aquifers, Bull. Int. Assoc. Sci. Hydrol. 13, 40-60. Henry, H. R., 1959, Saltwater intrusion into freshwater aquifers, J. Geophys. Res. 64, 1911-1919. Henry, H. R., 1960, Saltwater intrusion into coastal aquifer, Int. Ass. Sci. Hydrol. Publ. 52, 478-487. Henry, H. R., 1962, Transitory movements of the saltwater front in an extensive artesian aquifer, U.S.G.S. Prof. Paper. 450-b: 87-88. Henry, H. R., 1964, Interface between salt water and fresh water in coastal aquifers, sea water coastal aquifers. U.S.G.S. Water Supply Paper, 1613 C. Herzberg, B., 1901, Die WasserversQrgung einiger Nordseebader, J. Gasbeleuchtung und Wasserversorgung. 44, 815-819, 842-844. Hubbert, M. K., 1940, The theory of ground water motion, J. Geol. 48, 785-944. Huyakorn, R S., Anderson, E 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, 293312. Kashef, A. I. 1968, Fresh-salt-water interface in coastal groundwater basins, Proc. Internat. Assoc. Hydrogeologists. 8, 369-375. Kashef, A. I. 1986, Groundwater Engineering, McGraw-Hill, New York. Katwatni, T., 1980, Behaviour of saltwater intrusion in layered coastal aquifers, Proc. 3rd Int. Conf.

36

A.G. BOBBA

Finite Element Water Resour., Univ. of Mississippi, Oxford, Miss. Vol 2, pp. 199-208. Kohout, F. A., 1960, Cyclic flow of saltwater in the Biscayne aquifer of southeastern Florida, J. Geophys. Res. 65, 2133-2141. Kozeny, J., 1953. Hydraulik, Springer, Vienna. Lee, C. H., and Cheng, R. T., 1974, On seawater encroachment in coastal aquifers, Water Resour. Res. 10, 1039-1043. Mercer, J., Larson, S., and Faust, C., 1980, Finite-Diifference model to simulate the areal flow of saltwater and freshwater separated by an interface, U.S. Geol. Surv., Open-File Rep. 80-407. Mualem, Y., 1973, Interface refraction at the boundary between two porous media, Water Resour. Res. 10, 1207-1215. Mualem, Y. and Bear, J. 1974, The shape of the interface in steady flow in a stratified aquifer, Water Resour. Res. 10, 1207-1215. Muskat, M., 1946, The Flow of Homogeneous Fluids Through Porous Media, McGraw-Hill, New York. Pandit, A. and Anand, S. C., 1984, Groundwater flow and mass transport by finite elements - A parametric study, Proc. 5th Int. Conf. Finite Elements Water Resour., Univ. of Vermont. Pinder. G. E and Cooper, H. H. Jr., 1970, A numerical technique for calculating the transient position of the saltwater front., Water Resour. Res. 6, 875-882. Pinder, G. F. and Frind. E. O., 1972, Application of Galerkin's procedure to aquifer analysis, Water Resour. Res. 8, 108-120. Pinder, G. F. and Gray, W. G., 1977, Finite Element Simulation in Surface and Subsurface Hydrology, Academic Press, New York. Polubarinova-Kochina, E YA, 1962, Theory of Groundwater Movement (trans. J. M. R. De Wiest), Princeton University Press, Princeton. Reddel, D. I. and Sunada, D. K., 1970, Numerical simulation of dispersion in groundwater aquifers, Hydrol. Paper 41, Colorado State University, Fort Collins. Reilly, T. E. and Goodman, A. S., 1985, Quantitative analysis of saltwater-freshwater relationships in groundwater systems - a historical perspective, Y. Hydrol. 80, 125-160. Reilly, T. E. and Goodman, A. S., 1986, Analysis of saltwater upconing beneath a pumping well, aT. Hydrol. 89, 169-204. Rumer, R. R. and Shiau, J. C., 1968, Saltwater interface in layered coastal aquifer, Water Resour. Res. 4, 1235-1247. Sa da Costa, A. A. G. and Wilson, J. L., 1980, Coastal seawater intrusion: A moving boundary problem, Proc. 3rd Int. Conf. Finite Element Water Resour., Univ. of Mississippi, Oxford, Miss., Vol.2, pp. 209-218. Sahni, B. M., 1973, Physics of brine coning beneath skimming wells, Ground Water 11, 19-24. Schmorak, S. and Mercado, A., 1969, Upconing of fresh water-sea water interface below pumping wells, field study, Water Resour. Res. 5, 1290-1311. Segol, G., Pinder, G. F., and Gray, W. G., 1975, A Galerkin finite element technique for calculating the transition position of saltwater front, Water Resour. Res. 11,343-347. Segol, G. and Pinder, G. F., 1976, Transient simulation of saltwater intrusion in Southeastern Florida, Water Resour. Res. 12, 65-70. Shamit, U. and Dagan, G., 1971, Motion of the sea water interface in coastal aquifers: A numerical solution, Water Resour. Res. 7, 644-657. Shamir, U. and Harleman, D.R.F., 1967, Numerical solution for dispersion in porous mediums, Water Resour. Res. 3, 557-581. Sherif, M. M., Singh, V.E, and Amer, A.M., 1988, A two-dimensional finite element model for dispersion (2d-FED) in coastal aquifers, J. Hydrol. 103, 11-36. Shima, S., 1969, Transient characteristics of a saltwater wedge, Proc. 13th Congr. IAHR, Vol. 4, pp. 433-440. Souza, W, R. and Voss, C. I., 1987, Analysis of an anisotropic coastal aquifer system using variabledensity flow and solute transport simulation, aT.Hydrol..92, 17-41. Taylor, C. and Huyakorn, ES., 1978, Finite element analysis of three-dimensional groundwater flow with convection dispersion, Proe. 2nd Int. Symp. Finite. Element Methods Flow Problems, Italy. Todd, D. K., 1980., Groundwater Hydrology, Wiley, New York. Vappicha, V. N. and Nagaraja, S. H., 1976, An approximate solution for the transient interface in a coastal aquifer, J. Hydrol. 31, 161-173.

SALTWATER INTRUSION IN COASTAL AQUIFERS

37

Voss, C. I., 1984, SUTRA - A finite element simulation model for saturated-unsaturated fluid density dependent groundwater flow with energy transport or chemically-reactive single species solute transport, U.S. Geol. Surv., Water Resour. Invest. Rep., 84-4269. Voss, C. I. and Souza, W. R., 1987, Variable density flow and solute transport simulation of regional aquifers containing a narrow freshwater-saltwater transition zone, Water Resour. Res. 23, 1851-1866. Wikramaratna, R. S. and Wood, W. L., 1982, On the coupled equations of the groundwater quality problem, Proc. 4th. Conf. Finite Elements Water Resour., Hannover. Wilson, J. and Sa da Costa, A., 1982, Finite element simulation of a saltwater/freshwater interface with indirect toe tracking, Water Resour. Res. 18, 1069-1080.