ISSN 0097-8078, Water Resources, 2008, Vol. 35, No. 3, pp. 274–286. © Pleiades Publishing, Ltd., 2008. Original Russian Text © M.G. Khublaryan, A.P. Frolov, I.O. Yushmanov, 2008, published in Vodnye Resursy, 2008, Vol. 35, No. 3, pp. 288–301.

HYDROPHYSICAL PROCESSES

Seawater Intrusion into Coastal Aquifers M. G. Khublaryan, A. P. Frolov, and I. O. Yushmanov Water Problems Institute, Russian Academy of Sciences, ul. Gubkina 3, Moscow, 119333 Russia Received April 20, 2007

Abstract—The hydrological, hydrophysical, and hydrochemical aspects of the interaction between marine waters and groundwater in the adjacent land territory as a manifestation of the general process of interaction between surface and subsurface waters. Most attention is concentrated on the marine water intrusion into aquifers. Recent achievements in the field of modeling seawater intrusion into coastal aquifers are reviewed. DOI: 10.1134/S0097807808030032

INTRODUCTION The interaction between seawater and groundwater of the adjacent land areas is a manifestation of the general process of interaction between surface and subsurface waters. A component of this process, in addition to submarine groundwater discharge, is the intrusion of seawater into aquifers. Under natural conditions, seawater intrusion takes place when the density of groundwater is less than that of seawater with which it hydraulically interacts. The difference between the densities of seawater and fresh groundwater determines the character and rate of seawater intrusion. The denser seawater moves along the floor of the bed, while the interface between the seawater and groundwater acquires a complicated shape. Some results of theoretical and experimental studies of the motion and interaction of liquids with different densities in aquifers are given in [1, 7]. The process of intrusion can become much more intense if large amounts of groundwater are withdrawn for domestic needs. Thus, ~600000 m3/day of fresh groundwater are withdrawn from confined aquifers in the eastern coastal areas of the Baltic Sea. The deep depression cones (40–50 m) with a radius of up to 100 km that have formed as a result extend under the bed of the Baltic Sea and create the hydraulic conditions favorable for the intensification of seawater intrusion [4, 5, 8, 14]. A transitional zone with water salinity changing from freshwater to seawater level forms in the zone of contact between fresh groundwater and salt seawater. The location and size of this transitional zone depend on the density of seawater, the rate of submarine groundwater discharge, and other characteristics. The creation of subsurface water intakes in coastal areas and the intense pumping out of fresh groundwater can intensify the process of seawater intrusion into fresh aquifers, thus creating pollution of subsurface water sources, which is difficult to eliminate. The development of a strategy for the optimal use of coastal aquifers for water supply requires preliminary

studying the regularities of seawater intrusion into coastal areas, as functions of variations in water balance elements, including pumping fresh or salt waters into or out of the aquifer. Changes in the regime and volume of pumping in or out of the aquifer is associated with some transient process of deformation and propagation of the freshwater–seawater interface. The studies involving the mathematical modeling of seawater intrusion can be divided into two groups, depending on whether a sharp interface is assumed to exist between the fresh and salt water or these waters are assumed to mix. The choice of one or another approach to modeling depends on the specific hydrogeological conditions. The recent achievements in modeling seawater intrusion into coastal aquifers are reviewed in [16]. HYDROLOGICAL, HYDROPHYSICAL, AND HYDROCHEMICAL ASPECTS OF SEAWATER INTRUSION The interaction between aquifers and the sea under natural conditions is accompanied by the discharge of fresh groundwater into the sea and simultaneous penetration of seawater into fresh aquifers. With the growing water withdrawal from subsurface water intakes in coastal aquifers, the rate of submarine freshwater discharge decreases, which can be the cause of a significant intensification of salt seawater intrusion into fresh aquifers. Thus, the intensification of groundwater development in the Southern California since 1940 (because of the increasing aridity of climate) caused a drop in groundwater head below the sea surface in many coastal aquifers. The result was an intrusion of seawater into fresh coastal aquifers and considerable economic losses [21]. The penetration of seawater into freshwater aquifers affects their water quality and makes them unusable for drinking water supply. Moreover, the water-bearing rocks become saturated with salts; as a result, the deterioration of water quality in subsurface water intakes can persist over many years,

274

SEAWATER INTRUSION INTO COASTAL AQUIFERS

even after the intrusion has been eliminated and the natural process of freshwater discharge has recovered. In addition to water withdrawal, the intrusion of marine (salt) waters into freshwater aquifers can be facilitated by engineering activity in the coastal area, such as coastal land reclamation, development of quarries, etc. [31]. A good illustration is the intrusion of saltwater into a freshwater aquifer near Silver Bluff, Florida [23]. The construction of drainage systems here caused a drop in the groundwater table by several feet and thus intensified the landward intrusion of seawater by 2–3 km [38]. Under natural conditions, saltwater intrusion is governed by the density difference between salt and fresh waters and can be intensified by tidal and setup phenomena. The active development of groundwater causes a considerable drop in its head or level, which is an important cause of saltwater intrusion toward water intakes and the salinization of the withdrawn water. The intrusion depth of salt water can exceed several kilometers [3]. The growing demand for freshwater and the further forced intensification of groundwater withdrawal in marine coastal areas enhance the urgency of the theoretical and experimental studies in this field [32]. Studies of the seawater intrusion problem are being carried out along the following lines: analogous, physical, and mathematical modeling of intrusion processes; groundwater regime in coastal zones; seawater intrusion into stratified beds; optimal fresh groundwater development in the presence of its contact with seawater and the development of special engineering measures aimed to prevent the intrusion. The richest experience in studying intrusion in hydrogeological examinations of coastal aquifers has been gained in the United States, Japan, Israel, and the Netherlands. In the countries of the former USSR, the problem of seawater intrusion is of importance for some areas in the Baltic region, Kamchatka, and the coasts of the Black and Caspian seas, where seawater shifts toward coastal water intake structures [3, 4, 6]. In studying the penetration of salt seawater into freshwater aquifers, particular attention is paid to the location of the fresh–salt water interface and the prediction of water quality dynamics in water intakes operated in marine coastal zones. The intrusion of seawater into aquifers is a complex hydrodynamic process of joint motion of fresh and salt waters, which have different density and other physical properties. The processes taking place in the transition zone from fresh water to salt seawater include dispersion, as well as liquid mixing due to diffusion and advection. The place, shape, and length of the dispersion zone depend on many factors, including the density ratio of salt and fresh waters, the freshwater discharge or head, dispersion characteristics of the aquifer, etc. A critical role in the quantitative description of the interaction processes between fresh continental waters WATER RESOURCES

Vol. 35

No. 3

2008

275

and salt seawaters belongs to methods of mathematical modeling. The theoretical principles of the mathematical modeling of hydropysical processes in aquifers are represented in a number of monographs, among which the basic works of P.Ya. Polubarinova-Kochina [11] and G. Bear [22] are worth special mentioning. The general class of groundwater systems considered in intrusion problems consists of a saturated porous medium containing a liquid with varying density, depending on the concentration of dissolved salts. Systems of this type can be met in both coastal and inner continental areas. In relatively homogeneous porous media, the denser salt water commonly remains isolated from the overlying freshwater. However, a mixing or dispersion zone will commonly form between the two liquids. In this zone, part of salt water mixes with fresh water and moves toward the sea, thus causing the motion of seawater toward the mixing zone. The result is the circulation of salt water in coastal aquifers of large thickness. The thickness of the mixing zone can vary within wide limits. As early as the beginning of the past century, W. Badon-Ghyben [20] and A. Herzberg [34] independently established that salt seawater can penetrate far into freshwater aquifers and the occurrence depth salt– fresh water interface is ~40 times greater than the freshwater level in an unconfined aquifer. They used the condition of hydrostatic equilibrium Ä = Ç in points Ä and Ç, where Ä = ρcgH and Ç = ρfg(H + h) are overall hydrostatic pressures to derive the so-called Ghyben– Herzberg condition H = ρf h/(ρc – ρf), where h is the elevation of freshwater level above sea level; H is the depth of the interface between waters below sea level; ρf, ρc are the densities of fresh and salt waters, respectively. Assuming ρf = 1.00, ρc = 1.025, we obtain H = 40h. The hydrostatic theory of Ghyben–Herzberg does not take into account groundwater motion; therefore its level and the interface should have zero slope. However, under natural conditions, because of recharge, evaporation, and discharge, groundwater is in permanent motion, and fresh groundwater discharges (seeps) below sea level as well [4, 5, 7, 45]. Because of this, the position and the shape of the fresh–salt water interface should be governed by the conditions of their dynamic interaction. Nevertheless, Ghyben–Herzberg condition can be successfully applied for quasihydrostatic conditions at a sufficient distance from the shore where freshwater flow is almost horizontal. Detailed studies showed that Ghyben–Herzberg formula, along with the Dupuit hypothesis regarding the horizontal flow, is applicable to the calculation of a salt wedge under a dimensionless criterion B1 = π(H + h)Kf(ρc – ρf)/Q0ρf > 8, where Q0 is freshwater discharge per unit width of the bed, Kf is the hydraulic conductivity in the freshwater

276

KHUBLARYAN et al.

Table 1. Classification of salt waters by dissolved salt concentrations Waters

Concentration, ppt

Slightly saline Moderately saline Very saline Brines

1–3 3–10 10–35 >35

Table 2. Seawater properties at different salt concentrations DisDynamic solved Salinity, Chloridity, Density, viscosity, fraction ‰ ‰, at 20°C kg/l, at 20°C centipoise mass, % 0.00 0.50 1.00 2.00 3.00 3.50 4.00 5.00 10.00

0.00 4.94 9.92 19.89 29.89 34.84 39.82 49.79 99.83

0.00 2.72 5.48 11.00 16.53 19.29 22.05 27.53 55.12

0.9982 1.0019 1.0057 1.0132 1.0207 1.0245 1.0283 1.0358 1.0738

1.002 1.010 1.018 1.036 1.057 1.070 1.080 1.103 –

TRANSIENT INTRUSION PROBLEM IN HETEROGENEOUS BEDS IN THE REGIME OF ADVECTIVE DIFFUSION OF SALTS

Table 3. Major ions contained in seawater Ion Chlorides Na Sulfates Mg Ca K Bicarbonates Br Bromic acid Sr Others Total

such effect is the intrusion of salt seawater into freshwater coastal aquifers. Thus, the density of natural water under natural conditions varies within wide limits from 0.9982 at 20°ë for clear fresh water to >1.345 kg/l for brines [7, 24]. The classification of water from slightly saline to brines proposed by R. Krieger [39] divides salt waters into four major groups, depending on the concentration of dissolved solutes (Table 1). Water viscosity, as well as its density, is strongly dependent on the type and the amount of solids dissolve in it. Variations in some physico–chemical properties of natural waters as functions of their salinity are given in Table 2. The most common density of the surface seawater layers amounts to 1.022–1.028 g/l; these values depend on both water temperature and the concentration of dissolved solids. The kinematic viscosity of such seawater amounts to 1.826 at 0°ë and 1.049 centistokes at 20°ë. Commonly seawater contains ~34.5 ppt of dissolved solids. The major ions that occur in water are given in Table 3. Subsurface and lacustrine natural salt waters can differ from seawater in terms of chemistry, and their density and viscosity also are strongly dependent on the concentration of dissolved chemical components.

Ion content %

ppm

55.04 30.61 7.68 3.69 1.16 1.10 0.41 0.19 0.07 0.04 0.01 100.00

18980 10556 2649 1727 400 380 140 65 26 13 2 34483

zone [22, 48]. The respective error in determining the depth of the interface does not exceed 5%. This deviation increases with approaching seashore. The penetration of marine salts into the coastal hydrogeological system changes the density and other physico–chemical properties of groundwater and can have a significant effect on its dynamics. An example of

Depending on the hydrogeological conditions, the transient response of groundwater to changes in the natural regime of subaqueous discharge can last for decades. The reliable assessment of the economic consequences and the environmental effect of a large-scale intrusion of seawater into fresh aquifers requires detailed data on the duration of the transient state of the hydrological system, as well as the characteristics of the stationary dynamic equilibrium to which the system tends. The first solution of a transient problem of seawater intrusion taking into account the advective dispersion of salt was obtained by G. Pinder and H. Cooper [44] with the help of the method of characteristics. The transient state was shown to tend to the steady-state solution of H. Henry [33] at two different initial conditions, though this steady-state solution was not reached because of too large computation time. The assumption about the constancy of dispersion coefficients made in [44] is approximate because of the dynamic character of seawater intrusion into a confined aquifer. A typical feature of this process is that the fresh and salt waters in the lower part of the aquifer move toward one another and under dynamic equilibrium they form domains with low (zero) velocity (stagnant zones). Since the only mechanism of dispersion at the absence of advection of the fluid is molecular diffusion, the dependence WATER RESOURCES

Vol. 35

No. 3

2008

SEAWATER INTRUSION INTO COASTAL AQUIFERS

of dispersion coefficients on the velocity of liquid motion should be taken into account[51]. The first transient solution to intrusion problem on the basis of the dispersion coefficient depending on the flow velocity was obtained in 1975 by G. Segol [47] with the use of the finite element technique. The characteristics simultaneously calculated in each node of a finite-element grid included the pressure and two components of velocity; this allowed the author to obtain a continuous velocity field (the so-called three-equation system). The transient solution obtained by Segol also tends to the steady-state solution of Henry but does not reach it. A more efficient numerical method for the solution of solute transport problem in variable-density liquids was developed in [30] as applied to confined isotropic aquifers. The governing equations for the liquid transport equations are Darcy law; liquid continuity equation; dissolved salt conservation equation; and the equation of state relating liquid density and solute concentration. On the other hand, field and laboratory studies show that the depth of seawater intrusion into the shore is strongly dependent on the hydrological characteristics (permeability) of water-bearing rocks, which are commonly anisotropic and heterogeneous. Thus, the length of intrusion is negligible in thin low-permeability beds and is large (up to several kilometers) in thick, highly permeable beds [27, 31, 32]. Zones of intrusion in multilayer bed systems are very diverse. In stratified two-layer beds with a more permeable top layer, salt water penetrates further landward, in contrast to a bed with the same thickness but with a more permeable bottom layer. In multilayer systems that can be represented as several water-bearing layers, separated by intermediate low-permeability layers, an intrusion tongue will form in each layer, but the hydrodynamic analysis of the situation becomes more complicated if water exchange is possible between permeable layers by seepage through semipermeable layers. Depending on the location of each layer, the specific discharge of freshwater in it, and the thickness and permeability of water-bearing rocks, the extent of intrusion in individual layers will vary. A characteristic example of an intrusion zone with such structure is Long Island [41]. A multiwedge system was found to exist in coastal aquifers in Israel [46]. An attempt of theoretical consideration of multiwedge intrusion zones was made in [50] where intermediate layers were assumed to be completely impermeable. An intrusion zone of seawater into a two-layer aquifer, divided by a semipermeable interlayer, is considered in [25] in a hydrostatic approximation. Aquifers can contain isolated geological formations (structures) with different permeability. In such cases, the hydraulic conductivity in Darcy’s law varies in the horizontal direction as well. Such beds are commonly referred to as heterogeneous. WATER RESOURCES

Vol. 35

No. 3

2008

277

The most general mathematical model for studying the process of seawater intrusion into coastal freshwater aquifers is a system of groundwater flow equations and advection diffusion of dissolved salts [13, 17–19]. The model of intrusion discussed below is based on a generalized law of liquid flow in a porous medium, with liquid density depending on the concentration of dissolved solids. The proposed model can describe not only confined, but also unconfined groundwater flows, as well as heterogeneous and anisotropic aquifers. Let us introduce freshwater head ϕ(x, y, t) as P ( x, y, t ) ϕ ( x, y, t ) = ---------------------- + y, ρ0 g

(1)

where ρ0 is freshwater density, g is acceleration due to gravity, y is the vertical coordinate relative the reference line, p(x, y, t) is the pressure in point (x, y), x is the horizontal coordinate. Generalized Darcy’s law for groundwater flow of a variable-density liquid in an anisotropic bed has the form TK q = – ----------- ( ∇ p + ρg ), µ(ρ)

(2)

where q (qx, qy) are components of liquid flow in directions éï, éì, µ(ρ) is viscosity, ρ is liquid density, ∂ ∂ ∇ = ⎛ , ⎞, ⎝ ∂ x ∂ y⎠ ⎛ K K ⎞ TK = ⎜ xx xy ⎟ ⎝ K yx K yy ⎠

(3)

is the permeability tensor of water-bearing rocks, which in the case of coinciding principal axes of anisotropy with the axes éï, éY of the coordinate system will take the form ⎛ K O ⎞ TK = ⎜ x ⎟. ⎝ O Ky ⎠

(4)

Substituting relationship (1) into Darcy’s law (2), we transform it to Tkρ 0 g - [ ∇ϕ + ρ r j ], q = – -------------µ(ρ)

(5)

ρ where j is the unit vector along OY axis, ρr = ----- – 1. ρ0 The continuity equation for a variable-density liquid has the form [22] ∂p ∇ ( ρq ) = – ρ ( α + nβ ) ------ , ∂t

(6)

where α and β are the compressibility factors of the porous medium and the liquid, respectively; n is poros-

278

KHUBLARYAN et al.

ity. Since in the case of seawater intrusion, groundwater has a direct hydraulic connection with the water of an open water body, the compressibility of the porous medium and the liquid can be neglected; now the continuity equation will take the form ∇ ( ρq ) = 0.

(7)

The conservation equation of the dissolved substance (sea salt) can be written in the reduced concentration c = C/CM ∂c ∇ ( TD∇c ) – ∇ ( V c ) = -----, ∂t

(8)

where ë is salt concentration in the mixture of the fresh and sea water, CM is salt concentration in seawater, TD is the tensor of hydrodynamic dispersion ⎛ D D ⎞ TD = ⎜ xx xy ⎟ , ⎝ D yx D yy ⎠

(9)

whose components are evaluated by the following formulas from [19] D xx = α L V x / V + α T V y / V + D M T xx , 2

2

D yy = α T V x / V + α L V y / V + D M T yy , 2

2

(10)

D xy = D xy = ( α L – α T )V x V y / V , where αL and αT are the longitudinal and transverse dispersion constants, DM is molecular diffusion, Txx, Tyy are the principal components of the tortuosity tensor, |V| is flow velocity modulus. According to the data of laboratory studies presented in [18], the relationship between the density and viscosity of a mixture of fresh and marine water, with the relative concentration of sea salts varying within the range of interest, can be approximated by the linear relationships (11) ρ = ρ0(1 + ε), µ = µ0(1 + 2.8εc),

(12)

where ε = (ρM – ρ0)/ρ0, µ0 is freshwater viscosity and ρM is seawater density. Let us introduce dimensionless variables X = x/L, TK = TK/K 0 , µ = µ/µ 0 ,

Y = y/L,

Q = Lq/Q f , 2

t = tD 0 /L , Ψ = ψ/L,

ρ = ρ/ρ 0 ,

(13)

TD = TD/D 0 ,

where Qf is freshwater discharge per unit length of the shoreline, D0 is the maximum value of the coefficients of the hydrodynamic dispersion tensor in the intrusion zone, K0 is the maximum value of the permeability coefficient, L is bed thickness. The system of governing equations (5), (7), (8), (11), (12) transforms into

∇ ( ρQ ) = 0,

(14)

TK LK 0 ρ 0 g⎞ - ∇ϕ + cj , Q = – ------- ⎛ -----------------µ ⎝ µ0 Q f ⎠

(15)

Qf ⎞ ∂c = -----, ∇ ⎛ TD∇C⎞ – ∇ ( Qc ) ⎛ --------⎝ nD 0⎠ ⎝ ⎠ ∂t

(16)

ρ = 1 + εc,

(17)

µ = 1 + 2.8εc,

(18)

where ∇ the gradient in the new coordinate system. The use of the relative concentration in Darcy’s equation instead of the density improves the calculation accuracy, since the relative concentration can vary within the limits of 0 ≤ c ≤ 1.00, and the liquid density for the conditions of seawater intrusion varies within the limits of 1.000 ≤ ρ ≤ 1.033. Substituting relationship (18) into equation (14), we obtain ( 1 + εc )∇Q + ε∇cQ = 0.

(19)

Since ε for seawater commonly does not exceed the values of 0.035–0.040, the terms containing ε can be neglected in equation (16), and the continuity equation takes the form ∇Q = 0,

(20)

where Q is determined by (15). On the other hand, according to the results of field studies [40], the value M = K0εLρ0g/µ0Qf in coastal aquifers in the zone of intrusion can have an order of O(1), because of which the second term in the righthand part of (15) cannot be neglected in studies of intrusion processes. Equations (14)–(16) are related by the concentration of dissolved salts and liquid velocity and hence cannot be solved independently. Let us formulate the boundary conditions for this system of equations in head ϕ(x, y, t) and in sea salt concentration c(x, y, t). The boundary conditions for equation (14) can be of the first (Dirichlet) or the second (Neumann) type. The first-type boundary condition has the form (21) ϕ = ϕ0 on Γ1, where ϕ0 is the specified head. The second-type condition, expressed in terms of flow qf in the direction of the normal to the boundary of the flow domain has the form (22) qf = q0 on Γ2, where q0 is the specified liquid flow, and Γ1 + Γ2 = Γ is the boundary of the flow domain. If the aquifer is limited from above by a permeable aquiclude, the boundary condition on it will be characterized by water flow proportional to the head differenWATER RESOURCES

Vol. 35

No. 3

2008

SEAWATER INTRUSION INTO COASTAL AQUIFERS

tial in this aquiclude at the stationary flow of liquid in it. According to [30], the flow through a permeable aquiclude can be written as ϕ' – ϕ q f K 'y ⎛ -------------- + ρ 'r⎞ on Γ 2 , ⎝ b' ⎠

(23)

where K' is the vertical hydraulic conductivity in the low-permeability aquiclude, b' is its thickness, ϕ' is the known equivalent freshwater head on the upper boundary of the aquiclude, ϕ is the unknown head at the upper boundary of the aquifer, ρ r' is the relative density of the liquid in the low-permeability aquiclude. Equation (23) can be considered as a third-type boundary condition. The boundary conditions for equation (16) can have either first or third type. The first-type condition has the form (24) Ò = Ò1 on Γ3, where Ò1 is specified concentration. The third type of the boundary condition has the form q ∂c ⎞ ⎛ v c – D ------m i = -----0 c 0 on Γ 4 , i ij ⎝ ∂ x y⎠ n

(25)

where mi is the ith component of the external normal to boundary Γ4, vi is the projection of the liquid velocity on xi-axis, x1 = x, x2 = y, q0 is the specified flow in the direction of normal, c0 is the specified concentration, Γ3 + Γ4 = Γ is the boundary of the flow domain. Equation (25) is used for the boundary with the external inflow, when the flow of liquid with the specified solute concentration is known. Equation (25) can also be applied to the boundary with the outward flow, though in this case, q0 and c0 are unknown functions. For obtaining the boundary condition required for the solution of the boundary problem, the advective fluxes along the normal to the boundary are often assumed to be equal on both sides of the boundary. In such cases, boundary condition (25) degenerates into the secondtype boundary condition ∂c --------- = 0 on Γ 4 , ∂m 1

(26)

where m1 is the direction of the external normal. This boundary condition is also applicable in the case of an impermeable boundary. The computational algorithm and the computer code for calculations using finite-difference technique were developed for a rectangular domain (0 < x < lx, y = 1) for confined and unconfined groundwater flow [19]. Let us write the appropriate boundary conditions for such problems in more detail. On the aquiclude (floor of the bed), y = 0, 0 < x < lx, flow qy = 0 or, in accordance with (15), ∂ϕ ------ + εc = 0. ∂y WATER RESOURCES

Vol. 35

(27) No. 3

2008

279

A nonzero flow can also be specified on this boundary qy = f1(x, t). (28) On the roof of the bed y = 1, 0 < x < lx, the flow is equal either to zero ∂ϕ q y = ------ + εc = 0, ∂y

(29)

or a specified function qy = f2(x, t). On the free surface of the unconfined groundwater flow, pressure ê = 0, i.e., ϕ = y = f. On the vertical segments of the landward boundary, we specify either the head ϕ = f3(y, t), or discharge q = f4(y, t). On the seaward boundary below the water level, dm (x = 0, 0 < y < dm), the hydrostatic pressure distribution was specified P = ρmg(dm – y), ρ ϕ = -----m- d m – εy. ρ0

(30)

For an unconfined flow in the seepage zone x = 0, dm < y < f(0, t), pressure p = 0, i.e., ϕ = y. The boundary conditions for sea salt transport equation (16) are as follows x = lx,

0 < y < f ( l x, t ) or 0 < y < 1,

c = 0 ( for confined aquifer ); x = 0,

0 < y < R,

( q x > 0 ),

c = 1;

R < y < f ( 0, t ) or R < y < d m , ( q x < 0 ), 0 < x < lx,

∂c ----- = 0; ∂t y = 0,

(31)

y = 1,

∂c ----- = 0 ( for confined aquifer ). ∂y Here, R is the point on the interface between groundwater and seawater where the flow direction changes (the so-called return point). Below this point, salt seawater flow into the freshwater aquifer, while above this point, groundwater discharges into the sea. On the free surface, in the presence of infiltration recharge (for isotropic media) the condition has the form ∂c D ------- = E N ( c – c m ), ∂N

(32)

∂ is the derivative ∂N along the internal normal to the boundary, EN is the infiltration recharge per unit free surface, i.e., where D is diffusion coefficient,

E N = inf

∂ϕ 2 1 + ⎛ ------⎞ ⎝ ∂x ⎠

1 – --2

.

(33)

280

KHUBLARYAN et al.

The derivation of this condition and the condition above return point R are given in [18, 19]. The initial condition for the problem is the initial concentration distribution c(x, y, 0) = c0(x, y), and, in the case of unconfined beds, the position of free surface f(x, 0) = f0(x). The free surface moves in accordance with the equation ∂ϕ ∂ f ∂f ∂ϕ n ------ = – K y ⎛ ------ + εc⎞ + K x ------ ------ + in, ⎝ ⎠ ∂x ∂x ∂t ∂y where in(x, y, t) is the rate of infiltration recharge onto the free surface with salt concentration of ëin. The boundary problem formulated above was numerically solved on a computer. All partial differential equations were solved by using a splitting technique with respect to spatial variables. This procedure was used to calculate different variants in both confined aquifer and in aquifers with free surface [18]. The effect of hydraulic conductivity, the discharge value f4, and dispersion coefficient on the position of the zone of seawater penetration was evaluated. The results of model calculations were compared with field data [38]. Modeling of intrusion into a stratified confined aquifer was considered as an example. The position of isochlor ë = 0.5 under steady-state regime was presented in [18] for four cases: a homogeneous bed and an aquifer with a separating interlayer, where the hydraulic conductivity of the interlayer is less than that in the top and bottom layers by a factor of 10, 100, and 1000, respectively. A distinct trend toward a decrease in the mutual effect of intrusion processes in the top and bottom layers can be seen. When the interlayer is virtually impermeable, the position of the isochlor ë = 0.5 in the top and bottom layers is the same, i.e. the processes in these layers are independent. SEAWATER INTRUSION UNDER THE REGIME NONMIXING LIQUIDS The object of quantitative analysis in intrusion processes is the interaction of fresh and salt waters. The physical system consisting of a mixture of salt and fresh waters in a water-bearing porous medium and involving a three-dimensional natural flow and a gradual change in the sea salt concentration across the mixing zone is very complex and its full descriptions are few. Instead, various simplifying assumptions are commonly made, thus allowing a quantitative description of fresh—salt water interaction to be developed with sufficient detail [15, 32, 37]. The most important simplifying assumption regards the fresh–salt water mixing. Under concrete hydrological conditions, these liquids, which mix well, can be regarded as immiscible and separated by a sharp boundary. This approximation is widely used in intrusion problems and is helpful in their mathematical modeling.

The assumption that a sharp boundary exists between fresh and salt waters is especially realistic in large-scale studies that cover considerable land areas, including island territories [26, 42]. In the presence of a sharp boundary, the three-dimensional equations of salt and fresh water motion can be integrated along the vertical coordinate and reduced to a system of twodimensional equations. Let us consider the behavior of the salt–fresh water interface in a three-dimensional case in confined and unconfined aquifers under the following basic assumptions: the boundary line between the waters is sharp, no dispersion takes place; Dupuit assumption is valid for the motion of fresh and salt water, i.e., the hydraulic head does not vary along the depth; the hydraulic conductivity and specific yield are constant along the depth of the bed, i.e., stratified aquifers are not considered; the aquifer floor is impermeable. Under these assumptions, the main equations of motion of fresh and salt waters, integrated over the vertical coordinate, have the form [43] ρc ∂ϕ c ∂ϕ ϕf ⎛ ∂------------S – ρ *c --------c⎞ f b c -------- + n ⎝ ρ * f ∂t ∂t ∂t ⎠ ρf –

µ f ρ c ∂ϕ c⎞ ∂ ⎛ µ f ρ c ∂ϕ c⎞ ∂⎛ - -------- – - -------b K ---------b K ---------∂ x ⎝ c fx ρ f µ c ∂ x ⎠ ∂ y ⎝ c fy ρ f µ c ∂y ⎠

(34)

– Q wc δ ( x – x i, y – y i ) = 0, ∂ϕ ∂ϕ ∂ϕ ∂ϕ S f b f --------f – n ⎛ ρ *f --------f – ρ *c --------c⎞ + αn --------f ⎝ ∂t ∂t ∂t ⎠ ∂t –

∂ϕ ∂ϕ ∂ ∂⎛ b f K fx --------f ⎞ – ⎛ b f K fy --------f ⎞ ⎝ ⎠ ⎝ ∂x ∂y ∂y ⎠ ∂x

(35)

K' – Q wf δ ( x – x i, y – y i ) + αI – ( 1 – α ) ----- ( ϕ f – ϕ' ) = 0, b' where ρ *f = ρf /(ρf – ρc), ρ *c = ρc/(ρf – ρc), b is liquid layer thickness, ρ is liquid density, b' is confined-bed thickness, ϕ' is the head in the overlaying low-permeability layer, ϕ is hydraulic head, Kf is hydraulic conductivity for freshwater, K' is hydraulic conductivity in the confined aquifer, Qw is water intake rate, I is infiltration rate, Sf is specific yield in the freshwater zone, t is the time, z is interface level, δ is Dirac delta function, µ is dynamic viscosity, n is porosity, α = 0 is for confined bed, α = 1 is for unconfined bed, subscripts c and f refer to salt and fresh water, respectively. Equations (34) and (35) are written in the heads of the fresh ϕf and salt ϕc waters with the aim to enable the use of a constant-head boundary. The same equations can be written in ϕc and the level of the interface z or z and ϕf by using the relationship ϕ f ρ f – ϕc ρc -, z = ---------------------------ρ f – ρc WATER RESOURCES

Vol. 35

(36) No. 3

2008

SEAWATER INTRUSION INTO COASTAL AQUIFERS

which expresses the equality of the pressure of both liquids on the interface line. Thus, the resulting model (34)–(35) is quasi-three-dimensional. In accordance with (36), head changes in fresh and salt water can be mutually compensated, so that the position of the interface in this case can remain virtually unchanged. The presentation of the equation of motion of salt water in the form (34) creates additional numerical problem. In the presence of only one liquid in a part of the aquifer, the determination of the equation reduces to the equation of motion and an equation reflecting the fact that the time derivative of the interface position is zero. The use of an arithmetic mean or upstream weighing for the assessment of internodal coefficients may result in the appearance of nonzero discharges for both liquids, because of which the position of the interface can be above the upper boundary of the aquifer or below its lower boundary. As the result, the change in the interface position, which is required to counterbalance the discrepancies in the flow balance for a block can cause an artificial increase or decrease of the mass within the block [43]. Clearly the existence of this mass disbalance is an intrinsic feature of two-dimensional models with vertical averaging, and this fact should be taken into account in the analysis of the source of computation errors. The presence of wells in the aquifer makes the solution of the intrusion problem even more complicated. In the general case, a well can tap both the salt and fresh water layers, so that both liquids will be pumped out. The pumping rate of each liquid is partially determined by the position of the fresh–salt water interface in the well. The simplest approximation for the calculation of the volumes of liquids being pumped out is the linear proportionality of these volumes to the lengths of the segments occupied by the respective liquid in the well. These values can be explicitly specified in the model with the help of position of the interface in the well at the beginning of each time step. However, since the interface position varies within a time step, such explicit approach yields errors in mass balance, which can cause instability of calculation. In order to obtain a stable solution irrespective of the time step, we have to use an explicit scheme, which takes into account the volumes of pumped-out liquids as a function of the heads of fresh and salt waters [43] ∂Q Wf ∂Q Wf k+1 k - ∆ϕ f + ----------- ∆ϕ c , Q Wf = Q Wf + ----------∂ϕ f ∂ϕ c k

k

∂Q WC ∂Q WC - ∆ϕ f + ------------ ∆ϕ C , + -----------∂ϕ f ∂ϕ C k

k+1 Q WC

=

k Q WC

(37)

k

(38)

where the derivatives are calculated in the beginning of the time step, k is time step number, QWF and QWC are the pumped-out volumes of the fresh and salt water, respectively. WATER RESOURCES

Vol. 35

No. 3

2008

281

The use of coastal groundwater basins as reliable sources in water management systems requires the development of quick methods for the prediction of water mass behavior under different conditions of the intake of fresh and salt waters. One of the main problems is the long-term prediction of the direction and velocity of motion of salt seawater masses in an aquifer, filled with freshwater, under the effect of pumping out, artificial recharge, time-variable infiltration recharge, and other factors. For this purpose, analytical methods are used in some simplified hydrogeological situations to solve the problems of transient motion of the interface between fresh and salt waters in coastal aquifers. An analytical solution of the problem of motion of a rectilinear interface in a homogeneous isotropic confined aquifer with a constant thickness is constructed in [36]. Suppose that in moment t = 0, a vertical wall in point x = 0 separates salt (x < 0) and fresh (x > 0) water. The wall is removed at t = 0, and, because of the difference between the densities of liquids, the interface starts moving, while remaining rectilinear. The interface is described by the equation D x ζ ( x, t ) = ---- 1 + ---------- , 2 L(t )

(39)

∆ρkD where L(t) = A t, A = ⎛ ---------------⎞ , ∆ρ = ρc – ρf, L(t) is ⎝ πρ f ⎠ the intersection point of the interface and the bed floor. As can be seen from (39), the interface remains rectilinear and always passes through point M (x = 0, y = −D/2). Under Dupuit assumption that the hydraulic head is constant in the vertical direction in both fresh and salt water domains, the two-dimensional intrusion problem can be reduced to a one-dimensional problem, where all required variables depend only on x coordinate. This reduction of a two-dimensional to a onedimensional problem is the more accurate, the greater is the length of the intrusion zone relative to the aquifer thickness. Dupuit assumption commonly works well under natural hydrogeological conditions, and the errors in the evaluation of effective porosity and hydraulic conductivity of water-bearing rocks and the geometry of bed floor surpass the errors resulting from the use of Dupuit assumption. The equations of motion of fresh and salt water within the intrusion zone 0 ≤ x ≤ L(t) are combined in accordance with the pressure continuity condition on the interface and have the form [10, 12, 48, 49] ∂S ∂S ∂ζ ∂ αn n ----- + n c ----- – K ( ζ + αS ) ----- = R, ∂x ∂t ∂t ∂ x f ∆γ ∂ζ ∂ ∂ γf n c ----- + – ------ζ⎞ K c ( D – ζ ) ⎛ -----S γC ⎠ ∂t ∂ x ∂ x⎝ γ C For L(t) ≤ x ≤ B

= 0,

(40)

(41)

282

KHUBLARYAN et al.

∂S ∂S ∂ K è ( D + αS ) ----- = R, αn n ----- – ∂x ∂t ∂ x

(42)

where S(x, t) is the elevation of the phreatic surface (or the hydraulic head in the case of a confined bed) above sea level; ζ(x, t) is the distance from the sea level to the interface; nn is the effective porosity in the case of motion of the phreatic surface; nc is the effective porosity in the case of motion of the interface; Kf and KC are hydraulic conductivities in the fresh and salt water zones, respectively; D(x) is the thickness of the waterbearing layer from the sea level to the impermeable layer; R(x, t) is the resulting inflow into the aquifer from above (infiltration plus pumping in minus pumping out); ∆γ = γC – γè, γ is the specific weight of water; α is a parameter (α = 1 for unconfined aquifer, and α = 0 for confined bed). Equation (40) reflects the flow continuity condition in the freshwater domain above the interface, (41) and (42) reflect the continuity condition of fresh and salt water beyond the intrusion domain (x ≥ L(t)), respectively. The initial conditions are specified by the position of the interface ζ and the free (phreatic) surface S (or the distribution of the head for a confined aquifer) in some initial moment t = t0 ζ = ζ ( x, t 0 ), S = S ( x, t 0 ),

0 ≤ x ≤ L ( t 0 ),

(43)

0 ≤ x ≤ B.

On the sea shore at x = 0, Dupuit approximation requires that for all t > 0, S(0, t) = ζ(0, t) = 0. (44) However, under real conditions, there are some ooze zones here, i.e., S(0, t) and ζ(0, t) have some finite values. These values can be derived from approximate analytical solutions (e.g., by hodograph method), or laboratory experiments. Therefore, in the general case (45) S(0, t) = S0(t), t ≥ t0, ζ(0, t) = ζ0(t), t ≥ t0.

(46)

Boundary condition (46) can be replaced by the solution Q Of ( t ) -, ζ ( 0, t ) = ------------------------K f ( ∆γ /γ f )

(47)

where QOf(t) is the freshwater discharge into the sea per unit length of the shoreline. The condition on the lower end of the interface on the aquifer floor has the form ζ[L(t), t] = D[L(t)]. (48) On the boundary of the aquifer at x = B on the side of the land, freshwater discharge can be specified Q Bf

∂S = – K f ( D + αS ) ----∂x

,

(49)

x=B

or this boundary can be assumed to be impermeable

∂S ----∂x

= 0.

(50)

x=B

The boundary conditions (45)–(50) and the initial conditions (43) are sufficient for the correct solution of equations (40)–(42). The distribution R(x, t) is assumed to be known. PREDICTION OF SEAWATER INTRUSION INTO FRESHWATER AQUIFERS The development of a strategy of optimal use of coastal aquifers for water supply requires the preliminary examination of the regularities of seawater intrusion into the shore as a function of variations in water balance elements, including the pumping of fresh or salt waters into and out of the aquifer. Variations in the regime and volumes of such pumping correspond to a transient process of deformation and propagation of the fresh–salt water interface. As it was mentioned above, the general mathematical model for the description of seawater intrusion into freshwater coastal aquifers consists of a system of equations of groundwater flow and advective diffusion of dissolved salts. However, the application of such models is limited because of the computational problems and the virtual impossibility of detailed calibration [43]. When there is no mixing of liquids, the threedimensional equations of motion of salt water and fresh water can be integrated over the vertical coordinate to obtain system of quasi-three-dimensional equations of motion (40) and (41) in the zone of interaction of these waters. However, even with these simplifying assumptions, the transient problem of seawater intrusion into the shore is still nonlinear and even its numerical solution is difficult [43]. In the linear approximation, the solution of the replacement of a light liquid by a denser one in confined aquifers under different boundary conditions is proposed in [2]. A large-scale seawater intrusion into confined coastal aquifers in one-dimensional formulation at different regimes of water intake operation is considered below in detail by reducing a system of flow equations of salt and fresh waters to a nonlinear parabolic differential equation. A scheme of confined aquifer accepted in calculation is described in [10]. The flow equations of incompressible fresh and salt waters under the simplifying assumptions described in the previous section can be written in the form ∂h ∂η ∂⎛ ηK f --------f ⎞ – n ------ = 0, ∂t ∂x ⎠ ∂ x⎝

(51)

∂h ∂η ∂ ( m – η )K s -------s + n ------ = 0, ∂x ∂t ∂x

(52)

where m is aquifer thickness; t is time; η is the ordinate of the interface; n is effective porosity of water-bearing WATER RESOURCES

Vol. 35

No. 3

2008

SEAWATER INTRUSION INTO COASTAL AQUIFERS

rocks; Kf, Ks are hydraulic conductivities of fresh and salt waters, respectively. The hydraulic heads of fresh hf and salt hs waters are determined by the relationships Pf - – y, h f = -------ρfg

P h s = -------s- – y, ρs g

(53)

where Pf, ρf, Ps, ρs are the pressure and mass density of liquids in the flow domains of fresh and salt waters, respectively; g is acceleration due to gravity. The motion of freshwater in the zone x > L(t) takes place under elastic regime ∂h 1 ∂h ---------2-f = --- --------f , a ∂t ∂x 2

(54)

where a is the pressure conductivity factor. By using Lembke method to take into account the elastic flow regime in this zone, one can easily see that the freshwater head is a linear function of coordinate x, and the flow rate of the liquid, and hence the discharge of freshwater per unit width of the bed at x ≥ L depends only on time. The case of v(t) ≥ 0 is considered below. The conditions of pressure continuity on the interface between waters yield the following relationships (55) (1 + ε)hs = hf – εh, ε = (ρs – ρf)/ρf. The motion of the interface between fresh and salt waters in confined and unconfined coastal aquifers was examined by numerical methods in [48]. For confined aquifers, a finite-difference analogue of system of equations (51), (52), and (55) was used in the intrusion zone and elastic flow equation (54) was used outside this zone. The initial condition was specified by the position of the interface η at t = t0: η0 = η(x, t0), 0 ≤ x ≤ L(t0). On the seashore at x = 0, Dupuit approximation requires the condition η(0, t) = 0 for all t. The condition on the lower end of the interface from the landward side has the form η(L, t) = m. After each time step ∆t, a new value of the length of the intrusion zone L(t) is found by linear extrapolation of values η(x, t) in the last two points of the finite-difference grid until this line crosses the bed floor. Calculations [48] are in good agreement with experimental data obtained from Hele–Shaw model. Let us transform system of equations (51) and (52) in the zone of joint motion of fresh and salt waters to a single parabolic equation. Summing these equations and carrying out one-time integration with respect to x, we obtain the equation ∂h ∂h K f η --------f + ( m – η )K s -------s = B ( t ), ∂x ∂x

(56)

where B(t) is a constant of integration with respect to the space coordinate. From the condition of the equality of freshwater flows from the left and the right of x = L(t), it follows that B(t) = v(t) in equation (56). Flow v(t) is positive if WATER RESOURCES

Vol. 35

No. 3

2008

283

freshwater flows toward the sea. Using (55) we eliminate hs from equation (56) and obtain ∂h ( m – η ) ∂h ∂η K f η --------f + ------------------K s ⎛ --------f – ε ------⎞ = v ( t ). ⎝ ∂x ⎠ ∂x 1+ε ∂x

(57)

Now, solving equations (51) and (57) in derivative ∂hf /∂x, we finally receive a single equation in boundary η(x, t) K m – η ∂η v ( t ) ε ------s ⎛ -------------⎞ ------ + ---------⎝ ⎠ K Kf ∂η ∂ f 1 + ε ∂x Kf η ------------------------------------------------------ – n ------ = 0. ∂t ∂x K s ⎛ m – η⎞ η + ------ ------------Kf⎝ 1 + ε⎠

(58)

In problems of seawater intrusion, parameter ε usually does not exceed 0.035. Therefore, we can assume with a high reliability that in (58), Ks ----------------------≈ 1. K f (1 + ε)

(59)

Under this condition, equation (58) reduces to the parabolic equation Kf

∂η ∂ η ∂η v ( t )η εη ⎛ 1 – ----⎞ ------ + -------------- – n ------ = 0. ⎝ ∂t ∂x m⎠ ∂ x mK f

(60)

Introducing dimensionless variables X = x/m, ξ = η/m, K f tε -, V = v /mK f ε τ = ---------nm

(61)

we transform equation (60) to the form ∂ ∂ξ ∂ξ ξ ( 1 – ξ ) ------- + V ( τ )ξ = -----. ∂x ∂X ∂τ

(62)

An approximate solution of equation (62) was obtained in [28] by the method of characteristics and elimination from it derivative ∂2ξ/∂X2 only for jumplike variation in freshwater discharge in the aquifer ⎧ v 0 = const, V (τ) = ⎨ ⎩ v 1 = const,

τ ≤ 0, τ > 0.

(63)

This solution has the following form: in the case of v1 < v0 ⎧ 1 1 ⎪ --------- + τ ( v 0 – v 1 ), 0 ≤ τ ≤ ---------2 , 2v ⎪ 0 2v 0 ⎪ 1 1 ⎪ l ( τ ) = ⎨ 2τ – v 1 τ, ---------2 ≤ τ ≤ ---------2 , 2v 0 2v 1 ⎪ ⎪ 1 1 ⎪ --------- , -------- ≤ τ ≤ ∞, ⎪ 2v 1 2v 21 ⎩

(64)

284

KHUBLARYAN et al.

in the case of v1 > v0 1 1 ⎧ --------- + ( v 0 – v 1 )τ, 0 ≤ τ ≤ ---------------- , 2v 0 v 1 ⎪ 2v 0 l(τ) = ⎨ 1 1 ⎪ -------- , ---------------- ≤ τ ≤ ∞, ⎩ 2v 1 2v 0 v 1

⎧ v 0 = const, V (τ) = ⎨ ⎩ v 1 = const, (65)

where l = L/m. To obtain analytical solution of equation (62) at arbitrary time variations in flow rate V(t), let us average the derivative ∂ξ/∂t in this equation over coordinate X in domain 0 ≤ X ≤ l, assuming that [9] l

1 ∂ξ a 0 ( τ ) = --- ----- d X, l ∂t

∫

(66)

0

in this case, equation (62) reduces to an ordinary differential equation, which, after a single integration, yields ∂ξ ξ ( 1 – ξ ) ------- + V ( τ ) ( ξ – 1 ) = a 0 ( τ ) X + a 1 ( τ ), (67) ∂X where a1(τ) is an arbitrary function. Introducing u = ξ2/2, we rewrite equation (67) in the form a0 ( τ ) X + a1 ( τ ) du - + V ( τ ). ------- = ----------------------------------dX 1 – 2u

(68)

The solution of equation (68) at the boundary condition X = 0, u = u0 = 0 can be obtained by Picard-Lindelof technique x

uk + 1 ( X ) =

a0 ( τ ) X + a1 ( τ )

- d X + V ( τ ) X. ∫ ----------------------------------1 – 2u 0

(69)

k

Restricting ourselves to a first approximation in (69) and taking into account boundary condition X = l, u = 0.5, we obtain u = [ V ( τ ) + a 1 ( τ ) ] X + 0.5a 0 ( τ ) X ,

(70)

a0(τ) = l–2[1 – 2l(V + a1)].

(71)

2

To find a1(τ), we use the kinematic boundary condition on the floor of the bed at X = l in the form [2] dl dξ ------- = V ( τ ) + -----, dτ dX

1 dl a 1 = --- – ----- – 2V ( τ ). l dτ

(72)

The law of the propagation of the salt-water front l(τ) is determined from condition (66), from where, in the quasi-stationary case (∂ξ/∂τ ≈ 0), we can obtain the first-order ordinary differential equation in l(τ) dl 2l ----- + 2V ( τ )l – 1 = 0. dτ

(73)

When freshwater discharge in the aquifer varies step-wise

τ ≤ 0, τ > 0,

(74)

the integral of equation (73) has the form [30] 1 – 2V 1 l 0 1 - – 2V 1 ( l – l 0 ) , τ = ---------2 l 0 --------------------2V 1 1 – 2V 1 l

l 0 = l ( 0 ). (75)

The interface between the fresh and salt waters in this case is the parabola dl 2 η = 2 V ( τ ) + ----- X. dτ

(76)

Let us consider two regimes of water intake from the bed: pumping freshwater out of the bed—an increase in the intrusion zone (77) V1 = 0.5V0, V0 = 0.16; pumping freshwater into the bed (artificial freshwater recharge)—a decrease in the intrusion zone (78) V1 = 2V0, V0 = 0.08. Formula (75) provides the best approximation of the results of finite-difference modeling. By contrast to the results of [35], when the intrusion zone decreases (pumping freshwater into the bed), the propagation of the interface is more realistic and smooth. CONCLUSIONS The effect of seawater intrusion on fresh groundwater in coastal areas requires thorough examination, especially now, when its withdrawal for water supply to maritime territories increases. The developed mathematical models of salt seawater intrusion into coastal aquifers allow one to assess the extent of fresh groundwater contamination by marine salts in these regions. Depending on the concrete hydrogeological conditions and the hydrological regime of the coastal regions, the calculations of the penetration depth of saline seawater into fresh groundwater aquifers can be carried out by two models: in the regime of immiscible liquids and in the regime of advective diffusion (dispersion). The results of simulation of sea salts propagation in a coastal aquifer enables one to more rationally design the layout of groundwater intakes in coastal regions with the aim to ensure the water quality required for domestic water supply. REFERENCES 1. Bochever, F.M., Lapshin, N.N., and Oradovskaya, F.E., Zashchita podzemnykh vod ot zagryazneniya (Groundwater Protection against Pollution), Moscow: Nedra, 1979. 2. Verigin, N.N., Sarkisyan, V.S., and Shibanov, A.V., On the Determination of the Interface between Two ImmisWATER RESOURCES

Vol. 35

No. 3

2008

SEAWATER INTRUSION INTO COASTAL AQUIFERS

3. 4.

5. 6. 7.

8.

9. 10. 11. 12.

13. 14.

15.

16.

17.

18.

cible Liquids in Porous Medium, Izv. Akad. Nauk SSSR, Mekh. Zhidk. Gaza, 1973, no. 6, pp. 155–163. Gol’dberg, V.M., Seawater Intrusion into Fresh Aquifers, in Gidrologicheskie issledovaniya za rubezhom (Hydrological Studies Abroad), Moscow: Nedra, 1982. Gregorauskas, M.M., Mokrik, R.V., and Iokshas, K.K., Hydrochemical Aspects of Studying Fresh Groundwater Discharge into the Baltic Sea, Vodn. Resur., 1986, no. 4, pp. 13–22. Gregorauskas, M.M., Mokrik, R.V., and Meeris, S.A., On the Issue of Modeling Seawater Intrusion into Aquifers, Vodn. Resur., 1987, no. 4, pp. 21–30. Zektser, I.S., Present-Day State of the Problem of Studying Groundwater-Seawater Interaction, Vodn. Resur., 1983, no. 6, pp. 57–68. Zektser, I.S., Dzhamalov, R.G., and Meskheteli, A.V., Podzemnyi vodoobmen sushi i morya (Subsurface Water Exchange between the Continent and the Sea), Moscow: Gidrometeoizdat, 1984. Iodkazis, V.I. and Mokrik, R.V., Hydrogeology of the Submarine Part of the Pribaltic Artesian Basin, in Osadkoobrazovanie v Baltiiskom more (Sediment Formation in the Baltic Sea), Moscow: Nedra, 1981, pp. 35–44. Pirverdyan, A.M., Neftyanaya podzemnaya gidravlika (Petroleum Subsurface Hydraulics), Baku: Aznefteizdat, 1956. Podsechin, V.P. and Frolov, A.P., Transient Intrusion of Seawater into Coastal Confined Aquifers, Vodn. Resur., 1989, no. 2, pp. 78–84. Polubarinova-Kochina, P.Ya., Teoriya dvizheniya gruntovykh vod (Theory of Groundwater Flow), Moscow: Nauka, 1977. Frolov, A.P. and Podsechin, V.P., Modeling of Freshwater Aquifer Pollution by Seawater in Coastal Regions, in Matematicheskie problemy ekologii (Mathematical Problems of Ecology), Chita: Ekopress, 1990, pp. 42–43. Frolov, A.P. and Khublaryan, M.G., Salt Water Penetration into Fresh Coastal Aquifers, Vodn. Resur., 1986, no. 2, pp. 58–62. Frolov, A.P. and Yushmanov, I.O., Interaction between Fresh Groundwater and Salt Water at the Baltic Sea Coast, Vodn. Resur., 1998, vol. 25, no. 1, pp. 5–8 [Water Resour. (Engl. Transl.), vol. 25, no. 1, pp. 1–3]. Khublaryan, M.G. and Frolov, A.P., Studying the Problems of Saltwater Intrusion by Approximate Methods, in Tr. V Vsesoyuzn. gidrologicheskogo s”ezda (Proc. V AllUnion Hydrol. Congr.), Leningrad: Gidrometeoizdat, 1990, vol. 9, pp. 196–201. Khublaryan, M.G. and Frolov, A.P., Modelirovanie Protsessov Intruzii v Estuariyakh i Podzemnykh Vodonosnykh Gorizontakh (Modeling Processes of Intrusion in Estuaries and Aquifers), Moscow: Nauka, 1988. Khublaryan, M.G., Frolov, A.P., and Yushmanov, I.O., Modeling Submarine Discharge and Seawater Intrusion in Shelf Zones, in Tr. kongr. “Inzhenernaya geologiya shel’fa i kontinental’nogo sklona morei i okeanov mira” (Proc. Congr. “Engineering Geology of the Shelf and Continental Slope of Seas and Oceans of the World), Tbilisi: MAGI, 1988, pp. 148–149. Khublaryan, M.G, Frolov, A.P., and Yushmanov, I.O., Sea-Water Intrusion into Inhomogeneous Aquifers, WATER RESOURCES

Vol. 35

No. 3

2008

19.

20.

21.

22. 23.

24. 25.

26.

27.

28.

29.

30.

31.

32. 33.

34.

35.

36.

285

Vodn. Resur., 1996, vol. 23, no. 3, pp. 288–292 [Water Resour. (Engl. Transl.), vol. 23, no. 3, pp. 262–266]. Khublaryan, M.G., Churmaev, O.M., and Yushmanov, I.O., Examination of the Hydrodynamic Problem Flow and Advective Diffusion in Heterogeneous and Anisotropic Porous Media, Vodn. Resur., 1984, no. 3, pp. 23–29. Badon-Ghyben, W., Nota in Verband Met De Voorgenomen Putboring Nabij, Amsterdam: (Notes on the Probable Results of Well Drilling near Amsterdam) // Tijdschrift van het Koninklijk Instituut van Ingenieurs. Hague, 1888, pp. 8–22. Banks, H.O., Richter, R.C., and Harder, J., Sea Water Intrusion in California, J. Amer. Water Works Assoc., 1957, vol. 49, no. 1, pp. 71-88. Bear, J., Hydraulics of Groundwater, New York: McGrow-Hill, 1979. Brown, R.H. and Parker, G.G., Salt Water Encroachment in Limestone at Silver Bluff Miami, Florida, Econ. Geol., 1945, vol. 40, no. 4, pp. 235–262. Chow, V.T., Handbook of Applied Hydrology, New York: McGrow-Hill, 1986. Collins, M.A. and Gelhar, L.W., Seawater Intrusion in Layered Aquifers, Water Resour. Res., 1971, vol. 7, no. 4, pp. 971–979. Contractor, D.N., Numerical Modeling of Saltwater Intrusion in the Nothern Guam Lens, Water Res. Bul., 1983, vol. 19, no. 5, pp. 745–751. Das, F. and Datta, B., Optimization Based Solution of Density Dependent Seawater Intrusion in Coastal Aquifers, J. Hydrol. Eng. ASCE, 2000, vol. 5, pp. 82–89. Dei, H., Seawater Intrusion into Unconfined Coastal Aquifer, J. of Jap. Ass. of Groundwater Hydrol., 1975, vol. 17, no. 1, pp. 30–46. Essaid, H.I., A Comparison of the Coupled Fresh Water– Salt Water Flow and Thee Ghyben-Herzberg Sharp Interface Approaches to Modeling of Transient Behavior in Coastal Aquifer Systems, J. of Hydrology, 1986, vol. 86, nos. 1-2, pp. 169–193. Frind, E.O., Simulation of Long-Term Transient Density-Dependent Transport in Groundwater, Adv. Water Resour., 1982, vol. 5, no. 2, pp. 73–88. Gao Xueping, Tu Xiangyang, Yu Lili. Research on the Groundwater Contamination Caused by Using Seawater as Landscape Water in Channel Coastal, Transactions of the CSAE, 2006, vol. 22, pp. 132–137. Groundwater Problems in Coastal Areas, Custodio, E., Ed., Paris: UNESCO, 1987. Henry, H.R., Salt Water Intrusion into Coastal Aquifers, Int. Assoc. Sci. Hydrol. Publ., 1960, vol. 52, pp. 478– 487. Herzberg, A., Die Wasserversorgung Einiger Nordsee Baden, J. Gasbeleucht. Verw. Beleuchtungsraten Wasservergsorg, 1901, vol. 44, no. 6, pp. 815–819. Isaacs, V.T. and Hunt, B., A Simple Approximation for a Moving Interface in a Coastal Aquifer, J. of Hydrology, 1986, vol. 83, nos. 1-2, pp. 29–43. Keulegan, H.G., An Example Report on Model Laws for Density Currents, Gaithersburg, U.S. National Bureau of Standards, 1954.

286

KHUBLARYAN et al.

37. Khublarian M.G. and Frolov A.P., Seawater Intrusion in Estuaries and Aquifers, Publication of the Water and Environment Research Institute, Helsinki, 1989. 38. Kohout, F.A., Flow Pattern of Fresh Water and Salt Water in the Biscayne Aquifer of the Maiami Area, Florida, Int. Assoc. Sci. Hydrol. Publ., 1960, vol. 52, pp. 440–448. 39. Krieger, R.A., Hatchett, J.L., and Poole, J.F., Preliminary Survey of the Saline-Water Resources of the United States, Water Supply Paper, 1957, no. 1374, pp. 172– 184. 40. Lee, C.-H. and Cheng, R.T., S. On Seawater Encroachment in Coastal Aquifers, Water Resour. Res., 1974, vol. 10, no. 5, pp. 1039–1043. 41. Lusczynski, N.J. and Swarzenski, W.V., Fresh and Salty Groundwater in Long Island, New York, J. of Hydraul. Div. Proc. ASCE, 1962, vol. 88, no. 4, pp. 173–194. 42. Meinardi, C.R., Fresh and Brackish Groundwater under Coastal Areas and Islands, Geojornal, 1983, vol. 7, no. 5, pp. 413–425. 43. Mercer, J.W., Larson, S.P., and Faust, C.R., Simulation of Salt Water Interface Motion, Ground Water, 1980, vol. 18, no. 4, pp. 374–385. 44. Pinder, G.F. and Cooper, H.H., A Numerical Technique for Calculating the Transient Position of the Saltwater Front, Water Resour. Res., 1970, vol. 6, no. 3, pp. 875– 882.

45. Reilly T.E., Goodman A.S. Quantitative Analysis of Saltwater–Freshwater Relationship in Groundwater Systems—a Hystorical Perspective, J. of Hydrology, 1985, vol. 80, no. 1-2, pp. 125–160. 46. Schmorak, S. and Mercado, A., Upconing of Fresh Water–Sea Water Interface below Pumping Wells, Field Study, Water Resour. Res., 1969, vol. 5, no. 6, pp. 1240– 1311. 47. Segol, G., Pinder, G.F., and Gray, W.G., A Galerkin Finite Element Technique for Calculating the Transient Position of the Saltwater Front, Water Resour. Res., 1975, vol. 11, no. 2, pp. 343–347. 48. Shamir, U. and Dagan, G., Motion of the Seawater Interface in Coastal Aquifers: Numerical Solution, Water Resour. Res., 1971, vol. 7, no. 3, pp. 644–657. 49. Tamai, N. and Shima, S., Saltwater Wedge in Unconfined Coastal Aquifers, Transaction of Japan Soc. Civil Eng., 1967, vol. 139, no. 1, pp. 31–38. 50. Todd, O.K. and Huisman, L., Groundwater Flow in the Netherlands Coastal Dunes, J. of Hydraul. Div. Proc. ASCE, 1959, vol. 85, no. 7, pp. 63-81. 51. Volker, R.E. and Rushton, K.R., An Assessment of the Importance of Some Parameters for Seawater Intrusion in Aquifers and a Comparison of Dispersive and Sharp Interface Modeling Approaches, J. of Hydrol., 1982, vol. 56, nos. 3−4, pp. 293–250.

WATER RESOURCES

Vol. 35

No. 3

2008