Modelling climate change effects on a Dutch coastal groundwater system using airborne electromagnetic measurements

The forecast of climate change effects on the groundwater system in coastal areas is of key importance for policy makers. The Dutch water system has been deeply studied because of its complex system of low-lying areas, dunes, land won to the sea and dikes, but nowadays large efforts are still being done to find out the best techniques to describe complex fresh-brackish-saline groundwater dynamic systems. In this paper, we describe a methodology consisting of high-resolution airborne electromagnetic (EM) measurements used in a 3-D variable-density transient groundwater model for a coastal area in the Netherlands. We used the airborne EM measurements in combination with boreholelogging data, electrical conductivity cone penetration tests and groundwater samples to create a 3-D fresh-brackishsaline groundwater distribution of the study area. The EM measurements proved to be an improvement compared to older techniques and provided quality input for the model. With the help of the built 3-D variable-density groundwater model, we removed the remaining inaccuracies of the 3-D chloride field and predicted the effects of three climate scenarios on the groundwater and surface water system. Results showed significant changes in the groundwater system, and gave direction for future water policy. Future research should provide more insight in the improvement of data collection for fresh-brackish-saline groundwater systems as it is of high importance to further improve the quality of the model.


Introduction
Fresh water is an essential resource for all types of human activities.Worldwide, agriculture is the sector consuming about 70 % of all available fresh water (UNESCOPRESS, 14 March 2012).This makes agriculture a very vulnerable sector to changes in the water system.Densely populated areas and intensive agricultural activities in the coastal zone demand large quantities of fresh water to maintain their economy.However, fresh water in this zone is not abundant.In addition, it is threatened by salinisation processes such as lateral salt water intrusion as well as upward saline seepage (Custodio and Bruggeman, 1987;Jelgersma et al., 1993;Oude Essink, 1996;FAO, 1997;Post and Abarca, 2010;Oude Essink et al., 2010).The anticipated sea level rise and associated changes in recharge and evapotranspiration patterns will also intensify the pressure on the coastal groundwater system.
Salt water intrusion into coastal aquifers, on one hand, is a well-studied process worldwide which has been topic of many research works in the last years (Custodio, 2010;Barlow and Reichard, 2010;SWIM proceedings).Saline seepage, on the other hand, is so far only studied in the Netherlands (Oude Essink, 1996, 1999, 2001a; de Louw et al., 2010de Louw et al., , 2011)).Saline seepage is a process which takes place in lowlying areas with piezometric heads lower than its surroundings (e.g.mean sea level).In low-lying areas with a saline to brackish subsurface, seepage is a serious threat to fresh water resources, leading to a reduction of the fresh groundwater volume in the shallow subsurface.In coastal areas these shallow fresh groundwater volumes, known as rainwater lenses, are an important fresh groundwater resource for the agriculture.For this reason, nowadays different research studies focus on the dynamics of these shallow rainwater lenses (e.g. de Louw et al., 2011).
To ensure agricultural sustainability and fresh water supply of a coastal area, both policy makers and farmers need to know the future availability of fresh groundwater and the dynamics of saline groundwater in the deep and shallow subsurface.To implement appropriate countermeasures to protect coastal fresh water resources from salt water intrusion, we need to better understand the relevant salinisation processes.The possibility to more accurately predict the dynamic distribution of fresh, brackish and saline groundwater has opened the doors to hydro(geo)logists and policy makers to design strategic countermeasures (Klein et al., 1998;Oude Essink, 2001b).The challenge is to produce reliable quantitative hydrogeological information, good enough for these purposes.
In this paper we describe a model instrument which is able to assess the effects of the climate change on a coastal zone of the northwestern part of the Netherlands.Seasonal dynamics on the groundwater system and seasonal effects on a regional scale of three different climate scenarios are predicted.A 3-D variable-density transient groundwater model is built using existing data from the previous model MIPWA (Snepvangers et al., 2007) and new data acquired with airborne electromagnetic (AEM) surveys.AEM has already been successfully used by e.g.Sengpiel and Meiser (1981) or Fitterman and Deszcz-Pan (1998) in fresh-brackish-saline groundwater environments since the 1980s.However, it is within this framework of the Interreg IV-B project CLIWAT (www.cliwat.eu, a transnational project in the North Sea region) that Deltares, TNO, Province of Fryslân, Wetterskip of Fryslân (all the Netherlands), Federal Institute for Geosciences and Natural Resources (BGR, Germany), and Aarhus Geophysics (Denmark) worked together to make these AEM methods suitable and accessible for mapping fresh groundwater resources over large areas, and to use this AEM data in numerical modelling tools for the prediction of climate change effects on groundwater systems in northern Fryslân.In addition to this study there are other initiatives that combined AEM data and groundwater modelling (Sulzbacher et al., 2012).We built a new geological model and a new 3-D fresh-brackish-saline groundwater distribution with the data obtained through different geophysical AEM campaigns using frequency domain (Resolve) and time domain (SkyTEM) helicopter-borne systems.The final result is a calibrated groundwater model containing detailed underground data.The changes to the groundwater system caused by different climate scenarios are shown.Furthermore, the outcome of the model was processed with the objective that policy makers could directly employ it in the making of water management plans.

Study area
The study area (30 × 30 km 2 ) which is located in the northwestern part of the Netherlands in the province of Fryslân (Fig. 1), borders the shallow Wadden Sea to the north-west.The landscape consists of polders and other anthropogenic modified landscape elements such as dikes and dwelling hills.The main city and capital of the province is Leeuwarden, home of about 100 000 inhabitants.The study area has  an average surface elevation at mean sea level (m.s.l.) having just some dikes and cities outstanding at a higher altitude of maximum +11 m m.s.l. and some low nature areas lying at −2.1 m m.s.l.The area is drained by an intricate network of ditches to make the land fit for agriculture, an important component of the Friesian economy.The excess of water is pumped out into a higher elevated water reservoir and transported to the sea.The surface water levels are kept at a constant summer and winter level.The artificially controlled areas below sea level are called polders (van der Ven, 1993).

Geology
In the study area fluvial sediments, mainly consisting of sand and occasionally clay, were deposited in the early and middle Pleistocene (Peize-Waalre -PZWA, Appelscha -APand Urk -UR -formations).In the Elsterian ice age, tunnelvalleys formed and were filled with sand, sometimes capped by clay (Peelo formation -PE).After the Elsterian, fluvial deposits were again dominant (Urk-Tynje -URTY -formation) until the land-ice reached the area in the Saalian iceage and a till layer at about 10-30 m depth was deposited (Drente formation -DR).The till covers the entire area and dips slightly to the north-west.Its thickness is erratic and can vary considerably over short distances.The till is considered to be an important layer in the groundwater flow, due to its low hydraulic conductivity.After the land-ice retreated, a marine transgression caused the sedimentation of sand and clay (Eem formation -EE).Another ice-age (the Weichselian) did not produce land-ice in the area but widespread aeolian cover sands: the Boxtel formation -BX.The latest (Holocene -HL) transgression produced a sequence of alternating sand and clay deposits (with locally some peat), ranging from more than 10 m thickness in the northwest to less than 2 m in the southeast.Figure 2 shows the general geological setting in a cross-section.Due to the geological history of the area, the hydrogeological structure of aquifers and aquitards can be clearly differentiated.From bottom to top we find the hydrogeological base located at 270 m below mean sea level.Overlying this base we find the deepest aquifer which has a thickness of about 200 m and a transmissivity of 800 m 2 day −1 and is formed by the PZWA, AP, UR, PE, and URTY formations.This aquifer is confined by a layer 10 to 30 m thick (DR formation).Right on top we find the shallowest aquifer (BX formation) which is thin compared to the deepest aquifer and has a transmissivity between 30 and 300 m 2 day −1 .Overlying the shallowest aquifer, we find the Holocene layer (HL) which partly works as a confining layer.
In large parts of the study area, the groundwater is saline with the boundary between brackish and saline lying at less than 5 m below surface level.Only in the south-eastern part of the area the groundwater is fresh.The rich geological history and the saline groundwater offer an ideal field to test the application of airborne EM within variable-density groundwater models.Moreover, changes in climate and sea level are expected to have relevant effects on this coastal groundwater system and may have influence on the sustainability of the present water resources management.

Drillings
The study area is covered by more than six thousand drillings, 95 % of them only reaching a depth of 5-10 m, covering the Holocene and the upper part of the Pleistocene deposits.We extracted all the drillings in the area from the DINO database, the national database for subsurface data in the Netherlands (www.dinoloket.nl), in order to model the Holocene sequence in terms of lithological content.

Geophysics
Several geophysical techniques were used, both airborne and ground-based, to map the electrical conductivity of the area providing both information on lithology and groundwater.In the following paragraphs, we give a brief summary of the geophysics, but for a more detailed description of the techniques see Gunnink et al. (2012).

Ground-based geophysics
We carried out 71 ECPTs (electrical conductivity cone penetration tests) in order to measure the cone resistance and sleeve friction (from which the lithology is determined), as well as the electrical conductivity (EC) of the sediment and the groundwater together, which we call EC(total).In Fig. 3 the location of the ECPTs is depicted.

Airborne EM
We used airborne EM to map the electrical conductivity of the subsurface by means of two helicopter surveys.A frequency domain EM survey (HEM) was conducted by the German Federal Institute for Geosciences and Natural Resources (BGR), and a time domain EM survey (SkyTEM) was carried out by the University of Aarhus, Denmark.The surveys showed a small areal overlap (see Fig. 3).A radar station for air traffic control strongly disturbed the airborne HEM measurements; hence, an area of about 7 by 10 km in the centre of the planned survey could not be covered.The data of both airborne surveys acquired in summer 2009 were processed and inverted to smooth multi-layer electrical conductivity models.Siemon et al. (2009) presented a general description of both AEM methods.Details of the data collection and processing are given in Gunnink et al. (2012).Figure 4 shows two depth-slices and two cross sections which were interpolated from the AEM models.The slices show the spatial distribution of the EC, which formed a valuable source for establishing the salinity of the groundwater.The AEM profiles of EC were also used to determine the presence and thickness of the till in the area.We employed an artificial neural network (ANN) to use the distinct pattern of EC in the till (which is less conductive compared to the sediments above and below the till) to map the till.Results of the ANN,    in terms of probability that a certain voxel (volume element in a 3-D grid) contains till, are given in Fig. 5.For more details about the ANN method, see Gunnink et al. (2012).

Salinity data
In the study area, only limited information is available about the chloride concentration in groundwater.In general, large parts of the area have saline groundwater conditions with a chloride concentration higher than 1000 mg l −1 .To be able to convert the EC (total), -sediment + groundwater -to EC(groundwater), the influence of the lithology on the EC (total) needs to be determined by establishing the formation.
We measured the EC (total) and at the same location we collected groundwater samples from which we measured EC (groundwater).From these two sources we could obtain the formation factor (Sect. 2.3.4).
To obtain the current salinity of the groundwater, we used the following datasets:    -EC profiles from borehole-logs, converted to Cl − (1333 samples at 26 locations, measured once, from the DINO database) -EC profiles from the airborne EM, converted to Cl − (> 50 000 locations, measured once) -EC profiles from the ECPTs, converted to Cl − (71 locations, measured once).

Numerical code MOCDENS3D
The 3-D variable-density transient groundwater flow and coupled salt transport is simulated with the model code MOCDENS3D (Oude Essink, 2001a;Bakker et al., 2004;Vandenbohede, 2008;Vandenbohede et al., 2009).It is based on MODFLOW (McDonald and Harbaugh, 1988) and MOC3D (Konikow et al., 1996), which has been adapted for density differences.for groundwater flow in porous media.Advective and hydrodynamic dispersive solute transport processes through porous media are modelled by a particle tracking technique in combination with the finite difference method (Konikow et al., 1996).In Dutch coastal groundwaters, chloride is the major conservative negative ion.As such, the discussion about salinisation is focused on that predominant solute.Under the given circumstances in the Dutch coastal aquifers, the Oberbeck-Boussinesq approximation is valid as it is suggested that the density variations (due to concentration changes) remain small to moderate in comparison with the reference density (ρ) throughout the considered hydrogeologic system (Holzbecher, 1998;Nield and Bejan, 1992).As such, a substantial simplification of the governing differential equations can be derived.For additional information on the mathematical formulation, see e.g.Oude Essink et al. ( 2010).

Model discretization and boundaries
The 3-D model was built with cells of 100 by 100 m.This resolution was chosen based on the available data and a good compromise between the details of the pursued knowledge and the model runtime.Furthermore, this scale provided an optimal balance to model the regional hydrological processes and local details.We divided the area in 300 columns, 300 rows and 50 model layers.Every model layer has a constant thickness, but the thickness of the different model layers differs and increases in depth.As the model output is focused at the salt load from groundwater to the surface having a high spatial variation, we used a higher resolution for the shallow part of the model.−270 m m.s.l. and is the no flow bottom boundary of the model.The lateral boundaries are constant flow boundaries derived from a much larger existing regional and calibrated groundwater model that includes our study area (MIPWA model, Snepvangers et al., 2007).However, this larger model does not take into account the density-dependent flow.As the present model requires fresh water heads and the study area has saline groundwater, we corrected the heads to fresh water heads using the following equation (Oude Essink, 2001c;Post et al., 2007): where density of the groundwater (kg m −3 ); C = chloride concentration of the water (mg Cl − l −1 ); ρ f = reference density of the fresh water at mean ground temperature (equals 1000 kg m −3 ).
In addition, a linear equation of state couples groundwater flow and solute transport:   where β C is the volumetric concentration expansion gradient (m 3 kg −1 ), in this case 1.34 × 10 −6 l mg −1 Cl − .The concentration used for the constant flow boundaries is the one introduced as initial 3-D chloride field in the model (see Sect. 2.3.4).
The time horizon for which we wanted to analyze the changes of the groundwater system was the year 2100.Therefore, we created a model that starts in 2005 and ends in 2100.In order to reduce calculation time, we simulated the first 95 yr as a steady state with stress periods of one year with a yearly average recharge.After these 95 yr, a transient simulation was carried out in order to simulate the difference between a summer and a winter period.The transient simulation had time steps of 10 days.Test results indicated that 2-yr period of transient simulation was long enough to obtain a representative summer and winter situation.In this way, we rapidly achieved the situation of 2100 while still being able to simulate the winter and summer dynamics.

Hydraulic conductivity field generation
Hydraulic conductivity is an important parameter in hydrological modelling and is closely linked to the geological history and the sediments that have been deposited.To generate a hydraulic conductivity field for the study area, we started by constructing a geological/lithological model.We paid special attention to the Holocene sediments, because these constitute the upper part of the model in which upward seepage takes place and canals and ditches are located.
The lithological model was constructed on a 3-D grid of 100 × 100 × 0.5 m, a so-called "voxel" grid.For every voxel, the lithology was estimated from surrounding drillings and geological knowledge about the geological history.The hydrological model uses grid cells that have varying thickness with depth, and therefore the hydraulic conductivity from the original 100 × 100 × 0.5 m needs to be upscaled to the size of the voxel of the hydrological model.

Holocene deposits
The Holocene deposits in the area consist of alternating clay and sand with occasionally peat.The depositional circumstances resembled a "Wadden-type" environment, in which sand and clay occur in a fast alternating sequence.This knowledge is derived from the analysis of the large amount of shallow drillings in the area, which allowed us to use a stochastic simulation technique to model the lithologies in the Holocene deposits.We used sequential indicator dimulation to obtain stochastic realizations of lithological contentpeat, clay, sandy clay, fine sand and medium to coarse sand -which were used later to obtain a spatial distribution of the hydraulic conductivity.For details about the simulation technique, see Stafleu et al. (2011).We calculated 50 equiprobable realizations of the lithology.In Fig. 8 an example is given of one of the realizations.
The horizontal and vertical conductivity values for the Holocene deposits were assigned to the 5 different lithology classes we distinguished.In Table 1 these values are given.They are based on data that were collected in the past from a wide range of sources: laboratory tests, field tests and experience from users.The data collection for the hydraulic conductivity was carried out for the entire Netherlands and was summarized in the values given in Table 1.As there was never a special interest in the aquifer properties for drinking water because the groundwater is saline in this area, there are no pumping tests available in the study area.
For each one of the 50 realizations of the lithology, the natural logarithm (Ln) of the hydraulic conductivity for that lithology class was assigned to each voxel.We took the average of the 50 realizations of the hydraulic conductivity to get a representative Ln-hydraulic conductivity for every voxel in     the Holocene deposits.After taking the exponent of the Ln of the average hydraulic conductivity, we obtained the horizontal and vertical hydraulic conductivity for every voxel in the Holocene sequence.

Glacial till
The glacial till is an important layer in the groundwater model.Due to the clay content in the till and the fact that it was mechanically consolidated by the ice sheet it was deposited from, it acts as an impermeable layer.The till was mapped using the airborne EM models that were converted into the probability of having till using an artificial neural network (ANN) technique (Gunnink et al., 2012).The ANN procedure delivered a probability of glacial till for every meter in the vertical direction.Going from surface-level downwards, the first interval with a probability > 0.8 was chosen as the top of the till.As soon as the probability became less than 0.8, we assumed that the bottom of the till was reached.
In this way, every AEM model was converted into a top and bottom of the till.Together with the drillings and the current regional till model, the top and bottom of the till were calculated in the study area.

Pleistocene units
The Pleistocene units were taken from the REGIS II.1 model, which is a regional 2-D grid model containing top and base of aquifers and aquitards and its hydraulic conductivity.The glacial till in the current REGIS model was replaced by the till as it was modelled for this study (see above).
The 2-D grid model of top and base of the geohydrological units from REGIS and their corresponding hydraulic conductivity were converted to a 3-D voxel model of 100 × 100 × 0.5 m.For every voxel we determined the geohydrological unit to which the voxel belongs, based on the top and base of the unit.The voxel was then assigned the hydraulic conductivity corresponding to that unit.The 3-D voxel model was completed by adding the Holocene units on top of the Pleistocene units.
This 3-D voxel model consists of more than 200 layers, causing extremely large calculation times for density-driven groundwater flow.Therefore, in the 3-D variable-density groundwater model, the thickness of the model layers was increased with depth.For every cell in the groundwater numerical model, we calculated the hydraulic conductivity (k) based on the average of the Ln(k) values that corresponded to the initial voxel model.In this way, the hydraulic conductivity was upscaled correctly, since averaging the logarithmic values of hydraulic conductivity amounts to the geometric mean.The exponent of the average Ln(k) was calculated to obtain the hydraulic conductivity for the modelling.

Deriving a 3-D chloride field
The EC data derived from geophysical measurements (borehole-logging, AEM, ECPTs) were interpolated to the 3-D voxel grid using a kriging interpolation routine.This resulted in a 3-D model of EC(total) for every voxel of Hydrol.Earth Syst.Sci., 16, 4499-4516, 2012 www.hydrol-earth-syst-sci.net/16/4499/2012/ 100 × 100 × 0.5 m.To convert the EC (total) to EC (groundwater), we needed to know the intrinsic formation factor (FF i ).The FF i (Archie, 1942) is calculated as the ratio between the total resistivity (or EC) and the resistivity of the groundwater and is valid for clay-free sediments.Sulzbacher et al. (2012) used EC (groundwater) measured in boreholes and EC(total) from AEM surveys to compute the formation factor for the island of Borkum.We used boreholes to measure the EC(total) and to sample groundwater and measure the EC (groundwater).We also recorded the lithology at the depth of measurements.Unfortunately, not all boreholes had EC (total) measured, and therefore the EC (total) for some of the boreholes was extracted from the 3-D model of interpolated EC.This dataset was supplemented by 6 locations, at which an ECPT was conducted (resulting in EC (total)) together with groundwater samples that were collected to measure EC (groundwater).The estimation of the lithology was carried out by classifying the ECPT data in lithology classes (sand, clay and sandy clay).The amount of data for establishing the formation factor is rather limited.Only samples in which both EC (total) and EC (groundwater) are measured and from which the lithology is known can be used for such a purpose.Only the first two data sources mentioned in Sect.2.2.3 have measured EC (total) and EC (groundwater), and for the data from the first source the EC (groundwater) was often lacking.
From this limited dataset of EC (total) and EC (groundwater) data, we could establish a relation for the sandy sediments between the total EC and the EC of the groundwater, the FF i .Since our data are not clay-free, we applied a correction following the method as described in Soupios et al. (2007) to obtain the apparent formation factor, FF a .The relation between the FF i (as measured for each sample) and the FF a is given by where 1/FF i is the intercept of the straight line with the yaxis, and BQv/FF i represents the gradient.Figure 9 shows the relation between groundwater resistivity and the reciprocal of the apparent formation factor.BQv is related to the effects of surface conduction, mainly due to the clay particles.By plotting the 1/FF a vs. the fluid resistivity, we obtained the FF i .
Figure 9 also shows a large variability in FF a .This variability could be caused by an undetected variability in the lithology (especially clay content) from the interpreted ECPT.It could also be that the EC(total) for some boreholes, which was determined from the 3-D model of interpolated EC, is not representative for the location of the borehole.Another option could be that the high EC of the groundwater obscured the signal from the sediment.The intrinsic formation factor for sandy sediments was nevertheless determined from Fig. 9.It has a value of 5.2, which is in line with values found in sandy sediments in the Netherlands ( de Louw et al., 2011).Unfortunately, there were not enough samples to calculate   an intrinsic formation factor for the non-sandy sediments.Therefore, we used the average FF i for the clay and sandy clay from de Louw et al. (2011), which is about 3. From the geological voxel model (see Sect. 2.3.3),we took the value of 0.05 m day −1 to discriminate between clay and sand.For every voxel with hydraulic conductivity < 0.05 m day −1 , we applied the FF i of 3, while for voxels with a hydraulic conductivity > 0.05 m day −1 , we applied a FF i of 5.
All EC collected data were corrected for temperature assuming that the temperature of the groundwater at measuring was 10 • C. EC is standardized at 25 • C, and therefore the EC from the field measurements was corrected to reflect the EC at 25 • C.
The groundwater samples, from which both EC groundwater and Cl − were determined, were used to construct a linear regression for converting EC groundwater to Cl − .The fit of 15 samples, which have both Cl − and EC(groundwater), was good with an R 2 of 0.97.The linear equation is Cl g l −1 = EC (water) S/m • 3.66 − 0.42 (4) which is similar to the regression found by de Louw et al. (2011).In Fig. 10 the Cl − distribution is shown for the upper 40 m.At this point the 3-D chloride concentration distribution for voxels of 100 × 100 × 0.5 m was complete.To bring this distribution to the 3-D groundwater numerical model, with increasing model layer thickness (see Sect. 2.3.2),we averaged the ln(Cl − ) from the detailed voxel model to the groundwater numerical model and exponentiated the average chloride concentration to obtain the chloride concentration for the cells of the 3-D variable-density groundwater model.In Fig. 11 the workflow for obtaining kh, kv and Cl − is summarized.
Although this 3-D chloride field was the result of the interpolation of high-resolution AEM measurements, the chloride field was still not in balance with the groundwater system.We reached this conclusion after running the model for the first time for 50 yr.Salinisation and freshening of groundwater are long term and slow processes that can take many decades unless extreme events occur that unbalance the system (such as tsunamis or flooding) (e.g.Oude Essink and Van der Linden, 2005).The Dutch system is not in a total state of dynamic equilibrium because of the sea transgressions and creation of low-lying polders of the last centuries (Schultz, 1992;Oude Essink, 1996;Post, 2003); nevertheless, no extreme events have unbalanced the system lately.Therefore, when the autonomic salinisation processes are simulated, we do not expect substantial changes in concentrations.When we first run the model with the initial 3-D chloride field, we calculated the volume of fresh to brackish groundwater with a chloride concentration less than 1.5 g l −1 for every 5 yr. Figure 12 shows some of the results of this analysis.The results show that during the first 15 yr the number of fresh-brackish water cells changes rapidly compared to the years after, indicating that the system was not yet in equilibrium with the model boundaries and stress terms.The reasons for this non-equilibrium could depend on the transformation of the AEM data to chloride concentration values.The fact that the distinction between sandy and clayey sediments is set at a vertical hydraulic conductivity of 0.05 m day −1 might filter out the small-scale variability that might have influenced the airborne resistivity measurements.Another reason could be that the variation in thickness of the model grid cells causes a loss of detail in the chloride concentration distribution.These inaccuracies show that there is more research required regarding the incorporation of airborne EM techniques in numerical modelling.

Hydrol
Based on this analysis, and in order to eliminate irregularities, we decided to use, as initial 3-D chloride concentration field, the one generated after 15 yr of simulation time.

Surface water and drainage network
The complex surface water and drainage network of Fryslân is also characteristic for the rest of the Netherlands.In our model, we added the surface water and the drainage network by means of the river and drainage MODFLOW packages.The data we entered were originals from the earlier mentioned calibrated model MIPWA.This model is based on data provided by the water board of the region about the location, morphology and summer and winter levels of the surface water, and location and height of the drains.Nevertheless, the model schematization of MIPWA follows the geological layers while our model layers have a constant thickness (Fig. 6).For that reason, we adapted all the parameters corresponding to these packages for our model schematization.

Recharge
The recharge we put into our numerical model is calculated by the NHI (Netherlands Hydrological Instrument).The basis of the NHI is a state-of-the-art coupling of the groundwater (MODFLOW), the unsaturated zone (metaSWAP) and the surface water (MOZART-DM) models.The resolution of the NHI is 250 by 250 m, and the groundwater flow is computed on a daily basis.One of the interim results of the NHI is the groundwater recharge.The recharge is calculated as precipitation minus evapotranspiration plus irrigation.
For the reference (current) situation we calculated the mean recharge for the summer (1 April-30 September) and for the winter (1 October-31 March) for the period 1996-2005 (Fig. 13).The KNMI climate scenario W+ (Fig. 13) predicts a decrease of the precipitation of 13 % and an increase of the evapotranspiration of 13 % for the summer.For the winter, an increase of the precipitation of 9 % and an increase of the evapotranspiration of 5 % are forecasted (Klein Tank and Lenderink, 2009).The recharge for this climate scenario was also calculated with the NHI.

Model calibration
The model calibration was done by hand and consisted of adjusting the hydraulic properties of the glacial till until a satisfactory fit between calculated and observed fresh water heads was achieved (90 % within an absolute error of 0.46 m).We selected to adjust only the hydraulic properties of the glacial till, because the confinement of the second aquifer and its isolation from the first aquifer, caused by the very low permeability of the till (DR formation), is a clear structure that has a high influence in the flow patterns.The value giving the best approximation was a hydraulic conductivity of 0.003 m day −1 for the till.
The 721 piezometers in the area were not all equally suitable for the calibration.We made a selection based on their location (piezometers close to the model boundaries and to groundwater extractions were discarded), and based on the Figure 14 Plot showing the modelled fresh water heads versus the measured fresh water heads.amount of available data (piezometers with less than 20 measurements were also discarded).Finally, we calibrated the model with the data from 234 piezometers (Fig. 14).The calibration resulted in an absolute fresh water head difference of 0.18 m for 50 % of the piezometer and 0.46 m for 90 % of the piezometers.We considered these results satisfactory and used the calibrated model to simulate the different scenarios.
Besides the quantitative calibration, we made a qualitative validation.We used regional knowledge of the stakeholders to identify the fit between infiltration and seepage areas, rain water lens thickness, and regional groundwater flow.The main regional features of the system matched the results given by the model.

Model scenarios
The complexity of the climate change effects on the groundwater model is difficult to summarize in one scenario.Every process part of the climate change, such as an increase of the precipitation in the winter, has specific effects that can be best analyzed when modelled individually.Consequently, for the calibrated groundwater model we developed three climate scenarios.Besides the three climate scenarios, we also simulated the changes in the system due to the present situation, i.e. the autonomic processes.The four simulated scenarios are the following: 1. Autonomic scenario: we ran the model from 2005 to 2100 with the same input variables, i.e. no change in the boundary conditions and the stress terms.In such a way, we could analyze what are the changes that the system undergoes due to the still present disequilibrium.Because of this disequilibrium, saline groundwater flows upwards.This is what we call autonomic salinisation.By modelling this scenario and comparing it to   3. Sea level rise scenario for 2100 AD: this scenario illustrates a sea level rise of 0.85 m for 2100 AD estimated by the KNMI (van den Hurk et al., 2006).A sea level rise will cause an increase of the hydraulic head in the aquifers which will result in an increase of the seepage flux into the confining layer.However, the effect of sea level rise decreases rapidly with the distance from the Dutch coast as concluded by Oude Essink et al. ( 2010).
We implemented the sea level rise as a linear process.
4. Climate scenario W+ combined with sea lever rise scenario for 2100 AD: this scenario is the sum of scenarios two and three.

Model results
The model results are first described for the present situation, and then we compare the three KNMI-scenarios with the autonomous situation in 2100 AD in order to only show effects of climate change and sea level rise.From the model results we determined fresh water heads, seepage and infiltration fluxes, and salt loads to the surface water system for a  winter and a summer situation.We define seepage as the upward flow from the first aquifer into the top confining layer, and infiltration as the opposite flux.Salt loads are calculated by multiplying the seepage flux with the chloride concentration at the bottom of the confining layer.Additionally, we calculated the thickness of rainwater lenses, defined as the groundwater body with a chloride concentration less than 1.5 g l −1 .

Present situation: 2005 AD
The study area is characterized by a shallow groundwater level and a high chloride concentration of the groundwater.The average groundwater level lies at 0.5 m below the surface level in the winter and at 1 m below the surface in the summer (Fig. 15).The salinity of the groundwater at about 5 m depth is quite high at the coastal zone where it reaches values of 10 000 mg l −1 , whereas in the east of the region groundwater is fresh (Fig. 16).The thickest lenses (> 20 m) are found in the southeast of the area (Fig. 16), although in 50 % of the area rainwater lenses are thinner than 2 m or inexistent.Seepage and infiltration figures (Fig. 17) show that in the winter patterns are less defined due to the local systems, while in the summer infiltration and seepage areas can be better distinguished.Saline seepage is an important process occurring in almost 60 % of the area which can discharge up to 150 000 kg ha −1 yr −1 to the surface water system (Fig. 18).Highest salt loads occur in the low-lying polders close to the coastline.

Autonomic scenario: 2100 AD
The natural evolution of the system when no changes are applied, i.e. the autonomic evolution, shows a change in the concentration of the groundwater.Some areas become fresher and others more saline (Fig. 19).In general, the chloride concentration increases in the seepage areas because deeper and more saline groundwater is flowing to the surface.Consequently, also salt loads will increase in the future (Fig. 20), and rainwater lenses will become even thinner.The process is called autonomous salinization.Fresh rainwater lenses become thicker in most infiltration areas, because rainwater is infiltrating to greater depths and replaces the saline groundwater.The average variation of the freshwater lenses is about ± 0.25 m with a maximum variation of about ± 2 m.The effect of the autonomic processes in summer and winter    is practically the same; thus the described changes can be applied to both seasons.

Climate scenario W+: 2100 AD
In this scenario the visible seasonal effects enclose the most interesting aspect of the results.Figure 21 shows how the groundwater level will rise an average of 0.10 m in the winter, whereas in the summer it will drop an average of 0.50 m and up to 1.5 m.This will result for most areas in an increase of the salt load (Fig. 22).

Sea level rise of 0.85 m: 2100 AD
The effects of the sea level rise are most pronounced for the freshwater-heads of the first aquifer.Heads increase about 0.10 m in the first kilometre from the coastline that locally can reach up to 0.25 m (Fig. 23).Farther inland the effects of the sea level rise decrease rapidly, and, at a maximum distance of 5 to 7.5 km from the coastline, effects are less than 0.05 m.The extensional reach of the effect in this case is much smaller than the effect found in other studies in the Netherlands (Oude Essink et al., 2010) where the effects can be seen until 10 km inland of the coastline and main rivers.The explanation can be found on the shallowwater bathymetry of the Wadden Sea.Since sediments give more resistance to the transmission of water pressure than it would do in deeper sea bottom areas, the result is a more limited extensional effect.Furthermore, the increase of the heads close to the coastline causes higher seepage fluxes and therefore, higher salt loads to surface water system.Figure 23 shows that the salt load discharge can increase up to 50 000 kg ha −1 yr −1 .The summer and winter variations in this case are the same as we could see for the present situation, and are not accentuated due to the sea level rise.

Combination climate W+ and sea level rise scenario: 2100 AD
The combination of both scenarios has substantial effects on the pressure field.Heads increase in the summer and the winter close to the coastline causing higher seepage fluxes.On the other hand, the phreatic groundwater table rises in the winter and drops in the summer.The combination of higher heads in the first aquifer close to the coastline and a lower groundwater table in the summer benefits the increase of seepage in the summer, causing in this way the worst situation of the year.Just as in the previous scenario, due to saline groundwater, the salt loads towards the surface water increase (Fig. 24).Moreover, the thickness of the fresh water lenses decreases about 1 m until 10 km inland in the summer and just a few centimetres less in the winter as a result of the saline seepage (Fig. 24).As most lenses are already thin (< 2 m), the increase of seepage would probably result in problems for agriculture due to scarce fresh groundwater.

Conclusions
In this study we presented the model results of different climate scenarios on the salinisation of shallow groundwater and surface water.We successfully combined airborne EM geophysical techniques with a 3-D variable-density groundwater model.The results of the 3-D chloride field created with the airborne EM measurements go beyond the results achieved by former techniques (such as borehole-logging or groundwater sampling) to create a 3-D chloride concentration distribution of the subsurface.The reason is that airborne EM offers a higher spatial resolution of the measurements.Besides, model calculations showed that just after approximately 15 yr the numerical inaccuracies of the 3-D chloride field were removed.The inaccuracies became evident through changes in chloride concentration of several mg l −1 in relatively short periods and short distances.We could identify them as inaccuracies, because salt transport with groundwater flow was occurring too fast while it is a slow and long-term process.Introducing an incorrect chloride field (a field with an unrealistic distribution of the chloride concentration and not just a few inaccuracies) in the model will demand the model to run a simulation of up to even hundreds of years before the system reaches the equilibrium.For example, freshwater lenses under dune areas can take hundreds of years to evolve.Therefore, if the initial chloride field of the model does not include the lens, the model will have to be run until it is formed.This explains that the use of a good initial chloride field is of high importance for the simulation of the current evolution of a coastal groundwater system.With the calibrated model (50 % of the calibrated points show an absolute error smaller than 0.18 m), we calculated how the groundwater system will evolve as a consequence of the different climate scenarios.We showed that the model area is characterized by near-to-the-surface groundwater levels, high groundwater concentrations at shallow depths, and high salt loads by groundwater seepage in almost 60 % of the area.Due to autonomous processes, existing shallow freshwater lenses in the seepage areas decrease in thickness while other lenses in infiltration areas grow.The worst of the three analyzed scenarios is the one combining sea level rise and a change in the groundwater recharge.Sea level rise leads to the rise of the hydraulic heads in the first aquifer, which causes an increase of the seepage fluxes.The effect is visible in the first 5 km inland where the salt loads increase up to 50 000 kg ha −1 yr −1 .The lower groundwater recharge in the summer causes a drop of the groundwater levels.
Stakeholders such as water managers and policy makers can use the quantitative model results to (1) understand the dominant hydrogeological processes, (2) pin-point      vulnerable fresh water resources, (3) address policy tipping points in their present water management, (4) assess the impact of climate and global change, and (5), if applicable, come up with feasible adaptation strategies.Particularly, the model results are of importance in the field of water policy and agriculture.In the Netherlands there is a ranking system for the distribution of fresh water during periods of drought.The national water board is responsible for the priorities in distribution of the available water.The regional water boards, among others Wetterskip Fryslân, have the possibility to further detail the priorities for distribution of the available water to the various users.The effects of climate change might imply that some water users will simply not receive the water they need according to the ranking system.Furthermore, due to sea level rise effects in Fryslân, farmers will encounter more difficulties in the future when using surface water for irrigation purposes.The main counteracting measure would be to allow extra fresh water from Ijssel Lake into Fryslân to flush the brackish surface water more frequently.However, the results of this study stress the need of other solutions, such as storage of fresh water.On top of the regional needs for fresh water for irrigation purposes, an increasing number of farmers will encounter problems with salinisation on local scale, e.g.their own piece of land.According to the W+ scenario, the thickness of the fresh water lenses will diminish, or even completely disappear in the summer.This might be favoured by the traditional drainage systems that were designed to discharge excess water during the winter in order to allow the farmer to use machinery on his land.The results of this study show that not only water quantity aspects should be taken into account when designing drainage systems, but also other aspects such as the effects of diminishing fresh water lenses and salinisation should be incorporated.
Furthermore, this study sets an important step on the combination of airborne EM geophysical methods and groundwater modelling.It shows the importance of developing a good initial 3-D chloride field for long-term fresh-brackishsaline groundwater modelling.However, there are still inaccuracies and limitations to the method that necessarily require more research on incorporating airborne EM techniques in numerical modelling.A line of research could be to study the possibilities of increasing the depth and resolution of the EC measurements; to transform more accurately the EC measurements to chloride concentration, and to improve the initial 3-D chloride field making use of groundwater fluxes derived from numerical variable-density models.

Fig. 1 .
Fig. 1.The position of the study area within the Netherlands showing the surface elevation (left panel), and a topographic map showing the surface water, the urban areas and the Wadden Sea (right panel).

Figure 2
Figure 2 Geological cross-profile showing the main geological units.For explanation of the abbreviation of the geological units see the main text.For the location of the profile see Fig.3.

Fig. 2 .
Fig. 2. Geological cross-profile showing the main geological units.For explanation of the abbreviation of the geological units, see the main text.For the location of the profile, see Fig. 3.

Figure 3
Figure 3 Study area with location of Cl -measurements and HEM and S the location of the geological cross-profile is also indicated.

Figure 4 Figure 3
Figure 4 Depth slices and cross-sections of EC(total) from the AEM mod

Figure 4
Figure 4 Depth slices and cross-sections of EC(total) from the AEM models

Fig. 3 .
Fig. 3. Study area with location of Cl − measurements and HEM and SkyTEM survey lines; the location of the geological cross-profile is also indicated.

-Figure 3
Figure 3 Study area with location of Cl -measurements and HEM and SkyTEM survey lines; the location of the geological cross-profile is also indicated.

Figure 4
Figure 4 Depth slices and cross-sections of EC(total) from the AEM models

Fig. 4 .
Fig. 4. Depth slices and cross-sections of EC (total) from the AEM models.

Figure 5 Fig. 5 .
Figure 5 Probability of glacial till, as derived from AEM models (in this case the results of the HEM are shown) Fig. 5. Probability of glacial till, as derived from AEM models (in this case the results of the HEM are shown).

Figure 5 Figure 6
Figure 5 Probability of glacial till, as derived from AEM models (in this case the results of the HEM are shown)
The vertical discretization from top to bottom is as follows: 14 model layers of 0.5 m, 2 model layers of 1 m, 10 model layers of 2 m, 5 model layers of 5 m, 7 model layers of 10 m, and 4 model layers of 35 m (Fig. 6).The boundaries of the model are shown in Fig. 7.The drainage network and the recharge of groundwater define the top boundary.The geohydrological basis is situated at Figure 7 3-D view (left) and top view (right) showing the model boundaries.

Figure 8
Figure 8 Lithologies in the 3D voxel model of the Holocene sequence; one out of the 50 equiprobable realizations, as calculated with the stochastic simulation technique, is shown.

Figure 7
Figure 7 3-D view (left) and top view (right) showing the model boundaries.

Figure 8
Figure 8 Lithologies in the 3D voxel model of the Holocene sequence; one out of the 50 equiprobable realizations, as calculated with the stochastic simulation technique, is shown.

Fig. 8 .
Fig. 8. Lithologies in the 3-D voxel model of the Holocene sequence; one out of the 50 equi-probable realizations, as calculated with the stochastic simulation technique, is shown.

23Figure 9
Figure 9 Relation between groundwater resistivity and reciprocal of the apparent Formation Factor.

Figure 10
Figure 10 3-D chloride distribution field in the study area obtained after the interpolation of the measurements.

Fig. 9 .
Fig. 9. Relation between groundwater resistivity and reciprocal of the apparent formation factor.

Figure 9
Figure 9 Relation between groundwater resistivity and reciprocal of the apparent Formation Factor.

Fig. 11 .
Fig. 11.The flow of data and interactions for deriving vertical and horizontal hydraulic conductivities and chloride concentration for the groundwater model.

Fig. 12 .
Fig. 12. Change in the volume of fresh water in the model, consequence of the inaccuracies in the 3-D field.Every line represents one slice of the model, one cell thick.

Figure 13
Figure 13 Groundwater recharge for the summer and the winter of the present situation (left) and for the summer and the winter of the W+ scenario in the year 2050 (right).

Fig. 13 .
Fig. 13.Groundwater recharge for the summer and the winter of the present situation (left panels) and for the summer and the winter of the W+ scenario in the year 2050 (right panels).

25Figure 13
Figure13Groundwater recharge for the summer and the winter of the present situation (left) and for the summer and the winter of the W+ scenario in the year 2050 (right).

F
i g ur e 1 4 P l ot s h o w i n g t h e m o de l l e d f r e s h w a te r h e a d s v e r s u s th e m e a s u r e d f r e s h w a t e r heads.

Fig. 14 .
Fig. 14.Plot showing the modelled fresh water heads versus the measured fresh water heads.

Fig. 15 .
Fig. 15.Phreatic groundwater level in the winter (left panel) and in the summer (right panel), both in 2005.

Fig. 16 .
Fig. 16.Chloride concentration under the first confining layer in the winter and the summer of 2005 (left panel), and thickness calculated from the surface elevation of the fresh water lenses with Cl − < 1.5 g l −1 in the summer of 2005 (right panel).

Fig. 17 .
Fig. 17.Seepage and infiltration patterns under the first confining layer in the winter (left panel) and the summer (right panel) of 2005.

Figure 18
Figure 18 Salt load discharging from the groundwater into the surface water in the summer of 2005.

Fig. 18 .
Fig. 18.Salt load discharging from the groundwater into the surface water in the summer of 2005.

Fig. 19 .
Fig. 19.Left panel: change in the thickness of the fresh water lenses in the summer autonomic scenario (year 2100) and present case (year 2005).Red indicates that the thickness of the fresh water lenses decreases, and blue indicates that it increases.Right panel: change in concentration under the first confining layer between the autonomic scenario (year 2100) and the present case (year 2005).Red colours show salinisation and blue colours freshening of the groundwater.

Figure 20
Figure 20 Change in salt load between the autonomic scenario (year 2100) and the present case (year 2005) in the summer.Red colours indicate an increase of the salt load and blue colours a decrease.

Fig. 20 .
Fig. 20.Change in salt load between the autonomic scenario (year 2100) and the present case (year 2005) in the summer.Red colours indicate an increase of the salt load and blue colours a decrease.

Fig. 21 .
Fig. 21.Groundwater level rise in the winter in the climate scenario W+ compared to the autonomic scenario in 2100 (left panel), and groundwater level drop in the summer in the climate scenario W+ compared to the autonomic scenario in 2100 (right panel).

Figure 21
Figure 21 Groundwater level rise in the winter in the climate scenario W+ compared to the autonomic scenario in 2100 (left), and groundwater level drop in the summer in the climate scenario W+ compared to the autonomic scenario in 2100 (right).

Figure 22
Figure 22 Salt load decrease (blue) and increase (red) in the climate W+ scenario compared to the autonomic scenario in 2100.

Figure 23
Figure23Increase in hydraulic head in the sea level rise scenario compared to the autonomic scenario (left), and increase in salt load (red colours) in the sea level rise scenario compared to the autonomic scenario (right).

Fig. 22 .
Fig. 22. Salt load decrease (blue) and increase (red) in the climate W+ scenario compared to the autonomic scenario in 2100.

Fig. 23 .
Fig. 23.Increase in hydraulic head in the sea level rise scenario compared to the autonomic scenario (left panel), and increase in salt load (red colours) in the sea level rise scenario compared to the autonomic scenario (right panel).

Fig. 24 .
Fig. 24.Left panel: decrease (red) and increase (blue) of the fresh water lens thickness in the climate W+ and sea level rise scenario compared to the autonomic scenario.Right panel: salt load increase (red) and decrease (blue) due to climate W+ and sea level rise in comparison with the autonomic scenario.