Global Impact of Sea Level Rise on Coastal Fresh Groundwater Resources

Groundwater is the main freshwater source in many densely populated and industrialized coastal areas around the world. Growing future freshwater demand is likely to increase the water stress in these coastal areas, possibly leading to groundwater overexploitation and salinization. This situation will likely be aggravated by climate change and the associated projected sea level rise. Here, we assess the impact of sea level rise exclusively on coastal fresh groundwater resources worldwide (limited to areas with unconsolidated sedimentary systems) by estimating future decline in inland fresh groundwater volumes under three sea level rise scenarios following Representative Concentration Pathway (RCP) 2.6, 4.5, and 8.5. For that, 2D groundwater models in 1,200 coastal regions estimate the past, present and future groundwater salinity. Our results show that roughly 60 (range 16–96) million people living within 10 km from current coastline could lose more than 5% of their fresh groundwater resources by 2100 according to RCP 8.5 scenario compared to only 8 (range 0–50) million people based on RCP 2.6 scenario. We conclude that sea level rise will have severe consequences for many coastal populations heavily dependent on fresh groundwater.


Introduction
Low elevated coastal zones (LECZ) represent coastal areas with elevation lower than 10 m above sea level including areas at sea level and below sea level as in the Netherlands.At the beginning of the 21st century the LECZ were inhabited by 638 million people worldwide (Merkens et al., 2016).Following different Shared Socioeconomic Pathways (SSP) (O'Neill et al., 2014), this number is estimated to increase to more than 1 billion people by 2050 under all SSP scenarios considered, and slightly decrease (to 830-907 million) or increase (to 1.2 billion) depending on SSP scenario by 2100 (Merkens et al., 2016).Looking further landward, it is estimated that around 1.2 billion people lived within 100 km from the coastline and below elevation of 100 m at the beginning of 21st century (Small & Nicholls, 2003) indicating that the population density in the LECZ is multiple factors larger than the global average density.These coastal zones are often important commercial, agricultural and industrial hubs and as such are large freshwater consumers, frequently using groundwater as the main freshwater source (Carrard et al., 2020).Such high demand and pressure on fresh groundwater resources leads to overexploitation and can result in salt water intrusion (Custodio, 2002), see Figure 1.Additionally, changing boundary conditions due to indirect anthropogenic impacts pose another threat to the fresh groundwater resources in the LECZ.One such threat is climate change-induced sea level rise which is estimated by the latest Intergovernmental Panel on Climate Change (IPCC) reports (Fox-Kemper et al., 2023;Pörtner, et al., 2019) to reach anywhere between 0.4 and 0.7 m by 2100 and 0.8-3.7 m by 2300 comparing central values for Representative Concentration Pathways (RCPs) scenarios 2.6 and 8.5 (Van Vuuren et al., 2011) as likely sea level ranges of the indicative central estimates.While these ranges of future sea level rise urge humanity to limit further global warming, it is important to understand the scale of impending sea level rise impacts on fresh groundwater resources in the LECZ worldwide.It is projected that the immediate impacts will be in the form of increased permanent surface submergence or short-term flooding while iFoxn long-term the groundwater quality would be compromised by seawater intrusion (Nicholls & Cazenave, 2010).Given the projections of fast population and economic growth over the coming Abstract Groundwater is the main freshwater source in many densely populated and industrialized coastal areas around the world.Growing future freshwater demand is likely to increase the water stress in these coastal areas, possibly leading to groundwater overexploitation and salinization.This situation will likely be aggravated by climate change and the associated projected sea level rise.Here, we assess the impact of sea level rise exclusively on coastal fresh groundwater resources worldwide (limited to areas with unconsolidated sedimentary systems) by estimating future decline in inland fresh groundwater volumes under three sea level rise scenarios following Representative Concentration Pathway (RCP) 2.6, 4.5, and 8.5.For that, 2D groundwater models in 1,200 coastal regions estimate the past, present and future groundwater salinity.Our results show that roughly 60 (range 16-96) million people living within 10 km from current coastline could lose more than 5% of their fresh groundwater resources by 2100 according to RCP 8.5 scenario compared to only 8 (range 0-50) million people based on RCP 2.6 scenario.We conclude that sea level rise will have severe consequences for many coastal populations heavily dependent on fresh groundwater.
Plain Language Summary Sea level rise predictions for the upcoming centuries show that large strips of current coastal regions can be directly flooded.However, an often overlooked and hidden threat linked to sea level rise is the salinization of groundwater in those affected coastal regions.This can be caused by either direct lateral infiltration of saline water into the subsurface as well as by changes in groundwater pressure balance.Our modeling results show that different sea level rise magnitudes, based on future climate change scenarios, can lead to severe declines in fresh groundwater volumes in various coastal regions worldwide.
decades, LECZ areas in continents such as Africa and Asia are considered to be the most vulnerable (Nicholls & Cazenave, 2010).This combination of direct and indirect anthropogenic pressures on groundwater resources in coastal areas is called the coastal groundwater squeeze (Michael et al., 2000) and can already be observed in numerous regions around the world.Coastal flooding leading to erosion and actual disappearance of small islands around the world can already be observed (Ketabchi et al., 2014;Nicholls & Cazenave, 2010;Terry & Chui, 2012) and will lead to migration related to climate change (Warner et al., 2010).This threat is also looming over the Ganges-Brahmaputra-Meghna Delta which is already facing seawater intrusion effects on agricultural production and increased coastal flooding caused by sea level rise and exacerbated by mangrove forest removal (Faneca Sànchez et al., 2015;Khan et al., 2014;Michael & Voss, 2009).The effects of sea level rise can be further intensified by land subsidence which has been observed in many deltaic areas around the world (Herrera-García et al., 2021;Nicholls et al., 2021;Syvitski et al., 2009) for example, Ganges-Brahmaputra (Becker et al., 2020) and Mekong Deltas (Minderhoud et al., 2017).Land subsidence is caused by combination of overexploitation of coastal groundwater bodies, lower sediment deposition due to upstream dam construction and new infrastructure and construction leading to increased terrain compaction.Apart from these pressures, the need for freshwater to support human development in the LECZ is likely to continue given the projected population growth trends (Merkens et al., 2016) which calls for adaptation in water management strategies and policy development.(a) In a regional setting we recognize two sources of groundwater recharge-through surficial flow via infiltration and through potential flow from the mountains areas located inland.Sea water intrusion is often present due to human influences (e.g., pumping) or natural causes such as storm surges, we do not take these into account in this study.To quantify groundwater salinization due to future sea level rise we compare inland fresh groundwater volume change in the coastal zone stretching 10 km landwards.(b) Future sea level rise will presumably lead to increased salinization due to overtopping and seawater infiltration (1) or through increased pressure from the seaward direction (2).

10.1029/2023EF003581
3 of 24 In this study we focus on the effects of sea level rise on coastal groundwater resources and present the first global, spatially explicit quantitative analysis of sea level rise impact on coastal groundwater salinity.Our aim is to quantify the decline in fresh groundwater volumes in the vicinity of current coastline and thus identify regions worldwide that are the most endangered by effects of sea level rise on their present fresh groundwater resources.To this end we divide the global coast into 1,200 coastal sub-regions of limited length with representative coastal profiles (called sub-regional representative models (SRMs)-SRM) perpendicular to each coastal stretch.For reasons explained hereafter, we solely focus on regions with unconsolidated sediment formations.For each SRM, a two-dimensional variable-density groundwater flow and salt transport model (2D groundwater models hereafter) is built to simulate future groundwater salinization due to sea level rise projections following RCP scenarios 2.6, 4.5, and 8.0 (Pörtner et al., 2019) (See Chapter 2 and Texts S1-S4 for details).A paleo-reconstruction is conducted for each SRM by simulating the historical sea level rise since the Last Glacial Maximum (LGM) starting at 30,000 years before present (BP) (see Figure S4 in Supporting Information S1), starting with a constant sea level of −120 m below current sea level and continuing with a sharp sea level rise in the past 20,000 years.For each SRM we use three global digital elevation model (DEM) data sets (Kulp & Strauss, 2019;Weatherall et al., 2015;Yamazaki et al., 2017) to study their impact on estimated levels of groundwater salinization due to sea level rise.By including sub-regional hydrogeological uncertainty through the analysis of multiple realizations including three different DEMs, eight different geological scenarios and three RCP scenarios, our research is based on 86,400 individual model runs.It is worth mentioning that a previous investigation (Ferguson & Gleeson, 2012) studied the sensitivity of the interface between fresh and saline groundwater in the coastal zone to sea level rise.However, these results were not spatially explicit, based on simplified assumptions about groundwater salinity changes, coastal geology and groundwater dynamics, and assumed a steady-state sea level based on the maximum rise of 0.59 m predicted for 2090-2099 by the fourth IPCC assessment report.
A more recently conducted global analysis of sea level rise influence on coastal groundwater systems (Michael et al., 2013) identified two types of coastal environments based on their hydrological limiting factors.So-called recharge limited areas have higher topography that allows the groundwater table to rise in response to sea level rise and thus curb the ensuing groundwater salinization.On the other hand, topography limited areas do not have such leeway (due to a limited unsaturated zone) and thus have higher salinization rates compared to recharge limited areas.Even though such findings do neither provide local or regional quantifications nor directly help water management authorities, this classification helps to better understand the factors determining salinization risks in coastal areas.Thus, topography, as explained above (Michael et al., 2013a), is one of the main drivers influencing local and regional groundwater salinity patterns.This is especially true in low-lying coastal environments, where DEM input data set errors can have immense impact on the quality of groundwater models and sea level rise analysis (Minderhoud et al., 2019).Another important driver is local geological setting, particularly the type of rock formation in which groundwater is found and its heterogeneity (Zamrsky et al., 2018).Limiting ourselves to high-permeability hydrogeological systems that are most sensitive to salinity incursions due to sea level rise, we can distinguish two main types: karstic and porous rocks formations accommodating groundwater in complex fracture and fissure systems (Gingerich & Voss, 2005;Perriquet et al., 2014); and unconsolidated sedimentary formations composed of varying sized grained particles (Amir et al., 2013;Engelen et al., 2019).The composition of karstic systems often consists of highly variable fracture systems and is location dependent, therefore difficult to quantify on global scale without large amounts of local data.Unconsolidated sedimentary formations are a result of deposition and erosion processes governed by landward (rivers) and seaward forces.It is therefore possible to estimate and quantify geological heterogeneity characteristics in such formations on a global scale by combining available global data sets with empirical knowledge (Zamrsky et al., 2020).Thus, our results pertain to the LECZ with unconsolidated coastal sedimentary groundwater systems only, which constitute 25% of the global coast (Zamrsky et al., 2018) and currently accommodate 224 million people.

Methodology
Our global scale analysis is based on subdividing the global coastline into so-called sub-regional coastal stretches and building SRMs of variable-density groundwater flow and salt transport perpendicular to the coastline.This allows us to both implement local characteristics and increase the computational efficiency of 2D groundwater models.Sea level rise effects on local groundwater salinization are studied in more than 1,200 SRMs spanning over all continents.Their coastal extent (distance covered along the coastline) ranges from few tens to several hundred kilometers.The numerical model code SEAWAT (Langevin et al., 2008) is used to build 2D (cross-sectional) groundwater models perpendicular to the coastline for each SRM.The collected global data sets help us to define individual SRMs while taking into account local characteristics, illustrate the local hydrogeological conditions and define boundary conditions.Rapid changing sea levels in the past 20,000 years and associated climate variations had a significant impact on groundwater dynamics in coastal areas (Post et al., 2013).To properly pick up these dynamics, the 2D groundwater models follow past rapid sea level rise and climate evolution.The latter is implemented via estimated paleo groundwater recharge rates (see Text S1 in Supporting Information S1) which are calculated based on past potential evapotranspiration rates (calculated via temperature records), precipitation, land use and soil clay content.We also perform a sensitivity study to assess to what extent 2D groundwater model results are affected by spatial model resolution (Text S2 in Supporting Information S1) and the choice of the used DEM (Text S3 in Supporting Information S1).

Defining the SRMs
In our study we focus solely on coastal hydrogeological systems formed by unconsolidated sediments carried by rivers from higher elevated inland areas (sources) into these coastal lowlands (sinks).This sediment deposition has been continuously happening for millions of years and lead to a creation of coastal groundwater systems of varying thickness and geological complexity (heterogeneity).Sediment types and quantities deposited in the coastal areas depend on multiple regional factors, for example, inland lithology and elevation gradient.Our first task is therefore to divide the global coastal zone into regions whose uplands share similar geological characteristics.
To achieve this we use the so-called COSCAT regions that link the inland sources (Meybeck et al., 2006), that is, river basins, with the coastal sinks (Laruelle et al., 2013), that is, coastal plains and continental shelves.These COSCAT regions share the upland hydrogeological factors listed above, which provides an opportunity to use them to classify the coast into coastal stretches with similar upland hydrogeology resulting in a related similar coastal hydrogeology (Zamrsky et al., 2020).The average coastal hydrogeological system thickness (Zamrsky et al., 2018) is based on an approach that combines topographical profiles based on the General Bathymetric Chart of the Oceans (GEBCO) DEM (Weatherall et al., 2015) and lithological information (Hartmann & Moosdorf, 2012).
Coastal COSCAT regions can stretch over thousands of kilometers of coastline and using such large areas to derive an average representative elevation profile can potentially lead to over-simplifications.It is therefore desirable to further split the coastal COSCAT regions into smaller sub-regions.This procedure is based on establishing 2D cross-sectional profiles perpendicular to the coastline (as in Zamrsky et al., 2018Zamrsky et al., , 2020)).These coastal profiles are assigned three different coastal types based on their topographical (Weatherall et al., 2015) characteristics and potential location within a large deltaic area (in total 48 defined large deltas worldwide) (Tessler et al., 2015), see Figure 2. Coastal profiles with clearly determined onshore and offshore parts are assigned a coastal type Simple (S).If there are multiple stretches with above sea level elevations present in a coastal profile, then it is labeled as the coastal type Island (I).The Delta (D) coastal type is designated for coastal profiles that lie within a large deltaic region regardless of the topographical characteristics.
A set of rules is defined to create sub-regions in every coastal COSCAT region.First, a sub-region has to contain at least five coastal profiles placed equidistant (5 km apart) along the coastline to be considered part of our study.This rule is imposed to limit the total computation time and also to narrow down the focus on larger coastal areas.The second rule specifies the maximum gap (viz.two missing coastal profiles allowed) between individual coastal profiles in order to be considered part of the same sub-region.The third and last rule defines the maximum number (two) of coastal profiles of different coastal type allowed within a cluster of other coastal type profiles.Implementing such rule allows us to create larger sub-regions that share similar characteristics instead of generating multiple, yet almost identical, sub-regions.Figure 2 shows an example of different SRMs differentiated in a coastal COSCAT region.

Variable Density Groundwater Flow Coupled With Salt Transport Modeling
The main objective of this study is to assess the effects of sea level rise on fresh groundwater resources in coastal areas formed by unconsolidated sediments worldwide.This means evaluating the scale and rate of presumably increasing salinization with rising sea levels.To do so we used the SEAWAT numerical model code (Langevin et al., 2008) to model variable-density groundwater flow and coupled salt transport, in 2D cross-sectional profiles.The models are called 2D groundwater models throughout the following text.A 2D groundwater model is set up for each SRM (see Figure S1 in Supporting Information S1), this setup is based on previous research into offshore fresh groundwater volumes around the world which also used large-scale 2D groundwater models to estimate groundwater salinization over large geological time scales (Zamrsky et al., 2018(Zamrsky et al., , 2020)).
The chosen grid size for our groundwater models (100 m column width and 10 m layer thickness) provides almost identical estimates to finer grid size-based groundwater models while performing at much smaller runtime (Text S2 and Figure S12 in Supporting Information S1).This grid size also fits in the range of previous large-scale studies on coastal fresh groundwater salinization (Engelen et al., 2018;Feseker, 2007;Michael et al., 2013Michael et al., , 2016;;Pham, Van Geer, et al., 2019;Thomas et al., 2019;Van Camp et al., 2014;Vandenbohede & Lebbe, 2006;Vaeret et al., 2012;Xiao et al., 2018).The time resolution of our groundwater models consists of 41 stress periods (see Figure S4 in Supporting Information S1) spanning from 30,000 years BP to 300 years into the future.Including these long integration times, combined with changing paleo groundwater recharge rates, is necessary to correctly simulate the current groundwater salinity patterns, including past climate change and sea levels and particularly the early Holocene rapid sea level rise (Cohen et al., 2010;Delsman et al., 2014;Engelen et al., 2019;Gossel et al., 2010;Larsen et al., 2017;Meisler et al., 1984;Meyer et al., 2019).
Each SRM contains multiple coastal profiles based on coastal type and proximity.Once these SRMs are defined, the next step is to extract all necessary information to build the 2D groundwater models.This extraction along the coastal profiles is based on the methodology established in a previous study (Zamrsky et al., 2018) and consists of reading values from various data sets for equidistant cross-section points spaced along each individual coastal A sub-regional representative model (SRM) (SRM) will then be built for each of these 21 sub-regions.This procedure is implemented for all selected coastal COSCAT regions.
profile.An example of creating a single SRM from individual coastal profiles is given in Figure S1 in Supporting Information S1.
First, to define the SRM's extent, we extract topography and bathymetry values from the GEBCO data set (Weatherall et al., 2015) for each coastal profile within the given SRM.With its 30 arcmin resolution (approximately 1 km * 1 km at equator), the GEBCO data set is assumed to be fit for regional scale applications, and used in this study to define the SRM's mean topographical profile.Moreover, the GEBCO data set is the only global DEM that combines both onshore topography and offshore bathymetry into a single grid corrected to one vertical datum.This is important for a consistent definition of onshore and offshore domains in each coastal profile.
However, since the main focus of this study is on sea level rise impacts on onshore fresh groundwater resources, it is desirable to also include other DEM data sets with higher resolution, given that sea level rise is estimated to be in order of decimeters or few meters by 2300 (Pörtner, et al., 2019).In our study we chose to use the Multi-Error-Removed Improved-Terrain DEM (MERIT DEM) (Yamazaki et al., 2017) and CoastalDEM (Kulp & Strauss, 2019) data sets, both based on Shuttle Radar Topography Mission (SRTM) (Rodriguez et al., 2006) data but corrected for vertical errors (vegetation in MERIT DEM, vegetation and urban areas in CoastalDEM).Both MERIT DEM and CoastalDEM are corrected to mean sea level.Figure S2 in Supporting Information S1 illustrates a comparison of these three DEM data sets in an urban and a non-urban area.Even though major landscape features are captured by all the DEM data sets it is easy to spot the differences in elevation between the MERIT DEM data set (lower values) and other two DEM data sets in non-urban areas.Similarly, in an urban area, the CoastalDEM data set shows much lower elevations than the other two data sets.These differences in elevation have potentially a large effect on sea level rise impacts on fresh groundwater resources in low-lying coastal areas.

Hydrogeological Conditions
The hydrogeological schematizations used as input to the 2D groundwater models are based on a collection of global geological data sets and a randomization of unknown parameters.This approach is necessary due to lack of globally available lithological borehole data that could otherwise be used to create more accurate local hydrogeological representations.To create different hydrogeological schematizations, we use the geological heterogeneity estimation and simulation approach (Zamrsky et al., 2020).The basis in this approach is to fill the model domain with random combinations of permeable aquifers and low-permeable aquitards with varying thickness and random lithology, see Figure S1C in Supporting Information S1.
Several global geological data sets are used as direct input into this approach and help us constrain geological parameters.Amongst these are global thickness estimation of coastal unconsolidated sediment systems (Zamrsky et al., 2018), global horizontal hydraulic conductivity estimations data sets GLHYMPS (Gleeson et al., 2014) and GLHYMPS 2.0 (Huscroft et al., 2018) applied to aquitards and aquifers, presence of low permeable clay capping layer in continental shelf and slope domains (Dutkiewicz et al., 2015) and soil thickness and type SOILGRIDS data sets (Hengl, Mendes De Jesus, et al., 2014).For all models we assume a constant anisotropy of 0.1 between horizontal and vertical hydraulic conductivity.Mean hydraulic conductivity values (and standard deviation) are extracted on GLHYMPS and GLHYMPS 2.0 for each region.This mean and standard deviation is then used as input into a random function to define hydraulic conductivity of each cell in the 2D groundwater model (distinguishing between aquifer and aquitard cells), an example is shown in Figure S45 in Supporting Information S1.A previous study (Zamrsky et al., 2020) provides further information about the shape of individual aquifers and aquitards offshore and their respective volume fractions of the total estimated unconsolidated sediment volume.Other parameters are less well known and assigned a random value between reasonable boundaries, such as thickness of continental shelf and slope clay capping layer (between 10 and 30 m) and the number of aquifer and aquitard combinations in the model domain (between 2 and 5).The clay content of the aquitards is related to the sand/mud ration of the upstream sediment supply.Assuming that small size sediment particles such as clay are deposited during periods with higher sea levels, low permeable clay cells are randomly but preferentially assigned to the upper parts of each aquitard.Also, dependent on the upstream sediment supply and the subsidence rate that determine preservation potential, an erosion factor is determined and applied which results in removing a certain amount of these clay cells in random locations (Zamrsky et al., 2020).The (random) clay stacking factor (between 0.5 and 1) determines where clay cells are likely to be positioned within an aquitard in vertical direction.Higher clay stacking factor values mean that a clay cell is more likely placed in the upper part of an aquitard.In such cases, the aquitard will act as flow barrier compared to situations where clay cells are randomly distributed throughout the aquitard and potentially leaving space for permeable groundwater flow conduits.

Boundary Conditions in the 2D Groundwater Models
Once the spatial dimensions of a SRM are defined (see Figure S1 in Supporting Information S1), we assign hydrological boundaries to the 2D groundwater model domain.The bottom boundary is assumed to be impervious.Since our approach only takes into account unconsolidated sediment groundwater systems.This is probably an oversimplification in some cases since there might be consolidated sedimentary aquifers underneath which would mean a possible groundwater flow connection.However, simulating groundwater flow in such consolidated systems would require a large amount of local geological data (that are not available on a global scale), long simulation times and a different modeling approach all together.
A general head boundary (GHB) is assigned at the landward extremity of 2D groundwater model to simulate groundwater flow from the landward direction.The head and conductance of the GHB is chosen such that these cells are representing a (fixed) head elevation equal to the topographical elevation in given 2D groundwater model column and a fresh water concentration (= 0 total dissolved solids [TDS] g/l).The landward boundary location corresponds to a groundwater flow divide that is estimated based on the DEM elevation input.If the groundwater flow divide is not located within the first 50 km inland a GHB is placed 50 km landward from the coastline.The rest of the inland domain (where cell elevations are above sea level) receives fresh groundwater recharge and is also assigned a drainage system to capture exfiltrating groundwater in case the phreatic surface comes close to the surface.The drain elevation is equal to surface level and is implemented to avoid overestimating the forcing of fresh water into the model and thus artificially increase groundwater heads in the model domain.The upper most active cells in each offshore column (depending on fluctuating sea level) are assigned the GHB condition with head elevation equal to sea level and a seawater concentration (= 35 TDS g/l).The same conditions apply to the last active offshore column of the model domain.
The LGM sea level conditions are chosen to be the starting paleo conditions for the 2D groundwater models.This time period occurred approximately 20,000 years BP and the average global sea level is estimated to be −130m compared to current sea level (Lambeck et al., 2014).Initial salinity concentration (TDS g/l) and groundwater heads (m bsl.) profiles for an SRM is derived based on (i.e., interpolated from) the LGM conditions of an averaged representative profile 2D groundwater model of the corresponding COSCAT region (Zamrsky et al., 2021), see Figure S4 in Supporting Information S1.This profile follows from running the representative profile 2D groundwater model from the last high stand 125,000 years BP with completely saline initial conditions and is taken from Zamrsky et al. (2021).
The LGM period was followed by rapid sea level rise over the next 20,000 years until sea levels stabilized to current conditions.To simulate this process we divide this whole time period into separate stress periods each with different sea level approximated from the sea level rise estimates (Lambeck et al., 2014), see Figure S4 in Supporting Information S1.An extra stress period (SP0, 10 ka long, from 30,000 to 20,000 years BP) is implemented at LGM conditions to allow the 2D groundwater model to adapt the interpolated groundwater salinity initial conditions to the SRMs boundary conditions.The temporal resolution of following stress periods (SP1-SP18) is 1ka until sea level stabilized at current condition around 2000 years BP.From then onwards the stress periods are 100 years long (SP19-SP41), to take into account varying groundwater recharge, as is explained above.To finish the 2D groundwater model simulations, we implement three stress periods with estimated global mean sea level (GMSL) rise conditions (Pörtner, et al., 2019) based on three different RCP (RCP) scenarios (2.6, 4.5, 8.5 taking the GMSL 50th percentile values).The reason we decided to use the GMSL data set instead of the regional sea level rise predictions is its temporal span that provides estimates until 2300 compared to 2150 for the regional estimates (Fox-Kemper et al., 2023) for the three chosen scenarios.The implemented sea level rise levels in the 2D groundwater model simulations labeled as RCP 2.6, RCP 4.5 and RCP 8.5 correspond to indicative central estimates SSP1-2.6,SSP2-4.5 and SSP5-8.5 (and SSP3-7.0)through to 2100 respectively as indicated by the estimates of regional sea level rise (Fox-Kemper et al., 2023).The sea levels set for the final stress periods ending at years between years 2100 and 2300 are in accordance with probable regional sea level rise estimates for scenarios SSP1-2.6 and SSP5-8.5.The highest sea level rise considered in our models is 3.7 m in year 2300 for the RCP 8.5 scenario (see Figure S4 in Supporting Information S1).This value corresponds to the latest most likely sea level predictions (Fox-Kemper et al., 2023) not accounting for ice-sheet instability and close to SSPs 2 and 3.As such, the worst case sea level scenario (called RCP 8.5) considered in our 2D groundwater models can be interpreted as the most likely worst case scenario based on the most recent sea level rise estimates (Fox-Kemper et al., 2023).10.1029/2023EF003581 8 of 24

Inland Fresh Groundwater Volume Decline
Sea level rise effects on future fresh groundwater volumes of the LECZs are analyzed for each SRM by implementing a set of 2D groundwater models as described above.We quantify these effects by first estimating the current inland fresh groundwater volume (IFGV) from paleo groundwater modeling (see Section 2) and then comparing future IFGV projections among three RCP scenarios (Pörtner et al., 2019).We present the relative rate of change of IFGV in the future (expressed as decline in IFGV in year 2100 compared to situation in year 2000) rather than absolute IFGV change.This allows us to better identify and compare the threatened areas around the globe.The IFGV is calculated as the total volume of fresh groundwater in model cells that contain fresh groundwater and located in the area spanning up to 10 km from the current coastline position in the landward direction.Freshwater salinity concentration is defined as a value lower than and equal to 0.5 g/l TDS.The distance of 10 km is chosen arbitrarily and there is no depth limit in the vertical direction, we consider the whole model domain thickness for the IFGV change analysis.

Sea Level Rise Effects on Future Fresh Groundwater Volumes
Figure 3 shows for section SRM 1415_32 the impact of sea-level rise on coastal salinization.It shows that particularly for RCP 8.5 salinization of coastal groundwater is considerable by 2100 and for both RCP4.5 and RCP 8.5 from 2200 onward.Global maps of sea level rise impacts on coastal fresh groundwater volumes in the first 10 km from the coastline are provided in Figure 4 for each RCP scenario considered.Differences between RCPs are quite large; the estimated IFGV decline being severe for RCP 8.5 for several regions (e.g., East coast USA, Caribbean, Northern Java) that potentially lose more than 25% of their IFGV compared to the year 2000.Contrarily, future sea level rise according to the RCP 2.6 scenario shows only a very minor decline in IFGV by the end of the 21st century with the vast majority of SRMs only experiencing a IFGV decline lower than 5% compared to situation at the start of the century.End of 21st century estimates following RCP 4.5 show intermediate results with most SRMs experiencing a decline in IFGV lower than 5%, while several local hotspots (e.g., West coast Africa, Mekong Delta and the Netherlands) are visible where the decline in IFGV is larger than 25%.
Our global projections shows that a steady decline in IFGV worldwide will continue well after the year 2100 (see Figures S25-S36 in Supporting Information S1).The differences between the outcomes for the three RCP scenarios become more pronounced over time.Based on the RCP 2.6 scenario, sea level rise is projected to level off after year 2200 with only 0.3 m rise during the 22nd century and 0.1 m during the 23rd century (see Figure S4 in Supporting Information S1).These limited rates of sea level rise show little impact on the worldwide IFGV (see Figures S25-S28 in Supporting Information S1) as estimated by the 2D groundwater models as compared to the estimated situation in year 2100 shown in Figure 4a.By the end of the time period considered in our analysis in year 2300 (see Figure S28 in Supporting Information S1), the decline in IFGV worldwide under RCP 2.6 scenario is almost identical to the decline in IFGV projected to happen by 2100 under RCP 4.5 scenario (see Figure 4b).Projected sea level rise under RCP 4.5 follows an almost linear trend of 0.5 m per century until year 2300 (see Figure S4 in Supporting Information S1).This will translate into a steady decline in IFGV until year 2300 (see Figures S29-S32 in Supporting Information S1), at which IFGV decline to matches that projected in 2100 under RCP 8.5 (see Figure 4c).However, the majority of SRMs representing the global coastline would still experience only a minor decline in IFGV (<5%) in the coming centuries.
As anticipated, the largest decline in IFGV over the coming centuries is projected if sea level rises following the RCP 8.5 scenario, with a sea level i.e., 3.7 m higher than current level in 2300 and the rates of sea level rise increasing over the subsequent centuries: 0.7 m during the 21st, 1.4 m during the 22nd and 1.6 m during the 23rd century, see Figure S4 in Supporting Information S1).Such high rates of sea level rise create large hydraulic gradient changes and seawater overtopping in coastal areas which will lead to substantial declines in IFGV over the coming centuries (see Figures S33-S36 in Supporting Information S1).The projected situation at the end of 21st century, with several regions are already badly affected but with a majority of SRMs still showing small declines in IFGV, steadily worsens through the 22nd and 23rd centuries.By the end of 23rd century additional regions show up where declines in IFGV are larger than 25% compared to the year 2000.Among these regions are for example, the coast of West Africa, North-east Indian coast and the coastal areas in the Red Sea.These are densely populated and currently relatively impoverished areas and therefore probably not economically ready to 10.1029/2023EF003581 9 of 24 deal with the consequences of decreased fresh groundwater supplies (unless their economic conditions do not improve in future).
The salinization of coastal aquifers resulting from sea level rise can also be is observed by considering the landward shift of saline groundwater wedge (SWW) in Figure 5 (and Figures S37-S42 in Supporting Information S1).For the RCP 8.5 scenario, already by the year 2100 several coastal areas could experience groundwater salinities higher than the salinity level half of ocean water (17.5 g/l TDS) 10 km further inland than the situation in year 2000.Especially areas in the eastern USA, northern Europe and western Africa could experience significant aquifer salinization extents.Similar to the trend of decline in IFGV, the estimated shift in landward groundwater   salinity extent is significantly lower for the RCP 2.6 scenario and the difference between RCP 8.5 and 2.6 scenarios widens into the future.

Elevation Impacts on IFGV Estimates
The impact of sea level rises per RCP scenario for a given SRM is highly elevation dependent.Areas with elevated coastal elevation profile (Figures S13-S18 in Supporting Information S1) show almost no changes in IFGV for all RCP scenarios.This is due to the compensating effect of rising groundwater levels -which are possible due to the high surface level elevations-and associated increased fresh groundwater volumes following the Ghijben Herzberg principle, provided of course that recharge or freshwater influx from the inland hydrological boundary are sufficiently large.Such areas correspond to the earlier defined recharge-limited coastal areas and the rather negligible impact of sea level rise matches the findings of a previous study (Michael et al., 2013).On the other hand, as expected, the 2D groundwater model simulations show that coastal areas with low lying topography (i.e., topography-limited) experience a much more severe seawater intrusion due to sea level rise than recharge-limited areas.Naturally, large-scale overtopping (several kilometers up to tens of kilometers) by seawater due to sea level rise does not occur in topography-limited areas with elevation higher than projected sea level rise rates.In such cases the seawater intrusion tends to be limited to a first few kilometers at the coastline (see Figures S19-S21 in Supporting Information S1).This can be explained partly by low elevation of the first few hundred meters leading to small-scale overtopping but mostly by decrease in hydraulic gradient which leads to a landward shift of fresh-saline groundwater interface.If changes in the hydraulic gradient are only minor this shift can only be limited to few hundred meters depending on the RCP scenario.
The effects of sea level rise on fresh groundwater volumes are the most visible in topography-limited coastal areas with very low-lying elevation profiles (i.e., lower than projected sea level rise) stretching over a large portion of the inland coastal zone.The sea level rise magnitude clearly impacts the magnitude of these effects with large differences between RCPs.Where RCP 2.6 seem to have negligeable effects (see Figure S22 in Supporting Information S1), RCP 4.5 effects are clearly visible by year 2300 (see Figure S23 in Supporting Information S1) and RCP 8.5 clearly impacts the fresh groundwater volumes already during the 21st century (see Figure S24 in Supporting Information S1).The hydraulic gradient in these topography-limited coastal areas is already quite low under current conditions and therefore very sensitive to any variation due to sea level rise.In combination with seawater overtopping of the inland areas and the limited accommodation space to increase fresh groundwater levels without being drained by surface waters, this leads to a seawater intrusion spanning over several kilometers (see Figures S23-S24 in Supporting Information S1).Variation in projected seawater intrusion severity is also observed amongst 2D groundwater models with a different DEM input under the same RCP scenario (see Figures S23-S24 in Supporting Information S1).A few decimeters or meters change in elevation can already have significant impacts on estimated seawater intrusion in topography-limited coastal areas with very low elevation profile.This is in line with previous research conducted in low-lying areas (Minderhoud et al., 2019), and stresses the need to include the uncertainty in elevations when projecting sea level rise effects on groundwater salinity in coastal zones.

Impacts on Coastal Communities and Economies
The overall evolution of sea level rise impacts on IFGV decline in the coming centuries is presented in Figure 6.As already mentioned above, a clear difference in the severity of seawater intrusion (leading to IFGV decline) is observed between the RCP scenarios considered in this study.If the future sea level rise follows the RCP 2.6 scenario, less than 10% of the SRMs would experience a decline in IFGV by the year 2300 and for more than half of the affected areas the decline would be limited to 5%-10%.Under the RCP 2.6 scenario, human societies inhabiting the coastal areas would likely have enough time and resources to design adaptation measures to tackle the impacts of sea level rise on fresh groundwater volumes.Figure 6 shows that global adaptation efforts would have to be much larger under the RCP 4.5 and RCP 8.5 scenarios.For RCP 4.5 approximately 10% of the SRMs are projected to lose more than 5% of their IFGV by the end of 21st century and around 5% of the SRMs a IFGV decline larger than 10%.Under RCP 8.5 the fraction of SRMs losing more than 5% of the IFGV around 2100 is about 20% and more than 10% of the SRMs show a decline of more than 5%, with a non-negligible fraction of SRMs showing a decline larger than 25% reaching up to 50%.We note that these numbers are averaged over the results for three individual DEMs (Kulp & Strauss, 2019;Weatherall et al., 2015;Yamazaki et al., 2017).We have analyzed the impact of DEM input on IFGV decline for given scenarios (text S3 and Figure S43 in Supporting Information S1) and found results between DEMs to be limited to a few percent.
Based on global gridded population count for year 2020 (CIESIN, 2017) around 224 million people currently live in close proximity to the current coastline (closer than 10 km) in the SRMs considered in this study, see Table 1.In the analysis below the estimates are the mean decline (across 8 hydrogeological and 3 DEM realizations-24 runs) in IFGV for each SRM.A best and worst case scenario is also calculated by taking into account the standard deviation (across 24 runs) of the calculated decline in IFGV.This means that in the best case scenario this standard deviation is subtracted from the mean decline in IFGV and vice versa for the worst case scenario.By the end of the 21st century, taking the uncertainty in hydrogeology and elevations into account, we estimate that almost 60 (16-96) million people could be directly affected by limited access to fresh groundwater (decline in IFGV  between 5% and 50%) under the RCP 8.5 scenario sea level rise projections.Approximately one sixth of the 60 (19-96) million people is projected to live in severely impacted areas with decline in IFGV larger than 25%.In contrast, less than 8 (0-50) million people would experience a decline in IFGV ranging between 5% and 10% if sea level rise follows the RCP 2.6 scenario projections.The differences between RCP 2.6 and RCP 8.5 increase further when projecting further into the future with 16 (1-76) million and 120 (82-151) million people affected by 2300 respectively.The number of people living in severely affected areas by 2300 (decline in IFGV larger than 25%) is estimated to be only 0.1 (0.1-13) million for the RCP 2.6 scenario but more than 45 million (20-72) (one third of total population affected) for the RCP 8.5 scenario.
As a measure of economic costs we calculated the total gross domestic product (GDP) for year 2015 in constant international 2011 USD (Kummu et al., 2018) generated in the affected coastal regions (closer than 10 km to current coastline).Overall, the total yearly GDP produced in these areas is around 4,625 billion USD per year, see Table S1 in Supporting Information S1.Our estimates show that under RCP 2.6 in year 2100, an area with GDP of around 41 (0-785) billion USD per year will be affected in coastal areas with a IFGV decline of 5% or larger.None of these coastal areas will experience a decline in IFGV above 25%.By 2300, under the same scenario, the total GDP in affected areas with an IFGV decline larger than 5% rises to 214 (11-1,334) billion USD per year.However, only 3.3 (0-100) billion USD in yearly GDP would be located in severely impacted coastal areas where decline in IFGV would surpass 25%.Economic losses are expected to be far higher under more severe future sea level rise projected by the RCP 8.5 scenario.The total yearly GDP in affected coastal areas experiencing a decline in IFGV larger than 5% would reach 1,123 billion (201-1,733) USD by year 2100, out of which 75 (28-600) billion USD is located in coastal areas with decline in IFGV above 25%.By 2300, total yearly GDP in affected coastal areas with an IFGV decline larger than 5% than doubles compared to 2100, reaching 2352 (1,600-3,145) billion USD per year, of which 1,016 (374-1,530) billion USD is located in coastal areas with more than 25% decline in IFGV.

Validation and Impact of Dimensionality
We compared our regional modeled groundwater salinity estimates with a collection of measurements derived from local studies in order to evaluate our model results.The maximum inland saline groundwater extent is chosen as a metric since the collected local studies often provide limited information on other specifications (e.g., depth of saline groundwater).Table S2 in Supporting Information S1 and Figure 7 show the observed inland extent and the mean simulated extent (standard deviation in brackets in Table S2 in Supporting Information S1).
Given that we are comparing the mean (a representative estimate) of 24 alternate realizations of salinity distributions with a local study, we can only expect to have an order of magnitude resemblance this way.Based on the mean profiles we get a mean absolute error (MAE) of 4.5 km for non-delta coastal regions with coefficient of correlation (CC) equal to 0.46 and coefficient of determination (CD) of 0.21.When we look at the 95th percentile the results somewhat improve to MAE of 3.1 km, CC of 0.55 and CD equal to 0.3.Expectedly, our estimates in delta regions are worse with an overall MAE of 27.3 km, CC of 0.055 and CD of 0.003.This can be explained by quite a large number of outliers (see Figure 7).When taking into account the 95th percentile of our estimates we can see an improvement in MAE (23.8 km), CC (0.3) and CD (0.09).
The outliers indeed show large discrepancies (see Figure 7) and could be both caused by our 2D models not being fully fit for simulating groundwater salinities in complex environments or by not included emplacement of saline water during sea transgressions or more recent flooding by saline water due to storms or tsunamis.Also, we did not include groundwater pumping, which could increase the extent of salinity by upconing and reducing recharge to the sea.However, even though the models have difficulty in predicting the absolute salinity distribution at a given location, we believe that they contain sufficient regional (spatial) specificity to perform regional (spatial) sensitivity analyses and regional relative changes in fresh groundwater volumes.Particularly because differences in geological settings (subsidence rate, sediment influx, sand/mud ratio) and topography are very different between regions (see Zamrsky et al., 2020).
Furthermore, we built a regional three-dimensional groundwater model (where we also included groundwater extractions) to examine the differences in estimated decline in IFGV between our 2D groundwater models and a full 3D simulation (see Text S4 in Supporting Information S1).The results suggest that the differences in IFGV decline are within order of single percent units at most while the groundwater extraction had a very limited effect on the estimated decline in IFGV (see Tables S3-S6 and Figures S46-S43 in Supporting Information S1).

Discussion
The objective of this study is to analyze the impacts of sea level rise on fresh groundwater resources in regions with unconsolidated sedimentary systems along the global coastline.The novelty of this study lies in building a large set of regional scale 2D groundwater models based on state-of-the-art global data sets, geological heterogeneity simulations as well as conducting a paleo-reconstruction to simulate past sea level changes.In the presented analysis we opted to determine these impacts by calculating the decline in IFGV rather than providing quantitative IFGV estimates.We stress that this is a global study using global data sets and parametrizations that, albeit specific for the coastal regions identified, may deviate from local salinity measurements or local groundwater models with higher complexity, local input data input and grid resolution.Thus, expressing the future trend in IFGV decline in percentage change compared to the estimated situation in year 2000 ensures that we can identify the potentially most threatened regions by sea level rise (and its impact on fresh groundwater volume) in the context of limited local data and information.Our results show that the future state of fresh groundwater resources in coastal zones considered in this study varies dramatically between different RCP scenarios.The worst-case scenario RCP 8.5 would lead to around 20% of all SRMs experiencing reduced IFGV volumes by 2100, rising to almost 50% by year 2300.In comparison, the fraction affected SRMs under RCP 2.6 would be only around 1% by 2100% and 7% by year 2300.Coastal areas with low topography will be the most threatened, corroborating results of a previous global sensitivity study (Michael et al., 2013).These coastal areas with low topographical gradients also have low hydraulic gradients and are shown to be more prone to increased salinization as result of sea level rise (Ferguson & Gleeson, 2012).In areas with higher hydraulic gradients, defined as 1 m groundwater vertical head difference per 1,000 m profile length (Ferguson & Gleeson, 2012), sea level is deemed to have only very limited impacts on IFGV.In our study we observe the same trend with relatively low IFGV decrease on global scale by the end of 21st century.However, if sea level rise continues and even accelerates over the next centuries (RCP 8.5) the negative impacts will be much more severe and will affect even areas with higher hydraulic gradients.S2 in Supporting Information S1 lists all the collected studies).The marker size represents the frequency of discrepancies between observed and estimated values as some studies show the same extent of salinity inland while our estimates in these sub-regional representative models can also show identical values.

Groundwater Model Simplifications and Their Implications
For simplicity, we used a fixed but spatial varying future groundwater recharge (Figure S11 in Supporting Information S1) for all RCP sea level scenarios.The groundwater recharge used is on the high side (globally 20,963 km³/ yr) when compared to other global studies, with values ranging from 12,666 to 17,000 km³/yr (Sutanudjaja et al., 2018).On the other hand, the estimated average yearly groundwater recharge rate (142,8 mm/yr) is only slightly different or lower from recent global studies with global averages of 134 mm/yr (Mohan et al., 2018) and 218 mm/yr (Berghuijs et al., 2022).However, our groundwater recharge estimates seem to underestimate mean global groundwater recharge by almost 50% when compared to a study based on collected global groundwater recharge measurements (Moeck et al., 2020).The used paleo groundwater recharge reflects the climate changes that occurred since the LGM (Lone et al., 2018), showing a wetter period followed by a dry and arid period between 6,000 years BP and 3,000 years BP, as well as regional studies dealing with paleo groundwater recharge (Mabrouk et al., 2018;Wu et al., 2020).
The rapid sea level rise during the LGM also potentially involved deposition of saline water inland in some regions during the Holocene transgression which we did not take into account.In this way, our regional groundwater salinity estimates are probably overestimating the fresh groundwater volume in these areas (see Figure 7).However, since we focus on decline in IFGV compared to year 2000 it only mildly affect our results.Moreover, our models do not take into account gradual sedimentation processes that happened since the LGM as the geological conditions in our models are constant over the whole duration of the simulation.Implementing these processes could lead to more accurate regional estimates of groundwater salinity.However, our study only focuses on comparing the current and future rates of IFGV and not actually presenting the IFGV itself and therefore, similarly, as with inland flowing during the transgression, it can be expected that these past sedimentation processes wouldn't have substantial impacts on the estimated decline in IFGV.Due to our attention to the effect of sea level rise while limiting the number of climate runs, we did not consider future changes in groundwater recharge.However, as changes in future groundwater recharge may be considerable between different regions (Wu et al., 2020), it will be of interest to include it in future impact studies.
A recent literature review dealing with sea level rise effects on seawater intrusion in coastal aquifers concluded that it is crucial to move away from hypothetical or highly simplified coastal aquifer representations toward more complex hydrogeological representations (Ketabchi et al., 2016).In our study we navigate around missing local geological data (e.g., boreholes) by implementing a semi-randomized generation of hydrogeological schematizations (reflecting sub-regional coastal geological settings).This approach is based on several key global geological data sets (Gleeson et al., 2014;Hartmann & Moosdorf, 2012;Huscroft et al., 2018;Montzka et al., 2017) and an estimation of regional geological heterogeneity conditions (Zamrsky et al., 2018(Zamrsky et al., , 2020)); further details are provided in Chapter 2.3.
In accordance with aforementioned review (Ketabchi et al., 2016), the 2D groundwater modeling approach also includes past sea level rise (over half a glacial period) and topographic slopes based on three different global DEM data sets.Although this approach brings our large-scale analysis closer to actual regional conditions, there is still a gap to close to arrive at accurate sub-regional to local estimates.This gap could be closed by including local geological data (once available on global scale) into our hydrogeological simulations, as well as other local drivers such as pumping wells and river systems.These drivers were not included in our current 2D groundwater models due to both lack of information on global scale (e.g., pumping station locations and pumping rates) and because they require three-dimensional simulations that are as yet too expensive to apply globally.
We tested this by implementing groundwater extraction into a regional 3D groundwater model (see Text S4 and Figures S46-S53 in Supporting Information S1) and showed almost negligible impacts on estimated decline in IFGV in this specific region.However, this doesn't mean that groundwater extraction wouldn't have large effects in other regions with heavy reliance on groundwater as primary source of fresh water (e.g., the Nile Delta (Mabrouk et al., 2018) or the Mekong Delta (Minderhoud et al., 2019)).
Lastly, sea level rise implemented in the 2D groundwater models is based on GMSL rise estimates (Pörtner et al., 2019), with values chosen for our future sea level rise scenarios corresponding to the more recent regional sea level rise indicative central estimates (Fox-Kemper et al., 2023).Regional differences in sea level rise probably wouldn't have a large impact on our groundwater salinity estimates since our 2D groundwater models represent an approximation of sub-regional conditions and elevation (see Chapter 2).However, in future large-scale analyses using 3D groundwater models such regional changes in sea level rise should be taken into account since in such cases local differences in elevation can be captured.

Groundwater Model Uncertainty
It is important to mention that there is a rather large degree of uncertainty in our 2D groundwater model results, due to both intrinsic uncertainty of the input data sets (e.g., estimated groundwater recharge, hydraulic conductivity values etc.) and simplified geological conditions.This is to be expected in a set of regional scale models based on global data sets, but it is necessary to stress that the 2D groundwater models are unfit to assess groundwater salinity conditions at a specific location.However, we believe that our analysis can provide a glimpse into future effects of sea level rise on coastal groundwater resources by showing the relative change in IFGV at the regional scale.This, because the larger differences in hydrogeological context and topography between regions can be picked up by our models, despite the large remaining uncertainty.Furthermore, the large number of 2D groundwater model runs performed in this study can be replicated in future to study other future impacts on fresh coastal groundwater reserves.When location specific information is available, they can be improved locally provide information on local-scale impacts.The technical implementation of these 2D groundwater models can be accessed on GitHub.
Additionally, rising sea levels will also lead to increasing coastal flooding hazards (Muis et al., 2016;Vousdoukas et al., 2016;Wahl, 2017) which we do not take into account.Coastal flooding due to seawater transgressions would lead to increased frequency of short groundwater salinization events.These events could however have large effects on estimated groundwater salinity conditions.Therefore, our estimates likely underestimate the decline in IFGV by omitting these coastal flooding effects.
In this sense, our study can be interpreted as a "best case scenario" in relation to future fresh groundwater decline in coastal zones.In reality, this decline would most probably be further exacerbated by groundwater pumping, seawater flooding (e.g., by storm surges, tsunamis) and declines in groundwater recharge induced by both human activity (e.g., sealing of coastal aquifers due to urbanization) or natural causes (e.g., droughts, increased aridity, or higher runoff fractions due to higher intensity rainfall).Furthermore, sea level rise will also influence the magnitude of tides which can further exacerbate the decline in IFGV in areas where tides would get higher (Pickering et al., 2017).
Even though the 3D groundwater model shows only limited effects of groundwater pumping on decline in IFGV, it can be assumed that groundwater pumping over long periods of time (i.e., decades) could have a large negative impact on the estimated fresh groundwater volume in coastal zones with high groundwater consumption (Ferguson & Gleeson, 2012;Pauw et al., 2015;Van Camp et al., 2014).Therefore, future studies should involve three-dimensional groundwater models and include groundwater pumping effects to provide more accurate current and future groundwater salinization estimates.
DEM data sets (and their vertical accuracy) play a crucial role when estimating groundwater flow and salinization patterns, as previously shown in a study in the Mekong Delta (Minderhoud et al., 2019).Therefore, we evaluate the impact of three different global DEM data sets on coastal groundwater salinization caused by sea level rise.We observe limited yet significant differences between the 2D groundwater model outcomes for the three DEM data sets.The recently developed Coastal DEM data set (Kulp & Strauss, 2019) is focused on low lying areas and enhances the vertical correction of original SRTM input (Rodriguez et al., 2006) in urban areas (defined as an area with more than 100,000 inhabitants and "urban" land use type).The 2D groundwater models using the Coastal DEM as input show higher number of SRMs in urban areas affected, compared to 2D groundwater models using the other two DEMs (Figure S44 in Supporting Information S1).These differences confirm that DEM data sets play an important role in coastal groundwater modeling and sea level rise impacts analyses.

Sea Level Rise Effects on Future Coastal Fresh Groundwater Resources
We can compare our projection of total number of people living in the coastal areas that will be affected by IFGV decline to similar projections in global coastal flooding analyses (Kirezci et al., 2020;Kulp & Strauss, 2019).These studies project that between 176 and 287 million people (Kirezci et al., 2020) to 630 million people (Kulp & Strauss, 2019) will be living in areas affected by coastal flooding by year 2100 under RCP 8.5.Our projections suggest that approximately 225 (150-453) million people that live in close proximity to current coastline (up to 10 km) will be threatened by a decline in IFGV (albeit lower than 5%) due to sea level rise by 2100, with 60 (19-96) million experiencing a decline larger than 5%, which is similar in magnitude as estimated for coastal flooding.The decrease in IFGV will negatively impact the freshwater availability in affected areas and could lead to freshwater shortages for domestic, industrial and agricultural use.Increased groundwater salinity and flooding by rising sea levels can also lead to soil salinization further increasing the stress on agricultural production in affected areas (Daliakopoulos et al., 2016;Herbert et al., 2015;Pitman & Läuchli, 2002;Qadir et al., 2014) and potentially also increasing groundwater extractions for agricultural and domestic purposes.
The economic losses by coastal floods could potentially be devastating for human communities living in coastal areas; it is estimated that coastal flooding could expose assets totaling between 8,813 and 14,178 billion USD in GDP by the end of the 21st century (Kirezci et al., 2020).In comparison, we estimate that about 4,625 (3,408) billion USD in GDP worth of assets would be threatened by IFGV decline in year 2100 under RCP 8.5.This alone is an alarming number and if we take into account the combined effects of both sea level rise and coastal flooding, the costs (human, environmental and economic) for coastal areas in the coming century could be immense.Avoiding the path of the RCP 8.5 scenario should be prioritized as a global challenge for humanity in this century, along with adaptation and mitigation plans for limiting the impacts of sea level rise (and associated effects) on fresh groundwater volumes and thus also securing freshwater availability in coastal areas worldwide.
It is important to mention that our study only focuses on coastal regions with unconsolidated sedimentary aquifer (and aquitard) systems.Therefore, it is possible that low-lying areas with sedimentary rock aquifer systems would be similarly affected by sea level rise and experience decline in IFGV.In this sense, the impacts of sea level rise on IFGV presented in this study can be interpreted as slight underestimation of the likely decline in IFGV due to sea level rise.Finally, as in many other global hazard assessments (Tiggeloven et al., 2020), we use people-and GDP-affected as a measure of possible damage.An important next step to estimate actual damages and arrive at a risk-based approach, is the assessment of damage curves; this entails an assessment of the relationship between the depth of the fresh-saline groundwater interface and the costs of not pumping groundwater, the costs of pumping too saline groundwater or exposure to saline groundwater.These coasts would be in terms of reduced productivity linked to salt damage to crops (FAO, 2021), but also in terms of the costs to below ground infrastructure etc.In addition, increasing groundwater salinities will also, affect freshwater ecosystem and cause health issues related to increased salinity in drinking water from groundwater source (cardiovascular disease due to raised blood pressure (He & MacGregor, 2009)).Making such relationships and collecting the exposure data to use these is however a formidable task.
The comparison of inland saline groundwater extent estimated by the 2D groundwater models and collected from local studies shows that our estimates match the local information reasonably well (i.e., same order of magnitude) in most cases.However, in several cases we can observe large errors, mostly in the deltaic regions (see Figure 7).This can be due to the two-dimensional modeling approach considered in our study as well as the regional scale of the 2D groundwater models compared to local scale (and measurements) in the collected local studies.Especially in the delta-areas, a two-dimensional modeling approach may underestimate saline groundwater occurrence for example, due to the neglect of salinity in the surface water system in deltas that is influenced by tidal excursions.Furthermore, we did not consider human influence such as groundwater extraction (and the associated land subsidence) on current groundwater conditions in the coastal zones which can play a large role in current saline groundwater extent.These effects can only be included using full three-dimensional groundwater models.Regardless, we consider the match with local studies satisfactory especially since our study focuses on relative change in decline in IFGV due to sea level rise rather than on absolute fresh groundwater volume quantification.

Conclusion
To summarize, this study shows that sea level rise can have substantial negative impacts on fresh groundwater resources in many coastal regions.This impact is simulated as decline in fresh groundwater volume as compared to the situation in year 2000.The approach presented in this study consists of 2D groundwater models that are built to represent regional conditions (e.g., in topography, geology) and as such are unfit to quantify actual local-scale volumes of groundwater.Thus, the results are presented as relative changes (regional sensitivities) rather than quantified as volumes.We estimate that hundreds of millions of people living in coastal regions worldwide will experience decline in fresh groundwater by year 2100, along large economical losses in order of billions USD.The severity of these negative impacts differs largely between RCP 2.6 and RCP 8.5 scenarios.
Although the 2D groundwater models are meant to represent regional-scale variability a large degree of uncertainty remains, stemming both from the use of global input data sets and from the methodology applied (using representative regional models in a probabilistic framework).As such, they only provide a limited large-scale perspective into the effects of sea level rise on future fresh groundwater resources in coastal zones.To gain new and more detailed insights a shift toward local data and 3D groundwater models needs to be made in future.

Figure 1 .
Figure1.Schematization of groundwater salinization in a simplified cross-section of a coastal region with conceptual geological conditions where aquifer and very low permeable aquitard layers (i.e., clay layers) are continuous and stacked on top of each other.(a) In a regional setting we recognize two sources of groundwater recharge-through surficial flow via infiltration and through potential flow from the mountains areas located inland.Sea water intrusion is often present due to human influences (e.g., pumping) or natural causes such as storm surges, we do not take these into account in this study.To quantify groundwater salinization due to future sea level rise we compare inland fresh groundwater volume change in the coastal zone stretching 10 km landwards.(b) Future sea level rise will presumably lead to increased salinization due to overtopping and seawater infiltration (1) or through increased pressure from the seaward direction (2).

Figure 2 .
Figure 2. Example of a coastal COSCAT region and subdivision into smaller sub-regions.(a) The coastal COSCAT region is located in equatorial West Africa., around the Niger river delta.(b) It covers two large deltaic regions -Niger river delta in the east and Volta river delta in the west.Three different coastal profile types are distinguished-Simple (S), Island (I) and Delta (D)-based on their topographical profile characteristics or location within a deltaic area.(c) Individual coastal profiles are clustered together based on their proximity and coastal type resulting in this case into 21 sub-regions.A sub-regional representative model (SRM) (SRM) will then be built for each of these 21 sub-regions.This procedure is implemented for all selected coastal COSCAT regions.

Figure 3 .
Figure 3.An example of estimated groundwater salinization (the coastal segment of cross-sectional profile sub-regional representative model 1415_32) due to sea level rise under three different Representative Concentration Pathway scenarios for seven time steps.The blue rectangle delineates the boundary of the inland fresh groundwater volume (IFGV) assessment area, while the numbers in the lower left corner of each coastal segment show the percentage of fresh groundwater in the IFGV assessment area.

Figure 4 .
Figure 4. Difference in inland fresh groundwater volume (IFGV) in year 2100 expressed as percentage IFGV compared to situation in year 2000.Results are averaged over the three different digital elevation model inputs used in our modeling study for each Representative Concentration Pathway (RCP) scenario-(a) RCP 2.5, (b) RCP 4.5 and (c) RCP 8.5.

Figure 5 .
Figure 5. Increased inland salinization extent estimated in year 2100 expressed as inland shift (km) of saline water wedge (SWW, 17.5 g/l TDS considered as boundary-half of ocean water salinity) compared to situation in year 2000.Results are averaged over the three different digital elevation model inputs used in our modeling study for each Representative Concentration Pathway (RCP) scenario-(a) RCP 2.5, (b) RCP 4.5 and (c) RCP 8.5.

Figure 6 .
Figure 6.Proportional schematization of affected sub-regional representative model areas per Representative Concentration Pathway scenario quantified by decline in inland fresh groundwater volume compared to situation in year 2000.

Figure 7 .
Figure 7.Comparison of observed and estimated groundwater intrusion inland extent (TableS2in Supporting Information S1 lists all the collected studies).The marker size represents the frequency of discrepancies between observed and estimated values as some studies show the same extent of salinity inland while our estimates in these sub-regional representative models can also show identical values.

Table 1
Total Number of People Living in the Sub-Regional Representative Models (Up to 10 km From Current Coastline) Summarized for Each Representative Concentration Pathway Scenario and Time Step