Dynamics of negative hydraulic barriers to prevent seawater intrusion

María Pool & Jesús Carrera

Abstract Negative hydraulic barriers that intercept Introduction

inflowing saltwater by pumping near the coast have been
proposed as a corrective measure for seawater intrusion in Seawater intrusion is a common problem in coastal
cases where low heads must be maintained. The main aquifers. Density differences between freshwater and
disadvantage of these barriers is that they pump a seawater together with intense pumping cause seawater
significant proportion of freshwater, leading to contami- to intrude inland along the aquifer bottom. The problem is
nation with saltwater at the well. To minimize such especially severe in arid and semi-arid zones where
mixing, a double pumping barrier system with two alternative water resources are scarce. In such cases,
extraction wells is proposed: an inland well to pump corrective measures are needed to restore groundwater
freshwater and a seawards well to pump saltwater. A quality and optimally manage the aquifer. However, it is
three-dimensional variable density flow model is used to not easy to select the optimal corrective measures to
study the dynamics of the system. The system performs control seawater intrusion because their feasibility and
very efficiently as a remediation option in the early stages. impact depend on the hydraulic and geometric character-
Long-term performance requires a well-balanced design. istics of the aquifer and on the management structure (Post
If the pumping rate is high, drawdowns cause saltwater to 2005). Consequently, mitigating seawater intrusion continues
flow along the aquifer bottom around the seawater well, to be a challenge.
contaminating the freshwater well. A low pumping rate at The most simple and cost-effective measures for
the seawards well leads to insufficient desalinization at the seawater intrusion are designed to improve the ground-
freshwater well. A critical pumping rate at the seawater water balance in an effort to keep heads high. This may be
well is defined as that which produces optimal desalin- achieved by reducing pumping rates or relocating pump-
ization at the freshwater well. Empirical expressions for ing fields to zones that are less susceptible to intrusion
the critical pumping rate and salt mass fraction are (Sherif 1999). This can also be achieved by artificial
proposed. Although pumping with partially penetrating recharge, adding water to an aquifer through man-made
wells improves efficiency, the critical pumping rates systems including infiltration basins, canals and recharge
remain unchanged. ponds (Li et al. 1987; Reichard 1995). One special case is
the aquifer storage and recovery (ASR) approach, which
Keywords Seawater intrusion . Double pumping barrier . consists of injecting freshwater into the aquifer to recover
Critical pumping rate . Groundwater management . essentially the same water without the need for additional
Salt-water/fresh-water relations treatment. This approach could be used to set up a
hydraulic saltwater intrusion barrier (Misut and Voss
2007). Limiting pumping to the shallowest portions of the
Received: 6 February 2009 / Accepted: 1 August 2009
Published online: 3 October 2009 aquifer is a good method to redress the hydraulic balance
and minimize drawdowns. Although these measures are the
* Springer-Verlag 2009 most simple, they are not always easy to implement
because (1) wells are already drilled and drilling new
M. Pool ())
(shallower) ones is expensive, (2) other water sources or
Department of Geotechnical Engineering and Geosciences, space are not available and (3) implementation of these
School of Civil Engineering, measures requires a strong management infrastructure.
Technical University of Catalonia, One attractive solution is to design seawater barriers
08034, Barcelona, Spain that prevent seawater from flowing inland, thereby
e-mail: maria.pol@idaea.csic.es
Tel.: +34-93-4017244 protecting groundwater pumping zones. Different types
Fax: +34-93-4017251 of barrier designs can be considered. Generally speaking,
“positive” indicates a source of fluid to the aquifer. This is
J. Carrera
Institute of Environmental Assessment and true conceptually (heads rise and water resources are
Water Research (IDAEA-CSIC), increased) and in models and equations (a negative flow
08028, Barcelona, Spain often represents an outflow). This terminology is therefore

Hydrogeology Journal (2010) 18: 95–105 DOI 10.1007/s10040-009-0516-1

96
proposed to classify the barriers into four types: low wells and to a reduction in freshwater resources. Never-
permeability subsurface, positive, negative and mixed. theless, negative barriers may be the only possible
corrective measures to prevent seawater intrusion in
– Low permeability subsurface barriers (Fig. 1a) consist of aquifers where the water level cannot be raised. They
vertical walls (slurry walls, steel or concrete sheet piles) may also be used to feed desalinization plants.
placed inland to block seawater intrusion, (Sugio et al. – Mixed barriers (Fig. 1d) are injection-pumping systems.
1987). This system requires considerable engineering Such a system may be used to inject freshwater at
and investment. Moreover, it may be counterproductive shallow depths while pumping saltwater at depth, which
if pumping stops or sources of contamination exist makes sense only in thick aquifers and/or in vertically
(Bolster et al. 2007). anisotropic systems. Corrective measures may also be
– Positive hydraulic barriers (Fig. 1b) inject water into used to pump at locations of the aquifer that are well
the aquifer, which raises the piezometric head, thus connected to the sea while injecting at locations that are
preventing saltwater from flowing inland (Bray and well connected to the land. Therefore, a well-designed
Yeh 2008). It is often argued that positive barriers are mixed barrier is much more sensitive to a good hydraulic
inefficient because they inject freshwater into the characterization than other barrier systems.
ocean. This is not necessarily the case. Abarca et al.
(2006) demonstrated that their efficiency is greater Specific methods have been developed for a range of
than that belief (i.e., the allowable increase in inland specific problems. In this regard, Van Dam (1999)
pumping exceeds the recharged flow rate) because proposed scavenger wells as a technique to control the
they not only increase available resources but also position of the interface and prevent upconing when
protect inland wells. The main disadvantages of this simultaneously pumping freshwater and seawater. Scav-
corrective measure are the need for high-quality water enger wells consist of wells near freshwater pumping
whenever injection is undertaken through wells and wells with screens open in the saline groundwater at some
the need for space if recharge through basins is depth below the freshwater and seawater interface. This
necessary. These corrective barriers also require consid- author evaluated the depth of the interface in terms of the
erable maintenance to control clogging. pumping rates. Scavenger wells require careful monitoring
– Negative hydraulic barriers (Fig. 1c) pump near the and are only effective in very thick aquifers.
shore, thus intercepting inflowing saltwater (Todd The study reported here was prompted by seawater
1980). Seawater intrusion is impeded as long as the intrusion in the Mar del Plata aquifer, Argentina (Bocanegra
barrier pumps. In practice, this system is often applied et al. 2001), where deep pumping is necessary for ground-
with minimum or no planning. It is common knowledge water supplies because of contaminated sources. Moreover,
that inland wells are protected by coastal wells. Inland water levels cannot be raised given the susceptibility of
well owners are prepared to go to great lengths to keep underground urban infrastructures to flooding. This problem
coastal wells pumping. However, these barriers often occurs in many cities (Vázquez-Suñé et al. 2005). Negative
end up pumping much more freshwater than saltwater, barriers are therefore the only solution to prevent seawater
leading to a mixing of freshwater with seawater at the intrusion.

a) Low permeability subsurface barrier b) Positive hydraulic barrier

Freshwater
injection well

Seawater Freshwater Seawater Freshwater

c) Negative hydraulic barrier d) Mixed barrier

Seawater Seawater Freshwater
pumping well pumping well injection well

Seawater Freshwater Seawater Freshwater

Fig. 1 Types of barriers to control seawater intrusion: a low permeability subsurface barrier, b positive hydraulic barrier, c negative
hydraulic barrier, d mixed barrier (injection-pumping system)

Hydrogeology Journal (2010) 18: 95–105 DOI 10.1007/s10040-009-0516-1

 97
There are few papers on negative barriers. Sherif and Problem description
Hamza (2001) demonstrated qualitatively that seawater A homogeneous confined coastal aquifer with constant
intrusion may be controlled by pumping brackish water thickness b is considered. The system consists of lines of
from the dispersion zone. These authors used an ideal fully penetrating freshwater-seawater wells along the coast
two-dimensional (2D) leaky aquifer to evaluate the effect separated by a distance of 2Lc (Fig. 3). Because of the
of pumping brackish water on the width of the mixing symmetry, only a coastal portion of length Lc is simulated
zone. However, the loss of freshwater in the well was (Fig. 3). A constant flow rate (Qp =648m3/d) of freshwater
considerable. (with salt mass fraction equal to ω=0 kg/kg) is imposed
The present paper addresses the loss of freshwater along the inland boundary. Pressure p is prescribed along
through mixing at the negative barrier. To minimize this the seaside boundary (p=ρs gz, where ρs is density of the
loss, a double pumping barrier system with two extraction seawater, g is the magnitude of the gravitational accel-
wells (Fig. 2) is proposed. One of the wells, the “freshwater eration vector and z is the elevation).
well”, is designed to pump freshwater while the other, the A non-dispersion boundary condition is adopted for trans-
“seawater well”, pumps mostly seawater. A vertical cross port along the coastal boundary. This implies that salt mass
section (see Fig. 2) suggests little flow between the two fraction equals either that of seawater (ωs =0.0386 kg/kg) for
wells. In fact, Kacimov et al. (2008) used vertical two- inflowing portions of the boundary or that of the resident mass
dimensional cross-section simulations to study the possi- fraction for outflowing portions of the boundary (Voss and
bility of restoring groundwater quality in coastal aquifers by Souza 1987; Frind 1982). The remaining boundaries are
pumping freshwater and seawater from the vicinity of the closed to flow and solute transport.
shoreline. These authors demonstrated that a zone of Flow and transport parameters used for the selected
reduced velocity was created between the pumping strips model are as follows: freshwater hydraulic conductivity
and that the intrusion was mainly restricted to the area k=10.6m/d, porosity f=0.25, dimensionless


density
 dif-
between the sea boundary and the saline groundwater ference ε=0.027 (ε is given by " ¼ s  f = f with ρf
pumping field protecting the freshwater pumping. Unfortu- being density of the freshwater), longitudinal dispersivity
nately, the vertical cross section approach is only valid when αL =10 m, transverse dispersivity αT =1 m, molecular
pumping wells are closely spaced with respect to their diffusion coefficient Dm =10−9 m2/s, freshwater viscosity
distance from the coast. Otherwise, the flow system is μ=0.001 kg/ms, aquifer thickness b=50 m, the width of
complex and fully three dimensional. model domain Lc =500 m and the freshwater well, located
The aim of the present work is to elucidate the dynamics at a distance of Lf =300 m from the sea , pumps Qf =
of double pumping barriers and to assess to what extent the
three dimensionality of the flux compromises the efficiency
of the system. Efficiency is assessed here by the pumping A′ B′
rate at the seawater well needed to reach a minimum Inland
salinity at the freshwater well. The evolution of seawater
intrusion in an ideal system is analyzed in order to identify
the optimal location and pumping rates in the system
applied, protect freshwater wells and control seawater
intrusion. Using this approach the problem is formulated
in a dimensionless form and empirical expressions for the Freshwater Well
controlling variables are obtained. Sea Seawater Well
A B
Lc
B B′
Dynamics of double pumping barriers A A′
A methodology based on numerical three-dimensional
B′
simulations is employed to evaluate the dynamics of A′
double pumping barriers.

(kg/kg) Freshwater Well

Seawater well Freshwater well
0.0386 Qf
b
Seawater Well
0.015 Qs
Δ z
h=0 y x
Lf Ls Lc
0.00
Seawater Freshwater
Fig. 3 Plan view of the double negative barrier system, vertical
Fig. 2 Schematic description of a double pumping barrier. Ideally, cross-sections and model domain. Note that only the AB–A’B’
groundwater flux is low between the two wells region needs to be modeled because of symmetry

Hydrogeology Journal (2010) 18: 95–105 DOI 10.1007/s10040-009-0516-1

98
400 m3/d (here and below pumping rates refer to what is bottom because freshwater floats on top of saltwater
pumped from the model domain, which is effectively half (the well is fully penetrating). When the seawater well is
of the total pumping rate). Modeling seawater intrusion activated, groundwater fluxes around the freshwater well
requires the flow and transport equations to be coupled drift seawards as does the seawater at the bottom of the
through water density, giving a set of two coupled non- freshwater well, which becomes significantly desalinized
linear equations. Computer simulations were performed (Fig. 5, stage 2). However, as shown in Fig. 5, stage 3, the
with SUTRA (Voss and Provost 2002). drawdown generated by both wells causes the toe to
penetrate further inland along the boundary (i.e. between
seawater wells, B–B’, as displayed in Fig. 3). As a result,
Modeling methodology and results seawater flows laterally towards the freshwater well along
The simulations are carried out following a two-step the bottom of the aquifer. It should be noted that this
sequence. The first step simulates a typical coastal well mechanism is a 3D effect, resulting from buoyancy forces,
salinization. A well pumps freshwater (Qf) at a distance Lf which would not show up in vertically integrated (2D
from the sea. Salinity of pumped water rises gradually to areal) fluxes. At any rate, the higher the pumping rate at
steady state I. The second step constitutes a double pumping the seawater well, the greater the desalinization in stage 2
barrier solution to the problem. Using the end of step 1 as and the greater the lateral flow from the sea to the
the initial condition, a seawater well starts pumping (Qs) at a freshwater well. The salt mass fraction eventually
distance Ls from the sea, until steady state II is reached. increases at the freshwater well when the pumping rate
The resulting evolution of salinity at the freshwater at the seawater well intensifies. Therefore, the critical
well (ωf/ωs) is displayed in Fig. 4. At first, the system pumping rate at the seawater well (Qc/Qf) is defined as
appears to be effective. Pumping the seawater well causes that which yields the lowest salt mass fraction at the
salinity to fall sharply at the freshwater well. The freshwater well (ωm/ωs).
efficiency of the system remains high for a short period The location of the seawater well considerably
(several years with the parameters adopted here), but the affects the critical pumping rate value (Fig. 6). As the
situation is transient. After a relatively short period distance of the seawater well from the sea decreases
(3–5 years), salinity at the freshwater well increases again (from 150 m to 20 m), both the value of the critical
to reach steady state II. Figure 4 shows the two most pumping rate and the desalinization level increase. Note
remarkable features of the double negative barrier: first, an that the landward relocation of the seawater well at
immediate fall in salinity after installation of the seawater some distance from the shoreline allows one to reduce
well is followed by recovery with the result that the long- the pumping rate at the seawater well despite diminish-
term efficiency of the system is poor; second, the long- ing desalinization efficiency.
term behavior is sensitive to the barrier flow rate. These
two features are discussed in the following.
The evolution of salinization is displayed in Fig. 5. In Sensitivity analysis
stage 1 the freshwater well becomes salinized in the
absence of the seawater well. Salinization occurs at the A sensitivity analysis is carried out to describe the
behavior of the critical pumping rate at the seawater well
and the minimum value of the salt mass fraction at the
Qs /Qf =1.0 freshwater well. The system is defined in terms of
0.08 1.5 dimensionless numbers.
2.0
Steady State I 6.0
0.06 Dimensionless form of the governing equations
s

The governing equations are recast in a dimensionless form

/

Steady State II to obtain the dimensionless numbers that describe the

f

0.04 behavior of the system. A Cartesian x,y,z coordinate system

is adopted, where the z axis points vertically upwards. The
0.02
dimensionless coordinates x′, y′ and z′ are defined in terms
of the width of model domain Lc (see Fig. 3).
x 0 y z
0 x0 ¼ ; y ¼ ; z0 ¼ ; ð1Þ
0 50 Time(y) 100 Lc Lc Lc
Freshwater Well Double Barrier The dimensionless form of the distances of the freshwater
Fig. 4 Temporal evolution of seawater fraction at the freshwater well (Lf) and the seawater well (Ls) from the sea, and the
well for different pumping rates at the seawater well located at 90 m thickness of the aquifer become
from the coastal boundary. The first steady state represents a typical
coastal well salinization. The seawater well starts pumping Qs after 0 Ls 0 Lf b
this steady state, which causes an immediate drop in salinity but Ls ¼ ; Lf ¼ ; b0 ¼ ð2Þ
eventually leads to steady-state II Lc Lc Lc

Hydrogeology Journal (2010) 18: 95–105 DOI 10.1007/s10040-009-0516-1

 99
SALINIZATION (Steady State I)
0.08
Q s /Q f =2.0 Top Bottom
6.0
0.06 (1) (1) (1)

s
0.04 (3)

/ f
0.02
(2)

0
0 50 Time(y) 100
Freshwater Well Double Barrier

DOUBLE PUMPING BARRIER

INITIAL FINAL (Steady State II)
Bottom Top Bottom
(2) (3) (3)
Qs /Q f = 2

(2) (3) (3)

Q s /Q f = 6

Fig. 5 Three-dimensional dynamics of the double negative barrier. Concentration maps (same color scale as Fig. 3) and velocity vectors at
the top and bottom of the aquifer for two different pumping rates at the seawater well. A fully penetrating well is salinized from the aquifer
bottom (stage 1, above). When the seawater well is activated, salinity falls sharply at the freshwater well, especially when the barrier
pumping rate is high (stage 2). However, drawdown generated by both wells causes seawater to flow towards the freshwater well at depth
(stage 3). Note that the seawater well still draws some freshwater from the top portion of the aquifer at this stage

Darcy’s velocity, density, salt mass fraction and freshwater

head are written in a dimensionless form with respect to the
constant freshwater flux from inland (qp), the seawater
(Q c /Q f , m/ s)
density (ρs), the seawater salt mass fraction (ωs) and a
characteristic equivalent freshwater head, respectively
0.06
q 0  w 0 hky
q0 ¼ ;  ¼ ; w0 ¼ ;h ¼ ð3Þ
qp s ws Lc qp
Ls / Lc = 0.30
s

Note that ω′, the dimensionless salt mass fraction, can be

/

viewed as the proportion of seawater. The flow dimension-

f

0.18
0.04
0.14 less numbers are defined as follows

0.08 qd Kx Kz
a¼ ; rk1 ¼ ; rk2 ¼ ð4Þ
0.04
"Kz Ky Ky
0.02

0 5 10 15 where ε is given by " ¼ s  f =f with ρf and ρs the
Q s /Q f freshwater and seawater densities, respectively.
Fig. 6 Salt mass fraction at the freshwater well during steady state Here, a is the ratio of the freshwater flux to the
II versus pumping rate at the seawater well for five distances of the characteristic buoyancy flux, and rk1 and rk2 are the
seawater well from the coastal boundary (Ls/Lc =0.3, 0.18, 0.14,
0.08 and 0.04 m). The critical pumping rate at the seawater well hydraulic conductivity anisotropy ratios.
(Qc/Qf) is the one that causes minimum salinity at the freshwater According to these definitions, the dimensionless form
well (ωm/ωs) of the fluid mass balance equation in steady state (e.g.,

Hydrogeology Journal (2010) 18: 95–105 DOI 10.1007/s10040-009-0516-1

100
Bear 1972) including the pumping rates of freshwater and where rα denotes the ratio between longitudinal and
seawater wells reads as transverse dispersivities (ra ¼ aL =aT ) and qi and qj are the
  components of the flux.
0 @ 2 h0 @ 2 h0 @ 2 h0 1 @w0 The transport equation is subject to a zero dimension-
  rk1 0 2 þ 0 2 þ rk2 0 2 þ  q0 r0 0 less solute flux along all the boundaries except for the
@x @y @z a @z0
coastal boundary where the dimensionless salt mass flux is
0
 0 0 0
¼ 0 *Qf d 0 x  xf ; 0 ff ðzÞ þ 0 *Qs d 0 ðx  xs ; 0Þfs ðzÞ defined as follows
(
ð5Þ  p0 q0  n if q0 n < 0
where ρ′* is the dimensionless density of the incoming water ð0 q0 w0  0 D0 r0 w0 Þy0 ¼0  n ¼ 
0 q0 w0 y0 ¼0  n if q0 n > 0
through the source term, ∇′ expresses the dimensionless form
of the del operator, the Dirac’s delta δ′ represents the ð11Þ
dimensionless
0
location
0
of the wells and the dimensionless
functions ff ðzÞ and fs ðzÞ express the vertical distribution of where n is the unit vector normal to the boundary and
extraction, assumed to be uniform here. They are zero, pointing outwards.
except in the portion of the aquifer thickness where the wells
are pumping. The dimensionless form of the pumping rates
0 0
Qs and Qs are defined as the ratio of the pumping rates of the Critical pumping rate
wells to the discharge flowing into the sea Different sets of simulations were carried out by varying
independently the distances of the two wells from the sea
0 Qf 0 Qs
Qf ¼ ; Qs ¼ ð6Þ (Lf =200–400 m and Ls =20–150 m), the pumping rate at
qp Lc b qp Lc b the freshwater well (Qf =100–600 m3/d), the aquifer
thickness (b=9–50 m), the hydraulic conductivity (kx =
The dimensionless form of the boundary conditions along 10–100 m/d and kz =1–10 m/d), the longitudinal and
inland and coastal boundaries becomes transverse dispersivities (αL =10–50 m and αT =1–10 m)
 and the width of the model domain (2Lc =500–1,000 m).
@h0  1 0
0  ¼ 1 h0 jx0 ¼0 ¼ z ð7Þ For each problem, the pumping rate at the seawater well
@y Ly a rk2 was modified (382 cases) until the critical pumping rate
y0 ¼
Lc
was found. Table (1) shows the values of the dimension-
As regards the transport equation, Peclet numbers were chosen less numbers used for simulations.
as dimensionless parameters. They describe the relative The algorithm of Furnival and Wilson (1974) was used
importance of advective and diffusive transport mechanisms. to seek a simple empirical expression for the critical
pumping rate from the dimensionless numbers. The
Dm f aT algorithm finds the best subset regressions for a regression
bm ¼ ; bT ¼ ð8Þ model with candidate independent variables. The routine
qd Lc Lc eliminates some subsets of candidate variables by obtain-
ing lower bounds on the error sum of squares from fitting
where f is porosity, Dm the molecular diffusion coefficient and larger models.
αT is the transverse dispersivity [L].
The dispersivity anisotropy ratio is also proposed as a
dimensionless parameter. Table 1 Description of the simulations considered in the sensitivity analysis
This leads to the dimensionless form of the salt mass
Dimensionless Range Description
conservation in steady state as variables
r0 ð0 w0 q0 Þ  r0 f0 ðbT D0 þ bm IÞr0 w0 g rk1 =kx/ky 1.0–10 Hydraulic conductivity
rk2 =kz/ky 0.1–1.0 anisotropy ratios

 0 rL =Ls/Lf 0.06–0.75 Geometric shape factor
¼ þ0 *w0 *Q0f d 0 x  xf ; 0 ff ðzÞ b′=b/Lc 0.018–0.2 Aquifer thickness
0 0
a=qd/εkz 0.09–0.9 Characteristic
þ 0 *w0 *Qs d 0 ðx  xs ; 0Þfs ðzÞ ð9Þ 0
buoyancy flux
Q0f ¼ Qf =Lc qd b 0.6–0.9 Pumping rates
Qs ¼ Qs =Lc qd b 0.3–13
where 5′* is the dimensionless salt mass fraction of the bm ¼ ðDm fÞ=qd Lc 1.6∙10−6– Diffusive Peclet
incoming water through source term, I is the identity matrix 3.3∙10−6 number
and D′ is the dimensionless dispersion tensor, defined by bT =αT/Ls 0.013–0.1 Dispersive Peclet
" # number
0 0
qi qj ra ¼ aL =aT 1.0–25 Dispersivity
0 Dij 0  
D ¼ ¼ d ij jq j þ ðra  1Þ 0 ; i; j ¼ x; y; z 0
anisotropy ratio
aT qd jq j LAb ¼ L1f  L1c LGH −0.25–0.73 Toe penetration for
dispersive Henry
ð10Þ  FLDS
Lf problem

Hydrogeology Journal (2010) 18: 95–105 DOI 10.1007/s10040-009-0516-1

 101
Basically, the dimensionless parameters defined in the Note that the exponent in Eq. (15) depends on the
previous section were proposed as independent variables. dispersive case correction to the Ghyben-Herzberg toe
Some derived dimensionless numbers were also employed. approximation (Eq. 12) and on the ratio of the freshwater
In addition to these, the dimensionless toe penetration for the pumping to the freshwater inflow from inland, i.e. the
dispersive Henry problem (Abarca et al. 2007) was also exponent expresses the seawater invasion distance as
used. Abarca’s correction to the Ghyben-Herzberg toe evidenced by the reduction in freshwater discharge to the
penetration was modified to read sea because of pumping.
The empirical expression approximates the 3D dynam-
  ics of the system satisfactorily. However, it is only aimed
0 1 1 FLDS
LAb ¼  LGH  ð12Þ at preliminary assessments. The empirical approximation
Lf Lc Lf for the critical pumping rate is valid for the range of
values of the variables tested here. An exhaustive study is
where required to fine tune this kind of corrective measure to the
individual aquifer.
 pffiffiffiffiffiffiffiffiffiffiffi 0:724  pffiffiffiffiffiffiffiffiffiffiffi 0:362 Equation (15) facilitates the sensitivity analysis
aT aL a T aL (Fig. 8). Perhaps, the most surprising feature is the
FLDS ¼ 0:136 2 pffiffiffiffiffiffi þ 0:69 2 pffiffiffiffiffiffi 0
reduction in critical pumping rate for increasing Qf . The
a rk2 a rk2
reason for this is that the toe along the boundary BB’
ð13Þ 0
(Fig. 3) proceeds inland when Qf increases. Thus an
and LGH is the toe penetration without pumping under the increase in freshwater pumping causes the critical pump-
Ghyben-Herzberg approximation ing rate to decrease, leading to a reduction in the long term
lateral flux of seawater. The effect is more marked when
b " Ky the seawater well is close to the sea (Fig. 8a).
LGH ¼ ð14Þ When transverse dispersivity is increased with
2qd respect to the longitudinal one, the interface retreats
seawards and the salt mass fraction at the freshwater
An empirical expression for the critical pumping rate Fc was well is reduced. However, the efficiency of mixing
obtained from the regression model. The equation that best increases and the mixing zone is broadened. This
represents the numerical results is (see Fig. 7) promotes salinization of the freshwater well when the
pumping rate in the seawater well intensifies. Conse-
! quently, increasing transverse dispersivity decreases the
Qc 1 pffiffi 0 0
 Fc  1:6 e½ pðLAb Qf Þ ð15Þ critical pumping rate.
Qf rL ðrk1 Þ2=3 Buoyancy effects increase in importance with aquifer
thickness. Thus, an increase in aquifer thickness would
tend to favor lateral fluxes around the seawater well.
where rL is defined as the geometric shape factor of the two Intercepting these fluxes requires increasing the pump-
wells (rL =Ls/Lf). ing rate at the seawater well. Therefore, the thicker the
aquifer, the higher the critical pumping rate. This effect
is more significant when the seawater well is near the
15
sea (Fig. 8b).
y=x An increase in kx reduces lateral gradients. In fact,
R2 = 0.997 it is similar to a reduction in the width of the model
domain, which implies a reduction in lateral fluxes. As
a result, the critical pumping rate decreases (Fig. 8c).
10 The distance of the two wells from the sea (Lf and Ls)
and the width of the model domain (Lc) also affect the
critical pumping rate. When the distance of the
Qc /Qf

seawater well from the sea is increased, drawdown

caused by the seawater well also increases (pressure
5
along the coastal boundary is prescribed). This would
tend to favor the lateral flux towards the pumping well.
Thus, the critical pumping rate needs to be reduced to
minimize this lateral flux (Fig. 8a). When the distance
0 of the two wells from the sea exceeds that of Lc
0
Fc
5 10 15
(width of the model domain, see Fig. 3), the flow field
Fig. 7 Critical pumping rate in the seawater well obtained
becomes increasingly 2D, which reduces the critical
from the numerical simulations (Qc/Qf) versus the results from pumping rate (Fig. 8d). Figure 8 also illustrates the fit
empirical expression, Fc (Eq. 15), obtained from the regression between the numerical and empirical results (Eq. 15),
model which is fairly good.

Hydrogeology Journal (2010) 18: 95–105 DOI 10.1007/s10040-009-0516-1

102

Minimum salt mass fraction in the freshwater well freshwater. Nevertheless, Eq. (16) was used as an
The dimensionless salt mass fraction (ωf/ωs) is the independent variable in the regression model to obtain
fraction of seawater at the freshwater well. If the the empirical expression (Gc) for the minimum seawater
seawater well pumped only seawater and if dispersion mass fraction at the freshwater well (ωm/ωs), i.e. the salt
mechanisms and buoyancy forces were neglected, the mass fraction at the freshwater well when the seawater
seawater fraction at the freshwater well would become well is pumping the critical rate. The best regression
the seawater inflow not captured by the seawater well obtained was (see Fig. 9).
divided by Qf, i.e.,  
wm QTs  Qs
wf QTs Qs  Gc ¼ 0:2
¼ ð16Þ ws Qf
ws Qf !
 
bT pffiffiffiffiffi b0 rL
where QTs is the total saltwater flow rate entering from þ 0:8 0:175 0 ra þ 0 þ 0:04 0
Ls Lf Qf
the seaside boundary. Using the method of images for
an infinite line of freshwater and seawater wells, the ð18Þ
total saltwater flow rate can be approximated as
The fit between numerical and empirical results is not as
ZLb good as the one for the critical pumping rate but may be
Q f Lf X
n
1 used to obtain an idea of the efficiency of desalinization.
QTs ¼ ðqp b  ð17Þ
p i¼1 ðx  2nLc Þ2 þ L2f Desalinization efficiency (reduction in salinity at the
0 freshwater well between steady state I and steady state II)
Q s Ls X n
1 is displayed versus distance of the seawater well from the
 Þdx
p i¼1 ðx  2nLc Þ2 þ L2B sea, and versus aquifer thickness in Figs. 10 and 11,
respectively. Decreasing the distance of the seawater well
from the sea leads to a reduction in the salt mass fraction
where Lb is the fraction of the coastal boundary where at the freshwater well. However, the critical pumping rate
seawater flows inland (0≤Lb ≤Lc). at the seawater well shows a marked increase (Fig. 6). For
Unfortunately, this approach neglects the most notable example, the desalinization efficiency increases from 19 to
features of seawater intrusion, namely buoyancy and 28% when the seawater well is brought closer to the sea
interface mixing. In reality, the seawater well pumps some (Ls =40 m), but the critical pumping rate at the seawater

(a) (b)
Empirical Empirical
8 b′ = 0.1, L ′f = 0.6 Q ′f = 0.9, L ′f = 0.6

6
Q f′ =0 .6 L′s = 0.08
Qc /Qf

4
L′s = 0.14

2
L′s = 0.30
Q f′ =0 .9
0
0 0.1 0.2 0.3 0 0.05 0.1 0.15
L ′s b′
(c) (d)
Empirical Empirical
8 Q′f = 0.6, b′ = 0.1 Q′f = 0.6, b′ = 0.1

rL = 0.23
6
Qc /Qf

L′s = 0.08
4
rL = 0.30
2 L′s = 0.18

rL = 0.50
0
0 2 4 6 8 10 0.5 0.75 1 1.25
rk1 L′f
Fig. 8 Sensitivity analysis. The critical pumping rate decreases when a distance from the coast increases, b aquifer thickness decreases, c
anisotropy increases, or d distance of the freshwater well from the coast increases. Note that the proposed empirical expression, Fc Eq. (15),
is fairly accurate for the range of values analyzed here

Hydrogeology Journal (2010) 18: 95–105 DOI 10.1007/s10040-009-0516-1

 103
0.15
y=x Empirical
R2 = 0.982 0.15 initial
25%
0.12 L′s = 0.3
0.12 L′s = 0.18
L′s = 0.14

m/ s
0.09 0.09
m/ s

58%
0.06 0.06
85%
0.03
0.03
drinking water standard
0
0 0.02 0.04 0.06 0.08 0.1 0.12
0
0 0.03 0.06 0.09 0.12 0.15
b′
Gc
Fig. 11 Fraction of seawater and desalinization percentage (%) at
Fig. 9 Fraction of seawater at the freshwater well from the the freshwater well versus aquifer thickness for . The results from
numerical results versus the empirical expression, Gc, Eq. (18), the empirical expression provide a satisfactory representation of the
obtained from the regression model numerical results. Note that the double pumping barrier is more
efficient for thin aquifers, and that drinking water is obtained from
well is four times higher when the seawater well is near the freshwater well with the double pumping system
the freshwater well (Ls =150 m). Desalinization efficiency
is significantly better in the case of a low freshwater is obtained from the freshwater well, and the pumping rate
pumping rate. at the seawater well is low (Qc/Qf ≈1), see Fig. 11.
One critical factor in the efficiency of the double
negative barrier is aquifer thickness. A decrease in aquifer
thickness leads to a high desalinization efficiency with the Partially penetrating wells
double pumping system (Fig. 11). Note that drinking As discussed in the preceding, the wells were assumed to
water is pumped from the freshwater well with the double be fully penetrating. However, salinity is stratified in
negative barrier when the dimensionless aquifer thickness coastal aquifers, with freshwater floating on top of
is 0.018. seawater. Therefore, one would expect the location and
A reduction in the aquifer thickness from 50 to 9 m, length of the well screen to be a critical factor in a double
while
0
increasing the pumping rate at the freshwater well pumping barrier system. The effect of partially penetrating
(Qf ¼ 0:9), gives rise to a seawater fraction of 7% (i.e., wells on the efficiency of double pumping barriers is
some 1,400 mg/L of chloride for a 1,900 mg/L seawater) considered.
at the freshwater well in steady state I. Pumping the The simulations with partially penetrating wells were
seawater well causes chloride to fall below 200 mg/L at carried out with the same flow and transport parameters as
the freshwater well (<1 % of seawater), i.e. drinking water those of aforementioned models, where hydraulic con-
ductivity was assumed to be isotropic. First, the base case
0.15 was simulated with a fully penetrating freshwater well
and a partially penetrating seawater well. The numerical
28% 19%
Q′f = 0.9
0.12
(Qc Qf, m s)
Empirical
0.06
m/ s

0.09 initial

0.06 FPw FPw

s

67% 50%
/

0.04 L′s = 0.30 L′s = 0.18

f

Q′f = 0.6
0.03 PPw PPw

0
0 0.1 0.2 0.3 0.4 0.02
0 5 10
L′s Qs / Qf
Fig. 10 Comparison between numerical and empirical results of Fig. 12 Fraction of seawater in the freshwater well versus
seawater fraction at the freshwater well versus the distance of the pumping rates at the seawater well at distances from the coast of
seawater well from the sea. Reduction in the distance of the seawater 150 and 90 m ( and 0.18 respectively ). The dashed lines denote the
well from the sea improves the desalinization efficiency of the barrier, results with fully penetrating wells (FPw) and the solid lines
as measured by the percentage (%) reduction in freshwater well salinity represent the results with partially penetrating wells (PPw). Note
with respect to steady-state I. Note, however, the marked increase in that partial penetration significantly improves desalinization effi-
the critical pumping rate (see Fig. 8c) ciency but does not change the critical flow rate

Hydrogeology Journal (2010) 18: 95–105 DOI 10.1007/s10040-009-0516-1

104
results confirm a decrease in salinity at the freshwater top of the aquifer and the seawater well pumps at the
well. This effect would be more evident if both wells were bottom of the aquifer. Moreover, the numerical results
partially penetrating wells. Different sets of simulations suggest that the critical pumping rate does not depend on
were carried out to study the effect of partially penetrating the location of the screen well.
wells. The freshwater well pumps in the top 15 m of the In summary, the double pumping barrier system is not
aquifer whereas the seawater well pumps in the bottom as efficient as it looks from a vertical cross section of the
15 m. The numerical results provide evidence that buoy- aquifer. Nevertheless, in cases where other options are not
ancy allows the freshwater well to pump freshwater, while viable because of the lack of space or the need to keep
the seawater well pumps seawater. Hence desalinization water levels low (which is the case of urban aquifers), this
efficiency is significantly enhanced with respect to the corrective measure is a viable alternative. In the case of
fully penetrating wells, Fig. 12. thin aquifers and/or closely spaced wells, the double
It should be noted that the critical pumping rate does not pumping barrier is fairly efficient. The system proposed is
depend on the location of well screen, i.e. whether the wells much more efficient than a simple negative barrier.
are fully or partially penetrating. In an anisotropic system,
where the vertical hydraulic conductivity is lower than the Acknowledgements The authors acknowledge the financial support
horizontal hydraulic conductivity, the salinity at the fresh- provided by the European Commission (ATRAPO project, contract
water well must be lower than that in the isotropic system. MEC CTM2007-66724-C02-01/TECNO) and the Spanish CICYT
However, it is very probable that the critical pumping rate (HEROS project). The first author gratefully acknowledges the
receipt of an FI award from the Autonomous Government of
occurred at exactly the same discharge given that the critical Catalonia for the period during which this work was carried out.
pumping rate depends on lateral fluxes, i.e. it depends on the The authors wish to thank Mark Bakker and Eric Reichard for their
horizontal components of the hydraulic conductivity. constructive comments on the manuscript.

Conclusions References

A double pumping barrier is proposed as a corrective Abarca E, Vázquez-Suñé E, Carrera J, Capino B, Gámez D, Batlle F
(2006) Optimal design of measures to correct seawater
measure to control seawater intrusion in coastal aquifers intrusion. Water Resour Res 42(9), W09415
where the water level cannot be raised or where artificial Abarca E, Carrera J, Sánchez-Vila X, Dentz M (2007) Anisotropic
recharge or restrictions on pumping are not feasible. The dispersive Henry problem. Adv Water Resour 30(4):913–926
system appears to be effective in the early stages (3– Bear J (1972) Dynamics of fluids in porous media. Elsevier,
Amsterdam, 764 pp
5 years). Pumping the seawater well causes salinity to fall Bocanegra E, Massone H, Martinez D, Civit E, Farenga M (2001)
sharply at the freshwater well. However, after a relatively Groundwater contamination: risk management and assessment for
short period, the drawdown generated by both wells landfills in Mar del Plata, Argentina. Environ Geol 40:732–741
causes saltwater to flow around the seawater well along Bolster D, Tartakovsky D, Dentz M (2007) Analitycal models of
the bottom of the aquifer, contaminating the freshwater contaminant transport in coastal aquifers. Adv Water Resour 30
(9):1962–1972
well. When the seawater well pumping rate is high, this Bray B, Yeh W (2008) Improving seawater barrier operation with
lateral flow is dominant. However, a low pumping rate at simulation optimization in southern California. J Water Resour
the seawater well would fail to intercept all seawater Plann Manage 134:171–180
intrusion. Therefore, there is a critical pumping rate at the Dam JC (1999) Seawater intrusion in coastal aquifers: concepts,
methods and practices. Chap. Exploitation, restoration and
seawater well which produces the minimum salt mass management. Kluwer, Norwell, MA, pp 73–125
fraction at the freshwater well. Paradoxically, a too Frind E (1982) Simulation of long-term transient density-dependent
energetic pumping at the seawater well would consid- transport in groundwater. Adv Water Resour 5:73–88
erably increase salinity at the freshwater well. Furnival G, Wilson R Jr (1974) Regression by leaps and bounds.
Technometrics 16:499–512
The location of the seawater well and aquifer thickness Kacimov A, Sherif M, Perret JS, Al-Mushikhi A (2008) Control of
significantly affect the efficiency of the system. Efficiency sea-water intrusion by salt-water pumping: Coast of Oman.
increases as the seawater well is brought closer to the sea but Hydrogeol J 17(3):541–558
at the expense of a higher pumping rate at the seawater well. Li C, Bahr J, Reichard E, Butler J, Remson I (1987) Optimal siting
A decrease in aquifer thickness leads to a high desalinization of artificial recharge: an analysis of objective functions. Ground
Water 25:141–150
at the freshwater well. Therefore, the highest efficiency of a Misut PE, Voss CI (2007) Freshwater-saltwater transition zone move-
double pumping barrier is reached in thin aquifers. ment during aquifer storage and recovery cycles in Brooklyn and
Empirical expressions were derived for the critical Queens, New York City, USA. J Hydrol 337(1–2):87–103
pumping rate and for the corresponding (maximum Post V (2005) Fresh and saline groundwater interaction in coastal
aquifers: is our technology ready for the problems ahead?
achievable) desalinization at the freshwater well. These Hydrogeol J 13(1):120–123
expressions capture the 3D dynamics of the system Reichard EG (1995) Groundwater-surface water management with
satisfactorily. However, they should only be used for stochastic surface water supplies: a simulation-optimization
preliminary assessments. An exhaustive study is required approach. Water Resour Res 31(11):2845–2865
Sherif M (1999) Seawater Intrusion in Coastal aquifers: concepts,
to tailor this type of corrective measure to the individual methods and practices, vol. 14. Kluwer, Norwell, MA, pp 559–590
aquifer. Finally, the efficiency of double-negative barriers Sherif MM, Hamza KI (2001) Mitigation of seawater intrusion by
is enhanced if the freshwater well is screened only at the pumping brackish water. Transp Porous Media 43:29–44

Hydrogeology Journal (2010) 18: 95–105 DOI 10.1007/s10040-009-0516-1

 105
Sugio S, Nakada K, Urish DW (1987) Subsurface seawater intrusion Voss CI, Provost A (2002) SUTRA, a model for saturated-
barriers analysis. J Hydraul Eng 113:767–779 unsaturated variable-density ground-water flow with solute or
Todd D (1980) Groundwater hydrology, chap. 14. Wiley, Chichester, UK energy transport. US Geol Surv Sci Invest Rep 02-4231
Vázquez-Suñé E et al (2005) Groundwater modelling as a tool for Voss CI, Souza WR (1987) Variable density flow and transport
the European Water Framework Directive (WFD) application: simulation of regional aquifers containing a narrow freshwater-
the Llobregat Case. Phys Chem Earth 31:1015–1029 saltwater transition zone, Water Resour Res 26:2097–2106

Hydrogeology Journal (2010) 18: 95–105 DOI 10.1007/s10040-009-0516-1