Hydrogeology Journal DOI 10.1007/s10040-014-1214-1

A simulation/optimization study to assess seawater intrusion management strategies for the Gaza Strip coastal aquifer (Palestine) Marta Dentoni & Roberto Deidda & Claudio Paniconi & Khalid Qahman & Giuditta Lecca Abstract Seawater intrusion is one of the major threats to freshwater resources in coastal areas, often exacerbated by groundwater overexploitation. Mitigation measures are needed to properly manage aquifers, and to restore groundwater quality. This study integrates three computational tools into a unified framework to investigate seawater intrusion in coastal areas and to assess strategies for managing groundwater resources under natural and human-induced stresses. The three components are a three-dimensional hydrogeological model for densitydependent variably saturated flow and miscible salt transport, an automatic calibration procedure that uses state variable outputs from the model to estimate selected model parameters, and an optimization module that couples a genetic algorithm with the simulation model. The computational system is used to rank alternative strategies for mitigation of seawater intrusion, taking into Received: 13 December 2013 / Accepted: 9 November 2014 * Springer-Verlag Berlin Heidelberg 2014 M. Dentoni ()) : R. Deidda University of Cagliari, Cagliari, Italy e-mail: [email protected] Tel.: +39 070 5300 R. Deidda e-mail: [email protected] C. Paniconi Institut National de la Recherche Scientifique, Centre Eau Terre Environnement (INRS-ETE), Quebec City, Canada C. Paniconi e-mail: [email protected] M. Dentoni : G. Lecca Center for Advanced Studies, Research and Development in Sardina (CRS4), Loc. Piscina Manna, Building 1, 09010, Pula, Cagliari, Italy G. Lecca e-mail: [email protected] K. Qahman Environment Quality Authority, Al Bireh, Palestine K. Qahman e-mail: [email protected]

account conflicting objectives and problem constraints. It is applied to the Gaza Strip (Palestine) coastal aquifer to identify a feasible groundwater management strategy for the period 2011–2020. The optimized solution is able to: (1) keep overall future abstraction from municipal groundwater wells close to the user-defined maximum level, (2) increase the average groundwater heads, and (3) lower both the total mass of salt extracted and the extent of the areas affected by seawater intrusion. Keywords Coastal aquifers . Palestine . Numerical modeling . Automatic calibration . Simulation/ optimization (S/O) technique

Introduction The delicate equilibrium of the interface between freshwater and seawater in coastal aquifers can be perturbed by natural and anthropogenic forcing, such as sea level rise, precipitation regimes changes, and groundwater exploitation. This can lead to inland encroachment and/or vertical upconing of seawater, causing contamination of freshwater resources. Seawater intrusion (SWI) is a major problem in many coastal aquifers worldwide (Werner and Gallagher 2006; Sanford and Pope 2010), and can be exacerbated if water resources are not properly managed. In Europe, the Mediterranean basin has been recognized as a particularly sensitive area (Antonellini et al. 2008; Custodio 2010; Zghibi et al. 2011). A good knowledge of the aquifer system and the ability to forecast its future behavior under different scenarios of natural forcing and human interactions are necessary conditions to properly manage SWI problems (Werner et al. 2012). Numerical simulation is a commonly used tool to investigate SWI in coastal aquifers and to assess alternative groundwater resource-management scenarios (Bear et al. 1999). Amongst available modeling approaches, SWI can be represented as a variable density flow and miscible transport process in a three-dimensional (3D) domain solved using numerically efficient computational codes. Notwithstanding efficient codes, a 3D SWI model is not simple to apply, due to the complexity of the involved physical processes and the basin geometry, the scarcity of adequate hydrogeological field data, and the heterogeneity

of hydrogeological characteristics (Carrera et al. 2009). The numerical simulation itself cannot guarantee a successful modeling if the calibration procedure is not properly carried out (Carrera et al. 2005). Calibration should thus be considered an integral part of the process of modeling and understanding a hydrogeological system (Poeter and Hill 1997). Calibration techniques can be classified as local or global. Local methods follow deterministic rules, searching the single parameter set by using information on the objective function (direct methods) or by including derivative information (gradient-based methods; Carrera and Neuman 1986). Examples of global techniques include genetic algorithms (GA), shuffled complex evolution (SCE) schemes (Duan et al. 1992), and simulated annealing (SA) search strategies (Sumner et al. 1997). A calibrated simulation model can help us understand the behaviour of the modeled system and can be applied for different purposes, including the evaluation of aquifer response to different groundwater management strategies. These kinds of evaluations are fundamental for coastal aquifers where the extent of SWI is particularly critical (Zhou et al. 2003). Several authors (e.g., Cheng et al. 2000; Shamir et al. 1984; Willis and Finney 1988) have shown how an optimization approach can be useful in the solution of SWI management problems. A more efficient approach can be achieved by linking groundwater simulation models with optimization techniques in a single framework (Dhar and Datta 2009). For given groundwater abstraction needs, a simulation/optimization (S/O) coupled system can determine the optimal management strategy, under user-defined constraints such as minimum groundwater head levels or maximum groundwater salt concentrations permitted. The coupling of simulation and optimization models can be achieved with different techniques: by binding the discretized governing equations into the constraints of the optimization model (embedding technique, Das and Datta 1999), by using a response matrix (Ahlfeld and Mulligan 2000; Gorelick 1983; Ndambuki et al. 2000; Yang et al. 2001), or by building an external linkage of simulation and optimization model (Das and Datta 2001). Embedding and response matrix approaches are not very feasible for highly complex and nonlinear 3D numerical models of SWI processes (Werner et al. 2012). The third alternative, externally linking the SWI code with a general-purpose optimization-based management tool, is more feasible for this class of simulation model (Barlow 2005; Bhattacharjya and Datta 2005; Dhar and Datta 2009). As with calibration methods, optimization techniques can also be global, and these include genetic algorithms (GA; Carroll 1996), which are a family of combinatorial methods that search for solutions of complex problems using an analogy between optimization and natural selection (Goldberg 1989). The selection is made according to solution fitness (objective function): the more suitable a solution is, the more chances it has to reproduce, and the population thus evolves toward the optimal solution. Qahman et al. (2009) linked Carroll’s Hydrogeology Journal

FORTRAN GA Driver (Carroll 1996) and the densitydependent variably saturated groundwater flow and salt transport model CODESA-3D (Gambolati et al. 1999; Lecca 2000) and applied the resulting S/O framework to the Gaza Strip coastal aquifer in Palestine. The decision variables, tested for optimality, were pumping rates and the GA objective function and constraints were based on potential head and salt concentrations at the wells, respectively. The parameters of the CODESA-3D model were set up according to available literature, and steady-state conditions were examined for a small square area (4 km2) of the Gaza Strip aquifer system subjected to SWI. In the present study, a computational framework for assessing management measures to restore groundwater quality in coastal aquifers is presented. The framework consists of the 3D groundwater model CODESA-3D, the automatic calibration procedure PEST (Doherty 2002), and the GA optimization driver proposed by Carroll (1996). The developed framework is then applied to the 365 km2 Gaza coastal aquifer (GCA) to: (1) reconstruct how aquifer pumping and other practices have impacted groundwater and salinity levels over the past decades; and (2) support the planning and assessment of future management strategies for the GCA. This report builds on previous SWI modeling studies in the Gaza region (e.g., Moe et al. 2001; Qahman and Larabi 2006; Qahman et al. 2009) by: (1) focusing specifically on the GCA rather than a much smaller (field) unit; (2) modeling for the first time the full complexity of the 3D geometry; and (3) performing an integrated calibration, validation, and optimization analysis for a period spanning 1935–2020.

Computational framework The developed framework integrates different computational tools to encompass all the fundamental steps needed to study coastal aquifer management options. The main components are described in the following paragraphs.

CODESA-3D Model The COupled DEnsity-dependent variably SAturated groundwater flow and miscible salt transport 3D model (CODESA-3D) is a distributed, fully 3D, variably saturated flow and miscible transport finite element model, accounting for spatial and temporal variability of model parameters and boundary conditions (Gambolati et al. 1999; Lecca 2000; Putti and Paniconi 1995). The flow and solute transport processes are coupled through the variable density of the filtrating mixture made up of water and dissolved matter (salt, pollutants). The flow module considers the case of variably saturated porous medium, applicable both to the unsaturated (soil) and the saturated (groundwater) zone, while the transport module assumes complete mixing between freshwater and saltwater bodies, giving rise to a variably dense filtrating fluid with a nonreacting solute (salt). CODESA-3D solves the system in DOI 10.1007/s10040-014-1214-1

terms of pressure heads and concentrations and, from these state variables, it is straightforward to derive watertable levels, groundwater velocities, and the extent of the saltwater–freshwater mixing zone. The CODESA-3D mathematical model is expressed in terms of two state variables. The first one is the equivalent freshwater pressure head ψðx; y; z; t Þ ¼ ρp0 g [L], where p [M L−1 T−2] is the pressure on point (x,y,z) [L] at time t [T], ρ0 [M L−3] is the freshwater density, and g [LT−2] is the gravitational constant. A derived variable is the equivalent freshwater hydraulic head h=ψ+z [L], where z [L] is the vertical coordinate directed upward. The second unknown is the normalized concentration of salt cðx; y; z; t Þ ¼ ec [/], defined as the ratio between actual (ec ) and ecmax maximum (ecmax ) absolute concentration of salt [M L−3] in the water solution. The value of ecmax for saltwater intrusion problems corresponds to the average salt concentration of seawater, which is usually in the range of 25–40 g/l (in terms of total dissolved salts). In CODESA-3D the variable density ρ is expressed by the linear function ρ=ρ0(1+εc) where ε=(ρS−ρ0)/ρ0 is the density difference ratio and ρS is the solution density at the maximum normalized concentration c=1. The coupled system of variably saturated groundwater flow and miscible salt transport equations is: 8 ∂ψ ∂c ρ > ¼ −∇⋅v−ϕS w ε þ q <σ ∂t ∂t ρ0 > : ϕ ∂ðS w cÞ ¼ −∇⋅ðcvÞ þ ∇⋅ðD∇cÞ þ qc þ f ∂t

ð1Þ

where σ(ψ,c) is the general storage term [L−1], ∇ is the gradient operator, v is the Darcy velocity vector [LT−1], ϕ is porosity [/], Sw(ψ) is the water saturation [/], q is the injected (positive)/extracted (negative) volumetric flow rate [T−1], D is the dispersion tensor [T−1], c* is the normalized concentration of salt in the injected/extracted fluid [/], and f is the volumetric rate of injected/extracted solute that does not affect the velocity field [T−1]. The general storage term is equal to (1+ε⋅c)(SwSs+ϕdSw/dψ), with Sw=ρ0g(α+ϕβ) [L−1] the specific elastic storage at the reference density, and α, β the medium and fluid compressibility, respectively. The Darcy velocity vector is expressed as v=−Κ[∇ψ+(1+ε⋅c)∇z] [L T−1] where Κ=Κs(1+ε⋅c)kr [L T−1] is the hydraulic conductivity tensor, Ks [L T−1] is the saturated hydraulic conductivity tensor at the reference density, and kr(ψ) [/] is the relative conductivity. The dispersion tensor [L2 T−1] is equal to e ¼ αT jvjδi j þ ðαL −αT Þvi v j =jvj þ ϕS w D0 τδi j , w i t h ϕS w D e is defined as in Bear and Verruijt (1987) with i,j=x,y,z. D αL and αT the longitudinal and transverse dispersivity qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi coefficients [L], respectively, jvj ¼ v2x þ v2y þ v2z , δij the Kronecker delta, D0 [L2 T−1] the molecular diffusion coefficient, and τ (=1) the tortuosity. To complete the mathematical formulation of the flow and transport problem, initial and Dirichlet, Neumann, or Cauchy boundary conditions are added. Hydrogeology Journal

For a comprehensive mathematical and numerical description of the CODESA-3D model, the reader is referred to Gambolati et al. (1999). The model has been applied to several coastal aquifers in the Mediterranean basin affected by SWI problems (Cau et al. 2002; Kerrou et al. 2007; Lecca et al. 2001; Paniconi et al. 2001; Qahman et al. 2009). CODESA-3D is the computational engine of the grid-enabled Web demonstrative hydrology application AQUAGRID (GRIDA3 2011; Lecca et al. 2009, 2012).

Automatic calibration with PEST From available automatic calibration models, PEST (Doherty 2002) was adopted, which can be used to estimate parameters for a given simulation model working with ASCII input/output files. In the coupling of PEST with CODESA-3D (Lecca and Cau 2006), PEST iteratively reads state variables (e.g., groundwater heads, h) from the output files produced by CODESA-3D and perturbs the parameter file to find optimal model parameters (e.g., hydraulic conductivities, K; see section ‘Model results’) by minimizing the objective function Φ: Φ¼

NC X

 2 wi 2 hi −hi

ð2Þ

i¼1

where Φ is the sum of squared weighted residuals between simulated (h i ) and field measured ( hi ) groundwater heads at i-th control point, i.e., water observation wells, where i=1,..,NC. Groundwater heads are easily related to model state variables through the hydrostatic equation h=ψ+z. The weight wi associated to the i-th control point depends on the relative importance (in terms of accuracy, for instance) of the i-th observation with respect to others. The minimization process is based on the Marquardt–Levenberg algorithm (MLA; Levenberg 1944; Marquardt 1963), whose strength is to require fewer model runs than other estimation methods (Skahill and Doherty 2006). For a nonlinear model such as CODESA-3D, MLA finds an optimal parameter set through an iterative process. The iterative procedure is stopped when a certain convergence criterion is reached (e.g., the objective function reduces to a minimum corresponding to a user-defined threshold) or alternatively after a certain number of iterations, if over the last (default 3) successive iterations the objective values are within a relative distance prescribed by a termination criteria control variable. In addition to the objective function, another measure of goodness of fit used in PEST is the correlation coefficient (R). Generally, an acceptable value of R should be above 0.9 for the fit between model outputs and observations (Hill 1998). To DOI 10.1007/s10040-014-1214-1

account for the relative importance of observations, the correlation coefficient is calculated as: NC   X wi hi −m ðwi hi −mo Þ

R¼"

i¼1

X NC

 wi hi −m 2

#0:5 "

i¼1

NC X

#0:5

ð3Þ

ðwi hi −mo Þ2

i¼1

where hi is the model-generated counterpart to the i-th observation value hi , m is the   mean value of weighted NC 1 observations m ¼ NC ∑ wi hi , mo is the mean value of i¼1

weighted model-generated counterparts to observations   NC and wi is the weight associated with 1 mo ¼ NC ∑ wi hi , i¼1

the i-th observation. Unbiased estimation of m and mo are constrained by the condition that the mean of wi must be 1. The final output from PEST reports the detailed record of the optimization, including the parameter covariance and the correlation coefficient matrices calculated during the overall process. Another important output is the parameter sensitivity, which measures composite changes in model outputs generated by variations in the value of a single parameter. Parameter sensitivity is useful to identify parameters that could degrade the performance of the optimization process through lack of sensitivity to outputs. Similarly, the observation sensitivity, which measures overall adjustable parameter change caused by changes in the value of a single observation, is useful to identify observations that are more important to the inversion process due to their information content. The analysis of both these sensitivities could help the modeller to refine the conceptual model of the studied system, to better exploit available data, and also to eventually plan field campaigns for further data acquisition (Sanz and Voss 2006).

Simulation/optimization with Carroll’s GA

The S/O model is the final module of the proposed framework and it enables the feasibility of different management options to be assessed. It provides an efficiency ranking by finding the set or a representative subset of Pareto-optimal solutions, and it applies userdefined preferences to choose the best compromise solution from the generated set. In this coupled module, CODESA-3D is the simulation model while the GA of Carroll (1996) is the optimization model. At the beginning of the genetic algorithm, a large population of random chromosomes is created. The initial population is formed by individuals, which represent possible solutions that are selected within the predetermined lower and upper bounds of each model parameter to be optimized. Solutions from Hydrogeology Journal

one population are taken and used to form a new population. This is motivated by the hope that the new population will be better than the old one. New generations of individuals are reproduced from older generations through random selection, crossover, and mutation based on certain probabilistic rules. Selection is made according to solution fitness (objective function): the more suitable a solution the more chance it has to reproduce. In this way, the population gradually evolves toward an optimal solution. Let us consider a set of n wells, where the design variables Qi (decision variables tested for optimality) and ci (modeled variables) represent respectively the discharge rate and the salt concentration at the i-th well. In this optimization study two conflicting objectives are considered: (1) maximizing pumping rates from the aquifer wells while (2) limiting the salinity of the water withdrawn, composed into a single scalar objective function (Z), made up of two weighted sums (Qahman et al. 2005): (  )  Xn Xn  Max Z Z ¼ r1 ð4Þ Qi −r2 Qi ci  i¼1 i¼1 The constraints to Eq. (4) are: 0≤Qi≤Qmax, with the maximum discharge given by a user-defined threshold (generally based on aquifer safe yield); and 0≤ ci≤cmax, with the maximum relative concentration usually given by the maximum possible concentration in the system, i.e. cmax=1. The relative weights r1 and r2 are determined by means of a parametric calculation, set r1 to 1 and run the S/O model in transient state several times with changing r2 values in the appropriate range (usually 1–200). The optimal value for the salinity weighting parameter r2 is then chosen within the adopted feasibility region (based on the case-study management goals), ensuring a relatively high value of total allowed pumping while minimizing the total salt mass extracted for a given total pumping amount—see also the discussion in section ‘Optimization (2011–2020)’. After initialization of the model, the simulation model calculates the state variables (namely groundwater heads and concentration), which are processed by the GA to calculate the fitness value. Following the selection, crossover, and mutation probabilistic rules, from the old generation a new population is created, representing the new pumping rates dataset, which is updated in the simulation model. The procedure stops when the maximum number of generations (GenMax) is achieved. The GA parameters are set as follows: population size (how many sets of possible solutions (chromosomes) are in the population in one generation) is set to 5, as suggested by Carroll (1996); crossover and mutation probabilities are set to 0.5 and 0.02, respectively, following the recommendation given by Carroll (1996); maximum number of generations is set to 200, according to a parametric analysis showing that the fitness value is not changed by further increasing the number of DOI 10.1007/s10040-014-1214-1

Fig. 1 The Gaza Strip study area and the regional aquifer (adapted from UNEP/DEWA/GRID-Geneva 2002)

generations. In this configuration, CODESA-3D is run 1,000 times (i.e., 5 possible solutions multiplied by 200 generations) to find the optimal solution for a given problem.

Gaza strip coastal aquifer The Gaza Strip is a semi-arid region located in the Mediterranean basin; it covers a coastal area of about 365 km2 between Egypt and Israel, forming a long and narrow rectangle; its length is approximately 45 km and its width ranges from 6 to 12 km (Fig. 1). Land

surface elevations gradually slope downwards from east to west, ranging from 110 m above mean sea level in the eastern areas to mean sea level in the west. The Gaza coastal aquifer is the main source of water for agriculture, domestic, and industrial purposes in the Gaza Strip. The 2010 population of the Gaza Strip is 1.6 million and the density is about 4,500 people/km2, making it one of the most overcrowded areas in the world. Due to continuing population growth, the total water demand in the Gaza Strip is strongly increasing. The current available resources do not satisfy the need of water, and this is causing a huge deficit between water demand and supply (Qahman and Larabi 2006). Aquifer overexploitation exacerbates the saltwater intrusion problem and corrective measures are urgently needed to restore groundwater quality and properly manage the aquifer. During the last decades several studies have been carried out to analyze SWI in the Gaza Strip (Yakirevich et al. 1998; Melloul and Collin 2000; Moe et al. 2001; Qahman and Larabi 2006), but the aquifer quality status is so critical (Shomar et al. 2010) that this problem is still a long way from being resolved.

Table 1 Hydraulic parameters from field campaigns for the Gaza coastal aquifer Fig. 2 Schematic geological cross-section through the study area. This cross-section could be located anywhere in the GCA, although the existence and position of the subaquifers, clay lenses, and aquitards vary from north to south. The arrow indicates the direction of undisturbed groundwater flow driven by slope of the water table. Adapted from Dan and Greitzer (1967) Hydrogeology Journal

Parameter

Value

Hydraulic conductivity (K) – sand (m/s) Hydraulic conductivity (K) – clay (m/s) Porosity – sand (%) Porosity – clay (%) Specific storage (m−1)

2.3×10−4–9.3×10−4 1.2×10−7 35 40 10−4

DOI 10.1007/s10040-014-1214-1

Table 2 Water balance of the GCA: reference mean values, in Mm3/y, for the period 2001–2010 Water balance Inflows Rainfall recharge (effective infiltration) Lateral inflow Leakage from water distribution system Agricultural return flow Wastewater return flow Loss of aquifer storage Seawater intrusion Total inflows Outflows Agricultural wells Municipal wells Discharge Total outflows Net imbalance

40–45 15–40 10–15 20 10–15 2–3 10–45 107–183

general westward direction towards the coast where fresh groundwater discharges to the Mediterranean Sea. Local flow patterns are strongly disturbed by overpumping. Large cones of depression are located in different areas within the Gaza Strip (e.g., near Gaza City in the north or Khuzaa and Rafah in the south), causing water levels below mean sea level and producing hydraulic gradients from the Mediterranean Sea towards the major pumping wells (usually municipal). About 5,000 wells are located in the area, but most of them are no longer suitable for

80 90 10–40 180–210 27–73

Geological and hydrogeological setting The GCA is part of a regional groundwater system covering areas of Palestine, Israel, and Egypt (Fig. 1). It is classified as a Pleistocene-age granular aquifer (Kurkar Group), formed by marine and eolian calcareous sandstone (kurkar), reddish silty sandstone (hamra), silts, unconsolidated sands, and conglomerates with intercalation of clay of marine origin and loam. Some of these formations are lenses which are randomly distributed in the area, causing local perched water conditions (Shomar et al. 2010). Others begin at the coast, extend 2–5 km inland, and separate the aquifer into various subaquifers, indicated as A, B1, B2, and C units (Fig. 2). Subaquifer A is phreatic, whereas B1, B2, and C become confined toward the sea. The aquifer overlies a marine clay formation of Neocene age, known as the Saqiya Clay aquitard. The average thickness of the aquifer body is about 150 m, reaching its maximum of about 200 m near Gaza City and its minimum in the eastern part (landward boundaries) and in the south, near the city of Rafah. Several pumping test campaigns have been carried out to determine hydraulic parameters (Qahman and Larabi 2006), which are summarized in Table 1. The undisturbed regional groundwater flow is mainly perpendicular to the coastline (Moe et al. 2001), with a

Fig. 3 GCA groundwater pumping rates for the period 1935–2020 Hydrogeology Journal

Fig. 4 The Gaza Strip hydrogeological model: a 3D and b 2D meshes DOI 10.1007/s10040-014-1214-1

Fig. 5 a 3D representation of the aquifer (units 1, 3, 5, and 7) and aquitard (units 2, 4, and 6) layers and b soil map

drinking purposes because the water quality is very degraded, exceeding WHO (World Health Organization) standards for both chlorides (250 mg/l), due mainly to SWI, and nitrates (50 mg/l), due to polluting activities.

that the average annual potential evapotranspiration for the Gaza Strip is around 1,500 mm/y. It is estimated that the total rainfall recharge or effective infiltration ranges between 40 and 45 Mm3/y. The water balance of the aquifer takes into account the following other inflows:

Water balance

1. Lateral inflow (LI), coming from the eastern part of the aquifer which is contact with the larger regional aquifer. Its value cannot be evaluated unequivocally as it depends on outer conditions relative to the regional aquifer, so it is derived from the hydrologic balance and is estimated to be in the range 15–40 Mm3/ y (Moe et al. 2001; Qahman and Larabi 2006). 2. Leakage from the water distribution system, which is usually estimated as a percentage of distributed water, and amounts to approximately 10–15 Mm3/y. 3. Agricultural return, which is a portion of the water used for irrigation that infiltrates into the aquifer. This value

The water resources crisis in the Gaza Strip can be expressed in terms of a net imbalance of about 27– 73 Mm3/y between inflows and outflows, referring to average climatic conditions, total abstractions, and return flows (Table 2). Rainfall (P) is the main source of freshwater for the aquifer, although a large amount is lost by evapotranspiration (ET). Rainfall occurs in the period from October to March, while the rest of the year is usually an uninterrupted dry period. The average annual rainfall varies from 400 mm/y in the north to 200 mm/y in the south. Both measurements and calculations indicate

Fig. 6 a Operational agricultural and municipal wells and b model clustered wells Hydrogeology Journal

DOI 10.1007/s10040-014-1214-1

Table 3 Statistics for the calibration simulation

3D SWI model setup

Statistic

Value

Correlation coefficient, R Mean value of residuals (m) Maximum value of residuals (m) Minimum value of residuals (m) Standard deviation of residuals (m) Number of observation points

0.97 −0.10 1.11 −1.28 0.68 43

This study introduces for the first time a detailed reconstruction of the 3D geometry of the land side of the Gaza aquifer system. The CODESA-3D mesh of the aquifer system contains 39,408 nodes and 200,445 tetrahedra (Fig. 4a) and was obtained by replicating the 2D surface mesh vertically over several layers. The 2D mesh is the result of a triangularization of the study

is estimated as a percentage of pumping water for irrigation purposes and, for the Gaza Strip, it is usually set at 25 % of the pumping at individual wells (Qahman and Larabi 2006). 4. Wastewater return flow, both from pipes and from septic systems. 5. Loss of aquifer storage, calculated as a fraction of stored water in the aquifer. 6. Seawater intrusion, as confirmed by a multipurpose groundwater-monitoring network installed throughout the region to measure groundwater levels and nitrate and chloride content. The seawater inflow is estimated to be in a range of 10–45 Mm3/y (Moe et al. 2001; Qahman and Larabi 2006). The total inflow is estimated to be in the range 107– 183 Mm3/y (Table 2). Water abstraction (Q) is the main outflow from the aquifer. Groundwater is the main source of freshwater in the Gaza Strip for agricultural, industrial, and drinking purposes. On the basis of the available data on population, population growth, number of wells, and their reported abstraction values over the last decades, the abstraction rate from pumping wells has been estimated for the period from 1935 to 2020 (Fig. 3, adapted from Qahman and Larabi 2006). In 2006 approximately 170 Mm3/y of water was pumped from about 4,600 wells for domestic and agricultural purpose (PWA 2007). Another outflow is the freshwater discharge to the sea, usually estimated (Moe et al. 2001; Qahman and Larabi 2006). The total amount of outflow is estimated to be in the range 180–210 Mm3/y.

Table 4 Calibrated values of hydraulic conductivity Layers Aquifer

Clay

K (m/s) Horizontal 1 2 3 4 5 6 7 8 9 10 11 5.00E-07

Hydrogeology Journal

2.32E-04 1.00E-04 3.15E-04 1.77E-04 1.00E-03 2.76E-04 4.28E-04 2.79E-04 1.00E-04 4.48E-04 1.00E-04 1.00E-07

Vertical 1.00E-04 1.00E-05 1.00E-04 1.91E-05 1.00E-05 7.68E-05 1.00E-04 1.00E-05 1.57E-05 1.00E-05 3.15E-05

Fig. 7 a Comparison between calibrated (dotted lines) and measured (continuous lines) interpolated groundwater levels; b scatterplot of observed heads versus simulated heads for the year 1935 DOI 10.1007/s10040-014-1214-1

Table 5 Statistics for the comparison between simulated and measured groundwater heads for 1970, 1990, 2000, and 2010 Statistics 2

Variance of residuals (m ) Standard deviation of residuals (m) Minimum of residuals (m) Maximum of residuals (m) Mean of residuals (m) Correlation coefficient (/) Number of observation points

Symbol

1970

1990

2000

2010

σ σ Min Δh Max Δh Mean Δh R –

0.50 0.70 −1.24 1.67 0.09 0.89 38

0.62 0.79 −1.55 1.60 −0.10 0.71 68

0.78 0.88 −1.73 2.75 0.51 0.82 86

1.39 1.18 −3.16 4.07 −0.04 0.94 67

2

domain (which corresponds to the administrative limits of the Gaza Strip) into 9,545 triangles and 4,926 nodes (Fig. 4b). Grid horizontal spacing ranges from 25 to 500 m, with the finest elements in proximity of major wells. The vertical discretization consists of 7 layers of varying thickness, corresponding to the main regional hydrogeological units (Fig. 5a): phreatic aquifer (unit 1, corresponding to subaquifer A in Fig. 2), aquifer units 3, 5, and 7 (corresponding to subaquifers B1, B2, and C in Fig. 2), and aquitard units 2, 4, and 6 (corresponding to the clay lenses in Fig. 2). The average thicknesses of the phreatic aquifer, underlying aquifer, and aquitard units are 50, 30, and 7.5 m, respectively. The top elevation of the first aquifer layer is spatially variable in correspondence with the land surface elevation. The domain bottom corresponds to the upper surface of the Saqiya aquitard and is treated as an impervious boundary. The dominant soil type in the GCA is sandy, with some intercalation of clay lenses. A very simple conceptual hydrogeological 3D model was first assumed, with only four hydraulic conductivity (K) parameters to characterize the horizontal and vertical components of the sand and clay units. Several attempts to calibrate these four hydraulic conductivity parameters failed to provide an acceptable configuration, mainly due to the inhomogeneous spatial distribution of the sand/clay strata. In order to keep the hydrogeological setting as simple as possible, and after exploring several parameterization strategies, not described here, some degrees of freedom were introduced in the spatial distribution of soil properties, rather than changing the vertical profiles from point to point. As a final hydrogeological setting, 11 spatial zones (derived from the soil map, Fig. 5b; Hamad et al. 2012) were identified for the aquifer units, where the vertical and horizontal conductivity components for sand (corresponding to a total of 22 parameters) were calibrated as described in the following section, while keeping common components for vertical and horizontal conductivity for clay for all spatial zones, corresponding to two additional parameters to be calibrated as well. The porosity and specific storage parameters are not treated as calibration variables and their values are as reported in Table 1. Dispersivities are normally obtained from laboratory tests or field experiments (tracer tests), but the main problem is that the dispersion term is scale-dependent. For small laboratory experiments, small values are obtained (few centimetres), while values up to several hundreds of meters are obtained in field studies. In this study constant values of longitudinal and transverse dispersivity Hydrogeology Journal

(αL=100 m and αT=20 m respectively) were chosen using grid-Peclet criterion to avoid numerical instabilities. To assess tradeoffs between accuracy and efficiency, some tests were performed halving the dispersivity size (with αL=50 m and αT=10 m) showing overall comparable simulation results with some localized numerical oscillations near the seaside. A finer grid resolution will be needed to examine in more detail the 3D behaviour of local processes and to relax the grid dispersivity constraints. The freshwater and seawater densities are set to 1,000 and 1,025 kg/l, respectively (Lecca et al. 2001; Qahman and Larabi 2006). A zero flux Neumann boundary condition is prescribed on the northern and southern boundaries, assuming that the flow lines are perpendicular to the coastline under natural conditions. Another zero flux Neumann boundary condition is set at the bottom of the aquifer system, as mentioned previously. As the aquifer is only studied in the land side of the Gaza Strip, a vertical boundary is set along the coastline, with stationary Dirichlet boundary

Fig. 8 Location of 14 representative control wells for groundwater head operational during the period 1970–2010 DOI 10.1007/s10040-014-1214-1

Table 6 Measured and simulated groundwater heads (h) for 1970, 1990, 2000, and 2010 at 14 control wells Well ID

1 2 3 4 5 6 7 8 9 10 11 12 13 14 Mean

Measured h values (m a.m.s.l.) 1970 1990 2000 3.05 0.19 2.81 0.43 3.68 1.99 2.41 0.25 2.59 3.07 1.52 0.74 0.34 0.88 1.71

−1.35 −0.82 −1.19 −0.05 −0.22 0.11 −0.52 0.19 0.36 0.17 0.21 −0.90 −0.67 −0.69 −0.38

−3.40 −0.66 −1.86 0.65 −2.21 −0.13 −0.41 −0.11 0.88 1.20 0.23 −3.01 −0.20 −0.90 −0.71

2010

Difference 2010–1970

−10.10 −1.43 −8.00 −0.44 −7.91 −1.11 −2.66 −2.52 −1.29 −0.73 −1.19 −4.32 −1.65 −1.91 −3.23

−13.15 −1.62 −10.81 −0.86 −11.59 −3.11 −5.07 −2.78 −3.88 −3.80 −2.71 −5.06 −1.99 −2.79 −4.94

conditions of constant head derived from a hydrostatic pressure along this seaside boundary and a constant concentration for seawater (in this study set, for simplicity, equal to 20 g/l of chlorides, representing the maximum relative concentration cmax=1). An assigned inflow of freshwater is prescribed at the eastern border (inland lateral flow, LI) as derived from the hydrologic balance. The inland lateral flow rate is a model boundary condition that is adjusted during the model calibration. Surface nodes are subjected to atmospheric forcing, with the net recharge rate (Re) calculated as: Re ¼ ðP−ETÞ⋅C rech ¼ Pnet⋅ C rech

ð5Þ

where P is the total rainfall, assessed according to the recorded rainfall data from gauging stations located in the area during the period 1973–2010; ET is evapotranspiration, derived from potential evapotranspiration (ET0)

Simulated h values (m a.m.s.l.) 1970 1990 2000 1.93 0.51 1.97 1.67 2.51 2.64 1.83 0.87 2.52 2.67 2.25 0.22 1.33 1.54 1.75

−2.94 −0.39 −1.12 0.94 −1.44 −0.03 −0.73 −0.72 0.44 0.84 0.76 −0.38 0.46 0.83 −0.25

−4.87 −1.57 −2.79 −0.01 −3.10 −0.50 −1.09 −0.72 0.32 0.62 0.47 −2.17 −0.88 −1.14 −1.24

2010

Difference 2010–1970

−13.17 −2.71 −7.76 −1.53 −9.74 −3.41 −3.34 −1.53 −0.97 −0.44 −0.39 −3.88 −2.70 −3.09 −3.90

−15.10 −3.22 −9.73 −3.20 −12.25 −6.05 −5.16 −2.40 −3.49 −3.11 −2.64 −4.10 −4.03 −4.63 −5.65

estimated with the Penman-Monteith equation; Pnet is the net rainfall (P−ET); and Crech is the recharge coefficient, representing the fraction of net rainfall that actually infiltrates into the aquifer, set according to soil characteristics (Fig. 5b; Hamad et al. 2012). The water retention characteristics are parameterized according to the Huyakorn et al. (1984) semi-empirical constitutive relationships. The effective saturation is given 1 by S e ¼ ð1þΛ if ψ<ψa and Se=1 if ψ≥ψa, where Þγ β β Λ=α (ψα−ψ) , ψa is the air entry pressure, and α, β, and γ are constants. The water saturation is expressed as Sw=(1−Swr)Se(ψ)+Swr, where Swr is the residual water saturation. The relative conductivity is expressed as Kr(ψ)=Sen where n is a constant. The values set for these various parameters are α=0.015, β=2, γ=3, ψa=−10 m, Swr=0.01, and n=2, and are representative of a sandy soil. Well production rates in the study region have been estimated on the basis of information on water

Fig. 9 Groundwater level contours for year 2010: interpolation of a measured and b simulated values Hydrogeology Journal

DOI 10.1007/s10040-014-1214-1

Fig. 10 Interpolated 10 g/l salt concentration isolines for a 2000 and b 2010: measured (continuous lines) and simulated (dotted lines). The 11 control wells used in Table 7 are also shown

consumption since 1935 (Qahman and Larabi 2006) and the available groundwater abstraction rates from the Gaza aquifer since the late 1980s (Fig. 3). There are approximately 5,000 operational wells but, since they are concentrated in small areas, only 1,600 clustered wells (of which 139 municipal wells) are implemented in the model (Fig. 6a,b). Clustered wells account for the average depth and the sum of pumping rates of individual wells belonging to the cluster.

Model results Calibration The GCA model has been calibrated under steady-state flow conditions against measured groundwater heads (control variables) relative to a 1935 dataset, considered a pre-development scenario (no pumping). The initial setup of the PEST control file is based on parameters assumed by Lecca and Cau (2006). As already discussed, different conceptualizations of the

domain were tested for the calibration procedure, but the best results were obtained when applying the calibration procedure with 24 adjustable K parameters (i.e., horizontal and vertical conductivity components for sand aquifers in 11 spatial zones and common components for clay lenses). The chosen zonation takes into account “potential areas” with corresponding adjustable parameters to be optimized, in order to maintain both the flexibility of the computational model and the K value in the range of the measured values (Table 1). Thus, the 24K values (described in section ‘3D SWI model setup’) are to be optimized within a prior defined range (10−6–10−3 m/s for the aquifer layers and 10−8–10−6 m/s for the clay lenses). The observations are identified by 43 groundwater head values (h ) measured in the field, equally weighted (wi=1 in Eq. 2). For the initial calibration setup no prior information has been added. Evaluation of the calibration effectiveness is based on examination of the residuals (Δh ¼ h−h ), the correlation coefficient (R, Eq. 3), and the standard deviation of residuals (σ).

Table 7 Measured and simulated normalized salt concentrations for 2000 and 2010 at 11 control wells Control well ID

Measured normalized salt concentration 2000 2010 Difference 2010–2000

1 2 3 4 5 6 7 8 9 10 11

0.0053 0.0192 0.0044 0.0064 0.0129 0.0052 0.0161 0.0193 0.0198 0.0139 0.0015

Hydrogeology Journal

0.0048 0.0268 0.0415 0.0185 0.0206 0.0156 0.0199 0.0205 0.0372 0.1870 0.0026

−0.0005 0.0076 0.0371 0.0121 0.0078 0.0105 0.0038 0.0012 0.0175 0.1731 0.0010

Simulated normalized salt concentration 2000 2010 Difference 2010–2000 0.0001 0.0023 0.0074 0.2071 0.0000 0.0630 0.1572 0.0207 0.0375 0.1871 0.0002

0.0671 0.0086 0.0444 0.2800 0.0000 0.0897 0.2174 0.0641 0.0939 0.3508 0.0038

0.0671 0.0063 0.0370 0.0730 0.0000 0.0267 0.0602 0.0434 0.0565 0.1637 0.0036

DOI 10.1007/s10040-014-1214-1

indicates that additional measures in the northeastern portion of GCA could be very useful in improving the model calibration. It should be remarked that, even though the smaller spatial zones were also found to have an effect on the calibration results, the available observation data probably does not provide sufficient coverage to constrain all 11 zones. The spatial discretization and parameterization of the model, based in this study on the soil map data, should in future work also be guided by available observation data on the state variables being simulated.

Validation (1935–2010) Fig. 11 Trade-off curve for the S/O model

Based on the evidence provided by several calibration results, which were showing highly negative residuals in the northern areas and highly positive residuals in the southern areas, it was decided to treat the lateral inflow (LI) as a non-uniform flux distributed along the eastern side of the aquifer, with higher values in the northern part and lower values in the southern area. The LI flux spatial distribution has been established according to the spatial distribution of average yearly rainfall registered at the gauging stations located in the eastern part of the aquifer. Different total values of LI were tested to identify the best one, which resulted to be about 35 Mm3/y. The main outputs of the calibration procedure are reported in Tables 3 and 4, and Fig. 7. The correlation coefficient R is very high (0.97), with a mean value for residuals of −0.10 m and a standard deviation of 0.68 m. The calibration procedure attained reasonable results with residuals having positive and negative values throughout the aquifer. Comparison of simulated and measured values for the year 1935 shows that the two water-table fields are broadly consistent, although there are some discrepancies in the southeastern part of the domain. The calibrated K values are in the same range of field measurements (Table 1) both for sandy and for clayey areas. The sensitivity results from the PEST calibration show that all parameters have quite moderate sensitivity values, with the most sensitive parameters being those relative to horizontal conductivity, as would be expected. Likewise, most of the observations have moderate sensitivity values, with the weakest sensitivity attributed to the observations closest to the coast, owing to the effect of the imposed Dirichlet boundary condition (constant head, h), and the strongest sensitivity found in proximity to the northeastern boundary, where the measured groundwater heads and imposed lateral inflows are highest. This

The GCA model has been validated with an independent control dataset in a transient-state simulation for the period 1935–2010, with initial flow conditions taken from the simulated 1935 steady-state water levels (Fig. 7a) and with the calibrated K values provided in Table 4. The initial relative salt concentration in the aquifer is assumed equal to 0 for the entire domain, as the chloride dataset for 1935 appears to represent not only SWI but also dissolved salts from agricultural return flow. The validation procedure has been split into two periods, 1935–2000 (P1) and 2001–2010 (P2). For period P1, rainfall and lateral inflow have the same mean values as for the calibration procedure; for the northern part of the aquifer, during the years from 1980 to 2000 an additional freshwater recharge of 5000– 10,000 m3/day has been added to account for the discharge of partially treated wastewater in the sandy dunes area around the Beit Lahia wastewater treatment plant. Pumping rates, displayed in Fig. 3, were reduced by 20–30 % to account for the return flow coming from irrigation, sewage infiltration, and leakage from water networks. Due to uncertainty regarding the spatial distribution of this return flow, it was simply assumed to occur close to the pumping sites and attributed proportionally to the pumping rates. In the southern part of the aquifer, the pumping rates were increased for the last part of the simulation period to account for Israeli settlements, though there is considerable uncertainty concerning these abstraction rates. The values of groundwater head and solute concentration at the end of the period P1 (1935– 2000) simulation were used as initial conditions for the subsequent period P2 (2001–2010). For period P2, the net recharge rate is calculated as a percentage of the net recharging rainfall by considering the recorded daily rainfall data from gauging stations and evaluating evapotranspiration according to daily and monthly meteorological data of the period 2001–2010. Lateral inflow is again kept constant, as there is little

Table 8 Overall performance of non-optimized and optimized pumping strategies, and relative benefits at the end of 2020 Total pumpings (m3/y) Average head (m a.m.s.l.) Total extracted salt mass (kg/y) Total area affected by SWI (m2) Hydrogeology Journal

Non-optimized

Optimized

Δ

%

120,000,000 −7.63 790,322,398 11,662,932

120,566,965 −7.29 611,660,283 11,081,997

+566,965 +0.34 −178,662,115 −580,935

+0.5 +4.5 −22.6 −5.0

DOI 10.1007/s10040-014-1214-1

Fig. 12 Pumping magnitudes for the 139 clustered municipal wells for the a non-optimized and b optimized solutions, at the end of 2020

information about groundwater flow along the western boundary and outside the Gaza Strip. Additional freshwater recharge was added over the entire domain (but mainly in the northern region) to account for wastewater treatment plant locations and other artificial recharge points. Pumping rates of municipal and agricultural wells recorded for the same period (Fig. 3) were again reduced by an amount between 20 and 30 % for each pumping well to account for the irrigation return flow. To evaluate the performance of the numerical model for the flow problem, averaged simulated and measured groundwater-head (h) values are compared for the years 1970, 1990, 2000, and 2010. The validation results are summarized in Table 5 in terms of residual statistics for all available control points for each year. As the location of

some control wells changed during the period, a selected group of 14 control wells whose locations never changed during the period 1970–2010 (Fig. 8) is assumed as representative of the groundwater system, in order to also quantify and compare the evolution of GCA over the simulation period (Table 6). Figure 9 shows the comparison between measured and simulated water levels for the year 2010. It can be seen that there is a reasonable match between simulated and measured values for 2010, although there are, again, some discrepancies in the south-eastern part of the domain, which is difficult to model due to the various sources of uncertainty mentioned previously. For the 14 control wells, the mean difference between average measured heads from 1970 to 2010 is – 4.94 m, while the mean difference between average

Fig. 13 Interpolated groundwater head contours in m a.m.s.l. for the a non-optimized and b optimized models, at the end of 2020 Hydrogeology Journal

DOI 10.1007/s10040-014-1214-1

simulated heads from 1970 to 2010 is –5.65 m (Table 6). The model is thus simulating an excess of about 0.7 m of water level decrease in the period 1970–2010, or around 2 cm/y. To assess the sensitivity of the model to storage values, the model was run in the period 1935–2000 for Ss values of 10−3, 10−4 (the reference value as reported in Table 1), and 10−5. The difference in average groundwater head between the reference case and each of the other two Ss cases for the observation points of Table 5 and for the 14 representative wells of Fig. 8 at the end of the 65-year simulation was found to be less than 1 cm/y, compared to an average water table lowering of about 7 cm/y over this same period. The model is thus relatively insensitive to this storage parameter. For the solute problem, simulated and observed salt concentration isolines for 2000 and 2010 are compared, with Fig. 10 showing the 10 g/l isolines. In addition, point concentrations for the same years are compared at 11 control wells (Table 7; the location of the control wells is shown in Fig. 10). Both measured and simulated fields show a general encroachment of seawater, stronger in the northern zone where high pumping rates have significantly altered the natural SWI patterns.

Optimization (2011–2020) After the calibration and validation phases, the simulation/ optimization module of the computational framework was set up for the Gaza coastal aquifer for the period 2011– 2020. The 2010 validation results (i.e. groundwater heads and salt concentrations) are used to set the initial conditions for the S/O problem. The same lateral inflow and rainfall rates of the previous period (2001–2010) are considered. Pumping rates are estimated on the basis of current projections of agricultural and domestic water demand until 2020 (Fig. 3). Considering that 5,000 wells are operating in the area (Fig. 6a), in this study it is assumed that no further well will be installed until 2020; thus, the location of wells is kept constant during the optimization phase. The agricultural well pumping rates are kept constant (at around 77 Mm3/y), and only the pumping rates at the 139 clustered municipal wells are treated as decision variables, the reason being that in an operational context these can be more easily managed. The maximum allowed discharge value (Qmax) is based on the estimated mean maximum operational rate of overall wells, while the maximum allowed relative concentration is set as cmax=1, i.e., the maximum possible concentration in the system. After model setup and initialization, the transient management scenario (with eight different values of salinity weight r2, ranging between 1 and 100) is simulated for the 10-year period. In this study, the reasonable limits of the feasible region (Fig. 11) are: (1) for the pumping rates, between 95 and 105 % of the prescribed total pumping amount, corresponding to the estimated municipal abstraction at the end of 2020 (around 120 Mm3/y); and (2) for the total extracted salts, Hydrogeology Journal

the allowed upper level identified by the non-optimized value (around 790×106 kg/y). The optimal value of the salinity weighting parameter r2 is 15, ensuring in this way around 100 % of the prescribed total pumping amount, while keeping the total salt mass extracted at wells under the reasonable upper limit of the feasibility region. The results (Table 8 and Figs. 12, 13, and 14) show that the optimized model has, compared to the nonoptimized solution, a different distribution of municipal pumping rates but a total amount (around 120.5 Mm3/y) close to the prescribed upper limit, while reducing total salt extracted by about 23 %. At the same time, averaged groundwater heads show an increase of about 0.34 m (4.5 %) and the total area affected by SWI shows a reduction of about 5 %. Compared to the non-optimized solution, the major cones of depression in the northern area (around Gaza City) and in the southern area (Rafah and Khan Younis) for the optimized solution have much-reduced head gradients and head drawdowns (Fig. 13). For the northern area, the minimum groundwater level increases from less than 8 m below a.m.s.l. to 6 m below a.m.s.l., while for the southern part, the minimum groundwater level increases from less than 20 m below a.m.s.l. to 12 m below a.m.s.l. In line with this result, in the two corresponding coastal areas, the optimized solution is able to contain SWI. In the northern area in particular, where the calibrated and validated model performed best, the intrusion front for the non-optimized model is up to 300 m further inland compared to the optimized solution (Fig. 14). As is apparent from Fig. 12, the optimized solution compensates reduced pumping in the critical north and south regions of the GCA with increased pumping in the less vulnerable (lower population, fewer aquifer stress points) central part of the domain. In this area, pumping rates in some

Fig. 14 Interpolated normalized salt concentration contours for the non-optimized and optimized model at the end of 2020 DOI 10.1007/s10040-014-1214-1

wells are more than doubled with respect to the nonoptimized situation. As a result, the eastern flank of this area shows a moderate decrease in groundwater levels for the optimized solution (Fig. 13) and the central coastal plain shows some SWI impact as well (Fig. 14). Overall, the optimized solution provides for a general smoothing of pumping rates, lowering the higher values and increasing the lower values. No consideration is given here to additional pumping or to pipeline infrastructures; implementing the optimization plan would require that water be piped overland, and the feasibility of this is not addressed in the current study.

Conclusions The hydrogeological model of the Gaza coastal aquifer was set up using field data covering the period 1935– 2010. The model was calibrated under steady-state conditions using 1935 groundwater levels and considering average climate conditions and a natural (no groundwater pumping) scenario. The adjustable parameters in the calibration procedure were heterogeneous and anisotropic hydraulic conductivity values and the lateral inflow boundary condition for the eastern boundary of the GCA. Validation simulations were then carried out in transient mode for two time periods (1935–2000 and 2001–2010) using spatially and temporally variable recharge and pumping rates. Reasonable agreement was obtained for groundwater levels and coastal salt concentrations, with the largest discrepancies noted in the southern portion of the model domain, a region where field data is also highly uncertain. In the final phase of the work, the simulation model was linked to the optimization tool and applied over the period 2011–2020. The main goal of this part of the study was to find possible alternative management options considering: the same municipal wells currently operating in the area; the same location for all operational wells; the maximum allowed pumping rate at each well estimated on the basis of current mean upper operational rates. The S/O analysis identified an optimal scheme for the spatial distribution of water abstraction from the municipal wells for the 2011–2020 period that, compared to the nonoptimized solution (continuation of groundwater exploitation as per previous periods), keeps overall abstraction close to the user defined total amount while significantly lowering the total mass of salt extracted (by 23 %). Moreover, the optimized solution increases the average groundwater heads (by 4.5 %) and decreases the extent of the areas affected by SWI (by 5 %) compared to the nonoptimized situation. These benefits are most evident in the northern and southern areas of the GCA, which are currently also the regions of the Gaza Strip aquifer most affected by aquifer overexploitation and seawater contamination. Ongoing research on modeling the GCA will investigate management schemes under climate change scenarios provided by multimodel climatic audits for the case study (Sulis et al. 2012; Deidda et al. 2013). Hydrogeology Journal

Acknowledgements This study has been partially funded by the FP7-ENV-2009-1 project “Climate induced Changes on the Hydrology of Mediterranean Basins (CLIMB)” (GA 244151) and the Sardinian Regional Authority. The authors gratefully acknowledge the collaboration of Prof. Samir Afifi (Gaza Islamic University of Gaza) as the Palestinian CLIMB project coordinator.

References Ahlfeld DP, Mulligan AE (2000) Optimal management of flow in groundwater systems. Academic, San Diego Antonellini M, Mollema P, Giambastiani B, Bishop K, Caruso L, Minchio A, Pellegrini L, Sabia M, Ulazzi E, Gabbianelli G (2008) Salt water intrusion in the coastal aquifer of the southern Po Plain, Italy. Hydrogeol J 16:1541–1556. doi:10.1007/ s10040-008-0319-9 Barlow P (2005) Use of simulation-optimization modeling to assess regional ground-water systems. US Geol Surv Fact Sheet 2005– 3095 Bear J, Verruijt A (1987) Modeling groundwater flow and pollution. Reidel, Dordrecht, The Netherlands Bear J, Cheng AHD, Sorek S, Ouazar D, Herrera I (1999) Seawater intrusion in coastal aquifers: concepts, methods, and practices. Kluwer, Dordrecht, The Netherlands Bhattacharjya RK, Datta B (2005) Optimal management of coastal aquifers using linked simulation optimization approach. Water Resour Manag 19(3):295–320. doi:10.1007/s11269-005-3180-9 Carrera J, Neuman SP (1986) Estimation of aquifer parameters under transient and steady-state conditions: 2. uniqueness, stability and solution algorithms. Water Resour Res 22(2):211–227 Carrera J, Alcolea A, Medina A, Hidalgo J, Slooten LJ (2005) Inverse problem in hydrogeology. Hydrogeol J 13:206–222. doi:10.1007/s10040-004-0404-7 Carrera J, Hidalgo JJ, Slooten LJ, Vázquez-Suñé E (2009) Computational and conceptual issues in the calibration of seawater intrusion models. Hydrogeol J 18:131–145. doi:10.1007/s10040-009-0524-1 Carroll DL (1996) Genetic algorithms and optimizing chemical oxygen-iodine lasers. In: Developments in theoretical and applied mechanics, vol XVIII. The University of Alabama, Tuscaloosa, AL, pp 411–424 Cau P, Lecca G, Muscas L, Barrocu G, Uras G (2002) Saltwater intrusion in the Oristano plain (Sardinia). In: Proceedings of the Salt Water Intrusion Meeting (SWIM), Delft, The Netherlands, 6–10 May 2002 Cheng AHD, Halhal D, Naji A, Ouazar D (2000) Pumping optimization in saltwater intruded coastal aquifers. Water Resour Res 36(8):2155–2165 Custodio E (2010) Coastal aquifers of Europe: an overview. Hydrogeol J 18(1):269–280. doi:10.1007/s10040-009-0496-1 Dan J, Greitzer Y (1967) The effect of soil landscape and Quaternary geology on the distribution of saline and fresh water aquifers in the Coastal Plain of Israel. Publ. 670, TAHAL, Tel Aviv Das A, Datta B (1999) Development of multiobjective management models for coastal aquifers. J Water Resour Plan Manag 125:76–87 Das A, Datta B (2001) Application of optimisation techniques in groundwater quantity and quality management. Sadhana 26(4):293–316. doi:10.1007/BF02703402 Deidda R, Marrocu M, Caroletti G, Pusceddu G, Langousis A, Lucarini V, Puliga M, Speranza A (2013) Climate model validation and selection for hydrological applications in representative Mediterranean catchments. Hydrol Earth Syst Sci Discuss 10:9105–9145. doi:10.5194/hessd-10-9105-2013 Dhar A, Datta B (2009) Saltwater intrusion management of coastal aquifers. I: linked simulation–optimization. J Hydrol Eng 14:1263–1272 DOI 10.1007/s10040-014-1214-1

Doherty J (2002) PEST model-independent parameter estimation. Watermark, Brisbane, Australia Duan Q, Sorooshian S, Gupta VK (1992) Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resour Res 28(4):1015–1031 Gambolati G, Putti M, Paniconi C (1999) Three-dimensional model of coupled density-dependent flow and miscible salt transport. In: Seawater intrusion in coastal aquifers: concepts, methods, and practices. Kluwer, Dordrecht, The Netherlands, pp 315–362 Goldberg DE (1989) Genetic algorithms in search optimization and machine learning. Addison Wesley, Reading, UK Gorelick SM (1983) A review of distributed parameter groundwater management modeling methods. Water Resour Res 19(2):305– 319. doi:10.1029/WR019i002p00305 GRIDA3 (2011) GRIDA3: Gestore di RIsorse condivise per Analisi di dati e Applicazioni Ambientali [ GRIDA3: a shared resources manager for environmental data analysis and applications]. http://grida3.crs4.it. 6 November 2014 Hamad JT, Eshtawi TA, Abushaban AM, Habboub MO (2012) Modeling the impact of land-use change on water budget of Gaza Strip. J Water Res Prot 4:325–333 Hill MC (1998) Methods and guidelines for effective model calibration. US Geol Surv Water Resour Invest Rep 98-4005 Huyakorn PS, Thomson SD, Thompson BM (1984) Techniques for making finite elements competitive in modeling flow in variably saturated porous media. Water Resour Res 20(8):1099–1115 Kerrou J, Lecca G, Murgia F, Renard P (2007) Grid-enabled simulation of the impact of exploitation uncertainty on the seawater intrusion of the Korba aquifer (Tunisia). Proceedings of IST-Africa Conference, Maputo, Mozambique, May 2007 Lecca G (2000) Implementation and testing of the CODESA-3D model for density dependent flow and transport problems in porous media. CRS4-TECH-REP-00/40, Center for Advanced Studies, Research and Development (CRS4), Cagliari, Italy Lecca G, Cau P (2006) Automatic calibration of a 3D groundwater model applied to the Muravera-Flumendosa coastal aquifer (SE Sardinia, Italy). In: Proceedings of the XVI International Conference on Computational Methods in Water Resources, Copenhagen, Denmark, June 2006 Lecca G, Berjamy B, Paniconi C, El Hebil A (2001) Numerical modeling of seawater intrusion in the Sahel region of the Atlantic coast of Morocco. Proceedings of the First International Conference on Saltwater Intrusion and Coastal Aquifers – Monitoring, Modeling, and Management, Essaouira, Morocco, 23–25 April 2001 Lecca G, Lai C, Murgia F, Biddau R, Fanfani L, Maggi P (2009) AQUAGRID: an extensible platform for collaborative problem solving in groundwater protection. Earth Sci Inform 2(1–2):83– 95 Lecca G, Petitdidier M, Hluchy L, Ivanovic M, Kussul N, Ray N, Thieron V (2012) Grid computing technology for hydrological applications. J Hydrol 403(1–2):186–199 Levenberg K (1944) A method for the solution of certain problems in least squares. Q Appl Math 2:164–168 Marquardt D (1963) An algorithm for least-squares estimation of nonlinear parameters. SIAM J Appl Math 11:431–441. doi:10.1137/0111030 Melloul A, Collin M (2000) Sustainable groundwater management of the stressed Coastal aquifer in the Gaza region. Hydrol Sci J 45(1):147–159 Moe H, Hossain R, Fitzgerald R, Banna M, Mushtaha A, Yaqubi A (2001) Application of a 3-dimensional coupled flow and transport model in the Gaza Strip. In: First International Conference on Saltwater Intrusion and Coastal Aquifers – Monitoring, Modeling, and Management, 23–25 April 2001, Essaouira, Morocco Ndambuki JM, Otieno FAO, Stroet CBM, Veling EJM (2000) Groundwater management under uncertainty: a multi-objective approach. Water SA 26(1):35–42 Palestinian Water Authority (PWA) (2007) Agricultural and municipal water demand in Gaza governorates for 2006. PWA, Ramallah, Palestine Hydrogeology Journal

Paniconi C, Khlaifi I, Lecca G, Giacomelli A, Tarhouni J (2001) Modeling and analysis of seawater intrusion in the coastal aquifer of eastern Cap-Bon, Tunisia. Transp Porous Media 43(1):3–28 Poeter EP, Hill MC (1997) Inverse models: a necessary next step in groundwater modeling. Ground Water 35(2):250–260 Putti M, Paniconi C (1995) Picard and Newton linearization for the coupled model of saltwater intrusion in aquifers. Adv Water Resour 18(3):159–170 Qahman K, Larabi A (2006) Evaluation and numerical modeling of seawater intrusion in the Gaza aquifer (Palestine). Hydrogeol J 14:713–728. doi:10.1007/s10040-005-003-2 Qahman K, Larabi A, Ouazar D, Naji A, Cheng AHD (2005) Optimal and sustainable extraction of groundwater in coastal aquifers. Stoch Environ Res Risk Assess 19(2):99–110 Qahman K, Larabi A, Ouzar D, Naji A, Cheng HD (2009) Optimal extraction of groundwater in Gaza coastal aquifer. J Water Res Prot 4:249–259. doi:10.4236/jwarp.2009.14030 Sanford WE, Pope JP (2010) Current challenges using models to forecast seawater intrusion: lessons from the Eastern Shore of Virginia, USA. Hydrogeol J 18:73–93. doi:10.1007/s10040009-0513-4 Sanz E, Voss CI (2006) Inverse modeling for seawater intrusion in coastal aquifers: insights about parameter sensitivities, variances, correlations and estimation procedures derived from the Henry problem. Adv Water Resour 29:439–457 Shamir U, Bear J, Gamliel A (1984) Optimal annual operation of a coastal aquifer. Water Resour Res 20(4):435–444 Shomar B, Abu Fakher S, Yahya A (2010) Assessment of groundwater quality in the Gaza Strip, Palestine using GIS mapping. J Water Res Prot 2:93–104. doi:10.4236/ jwarp.2010.22011 Skahill BE, Doherty J (2006) Efficient accommodation of local minima in watershed model calibration. J Hydrol 329:122–139. doi:10.1016/j.jhydrol.2006.02.005 Sulis M, Paniconi C, Marrocu M, Huard D, Chaumont D (2012) Hydrologic response to multimodel climate output using a physically based model of groundwater/surface water interactions. Water Resour Res 48:W12510. doi:10.1029/ 2012WR012304 Sumner NR, Fleming PM, Bates BC (1997) Calibration of a modified SFB model for twenty-five Australian catchments using simulated annealing. J Hydrol 197(1–4):166–188. doi:10.1016/S0022-1694(96)03277-5 UNEP/DEWA/GRID-Geneva (2002) United Nations Environment Program, Division for Environmental Warning and Assessment, Global Resource Information Database (GRID) Network, December 2002. http://www.grid.unep.ch. 6 November 2014 Werner AD, Gallagher MR (2006) Characterization of sea-water intrusion in the Pioneer Valley, Australia using hydrochemistry and three-dimensional numerical modeling. Hydrogeol J 14:1452–1469. doi:10.1007/s10040-006-0059-7 Werner AD, Bakker M, Post VEA, Vandenbohede A, Lu C, AtaieAshtiani B, Simmons CT, Barry DA (2012) Seawater intrusion processes, investigation and management: recent advances and future challenges. Adv Water Resour 51:3–26. doi:10.1016/ j.advwatres.2012.03.004 Willis R, Finney BA (1988) Planning model for optimal control of saltwater intrusion. J Water Resour Plan Manag 114(2):163–178 Yakirevich A, Melloul A, Sorek S, Shaath S, Borisov V (1998) Simulation of seawater intrusion into the Khan Yunis area of the Gaza Strip coastal aquifer. Hydrogeol J 6:549–559 Yang R, Kalin M, Zhang Y, Lin X, Zou L (2001) Multi-objective optimization for sustainable groundwater resource management in a semiarid catchment. Hydrol Sci J 46(1):55–72. doi:10.1080/ 02626660109492800 Zghibi A, Zouhri L, Tarhouni J (2011) Groundwater modeling and marine intrusion in the semi-arid systems (Cap-Bon, Tunisia). Hydrol Process 25:1822–1836. doi:10.1002/hyp.7948 Zhou X, Chen M, Liang C (2003) Optimal schemes of groundwater exploitation for prevention of seawater intrusion in the Leizhou Peninsula in southern China. Environ Geol 43:978–985 DOI 10.1007/s10040-014-1214-1