Modeling saltwater intrusion scenarios for a coastal aquifer at the German North Sea

A 3d regional density-driven flow model of a heterogeneous aquifer system at the German North Sea Coast is set up within the joint project NAWAK (“Development of sustainable adaption strategies for the water supply and distribution infrastructure on condition of climatic and demographic change”). The development of the freshwater-saltwater interface is simulated for three climate and demographic scenarios. Groundwater flow simulations are performed with the finite volume code d3f++ (distributed density driven flow) that has been developed with a view to the modelling of large, complex, strongly density-influenced aquifer systems over long time periods. INTRODUCTION The project NAWAK started in 2013 as an associated project of eight German institutions including research institutes as well as local water supplying companies. It was funded by the German Federal Ministry of Education and Research (033W007A). The main objective was the implementation of a planning tool that enables the local stakeholders to know the possible range of impacts caused by climate and demographic change and to identify sustainable adaption strategies. In the BMBF-funded project “Development of sustainable adaptation strategies of water management in coastal areas under the conditions of climatic and demographic change” (NAWAK, was applied to coastal aquifer systems near the German North Sea. Three regional 3d density-driven flow models have been set up using the code d3f++. The aim of this works was forecasting the impact of different climatic and demographic scenarios on the freshwater availability. The largest of these models will be presented here, the Sandelermöns region, covering an area of about 1000 km2 and based on detailed geological and geochemical data, a subarea related recharge model (PANTA RHEI, Meon et al. 2012) as well as data from geophysics. The model contains groundwater recharge, river discharge and the well fields of three waterworks. Based on the calibrated model three scenarios regarding the impact of sea-level rise and changing rechargeand well pumping conditions were simulated. The results were provided to the local and regional stakeholders within a special developed planning tool. Furthermore, the same modeling approach was applied to two other coastal areas near the estuary of the river Elbe. © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/). E3S Web of Conferences 54, 00031 (2018) https://doi.org/10.1051/e3sconf/20185400031 SWIM 2018


INTRODUCTION
The project NAWAK started in 2013 as an associated project of eight German institutions including research institutes as well as local water supplying companies. It was funded by the German Federal Ministry of Education and Research (033W007A). The main objective was the implementation of a planning tool that enables the local stakeholders to know the possible range of impacts caused by climate and demographic change and to identify sustainable adaption strategies.
In the BMBF-funded project "Development of sustainable adaptation strategies of water management in coastal areas under the conditions of climatic and demographic change" (NAWAK, was applied to coastal aquifer systems near the German North Sea. Three regional 3d density-driven flow models have been set up using the code d³f++. The aim of this works was forecasting the impact of different climatic and demographic scenarios on the freshwater availability. The largest of these models will be presented here, the Sandelermöns region, covering an area of about 1000 km² and based on detailed geological and geochemical data, a subarea related recharge model (PANTA RHEI, Meon et al. 2012) as well as data from geophysics. The model contains groundwater recharge, river discharge and the well fields of three waterworks.
Based on the calibrated model three scenarios regarding the impact of sea-level rise and changing recharge-and well pumping conditions were simulated. The results were provided to the local and regional stakeholders within a special developed planning tool. Furthermore, the same modeling approach was applied to two other coastal areas near the estuary of the river Elbe.

HYDROGEOLOGICAL MODELING
The Sandelermöns model covers an area of about 1,000 km² situated at the German North sea coast as shown in Figure 1. Almost half of the model region is low-lying and formed part of the North Sea some hundred years ago so that the aquifers lead saline groundwaters. The south-western part consists of moraines with notable recharge rates and freshwater aquifers below. The three waterworks Sandelermöns, Feldhausen and Klein Horsten are situated in these areas, competing for pumping concessions and worrying about the close-by saline waters. The hydrogeological structure model was set up by the Oldenburg-Ostfriesian Water Company (OOWV). Six formations are distinguished, starting with a thick layer of fluvial sands at the bottom, followed by small fields of Tergast clay and a continuous layer of meltwater sands. The top of the model consists of thin layers of Lauenburg clay, dune sands and silty materials. The base surface of each hydrogeological layer was read-in into the ProMesh tool (Reiter 2016) and together with ground surface data composed to a 3d model.
The south-western boundary of the model is located on a watershed and therefore assumed to be impermeable. The salt concentration is set to zero. The north-western and south-eastern boundaries are choosen to be perpendicular to the water table isohypses and therefore also regarded as impermeable for the flow. The coastal boundary is equipped with a hydrostatic pressure for seawater and a salt concentration of 34 kg/kg (35 g/l). The bottom of the model is also assumed to be impermeable because geological knowledge suggests almost impermeable clayey formations here. On the upper boundary, a Dirichlet boundary condition p = 0 has to be choosen for the pressure prescribed by the levelset method. The concentration is also set to zero.
The groundwater recharge was simulated by using the hydrologic model PANTA RHEI (LWI), a conceptual model with partly physic-based modules, especially for the soil water processes. 147 polygonal recharge zones were distinguished (see figure 2). Time-dependent recharge rates for the three climate scenarios WETTREG2010 B1 (moderate), WETTREG2010 A1B (long-term dry) and WETTREG2010 A2 (dry) from 2011 to 2050 were computed and coupled with the density-driven flow model d³f++. A sea-level rise of 0.25 m was assumed for this period.
The north-eastern region of the model domain is characterized by a dense net of small draining ditches and rivers conducting water to the coastal pumping stations in the dikes to keep the groundwater level below land surface. This drainage plays a crucial role in the hydraulic regime and had to be incorporated in the model. Because an explicit mapping of all ditches in a regional model is impossible, only the rivers of first and second orders were integrated, and their influence was smeared over a user defined range of surface elements.
Regrettably there exists only little information about the water levels of the receiving streams in the Sandelermöns region as well as the pumping rates of the coastal stations. Missing information was replaced by reasonable assumptions to meet the natural hydraulic regime by a careful calibration process. Furthermore, the 51 pumping wells of the waterworks were included into the model. Additionally, 36 private wells are regarded. For 2050 three different socio-economic scenarios influencing the water demand were investigated by the OOWV, a basic scenario, a green world scenario and a growth scenario. The resulting time-dependent pumping rates were used in the scenario-simulations by d³f++. To find an appropriate initial condition for the free groundwater table the data of 284 gauge wells were averaged over the year 2011 and interpolated. Unfortunately, only few gauges exist near the coastal line, and the model reacts very sensitive on changes in these settings. The initial condition for the salt concentration is based on geoelectromagnetic data provided by the Leibniz Institute for Applied Geology.

METHODS
For the density-driven flow model the finite volume code d³f++ is used (Fein 1999, Schneider 2016). The code may be used for groundwater flow and transport modelling in porous media as well as in fractured rock or mudstone including options for modelling saltand heat transport as well as a free groundwater surface.
In d³f++, the equation system describing thermohaline flow is solved, consisting of the mass conservation of the fluid, the mass conservation of the brine, and energy conservation. The flow velocity follows Darcy's law. Fluid density and viscosity are depending on salt mass fraction and temperature. It should be mentioned that d³f++ solves the complete, nonlinear coupled equation system without simplifications such as the Boussinesq approximation.
d³f++ is based on the UG toolbox (Vogel et al. 2014). The use of modern numerical methods such as geometric and algebraic multigrid methods and their parallelisation enables simulations over long time periods with feasible computational effort.
In d³f++ the time dependent position of the groundwater table ( ) in a fixed domain is described as the zero set of a level set function (Frolkovič 2012). Using level set functions implicates some restrictions to these part of the model boundary: Regarding the pressure, the boundary condition on ( ) at any fixed time is set directly by numerical discretization and cannot be changed. That means groundwater recharge as well as discharge may not be treated as boundary conditions, both effects have to be modelled as factors directly influencing the normal velocity on ( ).

RESULTS
The main objective of this works was predicting the situation of the freshwater-/saltwater interface up to the year 2050. Simulations started with the calibrated model in 2011 and the initial and boundary conditions as described above.
Firstly, the groundwater model had to be calibrated, based on mean values of the recharge and pumping rates in the base year of 2011. It was observed that river drainage relatively to the recharge rates is much more influencing the results than well pumping, which plays a less important role. These sensitivities form a huge problem because of the low number of gauges und the poor knowledge about the pumping rates at the overflows of the dikes. Another problem is that the results for the fluid volume may completely fail without an adequate grid refinement. Subject to the uncertainties caused by lack of data, a model calibration could be reached.
Based on the calibrated model 9 scenarios were simulated, combining both climatic and socio-economic scenarios. One difficulty was that the current distribution of saline waters con not be explained only by the hydraulic situation. The influence of the hydrostatic pressure boundary condition at the coastal line is restricted to a zone near the coastal boundary. Saline waters detected by measures outside this zone must result from former inundations hundreds of years ago, before the areas were protected by dikes. That means the simulated salt distribution is mainly determined of the initial condition. Using the airborne electromagnetic data from LIAG was therefore essential for the simulation of the scenarios.
The results of the combination of the WETTREG B1 climate scenario and the socioeconomic scenario basic are shown in figure 3 and figure 4. Figure 4 shows the movement of the freshwater-saltwater interface within 39 years of model time in a depth of 65 m, corresponding to the mean screen of the well field next to the saline areas. distribution and the influence of the drainage ditches. Since the initial salt distribution cannot be derived from the hydraulic regime alone, a good data base for both of these factors is essential to obtain confidential results. To be able to quantify the exfiltration by the drainage ditches in the model it is important to have more gauches and, at the other hand, that the pumping stations measure the run-off to the sea at the overflow of dikes.
Furthermore it was detected that a minimum grid refinement is crucial for getting correct results for the fluid volume or the position of the groundwater surface, respectively. Thereby it turned out again that grid convergence is an indispensable precondition for confidential simulation results, especially regarding forecasts.
In consideration of the lack of data and the still coarse grid the calibration results are rather good and plausible. Nevertheless, further grid refinements are necessary for reliable simulation results. It is desirable to achieve grid convergence, that means at least grid levels 3 or 4 (1.1 or 4.6 millions of nodes) have to be used. In the framework of the project "Implementation of strategic development goals in coastal aquifer management" (go-CAM, 02WGR1427B) the described work is continued.