Intensively exploited Mediterranean aquifers: resilience to seawater intrusion and proximity to critical thresholds

. We investigate seawater intrusion in three promi-nent Mediterranean aquifers that are subject to intensive exploitation and modiﬁed hydrologic regimes by human activities: the Nile Delta, Israel Coastal and Cyprus Akrotiri aquifers. Using a generalized analytical sharp interface model, we review the salinization history and current status of these aquifers, and quantify their resilience/vulnerability to current and future seawater intrusion forcings. We iden-tify two different critical limits of seawater intrusion under groundwater exploitation and/or climatic stress: a limit of well intrusion , at which intruded seawater reaches key locations of groundwater pumping, and a tipping point of complete seawater intrusion up to the prevailing groundwater divide of a coastal aquifer. Either limit can be reached, and ul-timately crossed, under intensive aquifer exploitation and/or climate-driven change. We show that seawater intrusion vulnerability for different aquifer cases can be directly com-pared in terms of normalized intrusion performance curves . The site-speciﬁc assessments show that (a) the intruding seawater currently seriously threatens the Nile Delta aquifer, (b) in the Israel Coastal aquifer the sharp interface toe ap-proaches the well location and (c) the Cyprus Akrotiri aquifer is currently somewhat less threatened by increased seawater intrusion.


Abstract.
We investigate seawater intrusion in three prominent Mediterranean aquifers that are subject to intensive exploitation and modified hydrologic regimes by human activities: the Nile Delta, Israel Coastal and Cyprus Akrotiri aquifers. Using a generalized analytical sharp interface model, we review the salinization history and current status of these aquifers, and quantify their resilience/vulnerability to current and future seawater intrusion forcings. We identify two different critical limits of seawater intrusion under groundwater exploitation and/or climatic stress: a limit of well intrusion, at which intruded seawater reaches key locations of groundwater pumping, and a tipping point of complete seawater intrusion up to the prevailing groundwater divide of a coastal aquifer. Either limit can be reached, and ultimately crossed, under intensive aquifer exploitation and/or climate-driven change. We show that seawater intrusion vulnerability for different aquifer cases can be directly compared in terms of normalized intrusion performance curves. The site-specific assessments show that (a) the intruding seawater currently seriously threatens the Nile Delta aquifer, (b) in the Israel Coastal aquifer the sharp interface toe approaches the well location and (c) the Cyprus Akrotiri aquifer is currently somewhat less threatened by increased seawater intrusion.

Introduction
The generally slow rate of water movement and storage are innate advantages of groundwater resources in terms of resilience to climate variations and relative to surface water resources. Groundwater is often of potable quality, and thus does not require expensive treatment, permitting scaled development upon demand and use of infrastructure that is normally of lower cost than that of surface water supplies (Taylor et al., 2009). However, aquifers in general and especially coastal aquifers are often exploited so intensively that their natural hydrological regime is strongly disturbed and may be thrown out of balance (Custodio, 2010). In particular, changes of seawater intrusion are highly non-linear and exhibit important thresholds, or tipping points, beyond which full seawater intrusion into a coastal aquifer may occur in response to even small sea level and/or groundwater management changes (Mazi et al., 2013). The possibility of such acceleration into full seawater intrusion is a major concern for the resilience and sustainability of coastal populations that depend on groundwater for their water supply.
Several aquifers along the densely populated Mediterranean coasts are already suffering seawater intrusion. Acceleration of this phenomenon could also be particularly large here, because the Mediterranean region, and especially its semi-arid areas, is likely to be seriously affected by decline in water resources (Kundzewicz and Döll, 2009). Human vulnerability to the decrease of renewable groundwater resources induced by climate change, as projected by climate models up to 2055, is expected to be higher for the North African rim of the Mediterranean (Döll, 2009).
In this paper, we study three coastal aquifers in the southeastern Mediterranean region that are already suffering seawater intrusion, with the aim to determine how close these aquifers are to reaching critical seawater intrusion extents in relation to different possible groundwater pumping locations (Koussis et al., 2010a, b) and/or to tipping points for complete seawater intrusion in the aquifer (Mazi et al., 2013). The studied aquifers ( Fig. 1) are (1) that underlying the Nile Delta (the Nile Delta aquifer, referred to as NDA), which is of profound importance to Egypt; (2) the Israel Coastal aquifer (ICA), which is considered to be the most important freshwater source in Israel; and (3) the Cyprus Akrotiri aquifer (CAA), which is a major source of potable and irrigation water in Cyprus. These aquifers have been subject to extensive exploitation and modified hydrologic regime by human activities and constructions (e.g., dams), and are located in a sensitive region of the world, both in climatic and political respects.
We here review the salinization history and current status of these aquifers, and quantify their resilience/vulnerability to current and future exploitation schemes, in terms of pumping location and rate of extraction, and/or hydroclimatically changed recharge conditions, in view of their respective tipping points and other critical extents of seawater intrusion.  (2) east Nile Delta (adapted from ; (b) Israel Coastal aquifer, light grey area (left) (adapted from Yechieli and Sivan, 2011); and (c) Cyprus Akrotiri aquifer.

Site descriptions
In this section we summarize the physical characteristics of the three regional coastal aquifers ( Fig. 1) and present the manner of their exploitation and the resulting current and possible future development conditions.

The Nile Delta aquifer (NDA)
The Nile Delta area (22 000 km 2 ) is among the largest deltas in the world. This area is important as the outflow area of one of the largest rivers of the world, and because it concentrates 45 % of the total population, 45 % of the arable land, 50 % of the industrial production, 40 % of the agricultural production and 60 % of the fish catch in Egypt (Mikhailova, 2001). The Nile Delta has its apex near Cairo, approximately 200 km from the Mediterranean Sea, and its base at the Mediterranean coastline extending from Alexandria to Port Said, a distance of about 240 km (Fig. 2a). The area of the Nile Delta is bounded to the east by the Ismailiya and the Suez canals, to the west by the El Nubariya Canal, and to the north by the Mediterranean Sea. The Damietta branch of the Nile, flowing to the northeast, and the Rosetta branch, flowing to the northwest, divide the NDA into a middle, an eastern, and a western part.
The average annual rainfall in the Nile Delta ranges from 180 mm at the coast to ca. 26 mm in the Cairo area, occurring mostly in the winter (Sherif and Singh, 1999 (Mikhailova, 2001). The aquifer is recharged by infiltration from a network of irrigation conduits, from excess irrigation water and limited precipitation, and by flows from the Nile Valley aquifer. The main aquifer consists of unconsolidated coarse sands and gravels (with occasional clay lenses), having an estimated porosity varying between 25 and 40 % and thinning towards the east and west sides of the delta. The near-surface deposits change in the seaward direction into impervious clays and silt. The aquifer is mostly unconfined, but it is considered leaky and partially semi-confined close to the sea. The aquifer base is a clay aquiclude that is inclined towards the north with a slope of 0.3-0.4 %. The aquifer depth at the Mediterranean coast increases from west to east from ca. 600 to ca. 1000 m, near Cairo it is about 200 m and at its apex drops to under 100 m Sherif and Singh, 1999). Groundwater flows from the apex of the delta to the sea in a roughly radial pattern (in plan view) and its quality deteriorates towards the north mainly as a result of seawater intrusion due to over-pumping, but also due to leakage of domestic wastewater (septic tanks), discharges of untreated industrial wastewater (from industries located in the greater Cairo area) and fertilizers used in agriculture. The piezometric head field in the aquifer fluctuates mildly (±1 m) in an annual cycle (Kashef, 1983). A total volume of about 2.4 km 3 (Mikhailova, 2001) is extracted annually from the aquifer for all uses. An accurate balance for the Nile Delta cannot be calculated, due to unreliable or missing data (Mikhailova, 2001), and only a few components can be estimated with some certainty. Sherif et al. (2012) report that, today, seawater intrusion has reached 100 km from the coast.

The Israel Coastal aquifer (ICA)
The ICA (Fig. 2b) is considered the most important freshwater source among Israel's major water reservoirs (the next two are the Lake of Galilee and the Mountain aquifer), covering an area of 1900 km 2 , with a length of 126 km and an average width of 15 km (Assouline and Shavit, 2004). We refer the reader to the description of the ICA by Yechieli et al. (2010), here only outlining its major characteristics. The ICA slopes seaward at an inclination of 0.01 (1 % slope), its thickness reaches 200 m at the coast, diminishing to a few meters at its land-boundary to the east. The aquifer belongs to the Kurkar Group, consisting of inter-layered sandstone, calcareous sandstone, siltstone and red loam, alternating with continental and marine clays that overlie impervious marine clays of the Saqiye Group of the Pliocene age. To a distance of 5-8 km east of the shoreline, clay inter-layers subdivide the aquifer into four sub-aquifers, the upper ones of which are unconfined, while certain of the lower ones are confined. The general flow direction in the aquifer is from the east towards the Mediterranean Sea to the west. The average annual precipitation ranges from 550 mm yr −1 in the north, to 300 mm yr −1 in the south, while the infiltrated natu-ral recharge is about one-third to one-half of the precipitation (Aberbach and Sellinger, 1967).
The chlorinity of groundwater in the aquifer did not exceed 100 mg L −1 Cl − during the 1930s (Kass et al., 2005). Since then, however, the aquifer has been under extensive exploitation, and during the 1960s groundwater withdrawal was up to 400 × 10 6 m 3 yr −1 (or twice the natural recharge of the aquifer); in 1962, about 480 × 10 6 m 3 yr −1 of groundwater was extracted from the aquifer (Aberbach and Sellinger, 1967). In some locations, over-exploitation of groundwater has caused the water level to drop by 15-20 m and seawater to intrude 2.5 km from the coast (Aberbach and Sellinger, 1967). Assouline and Shavit (2004) report monitoring data for groundwater quality that indicate a steadily increasing mean salinity; in 2004, groundwater salinity was around 200 ppm Cl − and according to some predictions may reach 300 ppm Cl − by 2020.
To enhance the water balance of the aquifer, Israel practices artificial recharge, mostly in the winter, through wells and spreading grounds, with water conveyed annually by the National Water Carrier (NWC) from the Lake of Galilee (earlier also floodwater). Water brought by the NWC in 1973 was reported to be 80 × 10 6 m 3 yr −1 (Sellinger and Aberbach, 1973), corresponding to 40 mm yr −1 on average over the whole aquifer area; other reports estimate this contribution to be 100-150 × 10 6 m 3 yr −1 (Melloul and Zeitoun, 1999), corresponding to 50-80 mm yr −1 . Furthermore, 300 × 10 6 m 3 yr −1 of water is conveyed annually by the NWC for desalination and injection in ICA. After the 1990s, 35 × 10 6 m 3 yr −1 are artificially recharged (Assouline and Shavit, 2004) while natural recharge is calculated to be 275-300 × 10 6 m 3 yr −1 (Sellinger and Aberbach, 1973), or 200 mm yr −1 (Yechieli et al., 2010).

The Cyprus Akrotiri aquifer (CAA)
The CAA (Fig. 2c) extends over 40 km 2 in the southern edge of the Akrotiri Basin (78 km 2 ), which forms the southernmost peninsula of Cyprus. The Akrotiri Basin has a smooth and hilly relief, with average altitude over 200 m a.m.s.l. in the northern part of the basin. The Salt Lake, located in the middle of the peninsula, is a topographic low, with mean water level below sea level, receiving drainage from the Akrotiri Basin (Mazi, 2000).

K. Mazi et al.: Intensively exploited Mediterranean aquifers
The Akrotiri Basin and aquifer geology, hydrological budget, and water management scenarios and conditions have been studied extensively by Mazi (2000), Koussis (2001), Mazi et al. (2004a, b), and Koussis et al. (2010bKoussis et al. ( , 2012. The alluvial aquifer consists of river deposits, gravels with layers of marls, and sands with boulders interbedded with lenses of silts and clays. The impermeable base of the aquifer (alternating marls, chalks, chalky marls and marly chalks) slopes towards the south at about 1.7 % and the saturated aquifer thickness ranges from 10 m (water table at 50 m a.m.s.l.) at its northern edge to more than 100 m near the salt lake in the south. The mean annual precipitation  over the basin ranges from over 500 mm at its northern end to 420 mm near the salt lake.
Originally, the CAA was replenished by rainfall, leakage from the Kouris River in the adjacent basin, inflows from its northern boundary and return irrigation flows. The closure of the Kouris Dam in 1988 caused a 66 % decrease in the replenishment of the aquifer by Kouris flows (Balasha and Phedonos, 1992). In response to declining groundwater levels and attendant seawater intrusion in the aquifer, the Water Development Department of the Republic of Cyprus initiated artificial recharge. Artificial recharge was applied to spreading grounds and ponds in the main aquifer area, when water from external sources was available. Nowadays, the CAA is replenished by various supplementary sources: releases from the Kouris Dam and Kouris reservoir losses, and artificial recharge (effluents from Limassol's wastewater treatment plant and water from the Garyllis and Yermasogia reservoirs, located outside the Akrotiri Basin). Almost no artificial recharge took place during the drought years 1998 and 1999, however, and groundwater pumping from the CAA was then also reduced, falling to 10 × 10 6 m 3 yr −1 in 2000 (Mazi et al., 2004), well below the 14.5 × 10 6 m 3 yr −1 mean annual extraction rate during the period 1967-1977(Jacovides et al., 1982. The data of Mazi et al. (2004) are fairly compatible with those used by Milnes (2011) to assess the salinization risk of the Akrotiri aquifer.
Our study focuses on the area of Zakaki, west of Limassol (the second largest city of Cyprus). During past dry years, groundwater from this eastern part of the CAA was used as an additional potable water source, with groundwater salinity at 500 m from the coast reaching several thousand ppm of total dissolved solids (TDS) in the early 1990s. As a consequence, after 1997, groundwater withdrawals were reduced by 25 %, to counteract deterioration of groundwater quality, with salinity then stabilizing at 1000 ppm TDS in boreholes at about 1000 m from the coast.

Aquifer similarities and differences
The aquifers studied here are all important freshwater sources for the local populations and for the sustainability and prosperity of their economies. The natural regimes of the aquifers, replenishment from rainfall, and inflows from neighboring freshwater sources, have all been altered by human activities during the last 50 years, in the NDA through the Assuan Dams, in the CAA through the Kouris Dam, and in the ICA and CAA also through artificial recharge. All three are further unconfined sloping aquifers, from which large quantities of fresh groundwater are extracted annually. In the NDA, groundwater pumping is distributed over the area of the whole delta; for the ICA and the CAA, more concentrated pumping zones can be identified.
Among the three aquifers, the NDA area is the largest (22 000 km 2 ), its depth exposure to the sea is greatest (1000 m), the inclination of its base is lowest (0.3-0.4 %) and it is much more permeable (mean hydraulic conductivity K ≈ 100 m d −1 ) than the other two. The intruding seawater front has here advanced many tens of kilometers inland from the coast in the last 50 years, and this movement has accelerated in the last 20 years, reaching today a distance of ∼ 100 km from the coast. Another important difference from the other two aquifers is the absence of a serious attempt to manage the NDA, which is in contrast to the importance of this aquifer for the regional water supply; the lack of management is also manifested in the absence of credible data for calculation of an overall water balance for NDA.
For direct aquifer comparison, the ICA area is about 1900 km 2 , its depth exposure to the sea is about 200 m, the inclination of its base is 1 % and its mean hydraulic conductivity K = 30 m d −1 . The aquifer is recharged by precipitation and artificial recharge applied through wells and spreading grounds. No inflows occur through the land-boundary. An extensive monitoring network is in place and measurements indicate that the seawater front reaches today up to 2.5 km from the coast. Several studies have been performed to estimate the permissible withdrawal from the ICA.
Furthermore, the CAA is relatively small (40 km 2 ); its exposure to the sea is 50 m, the inclination of its base is 1.7 % and its mean hydraulic conductivity K = 28 m d −1 . The aquifer is recharged by precipitation and by freshwater inflows from its northern boundary, as well as by artificial recharge. In CAA there is an established network for monitoring the groundwater level and quality.

Modeling approach
This study focuses on a first-order vulnerability assessment of seawater intrusion in the three studied aquifers. These aquifers have also been studied previously by several scientists using variable-density models such as 2D-FED: NDA (Sherif and Singh, 1999), FEFLOW: NDA (Sherif et al., 2012), ICA (Yechieli et al., 2010), SUTRA: CAA (Prieto, 2005;Prieto et al., 2006). However, the present study is the first that addresses their vulnerability in direct relation to the high non-linearity of the seawater intrusion process with associated tipping points for complete aquifer intrusion (Mazi et al., 2013), in conjunction with consideration of critical seawater intrusion to different possible groundwater pumping locations (Koussis et al., 2010a, b), and comparatively so across the different aquifers.
Present-day personal computers are capable of solving the coupled non-linear equations governing variable-density flow and transport. Indeed, this is often done for specific studies at particular sites, like those of the present aquifer cases. However, the high computational demands of detailed simulations are an obstacle, and the usefulness of such simulations has also been questioned on the grounds that their results often lack reliability due to the common paucity of data for adequately characterizing the detailed hydrogeology (hydraulic parameters) and concentration field (Sanford and Pope, 2010), and particularly the aquifer's dispersive properties. In addition, the boundary conditions and hydrologic forcing (recharge and pumping) are often not sufficiently known over time. Thus, in the absence of credible data, the benefits of using detailed and complex numerical variabledensity models are undermined by the guesswork (or nonunique fitting) and simplifications needed for the estimation of multiple aquifer parameters and their variability in space and time.
In comparison, analytical solutions of steady-state sharpinterface flow of one flowing fluid are evaluated with efficiency and are suitable tools for screening-level, regionalscale assessments of coastal aquifer vulnerability to seawater intrusion. They are also useful for screening the impact of various groundwater pumping alternatives in management scenarios (Koussis et al., 2012;Mantoglou, 2003). After developing a sound conceptual model, properly schematizing aquifer geometry, and carefully estimating dominant hydraulic parameters and forcing quantities, analytical solutions can be used for systematic analyses, elucidating the way in which various physical and management parameters impact the seawater intrusion behavior of coastal aquifers, and determining likely intrusion responses to various forcing conditions (Werner and Simmons, 2009;Ferguson and Gleason, 2012;Mazi et al., 2013).
From this perspective, we use here the generalized analytical model of Koussis et al. (2012) (see Appendix A for the mathematical details) for steady interface flow in sloping phreatic aquifers, recently enhanced to allow for a gap in the sharp interface, through which submarine groundwater discharges (Koussis et al., 2014). That enhanced model retains the simplicity of classical Dupuit-type sharp-interface solutions, which however ignore that gap, assuming outflow through an idealized singular point. The solution of Koussis et al. (2012) accounts for the generally present and usually hydraulically significant aquifer slope that has previously been ignored in analytical sharp-interface solutions of seawater intrusion (e.g., Strack, 1976). This discharge-potential model approximates the gravity-driven flow component and represents the aquifer geometry in a schematized yet realis-tic manner. To compare our results to those obtained with variable-density models, we assigned to the sharp interface the 50 % salinity line (17 500 ppm) of the transition zone, as suggested by Reilly and Goodman (1987) and generally adopted. With the use of this generalized model and approach, Koussis et al. (2012) have identified important criteria useful in the management of an aquifer, and Mazi et al. (2013) have shown the aquifer's highly non-linear response to seawater intrusion and associated tipping points that should not be crossed in order to avoid rapid loss of control and complete seawater intrusion.
The model considers an inclined (sinϕ) unconfined homogeneous aquifer with depth H sea at the coast, with length L and hydraulic conductivity K. The aquifer is further recharged uniformly at rate r, in addition to receiving an inflow, q b (for prescribed-flux boundary condition, FB), or having a constant-head boundary, h L , that determines inflow through its land boundary (for prescribed-head boundary condition, HB). Furthermore, in the case of inland FB condition, the model solution can account for groundwater pumping in a collector trough (line sink, or well gallery), at a distinct distance l w from the coastline, that penetrates the aquifer completely and draws groundwater at the rate q w . In the case of inland HB condition, distributed groundwater pumping can be accounted for by reducing the recharge rate r with the pumped groundwater amount. Both types of boundary condition and their respective solutions enable calculation of the location of the seawater intrusion interface toe, l T , and the important (Destouni et al., 2008) but hard to determine in the field (Prieto and Destouni, 2011) submarine discharge, q SD . See further Fig. A1 in Appendix A for more details on all key parameters, boundary conditions and solutions discussed here.

Analysis approach
For the present analysis of aquifer cases with inland FB condition, and with a distinct groundwater pumping location, l w , and pumping rate q w , we introduce the term q norm = −(q o +q w )/(KL), which is the normalized remaining groundwater flow from l > l w after the pumping of q w , where q o = −r(L− l w ) + q b is the total groundwater flow from the area between the land boundary and the pumping location l w (and with q b , q o and all other groundwater flow terms defined as negative in the seaward direction). We note that, due to considering a vertical plane (see Appendix A), all flows are per unit width perpendicular to the aquifer plane (or parallel to the coastline), i.e., the flow units are m 3 m −1 per unit time (day or year). In order to make calculation results for different such case studies directly comparable with each other, they are shown and discussed in the following section in terms of the normalized position of the sharp interface toe, l T /L, as a function of q norm . For the analysis of aquifer cases with inland HB condition, and distributed groundwater pumping, resulting l T /L is instead shown as a function of the consistently normalized submarine discharge q SD /(KL), with q SD = q b -rL + q w = q b -r n L, where r n = (r + q w /L) is apparent net aquifer recharge after distributed pumping of groundwater at rate q w . The resulting curves of l T /L versus q SD /(KL) or q norm are then referred to as the performance curves of seawater intrusion (Koussis et al., 2012).
For each case study, we have further selected for the analysis typical cross sections that are parallel to the groundwater flow in the different aquifers. For the NDA, where the flow is radial, we have selected two such cross sections, one along the middle Nile Delta aquifer (mNDA), in the north-south direction, and one in the east Nile Delta aquifer (eNDA), following a northeast direction (Fig. 2a). For ICA, the cross section has an east-west orientation (Fig. 2b), while in CAA the orientation is northwest to southeast (Fig. 2c). Parameters used in the cross-section modeling of each case study are summarized in Table 1, and the cross-section conceptualizations are shown in Fig. 3 (a, b for NDA, c for ICA, d for CAA).
For the NDA, the lengths of the two cross sections are as follows: the mNDA section has a length of 175 km and the eNDA section has a length of 180 km (Fig. 3a, b). Using the surface of the Mediterranean Sea as datum, the piezometric head at the inland boundary is at 14 m for both sections. The hydraulic conductivities of the two sections were found by interpolating the curves in Fig. 7a, b of Sherif et al. (2012); for the mNDA profile section K = 100 m d −1 , and for the eNDA profile section K = 120 m d −1 . The slope of the aquifer base and the aquifer depth at the coast were determined for the mNDA profile-section as 0.3 % and 740 m, respectively, and for the eNDA profile-section as 0.4 % and 910 m, respectively (Fig. 3a, b).
In their recent paper, Sherif et al. (2012) performed simulations with FEFLOW for five horizontal layers (areal sections) of the NDA, reaching to a depth of 400 m; without explicitly reporting the aquifer recharge, they report the seawater intrusion to exceed 100 km (1000 ppm TDS) from the coast. As no information is openly available for the NDA recharge, and in order to use the latest published data on the Nile Delta seawater intrusion, we have used the present analytical model to calculate the likely current location of the interface toe, l T , by calibrating the net aquifer recharge, r n , in the definition of submarine groundwater discharge q SD . Through this calibration, the sharp interface toe for the mNDA was positioned at l T = 100 km, according to the intersection of the constant slope aquifer base with the 50 % salinity line of profile 1 in Fig. 10 of Sherif et al. (2012) (the 50 % salinity line was interpolated from the 1000 and 35 000 ppm lines in the salinity transition zone). The thus calibrated net recharge value r n was found to be about 10 mm yr −1 , incorporating the aquifer recharge from all possible (natural and artificial) sources minus the known magnitude of groundwater extractions from the mNDA. Because the extraction locations are distributed in the NDA, there is no single representative pumping location, l w , that we could use in the   (km) 143.5 153.6 20 Submarine discharge, q SD (m 3 m −1 day −1 ) 3.9 4.2 4.9 0.9 L/H sea : aquifer length to depth at the sea 237 198 100 60 Normalised difference between the densities of 1/40 1/40 1/40 1/40 sea-and freshwater ρ s and ρ f , respectively: δ = (ρ s -ρ f )/ρ f * Water withdrawals have been subtracted from the aquifer recharge rate, to yield net recharge rate r n . modeling and the illustration of resulting intrusion performance curves for this aquifer case; performance curves for the NDA sections are therefore shown in terms of l T /L as a function of q SD /(KL). Table 1 further summarizes the modeling parameters for the NDA sections, with r n for the eNDA being assumed approximately equal to that calibrated for the mNDA based on results reported by Sherif et al. (2012).
For the distributed present pumping in the NDA cross sections, we then calculated first the interface toe position l T for various exploitation conditions of the aquifer (at various times in the past), assuming that the head at the inland boundary remained unaltered (at today's levels) at 14 m. We also performed a sensitivity analysis of the combined effects of calculated net recharge changes from former to current pumping rates, in conjunction with possible groundwater head changes at the inland boundary. The latter could be brought about by changing management conditions in the Nile River valley aquifer, or, more generally, by a change in the apportionment of the Nile waters among neighboring countries, and/or by climate change.
For the ICA, our cross-section conceptualization is based approximately on the work of Yechieli et al. (2010), but considers a 20 km-long cross section ending in a perfect cliff, with an aquifer depth at the coastline of 200 m, slope of the impermeable aquifer base of 1 %, hydraulic conductivity K = 30 m d −1 and no boundary inflows (Fig. 3c). The total aquifer recharge considered here is 240 mm yr −1 , with 200 mm yr −1 derived from natural recharge and 40 mm yr −1 from artificial recharge. We further consider a pumping rate of 285 × 10 6 m 3 yr −1 (3000 m 3 m −1 yr −1 ) through a fully penetrating trough at 3000 m from the coast. Our analytical modeling applied to current conditions yields then the sharp-interface toe position at 2.6 km, in agreement with the intrusion stated by Aberbach and Sellinger (1967). The complete set of the parameters used for the ICA modeling is further summarized in Table 1.
For the CAA, the one-layer cross-sectional conceptualization is shown in Fig. 3d. The length of the aquifer is L = 3000 m, its hydraulic conductivity K = 28 m day −1 , the slope of the aquifer base is sinϕ = 0.017, and the depth of the sea to the aquifer base at the coastline is H sea = 50 m; a representative pumped collector trough is located at l w = 1000 m. Conditions are typical of the eastern area of the unconfined aquifer in the Akrotiri peninsula, the Zakaki area bordering the city of Lemessos. The natural recharge, r, and the inflows through the aquifer's landside boundary, q b , were estimated by Mazi et al. (2004) from the hydrological balance of the Akrotiri Basin for the period 1972-1992 at the following mean values: r = 92 mm yr −1 and q b = −549 m 2 yr −1 (m 3 per meter width). For the same period, Koussis et al. (2012) translated the annual groundwater withdrawals in the Zakaki area of 2 × 10 6 m 3 to the 2-D-equivalent pumping rate at a representative well gallery of q w = 500 m 3 m −1 yr −1 for use in the profile model of the aquifer. The full parameter set used in the CAA modeling is further summarized in Table 1.
For both the ICA and CAA cross sections, with their more well-defined average pumping locations, l w , we have calculated their respective intrusion performance curves in terms of l T /L as a function of q norm = −(q o + q w )/(KL), in other words, the normalized remaining groundwater flow from l > l w after the pumping of q w . We further analyzed the effects of different possible pumping locations on the intrusion performance curves. The results of these calculations are illustrated in Fig. 6 below.

Model result interpretation of critical intrusion points
In the ICA and CAA cases with inland FB condition (prescribed boundary inflow q b ) and distinct pumping location l w , the associated analytical solution allows for evaluating whether a hydraulic divide (maximum hydraulic head) is formed at a location l div between the well and the coastline, 0 < l div < l w . The ratio l T /l div , which expresses the relation between the interface toe location and the divide location, then depends on and is inversely proportional to the submarine discharge, q SD = q b -rL + q w . As pumping increases, q SD decreases and l div moves seawards, that is, decreases. If a hydraulic divide is formed in the region 0 < l < l w (i.e., if l div /l w ≤ 1) then the pumping will also draw groundwater from the region (l wl div ) between the coast and l w . A critical tipping point of intruding seawater that is just on the verge of invading the whole coastal aquifer will be reached here if the interface toe reaches the groundwater divide (i.e., if l T /l div = 1), as identified by Mazi et al. (2013). If no groundwater divide is formed between the coastline and the pumping location (i.e., if l div /l w ≥ 1) then the condition l T /l w = 1 constitutes another point of critical seawater intrusion. This condition means that the invading seawater reaches the distinct well location, l w , so that the well is lost before the intruding seawater reaches the tipping point of complete aquifer intrusion. The critical intrusion point l T /l w = 1 is reached if the total groundwater flow q o = −r(Ll w ) + q b from the area between the land boundary and the pumping location l w is too small to counteract seawater intrusion to l w .
The pumping of q w may further draw groundwater only from l > l w , or from both l > l w and from the region 0 < l < l w . When the value of the normalized remaining groundwater flow after pumping, q norm = −(q o + q w )/(KL), is positive, the pumped water is only from the region between the well and the land boundary; when q norm is negative, some of the pumped water is from the region between the well and the coastline. The flow q norm cannot become negative if a divide is not formed in the region 0 < l < l w . If the pumping location is moved further inland, the pumping of q w may draw water from both sides of the well gallery, such that a divide can form between the sea and l w if sufficient submarine groundwater discharge q SD = q b −rL+ q w remains after that pumping. However, just moving a pumping well farther inland does not by itself ensure that the well will not be invaded by intruding seawater; such protection can only be achieved if a hydraulic divide indeed is formed between the well location and the coastline, and if the tipping point l T /l div = 1 of complete aquifer invasion is not reached by the toe extending to that divide.
In the case of the inland HB condition (prescribed head at the boundary), the solution allows for evaluating whether a hydraulic divide is formed at a location l div < L (the location of the inland boundary). If a hydraulic divide exists in the aquifer and the net recharge r n decreases (e.g., due to increased distributed groundwater pumping), the submarine groundwater discharge q SD will decrease and the groundwater divide location l div will move landward, such that it may then reach the end of the aquifer, so that l div /L ≥ 1. This means that a tipping point of the type l T /l div = 1 for complete aquifer intrusion will not be reached in this case. However, if the hydraulic head at the inland boundary decreases (e.g., due to groundwater pumping, or climatedriven recharge decrease in the hydrological catchment upgradient of the inland aquifer boundary), the divide location l div will move seaward, allowing for the interface toe of seawater intrusion to possibly meet it, and the critical tipping point of l T /l div = 1 to be reached such that complete seawater intrusion occurs.
In both cases of possible inland boundary conditions (FB and HB) there exists thus the possibility of reaching critical tipping points of complete seawater intrusion in the coastal aquifer l T /l div = 1. In the flux boundary (FB) case, with distinct pumping location l w , the other type of critical seawater intrusion point, in which the intruded wedge reaches the well gallery (i.e., l T /l w = 1) may also occur. Which of the two types of critical seawater intrusion points may be reached first, and thus be limiting for aquifer resilience to seawater intrusion, depends on actual aquifer conditions and their changes, as quantified further in the following section.

Aquifer-specific calculation results
According to our calculations, prior to the 1950s, when pumping from the NDA was negligible, the interface toe location was at 14 km in the mNDA and at 23 km in the eNDA. In the 1990s (Sherif and Al-Rashed, 2001) the pumping from the aquifer was 1.92 billion m 3 yr −1 and, according to our calculations, the toe moved then to 26 km in the mNDA and to 40 km in the eNDA. Increased pumping after the year 2000 to 2.4 billion m 3 yr −1 (Mikhailova, 2001) caused the interface to move to its current position, that is, 100 km from the coast in mNDA and 120 km from the coast in eNDA. The strong non-linearity of the intrusion performance curves implies that the increased pumping has brought the NDA to a marginal situation (Fig. 4).
More specifically, the submarine discharge q SD is currently 4.2 m 3 m −1 day −1 for the eNDA, and 3.9 m 3 m −1 day −1 for the mNDA, with the associated groundwater divide and interface toe locations then being at ∼ 154 km and 116 km, respectively, in the eNDA, and at ∼ 143 km and 98 km, respectively, in the mNDA. Thus, the ratio of the toe to the divide location, l T /l div , which is a key parameter for aquifer resilience/vulnerability to seawater intrusion, is 0.755 for the eNDA and 0.684 for the mNDA, which explains the higher vulnerability of eNDA. As r n decreases, the hydraulic divide will reach the end of the aquifer, l div = L, for r n = 6.5 mm yr −1 in the eNDA and for r n = 5.7 mm yr −1 in the mNDA. The ratio l T /l div will then be 0.69 and 0.64 in the eNDA and mNDA, respectively, showing that l div moves faster than l T under decreasing r n . Because of the landside HB condition, with piezometric head at 14 m, the aquifer will not be completely invaded by seawater even if net recharge r n diminishes to zero; under those conditions l T ≈ 140 km in the eNDA and l T ≈ 131 km in the mNDA (Eq. 2, Appendix A). In other words, even for r n = 0, the interface toe l T cannot here intrude to the inland boundary at L = 180 km in the eNDA and L = 175 km in the mNDA, if the prescribed hydraulic head and the implied inflow there are sustained at current levels (Fig. 4). However, if water management practices and/or climate change in the hydrological catchment upstream of inland aquifer boundary decreased the hydraulic head at that boundary, the tipping point l T /l div = 1 of full seawater intrusion could be reached, even if the net recharge were maintained at current levels (inset, Fig. 4).
The sensitivity analysis of the combined effect of recharge and boundary head change for the eNDA and mNDA shows that, if net aquifer recharge and associated submarine groundwater discharge q SD had remained at their former greater values, the head at the inland boundary would not play such a big role for the vulnerability of the aquifer to seawater intrusion. Specifically, a change of ±10 % in the boundary head would only cause the interface toe to fluctuate by ±350 m at the most in the mNDA section and by  −860 m to +420 m in the eNDA section ( Fig. 5a and Table 2). However, under current conditions of net recharge limited to about 10 mm yr −1 , a +10 % head change causes the interface to retreat by 6.5 km in both sections, while a −10 % head change causes the interface to advance by 5.3 km in the mNDA and by 7.7 km in the eNDA ( Fig. 5b and Table 2). In particular, a decrease in the inland boundary head of 22 % in the mNDA and of 17 % in the eNDA can bring the NDA to its tipping point of complete aquifer intrusion (Fig. 4, inset). Therefore, under current recharge conditions, it is critical for the resilience and sustainability of the Nile Delta aquifer to control the up-gradient catchment hydrological conditions so that sufficient piezometric groundwater head is maintained at the inland aquifer boundary. For the ICA and CAA cross sections, resulting seawater intrusion performance curves are shown in terms of normalized intrusion toe position l T /L as a function of the normalized remaining groundwater flow q norm = −(q o + q w )/(KL) (Fig. 6). In general, depending on which condition is the limiting one for each case, the resulting performance curves terminate either at the tipping point of complete aquifer intrusion l T /l div = 1, or at the limit of well intrusion l T /l w = 1.
If the tipping point condition l T /l div = 1 is the limiting one, Fig. 6 shows also the corresponding value of l T /l w < 1. If the well intrusion condition l T /l w = 1 is instead limiting, the performance curves are continued (simulating the condition of pumping seawater) as dotted lines from that limiting point until the tipping point of complete seawater intrusion in the aquifer.
For current pumping rates (corresponding to the q norm = −(q o + q w )/(KL) value of the red dashed lines, Fig. 6; with q w = 3000 m 3 m −1 yr −1 in ICA and q w = 500 m 3 m −1 yr −1 in CAA) and pumping locations (black curves; l w /L = 0.15 (l w = 3 km) in ICA, Fig. 6a; l w /L = 0.33 (l w = 1 km) in CAA, Fig. 6b), the well intrusion condition l T /l w = 1 (black filled circle) is the limiting one for both aquifers. The ICA (Fig. 6a) is then closer to that limit than the CAA (Fig. 6b). Increasing the pumping rate q w at the same well location will decrease q norm = −(q o + q w )/(KL) and move the aquifer closer toward the well intrusion limit l T /l w = 1 (along the black curves, Fig. 6). The risk of crossing that limit can be reduced by moving the pumping location further inland (shifting from the black to, e.g., the green or blue curves, Fig. 6) but the tipping point of complete aquifer intrusion l T /l div = 1 (open square symbol, Fig. 6) may then become the limiting condition instead of the well intrusion limit. Table 3 summarizes for the different well location choices illustrated in Fig. 6 the highest possible groundwater pumping rate, (q w ) max (and other associated variable values), until a critical seawater intrusion limit is crossed.
In the ICA simulations for pumping locations l w = 6000 m and l w = 8000 m, the maximum pumping rate (q w ) max = 3603 m 3 m −1 yr −1 is limited by the condition l T = l div = 4987 m. In the CAA simulations for pumping locations l w = 1400 m and l w = 2000 m, the maximum pumping rate (q w ) max = 697.6 m 3 m −1 yr −1 is limited by the condition l T = l div = 1385 m. We note then that moving pumping locations more "up-gradient" does not, without further specification, automatically imply well protection from seawater intrusion and increase of exploitable water volume. What matters is the specific location of the hydraulic divide to be between the coastline and the pumping location (critical point of seawater intrusion to well) and the interface toe not to reach the divide (tipping point of complete intrusion). In the specific ICA and CAA cases, there is then no gain at all, in either exploitable water volume or aquifer protection, of moving pumping locations further inland than l w = 6000 m in ICA and l w = 1400 m in CAA.

Conclusions
Directly comparable performance curves have here been developed for analyzing the resilience/vulnerability of different coastal aquifers to seawater intrusion. Two different types of critical seawater intrusion limits are then identified: the Table 3. Variable quantification related to the results of Fig. 6 for the Cyprus Akrotiri aquifer and the Israel Coastal aquifer. Quantified variables include groundwater exploitation (current pumping rate q w , possible maximum pumping rate (q w ) max ), pumping l w in relation to aquifer length L, intrusion toe position l T /L and (l T /L) max for q w and (q w ) max pumping, respectively, and relative location l div /L of a possible groundwater divide formed between the well and the coastline.
Curve, Position q w (q w ) max Limiting Fig. 6 of well, l w l w /L l w /H sea (m 3 m −1 yr −1 ) l T /L l div /L (m 3 m −1 yr −1 ) (l T /L) max Criteria well intrusion limit l T /l w = 1, and the complete intrusion limit l T /l div = 1. The risk of crossing the first limit can be reduced by moving pumping locations further inland, but then the tipping point of complete aquifer intrusion may be crossed instead. The measure of moving pumping locations more inland beyond the limiting location of the prevailing groundwater divide in a coastal aquifer will offer no gain in more exploitable water volume or improved aquifer protection from seawater intrusion. The site-specific vulnerability assessments show that the advance of seawater currently seriously threatens the Nile Delta aquifer (NDA). The inland boundary head must here be sustained at any cost, as even a 10 % decline will cause seawater intrusion advancement by around 5-8 km and a 20 % decline may bring the NDA to its tipping point of complete aquifer intrusion. These indications call for careful further study for establishing reliable hydrologic conditions for the NDA.
In the Israel Coastal aquifer (ICA), the current pumping location does not allow for increased groundwater abstractions. The maximum pumping could here be increased by about 20 % if the pumping location were moved more inland, at least at l w = 5 km from the coast, thereby allowing a groundwater divide to form between the well and the coastline. In the Cyprus Akrotiri aquifer (CAA), current forcing conditions do not directly threaten groundwater sustainability by increased seawater intrusion. Maximum pumping could here be raised by up to around 28 % above the present rate. If pumping locations here were moved more inland (to about l w = 1.4 km from the coast), the maximum pumping rate would increase by up to 40 % over the present rate. In aquifers with flux-boundary and pumping from a welldefined location l w , if a hydraulic divide is formed between the interface toe and the pumping location, such that l T ≤ l div ≤ l w , the divide location can be calculated from while the head at the divide is calculated from with the physically meaningful solution being the root with the positive square root.
In aquifers with prescribed-head boundary and without a distinct pumping location (but still possible distributed pumping that may be represented by reduced recharge r), if a hydraulic divide is formed inside the aquifer, its location can be calculated (Mazi et al., 2013) from where z b denotes bed elevation and the hydraulic head is given by Finally, for the case of a shallow, horizontal coastal aquifer, under exploitation and negative recharge (evaporation), we refer the reader to the analytical solution of Kacimov et al. (2009), who also demonstrated that, under extreme conditions, even a gap can form separating the saturated fresh and saline groundwater zones.