# Coastal aquifer recharge and groundwater-seawater exchanges using downscaled GRACE data: case study of the Djeffara plain (Libya–Tunisia)

Julio Gonçalvès, Pierre Séraphin, Thomas Stieglitz, Amine Chekireb, Bruno Hamelin, Pierre Deschamps

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Julio Gonçalvès, Pierre Séraphin, Thomas Stieglitz, Amine Chekireb, Bruno Hamelin, et al.. Coastal aquifer recharge and groundwater-seawater exchanges using downscaled GRACE data: case study of the Djeffara plain (Libya–Tunisia). Comptes Rendus. Géoscience, 2021, 353 (1), pp.297-318. 10.5802/crgeos.74. hal-03361219

HAL Id: hal-03361219

https://amu.hal.science/hal-03361219

Submitted on 20 Oct 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Géoscience
### Sciences de la Planète

Julio Gonçalvès, Pierre Séraphin, Thomas Stieglitz, Amine Chekireb, Bruno Hamelin and Pierre Deschamps

**Coastal aquifer recharge and groundwater-seawater exchanges using downscaled GRACE data: case study of the Djeffara plain (Libya-Tunisia)**

Volume 353, issue 1 (2021), p. 297-318

<https://doi.org/10.5802/crgeos.74>

© Académie des sciences, Paris and the authors, 2021.

*Some rights reserved.*

<figure><img src="image_3.png" /></figure>

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.
http://creativecommons.org/licenses/by/4.0/
Original Article — Hydrology, Environment

# Coastal aquifer recharge and groundwater-seawater exchanges using downscaled GRACE data: case study of the Djeffara plain (Libya–Tunisia)

Julio Gonçalvès<sup>a, *</sup>, <sup>a</sup> Pierre Séraphin<sup>a</sup>, Thomas Stieglitz<sup>a</sup>, <sup>a</sup> Amine Chekireb<sup>a</sup>, Bruno Hamelin<sup>a</sup> and Pierre Deschamps<sup>b, a</sup>

<sup>a</sup> Aix Marseille Univ, CNRS, IRD, INRAE, Coll France, CEREGE, Aix-en-Provence, France

*E-mails*: goncalves@cerege.fr (J. Gonçalvès), seraphin@cerege.fr (P. Séraphin), stieglitz@cerege.fr (T. Stieglitz), chekireb@cerege.fr (A. Chekireb), hamelin@cerege.fr (B. Hamelin), deschamps@cerege.fr (P. Deschamps)

**Abstract.** The Djeffara plain, extending over Tunisia and Libya (48,000 km²), experiences unsustainable exploitation of groundwater, the severity of which needs to be determined. Here, we used data from Gravity Recovery and Climate Experiment (GRACE) and Global Land Data Assimilation System (GLDAS) to assess the groundwater storage variation (ΔGWS) for the period 2002–2016. ΔGWS values, downscaled over a 0.5° × 0.5° grid using pumping wells spatial distribution and rates, were introduced in a discrete set of groundwater and salt budget equations to ascertain the recharge and the coastal fluxes. An average recharge of 4.3 ± 2.7 mm·yr⁻¹ representing 2.2 ± 1.4% of the domain-averaged annual rainfall was obtained. A saline intrusion of 147.3 × 10⁶ and 46.9 × 10⁶ m³·yr⁻¹ was calculated for the Tripoli (Libya) and Djerba (Tunisia) areas, respectively. The overall Submarine Groundwater Discharge (SGD) of freshwater along the 450 km long Djeffara coast represents about 339.6 × 10⁶ m³·yr⁻¹ for the 2002–2016 period.

**Keywords.** Groundwater mass balance, Satellite-derived gravity, Recharge, Submarine groundwater discharge, Arid domain, Groundwater sustainability.

Manuscript received 25th March 2021, revised 9th June 2021 and 1st July 2021, accepted 7th July 2021.

## 1. Introduction

In the Djeffara plain, like in other areas of the Maghreb, demographic and economic growths as well as the development of irrigated agriculture contribute to massive and unsustainable groundwater exploitation. In the transboundary multilayer aquifers system underlying the Djeffara plain shared

by Tunisia and Libya, withdrawals increased from 0.03 to 1.0 km³·yr⁻¹ between 1950 and 2011 [Alfarrah and Walraevens, 2018, OSS, 2006]. About 90% of the water extracted from these aquifers is used primarily for irrigated agriculture. The very intensive mining of these resources since 1950 has caused a decline of groundwater reserves together with the drying up of naturally flowing springs which were historically exploited [OSS, 2006]. An accurate assessment of the recharge, whether natural or artificial, is pivotal

<footnote>* Corresponding author.</footnote>
to determine the degree of (un)sustainability of
groundwater exploitation and to guide management
policies. The few available values for the recharge in
the Djeffara plain are derived from hydrogeological
modeling [OSS, 2006]. With the exception of the stud-
ies by Allemmoz and Olive [1980] and Zouari et al.
[1985] in the regions of Tripoli and Gabès (see Fig-
ure 1), alternative approaches using geochemical or
hydrogeophysical methods, for instance, are missing
so far. While the estimation of a still unclear current
recharge is a common issue for all the circum-Sahara
regions, coastal aquifers such as the Djeffara system
raise specific problems. Heavy and increasing with-
drawals in coastal domains lead to salinization prob-
lems caused by seawater intrusion [Steyl and Dennis,
2010]. This has been clearly established in the Tripoli
area (Libya) with saline intrusion due to increasing
urbanization and agricultural activities [Alfarrah and
Walraevens, 2018, Sadeg and Karahanoglu, 2001].
The identification of groundwater-seawater fluxes,
termed Submarine Groundwater Discharge (SGD), is
another issue specific to coastal aquifers. The quan-
titative importance of SGD in comparison to surface
water fluxes for oceans mass balance is a subject of
controversy but their qualitative impact is largely
accepted [e.g. Luijendijk et al., 2020]. SGD is often
associated with a substantial supply of nutrients that
affects the ecological equilibrium of coastal areas
[e.g. Rodellas et al., 2015].

A commonly used approach to derive SGD fluxes
on local to regional scales is to close a radionu-
clide mass balance, using predominantly radon and
radium as SGD tracers. Local SGD fluxes (on a
bay/beach scale) were estimated using these trac-
ers at many locations along the northern (Euro-
pean) coast of the Mediterranean Sea, but there
is a conspicuous lack of any data-based SGD esti-
mates in the southern Mediterranean basin along
north African coastlines, including those adjacent to
the Djeffara plain. Equally important, such radionu-
clide tracer-derived estimates include all advective
fluxes across the sediment water interface in the term
SGD [e.g. Rodellas et al., 2015]. It has been widely
demonstrated that the major process contributing to
radionuclide fluxes (>80% in most cases) is the recir-
culation of seawater through coastal sediments, and
freshwater SGD generally are one order of magni-
tude lower than such total SGD fluxes [e.g. Stieglitz
et al., 2013]. While radionuclide-derived values are

useful from a geochemical perspective, they do not
yield estimates comparable to freshwater flow esti-
mates based on hydrogeological or geophysical ap-
proaches, unless combined with other tracers such
as salinity [e.g. Stieglitz et al., 2013, Weinstein et al.,
2007]. Similar to the quantification of recharge, no es-
timates other than those produced by hydrogeolog-
ical modeling are available for freshwater SGD [e.g.
Zektser et al., 2007] or saline intrusion.

Since 2002, a satellite monitoring of Earth's grav-
ity field (GRACE experiment) enables groundwater
mass storage (GWS) variation assessments over large
regional-scale aquifer systems. Satellite approaches
provide an integrated view of processes over large
areas, which would be much more time- and cost-
consuming to obtain from ground observations, if
not impossible. GRACE data have been widely used
for quantifying aquifer storage changes [e.g., Ahmed
et al., 2014, Bonsor et al., 2018, Carvalho Resende
et al., 2019, Leblanc et al., 2009, Ramillien et al.,
2014, Richey et al., 2015, Rodell et al., 2009]. Some
studies used a global mass balance to estimate the
domain-averaged recharge for large regional aquifers
of the Saharan belt [Gonçalvès et al., 2013, Mohamed
et al., 2017]. Similar to these studies, GRACE-derived
groundwater mass variations (ΔGWS) could be com-
bined with ground data to estimate hydrological vari-
ables such as the recharge or the coastal fluxes in-
dependently from any hydrogeological model. How-
ever, all these previous applications using GRACE
data have been carried out on large surface basins
(>100,000 km²). In fact, due to the relatively low
resolution of the GRACE data (>1°) and signal-to-
noise ratio issues, poor quality applications are ex-
pected for basins with a surface area of less than
about 100,000 km². One way to address this prob-
lem for more local-scale basins such as the Djeffara
aquifers system is to downscale the GRACE data us-
ing ground-based information or data assimilation
systems variables [Miro and Famiglietti, 2018, Sey-
oum et al., 2019].

In this context, GRACE-derived ΔGWS for the
Djeffara domain are first downscaled at the res-
olution of the gridded annual average rainfall
(AAR) data (0.5° × 0.5°) assuming a major control
of pumping wells on ΔGWS. The weights in the
downscaling approach proposed here are therefore
based on the distance to pumping well areas and
pumping rates. Downscaled GRACE-derived ΔGWS
<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) Study area showing the Djeffara plain geographical extension, and (b) Schematic view of the geological context reproduced from Ben Dhia [1987].</figcaption></figure>

values are compared to their independent calcu-
lations of groundwater mass variations based on
ground data, i.e., aquifer storativity and water level
drawdown values available for the Libyan Djeffara.

The classical groundwater mass balance equation
is discretized over the 0.5° × 0.5° grid introducing
the fluxes to be ascertained, i.e., the recharge and
the coastal fluxes as well as the downscaled ΔGWS
values. Then, the linear system of discrete mass balance equations obtained was solved for recharge and saline intrusion or SGD assuming a linear relation between average annual rainfall (AAR) and recharge. The results, which are the first independent estimates from any numerical models to our knowledge, are discussed in comparison with hydrogeological modeling outputs as well as with SGD flux estimates based on radium tracer approaches.

## 2. Material

### 2.1. Hydrogeological context

Located between the Mediterranean coast and the mountains of the Nafusa mountains to the south and the Dahar to the west, the Djeffara plain covers a surface area of about 48,000 km² over eastern Tunisia and northern Libya (Figure 1a). This region receives an annual average precipitation of about 200 mm for a potential evaporation of 2760 mm characterizing an arid climate [Trabelsi et al., 2009]. The mountains (jebels) correspond to anticlines that structurally separate the Djeffara Plain from the north-western Sahara Aquifers System (NWSAS, "Continental Intercalaire" CI and "Complexe Terminal" CT aquifers Figure 1b). The Djeffara basin is characterized by an almost continuous stratigraphic succession from the Permo-Triassic to the Quaternary [OSS, 2006]. However, in the southern Libyan part, the Triassic layers are unconformably overlain by Miocene formations (erosion surface). Therefore, from a (hydro)geological standpoint, the Djeffara plain can be subdivided into a northern domain (Tunisia) where the aquifer formations are mainly the Senonian and Mio-Pliо-Quaternary units, and a southern domain where the Triassic and Miocene formations are of primary importance in addition to the Plio-Quaternary deposits. Overall, the stratigraphic succession can be closely examined to define a multilayered aquifer system with lateral variations from the Tunisian to the Libyan Djeffara. In Tunisia, the Trias is composed of sandstones, dolomites, and, in the upper part, an intercalation of clay and evaporites which extends stratigraphically until the Jurassic. The Trias and Jurassic formations are overlain by the dolomitic Turonian and the calcareous Senonian formations which constitute the main exploited and productive

aquifer together with the Mio-pliocene-Quaternary sands. The Mio-pliocene-Quaternary aquifer is separated from the Senonian by more or less continuous marly gypsum and clayey formations. In Libya, the Lower Triassic consists of sandy clay then dolomite, marls, and clay layers overlain by sandstone. The Upper Triassic and Lower Jurassic are made of gypsum and anhydrite. The Jurassic and Cretaceous formations (only along Nafusa Mountains) are of small hydrogeological importance either because they are absent (erosion) or because they are not really exploited [Ben Lasmar et al., 2017]. From bottom to top, the Miocene begins with a clay layer followed by a limestone formation (sometimes referred to as the Lower Miocene Aquifer), then a sandy limestone (Middle Miocene Aquifer). Separated from the latter by a sandy clayey layer, the Mio-Plio-Quaternary is extremely heterogeneous (sand, alluvial deposits, limestones, clay) and intensively exploited in terms of groundwater.

The hydrostratigraphic analysis conducted by OSS [2006] leads to the definition of three aquifers which are in fact multilayer aquifers. The Upper Miocene, Pliocene, and Quaternary deposits form a continuous aquifer over the entire Djeffara (known as the phreatic or Upper aquifer). The Lower and Middle Miocene, also associated with the Jurassic and Upper Triassic when confining units are absent, form the Middle aquifer (deep coastal aquifer). The deep aquifer comprises the Lower and Middle Triassic formations.

Groundwater flow directions are generally from the reliefs of the jebels Nafusa and Dahar to the coast and therefore mostly from west to east in the Tunisian part and almost south to north in the Libyan Djeffara. All aquifers are characterized by high transmissivities between $1 \times 10^{-2}$ and $5 \times 10^{-2} \text{ m}^2 \cdot \text{s}^{-1}$ and the regional hydraulic head gradient is in the order of 0.1% [OSS, 2006]. Interestingly, introducing these values in a Darcy flow calculation for a coastal length of 450 km leads to a horizontal drainage of 4.5 to $22.5 \text{ m}^3 \cdot \text{s}^{-1}$ using the above transmissivity range which can be translated into a steady-state recharge of between 3.0 and $14.8 \text{ mm} \cdot \text{yr}^{-1}$. Recharge occurs locally, either focused or diffuse as well as from inflows at the boundaries. Recharge estimates (natural and irrigation return flow) were obtained from modeling approaches at $6.8 \text{ mm} \cdot \text{yr}^{-1}$ [$1.8 \text{ mm} \cdot \text{yr}^{-1}$ for focused recharge; OSS, 2006]. Using a water table
fluctuation method and soil moisture profiles, Allem-
moz and Olive [1980] found a recharge of 10 mm·yr⁻¹
in the area of Tripoli (Libya) exclusively from focused
recharge (wadi flood infiltration). Zouari et al. [1985]
used chloride and stable isotope of water profiles in
the vadose zone near Gabès (Tunisia) to estimate in-
filtration (5 mm·yr⁻¹) and evaporation (1 mm·yr⁻¹,
presumably from the water table) leading to a net dif-
fuse recharge of about 4 mm·yr⁻¹.

Groundwater inputs at the boundaries of the Djef-
fara aquifer system correspond to (i) a horizontal flow
from the south in the Trias at the Nafusa [OSS, 2006]
and (ii) an upward flow along faults which is one of
the main outflows of the deep confined CI aquifer of
the NWSAS [Gonçalvès et al., 2015]. Heavy ground-
water exploitation (see Section 2.3) causes continu-
ous seawater intrusion and salinization of groundwa-
ter. This was particularly reported in the most pop-
ulated Tripoli area [Alfarrah and Walraevens, 2018,
Sadeg and Karahanoglu, 2001]. Groundwater salin-
ization by dissolution of evaporites is also docu-
mented [Trabelsi et al., 2012].

## 2.2. GRACE and GLDAS data

The GRACE twin-satellite system launched in 2002
by NASA and the German Aerospace Center (DLR)
monitors the earth gravity on a monthly basis [Ta-
pley et al., 2004]. The mission is based on the as-
sumption that gravity variations are mostly due to
mass variations in water bodies (groundwater, soil
water, and surface water). Such data expressed in wa-
ter height can thus be used for hydrological mass bal-
ance calculations at the regional scale. After more
than one decade of gravity field processing using
spherical harmonics (SH), an alternative interpre-
tation using Mass Concentration blocks or Mas-
cons (M) of discrete spherical caps at Earth's surface
was proposed in recent years, GRACE-M [Watkins
et al., 2015]. The main disadvantage of the GRACE-
SH solution is a noisy signal which requires a fil-
tering post-processing in which an undetermined
fraction of the geophysical signal is lost [Landerer
and Swenson, 2012]. A fraction of these signal losses
can be recovered by using gain or scale factors ob-
tained, for instance, by minimizing filtering signal
degradation (undesired filter effect) using land sur-
face model (LSM) water mass outputs [Landerer and
Swenson, 2012]. Conversely, GRACE-M requires less

post-processing treatment due to lower geophysi-
cal signal losses. However, mascon gravity data in-
terpretation with a native resolution of 3° × 3° for
the Jet Propulsion Laboratory (JPL) data center so-
lution used here (see below) produces a smooth sig-
nal, which is a serious limitation for hydrological ap-
plications that require a higher resolution and thus
downscaling approaches [Wiese et al., 2016]. Two al-
ternative solutions can therefore be used: (1) GRACE-
SH, which is noisier but whose resolution (1° × 1°) is
closer to the target one (0.5° × 0.5°), and (2) GRACE-
M, which is less noisy but requires further downscal-
ing efforts. In this study, we used both solutions for
comparison. Monthly GRACE-derived Terrestrial Wa-
ter Storage (TWS) anomalies (difference between a
TWS value and its average value calculated for the
period January 2004–December 2009) were retrieved
from the JPL data center (http://grace.jpl.nasa.gov).
GRACE-SH data correspond to the RL05 solution
from the Centre for Space Research (CSR, University
of Texas, http://www.csr.utexas.edu/grace/) scaled
with the scaling factor of the NCAR-CLM4.0 model
[Landerer and Swenson, 2012]. GRACE-M data are
the RL06M solution from JPL. We retrieved the JPL
R06M filtered from the contribution of oceans [CRI
filter, Wiese et al., 2016], therefore giving the TWS
time series corresponding only to the land part of a
Mascon. The proposed scaling factors which corre-
spond in fact to downscaling factors [Scanlon et al.,
2016] were not used here since an alternative ap-
proach is proposed. These downscaling factors are
based on the mass distribution calculated by LSM
(soil–atmosphere models) while water height in the
Djeffara plain is expected to be mostly controlled by
groundwater mass variations. Therefore, we derive a
specific downscaling algorithm in this study. We used
170 monthly solutions for the period between April
2002 and December 2016 (Figure 2). This is the pe-
riod with the minimum amount of missing GRACE-
SH data and fairly close to the end of the period when
water level data are available. The trend of a TWS
time series can be analyzed and interpreted in terms
of water amount variation (e.g., in mm·yr⁻¹) over a
land surface and thus water fluxes over a domain of
interest.

As outlined above, TWS accounts inextricably for
groundwater, surface water, and soil moisture. There-
fore, and in the absence of substantial surface wa-
ter in the Djeffara plain, Soil Water Storage (SWS)
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Domain-averaged time series for the Terrestrial Water Storage (TWS) and the Soil Water Storage (SWS) using VIC and CLM models, and linear trends corresponding to TWS and SWS variations (ΔTWS and ΔSWS), (b) Geographical distribution of GRACE data and downscaling grid with in blue GRACE-SH 1° × 1° grid, in green GRACE-M 3° × 3° grid, in red the 0.5° × 0.5° grid corresponding to the CRU AAR data, in gray the 0.5° cells outside the Djeffara domain (NWSAS), and plain gray circles are the centroïds used for the downscaling in the NWSAS domain with the indication of the overall pumping rates, and (c) time series for the Terrestrial Water Storage (TWS) for the land areas of the two Mascons corresponding to the Djeffara plain with the linear trends ΔTWS.</figcaption></figure>
has to be integrated to identify Groundwater Storage (GWS). However, soil water storage variations (ΔSWS) are usually small due to almost constant and low soil moisture in arid areas, and may represent a small or even negligible component of TWS [Gonçalvès et al., 2013, Nimmo et al., 1994]. We can easily verify the relative importance of ΔSWS using GRACE-SH data whose resolution is the same as the SWS data and which correctly describes the Djeffara plain (see Figure 2b for GRACE-SH and GRACE-M grid cells location). For this purpose, soil storage data were obtained from Global Land Data Assimilation System (GLDAS), which provides soil moisture for the first meters by combining ground-based, satellite data and hydrological surface model results [Hogue et al., 2005, Rodell et al., 2004]. We retrieved monthly data from the Variable Infiltration Capacity (VIC) and the Community Land Model (CLM) versions with a 1° × 1° degree spatial resolution (http://disc.sci.gsfc.nasa.gov/). VIC and CLM were previously recognized as the best models for the Saharan region [Gonçalvès et al., 2013]. The time series for the SWS (anomalies from their average value calculated over the period 01/2004–12/2009 as TWS) and TWS averaged over the Djeffara plain are shown in Figure 2a which provides a ΔTWS value by far greater than ΔSWS. Considering GRACE-M data, Mascons #1 and #2 extend beyond the Djeffara plain domain, i.e., over the NWSAS domain (Figure 2b), making the comparison between ΔSWS and ΔTWS more difficult. The ΔTWS values for land surfaces of Mascons #1 and #2 are −6.1 ± 0.2 and −5.5 ± 0.2 mm·yr⁻¹, respectively (Figure 2c). These values are also substantially larger than the ΔSWS of less than 0.5 mm·yr⁻¹ discussed above for the Djeffara (Figure 2a) and reported for the NWSAS [Gonçalvès et al., 2013]. Therefore, for both GRACE-SH and GRACE-M data, the contribution of ΔSWS to ΔTWS is considered negligible.

## 2.3. Hydrological data

Prior to the 1940s, most of the water used was supplied by springs, notably in Tunisia. With the development of pumping wells and the continuous increase of withdrawals, the water fraction supplied by springs declined until their complete drying up [OSS, 2006]. Pumping rates are available for years 1950, 1972, 1978, and 2003 from OSS [2006], and for year 2011 (only global exploitation) from

Alfarrah and Walraevens [2018]. The most documented study from OSS [2006] presents relatively precise values for irrigation but highly approximate pumping rates for domestic water. Since most of the domestic water consumption occurs in the Tripoli area, we have considered the more accurate data by Sadeg and Ghazali [2008]. In this study, domestic water pumping (about 6% of the withdrawals) will be considered explicitly only in the Tripoli area where a salt mass balance calculation is made. For Tunisia, OSS [2006] provides overall pumping rates, i.e., domestic and irrigation uses. However, pumping in Tunisia represents only about 15% of the global withdrawals in the Djeffara plain as shown in Figure 3. This figure describes the temporal evolution of global withdrawals for the whole Djeffara showing a linear increase from almost zero in 1945 to a value of 1.0 km³·yr⁻¹ in 2011. Most of the extracted groundwater, i.e., 61%, comes from the Upper aquifer, 32% from the Middle aquifer, and the remaining fraction (7%) from the Triassic aquifer. The only area where the increase in pumping rates is not linear but seems to reach a plateau is the Tripoli region. This may be explained by a combination of salinization problems and the substitution of some part of local water demand by groundwater from aquifers in southern Libya (Nubian sandstones) conveyed by a pipeline, the Great Man-Made River (GMMR), sourcing water from outside the catchment. This is clearly true for domestic water supply until 1997 produced by local wells in the Tripoli area, and since then exclusively supplied by the GMMR at a rate of 146 × 10⁶·m³·yr⁻¹ [Alfarrah and Walraevens, 2018, Sadeg and Ghazali, 2008]. Therefore, intensive groundwater pumping for domestic water mainly concerned the Tripoli area and was interrupted in 1997 with the connection to the GMMR system.

Considering the recharging fluxes at the western and southern boundaries of the system, a constant value of 2.6 m³·s⁻¹ was estimated for the Triassic flux by means of a hydrogeological model [OSS, 2006]. Using an approach combining water stable isotopes, geostatistical and reservoir simulations, and a simple hydrological model based on piezometric levels, Gonçalvès et al. [2015, 2020] obtained a temporal variation for the recharging flux from the CI to the Djeffara aquifers showing a decline from 1.8 m³·s⁻¹ in 1950 to 0.6 m³·s⁻¹ in 2009 due to the overexploitation of the CI aquifer. In addition, global rainfall rates
<figure><img src="image_6.png" /></figure>

**Figure 3.** Temporal evolution of pumping rate. Symbols for data and solid lines are fitted regressions mostly linear except for the Tripoli area in the lower panel (exponential functions, see Section 3.3). Open symbols are used for values calculated using the total withdrawals in 2011 and pumping fractions for Tunisia and Libya. In the lower panel, domestic water is provided by local pumping wells until 1997, then wells in the Tripoli area were stopped and GMMR water was supplied; the gray line is the sum of irrigation and local domestic water consumption.

beginning in 1901 were retrieved from the CRU data center (Climatic Research Unit, http://wps-web1.ceda.ac.uk) at a resolution of 0.5° × 0.5°. These values indicate an annual domain-averaged rainfall rate of around 200 mm·yr⁻¹ for the Djeffara plain.

Piezometric information available from OSS [2006] is mostly made of cumulative drawdown Δh values in 2003, and drawdown (not absolute piezometric levels) time series but for a small fraction of piezometers. Moreover, the data is not homogeneously distributed over the domain with a higher density in Libya especially for the Upper aquifer (10 data in Tunisia against 33 in Libya). Therefore, ground-based ΔGWS estimates were focused on the Libyan Djeffara (cells 11–14, 19–21, Figure 2). Drawdown data from OSS [2006] do not fully cover the period of investigation here (2002–2016). Extending the applicability of such Δh data requires analyzing

longer time series. Piezometric levels time series until the 2010's could only be reconstructed for nine piezometers (see Figure 4a, b, and c) using punctual hydraulic head values [1972, 1978, 1993, and 2010; Gejam, 2018], and drawdown records [OSS, 2006, Saleh and Rashrash, 2015]. Among the available piezometers, well 1006 is of particular importance since it belongs to the reference piezometric network and is located at the center of the major piezometric depression affecting the Upper aquifer in Libya. This piezometer shows an almost continuous deepening of the piezometric depression. Continuous albeit moderate drawdown is also recorded in coastal piezometers (1020, 1057, and 1058; Figure 4a) while piezometer 1344 shows a stabilization which is inconsistent with the overall piezometric behavior. Piezometers in the Middle and Trias aquifers show a continuous lowering of the piezometric level.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Piezometric times series for (a) the Upper aquifer, (b) the Middle aquifer, and (c) the Trias; (d) comparison between long-term drawdown values (until 2003) from OSS [2006] and recent drawdown after 2002 in m·yr⁻¹.</figcaption></figure>

The comparison between recent drawdown (since 2002) and long-term values given by OSS [2006] shown in Figure 4d and expressed in meters per year for comparison suggest that the latter can be used for the recent period considered here (2002–2016) within an error of about 16% (regression slope uncertainty). Moreover, these drawdown values are the only independent and relevant ground-based estimates for comparison purposes with the GRACE-derived estimates. Ground-based GWS variations ΔGWS<sub>GWL</sub> can

be calculated by linear interpolation of the long-term drawdown values Δh [OSS, 2006] on a 0.1° × 0.1° grid, then an averaging on a 0.5° × 0.5° grid and by calculating ΔGWS<sub>GWL</sub> = S<sub>t</sub><sup>k</sup> × (Δh/Δt)<sup>k</sup>, with S<sub>t</sub><sup>k</sup> the stora-tivity (specific yield or storage coefficient) and Δt the time interval for the drawdown, for each aquifer k (Upper, Middle, and Trias) and by summing all these values. A specific yield at 7.5 × 10<sup>-2</sup> corresponding to well tests [Allemmoz and Olive, 1980] was used for the Upper aquifer, and values for the Middle and Trias aquifer were taken from the model calibration by OSS [2006] at 2 × 10<sup>-3</sup> and 1 × 10<sup>-3</sup>, respectively. The beginning of the drawdown was considered in 1950 for the Upper aquifer (Δt = 53 years) and 1970 (Δt = 33 years) for the Middle and Trias aquifer whose exploitation started later [Saleh and Rashrash, 2015]. The ΔGWS<sub>GWL</sub> values for the Tripoli area are given in Table 1.

### 3. Methodological approaches

#### 3.1. Downscaling GRACE-derived groundwater mass variations

Regional-scale groundwater budgets are used in this study to estimate different hydrological terms such as recharge, aquifer salinization, and SGD. The solution of such balance equations requires the groundwater mass variations ΔGWS to be estimated. GRACE data provide a valuable integrated mapping of water storage over large domains. This would be difficult to achieve using field data such as piezometric level monitoring due to discontinuous data acquisition in time and/or space. The first step is to isolate the groundwater storage anomalies (GWS) from the TWS anomalies expressed in water height to facilitate hydrological applications. In arid areas such as the Djeffara plain, surface water can generally be neglected leading to the following simple decomposition of the GRACE data:

$$TWS = GWS + SWS, \qquad (1)$$

and thus

$$\Delta TWS = \Delta GWS + \Delta SWS. \qquad (2)$$

Given the fact that ΔSWS ≪ ΔTWS (Section 2.2, Figure 2) in the case of the Djeffara plain, ΔGWS can be assimilated to ΔTWS allowing a direct interpretation of GRACE data. Due to the coarse resolution of the GRACE data, the Djeffara plain only corresponds
<table><caption>**Table 1.** Groundwater mass variations over the 2002–2016 period obtained by downscaling Grace-SH (ΔGWS<sub>SH</sub>), GRACE-M (ΔGWS<sub>M</sub>), and weighted averaged values ΔGWS<sub>A</sub> in mm·yr<sup>-1</sup></caption><thead><tr><th>Cell</th><th>ΔGWS<sub>SH</sub></th><th>ΔGWS<sub>M</sub></th><th>ΔGWS<sub>A</sub></th><th>Cell</th><th>ΔGWS<sub>SH</sub></th><th>ΔGWS<sub>M</sub></th><th>ΔGWS<sub>A</sub></th><th>ΔGWS<sub>GWL</sub></th></tr></thead><tbody><tr><td>1</td><td>-4.2 ± 0.4</td><td>-6.2 ± 0.2</td><td>-4.5 ± 0.3</td><td>11</td><td>-13.6 ± 1.3</td><td>-10.0 ± 0.3</td><td>-13.2 ± 1.2</td><td>-12.4</td></tr><tr><td>2</td><td>-2.8 ± 0.3</td><td>-8.1 ± 0.2</td><td>-3.4 ± 0.3</td><td>12</td><td>-18.1 ± 1.9</td><td>-14.8 ± 0.4</td><td>-17.8 ± 1.7</td><td>-15.3</td></tr><tr><td>3</td><td>-3.0 ± 0.3</td><td>-8.6 ± 0.27</td><td>-3.7 ± 0.3</td><td>13</td><td>-23.7 ± 2.9</td><td>-13.7 ± 0.4</td><td>-22.6 ± 2.6</td><td>-24.5</td></tr><tr><td>4</td><td>-2.4 ± 0.2</td><td>-6.8 ± 0.21</td><td>-2.9 ± 0.2</td><td>14</td><td>-15.7 ± 1.9</td><td>-8.8 ± 0.2</td><td>-14.9 ± 1.7</td><td>-14.0</td></tr><tr><td>5</td><td>-2.9 ± 0.3</td><td>-7.8 ± 0.2</td><td>-3.4 ± 0.2</td><td>15</td><td>-2.6 ± 0.2</td><td>-5.2 ± 0.2</td><td>-2.9 ± 0.2</td><td></td></tr><tr><td>6</td><td>-1.4 ± 0.1</td><td>-5.5 ± 0.2</td><td>-1.9 ± 0.1</td><td>16</td><td>-3.2 ± 0.2</td><td>-6.3 ± 0.2</td><td>-3.5 ± 0.2</td><td></td></tr><tr><td>7</td><td>-2.8 ± 0.2</td><td>-5.7 ± 0.2</td><td>-3.2 ± 0.2</td><td>17</td><td>-5.3 ± 0.5</td><td>-5.1 ± 0.1</td><td>-5.3 ± 0.5</td><td></td></tr><tr><td>8</td><td>-3.5 ± 0.3</td><td>-6.7 ± 0.2</td><td>-3.8 ± 0.2</td><td>18</td><td>-6.1 ± 0.6</td><td>-6.1 ± 0.2</td><td>-6.1 ± 0.5</td><td></td></tr><tr><td>9</td><td>-6.2 ± 0.5</td><td>-5.3 ± 0.1</td><td>-6.1 ± 0.5</td><td>19</td><td>-10.7 ± 1.0</td><td>-8.1 ± 0.2</td><td>-10.4 ± 0.9</td><td>-11.4</td></tr><tr><td>10</td><td>-7.7 ± 0.6</td><td>-6.8 ± 0.2</td><td>-7.6 ± 0.6</td><td>20</td><td>-13.2 ± 1.4</td><td>-10.5 ± 0.3</td><td>-12.9 ± 1.2</td><td>-12.2</td></tr><tr><td></td><td></td><td></td><td></td><td>21</td><td>-17.7 ± 2.1</td><td>-10.2 ± 0.3</td><td>-16.8 ± 1.9</td><td>-13.7</td></tr></tbody></table>

to seven 1° × 1° grid cells for GRACE-SH and two 3° × 3° Mascons for GRACE-M (blue and green squares, respectively in Figure 2b). Similar to what was done using the domain-averaged TWS as depicted in Figure 2, the interannual trend (slope of the linear regression on TWS time series) gives the required ΔGWS = ΔTWS values in mm·yr⁻¹ for each grid cell of GRACE-SH and GRACE-M solutions. However, to enable the integration of the CRU annual rainfall data which is at a higher resolution (0.5° × 0.5°) and to refine the water budget calculations, the ΔGWS values obtained for each of the coarse grid cells have to be downscaled. The downscaling here is based on the assumption that groundwater mass balance is mostly influenced by pumping. However, the groundwater budget within a domain described using GRACE grid cells is controlled by both internal fluxes (abstractions, recharge) and boundary net fluxes. The validity of the assumption adopted here which leads to neglect boundary fluxes is discussed in Section 5. Given that the drawdown caused by a pumping well is proportional to $Q_W \times \ln(r)$ [de Marsily, 1986] with $Q_W$ (m³·s⁻¹) the pumping rate and r the radial distance to the well, the lateral drainage (radial component) upon application of Darcy’s law may be dependent on $Q_w/r$. Since the lateral drainage and ΔGWS are related in the groundwater mass balance equation (see Section 3.2), pumping rates divided by the distances to pumping zones can be heuristically introduced as weighting factors in a downscaling approach. As illustrated in Figure 5a, for a given cell j

out of the $n_c$ 0.5° × 0.5° refined grid cells contained in a sub-domain i to be downscaled (e.g., $n_c = 4$ for a 1° × 1° grid) with a known groundwater storage variation ΔGWS<sub>i</sub>, the weight writes:

$$W_{j=1,...,n_c}^{i} = \frac{n_c \beta_{j}^{i}}{\sum_{k=1}^{n_c} \beta_{k}^{i}}, \quad (3)$$

with:

$$\beta_{k=1,...,n_c}^{i} = \sum_{l=1}^{n} \frac{Q_{W}^{l}}{r_{k}^{l}}, \quad (4)$$

where $Q_W^l$ is the total withdrawal (m³·s⁻¹) of pumping area l and $r_k^l$ is the average distance between the four corners of the refined grid cell k and the centroid of the pumping area l (Figure 5a) to avoid singularity issues ($r_k^l \to 0$). Note that the centers of the main pumping zones (centroid) were used here to avoid considering thousands of wells whose location is uncertain for a large number of them. The groundwater storage variation for a given refined grid j writes finally:

$$\Delta GWS_{j}^{i} = W_{i} \Delta GWS_{i}. \quad (5)$$

The closer the node j is to a pumping centroid and the larger the pumping rate, the greater the weight. Moreover, from (3) to (5), it can be easily verified that $\Delta GWS_i = (1/n_c) \sum_{j=1}^{n_c} \Delta GWS_j^i$. The downscaling approach proposed here is applied in two different manners depending on the GRACE solution used. For GRACE-SH, the application is straightforward using $n_c = 4$ for each 1° × 1° grid and the pumping rates and
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Methodological steps for the resolution of the linear system of groundwater mass balance equation solved for the natural recharge and the coastal fluxes. (a) Downscaling using the pumping wells centroids and overall pumping rates for GRACE-SH and GRACE-M data. (b) Discrete mass balance equations over the 0.5° × 0.5° grid with the blue arrows for the known boundary fluxes from the Trias ($F_T$) and from the CI ($F_{CI}$) and the red arrows are the unknown coastal fluxes. (c) The salt balance calculation in cell 13 corresponds to Tripoli.</figcaption></figure>

centroids within the Djeffara plain. For the downscaling of each Mascon land surface ΔGWS value, the weights are calculated for each 0.5° × 0.5° cell but using exclusively the pumping wells located in the subdomain (NWSAS or Djeffara) of the considered cell. Therefore, the main pumping areas and rates of the NWSAS within the retained Mascons had to be introduced [OSS, 2003, see Figure 2b].

For consistency reasons, we considered average pumping rates for the same period as ΔGWS values, i.e., 2002–2016. The pumping rates in 2009 considered as accurate estimates of the average values over 2002–2016 were obtained by extrapolating the reported withdrawals in 1972, 1978, and 2003 for each pumping zone (centroid). It was verified that the cumulated extrapolated pumping rates satisfy the global value in 2009 (Figure 3) within an acceptable error (here about 10%). No further refinement of the grid in the north-west of the domain (first row) was made to avoid making the linear system resolution unnecessarily more complex for an uncertain gain in precision (see discussion). The mass balance grid is therefore maintained at the resolution of the AAR data. The flux $F_{x,r}^3$ therefore represents the overall exchange with the sea in the area of Gabès-Djerba-Zarzis.

## 3.2. Mass balance approach

Groundwater storage variations ΔGWS provided by GRACE can be introduced in a regionalized mass balance equation to identify recharge and coastal exchange fluxes. For this purpose, we consider the groundwater mass balance equation for an aquifer system as:

$$ \frac{\partial}{\partial x} \left( T_x \frac{\partial h}{\partial x} \right) + \frac{\partial}{\partial y} \left( T_y \frac{\partial h}{\partial y} \right) = S_t \frac{\partial h}{\partial t} + Q_W - R_N - R_A, \quad (6) $$

where $T_x$ and $T_y$ are the transmissivities (m$^2$·s$^{-1}$) in the x and y direction, $h$ is the hydraulic head (or piezometric level in m), $S_t$ is the storativity (specific yield $S_y$ or storage coefficient S for unconfined or confined aquifers), $Q_W$ is the pumping rate in m·s$^{-1}$ (discharge in m$^3$·s$^{-1}$ divided by a cell grid surface area), and $R_N$ and $R_A$ are natural and artificial recharges, respectively (m·s$^{-1}$). The artificial recharge mostly corresponds to irrigation return flow to the aquifer which is accounted for using the relation $R_A = f_R \times Q_W$ with $f_R$ standing for the fraction of irrigation return flow (to be ascertained, see
Section 4). For the most populated region of Tripoli only, an additional term due to domestic water network leaks is considered (see Section 4). Generally, Equation (6) is solved for *h* and a single aquifer, and then the fluxes are calculated using Darcy's law. Here, identifying *h* is not the purpose since several aquifers are in fact present, and our goal is to obtain mass balance equations written in fluxes, applied to a multilayer aquifer where GRACE-derived ΔGWS values can be introduced. Upon discretization over the 0.5° × 0.5° grid considered here for the Djeffara plain (Section 3.1) using Δ*x* and Δ*y* as *x* and *y* spacing, and a time step Δ*t* of one year (6) becomes:

$$ \frac{Q_{x,l}^i - Q_{x,r}^i}{\Delta x \Delta y} + \frac{Q_{y,l}^i - Q_{y,u}^i}{\Delta x \Delta y} = S_t \frac{\Delta h_i}{\Delta t} + Q_{Wi} - R_{Ni} - R_{Ai}, \quad (7) $$

where $Q_{x,l}^i$ and $Q_{x,r}^i$ are the discharges (m³·yr⁻¹) in the *x* direction at the left and right edges of cell *i* positive if oriented in the *x* axis direction; $Q_{y,l}^i$ and $Q_{y,u}^i$ are the discharges (m³·yr⁻¹) in the *y* direction at the lower and upper edges of cell *i* positive if oriented in the *y* axis direction. Introducing the fluxes (in m·s⁻¹) corresponding to discharges (e.g., $F_{x,l}^i = Q_{x,l}^i / (\Delta x \Delta y)$) in the left-hand side of (7), and noting that $S_t \Delta h_i / \Delta t = \Delta GWS_i$, Equation (7) yields the mass balance equation applying to each of the 210.5° × 0.5° grid cells of the Djeffara which writes:

$$ F_{x,l}^i - F_{x,r}^i + F_{y,l}^i - F_{y,u}^i = \Delta GWS_i + Q_{Wi} - R_{Ni} - R_{Ai}. \quad (8) $$

In (8), $Q_{Wi}$, $R_{Ai}$, and $\Delta GWS_i$ are known at each grid cell, the left-hand side fluxes and the natural recharge are the unknowns. In the most general case, the system formed by the 21 equations (mass balance (8) for each grid of the Djeffara) is unsolvable since 76 unknown variables (fluxes and natural recharge values) are involved. Therefore, two simplifications are proposed to solve the set of equations. The first simplification is to consider a linear relation between natural recharge and AAR in arid areas [Scanlon et al., 2006], i.e., $R_{Ni} = \alpha AAR_i$ with $\alpha$ a dimensionless coefficient to identify. The time lag of several years to several decades [e.g. Scanlon et al., 2006] for the recharge to reach the water table makes the use of the annual rainfall inaccurate when the use of AAR is more adapted to account for the buffering effect of the generally thick vadose zones in arid domains. According to OSS (2006, and see Section 2.1), diffuse recharge characterized by large lag times could represent almost 75% of the overall recharge in the

Djeffara plain. The second simplification consists in considering only the main flow directions, i.e., west to east in the Tunisian Djeffara and south to north in the Libyan part (Figure 5b). This assumption relies on the general natural hydraulic gradients from reliefs to coastal zones, and on the alignment of pumping wells along the coast creating hydraulic divides also globally directed from the reliefs toward the coast. These two assumptions lead to a groundwater mass balance linear system of 21 equations and 22 unknowns (spatially variable fluxes) shown in Figure 5b which is therefore still undetermined. However, the mass balance equation system can be simplified to a stream tube form by summing groundwater balance equations of individual cells for the following blocks: 1–3; 4–6; 7–10 and 15–18; 11–19; 12–20; 13–21; 14. For all blocks (“tubes”) intermediate fluxes at meshes sides would simplify leaving only one unknown (besides $\alpha$) per equation: the coastal flux. Note that for the coastal cells, negative values of $F_{x,l}^i$ or $F_{y,u}^i$ represent seawater intrusion. Conversely, positive values correspond to groundwater discharge (SGD).

### 3.3. Salinity model for the Tripoli refined grid cell

In order to solve the indetermination problem of the groundwater mass balance equations system of Section 3.2, a simple domain-averaged salinity model for the cell corresponding to the Tripoli region (cell 13, Figure 5b and c) is considered. Given the heavy withdrawals in this area, fresh groundwater discharge is assumed negligible. We consider that the coastal flux mostly corresponds to seawater intrusion which occurs mainly in the Upper aquifer, therefore neglecting the seawater fluxes for deeper aquifers. The time evolution of the salinity in the Tripoli area is sufficiently well documented which enables the calculation of average salinity values for the corresponding 0.5° × 0.5° grid cell using the salinity maps available for 1950, 1960, and 2000 from OSS [2006] for the Upper aquifer. The simple salinity mass balance equation in the Upper aquifer for cell 13 writes:

$$ \frac{\partial(VC)}{\partial t} = -Q_W^T C + f Q_W^T C_1 + Q_M C_M + Q_S C_0, \quad (9) $$

where $V$ (m³) is the volume of the Upper aquifer in the Tripoli grid cell, $f$ is the fraction of the pumping in the Tripoli area $Q_W^T$ (m³·s⁻¹) which returns to the aquifer (irrigation return flow and/or domestic water
network leaks) with a salt concentration $C_1$ (kg·m⁻³), $Q_M$ and $Q_S$ (m³·s⁻¹) are the seawater entering the aquifer and the flux coming from the south, and $C_M$ and $C_0$ (kg·m⁻³) are the seawater salinity and the natural salinity (water rock interactions equilibrium concentration). Note that in (9) the natural recharge does not appear because it involves the concentration of rainwater (or evaporated rainwater) and represents a negligible term compared to the others. We assume that the time variation of the total pumping rate in the Tripoli area (Figure 3) can be described by:

$$Q_W^T = A(1 - \exp(-\beta t)), \quad (10)$$

with the fitted values $A = 13.3$ m³·s⁻¹, $\beta = 0.05$ yr⁻¹ and a time origin in 1948 when $Q_W^T = 0$. Assuming that seawater intrusion and southern fluxes closely follow the global pumping rate time evolution leads to consider similar expressions for $Q_M$ and $Q_S$ writing:

$$Q_{M,S} = A_{M,S}(1 - \exp(-\beta t)). \quad (11)$$

A simplified resolution of the salt balance equation (9) can be achieved by considering that the return flow is mainly related to irrigation (85% of the water pumped before 2000, i.e., the calibration period). The salt is extracted from aquifers at a rate $Q_W^T \times C$ but the fraction $(1-f)$ of the groundwater abstraction $Q_W^T$ is lost by evaporation leaving the salt in the top soil. Assuming a perfect salt leaching by the excess irrigation water, the concentration in the water returning to the aquifers is $C_1 = [(1-f)C+fC]/f = C/f$, which exactly corresponds to an evaporation effect. Introducing this expression for $C_1$ in the mass balance equation (9), the two first terms at the right-hand side vanish (total salt return). This simplification is therefore equivalent to considering that the salt balance is controlled by lateral fluxes ($Q_M$, $Q_S$) depending on pumping but is not directly affected by withdrawals due to irrigation practices. In addition, assuming small temporal variations for $V$, and introducing (11) into (9) yields, upon integration:

$$C = \left( \frac{A_M C_M}{V} + \frac{A_S C_0}{V} \right) \left[ t + \frac{1}{\beta} (\exp(-\beta t) - 1) \right] + C_0. \quad (12)$$

Equation (12) is used to reproduce the cell-averaged salinity values available at three measurement times by calibrating the initial (natural) salinity of the Upper aquifer in 1948 $C_0$, and $A_M$ and $A_S$. Once calibrated, Equation (11) can be used to calculate the seawater flux.

## 4. Results

GRACE-derived ΔGWS for the period 2002–2016 were downscaled based on the distance to the centroids of the pumping areas and on the pumping rates as reconstructed in 2009 (see Section 3.1, Figures 2b and 6a for locations and rates) to obtain the same resolution as the AAR data from the CRU (0.5° × 0.5°, Figure 6c). Downscaling produces realistic groundwater depletion especially in the region of Tripoli concentrating about one half of the Libyan withdrawals as shown in Table 1 and Figure 6b using GRACE-SH and GRACE-M data. A better agreement between downscaled GRACE-derived ΔGWS values and ground-based estimates ΔGWS<sub>GWL</sub> is obtained using GRACE-SH than GRACE-M as shown by the Root Mean Square error (RMSE) values in Figure 7. The discrepancy is mostly related to cell 13 (Tripoli) with a large underestimation. For groundwater budget calculations proposed hereafter, the average values ΔGWS<sub>A</sub> calculated as weighted mean values of ΔGWS<sub>SH</sub> and ΔGWS<sub>M</sub> using the MSE as weighting factors were adopted. The calculated domain-averaged value (using ΔGWS<sub>A</sub> in Table 1) for ΔGWS over the 0.5° square grid is about −8.0 ± 0.2 mm·yr⁻¹, which confirms that the depletion shown by the piezometric levels decline. GRACE data points to an overall aquifer depletion with ΔGWS values varying in space between −1.9 and −22.6 mm·yr⁻¹. The highest depletion rate is observed in the Libyan part (Tripoli area).

The transient salt balance for the 0.5° × 0.5° cell corresponding to Tripoli (cell 13 in Figure 5b, see Section 3.3) was calculated to obtain an independent estimate of the seawater flux at the shoreline. This flux is required to enable solving the linear system (8) of groundwater mass balance equations accounting for the downscaled ΔGWS values. Combining the mass balance equations for cells 13 and 21 (Figure 5b) yields:

$$F_T - F_{y,u}^{13} = \Delta GWS_{13} + \Delta GWS_{21} - \alpha(AAR_{13} + AAR_{21}) + (1 - f_R)Q_W^{13} - R_A^D, \quad (13)$$

where $R_A^D$ (mm·yr⁻¹) is the leaks from the domestic water system that recharges the Upper aquifer, which is considered only for cell 13 corresponding to Tripoli (near 1.5 × 10⁶ inhabitants). According to CEDARE [2014], these leaks represent 43 ± 4% of the 146 × 10⁶ m³·yr⁻¹ of domestic water supplied to the Tripoli area by the GMMR. This corresponds to an
<figure><img src="image_9.png" /><figcaption>**Figure 6.** (a) GRACE-derived ΔGWS at 1° × 1° resolution using GRACE-SH data over the Djeffara Plain, black diamonds indicate the location of the pumping areas centroids (number and rate in 2009) used for the downscaling within the Djeffara plain domain, (b) downscaled GRACE-derived ΔGWS using GRACE-SH data, and (c) AAR data from the CRU.</figcaption></figure>

artificial recharge of 24.3 ± 2.3 mm·yr⁻¹ considering the cell surface area of about 2580 km². Note that in (13) applied to the period 2002–2016, Q<sub>W</sub><sup>13</sup> only corresponds to irrigation pumping. Solving (13) requires identifying the fraction of irrigation return flow f<sub>R</sub>. For this purpose, we consider the available estimates

for the irrigation efficiency:

$$E_{\text{irr}} = \frac{Q_D}{Q_W}, \qquad (14)$$

for the Libyan Djeffara, where Q<sub>D</sub> is the plant water demand while Q<sub>W</sub> is the water effectively supplied
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Comparison between GRACE-derived (using GRACE-SH ΔGWS<sub>SH</sub> or GRACE-M ΔGWS<sub>M</sub>) and Ground-derived (ΔGWS<sub>GWL</sub>) groundwater storage variations in the Tripoli area with the Root Mean Square Error (RMSE) associated to each alternative GRACE data used.</figcaption></figure>

(pumping). According to Abdudayem and Scott [2014], $E_{\text{irr}}$ is about 0.64 which is in agreement with other estimates observed in southern parts of Libya [Shaki and Adeloye, 2006]. Given the fact that water supplied satisfies the evapotranspiration demand of cultivated plants, the water provided in excess infiltrates except a certain fraction which is evaporated from the soil surface. Quantifying this evaporation from the soil surface is not straightforward despite its pivotal importance for water balance calculations [Dewandel et al., 2008]. Here, we consider a simple calculation derived from the method proposed by Allen et al. [1998] in an irrigation context. Irrigation water saturates a soil surface layer whose thickness is the excess irrigation (in m) divided by the porosity. A fraction of this saturated layer corresponding to the drainage porosity $S_y$ is immediately removed by vertical drainage due to gravity, the remaining part evaporates. Assuming that evaporation can bring the soil moisture of top soils to almost zero (residual porosity is evaporated), our simple model writes:

$$f_R = (1 - E_{\text{irr}}) \frac{S_y}{\omega}, \quad (15)$$

where $\omega$ is the total porosity. More elaborated models are available, but require additional parameterization such as unsaturated flow coefficients [Dewandel et al., 2008, Or and Lehmann, 2019]. In the Tripoli area, Allemmoz and Olive [1980] reported a value for $S_y$ of between 5 and 10%. Considering this range as the 95% confidence interval leads to $S_y = 7.50 \pm 1.25\%$. Together with a value of 10% for $\omega$ [Sadegh and Karahanoglu, 2001], this $S_y$ estimate leads to $f_R = 27 \pm 4\%$ consistent with the value of 30% considered by Sadegh and Karahanoglu [2001] in their modeling study. The volume $V$ of groundwater in cell 13 for the Upper aquifer was calculated at $3.87 \times 10^{10} \text{ m}^3$ using $\omega = 10\%$, a mean thickness of 150 m [OSS, 2006], and a cell surface area of 2580 km$^2$. We considered an uncertainty of 20% for $V$ in particular to account for the error associated with the assumption of a time constant volume (see discussion). This value for $V$ was introduced in (12) to reproduce the average salinity evolution from 1948 to 2000 as shown in Figure 8. Then, upon calibration of $C_0$ (0.63 g·L$^{-1}$), $A_S$ (2 m$^3$·s$^{-1}$) and $A_M$ (1.1 m$^3$·s$^{-1}$), an average value of $1.0 \pm 0.2 \text{ m}^3$·s$^{-1}$ for the seawater intrusion in mesh 13 for the period 2002–2016 was calculated using (11). This value corresponds to $F_{y,u}^{13} = -13.2 \pm 2.8 \text{ mm} \cdot \text{yr}^{-1}$ (discharge divided by surface area of cell 13 and negative sign due to incoming flux orientation). Introducing the value for $F_{y,u}^{13}$ into (13) yields $\alpha = 2.2 \pm 1.4\%$ using the average value for $Q_W^{13}$ (irrigation only) of $332 \times 10^6 \text{ m}^3$·yr$^{-1}$ for the period 2002–2016 obtained using the exponential fit of Figure 3 (Equation (10)) with $A = 13.3 \text{ m}^3$·s$^{-1}$ and $\beta = 0.05 \text{ yr}^{-1}$. Now that $\alpha$ is known, the full linear system of groundwater mass balance equations can be solved for the coastal fluxes (red arrows in Figure 5b). The results are given in Table 2.

## 5. Discussion

Several simplifying hypotheses supporting our approach must first be discussed. The downscaling approach as developed here requires boundary fluxes for a given domain described by GRACE cells to be low in comparison with the internal withdrawals. This requirement is satisfied using GRACE-SH data, the resolution of which allows obtaining a domain close to the real geographical extent of the Djeffara characterized by low boundary fluxes. Using coarser GRACE-M grids which extend beyond the Djeffara
<table><caption>**Table 2.** Terms of the groundwater budget for the 21 cells in Million m³·yr⁻¹</caption><thead><tr><th>Cell i</th><th>F<sup>i</sup>(×10<sup>6</sup> m<sup>3</sup>·yr<sup>-1</sup>)</th><th>Q<sup>i</sup><sub>W</sub>(×10<sup>6</sup> m<sup>3</sup>·yr<sup>-1</sup>)</th><th>R<sup>i</sup><sub>N</sub>(×10<sup>6</sup> m<sup>3</sup>·yr<sup>-1</sup>)</th><th>R<sup>i</sup><sub>A</sub>(×10<sup>6</sup> m<sup>3</sup>·yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1</td><td>41.0</td><td>—</td><td>10.5</td><td>—</td></tr><tr><td>2</td><td>-3.7</td><td>87.0</td><td>10.0</td><td>23.5</td></tr><tr><td><b>3</b></td><td><b>-46.9</b></td><td><b>87.0</b></td><td><b>10.8</b></td><td><b>23.5</b></td></tr><tr><td>4</td><td>17.6</td><td>—</td><td>10.0</td><td>—</td></tr><tr><td>5</td><td>37.2</td><td>—</td><td>10.8</td><td>—</td></tr><tr><td><b>6</b></td><td><b>54.3</b></td><td><b>—</b></td><td><b>12.3</b></td><td><b>—</b></td></tr><tr><td>7</td><td>17.5</td><td>—</td><td>9.3</td><td>—</td></tr><tr><td>8</td><td>36.9</td><td>—</td><td>9.5</td><td>—</td></tr><tr><td>9</td><td>63.0</td><td>—</td><td>10.5</td><td>—</td></tr><tr><td><b>10</b></td><td><b>170.9</b></td><td><b>—</b></td><td><b>11.8</b></td><td><b>—</b></td></tr><tr><td><b>11</b></td><td><b>79.8</b></td><td><b>30.0</b></td><td><b>12.8</b></td><td><b>8.1</b></td></tr><tr><td><b>12</b></td><td><b>-113.3</b></td><td><b>324</b></td><td><b>14.1</b></td><td><b>87.5</b></td></tr><tr><td><b>13</b></td><td><b>-34.0</b></td><td><b>331.0</b></td><td><b>15.0</b></td><td><b>89.4(62.8)*</b></td></tr><tr><td><b>14</b></td><td><b>34.8</b></td><td><b>47.0</b></td><td><b>14.4</b></td><td><b>12.7</b></td></tr><tr><td>15</td><td>15.1</td><td>—</td><td>7.7</td><td>—</td></tr><tr><td>16</td><td>32.4</td><td>—</td><td>8.3</td><td>—</td></tr><tr><td>17</td><td>33.9</td><td>29.0</td><td>9.1</td><td>7.8</td></tr><tr><td>18</td><td>76.4</td><td>—</td><td>10.3</td><td>—</td></tr><tr><td>19</td><td>54.9</td><td>—</td><td>11.8</td><td>—</td></tr><tr><td>20</td><td>63.3</td><td>—</td><td>13.5</td><td>—</td></tr><tr><td>21</td><td>71.7</td><td>—</td><td>12.7</td><td>—</td></tr></tbody></table>

Lateral fluxes $F^i$ (8) are positive when oriented in the same direction than x or y axis (see Figure 5). * Value between parentheses is the domestic water leaks considered only for the Tripoli area. Coastal cells are indicated in bold character. The boundary fluxes $F_{CI}$ (cell 1) and $F_T$ (cells 18–21, 14; see Figure 5) are 18.9 × 10⁶ and 16.4 × 10⁶ m³·yr⁻¹, respectively.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Calculated temporal evolution of the salinity (blue plain line) for cell 13 corresponding to Tripoli and average salinity values (red diamonds) derived from salinity maps in 1950, 1960, and 2000 [OSS, 2006].</figcaption></figure>

limits and more precisely in the NWSAS domain raises the question of the boundary fluxes. According to the piezometric maps for the NWSAS [Baba-Sy, 2005, Besbes and Horriche, 2008], the net fluxes at the southern and western boundaries of Mascon#1 (Figure 2b) are almost zero for the CT (orthogonal hydraulic head contours). Regarding the CI, the piezometric map points to an entering flux at the western boundary of Mascon#1 corresponding to the outflow of the CI in the Djeffara plain estimated at 0.6 m³·s⁻¹ (Section 2.3). Therefore, the neglected net boundary flow for Mascon#1 represents no more than 0.2 mm·yr⁻¹ considering its surface area. For Mascon#2, the Nafusa represents a recharge area for both the CI and CT with a cumulated value at 1.7 m³·s⁻¹ calculated by Baba-Sy [2005]. This would
represent a neglected net boundary outflow of less
than −0.5 mm·yr⁻¹ for Mascon #2. For both GRACE
version (M and SH), considering small boundary
fluxes is thus a reasonable assumption. In case of
non-negligible boundary fluxes, the application of
the proposed approach is more complex and requires
additional steps: (i) identification of the net bound-
ary flux, (ii) subtraction of this flux from the GRACE-
derived ΔGWS signal to isolate the purely internal
ΔGWS, (iii) application of the proposed downscaling
approach and, (iv) addition of the net boundary flux
to downscaled ΔGWS values.

A salt balance provides coastal flux estimates, as-
suming a constant volume of groundwater in the
phreatic aquifer for cell 13 (Tripoli) along the sim-
ulation extending over the period 1948–2000. As ac-
tive recharge and withdrawals are at work, this may
represent a rough approximation. For the calcula-
tion period (1948–2000), the average pumping rate is
267 × 10⁶ m³·yr⁻¹, of which about 35% (average value
between domestic leaks and irrigation return flow, *f*)
returns to the aquifer (93 × 10⁶ m³·yr⁻¹, about 80%
comes from irrigation and 20% from domestic water
leaks). Recharge is approximately 15 × 10⁶ m³·yr⁻¹,
average seawater intrusion 22 × 10⁶ m³·yr⁻¹, and the
flux out of cell 21 is about 40 × 10⁶ m³·yr⁻¹ (temporal
average values of Q<sub>M</sub> and Q<sub>S</sub> using (11)). The net an-
nual average water withdrawal from cell 13 is there-
fore 97 × 10⁶ m³·yr⁻¹ which leads to a net volume ex-
tracted from the phreatic aquifer over the 52 years
of simulation of 5.0 km³ representing a groundwa-
ter storage variation of about 13%. This is consistent
with the 20% uncertainty considered for *V* in the cal-
culation of seawater intrusion flux. In this study, we
considered only principal flow directions, i.e., east to
west in Tunisia and south to north in Libya, which
makes it possible to solve the discrete linear sys-
tem of mass balance equations. This clear simplifi-
cation is however relatively well supported by piezo-
metric data [e.g. OSS, 2006, Sadeg and Karahanoglu,
2001, Trabelsi et al., 2012]. Although adapted to the
overall objective of fluxes assessment using parsimo-
nious modeling approaches and to the resolution of
the rainfall data, the mesh may appear coarse. How-
ever, noting that the nine pumping zones are mostly
aligned along the 450 km of the shoreline, the average
spacing of hydraulic divides associated with pump-
ing is about 50 km which corresponds to the grid res-
olution. This resolution enables to compensate for

the fact that the pumping wells are actually located
around the centroids which were conveniently used
here for downscaling and mass balance calculations.
Even though small portions of the Djeffara plain in
the south of the domain covering no more than half
a grid cell were excluded, we believe that the spa-
tial resolution does not significantly affect the results.
As mentioned in Section 3.1, additional downscal-
ing and thus grid refinement would have been es-
pecially beneficial for the area of Gabès (Figure 1).
However, Chulli et al. [2007] estimated the SGD in the
gulf of Gabès at 0.2 m³·s⁻¹ (0.1 m³·s⁻¹ for half of the
coast length) using Darcy's calculations for the Upper
aquifer. This relatively low value also obtained in a re-
cent modeling study by Vernoux et al. [2020] justifies
maintaining the grid resolution as adopted here.

Regionalized groundwater mass balance pro-
posed here and based on GRACE and ground data
provides valuable estimates for recharge and coastal
fluxes. Our approach leads to a natural recharge
representing 2.2 ± 1.4% of the rainfall which corre-
sponds to a mean value of 4.3 ± 2.7 mm·yr⁻¹ con-
sidering the average rainfall in the Djeffara plain
at 196 mm·yr⁻¹. This value is relatively consistent
with the estimates from hydrogeological model-
ing at 6.8 mm·yr⁻¹ [OSS, 2006]. In the context of
circum-Saharan arid domains, this rate of 2.2 ± 1.4%
is also in good agreement with some estimates for
the NWSAS [2.2 ± 1.4%, Gonçalvès et al., 2013], and
the Nubian Sandstones Aquifers System [1.4 ± 1.7%,
Mohamed et al., 2017] using global groundwater
mass balance based on GRACE data. The fluxes ob-
tained here enable to propose a global groundwater
budget for the Djeffara plain as shown in Table 3.
This groundwater budget is characterized by an out-
flow exceeding the inflow in accordance with the
overall groundwater storage depletion shown by
GRACE data. In addition to the Tripoli area, seawa-
ter intrusion also occurs in the Zarzis-Djerba area
(see Figure 1) which is consistent with geochemical
evidences from Trabelsi et al. [2012]. In Libyan Dj-
effara, a lower seawater intrusion is calculated for
cell 13 in comparison to cell 12 due to the domestic
water leaks considered only for the most populated
area (cell 13 Tripoli; see (13) and Table 2) while the
remaining terms of the groundwater budget (includ-
ing pumping) are rather similar. The overall seawa-
ter intrusion amounts to 193 × 10⁶ m³·yr⁻¹ which
is larger than the 34.7 × 10⁶ m³·yr⁻¹ estimated from
<table><caption>**Table 3.** Global groundwater budget for the Djeffara plain (2002–2016) obtained using GRACE-derived ΔGWS downscaled over a 0.5° × 0.5° grid</caption><thead><tr><th>Input (×10<sup>6</sup> m<sup>3</sup>·yr<sup>-1</sup>)</th><th colspan="2">Output (×10<sup>6</sup> m<sup>3</sup>·yr<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Natural Recharge</td><td><strong>236.3</strong> ± 154.7</td><td>Pumping</td><td>933.5</td></tr><tr><td>Irrigation return flow</td><td>252.5 ± 38</td><td>Coastal Groundwater discharge</td><td><strong>339.6</strong> ± 56.7</td></tr><tr><td>Domestic water leaks</td><td>62.8 ± 5.8</td><td></td><td></td></tr><tr><td>Seawater intrusion</td><td><strong>193.0</strong> ± 32.3</td><td></td><td></td></tr><tr><td>Trias discharge</td><td>82.0</td><td></td><td></td></tr><tr><td>CI discharge</td><td>18.9</td><td></td><td></td></tr><tr><td>Total</td><td><strong>845.5</strong> ± 162.6</td><td>Total</td><td><strong>1273.1</strong> ± 56.7</td></tr></tbody></table>

Uncertainty propagation using the error on ΔGWS<sub>A</sub>, α and f<sub>R</sub>.

the OSS model (2006). Total submarine groundwater discharge (SGD) is 339.6 × 10⁶ m³·yr⁻¹ which is also higher than the value at 186 × 10⁶ m³·yr⁻¹ found with the OSS model (2006). These discrepancies can be explained by a more important recharge deduced here (10.2 mm·yr⁻¹ including an artificial recharge) compared to the model by OSS (2006; 6.8 mm·yr⁻¹ only comprising the natural contribution). This higher recharge, if taken into account in the hydrogeological model, would lead to larger lateral drainage (SGD) and impact the model calibration probably affecting the calculated seawater intrusion fluxes.

Distributed along the 450 km long coastline, this corresponds to a shore-length-normalized SGD flux 7.5 × 10⁵ m³·km⁻¹·yr⁻¹. This is one of the first semi-experimentally derived SGD estimates for the north African coastline, a region that presents a notorious data gap in coastal hydrological data [e.g., El-Gamal et al., 2012; see Figure 1 in Rodellas et al., 2015]. This value compares well with model-derived estimates for freshwater SGD on the scale of the entire Mediterranean basin of 8.1 × 10⁵ m³·km⁻¹·yr⁻¹ [Zektser and Everett, 2004] and 1.1 × 10⁶ m³·km⁻¹·yr⁻¹ [Zektser et al., 2007], perhaps indicating that the aquifer system under investigation may be representative for the semi-arid Mediterranean basin. Indeed, the few studies carried out in the southern/eastern Mediterranean regions with similar semi-arid conditions also yield comparable results, e.g., 1.8 × 10⁶ m³·km⁻¹·yr⁻¹ [Weinstein et al., 2007]. This estimate obtained for Israel's Carmel Coast using a combination of radon and salinity as tracers is a factor 2 of the flux estimates obtained here for the Djeffara. Larger scale estimates obtained from radium isotope mass balances such as

those reported by Rodellas et al. [2015] ranging from 4.0 × 10⁵ m³·km⁻¹·yr⁻¹ to 3.7 × 10⁸ m³·km⁻¹·yr⁻¹ (average 3.7 × 10⁸ m³·km⁻¹·yr⁻¹) cannot be compared in a meaningful way to the above freshwater SGD fluxes because radium-derived estimates include both freshwater SGD and (saline) seawater recirculation flux, and are correspondingly significantly larger. Although a direct comparison between our results and the radium-derived estimates is difficult, several objective reasons can however support the hypothesis of lower SGD in the southern Mediterranean domain: more extensive groundwater exploitation, lower recharge rate (recharge to AAR ratio), and lower annual precipitation. All these factors contribute to a relatively modest lateral drainage in the aquifers and therefore to lower SGD values. The groundwater exploitation results in water being removed from the lateral natural groundwater flow to the coast. For the Djeffara plain, the net water withdrawn corresponds to the pumping minus the irrigation return flow, i.e., about 70% of the exploited volumes or about 0.7 km³·yr⁻¹. This represents 1.56 × 10⁶ m³·yr⁻¹·km⁻¹ suggesting a still low “natural” SGD value. A recharge ratio of between 23 and 52% (average 35%) for karstified carbonates [Allocca et al., 2014, Arfib and Charlier, 2016, Baudement et al., 2017, Martos-Rosillo et al., 2015, Polemio, 2016, Zagana et al., 2007], and between 10 and 21% (average 16%) for detritic aquifer formations [Séraphin et al., 2016, Yagbasan, 2016, Zagana et al., 2007] were reported in the north Mediterranean area (Spain to Turkey). Given that karstic aquifers represent almost 20% of the coastal areas in Mediterranean Europe [Goldscheider et al., 2020], an average value of about 20% for the recharge ratio, therefore tenfold higher than in the southern
part of the Mediterranean domain (2 to 3%), can be considered. In addition, the average annual rainfall is generally two to three times higher on average in the coastal areas of north Mediterranean Sea. Together with the recharge ratio, this could readily explain lateral groundwater drainage toward the coasts and thus SGD values one order of magnitude higher in the north Mediterranean domain in comparison to southern parts. To date, SGD studies are primary carried out on a local scale (beach/embayment). There is an increasing interest in upscaling these coastal fluxes to, e.g., basin scale [e.g. Rodellas et al., 2015]. These large-scale estimates have their own limitations, e.g., the choice of approach which in the case of Rodellas et al. [2015] does not allow a determination of the freshwater fraction of SGD which is of interest to hydrological studies. Based on the results of this study, we propose that the correct unit to report SGD fluxes may be the catchment/aquifer scale as investigated here. It is important to note that in our GRACE-derived estimate, no SGD flux occurs from 3 out of 7 coastal cells of the Djeffara aquifer investigated in this study (Table 2). This spatial variability has to be considered when upscaling, and we suggest that GRACE data could be used in combination with local “point measurements” with other approaches to provide more accurate upscaling estimates in future studies elsewhere.

## 6. Conclusions

In this study, we proposed a downscaling approach for the GRACE-derived groundwater mass storage variations based on the pumping wells spatial distribution and rates of exploitation. In arid zones, where groundwater is the main or the only water resource, pumping is often the dominant term of the hydrological budget. This approach leads to a realistic spatial distribution of the groundwater storage variations averaged over the period 2002–2016 with lower values in the areas where large drawdown of the piezometric levels are reported, such as the Tripoli area. Downscaled groundwater storage variations were then introduced in a set of 21 discrete groundwater balance equations over the finer 0.5° × 0.5° (2580 km²) grid to obtain new estimates of the recharge and coastal fluxes (outflowing or incoming) independently from any hydrogeological model. Due to indetermination issues, the saline intrusion rate in

the Tripoli region was identified separately using a local salt mass balance calculation. According to our calculations, the Tripoli and Djerba-Zarzis areas are subjected to a strong seawater intrusion. The remaining coastal zones produce a cumulated groundwater outflow (SGD) of $339.6 \times 10^6 \text{ m}^3\cdot\text{yr}^{-1}$ for the 2002–2016 period. This geophysically based value fills an important gap in SGD estimations for the southern Mediterranean coastal area which remains unquantified to date, including by radiochemical tracers widely used elsewhere. The average recharge over the Djeffara plain was estimated at $4.3 \pm 2.7 \text{ mm}\cdot\text{yr}^{-1}$ or $2.2 \pm 1.4\%$ which is consistent with the recharge to average annual rainfall ratio found for other Saharan aquifers. The GRACE-derived domain-averaged water budget which is obtained here, i.e., overall groundwater depletion at $-8.0 \pm 0.2 \text{ mm}\cdot\text{yr}^{-1}$, not surprisingly confirms the unsustainable mining of groundwater in the Djeffara plain. Available options to stop this unsustainable mining include pumping limitations to prevent additional salinization effects, use of underexploited resources such as the Jurassic aquifer at the border of the Djeffara, exploitation of remote but unsustainable resources (GMMR), irrigation efficiency increase, and seawater desalinization. These are all ongoing (or future) options under consideration to improve water uses and access to water. Here, we have taken advantage of the alignment of the main groundwater flow directions (toward the coast) with the GRACE grids. This favorable situation can be found for other Mediterranean coastal aquifers enabling a similar approach of distributed mass balance using satellite gravity.

## Acknowledgments

GRACE land data are available at http://grace.jpl.nasa.gov, supported by the NASA MEaSURES Program. This study has been supported by the Labex OT-Med project of the French program "Investissement d'Avenir." The two anonymous reviewers and the Editors are acknowledged for their relevant comments which largely improved the manuscript.

## References

Abdudayem, A. and Scott, A. H. S. (2014). Water infrastructure in Libya and the water situation in
agriculture in the Jefara region of Libya. *Afr. J. Econ. Sustain. Dev.*, 3, 33–64.

Ahmed, M., Sultan, M., Wahr, J., and Yan, E. (2014). The use of GRACE data to monitor natural and anthropogenic induced variations in water availability across Africa. *Earth Sci. Rev.*, 136, 289–300.

Alfarrah, N. and Walraevens, K. (2018). Groundwater overexploitation and seawater intrusion in coastal areas of arid and semi-arid regions. *Water*, 10, article no. 143.

Allemmoz, M. and Olive, P. (1980). Recharge of groundwaters in arid areas: case of the Djeffara Plain in Tripolitania, Libyan Arab Jamahiriya. In *Arid Zone Hydrology Investigations with Isotopic Techniques*, pages 181–191. IAEA, Vienna.

Allen, R. G., Pereira, L. S., Raes, D., and Smith, M. (1998). Crop evapotranspiration: Guidelines for computing crop requirements. Irrigation and Drainage, Paper No. 56, FAO, Rome, Italy.

Allocca, V., Manna, M., and De Vita, P. (2014). Estimating annual groundwater recharge coefficient for karst aquifers of the southern Apennines (Italy). *Hydrol. Earth Syst. Sci.*, 18, 803–817.

Arfib, B. and Charlier, J. B. (2016). Insights into saline intrusion and freshwater resources in coastal karstic aquifer using a lumped Rainfall-Discharge-Salinity model (the Port-Miou brackish spring, SE France). *J. Hydrol.*, 540, 148–161.

Baba-Sy, O. (2005). *Recharge et paléorecharge du système aquifère du Sahara septentrional*. PhD thesis, Univ. Tunis, Tunis, Tunisia.

Baudement, C., Arfib, B., Mazzilli, N., Jouves, J., Lamarque, T., and Guglielmi, Y. (2017). Groundwater management of a highly dynamic karst by assessing baseflow and quickflow with a rainfall-discharge model (Dardennes springs, SE France). *BSGF Earth Sci. Bull.*, 188, article no. 40.

Ben Dhia, H. (1987). Geothermal energy in Tunisia: Potential of the southern province. *Geothermics*, 16, 299–318.

Ben Lasmar, R., Guellala, R., Garrach, M., Mahroug, A., Naouali, B. S., and Inoubli, M. H. (2017). The Tunisian Jurassic aquifer in the North African Sahara aquifer system: information derived from two-dimensional seismic reflection and well logs. *Hydrogeol. J.*, 25, 2281–2301.

Besbes, M. and Horriche, F. J. (2008). Design of a piezometric monitoring network for the North

Western Sahara Aquifer System. *Secheresse*, 18, 13–22.

Bonsor, H. C., Shamsudduha, M., Marchant, B. P., MacDonald, A. M., and Taylor, R. G. (2018). Seasonal and decadal groundwater changes in African sedimentary aquifers estimated using GRACE products and LSMs. *Remote Sens.*, 10, article no. 904.

Carvalho Resende, T., Longuevergne, L., Gurdak, J. J., Leblanc, M., Favreau, G., Ansems, N., Van der Gun, J., Gaye, C. B., and Aureli, A. (2019). Assessment of the impacts of climate variability on total water storage across Africa: implications for groundwater resources management. *Hydrogeol. J.*, 27, 493–512.

CEDARE (2014). Libya water sector M&E rapid assessment report. Monitoring & Evaluation for Water In North Africa (MEWINA) Project, Water Resources Management Program, CEDARE, 83 p.

Chulli, A., Tizro, A. T., and Jabnoun, N. (2007). Modélisation de l'intrusion marine dans l'aquifère côtier de Gabès (sud tunisien). In *Proceedings of Symposium "A New Focus on Groundwater-Seawater Interactions"*, Perugia, Italy, volume 312. IAHS Publ., https://iahs.info/uploads/dms/13946. 32-219-226-312_07_Chulli.pdf. Assessed 27 May 2020.

de Marsily, G. (1986). *Quantitative Hydrogeology: Groundwater Hydrology for Engineers*. Academic Press, New-York, USA.

Dewandel, B., Gandolfi, J.-M., de Condappa, D., and Ahmed, M.-S. (2008). An efficient methodology for estimating irrigation return flow coefficients of irrigated crops at watershed and seasonal scale. *Hydrol. Process.*, 22, 1700–1712.

El-Gamal, A. A., Peterson, R. N., and Burnett, W. C. (2012). Detecting freshwater inputs via groundwater discharge to Marina Lagoon, Mediterranean Coast, Egypt. *Estuar. Coast.*, 35, 1486–1499.

Gejam, A. M. S. (2018). *Groundwater Potentiality in Jefara Plain*. PhD thesis, Tripoli University, Tripoli, Libya.

Goldscheider, N., Chen, Z., Auler, A. S., Bakalowicz, M., Broda, S., Drew, D., Hartmann, J., Jiang, G., Moosdorf, N., Stevanovic, Z., and Veni, G. (2020). Global distribution of carbonate rocks and karst water resources. *Hydrogeol. J.*, 28, 1661–1677.

Gonçalvès, J., Deschamps, P., Hamelin, B., Vallet-Coulomb, C., Petersen, J., and Chekireb, A. (2020). Revisiting recharge and sustainability of the North-
Western Sahara aquifers. *Reg. Environ. Change, 20*, article no. 47.

Gonçalvès, J., Petersen, J., Deschamps, P., Hamelin, B., and Baba-Sy, O. (2013). Quantifying the modern recharge of the “fossil” Sahara aquifers. *Geophys. Res. Lett.*, 40, 2673–2678.

Gonçalvès, J., Vallet-Coulomb, C., Petersen, J., Hamelin, B., and Deschamps, P. (2015). Declining water budget in a deep regional aquifer assessed by geostatistical simulations of stable isotopes: case study of the Saharan Continental Intercalaire. *Hydrol. J.*, 531, 821–829.

Hogue, T. S., Bastidas, L., Gupta, H., Sorooshian, S., Mitchell, K., and Emmerich, W. (2005). Evaluation and transferability of the Noah land surface model in semiarid environments. *J. Hydrometeorol.*, 6, 68–84.

Landerer, F. W. and Swenson, S. C. (2012). Accuracy of scaled GRACE terrestrial water storage estimates. *Water Resour. Res.*, 48, article no. W04531.

Leblanc, M. J., Tregoning, P., Ramillien, G., Tweed, S. O., and Fakes, A. (2009). Basin-scale, integrated observations of the early 21st century multiyear drought in southeast Australia. *Water Resour. Res.*, 45, article no. W04408.

Luijendijk, E., Gleeson, T., and Moosdorf, N. (2020). Fresh groundwater discharge insignificant for the world's oceans but important for coastal ecosystems. *Nat. Commun.*, 11, article no. 1260.

Martos-Rosillo, S., Gonzalez-Ramon, A., Jimenez-Gavilan, P., Andreo, B., Duran, J. J., and Mancera, E. (2015). Review on groundwater recharge in carbonate aquifers from SW Mediterranean (Betic Cordillera, S Spain). *Environ. Earth Sci.*, 74, 7571–7581.

Miro, M. E. and Famiglietti, J. S. (2018). Downscaling GRACE remote sensing datasets to high-resolution groundwater storage change maps of California's Central Valley. *Remote Sens.*, 10, article no. 143.

Mohamed, A., Sultan, M., Ahmed, M., Yan, E., and Ahmed, E. (2017). Aquifer recharge, depletion, and connectivity: Inferences from GRACE, land surface models, and geochemical and geophysical data. *Geol. Soc. Am. Bull.*, 129, 534–546.

Nimmo, J. R., Stonestrom, D. A., and Akstin, K. C. (1994). The feasibility of recharge measurements using the steady state centrifuge method. *Soil Sci. Soc. Am. J.*, 58, 49–56.

Or, D. and Lehmann, P. (2019). Surface evaporative

capacitance: How soil type and rainfall characteristics affect global scale surface evaporation. *Water Resour. Res.*, 55, 519–539.

OSS (2003). *Système aquifère du Sahara Septentrional, Hydrogéologie, Volume II. Technical report. OSS, Tunis, Tunisia*, 322 pp.

OSS (2006). *Etude hydrogéologique du système aquifère de la Jeffara Tuniso-Lybienne. Technical report. OSS, Tunis, Tunisia*, 210 pp.

Polemio, M. (2016). Monitoring and management of karstic coastal groundwater in a changing environment (Southern Italy): A review of a regional experience. *Water, 8*, article no. 148.

Ramillien, G., Frappart, F., and Seoane, L. (2014). Application of the regional water mass variations from GRACE satellite gravimetry to large-scale water management in Africa. *Remote Sens.*, 6, 7379–7405.

Richey, A. S., Thomas, B. F., Lo, M.-H., Reager, J. T., Famiglietti, J. S., Voss, K., Swenson, S., and Rodell, M. (2015). Quantifying renewable groundwater stress with GRACE. *Water Resour. Res.*, 51, 5217–5238.

Rodell, M., Houser, P. R., Jambor, U., Gottschalck, J., Mitchell, K., Meng, C.-J., Arsenault, K., Cosgrove, B., Radakovich, J., Bosilovich, M., Entin, J. K., Walker, J. P., Lohmann, D., and Toll, D. (2004.). The global land data assimilation system. *Bull. Am. Meteorol. Soc.*, 85, 381–394.

Rodell, M., Velicogna, I., and Famiglietti, J. S. (2009). Satellite-based estimates of groundwater depletion in India. *Nature*, 460, 999–1002.

Rodellas, V., Garcia-Orellana, J., Masqué, P., Feldman, M., and Weinstein, Y. (2015). Submarine groundwater discharge as a major source of nutrients to the Mediterranean Sea. *Proc. Natl. Acad. Sci. USA*, 112, 3926–3930.

Sadeg, S. and Ghazali, A. M. (2008). Integrated water supply initiative for the Northwest of Libya. In *Conference: Regional Process of the 5th World Forum: Regional Meeting on Water in the Mediterranean Basin, Lefkosa, North Cyprus, 09–11 October 2008*, page 12, Lefkosa, Turkey (North Cyprus). Near East University.

Sadeg, S. A. and Karahanoglu, N. (2001). Numerical assessment of seawater intrusion in the Tripoli region, Libya. *Environ. Geol.*, 40, 1151–1168.

Saleh, H. and Rashrash, S. (2015). Status study of the Shaybah aquifer in the middle of the part of the
Jifarah plain, Libya. *Int. Res. J. Nat. Appl. Sci.*, 2, 49-61.

Scanlon, B. R., Keese, K. E., Flint, A. L., Flint, L. E., Gaye, C. B., Edmunds, W. M., and Simmers, I. (2006). Global synthesis of groundwater recharge in semiarid and arid regions. *Hydrol. Process*, 20, 3335–3370.

Scanlon, B. R., Zhang, Z., Save, H., Wiese, D. N., Landerer, F. W., Long, D., Longuevergne, L., and Chen, J. (2016). Global evaluation of new GRACE mascon products for hydrologic applications. *Water Resour. Res.*, 52, 9412–9429.

Séraphin, P., Vallet-Coulomb, C., and Gonçalvès, J. (2016). Partitioning groundwater recharge between rainfall infiltration and irrigation return flow using stable isotopes: The Crau aquifer. *J. Hydrol.*, 542, 241–253.

Seyoum, W. M., Kwon, D., and Milewski, A. M. (2019). Downscaling GRACE TWSA data into high-resolution groundwater level anomaly using machine learning-based models in a glacial aquifer system. *Remote Sens.*, 11, article no. 824.

Shaki, A. A. and Adeloye, A. J. (2006). Evaluation of quantity and quality of irrigation water at Gadoowa irrigation project in Murzuq basin, southwest Libya. *Agric. Water Manag.*, 84, 193–201.

Steyl, G. and Dennis, I. (2010). Review of coastal-area aquifers in Africa. *Hydrogeol. J.*, 18, 217–225.

Stieglitz, T. C., van Beek, P., Souhaut, M., and Cook, P. G. (2013). Karstic groundwater discharge and seawater recirculation through sediments in shallow coastal Mediterranean lagoons, determined from water, salt and radon budgets. *Mar. Chem.*, 156, 73–84.

Tapley, B. D., Bettadpur, S., Watkins, M., and Reigber, C. (2004). The gravity recovery and climate experiment: mission overview and early results. *Geophys. Res. Lett.*, 31, article no. L09607.

Trabelsi, R., Abid, K., Zouari, K., and Yahyaoui, H. (2012). Groundwater salinization processes in shallow coastal aquifer of Djeffara plain of Medenine, Southeastern Tunisia. *Environ. Earth Sci.*, 66, 641–653.

Trabelsi, R., Kacem, A., Zouari, K., and Rozanski, K. (2009). Quantifying regional groundwater

flow between Continental Intercalaire and Djeffara aquifers in southern Tunisia using isotope methods. *Environ. Geol.*, 58, 171–183.

Vernoux, J. F., Horriche, F. J., and Ghoudi, R. (2020). Numerical groundwater flow modeling for managing the Gabes Jeffara aquifer system (Tunisia) in relation with oasis ecosystems. *Hydrogeol. J.*, 28, 1077–1090.

Watkins, M. M., Wiese, D. N., Yuan, D.-N., Boening, C., and Landerer, F. W. (2015). Improved methods for observing Earth's time variable mass distribution with GRACE using spherical cap mascons. *J. Geophys. Res. Solid Earth*, 120, 2648–2671.

Weinstein, Y., Burnett, W. C., Swarzenski, P. W., Shalem, Y., Yechieli, Y., and Herut, B. (2007). Role of aquifer heterogeneity in fresh groundwater discharge and seawater recycling: An example from the Carmel coast, Israel. *J. Geophys. Res.*, 112, article no. C12016.

Wiese, D. N., Landerer, F. W., and Watkins, M. M. (2016). Quantifying and reducing leakage errors in the JPL RL05M GRACE mascon solution. *Water Resour. Res.*, 52, 7490–7502.

Yagbasan, O. (2016). Impacts of climate change on groundwater recharge in Küçük Menderes River Basin in Western Turkey. *Geodin. Acta*, 28, 209–222.

Zagana, E., Kuells, C., Udluft, P., and Constantinou, C. (2007). Methods of groundwater recharge estimation in eastern Mediterranean—A water balance model application in Greece, Cyprus and Jordan. *Hydrol. Process.*, 21, 2405–2414.

Zektser, I. S. and Everett, L. G. (2004). Groundwater resources of the world and their use. IHP VI, Series on Groundwater no. 6, unesco.org/images/0013/001344/134433epdf, UNESCO, Paris, France.

Zektser, I. S., Everett, L. G., and Dzhamalov, R. G. (2007). *Submarine Groundwater*. CRC Press, Boca Raton, USA.

Zouari, K., Fontes, J. C., Aranyossy, J. F., and Mamou, A. (1985). Isotopic and geochemical study of the soil moisture movement in the unsaturated zone in semi-arid climate (Southern Tunisia). In *Stable and Radioactive Isotopes in the Study of the Unsaturated Soil Zone (IAEA-TECDOC-357)*, pages 121–143. International Atomic Agency, Vienna.