A three-dimensional palaeo-reconstruction of the groundwater salinity distribution in the Nile Delta Aquifer

10 The Nile Delta is an important agricultural area with a fast-growing population. Though traditionally irrigated with surface water, the delta increasingly relies on groundwater. However, saline groundwater extends far land inward, rendering groundwater close to the coastal zone useless for consumption or agriculture. To aid groundwater management decisions, hydrogeologists reconstructed this saline and brackish groundwater zone using variable-density groundwater models with very large dispersivities. However, this approach cannot explain the observed freshening of this zone as observed by 15 hydrogeochemists, who hypothesize that the coastal saline zone is the effect of the Holocene transgression. Here, we investigated physical plausibility of this hypothesis by conducting a palaeo-reconstruction of groundwater salinity for the last 32 ka with a complex 3D variable-density groundwater flow model, using state-of-the-art model code that allows for parallel computation. Several scenarios with different lithologies and hypersaline groundwater provenances were simulated, of which five were selected that showed the best match with the observations. Amongst these selections, total fresh water volumes 20 varied strongly, ranging from 1526 to 2659 km, mainly due to uncertainties in the lithology offshore and at larger depths. This range is smaller (1511-1989 km) when we consider the volumes of onshore fresh groundwater within 300 m depth. Regardless of the variance, in all cases the total volume of hypersaline groundwater exceeded that of sea water. We also show that during the last 32 ka, the total fresh groundwater volumes significantly declined, with a factor ranging from 1.9 to 5.4, due to the rising sea-level. Compared to a steady-state solution with present-day boundary conditions, the palaeo25 reconstruction improved our validation for the saline zone (5 g/L – 35 g/L TDS). Also, under highly permeable conditions the marine transgression simulated with the palaeo-reconstruction led to a steeper fresh-salt interface compared to its steadystate equivalent, while low permeable clay layers allowed for the preservation of volumes of fresh groundwater. This shows that long-term transient simulations are needed when estimating present-day fresh-salt groundwater distribution in large deltas. The insights of this study are also applicable to other major deltaic areas, given the wide-range of lithological model 30 scenarios used in this study and since many deltas also experienced a Holocene marine transgression. Hydrol. Earth Syst. Sci. Discuss., https://doi.org/10.5194/hess-2019-151 Manuscript under review for journal Hydrol. Earth Syst. Sci. Discussion started: 2 May 2019 c © Author(s) 2019. CC BY 4.0 License.


Introduction
The Nile Delta is vital for Egypt.It is an important bread basket, as it has been through history (Dermody et al., 2014), since it is the main area suitable for agriculture in an otherwise water-scarce region (WRI, 2008).Furthermore, the area is densely populated (Higgins, 2016), and the population is only expected to increase over the coming decades (World Bank, 2018).
Although the area is traditionally irrigated with surface water, agriculture is increasingly relying on groundwater (Barrocu and Dahab, 2010).For instance, a majority of interviewed farmers in the central delta indicated pumping groundwater almost continuously during the summer (El-Agha et al., 2017).This groundwater pumping will increase in the future, as large dams are planned or built upstream, e.g. the Grand Ethiopian Renaissance Dam (GERD), which will reduce the discharge of the Nile and will consequently hamper flushing of the irrigation system.This is expected to further deteriorate the quality of surface water in the irrigation system, which currently already is in a critical state; the water that currently reaches the sea is mostly saline and highly polluted (Stanley and Clemente, 2017).For example, the nitrate concentrations of the water inflow into the Manzala lagoon have increased by a factor 4.5 in the past 20-25 years, as an effect of the intensified agriculture and the reduced inflow of river water (Rasmussen et al., 2009).Given the increasing stress on the groundwater system, it is important to assess the status of the current and future fresh groundwater volumes.
The Nile Delta Aquifer (NDA) is very thick, especially near the coast, where thicknesses up to 1000 m are observed (Sestini, 1989).For reference, only 0.3% of the coastal aquifers in the world is estimated to exceed this thickness (Zamrsky et al., 2018).Despite its large groundwater volume, the amount of groundwater suitable for domestic, agricultural or industrial water supply is limited.A considerable fraction of the groundwater volume is saline (Kashef, 1983;Sefelnasr and Sherif, 2014) or even hypersaline (van Engelen et al., 2018), and unsuitable for most uses.Saline groundwater can be detrimental to Egypt's already critical food supply.For example, when increased irrigation caused the level of the brackish groundwater to rise in Tahrir, Cairo, wheat yields reduced by 41% within a four year period as soil salinity increased (Biswas, 1993).
The saline groundwater issue has led to a continuous line of research into fresh and saline groundwater occurrence in the NDA.These studies can be divided into hydrogeochemical and groundwater mechanical studies.Starting with the former, Geirnaert and Laeven (1992) provided the first conceptual palaeohydrological model (dating back to 20 ka (=20,000 years BP)).They used groundwater dating to show that the shallow groundwater was likely recharged as fresh river water around 3.5 ka.Barrocu and Dahab (2010) extended this conceptual model to up to 180 ka and concluded that the saline groundwater in the north is old connate groundwater trapped in the Early Holocene, since it is freshening.Geriesh et al. (2015) further supported this conceptual model with additional measurements.
As an example of a study based on groundwater mechanics, Kashef (1983) used the Ghyben-Herzberg relationship to show that the salt water wedge reached far inland, despite that the aquifer was gaining fresh water at that time.This was further confirmed by the 2D numerical model constructed by Sherif et al. (1988), who concluded that the width of the dispersion zone may be considerable.Sefelnasr and Sherif (2014) showed the sensitivity of the area to sea-level rise with a 3D numerical model, as the low topography of the Nile Delta allowed for a far reaching land surface inundation that can be Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2019-151Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 2 May 2019 c Author(s) 2019.CC BY 4.0 License.detrimental to fresh groundwater volumes (Ketabchi et al., 2016;Kooi et al., 2000).Mabrouk et al. (2018) showed that groundwater extraction is a larger threat to fresh groundwater volumes than sea-level rise, under the assumption that land surface inundation with seawater is prevented.Van Engelen et al. (2018) investigated the origins of hypersaline groundwater that is found starting from 400 m depth, which was overlooked in preceding studies.They tested four hypotheses of which two remained valid: 1) free convection of hypersaline groundwater formed in the Late-Pleistocene coastal sabkha deposits and 2) upward compaction-induced flow of hypersaline groundwater, formed during the Messinian Salinity Crisis.
Moreover, they showed the influence of the uncertain lithology on the groundwater salinity distribution.
When comparing the results of the hydrogeochemical studies with those of the groundwater mechanical studies, a discrepancy can be observed that is difficult to reconcile.Many of the groundwater-mechanical studies that use physicallybased variable-density groundwater flow models tend to neglect the influences of palaeohydrogeological conditions that were inferred by hydrogeochemists.Instead, they attribute the large brackish zone in the groundwater system (Fig. 1) to a large influence of hydrodynamic dispersion (Sefelnasr and Sherif, 2014;Sherif et al., 1988), with longitudinal dispersivities exceeding 70 m.The hydrogeochemical studies, however, attribute it to the Holocene transgression, since a large fraction of the brackish zone is presumably freshening (Geriesh et al., 2015;Laeven, 1991) and the extent of this brackish zone coincides with palaeo-shorelines (Stanley and Warne, 1993) (Fig. 1).Yet a variable-density groundwater flow model with present-day boundary conditions and a large amount of hydrodynamic dispersion cannot support the freshening of the brackish zone, since hydrodynamic dispersion would only lead to progressive salinization.In addition, currently no reliable field data exists that supports longitudinal dispersivities beyond 10 m (Zech et al., 2015).This discrepancy in conceptualization still exists, despite that the large influence of palaeohydrogeologcial conditions on salinities was previously shown by applying variable-density groundwater flow models for other cases (Kooi et al., 2000;Meisler et al., 1984).
The lack of including palaeohydrogeological reconstructions in variable-density groundwater flow and coupled salt transport models is understandable as it requires vast amounts of computational resources, especially when 3D reconstructions are considered.Yet these reconstructions can deliver important insights on the groundwater system, as follows from some recent examples.Gossel et al. (2010) created a large-scale 3D variable-density groundwater model of the Nubian Aquifer System and showed that sea water intrusion occurred already during the Pleistocene Lowstand towards the Qattara Depression, which lays roughly 200 km West from the Nile Delta.Later, Voss and Soliman (2013) showed with a parsimonious 3D model of the same groundwater system that groundwater tables are naturally lowering during the Holocene, since they receive limited recharge and are drained in oases or sabkhas.Moreover, the authors used an inventive validation method by comparing the position of discharge areas in the model with a dataset on oases or sabkha locations.Delsman et al. (2014) conducted a detailed palaeo-reconstruction over a cross-section in the Netherlands and showed that the system has never reached steady-state.They showed that the Holocene transgression caused substantial salt water intrusion, from which the system is still recovering.Larsen et al. (2017) showed with a combination of geophysics and 2D numerical models that during the Holocene transgression salt water preferentially intruded in former river branches in the Red River Delta, Given the stress on the groundwater system, the recent impetus of numerical palaeo-reconstructions, and the availability of a newly developed model code that allows for high-performance computing (Verkaik et al., 2017), we deem it necessary and possible to create the first numerical palaeohydrogeological-reconstruction of the Nile Delta Aquifer in 3D.Therefore, in this research we have constructed a variable-density groundwater flow and coupled salt transport model of the Nile Delta Aquifer and use it for a palaeo-reconstruction of Nile Delta Aquifer salinity over the last 32 ka.Specifically, we use the model to 1) investigate the physical plausibility of the Holocene transgression hypothesis for the Nile Delta; 2) investigate the influence of the uncertain geology; 3) provide estimates of the present-day fresh groundwater [FGW] volumes; 4) assess the importance of using palaeo-reconstructions compared to less expensive steady-state modelling.

Hydrogeology
The Nile Delta Aquifer reaches up to 1 km depth and consists of the Late-Pliocene El Wastani formation and the Pleistocene Mit Ghamr formation.These formations both contain mainly coarse sands with a few clay intercalations (Sestini, 1989).It is capped by the Holocene Bilqas formation, which has recently been mapped extensively (Pennington et al., 2017).The Pliocene Kafr El Sheikh formation is generally assumed to be the hydrogeological base as it is 1.5 km thick and consists mainly of marine clays, though it is possible that there is compaction-driven salt transport through this layer (van Engelen et al., 2018).Hydrogeological research dealing with this area generally assumed that the connection between aquifer and sea is completely open, but there exists ample evidence against this assumption from both seismics (Abdel-Fattah, 2014;Abdel Aal et al., 2000;Samuel et al., 2003) as well as from the fact that offshore borelogs contain more clay than sands (Salem et al., 2016).These observations can be explained by the fact that in the seaward direction marine clays are more commonly observed in boreholes.These marine clays can even form large vertical clay barriers when deposited during aggradational or retrogradational phases (Nichols, 2009) with a major effect on the groundwater system (van Engelen et al., 2018).

Groundwater salinity
Figure 1 shows all groundwater Total Dissolved Solids (TDS) measurements above 250 m depth that we could find in the literature (Geriesh et al., 2015;Laeven, 1991;Salem and El-Bayumy, 2016), combined with data of the Research Institute for Groundwater (RIGW) (Nofal et al., 2015).Other publications are not included as 1) these did not report measurement depths, 2) did not report measurement locations, or 3) only showed isohalines, from which the actual measurement locations are impossible to discern.The few available salinity measurements deeper than 250 m are not shown here as these are hypersaline and likely not influenced by the Holocene transgression (van Engelen et al., 2018).The different dates of the measurement campaigns in Fig. 1 matter little as groundwater ages exceeding 1000 years have been determined at similar depths at most locations (Geirnaert and Laeven, 1992).Despite this, the measurements still show considerable spatial variation, especially in the brackish zone.This variability in measured salinity can be explained with 1) the different measurement depths, 2) heterogeneity in the hydraulic conductivity of the subsoil resulting in heterogeneous salt transport, and 3) heterogeneous evapoconcentration.The latter is inferred from the observed salinities that exceed the salinity of sea water (35 g TDS/L), which presumably are pockets of evapoconcentrated Holocene groundwater (Diab et al., 1997).Despite the spatial variability in salinity, we observe a trend: the extent of the brackish zone seems to conform in the east to the coastline during the maximum transgression at 8 ka.Westwards, the extent of this zone can be explained with the location of the former Maryut lagoon, which had several periods where its salinities approached that of seawater (Flaux et al., 2013).
It is important to note that our dataset is biased towards areas with low soil salinities, as measurements have been preferentially taken in the most productive agricultural areas and thus of higher economic interest to monitor.This becomes evident from a comparison with soil salinity maps (Kubota et al., 2017).This is particularly evident from the data gathered by Salem et al. (2016b).These authors conducted research in a specific coastal area with low soil salinities, which is the area that has predominantly consisted of coastal dunes from 7.5 ka until present (Sestini, 1989;Stanley and Warne, 1993), thus providing the hydrogeological circumstances for the development of a freshwater lens.The other, former dune areas (Fig. 3) were either eroded naturally or removed by humans (El Banna, 2004;Malm and Esmailian, 2013;Stanley and Clemente, 2014) and have high soil salinities, presumably causing former freshwater lenses to salinize.Given this bias in our dataset towards fresh groundwater, the extent of the saline groundwater problem is likely underexpressed in Fig. 1.

Lithological model
A dataset of 159 borelogs was compiled of georeferenced data in different literature sources (Coleman et al., 1981;Nofal et al., 2016;Salem et al., 2016;Summerhayes et al., 1978).These borelogs were used to constrain a 3D lithological model using SKUA-GOCAD's implicit modelling engine (Paradigm, 2017).There are two main uncertainties in the lithological model thus constructed.The first uncertainty pertains to the question whether the continental slope is covered completely with low-permeable clayey sediments or not and the second whether the clay layers observed in the NDA are continuous, forming low permeable structures, or are disconnected with only a limited effect on regional groundwater flow.To account for these uncertainties, nine different lithological model scenarios were constructed (Fig. 2), where we varied the height of the clayey sediments on the continental slope and the hydraulic conductivity of the onshore-reaching clay layers.

Model code and computational resources
We applied the newly developed iMOD-SEAWAT code (Verkaik et al., 2017) which is based on SEAWAT (Langevin et al., 2008), but supports distributed memory parallelization that allows for a significant reduction in computation times.The SEAWAT code is the industry standard for solving variable-density groundwater and coupled solute transport problems, therefore the reader is referred to its' manuals for an extensive explanation (Guo and Langevin, 2002;Langevin et al., 2003Langevin et al., , 2008)).The main improvement of the iMOD-SEAWAT code is that it replaces the original solver packages for variabledensity groundwater flow and solute transport (respectively PCG and GCG) with the Parallel Krylov Solver package (PKS).
The PKS linear solver is largely based on the unstructured PCGU-solver for MODFLOW-USG (Panday et al., 2013) and solves the global linear system of equations with the additive Schwarz parallel preconditioner (Dolean et al., 2015;Smith et al., 1996) using the Message Passing Interface (The MPI Forum, 1993) to exchange data between subdomains, while keeping memory locally.The variable-density flow problem and the solute transport problem are solved in parallel using respectively the additive Schwarz preconditioned Conjugate Gradient and BiConjugate Gradient Stabilize linear accelerators (Barrett et al., 2006;Golub and Van Loan, 1996).Simulations were conducted on the Dutch national computational cluster Cartesius (Surfsara, 2014).

Model domain and numerical discretization
The study area encompasses the complete deltaic plain, the coastal shelf, coastal slope, and the Eastern desert and Western

Timeslices
The palaeo-reconstruction consisted of several consecutive time slices (Fig. 3), following Delsman et al. (2014), in which boundary conditions were kept constant.Each time slice consisted of a palaeogeographic map (Pennington et al., 2017;Stanley and Warne, 1993) that defines the location of geographical classes.Boundary conditions were assigned to five of these classes, namely sea, lagoon, sabkha, river, and dune/beach.These boundary conditions are explained more extensively in sections 3.4.2-3.4.5; a summary of their data sources can be found in table 2. The classes "clay" and "sand" indicate where the Holocene confining clay layer was respectively located and absent in each timeslice.Firstly, time slices 1 and 2 covered the Late-Pleistocene, where sea-levels were low and therefore the coastline was located ~70 km further to the North and the hydraulic gradient was high.The area consisted of an alluvial plain with braided rivers.In between these rivers, sabkha (salt flat) deposits were found in local depressions, which stayed fixed in location (Stanley and Warne, 1993).Next, time slice 3 covered the marine transgression, which was a period of rapid sea-level rise that announced the start of the Holocene.In time slice 4 the delta started prograding in the west and centre and lagoons were formed.Of particular interest is the large extent of the Maryut lagoon in the west, filling up what is now known as the Maryut depression (Warne and Stanley, 1993).In time slice 5, progradation started in the east, leading to the symmetrical shape in time slice 6.The eastwards progradation was caused by the decreased hydraulic gradient that led to finer sediments being deposited, which were then transported eastwards with the dominantly eastwards flowing sea currents.During time slice 7, humans converted most river branches to a system of irrigation channels, leaving only the Rosetta and Damietta river branches to persist.

Surface waters
The surface water systems, sea, rivers, and lagoons, were incorporated as a Cauchy boundary condition, requiring a specified head, salinity, and bed resistance, the latter being the resistance exerted by the bed sediments on the surface-groundwater interaction.Of these three inputs, the least is known about the bed resistance.A 100 day resistance was used for all surface water systems (Table 3), which is a common value for models of the Rhine-Meuse Delta (De Lange et al., 2014;Timmerman and Hemker, 1993).
Sea boundary cells were placed on the edge of the coastal shelf and slope, following the bathymetry (GEBCO, 2014) and, where possible, palaeogeographic maps that specified the coastline (Fig. 3).This meant that for the Late-Pleistocene (time slices 1 and 2) we had to recourse to the present-day bathymetry and the contemporary sea-level to approximate the coastline.We specified the head using palaeo sea-level curves.No specific sea-level curves are available for the Nile Delta (Pennington et al., 2017).Hence we resorted to eustatic sea-levels for the Pleistocene (Spratt and Lisiecki, 2016) and a more local curve for the Holocene (Sivan et al., 2001), as this was the best information available.The minimum sea-level was fixed at -100 m, as below that acquiring numerical convergence was cumbersome.This only slightly affected results as land The lagoon stage was set such that it was in hydrostatic equilibrium with the sea, that is pressure is corrected for salinity (Post et al., 2007).Lagoonal palaeo-salinities were estimated from the published strontium isotope ratios from the Maryut lagoon (Flaux et al., 2013).These salinities show strong variation through time, varying between 2-18 g TDS/L, as the inflow of the Nile varied through time.
River stages were specified by creating linear profiles from the apex to the coastline.The location of the delta apex and its river stage were fixed through time at the present-day conditions.This is a simplification, as in reality the apex' location varied through time.For instance, the delta apex was located up to 65 km south of its current location somewhere in the last 6000 years (Bunbury, 2013), but due to the lack of data on apex migration for the rest of the modelled time domain we chose this simplification.Since the coastline is irregular in shape and the location of rivers was not fixed through time, river stages were determined in the following manner: Firstly, lines were drawn from the apex to points at the coastline for each raster cell at the coastline.Along these lines, river stages declined linearly to the contemporary sea-level at the coastline.Secondly, these profiles were subsequently interpolated to a surface using Inverse Distance Weighting.Finally, the contemporary river branches were clipped out of this surface.Rivers were assigned a salinity of 0 g TDS/L, since the Nile Delta is not tide dominated (Galloway, 1975), meaning that saline water intrusion in river branches is small.

Dunes and beaches
Locations where dunes or beaches occurred were assigned a fixed recharge of 200 mm/year, equal to the present-day average precipitation near Alexandria (WMO, 2006) and also the current recharge along the Levant coast (Yechieli et al., 2010).This is reasonable, as the climate was mainly wetter throughout the Holocene than it currently is (Geirnaert and Laeven, 1992).This recharge allowed for the formation of freshwater lenses underneath dunes and beaches, which have already been observed in the area during ancient times (Post et al., 2018).

Extractions
Groundwater extractions mainly occur in the South-West.This area, near Wadi El Natrun, was appointed a reclamation area by the Egyptian Government in 1990 and since then has seen a rapid increase in extraction rates (King and Salem, 2012;Switzman et al., 2015).This was implemented in the numerical model by including extraction wells (locations and rates) that were in the RIGW database (Nofal et al., 2018) for the last 30 years of the simulation.

Hypersaline groundwater provenances
Below 400 m depth, hypersaline groundwater (HGw) has been observed in the Nile Delta Aquifer.Van Engelen et al. (2018) investigated potential origins of this hypersaline groundwater and identified two potential sources that are included in this model.Firstly, the Late-Pleistocene sabkhas could have been sources of the hypersaline groundwater.Therefore, we set the concentrations of these areas at 120 g/L for certain scenarios that are specified further in section 3.4.6,so that hypersaline groundwater could be formed in these areas.We refer to these as scenarios where HGw originates from the "top".
Secondly, another potential source of hypersaline groundwater could be seepage of hypersaline groundwater expulsed from the low-permeable Kafr el Sheikh formation due its compaction.This seepage flux was included in the model as fixed fluxes at the bottom.The magnitude and salt concentration of these fluxes were set at 6E-08 d -1 and 120 g TDS/L (van Engelen et al., 2018).These scenarios are referred to as those where HGw originates from the "bottom".

Initial salinity distributions
Initial 3D salinity distributions were set with a dynamic spin-up time for 160 ka using Late-Pleistocene boundary conditions.
Although the 160 ka exceeds a glaciation cycle, we observed that no significant changes in salinity concentrations occurred after 4 ka for the homogeneous and 80 ka for the heterogeneous cases, which is well within the range of a glaciation cycle.
For the majority of the model scenarios we started the spin-up with a completely fresh Nile Delta Aquifer, except for the model scenarios receiving HGw from the bottom, where the output of a simulation from previous research (model "NDA-c" in van Engelen et al., 2018) was used to set the initial conditions.This was a simulation of the effect of 2.5 million years of compaction-induced upward flow solute transport through the Kafr El Sheikh formation into the aquifer under interglacial conditions (low hydraulic gradient).

Model scenarios
Van Engelen et al. ( 2018) identified two sensitive inputs for hydrogeological models in predicting the distribution of hypersaline groundwater, which are the geological model and the source of hypersaline groundwater.Additionally, they found that when the hypersaline groundwater originates from compaction-induced upward salt transport, the bottom of the aquifer should be closed off.Even though there are 18 possible combinations between the nine lithological model scenarios and the two HGw provenances (top and bot), only 13 were calculated.Based on previous findings in Van Engelen (2018), three combinations could be excluded.These were the combinations where HGw was coming from the bottom and the sea boundary was completely open.The other two excluded scenarios are those with a closed sea boundary, horizontal clay layers (fluvial or marine), and HGw coming from the top.The free-convective plumes in these scenarios created entrapped volumes of fresh groundwater in between clay layers.Pressures in these zones increased fast, which, in combination with the heterogeneity, made it impossible to maintain numerical convergence for these two model scenarios, despite our best efforts.
Additionally to the palaeo-hydrological reconstruction, we formulated for each geology-HGw provenance combination, an equivalent steady-state model which used only the present day hydrological forcings until its salinity distribution reached a steady-state, comparable to previous 3D models made for the area (Mabrouk et al., 2018;Sefelnasr and Sherif, 2014).This all resulted in 26 scenarios in total.To distinguish different model scenarios, we introduced a coding system.For each feature, the corresponding letters in table 1 are converted to a code as follows: {sea}-{clayer}-{prov}-{temp}.For example, the palaeo-reconstruction of an aquifer with a half-open sea connection, fluvial horizontal clay layers, and HGw seeping in the Pleistocene sabkhas.The scenarios with a "Homogeneous" lithological model (Fig. 2) get the scenario codes "O-N".

Model evaluation
We compared our modelled salinity distributions to all the available TDS measurements (n=293), of which the majority is shown in Fig. 1.Comparing measurements with model salinities at the point location can be deceiving as sharp transition zones often occur in salinity distributions.Therefore, we assessed our model's ability to reproduce measured salinity patterns as follows.First, TDS measurements were binned into classes as shown in Fig. 1, similar to what is common in the validation of hydrogeophysical products (e.g.Delsman et al., 2018).Second, isosurfaces were determined in our model output for all bin edges with the Marching Cubes algorithm (van der Walt et al., 2014).If these classes did not correspond at the measurement location, the minimum displacement (Λ) to the observed class in the model output was determined by calculating the minimum Euclidean distance in 3D:

|
(1) where o i is the location of the observation in dimension i and   and   the locations of the isosurface vertices in dimension i for respectively the lower and upper bin edge.All locations were normalized to a range of [0,1] since the model domain was a rectangular cuboid, in other words not a perfect cube, and thus unequal in size across each dimension.Subsequently, Δ Λ was calculated to assess the improvement in validation of the palaeo-reconstructions over their steady-state equivalents: where p and s respectively denote the palaeo-reconstruction and its equivalent steady-state model.A positive Δ Λ thus means that incorporating a palaeo-reconstruction improved the fit to the observations.

Comparison palaeo-reconstruction with its equivalent steady-state model
The palaeo-reconstruction was compared with its equivalent steady-state model in two manners.Firstly, we calculated: where  is the 3D salinity matrix at the last timestep, and p and s respectively are the palaeo and steady-state model scenarios.The locations where   > 0, the palaeo-reconstruction is saltier than its steady-state equivalent and vice versa.Secondly, we calculated isohalines for a set of TDS values (3, 10, 20, 30 g/L) for both model scenarios and then looking at the distance between these isohalines in the y-dimension (North-South).As our coastline was irregular, these distances varied in the x-dimension (West-East), therefore the median (Md) was calculated across this dimension as a conservative measure of isosurface separation (ω), which for a given depth and concentration is defined as: where  are the isohaline locations in the y-dimension, and p and s respectively are the palaeo and steady-state model scenarios.

Model evaluation
Figure 4 shows the evaluation results for all palaeo-reconstructions.Note that the 4 classes above 5 g/L in Fig. 1 were grouped into 2 classes: "saline" (5-35 g TDS/L) and "hypersaline" (35-100 g TDS/L), as the number of TDS measurements was low for these classes.All model scenarios matched the location of fresh groundwater correctly (Λ = 0) at a majority of the observation points.However, their prediction quality worsens with increasing salinity.Despite having sometimes quite different salinity distributions (Fig. 6), the majority of the model scenarios predict the location of the brackish and saline zone with the same error.More striking differences are observed in the hypersaline zone, where major differences occur amongst model scenarios.Therefore, the latter class is the sole criterium in assessing model scenario behaviour.Model scenarios with Md|| < 0.07 were considered behavioural, which are the following five model scenarios: C-M-B-P, C-N-

T-P, H-M-T-P, H-F-T-P, H-N-T-P.
Figure 5 shows the improved skill of a palaeo-reconstruction compared to a steady-state model.For 11 out of 13 model scenarios, the palaeo-reconstruction improved the goodness of fit over the interquartile range of the saline measurements.
For the other two (C-M-B and H-F-T), the palaeo-reconstruction led to no clear improvement with a median Λ at zero.The palaeo-reconstruction had no clear effect on goodness of fit to fresh or brackish observation points.For four B-model scenarios with a bottom HGw provenance (except C-M-B), a steady-state simulation improved the fit, as the HGw had virtually unlimited time to flow upwards and get into steady-state with the present-day hydraulic gradient.This gradient is roughly an order of magnitude lower than the Pleistocene gradient, the starting situation of the palaeo-reconstruction.These four steady-state equivalent model scenarios required more than 60 ka to reach steady-state conditions, which is six times the extent of the Holocene period.In contrast, the C-M-B-P model scenario had a better fit to the observations than its steadystate equivalent, as the low permeable structures allowed preservation of HGw even under steep gradients.

Current spatial TDS distribution as predicted by behavioural model scenarios
The behavioural model scenarios show different salinity distributions, despite all having a nearly similar fit to the data (Fig. 6).The model scenarios without clay layers (N) or fluvial clay layers (F) have salinity distributions very similar to the O-N-T-P model scenario.In the model scenarios with marine clay layers (M), however, fresh groundwater is preserved in between low-permeable clay layers, especially in the centre.Regardless of the differences between scenarios, in all realizations there is far extending salt water intrusion visible towards Wadi El Natrun (Fig. 1).Next to the influence of this depression, more 3D patterns are visible in the salinities, especially near the (former) lagoons and the fresh groundwater lenses at (former) dune areas (Fig. 3).In general, the fresh-salt interface partly has a conical shape.

Fresh groundwater volume dynamics and sensitivity analyis
Figure 7a shows that fresh groundwater reserves have declined strongly throughout the Late-Pleistocene and Holocene.In The sensitivity analysis (Fig. 7b) shows a clear influence of horizontal clay layers.For the M-model scenarios, more fresh groundwater is maintained in the transient model scenarios than in the steady-state model scenarios.In contrast, the opposite is visible for the N-model scenarios: fresh groundwater volumes are lower in the transient simulations than in their steadystate counterparts.
In the model scenarios with marine clay layers there are still considerable fresh groundwater reserves available below 300 m depth and in one case even offshore [C-M-B-P] (Table 4).This table also shows that there is less spread amongst model scenarios when looking at the shallow onshore fresh groundwater volume.Disregarding potential deep and offshore fresh groundwater volumes decreases the spread from 74% to 32%.

Salt sources over time
Figure 8 shows the provenance of the different groundwater types as fraction of the total modelling domain.In all four Tmodel scenarios, the model domain starts initially mainly fresh, dominated by infiltrated Pleistocene river water, this water is then replaced by hypersaline groundwater.As sea-level rises, this hypersaline groundwater is in turn replaced with sea water.
The C-M-B-P model scenario takes a slightly different course, as the model is in steady-state during time slice 1, dominated by river water.As sea-level rises and the hydraulic gradient decreases, the hypersaline groundwater and infiltrated sea water volumes increase at the expense of river water volumes.Regardless of differences amongst model scenarios, the total volume of sea groundwater is lower than the total volume of hypersaline groundwater and river water in all model scenarios.
The total amount of dune water remains small for all five model scenarios over the entire modelled time.

Fresh-salt distribution of the palaeo-reconstruction in comparison with its equivalent steady-state distribution
Δ  (eq. 3 in sect.3.8) is generally negative in the marine clay (M) comparisons (snapshots 2 and 4 in Fig. 9), as the low permeable clays block vertical intrusion during the marine transgression, concurrent with previous research (Kooi et al., 2000;Post et al., 2013).In the no clay (N) comparison (snapshots 1, 3, and 6 in Fig. 9) however, a different observation is made: the shallow groundwater (< 200 m depth) is saltier, whereas the toe of the wedge of the steady-state equivalents lays further land inward.This pattern is the most clearly visible in the homogeneous model (Fig. 9, plot 1) and is obscured more as heterogeneity increases, as also visible in Fig. 10.In the C-N-T scenario, the palaeo-reconstruction is more saline nearly everywhere, as there is a larger volume of hypersaline groundwater that is influencing the toe of the salt water wedge (see also Fig. 6, snapshot 3).
In all comparisons, (former) lagoons are visible (Fig. 3 and Fig. 9).In the west, the former Maryut lagoon is clearly visible in all model scenarios as a red zone, stretching all the way towards the Wadi El Natrun depression.The central and eastern lagoons, however, are coloured blue.At these areas, the palaeo hydraulic gradients were the steepest as the lagoons lay more land inward up till 0.2 ka (Stanley and Warne, 1993), resulting in a quicker freshening of the groundwater at these locations.
Figure 11 shows the isosurface separation ω for the model scenarios where the previous described pattern was visible.An Sshaped trend in ω is visible with depth, where there is a maximum isosurface separation at around -130 m and minimum at around -350 m depth.The latter is not visible in C-N-T as the hypersaline groundwater influences the fresh-salt interface.
This trend means that, for these four model scenarios, using a steady-state model results in an underestimation of salt water intrusion length above -250 m depth and overestimation of the intrusion of the "toe" of the wedge.The median of this overand underestimation of the salt water wedge top and toe ranges up to 10 km.The isosurface seperation increases with concentration, because the denser surface has to rotate further to reach steady state.

Discussion
There is a large variance in estimated fresh groundwater volumes in the NDA, mainly located offshore and at larger depths.
This variance is mainly attributed to the limited availability of salinity observations in the saline zone.Although the number of observations may look adequate on a 2D map, it is still insufficient in 3D.Though it is generally assumed for the sake of simplicity that the northern part of the delta is completely saline, here we show that there might be just as well large overlooked quantities of offshore fresh groundwater.Still, the fresh groundwater volume estimations are on the large side.
For instance, Mabrouk et al. (2018) estimated a total amount of 1290 km 3 and Sefelnasr and Sherif (2014) 883 km 3 .We stress that the latter authors disregarded the 80 km of coastal shelf offshore and used a lower porosity value.Reducing the effective porosity from our 25% to their 17%, reduces our estimated FGw volume with 32%; disregarding the increased groundwater flow velocities this would result in, which in turn would slightly increase the extent of the fresh groundwater zone.The variance in the results should also affect management decisions.Three-dimensional groundwater model scenarios are becoming increasingly important for management decisions, as they are currently the only available tool to estimate the fresh groundwater volumes on a deltaic groundwater system scale (~200 km * 200 km * 1 km) (Delsman et al., 2014).As Egypt's reliance on groundwater increases (World Bank, 2018), more hydrogeological data should be obtained and made available to better constrain model results, which should aid management decisions.
Of all groundwater constituents, hypersaline groundwater in the NDA is the main contributor to total salt mass (Fig. 8).We can limit the amount of potential combinations between lithology and HGw provenance that were posed in van Engelen et al ( 2018), as the current study also incorporated a steep glacial hydraulic gradient, whereas van Engelen et al (2018) kept a constant low (interglacial) hydraulic gradient.Apparent from the model evaluation is that the only behavioural model where HGw originates from the bottom includes a closed off continental slope and extending, low-permeable horizontal clay layers, resulting in a compartmentalization of the bottom of the aquifer.When HGw originates from the top, also a closed off bottom of the coastal slope is necessary to preserve the observed amount of hypersaline groundwater since the last Glacial.
This implies that regardless of HGw provenance, the interaction aquifer-sea connection needs to be very limited at the lower half of the aquifer to explain the observed salinities, as otherwise all HGw would have flowed out to the sea under the steep Late-Pleistocene (32 -13.5 ka) hydraulic gradients.
Next to estimating FGw volumes, the three-dimensionality of our model showed us a vulnerable zone in the area: The Wadi El Natrun depression.This area might have started attracting, and eventually draining, sea water from the Maryut lagoon as sea-level rose above the depression height.The amount of evidence of this actually occurring is not overwhelming, but there are some observations in support of this hypothesis (Ibrahim Hussein et al., 2017).Nevertheless, it is concerning that the area with the most far reaching salt water intrusion in our model scenarios, is coincidentally in reality also the area with the most intense groundwater pumping.Even though the 30 years of groundwater extraction in our model scenarios did not seem to show a large influence on salt water intrusion, this influence can increase the coming century (Mabrouk et al., 2018), since the Nile Delta Aquifer is a slow responding groundwater system.In addition, the large cell sizes of our model can also be the upconing effects around wells (Pauw et al., 2015).In addition, it should be kept in mind for future hydrogeochemical analyses of the Wadi El Natrun area that the former Maryut lagoon used to extend a lot more land inward towards this depression than currently, and that the composition of the water in this lagoon fluctuated strongly over time (Flaux et al., 2013).
The efforts of conducting 3D palaeo-reconstructions were rewarded with new insights.We quantitatively show for the first time that a strong reduction in FGw volume presumably occurred in the Nile Delta Aquifer during the last 32 ka.As sealevel rose, FGw volumes were reduced with a factor ranging from 1.9 to 5.4 (Fig. 7).Just as in the Netherlands (Delsman et al., 2014) and the Mekong (Van Pham et al., 2019), the system never reached steady-state in the last 9000 years.Moreover, the variance in FGw volumes of the palaeo-reconstructions was larger than that of the steady-state model scenarios.This implies that conducting sensitivity and uncertainty analyses on steady-state model scenarios, can lead to grave underestimations of the uncertainty.In addition, the palaeo-reconstruction improved simulations of the location of the saline zone (Fig. 5).Furthermore, we have shown that it is physically possible that the influences of the Holocene marine transgression can be still observed in the delta.Either as fresh groundwater volumes in between clay layers in the northern part of the NDA, concurrent with earlier research (Kooi et al., 2000), or as a steeper fresh-salt interface than the steady-state interface, which to the authors' knowledge has never been shown before in this detail.The implications of this steeper interface are that using the steady-state approach can lead to an underestimation of the amount of salt groundwater above roughly 250 m depth and an overestimation of the toe length of the interface (Fig. 11).Likewise, it could explain the observed freshening of the shallow NDA as observed in hydrogeochemical studies (Barrocu and Dahab, 2010;Geriesh et al., 2015), since only the zone above 250 m depth has been sampled.It was however difficult to compare the observed freshening with model results directly, since the time scale over which this observed freshening has developed is often unknown (Stuyfzand, 2008).Which error is made in simulating salinity distributions of deltaic areas with steady-state scenarios largely depends on the lithology of their subsurface.Deltas with abundant low permeable clay layers are expected to possess larger quantities of deep fresh groundwater than would be approximated with steady-state scenarios.Examples of such deltas are the Chao Praya Delta, Thailand (Yamanaka et al., 2011), Red River Delta, Vietnam (Larsen et al., 2017), and the Pearl Delta, China (Wang and Jiao, 2012).In deltaic groundwater systems consisting of mainly coarse material, a steadystate approximation may underestimate the slope of the fresh-salt interface.Examples of these deltaic areas are the Niger Delta, Nigeria (Summerhayes et al., 1978), and the Tokar Delta, Sudan (Elkrail and Obied, 2013).Note that the Holocene transgression in the above-mentioned Asian deltas reached more land inward than it did in the Nile Delta (Sinsakul, 2000;Tanabe et al., 2006;Zong et al., 2009), thus their salt groundwater distributions are expected to more strongly deviate from from steady-state approximations.
Despite our efforts at increasing the realism of our model scenarios, some processes were only incorporated in a very simplistic manner.The most prominently simplified process in this arid area is salinization due to evapoconcentration and redissolution of precipitated salts.Though this does not currently seem to influence groundwater salinities, it is hypothesized to have been of influence in the past (Diab et al., 1997;Geirnaert and Laeven, 1992).In our T-scenarios, the Pleistocene sabkhas were assigned a fixed concentration (120 g/L), meaning that free convection here was unconstrained by available salt and water mass.Our T-scenario results indeed show a very rapid change in the first 5000 years (Fig. 7 and 8), after which they reach equilibrium conditions until the next time slice.This may seem very rapid, but experiments with smaller scale models which fully accounted for salt precipitation and dissolution, evapo-concentration, variable-density groundwater flow and the unsaturated zone, showed that higher salinities ( > 200 g/L) are possible under evaporation rates similar to those in the Nile Delta (Shen et al., 2018;Zhang et al., 2015).Given that the Late-Pleistocene sabkhas did not migrate (Stanley and Warne, 1993), we deem our T-model scenarios to be conservative.Likewise, rock dissolution is thought to influence salinities around the borders of the delta (Ibrahim Hussein et al., 2017).Either of these processes can explain the saline observations in the South-East (Fig. 1), which our model cannot.
Future research for this area would greatly benefit from an enhanced geological model, especially of the lithology of the delta shelf and slope.A lot of measurements were conducted for the petroleum industry (Sestini, 1989), but most of this data is still inaccessible for external researchers.Furthermore, the available publications on offshore geology are written mainly for the petroleum industry, focusing on the area >2.5 km depth, rendering them unsuitable for hydrogeological research of the Nile Delta Aquifer.Therefore, we think research for this area would benefit from a re-analysis of the available data with a hydrogeological perspective.In addition, a larger salinity dataset, especially in the saline zone would help validating groundwater models that are used more frequently in management decisions, e.g.Nofal et al. (2016).

Conclusions
A 3D variable-density groundwater flow and coupled salt transport model was used for a palaeo-reconstruction of the salinity distribution of the Nile Delta Aquifer.Compared to a steady-state analysis, the palaeo-reconstruction was better in reproducing salinities in the saline zone.Therefore, we deem the marine transgression as a valid hypothesis explaining the occurrence of extensive saline zones land inward.Nevertheless, the estimated FGw volumes were subject to considerable uncertainty, due to the lack of constraining data in this zone, mainly due to uncertainty in the permeability of the horizontal clay layers.Regardless, the palaeo-reconstruction provided several insights.Firstly, during the Late-Pleistocene the total FGw volumes declined strongly, with estimates ranging from a factor 1.9 to 5.4.Secondly, ignoring past boundary conditions lead to an underestimation of FGw uncertainty resulting from lack of knowledge on geological schematizations and HGw provenances.Thirdly, the differences between the fresh-salt distribution of the palaeo-reconstruction and the steady-state scenarios varied, depending on the geology.In case of low permeable clay layers, the palaeo-reconstruction resulted in considerable volumes of fresh groundwater stored underneath these clay layers.However, with more permeable or no clay layers, the fresh-salt interface was steeper than in the steady-state equivalents.Therefore, using steady-state boundary conditions is likely to result in erroneous fresh-salt distributions and an underestimation of the uncertainty in FGw volumes.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2019-151Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 2 May 2019 c Author(s) 2019.CC BY 4.0 License.Vietnam.Van Pham et al. (2019) showed that most of the fresh groundwater in the Mekong Delta Vietnam was likely recharged during the Pleistocene and preserved by the Holocene clay cap.Despite being in a humid climate, recharge to the deeper Mekong Delta groundwater system is very limited and fresh groundwater volumes are still declining naturally.
desert fringes (near the Suez Canal and Wadi El Natrun area).Its Western boundary is close to Alexandria, its Eastern the Suez Canal, its Southern near Cairo, and the Northern boundary laid 70 km offshore, resulting in a rectangle of 240 km meridional (W-E) distance and of nearly 260 km zonal (S-N) distance.This area was discretized into 1 by 1 km cells.To determine the vertical dimension, borelogs and bathymetrical/topographical data were used to respectively determine the bottom and top of the NDA.For the topography of our model, we used a combination of the Shuttle Radar Topography Mission (NASA, 2014) for the onshore part and the General Bathymetric Chart of the Oceans (GEBCO, 2014) for the offshore part.The top of the NDA was clipped off above 20 m AMSL, as the hills above this height have no important contribution to the groundwater flow.The model domain was discretized into 35 model layers, from top to bottom: 21 model layers of 20 m thickness, 10 model layers of 40 m thickness, and 4 model layers of 50 m thickness.In total, the model has more than 2 million cells.We simulated a time period of 32,000 years, covering the Late-Pleistocene and Holocene.This period was chosen as 1) no palaeogeographic maps were available of the preceding period, and 2) this research only focused on the potential effects of the Holocene transgression, not previous Pleistocene transgressions.Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2019-151Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 2 May 2019 c Author(s) 2019.CC BY 4.0 License.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2019-151Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 2 May 2019 c Author(s) 2019.CC BY 4.0 License.surface inundation at these sea-levels was the near equal due to the steep coastal slope.The salinity of the sea was set at a constant 35 g TDS/L.
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2019-151Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 2 May 2019 c Author(s) 2019.CC BY 4.0 License.from the bottom gets the following code: H-F-B-P.Furthermore, all model scenarios that share the same feature are noted as "{feature}-model scenarios", e.g."T-model scenarios" means all model scenarios where HGw originates from the top, i.e.
case of the T-model scenarios, first by free convection of hypersaline groundwater in timeslice 1, next, for all model scenarios, by sea-level rise (timeslice 2) and from 13.5 ka by the marine transgression.The total fresh groundwater volumes drop considerably with a factor ranging from 1.9 (C-M-B-P) up to 5.5 (C-N-T-P).There is also considerable variance amongst behavioural model scenario results, with the C-M-B-P having 74% more fresh groundwater than C-N-T-P.A peculiar observation for the model scenarios C-N-T-P, H-N-T-P, and H-F-T-P is that during the marine transgression (timeslice 3) the total fresh groundwater volume recovers after a quick drop.This is caused by the disappearance of the sabkhas, stopping inflow of HGw, whereas outflow of HGw continues, allowing the total FGw volume to recover (see video supplement).

Figure 1 :
Figure 1: Palaeogeography and groundwater salinity measurements up to 250 m depth of the Nile Delta.Palaeogeographical data from Pennington et al. (2017).

Figure 2 :
Figure 2: Snapshots of the lithological model scenarios that are used as input for the numerical variable-density groundwater flow model.The eastern half of the model is plotted here.The first word refers to the connection to the sea of the deeper part of the NDA, the second word to the assigned hydraulic conductivity of the onshore-reaching clay layers.

Figure 3 :
Figure 3: Timeslices with simplified palaeogeography used as boundary conditions.Figures modified after (Pennington et al., 2017; Stanley and Warne, 1993).The panels focus on the area around the present-day coastline, where upper boundary conditions varied the most.Actual model domain extended beyond the extents of these panels.In the bottom right: The median eustatic sealevel curve (Spratt and Lisiecki, 2016) and the selected sea-level for each time slice.

Figure 5 :
Figure 5: Difference in goodness of fit between palaeo-reconstruction and equivalent steady-state model scenarios.Codes indicate model scenarios.TDS values are binned in the [0, 1], [1, 5], [5, 35], [35, 100] g TDS/L for respectively "fresh", "brackish", "saline", "hypersaline".Diamonds indicate outliers, defined as values separated from the first or third quartile at 1.5 times the interquartile range.When no box is plotted for a scenario, 75% of the measurements equal zero, which consequently causes the interquartile

Figure 7 :
Figure 7: (a) the development of the total fresh groundwater volume through time for the five behavioural model scenarios.Timeslice numbers are indicated in bold with "ts".(b) End state total fresh groundwater volumes of all model scenarios, grouped for different inputs.Colours of the dots indicate a steady-state or palaeo-reconstruction result.Dots are jittered to better show overlapping dots.

Figure 8 :
Figure 8: Groundwater provenance of different water types as fraction of the total modelling domain for the five behavioural model scenarios.The dashed lines represent the C-M-B-P model.Shaded area indicates the range of the behavioural T-model scenarios(C-N-T-P, H-N-T-P, H-F-T-P, H-M-T-P), the thick line their median.Timeslice numbers are indicated in bold, preceded by "ts".

Figure 9 :
Figure 9: Comparison between the palaeo-reconstructions and the steady-state model scenarios.Everywhere red,  > , meaning that the palaeo-reconstruction is saltier than its steady-state equivalent.Vice versa for blue, where  < −.All values in between -1 and 1 are made transparent to focus purely on major differences.Groundwater below 450 m depth is not plotted, as the present-day fresh-salt interface was not located below these depths.Camera angle is the same as in the six bottom plots of Fig.5