Paleo-modeling of coastal saltwater intrusion during the Holocene : an application to the Netherlands

Coastal groundwater reserves often reflect a complex evolution of marine transgressions and regressions, and are only rarely in equilibrium with current boundary conditions. Understanding and managing the present-day distribution and future development of these reserves and their hydrochemical characteristics therefore requires insight into their complex evolution history. In this paper, we construct a paleo-hydrogeological model, together with groundwater age and origin calculations, to simulate, study and evaluate the evolution of groundwater salinity in the coastal area of the Netherlands throughout the last 8.5 kyr of the Holocene. While intended as a conceptual tool, confidence in our model results is warranted by a good correspondence with a hydrochemical characterization of groundwater origin. Throughout the modeled period, coastal groundwater distribution never reached equilibrium with contemporaneous boundary conditions. This result highlights the importance of historically changing boundary conditions in shaping the present-day distribution of groundwater and its chemical composition. As such, it acts as a warning against the common use of a steady-state situation given present-day boundary conditions to initialize groundwater transport modeling in complex coastal aquifers or, more general, against explaining existing groundwater composition patterns from the currently existing flow situation. The importance of historical boundary conditions not only holds true for the effects of the largescale marine transgression around 5 kyr BC that thoroughly reworked groundwater composition, but also for the more local effects of a temporary gaining river system still recognizable today. Model results further attest to the impact of groundwater density differences on coastal groundwater flow on millennial timescales and highlight their importance in shaping today’s groundwater salinity distribution. We found free convection to drive large-scale fingered infiltration of seawater to depths of 200 m within decades after a marine transgression, displacing the originally present groundwater upwards. Subsequent infiltration of fresh meteoric water was, in contrast, hampered by the existing density gradient. We observed discontinuous aquitards to exert a significant control on infiltration patterns and the resulting evolution of groundwater salinity. Finally, adding to a long-term scientific debate on the origins of groundwater salinity in Dutch coastal aquifers, our modeling results suggest a more significant role of pre-Holocene groundwater in the present-day groundwater salinity distribution in the Netherlands than previously recognized. Though conceptual, comprehensively modeling the Holocene evolution of groundwater salinity, age and origin offered a unique view on the complex processes shaping groundwater in coastal aquifers over millennial timescales. Published by Copernicus Publications on behalf of the European Geosciences Union. 3892 J. R. Delsman et al.: Paleo-modeling of coastal saltwater in the Netherlands


Introduction
While fresh groundwater reserves in coastal areas are a vital resource for millions of people, they are vulnerable to salinization, given both their proximity to the sea and the usually large demands on freshwater by the larger population densities in coastal areas (Barlow and Reichard, 2009;Custodio and Bruggeman, 1987;Ferguson and Gleeson, 2012;Post and Abarca, 2009;Werner et al., 2013).Reported impacts of salinizing coastal aquifers include the salinization of abstraction wells (Custodio, 2002;Stuyfzand, 1996), decrease of agricultural yield (Pitman and Läuchli, 2002), degrading quality of surface waters (De Louw et al., 2010), and adverse effects on vulnerable ecosystems (Mulholland et al., 1997), issues that will only intensify in the future, given the prospects of global change (Kundzewicz et al., 2008;Oude Essink et al., 2010;Ranjan et al., 2006).The above issues have sparked a surge in renewed scientific interest in the "classic" saltwater intrusion process, i.e., the development of a landward-protruding saline groundwater wedge under the influence of groundwater density differences, as reviewed by Werner et al. (2013).
Given their vulnerability, sustainable management of coastal fresh groundwater reserves is of paramount importance.A prerequisite is an accurate description of the presentday distribution of fresh groundwater reserves.That accurate description is, however, difficult to obtain: measurements are sparse, especially at greater depths, while salinity varies within short distances, driven by relatively minor head gradients that vary over time (De Louw et al., 2011).And although recent advances in airborne geophysics (Faneca Sànchez et al., 2012;Gunnink et al., 2012;Siemon et al., 2009;Sulzbacher et al., 2012) are promising, the availability of airborne data is still limited and its reliability decreases with depth.Variable density groundwater modeling may be used to assess coastal freshwater resources and management strategies (e.g., Nocchi and Salleolini, 2013;Oude Essink et al., 2010).However, as a result of the density feedback of solute concentration on groundwater flow, this requires an adequate description of the initial solute concentration: a vicious circle of having to know the salinity distribution to model the salinity distribution.A frequent workaround is the assumption of steady state, obtained by a spin-up period applying current boundary conditions (e.g., Souza and Voss, 1987;Vandenbohede et al., 2011;Vandenbohede and Lebbe, 2002).However, given the usually long timescales involved, coastal groundwater systems are rarely in equilibrium, often still reflecting events occurring thousands or even millions of years ago (e.g., Groen et al., 2000;Post et al., 2003;Stuyfzand, 1993).
Paleo-hydrogeologic modeling, or the transient modeling of the long-term co-evolution of landscape and groundwater flow, may provide a way out of this vicious circle.This involves starting a model run at a reference point in time where the salinity distribution is either more or less known, or is certain not to influence the present-day salinity distribution.Successful use of paleo-hydrogeologic modeling is difficult however, given the long timescales considered, the often limited availability of data on paleo-boundary conditions and the impossibility of validating past time frames (Van Loon et al., 2009), on top of the "normal" difficulties in hydrogeologic (transport) modeling (Konikow, 2010).Paleo-hydrogeologic modeling has been previously applied to study the influence of groundwater during glacial cycles (Bense and Person, 2008;Lemieux and Sudicky, 2009;Person et al., 2012;Piotrowski, 1997), to better explain the observed pattern in groundwater ages using carbon dating (Sanford and Buapeng, 1996), to study the degradation of fen areas in the Netherlands (Van Loon et al., 2009;Schot and Molenaar, 1992), and to relate archeological settlements to historic phreatic groundwater levels (Zwertvaegher et al., 2013).Applications of paleo-hydrogeologic modeling in variable-density flow situations are scarce however, and are limited to the evolution of fresh-and saltwater over the last century (Nienhuis et al., 2013;Oude Essink, 1996) or millennium (Lebbe et al., 2012;Vandevelde et al., 2012), using available historic information.
In this paper, we apply paleo-hydrogeologic modeling to study the processes controlling the Holocene evolution of groundwater salinity in a representative deltaic coastal aquifer: the coastal region of the Netherlands.The studied region (Sect.2) has a complex paleo-geographic history of marine trans-and regressions, peat accumulation and degradation, and more recently land reclamation, drainage and groundwater abstraction.The groundwater salinity distribution still reflects this complex history (Oude Essink et al., 2010;Post et al., 2003;Stuyfzand, 1993), and both the paleogeographic evolution (Vos et al., 2011) and the distribution of aquifer properties (Weerts et al., 2005) are relatively well known.As such, the region is well suited to a paleohydrogeologic modeling approach.In addition, societal interest in the region's groundwater salinity distribution is spurred by a deterioration of surface water quality through exfiltration of brackish groundwater, adversely affecting agriculture and vulnerable ecosystems (De Louw et al., 2010;Oude Essink et al., 2010;Van Rees Vellinga et al., 1981).While salinity is the prime focus of the present paper, the approach presented is considered relevant for the many other societally relevant groundwater constituents in coastal aquifers, like nutrients (Van Rees Vellinga et al., 1981;Stuyfzand, 1993) or arsenic (Harvey et al., 2006;Michael and Voss, 2009).

Study area
We studied an approximately west-to-east-oriented transect, located some 10 km south of the city of Amsterdam, the Netherlands (Fig. 1).The 65 km long transect is oriented perpendicular to the coastline and extends from 12 km offshore to the midpoint of an ice-pushed ridge, forming a regional groundwater divide.The transect is exemplary for this part of the coastal region of the Netherlands, intersecting coastal sand dunes, reclaimed lakes, managed fen areas and the aforementioned ice-pushed ridge.Elevations along the transect range from 5 m below mean sea level (b.s.l.) in the deep polder areas, to locally 35 and 30 m above mean sea level (m.s.l.) for the dune area and ice-pushed ridge, respectively.Present-day climate is categorized as moderate maritime, with temperatures that average 10 • C, an average annual precipitation total of 840 mm, and an average annual Makkink reference evapotranspiration total (Makkink, 1957) of 590 mm (KNMI, 2010).The hydrogeology of the area is characterized by 300 m thick deposits of predominantly Pleistocene marine, glacial and fluvial deposits, forming alternating sandy aquifers and clayey aquitards (Fig. 1b).The Maassluis formation comprises the oldest Pleistocene deposits and includes sandy and clayey sediments of marine origin, limited dated samples indicate remaining connate marine groundwater (Post et al., 2003).An aquiclude of Tertiary clays is present below these deposits (Dufour, 2000).Excluding the coastal dune area, Holocene deposits are generally no more than 10 m thick, thinning out in an easterly direction.A more elaborate description of these Holocene deposits and their genesis is presented below.Present-day groundwater flow is directed from the elevated dune and ice-pushed ridge areas towards the deep polder areas in the center of the transect.Water management in the central part is aimed at keeping groundwater levels at an optimal level for agriculture, within 1-2 m below the ground surface, and requires an extensive network of canals, ditches and subsurface drains to drain excess precipitation and exfiltrating groundwater.Flow direction reverses during summer, when freshwater from the river Rhine is redirected to compensate for precipitation deficits and salinity increases.

Holocene paleo-geographical development
An overview of the Holocene paleo-geographical development of the area is presented in Fig. 2. At the end of the Pleistocene, up to about 13 000 BC, the area was characterized by sandy plains with braided rivers, sloping gently from the ice-pushed ridge towards the contemporaneous coastline.Because of post-glacial sea level rise during the early Holocene, groundwater levels started to rise in the coastal Red dots and letters in (g) refer to the corresponding paleogeographical map (a-f).For reference, note that the extent of the paleo-geographical maps equals the extent of Fig. 1a.
zone and promoted the widespread formation of peat.The continuing sea level rise resulted around 6500 BC in the submersion of these peat deposits, when an open barrier system with barriers and a tidal basin formed to a maximum extent of about three-quarters of the studied transect (transgression phase).Around 3950 BC the Dutch coast became a closed system, when sediment availability had begun to match the decreased sea level rise rate (regression phase).
The coast now changed into a prograding system that extended into the North Sea until 2500 BC.The tidal areas silted up and freshened, stimulating large-scale peat development behind the coastal barriers.Peat development was at a maximum around AD 1000, reaching a maximum thickness of 6 m (elevation of 2 m m.s.l.).Subsequently, peat extraction and anthropogenic drainage resulted in rapid peat degradation, a lowering of the ground surface and the eventual formation of several freshwater lakes.Increased sand availability in the coastal zone around AD 900 led to the formation of an extended and higher coastal dune system (Jelgersma et al., 1970).Anthropogenic influence grew in importance from AD 1500 onwards, through land cultivation, improved agricultural drainage, river embankment and urban development.Large-scale land reclamation projects were carried out on most freshwater lakes in the 19th century, resulting in the deep polders of Haarlemmermeer (AD 1852) and Horstermeer (AD 1888).Groundwater abstraction started in the coastal dunes and the ice-pushed ridge in the mid-1800s.Subsequent salinization problems prompted the abandon-ment of most abstraction wells in the coastal dunes in favor of the current Rhine water infiltration scheme in use since 1957, whereby water is infiltrated in infiltration ponds and extracted using recovery canals and horizontal drains.

Hydrochemical facies analysis (HYFA)
In the 1980s, about 20 piezometer nests in the western part of the studied transect, each with four to fifteen 1 m long monitor well screens, were sampled and analyzed on main constituents, trace elements and environmental tracers (locations and depths in Fig. 6a).Stuyfzand (1993Stuyfzand ( , 1999) ) used the resulting data set (with many more data from monitor wells along the Dutch coast) to depict the spatial distribution of groundwater bodies with a specific origin (hydrosomes), and their hydrochemical facies (distinct hydrochemical zones within each hydrosome).Environmental tracers (Cl / Br ratio, 18 O, 3 H, 14 C, SO 4 and HCO 3 ) were used to discern the following hydrosomes, in order of increasing salinity: (i) fresh dune groundwater (rainwater infiltrated in coastal dunes) (in Fig. 6d), (ii) fresh, artificially recharged Rhine River water, (iii) slightly brackish polder water (a mix of rainwater, Rhine River water and exfiltrated Holocene transgression water), which after mixing infiltrated via canals and ditches on a higher topographical level than the deep polders from which the mix originated (P), (iv) two types of brackish groundwater, which infiltrated during the Holocene transgression (LC / Lm), (v) brackish to saline paleo-groundwater upconing from deep marine Tertiary aquitards (M), and (vi) intruding North Sea water (S) (Fig. 6a).Within each hydrosome a variety of hydrochemical facies was discerned by a combination of four aspects:

Paleo-hydrogeological modeling
To model the evolution of groundwater salinity throughout the Holocene, we used the variable density groundwater modeling software SEAWAT (Langevin and Guo, 2006) to set up a 2-D model for the described transect (A-B in Fig. 1a, conceptual outline in Fig. 3).We assumed a Dirichlet boundary condition (sea level) on the western side, and no-flow boundaries (groundwater divide and geohydrological base, taken as the top of Tertiary clays (below the Maassluis formation, fifty-one, 100 m wide model cells in the horizontal, and 102 layers in the vertical, whose thicknesses increase with depth (thickness 1 m in upper 60 m, increasing to 10 m at maximum depth) in the vertical.Cell-specific geohydrological properties were taken from national geohydrological databases REGIS (Vernes and Van Doorn, 2005) and GEOTOP (Van der Meulen et al., 2013;Stafleu et al., 2011) (both available from http://www.dinoloket.nl).GEOTOP provides detailed (100 × 100 × 0.5 m) estimates of lithology to a depth of 50 m, we applied REGIS-derived formation-specific hydraulic properties for the deeper subsurface.We assumed a homogeneous aquifer seaward of the present coastline, given the limited availability of geohydrological information.Information on present-day water management was obtained from the Netherlands Hydrological Instrument model (De Lange et al., 2014, available from http://www.nhi.nu).
We used chloride to represent salinity, as chloride is the dominant anion in Dutch coastal groundwater and density is linearly related to it within naturally occurring concentrations.To better understand the evolution of groundwater salinity, we included several fictitious inert tracers, representing the various inputs to the groundwater system, as additional mobile species in the simulation.These tracers were given a concentration of one when they entered the model domain, or were present during model start.We refer to these fictitious tracers as origin tracers in the remainder of this paper.Furthermore, we modeled direct groundwater age (Goode, 1996) by including an additional specie with a negative zero-th order decay term (Zheng, 2009).This specie is zero when entering the model domain, then gains one for every year spent inside the model.An overview of the different origin tracers and their relation to hydrosomes (Sect.2.3) is presented in Table 1.Longitudinal dispersivity was set to 1 m, the lower bound found for similar settings in experimental work reviewed by Gelhar et al. (1992), and similar to values used in comparable settings (Lebbe, 1999;Oude Essink et al., 2010).Horizontal and vertical transversal dispersivities were assumed 0.1 and 0.01 m, respectively (Zheng and Wang, 1999), and we assigned a uniform molecular diffusion coefficient of 10 −9 m 2 s −1 .We did not attempt to calibrate our model, recognizing that calibration would only be possible for the most recent periods, and a rigorous sensitivity analysis was impossible given the long calculation times.We regard our model therefore primarily as a conceptual tool.Still, we assessed the validity of the model by comparing model results to measured heads and chloride measurements, tritium-derived groundwater ages and a hydrochemical interpretation of groundwater origin (HYFA, see Sect.2.3).Available radiocarbon measurements were proven impossible to use for accurate dating in this area, due to the large contribution of heterogeneously aged sedimentary carbon sources to inorganic carbon dissolved in groundwater (Post, 2004).
The geographical changes throughout the Holocene were implemented using 10 successive time slices, with each time slice representing a distinct period in the paleo-geographical evolution (Fig. 2, Table 2).Model start was set at 6500 BC, marking the start of marine influence in the area.Conditions during a time slice were assumed constant, with the exception of the rapidly rising sea level (implemented using 10 stress periods) during the first 2 time slices (transgression phase).The model state (head and concentration) at the end of each time slice was used as the starting state for the subsequent time slice; model cells not present in a previous time slice were given the state of the previously uppermost model cell per column.Specific to each time slice were its sea level, surface elevation, near-surface geohydrological properties, drainage structure and groundwater abstractions, which were reconstructed based on the depositional history reflected in the near-surface geological record and various literature sources (Table 3).In addition to reconstructed largerscale drainage structures such as e.g. the river Vecht, we applied infinite drainage at the model surface to represent small streams and creeks.As little erosion has taken place since the start of modeling, and compaction of clay was not considered significant, the near-surface geological record provided a good approximation of the historical landscape.An   Sea level rise Beets et al. (2003), Denys and Baeteman (1995), Jelgersma (1961), Kiden (1995), Ludwig et al. (1981), Van de Plassche (1982) Geohydrological Reclaimed areas Dufour (2000), Schultz (1992) Groundwater abstractions Van Loon (2010), Oude Essink (1996) important exception is the build-up and subsequent degradation of peat domes; we derived model parameters for peat elevations and extent from a detailed reconstruction located in a similar setting just north of Amsterdam (Vos, 1998).Geohydrological properties for historical surface sediments were assumed to equal their current (buried) properties, except for uncompacted peat deposits, set in accordance to relevant literature values (Kechavarzi et al., 2010).
As no long-term precipitation record exist for the Netherlands, and annual temperatures have remained approximately constant over the past 7 kyr (Davis et al., 2003), we chose to apply a constant uniform recharge of 0.7 mm d −1 , equal to the current long-term average, for all time slices.We did not differentiate recharge amounts for vegetation types, given the lack of data on past vegetation patterns.North Sea chloride concentration was kept at a constant 16 g L −1 over the entire model time, the present average concentration.Initial chloride concentration (at 6500 BC) was set to 16 g L −1 below the area initially inundated by the sea, and zero throughout the remainder of the model, as the 100 kyr period of the Weichselian glacial stage preceding the modeled period is expected to have caused extensive freshening of the Pleistocene aquifers (Post et al., 2003).The only exception is the low permeable Maassluis formation at the bottom of the transect, where limited dated samples indicate only partial freshening of this connate marine groundwater (Post et al., 2003).The initial concentration in the Maassluis formation was therefore assumed to be 10 g L −1 , the approximate upper limit of measured concentrations (Stuyfzand, 1993).In addition to the described scenario, four additional sensitivity runs were performed to explore two main model uncertainties: dispersivity and the chloride concentration of water present in the Maassluis formation at model start.Dispersivities were decreased 10-fold, and the initial Maassluis concentration was set to 0, 5 and 15 g L −1 in these runs, respectively, all other parameters remaining unchanged.

Sensitivity runs
SEAWAT model run time simulating 8.5 kyr and including seven additional mobile species was approximately 3 days on a standard single cpu; no convergence problems were encountered, including at time-slice transitions.The performed sensitivity runs did not reveal a significant influence of dispersivity on the general shape of the present-day salinity distribution.A 10-fold decrease in dispersivity did, however, result in a narrower transition zone between fresh and saline groundwater, most importantly beneath the coastal dunes (100 to 25 m).Further inland, the width of the transition zone decreased from around 50 m to around 25 m.While measurements indicate a narrow transition zone beneath the coastal dunes, they suggest a wide transition zone further inland and are therefore inconclusive in regards to the "better" value.The Maassluis sensitivity runs revealed a clear dependence of the modeled historical trajectory and presentday location of Maassluis water on its initial chloride concentration.While the general processes acting on this water type remain the same, the extent to which Maassluis water is transported through the subsurface is negatively related to its chloride concentration.At 15 g L −1 , the density difference between infiltrating transgression water and resident Maassluis water clearly provides less incentive for flow than at 0 g L −1 .As a result, while a significant fraction of Maassluis water is still present at its original location in the 15 g L −1 run at model end, it is almost completely displaced in the 0 g L −1 run, with the 5 and 10 g L −1 scenarios in between those extremes.Comparison with the HYFA results of Stuyfzand (1993), although based on only limited samples at the relevant depths, suggests an initial Maassluis concentration of around 10 g L −1 , showing agreement in the relatively shallow occurrence of Maassluis at a depth of 100 m around x coordinate 112 km.

Model validation
We compared modeled heads with averaged heads measured in 382 piezometers located along the modeled transect (RMSE of 1.4 m, Fig. 4, locations in Fig. 1).Visual inspection revealed a concentration of absolute errors in the coastal dune area: the RMSE of modeled versus measured heads excluding the dune area is a mere 0.7 m.RMSE normalized to the range of observed values (NRMSE) is in both cases 11 %.Larger absolute deviations in the coastal dune area were expected given the large head variation over short distances caused by the varied relief, the concentration of well fields and the presence of an artificial recharge installation, factors implemented in the model only in a simplified way.Chloride measurements are available in the study area from AD 1891 onwards, which we compared to concomitant model results (RMSE, NRMSE of 2.7 g L −1 , 16 %, respectively, 1.5 g L −1 , 22 % excluding the coastal dunes, locations in Fig. 1).These relatively large RMSEs reflect the difficulty in obtaining good fits between measured and modeled values along a transect due to the large spatial variation in (modeled) chloride concentrations over short distances.Nevertheless, the measured groundwater chloride distribution is approximated quite well (Fig. 5), with the model capturing both the depth of the Badon Ghijben-Herzberg (BGH)-type lens below the coastal dunes and the upward movement of chloride below the inland deep polder areas, even including the occurrence of very localized brackish groundwater below the Horstermeer.Discrepancies around x coordinates 113 and 125 are likely due to the orthogonal projection of chloride measurements on the model transect, and caused by local head gradients not included in the chosen model transect.And while the several tens-of-meters-wide transition zone between fresh and saline water inland is well-captured by the model, the width of the transition zone beneath the BGH lens below the dunes is overestimated by the applied model-wide dispersion parameters.A comparison of available tritium measurements with the modeled 1952 age contour was hampered by the necessary orthogonal projection of the sparsely available tritium measurements on the modeled transect.Notwithstanding, tritium ages generally confirmed the modeled vertical extent of both post-1952 groundwater in infiltration areas around the deep polder Horstermeer, the coastal dunes and the ice-pushed ridge area, and pre-1952 groundwater beneath the exfiltrating deep polder Haarlemmermeer.Age calculation results were omitted from this paper for the sake of brevity, but are included in the film available as supplementary information.
We compared the present-day distribution of modeled origin tracers to results of a hydrochemical facies analysis (HYFA) (Sect. 2.3;Stuyfzand, 1993Stuyfzand, , 1999) ) along the model transect, approximately between x coordinates 95 and 110 km (Fig. 6).HYFA uses the hydrochemistry of groundwater to identify groundwater bodies (hydrosomes) and the hydrochemical zones within them, and thus provides clues to their respective histories.This comparison therefore provides a comprehensive, independent model test.Various discrepancies exist between the distribution of modeled origin tracers and measured distribution of hydrosomes.Discrepancies are largest at depths greater than 150 m b.s.l.beneath the coastal dunes, and below 80 m b.s.l. in the deep polder area.This is due to the depth limits of observation wells (differing in both zones), giving rise to uncertainties in both the hydrochemical, hydraulic and hydrogeological data underlying the HYFA analysis and in the structure of the model.The differences call for further model optimizations where the hydrochemical patterns are very reliable, for instance regarding the advance of the intruding North Sea water, shape of the fresh dune water lens and transition zone below it, and the presence of polder water in the central parts of the deep polder Haarlemmermeer.Overall, however, the comparison shows a clear correspondence between the position of modeled origin tracers and hydrosomes, in both relatively recent (seawater wedge, infiltrating dune water) and older water types (Maassluis water, and water infiltrated during transgression and after extensive peat formation).

Evolution of groundwater salinity
An overview of the modeled evolution of the groundwater chloride distribution is presented in Fig. 7; a film of the evolution of chloride, age and origin tracers is available as supplementary information.Before 4500 BC, the coastline shifted gradually landward to a maximum of about threequarters of the model transect (x coordinate 129 km), receding to x coordinate 125 km in 3300 BC.Saline water infiltrated below the zone of marine influence through free convection, showing classic fingering patterns (Elder, 1967).Infiltration of marine water was influenced significantly by the presence of aquitards between 50 and 100 m b.s.l. and below 150 m b.s.l.Infiltration in the absence of these aquitards (around x coordinate 105 km) was rapid, reaching a depth of 150 m within decades.However, where infiltration water encountered aquitards, the concentration gradient driving free convection and hence infiltration rates decreased.Infiltration water subsequently expanded horizontally (e.g., Fig. 7c, x coordinate 115 km), forcing the resident freshwater to flow upwards, resulting in an effective stop to infiltration in this region.Although salinization rates were significantly lower in low-permeable strata, the time available was enough to also completely salinize the aquitards between 50 and 100 m b.s.l.Infiltration of transgression water halted completely after the coastal barriers closed around 3300 BC and the inland area freshened.Salinization of deeper layers continued as density differences resulted in a further downward movement and horizontal spreading of the infiltrated water.A BGH-type freshwater lens developed underneath the coastal barriers within centuries and expanded when the elevation of the coastal dunes increased significantly around AD 900.The widespread gradual build-up of peat between 3300 BC and AD 1500 elevated the land between the coastal dunes and the ice-pushed ridge to approximately 3 m above contemporaneous sea level, resulting in extensive infiltration of meteoric freshwater.Maximum infiltration depth varied between 20 and 100 m and was reached after about 4 millennia of infiltration.Rates of forced convection of meteoric infiltration water were much lower than those observed for free convection of transgression water, especially where low-permeable strata and underlying saline groundwater impedes downward flow.Infiltration rates were considerably reduced after AD 1500, caused by peat degeneration and the resulting decrease in hydraulic gradients.After the development of the Vecht river system around AD 500 (Fig. 7e, x coordinate 130 km), the gaining river Vecht started to attract regional groundwater flow resulting in upconing of the deeper brackish groundwater to about 50 m b.s.l.Lake reclamations around AD 1850 caused large vertical upward hydraulic gradients, shifting groundwater flow patterns towards the deep polders.The reclamation of the Horstermeer transported the upconed water beneath the Vecht river eastward and further upwards, where it eventually exfiltrated.Effects of the reclamation of the Haarlemmermeer are most pronounced where the aquitards between 50 m and 100 m b.s.l. are absent, and brackish groundwater exfiltrates (between x coordinate 103 and 108 km).There, groundwater flows upward at a rate of approximately 0.5 m yr −1 , momentarily slowed down between 1900 and 1950 by the overextraction in the dune area.Groundwater below the remaining part of the Haarlemmermeer is still fresh up to a depth of 50 m b.s.l., owing to the lower upward flow rates due to the presence of aquitards and larger distance to the coastal dunes.The depth and shape of the BGH-type freshwater lens beneath the coastal dunes is not straightforwardly determined by the BGH approximation, but is influenced by both the occurrence of lowpermeable sediments below the dunes and groundwater abstraction.Moreover, the freshwater lens is displaced eastward as a result of the steep hydraulic gradient towards the Haarlemmermeer.

Total salt present
We calculated the total amount of chloride present in the model (SP), east of the present-day coastline, and the contribution of the different sources of chloride to SP (Fig. 8).SP of transgression origin was the dominant input of chloride (60 %), Maassluis origin contributed 24 %, while seawater intrusion accounted for 17 % of SP.SP predominantly increased before 4500 BC, owing to the infiltration of transgraded seawater.As the sea level rose and the coastline progressively moved landward, the passing of discontinuities in these aquitards caused a distinctive accelerating and decelerating pattern in the increase of SP in the model.Infiltration of transgression water reduced after 4500 BC, when the continued deposition of marine clays started to pose a considerable resistance to vertical flow.Infiltration of transgression water halted completely after the coastal barriers closed around 3300 BC, the continuing increase of SP until 2100 BC is due to a landward migration of coastal barriers.Subsequent infiltration of meteoric water slowly decreased SP until around AD 1500.Extensive drainage, land reclamation and groundwater abstraction resulted in the reversal of groundwater gradients and promoted the landward migration of seawater intrusion, increasing SP in the groundwater = L -young marsh type, M = Maassluis, P = polder, S = (actual) North Sea.See Sect.2.3 for a description of discerned hydrosomes, and Table 1 for the mapping of origin tracers to hydrosomes, dots and dashed lines in (a) denote locations of piezometers used in HYFA.system.The steep increase in SP of sea origin between 1900 and 1950 was caused by groundwater (over-) abstraction in the coastal dunes.

Origin tracers
The modeled origin tracers provided a comprehensive overview of the evolution of groundwater types of differ- ent origin (Fig. 9).As already noted, widespread infiltration of transgression water occurred through free convection between 6500 and 3300 BC.The infiltration had a pronounced effect on other tracers, as resident older water was pushed upwards around the infiltrating fingers.This mechanism is most readily apparent in the evolution of Maassluis water, mobilized from its original position below 170 m b.s.l.upwards to the ground surface (Fig. 9b).After marine infiltration stopped, the mobilized Maassluis water moved gradually downward, while the Maassluis water still present at greater depths was slowly displaced upwards by the more saline transgression water due to density differences.This of polder Haarlemmermeer, put forward by Stuyfzand (1993) to explain the present-day position of Maassluis water, was shown to be only of minor influence.The infiltration of meteoric infiltration water caused the submarine groundwater discharge of earlier transgression water seaward of the coastal dunes (Fig. 9d, x coordinate 95 km), owing to a regional seaward flow system underlying the coastal dunes.The influence of the classic seawater intrusion mechanism, the landward-intruding seawater wedge (origin tracer sea), was limited to about 5 km from the coastline, in general agreement with results of Stuyfzand (1993).Maassluis water constituted the third input of chloride, mainly contributing to chloride concentrations in the eastern part of the model transect.

Discussion
While the paleo-geographical development of the Netherlands is relatively well documented, based on numerous investigations of the near-surface geology, carbon dating and archeological evidence (Vos and Gerrets, 2005;Vos et al., 2011;Weerts et al., 2005), such reconstructions necessarily entail a significant degree of interpretation and hence uncertainty.This uncertainty is only enlarged in a paleohydrogeological model, both through model simplification and the introduction of additional uncertain parameters (e.g., climate, sea salinity, historical surface elevation, and hydraulic parameters of the subsurface).In addition, the employed 2-D approach neglects variation perpendicular to the model transect, and the use of time slices introduces discrete and significant jumps in the gradually changing geology and hydrology.Finally, the chosen model discretization prevents the inclusion of boils, highly localized conduits penetrating the low-permeable Holocene strata (De Louw et al., 2010) that contribute significantly to the exfiltration of salts in the polder Haarlemmermeer (Delsman et al., 2013).We did not attempt to quantify the uncertainty in our model, and therefore do not claim its validity other than as a conceptual tool.Nevertheless, judging from comparisons to measured heads, chloride concentrations, age patterns and, perhaps most assuring, hydrochemical facies analysis, the model appears to explain the present-day distribution of both groundwater salinity and origin quite well.We are confident, therefore, that the important processes occurring over the modeled period, responsible for the present-day salinity distribution on the scale considered, are well represented by the paleohydrogeological model.
The influence of variable density on groundwater flow patterns has been widely demonstrated in either idealized smallscale numerical or sandbox experiments (e.g., Simmons et al., 2001;Post and Kooi, 2003;Post and Simmons, 2009;Jakovovic et al., 2011) or in numerical studies describing present-day salinization patterns in real-world aquifers (e.g., Oude Essink et al., 2010;Nocchi and Salleolini, 2013;Cobaner et al., 2012).However, reports on the long-term effects of these processes in real-world aquifers over the timescales considered here have remained scarce.Post and Kooi (2003) investigated the ability of free convective infiltration to salinize high-permeable aquifers, and reported possible infiltration velocities of several meters per year.With infiltrating saltwater reaching 150 m depth within decades where aquitards are absent, our results confirm these findings.In addition, they show the importance of these "infiltration hotspots" for the salinization of regions underlying lowpermeable strata (Simmons et al., 2001).Our results further signify that aquitards at greater depths can effectively impede infiltration in overlying strata, with effects still visible in the present-day salinity distribution.However, we did not find aquitards that effectively resist salinization or freshening (Groen et al., 2000): the timescales considered are evidently long enough, and hydraulic gradients high enough in the considered setting to eventually lead to salinization or freshening.Mixing zones between salt-and freshwater have been reported to vary widely in geologic settings across the world (Werner et al., 2013), resulting from diffusive processes and kinetic mass transfer (Lu et al., 2009).Measurements indicate that mixing zones even vary widely within our modeled transect, from a narrow zone of some meters around the BGH lens beneath the coastal dunes, to a wide vertical mixing zone of several tens of meters inland.We attribute this variation to the difference between the relatively steady evolution of the BGH lens, resulting in a transition width controlled by molecular diffusion and pore-scale transversal dispersion (Paster and Dagan, 2007), versus the highly transient evolution history of groundwater salinity on the landward side.Wide mixing zones associated with vertical seawater intrusion result from the slow mixing of small-scale (meters) seawater fingers with the intermediary freshwater after the fingers reach the aquifer bottom and flow effectively stops (Kooi et al., 2000).Model-wide dispersivity values could not satisfactorily simulate both these extremes.A dual-domain approach (Lu et al., 2009) could perhaps yield better results.
The present-day salinity distribution in the coastal zone of the Netherlands has been widely recognized to result from free-convective infiltration during Holocene transgressions, widespread peat development and subsequent degradation, and the increasing anthropogenic influence in the more recent past (Post, 2004;Stuyfzand, 1993).Our modeling results clearly support this evolution history, but provide a more detailed overview of the processes involved, signifying the role of aquitards, the ice-pushed ridge flow system, and premodel Maassluis water.The role of (connate) salt present in the Maassluis formation in explaining the present-day salinity distribution has been the subject of some controversy over the past decades, with some authors attributing the presence of brackish groundwater in Pleistocene aquifers to upwards diffusion of salts from this formation (Meinardi, 1991;Volker, 1961).Our results indicate that (1) some salt must have been initially (i.e., at model start, 6500 BC) present in this formation to explain the present-day occurrence of this water type, (2) this salt was rapidly mobilized during the marine transgression in the Holocene, and (3) pre-model salt still contributes with a significant percentage to the presentday salinity distribution.The evolution of salinity in the Maassluis formation throughout the trans-and regression periods during the Quaternary clearly deserves further study.
Regional-scale numerical modeling studies to understand present-day or predict future variable density flow in realworld aquifers are faced with the problem of applying an initial salinity distribution, inherently unknown as measurements are sparse and salinity varies on short-length scales.An often-applied approach is the assumption of steady state, conditioned on present-day boundary conditions (e.g., Souza and Voss, 1987;Vandenbohede and Lebbe, 2002).Such an approach is clearly not warranted for coastal aquifers as the one considered here, as the present-day salinity distribution is almost entirely determined by the significant paleogeographical changes occurring throughout the Holocene.Moreover, neither the salinity distribution nor, consequently, the head distribution was, throughout the modeled period, ever in steady state.A second approach is the construction of the initial salinity distribution based on available measurements, generally by three-dimensional interpolation of available point measurements of salinity (Van der Meij and Minnema, 1999;Oude Essink et al., 2010).This approach necessarily disregards the many small-scale salinity variations not represented by measurements, resulting in an initial salinity distribution not in equilibrium with applied boundary conditions and, consequently, in difficult-to-detect model artifacts.Recent advances have been made in the application of airborne geophysical techniques (AEM) to map groundwater salinity variations over larger scales as input to numerical models (Faneca Sànchez et al., 2012).While this technique is very promising, adequate separation of the contribution of lithology and salinity in the acquired signal requires elaborate ground-truthing (Gunnink et al., 2012), but the resolution of AEM is too coarse to delineate small-scale features and its accuracy decreases with depth.The primarily conceptual paleo-geographical modeling approach presented in this paper cannot yet claim to provide an alternative to the above approaches.Ultimately however, the incorporation of the presented approach within a rigorous uncertainty framework, calibrated to an increasing amount of present-day salinity data supplemented with airborne techniques, may prove successful in adequately describing present-day salinity distributions.

Conclusions
We successfully modeled the effect of paleo-geographical changes throughout the Holocene on the intrusion and redistribution of salts in a representative coastal aquifer.This approach refined our current understanding of the evolution of the salinity distribution in the coastal region of the Netherlands, and yielded insights into the long-term, real-world effects of processes previously investigated in idealized experimental settings.Not once reaching steady state throughout the Holocene, our results attest to the long-term dynamics of salinity in coastal aquifers.The implications of our results extend beyond understanding the present-day distribution of salinity, as the proven complex history of coastal groundwater holds important clues for understanding the distribution of other societally relevant groundwater constituents like nutrients (Van Rees Vellinga et al., 1981;Stuyfzand, 1993) and arsenic (Harvey et al., 2006;Michael and Voss, 2009).
The Supplement related to this article is available online at doi:10.5194/hess-18-3891-2014-supplement.

Figure 1 .
Figure 1.Location of studied transect (A-B), elevation and main topographical features (a), and a lithological cross section along the transect (A'-B) (b); dashed line in (b) demarcates Pleistocene and Holocene deposits.

Figure 2 .
Figure 2. Overview of Holocene paleo-geographical development (a-f) and sea level rise (g, adapted from Van de Plassche, 1982).Red dots and letters in (g) refer to the corresponding paleogeographical map (a-f).For reference, note that the extent of the paleo-geographical maps equals the extent of Fig.1a.
(a) the redox level, as deduced from the concentrations of O 2 , NO 3 , SO 4 , Fe, Mn and CH 4 ; (b) the calcite saturation index; (c) a pollution index (POLIN) based on six equally weighted quality aspects; and (d) the Base EXchange index (BEX), defined as the sum of the cations Na, K, and Mg (in meq L −1 ), corrected for a contribution of sea salt.
Sea level rise, linearly from 22 to 8 m b.s.l.(10 stress periods).Maximum transgression extent reached.Tidal area develops over Pleistocene surface, "basal" peat deposits left mostly intact.Surface drainage.4500 BC-3300 BC Sea level rise, linearly from 8 to 5 m b.s.l.(10 stress periods).Open system with strong marine influence.Deposition of marine clay and sand.Peat extent expands.3300 BC-2100 BC Sea level at 3.5 m b.s.l.Closed system, freshening of hinterland.Peat development behind barriers, peat elevation 3 m b.s.l.2100 BC-700 BC Sea level at 2 m b.s.l.Peat development accelerates, peat domes rise to 1 m m.s.l.700 BC-AD 500 Sea level at 1 m b.s.l.Peat elevation 1.5 m m.s.l.River Vecht system develops (0.7 m b.s.l.).AD 500-AD 1500 Sea level at 1 m b.s.l.Maximal peat elevation: 2 m m.s.l."young dunes" develop, coastal dunes rise to 12 m m.s.l.AD 1500-AD 1850 Sea level at 0.3 m b.s.l.Rapid degradation of peat due to peat extraction and anthropogenic drainage (0 m m.s.l.).Freshwater lakes develop.Water level River Vecht 0.05 m b.s.l.AD 1850-AD 1900 Sea level at 0.1 m b.s.l.Reclamation of Haarlemmermeer (AD 1852) and Horstermeer (AD 1882).Anthropogenic drainage through canals and ditches.AD 1900-AD 1950 Sea level at 0.05 b.s.l.Subsurface drains introduced.Groundwater (over-) abstraction in coastal dunes.AD 1950-AD 2010 Present-day situation, sea level at 0 m.s.l.Groundwater abstraction in ice-pushed ridge, groundwater abstraction in coastal dunes decreased.

Figure 4 .
Figure 4. Comparison of measured versus modeled heads.Locations (Fig. 1) were selected within a trapezoidal buffer (0.5 km at the surface to 5 km at 300 m depth) around and projected orthogonally onto the modeled transect.Measurement values are the average of time series of head measurements from AD 1990 onwards containing at least 25 measurements.

Figure 5 .
Figure 5. Chloride measurements versus modeled chloride concentration at AD 2010.Measurements were selected within a trapezoidal buffer (2 km (b) at the surface to 5 km at 300 m depth) around and projected orthogonally onto the modeled transect.

Figure 7 .
Figure 7. Modeled evolution of groundwater chloride concentration (a-g).White lines are contours of the stream function, contour intervals are equal for all time slices.Except for (a) (starting concentration), transect times correspond to paleo-geographical maps in Fig. 2.

Figure 8 .
Figure 8. Contribution of transgression, sea, and Maassluis to total salt present (SP).Only the model domain east from x coordinate 95 km is considered.Vertical dashed lines denote time-slice transitions, discontinuities at transitions result from changing numbers of active model cells.Legend percentages are percentages at model end.

Figure 9 .
Figure 9. Modeled evolution of groundwater origin (a-g).White lines are contours of the stream function, contour intervals are equal for all time slices.Except for (a) (starting situation), transect times correspond to paleo-geographical maps in Fig. 2.

Table 1 .
Description of modeled origin tracers and relation to hydrosomes.

Table 2 .
Description of model time slices.

Table 3 .
References for paleo-model implementation.