Sustainability of fresh groundwater resources in fifteen major deltas around the world

Population growth, urbanization and intensification of irrigated agriculture in the world’s deltas boost the demand for fresh water, with extensive groundwater extraction as a result. This, in turn, leads to salt water intrusion and upconing, which poses a threat to freshwater and food security. Managing fresh groundwater resources in deltas requires accurate knowledge about the current status and behaviour of their fresh groundwater resources. However, this knowledge is scarcely present, especially for groundwater at larger depths. Here, we use three-dimensional variable-density groundwater model simulations over the last 125 ka to estimate the volume of fresh groundwater resources for 15 major deltas around the world. We estimate current volumes of onshore fresh groundwater resources for individual deltas to vary between 1010 m3 and 1012 m3. Offshore, the estimated volumes of fresh groundwater are generally smaller, though with a considerably higher variability. In 9 out of 15 simulated deltas, fresh groundwater volumes developed over thousands of years. Based on current groundwater extraction and recharge rates, we estimate the time until in-situ fresh groundwater resources are completely exhausted, partly leading to groundwater level decline and mostly replacement with river water or saline groundwater. This straightforward analysis shows that 4 out of 15 deltas risk complete exhaustion of fresh groundwater resources within 300 m depth in 200 years. These deltas also suffer from saline surface water which means their groundwater resources will progressively salinize. With a fourfold increase in extraction rates, seven deltas risk a complete exhaustion within 200 years. Of these seven deltas, six suffer from saline surface water. We stress that the groundwater of these six vulnerable deltas should be carefully managed, to avoid non-renewable groundwater use. The progressive exhaustion of fresh groundwater resources in these deltas will hamper their ability to withstand periods of water scarcity.


Introduction
River deltas are economic hotspots that are commonly densely populated and harbour lands with high agricultural productivity (Seto 2011, Neumann et al 2015. Their growing population in combination with urbanization and intensification of irrigated agriculture are a cause for rising freshwater shortages (Bucx et al 2010), which encourages increased groundwater extractions. These, in turn, lead to salt water intrusion and upconing , Michael et al 2017. Currently, salinization is already affecting many large deltas (Rahman et al 2019), leading to significant agricultural yield loss (Biswas 1993), which in turn decreases regional food security and creates new socio-economic challenges such as migration (Giosan et al 2014) and changing farming practices (Smajgl et al 2015). Moreover, a relatively high chloride content (>500 mg Cl l −1 ) in drinking water is associated with serious health risks (Khan et al 2014, Talukder et al 2017, Al Nahian et al 2018.
Groundwater management in deltas requires proper knowledge about the current status and behaviour of their fresh groundwater resources. This knowledge, however, is very limited, especially at large depths, which is caused by a lack of observations. Acquiring better insight into fresh groundwater resources in the world's deltas will pinpoint vulnerable areas and is an essential missing component to include in delta inter-comparison projects, of which several have been accomplished in the last two decades (e.g. Coleman and Huh 2003, Ericson et al 2006, Syvitski et al 2009, Bucx et al 2010, Tessler et al 2015, Van Driel et al 2015, Wolters and Kuenzer 2015. Here, we apply three-dimensional variabledensity groundwater model simulations spanning the last 125 ka to estimate fresh groundwater volumes of 15 major deltas. For every delta, these volumes are compared to estimates of present-day groundwater extraction rates, to acquire a sense of the deltas' vulnerability to groundwater extraction. The selected deltas differ greatly in size, climate, and population, but face similar issues. For instance, the thick, highly populated Ganges-Brahmaputra delta (India and Bangladesh) dwarfs the thin, scarcely populated Saloum delta (Senegal) in size and population, but both suffer from surface water salinization and are therefore highly depending on groundwater for drinking water and irrigation (Faye et al 2005, Ayers et al 2016. This work is an advancement of earlier work of the authors (van Engelen et al 2020), as sitespecific data is used as model input and to validate model results, in order to analyse real world problems, instead of being a pure synthetic study.

Background: genesis of onshore groundwater salinity in deltas
The present-day groundwater salinity distribution in most deltas presumably developed over thousands of years (Larsen et al 2017, van Engelen et al 2020, and is often not in a steady state with the groundwater systems' current forcings (Oude Essink 2001). This has been shown with model simulations (Delsman et al 2014, Larsen et al 2017, Meyer et al 2019, Van Pham et al 2019, van Engelen et al 2019, and is further illustrated by the fact that in many deltas across the globe, zones of saline groundwater are observed far landwards, which have been hypothesized or proven to be caused by a marine transgression (figure 1 and table 1), making this a global phenomenon. Nearly all of these deltas have experienced a marine transgression sometime between 8.5 ka to 6.5 ka (Stanley and Warne 1994), which serves as a reasonable hypothesis to, either fully or partly, explain the inland salinities of most deltas. The only exception to this is the Incomati delta, which did not experience a Holocene transgression but endured several Pleistocene transgressions (Salman and Abdula 1995), of which traces still can be found in its groundwater salinity today (Nogueira et al 2019). Figure 2 sketches a general picture of the evolution of the groundwater salinity distribution. During the last glacial period, the eustatic sea level was significantly lower, up to 125 m ( figure 2(a)). This led to a steep hydraulic gradients and shorelines being located further seawards than nowadays, which in turn exposed a larger area of the continental shelf to natural groundwater recharge, resulting in mostly fresh deltaic aquifers. After this period, sea levels rose rapidly leading to the Holocene transgression (Stanley and Warne 1994), during which sea water infiltrated and consequently the fresh groundwater volumes rapidly declined ( figure 2(b)). After the transgression reached its maximum extent, the fresh groundwater volumes partly recovered due to subsequent delta progradation, but not to their level during the Glacial Maximum (figure 2(c)). Finally, humans start to influence the groundwater salinity distribution (figure 2(d)), for example by extracting groundwater and land reclamation projects.

Model description
To estimate the present-day volumes of fresh groundwater in 15 major deltas, we used three-dimensional variable-density groundwater models ('paleohydrogeologic models' in short) with which we simulated the last 125 ka, encompassing a full glacial cycle. We applied the paleohydrogeologic models to these 15 deltas with an idealized geometry and lithology. With these models, groundwater salinity distributions were simulated, which aim to approach the present-day groundwater salinity distribution (figure 3). The models were constructed based on a set of 21 inputs, for which an extensive literature study was conducted to find representative values (tables A1 and A2). The inputs with the largest uncertainty and impact on fresh groundwater resources (aquifer horizontal conductivity (K h, aqf ) and aquitard vertical conductivity (K v, aqt ), as found previously in van Engelen et al 2020) were varied to account for uncertainty. We discretized the input range of K h, aqf and K v, aqt in three levels (minimum, logarithmic mean and maximum), and conducted a full factorial analysis for each delta. This resulted in nine simulations per delta.
The planar geometry of the delta was captured with a fan-shape, of which the surface elevation linearly decreased from delta apex to present day coastline. The hydrogeological base was conceptualized as a half-ellipsoid along the radial lines of the fan, with a linear decline in the coastal direction. The deterministic lithology featured a confining layer on top, underneath which a N aqt number of aquitards was positioned. Each aquitard was incised by a N pal number of paleochannels. A combination of the eustatic sea level curve (Spratt and Lisiecki 2016) and a Holocene transgression (van Engelen et al 2020) provided the boundary conditions, namely the change of the stage of, and border between, onshore and offshore surface water. These were conceptualized as assigned heads and concentrations of incoming water. Onshore surface water was conceptualized Groundwater in 32 deltas (possibly) affected by marine transgressions, classified into three discrete classes. 'Evidence' means that based on several proxies a marine transgression has been concluded in the case study in table 1. 'Hypothesis' means a marine transgression is hypothesized in the case study in table 1 to explain salinities far landwards. 'Potentially' means a large fraction of the present-day delta was covered with seawater during the Holocene transgression, which might still be reflected in groundwater salinities. The numbers refer to the entry number in table 1.
as a linear decreasing profile of assigned heads from the delta apex to the coastline. The elevation and location of the apex were constant in height and location through time, whereas those of the coastline varied through time, driven by the sea level curve and Holocene transgression. Salinity of the onshore surface water system was set with a linear salinity profile, emulating the effect of saline surface water intrusion. The models start initially completely saline, as it is assumed that the deltaic aquifers were completely  Overview of the paleohydrogeologic modelling process, in which the grey arrows indicate dataflow. A set of geometrical inputs is used to create a geometry (a). This geometry is then combined with a set of hydrogeological inputs to create a hydrogeological schematization (d). The geometry is combined with a discretized sea-level curve (c) and with data on the extent of the salinity intrusion in the surface water system and marine transgression to create boundary conditions (b). The geometry, hydrogeology, and boundary conditions are consequently used to compute a fresh-salt groundwater distribution with the computer code iMOD-WQ (e). This results in a fresh-salt groundwater distribution (f), which is halved in panel f for the sake of a clear visualization.
saline just after the marine transgression in the interglacial highstand 125 ka (Zamrsky et al 2020). This initially saline groundwater is traced with a separate tracer in the solute transport component of the model, to get a sense of the memory of the groundwater system , Delsman et al 2014. For further details, the reader is guided to the extended model description in appendix A. The simulations were computed with the iMOD-WQ code (Verkaik et al 2021), which allows for parallellization, on the Dutch national supercomputer (Surfsara 2014).

Metrics
Our simulations computed equivalent present-day fresh-salt groundwater distributions for 15 deltas under natural conditions, named 'end-state' in the rest of this paper. 'End-state' therefore refers to the results of the final timestep of the results, which represent the natural groundwater salinity distribution. In other words, the effects of groundwater extractions were not included in the model simulations. These effects, however, were corrected for in the analysis of the results, as described further. To characterize the volumes of fresh groundwater and compare between deltas we used the following metrics. Firstly, V fw,on and V fw,off which are respectively the total onshore and offshore volumes of natural fresh groundwater up to 300 m depth. Groundwater is only rarely extracted beyond this 300 m limit (Perrone and Jasechko 2017), as this is generally economically unfeasible (Wittmeyer et al 1996). Secondly, S init is the mass fraction of initial salt (initial condition of the paleohydrogeologic model 125 ka ago) still present in the entire model domain over the total salt mass. Finally, t d is the exhaustion time, which is the time until the onshore fresh groundwater is exhausted, when taking only the diffuse recharge of infiltrated rainfall into account. This was calculated as: where Q is the total yearly groundwater extraction flux of the delta and R is the yearly diffuse groundwater recharge of infiltrated rainfall in the delta. . When R exceeded Q, Q − R is set to 0, as we assumed that any excess recharge water is drained by the surface water system. The rationale behind t d is that it relates to the theoretical volume of 'clean' fresh groundwater that can be pumped. After this time, this volume is exhausted and a) partly taken out of storage, leading to head decline; b) mostly replaced by (often brackish or polluted) river water from concentrated recharge or by sea water. The definition of t d does therefore not account for enhanced recharge from rivers (Bredehoeft 2002). Moreover, t d does not include local-scale problems due to overpumping, such as saline groundwater upconing. Instead, the underlying assumption is that whenever a groundwater well salinizes, a new one is immediately created to pump in a region with fresh groundwater. To restrict understating the problem of groundwater salinization, we considered deltas with a logarithmic mean t d of less than 200 years at high risk. Deltas that cross the 200-year threshold only after increasing Q 4 times, equaling the predicted increase in groundwater extraction for the Nile Delta (Mabrouk et al 2018), were appointed a medium risk.

Validation
Validating our fresh groundwater volume estimates is challenging as data is limited. To utilize the available data as much as possible, three approaches were taken. Firstly, we compared our simulated volumes to three-dimensional distributions inferred from observations ( For each x, y location in the model, wells are dug until an aquifer is found that contains no salt in the vertical direction, up to a limit of 300 m. That is, if the upper aquifer is saline, the well filter will be placed in the first aquifer that is fresh below it, where we have used 20 m under its confining layer, as suggested in Wittmeyer et al (1996). In this way well distributions are created with depth that can be compared to the local data ( Total extraction rates were also compared between the three regional datasets and PCR-GLOBWB simulations. The total extraction rates estimated in Michael and Voss (2009a) for the Ganges-Brahmaputra were added to this comparison. Figure 4 shows a comparison between the observed and simulated depth distributions of groundwater wells. In general, most simulations were not able to capture the well depth distributions. In the Nile Delta, groundwater is mainly extracted at 50 m depth, whereas our simulated distributions seem to avoid this depth. We attribute these differences to our simplification of the confining layer. In reality, the thickness of this layer increases from a few meters , the smaller matrix of panels on the right-hand side represent the nine simulations. For each panel matrix, the columns indicate the vertical hydraulic conductivity value for the aquitards (Kv, aqt) and the rows the horizontal hydraulic conductivity of the aquifers (K h, aqf ). Table 2. A comparison between total groundwater extraction rates in regional datasets and simulated by PCR-GLOBWB. The 'factor' column states the division between these two datasets. . For the Rhine-Meuse delta, the models with a low aquifer hydraulic conductivity seemed to reproduce the observed well depth distribution the best. These simulations captured the low number of wells at 10 m depth and the peak in wells at 30 m depth.  The tail of the depth distribution was not simulated, presumably due to the simplifications in our lithological model. We also checked the total extracted fluxes in the available extraction data to the simulated PCR-GLOBWB data, shown in table 2. This table shows that considerable differences can exist. The total extraction rates in the Mekong delta were simulated quite well by PCR-GLOBWB but were underestimated in the Nile and the Ganges-Brahmaputra delta and overestimated in the Rhine-Meuse delta. Figure 5 shows example end-state fresh-salt distributions for 15 deltas. The wide variability in geometry is visible. For example, the groundwater system of the Pearl delta (No. 8 in figure 5) is significantly shallower than the thick Po delta (No. 9 in figure 5). It also shows that in some simulations, for example the Kelantan  2(b)). This delta has a lot of clay layers and paleochannels, following the conceptual model of (Griffith 2003), causing fingers of salt water to persist in the deep groundwater system after the Holocene transgression. The Kelantan delta (No. 4) clearly shows a saline upper aquifer and fresh lower aquifers just like in figure 2(c). In the Po delta (No. 9) a cone develops towards the onshore, following figure 2(d). Figure 6 presents the time series of the total volume of fresh groundwater onshore (V fw, on ) over 125 ka. Fresh groundwater volumes increase throughout the Pleistocene, after which they rapidly decline during the marine transgression that marks the start of the Holocene (note the logarithmic horizontal axis).

Fresh groundwater volumes through time
In some small deltas, namely the Burdekin and the Figure 7. End state results for four metrics. The line indicates the range of the 9 simulations conducted for each delta, the dot the logmean for V fw, on , V fw, off (the total onshore and offshore volumes of fresh groundwater in the delta), and t d (the groundwater exhaustion time) and the mean for S init (mass fraction of initial salt still present in the model domain over the total salt mass). Along the simulated V fw, on estimates are plotted the observed volumes ('observed') and estimates based on cross-sections or maps ('estimated'). The orange colour in t d indicate the effect of increasing the extraction rates with a factor four, the light brown colour decreasing these rates with a factor 2.5, similar to the uncertainties observed in table 2. The dotted line presents t d = 200 year. An infinite t d implies that the diffuse recharge of infiltrated rainfall exceeds the extraction rate.
Tokar, a full recovery from this transgression to the Pleistocene fresh groundwater volume is simulated as the delta prograded, but in most cases V fw,on does not recover. In the very thick groundwater systems (Ganges-Brahmaputra, Po, Nile, and Mekong), it seems that V fw has not reached its potential maximum value in one glacial cycle. Thus, for these systems V fw is a conservative estimate, because of our salt initial condition. Given more glacial cycles, which have relative long periods of low sea level, V fw would likely be larger in these large deltas. In most cases, the range in estimated V fw values is mainly caused by the uncertainty in the horizontal hydraulic conductivity of the aquifers, with the Nile delta as exception where the properties of the aquitards are much more uncertain (van Engelen et al 2019). Figure 7 presents the end-state metrics for the 15 deltas (see table 4 for the results in a table). Most deltas have an onshore fresh groundwater volume (V fw, on ) ranging from 10 10 to 10 12 m 3 , with the main exception being the Saloum delta in Senegal which has a very shallow groundwater system and a saline surface water system. In all cases where there are observations to validate V fw, on estimates based on simulations, The Mekong delta, the Rhine-Meuse delta and Nile delta, our the logmeans of the simulations are on the conservative side. The estimates based on reported maps and cross-sections (Ganges-Brahmaputra, Red River, Kelantan, and Saloum) are in agreement with this. The hydrogeology of Saloum delta is presumably not well captured in our model, as the simulated V fw, on has a high variance and differs a lot from the estimated V fw, on based on cross-sections. The variance in V fw, on estimates is small compared to the variance in the offshore fresh groundwater volumes (V fw, off ), which is controlled more by free-convection. Furthermore, the logarithmic means of these offshore volumes are usually at least an order of magnitude smaller than their onshore counterparts, except for the Burdekin and Kelantan delta. This is because the offshore domain of these two deltas is a lot larger than its onshore counterpart. It can also be seen that the residence times of salt in many deltas presumably exceeds that of a glaciation cycle, as S init often has a value higher than 0%. Finally, we observe that four of the simulated deltas have an expected exhaustion time (t d ) that is shorter than 200 years. The deltas with low extraction rates in the PCR-GLOBWB data (Tokar) or a high recharge (Mississippi, Po, Kelantan, Burdekin) have a large exhaustion time. Multiplying extractions by a factor four causes nine of the deltas to have a t d shorter than 200 years. The exhaustion time t d of the Saloum delta varies strongly because of its low extraction rate and strongly varying V fw, on .

Discussion and conclusion
We have conducted paleohydrogeologic reconstructions for 15 deltas across the world, varying widely in geometry, lithology, and boundary conditions. The results show that fresh groundwater resources of these areas have experienced large changes over time, especially at the start of the Holocene (9 ka), when the eustatic sea level increased rapidly. The fresh groundwater resources of two simulated smaller deltas, namely the Burdekin and Tokar, experience a complete recovery to Pleistocene freshwater volumes after the Holocene transgression. In the field, traces of the past marine transgressions are still found in these deltas (Fass et al 2007, Elkrail andObied 2013), thus the fresh groundwater volumes in these small deltas are presumably overestimated. A possible explanation for this overestimation is that a single porosity model, as used in this research, is not able to capture the tailing of a passing salt front. A dualporosity model manages to capture this phenomenon (Lu et al 2009), which would result in a slower flushing of salt water after the Holocene transgression. For the larger deltas, however, our simulations possibly underestimate the available onshore fresh groundwater resources (V fw,on ), as shown by The Red River, the Nile, Rhine-Meuse, and to a lesser degree the Mekong, Ganges-Brahmaputra and Kelantan delta in figure 7. Nearly all deltas possess onshore fresh groundwater volumes that range between 10 10 to 10 12 m 3 (2 [Red River] to 280 [Po] m3/m2) and offshore fresh groundwater volumes that are presumably one or more orders of magnitude smaller. Comparing these volumes to present-day recharge of infiltrated rainfall and extraction rates shows that four deltas have a fresh groundwater exhaustion time t d of less than 200 years and are presumably at risk, another four might be at risk with increasing extractions, and six deltas have a low risk since they receive ample recharge of infiltrated rainfall or have little groundwater extractions. The t d computed for the Saloum delta had a very high variance. Table 3 provides, based on literature, the reasons why groundwater is used rather than surface water in the considered deltas. We observe that groundwater is extracted for water quality reasons in most deltas, the Burdekin being the only exception. The most common issues with surface water are salinity, caused by either tidal differences or evapoconcentration, and nutrients, often caused by over fertilization or human waste effluents. Ten deltas experience saline surface water, and eight deltas experience an abundance of nutrients in their surface water. In only five deltas, groundwater was extracted for quantitative reasons because of surface water scarcity. We have listed whether deltas have saline surface water and their exhaustion risk (t d < 200 years) in table 4. Note that the Nile and Rhine-Meuse were appointed respectively a high and low risk, instead of a medium risk, as we corrected their extraction rates with the available field data, and thus respectively used 4.0 Q and 0.4 Q for these deltas (table 2). Furthermore, the Saloum was appointed an uncertain risk given the high variance in V fw,on . There are four deltas that have both a high fresh groundwater exhaustion risk and surface water salinity issues, namely the Ganges-Brahmaputra, the Nile, the Pearl, the Yangtze delta. The Mekong and Red River delta, also suffering from saline surface water, currently possess ample fresh groundwater resources, but are seriously under pressure with increasing extractions.
The definition of t d does not account for enhanced recharge from rivers (Bredehoeft 2002), which is likely the main source of groundwater recharge in case of intensive pumping in deltaic areas. We deliberately left this term out of the sustainability measure t d , because the usefulness of riverine recharge for freshwater supply is questionable. Groundwater is usually extracted in deltas because of an insufficient surface water quality (table 3). Riverine recharge contaminated with pollutants that break down or are adsorbed to sediments might be good enough for extraction after treatment. An example of this is the Rhine water that is first treated and then used to artificially recharge drinking water resources in the Netherlands (Stuyfzand 1997). However, when surface water is saline, riverine recharge will progressively salinize groundwater resources, gradually rendering these resources unusable. This is the case for 10 out of the studied 15 deltas (tables 3 and 4), so for these at least t d is a useful metric. A combination of a low t d and saline recharge leads to progressively saline groundwater resources on a human timescale, which means the delta's capacity to withstand periods of water scarcity decreases. This combination exists in four deltas (Ganges-Brahmaputra, Nile, Pearl, and Yangtze), and arises in two more deltas with increased extraction rates (Mekong and Red River) (table 4). Some further remarks should be made on t d . A t d of 200 years does not mean that the fresh groundwater volumes are safeguarded for 200 years. Before t d is reached, deltas already face more local-scale problems from groundwater extraction, such as enhanced salt water intrusion (Michael et al 2017) and upconing , Pauw et al 2015. In the latter case, pumping results in a rise of saline groundwater (Jakovovic et al 2016) and furthermore the induced vertical groundwater velocities will increase mechanical dispersion and cause mixing of fresh and saline groundwater (Zhou et al 2005). The latter process increases the salinity of the pumped groundwater becoming saline, making it unusable. Underneath abandoned wells, a thick brackish zone can persist for tens of years, significantly hindering the restart of these wells (Zhou et al 2005). An example of these local-scale problems can be found in the Rhine-Meuse delta, where issues with saline groundwater do exist near the coast despite it having a low risk in table 4 (Delsman et al 2018).
Another remark is that t d will presumably decrease, since extraction rates are expected to increase in most areas, as can be seen for a fourfold increase in extraction rates in figure 7. Reducing extractions, on the other hand, can increase t d by an order of magnitude. The fourfold increase of Q in figure 7 may seem large, but mind that similar increases are expected for deltas with a booming population, such as the Nile Delta (Mabrouk et al 2018). In addition, PCR-GLOBWB underestimated total extraction rates in the Nile Delta and Ganges-Brahmaputra from local data by a factor four. PCR-GLOBWB's water allocation scheme has a tendency to underestimate extraction rates in deltas, as the current version of the model does not incorporate the effects of surface water quality (Sutanudjaja et al 2018). Since in terms of water quantity, ample surface water is usually available in most deltas, groundwater is only extracted mainly because of an insufficient surface water quality (table 3). On the other hand, Table 4. End state results. V fw here stands for the onshore fresh water volume (m 3 ), the min, avg, max, obs, est suffixes respectively stand for minimum, logmean, maximum, observed, and estimated. Q and R is the extracted groundwater and recharge flux (m 3 /d). t d is the exhaution time (y), the q5 suffix indicates the exhaustion time computed using a Q multiplied with a factor 4 and 'Risk' is the risk of groundwater exhaustion.
the higher extraction rates of PCR-GLOBWB in the Rhine-Meuse delta are presumably caused by the fact that our model area does not include the ice-pushed sandy moraines. Many extractions in the Netherlands are located in these moraines, because of its good water quality (Stuyfzand andStuurman 2006, De Lange et al 2014). The uncertainty propagated by Q affected t d more severely than the uncertainty propagated by V fw,on . The Rhine-Meuse and Nile delta both showed strong error in extracted groundwater (table 2), which affected t d more than the uncertainty in V fw,on . V fw,on was estimated with an acceptable uncertainty for this study (figure 7), despite that most of the simulations did not simulate a groundwater salinity distribution which corresponds to the observed well depth distributions (figure 4). We therefore deem these simulations acceptable for this broad, general analysis, but extra work would be required to analyze smaller scale groundwater salinity issues.
Finally, we would like to stress that for 9 out of 15 simulated deltas (figures 6 and 7), fresh groundwater volumes were formed over very long timescales. The fact that mankind currently already has the capacity to exhaust the fresh groundwater volumes for six of these (Nile, Pearl, Yangtze, Yellow River, Mekong, and Ganges-Brahmaputra), undoing thousands of years of natural development within three human life spans, is concerning. The fresh groundwater resources should be carefully managed in the deltas that have poor quality surface water, as replacing good quality groundwater with poor quality surface water reduces a delta's defenses from (future) water scarcity. Reducing groundwater extractions can pay off, as halving the total extracted groundwater often increases t d by an order of magnitude (figure 7).

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: 10.5281/ zenodo.6350503.

Acknowledgments
We thank Jaap Nienhuis for helping with the estimation of tidal factors and tidal intrusion length. Furthermore, we thank Philip Minderhoud and Liduin Burgering for providing the well extraction data of the Mekong and Rhine-Meuse delta respectively. In addition, we thank Joost Delsman for providing the salinity data of the Rhine-Meuse delta. Finally, we like to thank three anonymous reviewers for their aid in the improvement of this manuscript delta_aquifer (van Engelen 2022). This research has been supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) under the New Delta programme (Grant No. 869.15.013)

Appendix A. Model description
Idealized models of 15 major deltas were created based on a set of 21 inputs, in order to conceptualize the deltas' geometry, lithology, hydrogeology, and boundary conditions. The model conceptualization was based on earlier work in van Engelen et al (2020), with two main differences. First, a full glacial cycle (125 ka) was simulated, starting from an initially saline groundwater system, instead of half a glacial cycle (40 ka) with an initially fresh groundwater system. Second, in this research, inputs were selected to create nine representative models (section A.4) for each delta, instead of the purely synthetic work in van Engelen et al (2020), which was required for the global sensitivity analysis in that study. The representative inputs are based on a literature study (section A.4), for which the values are listed in table A1. The next subsections will describe the creation of the models' components. A visual summary for all inputs that set spatial aspects of the model are summarized in figure A1.

A.1. Geometry
First, a fan shaped delta was created by setting a total length (L) and a sector angle (φ f ) in the horizontal plane. L is subdivided in an onshore length (L a ), i.e. the distance between delta apex to coastline, and an offshore length (L b ), i.e. the distance between the coastline and the foot of the coastal slope. Therefore: L = L a + L b . The inputs α, β, and γ specify respectively the slope of the onshore part, the coastal shelf and coastal slope. Finally, H a and H b respectively set the depth of the aquifer at the delta apex and the center of the coastline. The hydrogeological base is conceptualized as a half-ellipsoid across the φ − z plane.

A.2. Lithology and hydrogeology
The lithology was conceptualized deterministically. On top a confining layer was positioned, which from the coastline, reached up to l conf the distance to the apex. Underneath the confining layer, a number N aqt of aquitards were positioned. Each aquitard was incised by a number N pal of straight paleochannels, at a relative distance s pal from paleochannels in an overlying aquitard. f aqt sets the relative thickness of clay in the sediment column. The hydrogeology was specified by setting the porosity (n), the horizontal hydraulic conductivity of the aquifers (K h,aqf ), the vertical hydraulic conductivity of the aquitards (K v,aqt ), as these were the inputs usually provided in the literature, adhering to the main groundwater flow directions (Domenico and Schwartz 1990). These were converted to a K v,aqf and a K h,aqt with the anisotropy (K h /K v ). Deposition of the confining layer was simulated by converting the hydraulic conductivity of the cells in the confining layer from K h,aqf to K h,aqt during the marine regression. The longitudinal dispersivity a l was fixed to 2 m, a representative value for regional models (Zech et al 2015). These were converted to a transversal dispersivity and vertical dispersivity by multiplying a l with a factor 0.1 and 0.01 respectively (Zech et al 2015).

A.3. Boundary conditions
We used the median out of the ensemble eustatic sea level derived by Spratt and Lisiecki (2016) of the last 125 ka to set the heads offshore and determine the coastline. The code used (iMOD-WQ), requires constant sea levels in stress periods, so each stress period was set to its mean sea level. We used stressperiods of 8000 years during the Pleistocene, and refined them starting from 11 ka, as during this period the sea started to rise very rapidly (figure 3(c)). A marine transgression was forced by t tra , which is the time when the transgression reached its maximum extent, and l tra , which is the part of the onshore length L a covered by the transgression at maximum extent. Onshore surface water was conceptualized as a linear decreasing profile from the delta apex with constant height towards the coastline, which varied in height and location along with the sea level curve. Saline surface water intrusion in the surface water system was included as a linear salinity profile, which started at the concentration of sea water at the coastline, Figure A1. Visual description of symbols used for the inputs described in section A.4, for their values see table A1. The panel at the bottom right shows the boundary conditions in more detail. The 'Recharge' area is the area outside of the confining layer where recharge takes place. The 'Head, fresh' area is the onshore area with onshore surface water, the 'Head, brackish' is the area is the onshore surface water with brackish water (due to salinity intrusion in the surface water), and the 'Head, saline' is the offshore area. 'Extent transgression' shows up to where the Holocene marine transgression reaches. and linearly decreased landwards to zero at l sal . These boundary conditions together with the boundary conditions offshore consisted of assigned heads and assigned salinity of incoming fluxes. The area of the delta not covered by the confining layer received a natural groundwater recharge flux of infiltrated rainfall R, of which the magnitude varied through time. In reality, the groundwater recharge was very dynamic (Gossel et al 2010), so the constant recharge is a firstorder approximation. We deemed this as sufficient, as a global sensitivity analysis published in van Engelen et al (2020), showed only a very limited sensitivity of V fw to R.

A.4. Input data and parameterization
Our idealized models require a set of 21 inputs to create a geometry, lithology, hydrogeology, and boundary conditions. We conducted a literature study to obtain representative inputs for 15 deltas, of which the values can be found in table A1 and the literature sources can be found in table A2. Parameters to specify lateral geometry, the topography, and bathymetry were determined from global datasets (Syvitski and Saito 2007, Weatherall et al 2015, Kulp and Strauss 2019. Groundwater system thicknesses were, if available, retrieved from case studies, or, if not, from the global aquifer system thickness estimation dataset (Zamrsky et al 2018). The lithological inputs were retrieved purely from local studies (table A2), as no global dataset of groundwater system interiors exist. No lithological data could be found for the Tokar delta (Sudan), so here we used the median value of each lithological input across all deltas in our dataset. The boundary condition inputs (salinity intrusion length and marine transgression) were also retrieved from local studies, except for the cases where no data could be found for the salinity intrusion length in the surface water system (The Tokar, Kelantan, and Burdekin delta). A comparison with a global dataset of the tidal factor and tidal intrusion length (Nienhuis et al 2020) confirmed that saline water intrusion is presumably not occurring in these deltas and could be set to zero. For the hydrogeologic inputs, input ranges were usually provided in the local studies. Based on a global sensitivity analysis (van Engelen et al 2020), we varied the inputs that  create the most output uncertainty, namely the horizontal hydraulic conductivity in the aquifers (K h, aqf ) and the vertical hydraulic conductivity of the aquitard (K v, aqt ). If no input range was available for K h, aqf in a delta, a data range was taken from a global permeability dataset (Huscroft et al 2018), which provided a wide range that varies roughly three orders of magnitude. If no input range could be found for K v, aqt in a delta, an input range was created by taking the logarithmic mean of the minimum and maximum values for K v, aqt across all deltas. We discretized the input range of K h, aqf and K h, aqt in three levels (minimum, logarithmic mean and maximum), and conducted a full factorial analysis for each delta. This resulted in nine simulations per delta. The other hydrogeologic inputs were averaged over their reported range, which were the porosity, infiltration recharge, and the formation anisotropy. If data was missing for one of these four inputs for a delta, the median value of all deltas was assigned.