Introduction

In the southwestern delta of the Netherlands large areas lie at or below mean sea level (MSL). Natural dunes and man-made dikes prevent these areas from flooding, and the surface-water level in the ditches is artificially controlled using weirs, dams, and mechanical pumps. In this landscape, sandy tidal creek deposits lie up to 2 m higher than the surroundings. These geomorphological features are known as ‘creek ridges’. Below the creek ridges, 5–30-m-thick freshwater lenses are present. These lenses are an important source for irrigation, because the surrounding saline groundwater is predominantly saline and because in summer most of the surface waters are brackish to saline and there is no supply of fresh surface water from other areas (De Louw et al. 2011). The amount and rate of groundwater extraction from freshwater lenses is regulated by the local water board to prevent excessive saltwater upconing and drawdown of the water table (Province of Zeeland 2002). Currently, the irrigation demand already exceeds the permissible groundwater extraction, which leads to crop damage. As the anticipated sea-level rise, climate change, and land subsidence will further jeopardize fresh groundwater resources in the future (Oude Essink et al. 2010; Pauw et al. 2012), measures are necessary in this delta to assure a sustainable freshwater supply.

In this paper a new system to increase the thickness of freshwater lenses below creek ridges is presented. This controlled artificial recharge and drainage (CARD) system aims at increasing the freshwater lens by raising the water table during winter, when there is a precipitation excess. In view of the classical Badon Ghijben-Herzberg relationship, a higher water table leads to a thicker freshwater lens. Evidently, from thicker freshwater lenses, more fresh groundwater can be extracted (Dagan and Bear 1968; Volker et al. 1985; White and Falkland 2009). The CARD system differs from other artificial recharge techniques (Bouwer 2002) that are aimed at increasing freshwater supply in coastal areas such as spreading basins (Stuyfzand 1993; Tredoux et al. 2003; Pyne 2005) and aquifer storage and recovery (ASR) using wells (Miotlinski et al. 2013; Zuurbier et al. 2014a, 2014b), in that agricultural drain tiles are used to raise the water table.

The CARD system was investigated in the field on a creek ridge on Walcheren, a peninsula in the southwestern part of the Netherlands (Fig. 1). The development of the freshwater lens below the creek ridge was monitored using various types of field measurements. The CARD system at the field site was further investigated using numerical models. This paper describes the hydrological aspects of this case study; hydrochemical and economic aspects of the CARD system were not considered. The important hydrogeological factors for applying the CARD system at other locations are qualitatively discussed.

Fig. 1
figure 1

Location of the study area (red box) on Walcheren, a peninsula in the southwestern part of the Netherlands. Walcheren is surrounded by three surface waters: the North Sea, the Western Scheldt estuary, and Lake Veere, which are indicated in blue

Material and methods

Theory and definitions

In this paper, the classification of the salinity of surface water and groundwater into ‘fresh’, ‘brackish’, or ‘saline’ is based on the electrical conductivity of the water referenced to a temperature of 20 °C (ECw20). The water is fresh if the ECw20 ≤ 2.0 mS/cm, the water is brackish if the ECw20 is between 2.0 and 5.0 mS/cm, and the water is saline water if the ECw20 ≥ 5.0 mS/cm. This classification is based on guidelines that freshwater users in the area use for irrigation; freshwater can be used safely for irrigation, brackish water is only used occasionally in times of severe water scarcity, and saline water is unfit for irrigation.

The thickness of the freshwater lens is defined as the depth at which the ECw20 or ECbulk (the electrical conductivity of the water and porous medium) is equal to 50 % of the difference between their maximum and minimum value below the creek ridge. This value is referred to as D mix (modified from De Louw et al. 2011). D mix approximately corresponds to the centre of the mixing zone of the freshwater lens. The reason for using D mix rather than depth of the ‘fresh–brackish’ interface for the thickness of the freshwater lens is that, from the different geophysical measurement techniques that were used, it was difficult to determine the fresh–brackish interface. The main reason for this was that data about the formation factor for the conversion of the ECbulk to ECw20 using Archie’s law (Archie 1942) were not collected.

The CARD system aims at increasing the freshwater lens raising the water table below the creek. A classical principle of the relationship between the water table and the thickness of the freshwater lens is the Badon Ghijben-Herzberg principle (Du Commun 1828; Drabbe and Badon Ghijben 1888; Herzberg 1901):

$$ h=\delta H $$
(1)

h is the elevation of the water table and H is the thickness of the lens (both H and h are relative to a reference level). δ is the dimensionless density difference:

$$ \delta =\frac{\rho_{\mathrm{s}}-{\rho}_{\mathrm{f}}}{\rho_{\mathrm{f}}} $$
(2)

in which ρ s is the density of the saline water and ρ f is the density of the freshwater. The two most important assumptions of the Badon Ghijben-Herzberg principle are (1) no mixing zone is present between the fresh and saline groundwater and (2) the vertical pressure distribution is hydrostatic. In reality, these conditions are never exactly met. Amongst others, water-table fluctuations will never result in immediate changes of the thickness of the freshwater lens according to Eq. (1), because hydraulic heads in the subsurface react much faster to changing hydraulic conditions than the groundwater salinity distribution does (Oude Essink 1996; Post 2003). Nevertheless, in the long term, an increase of the freshwater lens can be expected when the time-averaged water table is raised (Eeman et al. 2012), which is the basic idea behind the CARD system.

Study area and design of the CARD system

The study area is located on Walcheren (Fig. 1). In the low-lying areas of Walcheren, the upper hydrogeological unit is a semi-confining layer up to 10 m thick, which consists of Holocene peat, clay, and silt deposits. Below this unit, sandy deposits from the Quarternary and Neogene periods form a 10–80-m-thick semi-confined aquifer. In the semi-confined aquifer, low-permeability layers are present, but they do not form a coherent confining unit. The semi-confined aquifer is bounded below by Oligocene marine clay, which is considered as the hydrogeological base (TNO 1997; Goes et al. 2009; De Louw et al. 2011; Stafleu et al. 2011).

At many locations on Walcheren, Holocene sandy tidal creek deposits fully incise the semi-confining layer and partly incise the underlying aquifer up to a depth of −35 m MSL. The upper 2 m of the sandy tidal creek deposits show a fining up sequence from medium sand to clayey silt, which is related to the silting-up of the tidal creeks (van Meerten 1986). These tidal creeks were active between 700 BC and AD 750. Around AD 1,200, natural and man-made drainage led to compaction of peat and clay deposits of the tidal flat, which initiated an inversion of the landscape (Fig. 2; Ervynck et al. 1999; Vandenbohede and Lebbe 2002); it is in this way that the creek ridges developed.

Fig. 2
figure 2

a Concept of the landscape before land reclamation. b Concept of the landscape after the silting up of the tidal creek and drainage of the peat and clay deposits in the semi-confining layer, which led to an inversion of the landscape

Within the study area, the focus is on the fields of two farmers (Figs. 1, 3, and 4). The farmers use the land for crop rotation of cichorium (chicory), potatoes, cauliflower, and fennel. In the period 1971–2010, the average precipitation excess in autumn and winter was 261 mm, whereas in spring and summer the average potential evapotranspiration was 135 mm higher than the precipitation (KNMI 2015, weather station Vlissingen, ~11 km from the study site). In the north and east, the study area is bordered by a 10-m-wide ditch (Fig. 4). The ECw20 of this ditch varies between 3 (winter) and 9 (summer) mS/cm due to seepage of saline groundwater, making the salinity of this water too high for irrigation. In summer, the water level in the ditch is maintained at −1.0 m NAP (the Dutch ordnance datum, approximately equal to MSL) and in winter at −1.2 m NAP. In between the fields of the two farmers, a 1-m-wide ditch is present.

Fig. 3
figure 3

a Concept of the CARD system in vertical cross section. b Overview of the CARD system in plan view (part of the study area, see Figs. 1 and 4). Close and parallel to the cross section (a), a ditch is present that separates the fields of two farmers. As this ditch was draining fresh groundwater, a dam was placed in conjunction with the installation of the CARD system to prevent undesirable drainage of fresh groundwater

Fig. 4
figure 4

a The study area (also shown in Fig. 1) and the horizontal extent of the numerical model in plan view, including the surface elevation. b Overview of the horizontal extent of the CARD system, the measurement locations (MLs) and continuous vertical electrical sound (CVES) transects. Table 1 gives an overview of the different types of field measurements per ML. c Cross-sectional transect of the horizontal hydraulic conductivity (K h) in the numerical model. This transect is indicated by the black dashed line in a. The purple solid line indicates the cross-sectional transect of the numerical model results shown in Fig. 11. CVES transect B–B′ is partly located along this transect, as well as partly along the cross-sectional transect of the horizontal hydraulic conductivity (c)

In April 2013, the CARD system was installed on the two fields and a dam in the small ditch between the fields was placed to prevent drainage of the fresh groundwater (Fig. 3). The drain tiles of the CARD system are composed of plastic perforated tubes of 60 mm with an outer envelope of polypropylene (PP450). The depth of the drain tiles varies between 1.1 m (higher parts of the creek ridge) and 0.8 m (lower parts of the creek ridge) below ground level (BGL). At the higher parts of the creek ridge the distance between the drain tiles is 7 m. Artificial recharge using fresh surface water takes place here. The fresh surface water is collected downstream of a weir just north of the CARD system (Fig. 3). The weir controls the discharge of fresh surface water into a brackish/saline ditch. Discharge mainly occurs in late autumn, winter, and early spring, as a result of the precipitation excess, so that, in this period, water is available for artificial recharge. The collected freshwater is passed through a 0.2-mm filter (Meeuwse Handelsonderneming), is transported over the brackish/saline ditch, and is then pumped into the CARD system at an intake point (Fig. 3). When the water table is higher than the highest desired drainage level (about 1 m BGL at the highest parts of the creek ridge), the pump automatically turns off. During the other period of the year, artificial recharge is generally not possible, as the weir prevents the discharge of freshwater into the brackish/saline ditch.

At the lower parts of the creek ridge, artificial recharge does not take place, but the drainage level is kept as high as possible (between 0.6 and 0.8 m BGL) to reduce the groundwater flow from the higher parts of the creek ridge. A water table higher than 0.6 m BGL at the lower parts of the creek ridge is undesired as it leads to crop damage. In Fig. 3 it is shown how the drainage level is controlled. Individual drain tiles at 6-m spacing are oriented perpendicular to the plane of Fig. 3. These drain tiles are connected to a collector drain. The water in this collector drain discharges into a collector well, where the drainage level can be controlled at the desired level. Additional illustrative pictures of the CARD system are given in the electronic supplementary material (ESM) of this paper.

Field measurements

In Fig. 4b and Table 1 an overview is shown of the measurement locations (ML) and the different field measurements that were used to study the effect of the CARD system on the thickness of the freshwater lens. Prior to the installation of the CARD system, field data of the lithology, ECbulk, and ECw20 was collected for the design of the system and for the calibration of the numerical model (Fig. 4a). Cone penetration tests (CPTs) were used to determine the vertical distribution of the lithology and ECbulk at five MLs. For the lithology, a friction ratio (R f) > 2 was used to indicate deposits with a low permeability such as clay or peat. The friction ratio is the ratio of the pressure on the tip of the cone and on the sleeve (i.e., the side) of the cone (Lunne et al. 1997). These pressures were recorded together with ECbulk at a 2-cm interval as the cone was pushed downwards. As the ECbulk is also influenced by lithology (e.g., clay), the influence of the lithology was qualitatively considered using R f.

Table 1 Overview of the different types of field measurements and their measurement location (ML)

Continuous vertical electrical sounding (CVES) measurements were carried out at three transects to estimate the distribution of ECbulk in a vertical cross section (Fig. 4b). The ABEM Lund Imaging System, comprising the ABEM SAS 4000 terrameter and the ES10-64 electrode selector, was used for this purpose. The Wenner method was used with an electrode distance of 2 m to measure the ‘apparent’ ECbulk distribution (Reynolds 1997). The RES2DINV inversion program (Loke 2006) was used to obtain a model of ECbulk.

Monitoring

After the installation of the CARD system, the water table, ECbulk and ECw20 were monitored at MLs 1, 2, 4, 5, and 6 (Fig. 4b), for the period May 1, 2013–May 1, 2014. Measuring location 6 served as a reference measurement, as the CARD system was not installed here.

At ML 2, a vertical electrical resistivity monitoring system (subsurface monitoring device; SMD) was installed. The borehole for this system was made using pulse drilling with a diameter of about 15 cm. The SMD was used to measure the vertical distribution of ECbulk at a daily interval, using multiple electrodes located at various depths. The electrode spacing was 0.1 m from −6.5 to −10.3 m NAP. Above and below this section the electrode distance gradually increased to 0.7 m. The Wenner configuration was used to measure ECbulk.

At MLs 1 and 4, six 2.5-cm-diameter piezometers with a 40-cm-long screen were installed at depths where the mixing zone was present before the installation of the CARD system. The borehole was made using rotary flush drilling, with a diameter of about 20 cm. Bentonite plugs in between the screen were used to prevent short-circuiting of the groundwater (Stuyfzand 1993). About two times the volume of the piezometer was pumped to determine the ECw20, which was sufficient to obtain a stable value. During the period September 2013–May 2014, the ECw20 was sampled approximately bi-monthly.

In December 2013 and February 2014, electromagnetic (EM) borehole logging using a thin probe (SLIMFLEX) was conducted at MLs 1, 5, and 6. In the probe, an alternating electrical current is generated in a ‘transmitting’ coil to produce a primary EM field. This primary EM field produces electrical currents in the subsurface which lead to a secondary EM field. Both fields are measured in another coil at 50-cm distance from the transmitting coil. From the difference in amplitude and phase lag of the EM fields, ECbulk can be determined (McNeill 1980). The coils are situated such that the influence of the water in the borehole is negligible.

At MLs 1, 2, 4, and 6, groundwater levels were recorded at a 10-min interval using pressure transducers installed in piezometers. The piezometers had a 40-cm screen at the bottom, and the screens were located at 4 m from the surface. The borehole for the piezometer was made using hand pulse drilling equipment. The pressure transducers were read out, compared with manual measurements and re-installed approximately bi-monthly. The maximum difference between the manual measurement and the recorded pressure was 4 cm.

Numerical model

Set-up of the initialization model

SEAWAT version 4 (Langevin et al. 2007) was used to simulate variable density groundwater flow and solute (salt) transport in the study area (Fig. 4a). The horizontal extent of the model is 1,600 × 1,600 m and discretized using cells of 10 × 10 m. The vertical extent of the model is from +1.5 m NAP to −87.0 m NAP. Up to a depth of −25 m NAP, the model layers are 0.5 m thick. The hydraulic conductivity distribution in this section was based on the hydrogeological model GEOTOP (Stafleu et al. 2011). Below −25 m NAP, the thickness of the model layers increases to up to 5 m. The hydraulic conductivity distribution in this section was based on the regional hydrogeological model REGIS II (REGIS 2005). The vertical anisotropy ratio—the horizontal hydraulic conductivity (K h) divided by the vertical hydraulic conductivity (K v)— in the numerical model is 1.4. In Fig. 4c, a cross section of K h in the numerical groundwater model is shown, which indicates the incision of the semi-confining layer by the more permeable deposits below the creek ridge. This is in line with the concepts shown in Figs. 2 and 3. The effective porosity (0.3), confined storage coefficient (1E-5 1/m), specific yield (0.15), longitudinal dispersivity (0.1 m), transverse dispersivity (0.01 m), and molecular diffusion coefficient (8.64E-5 m2/d) were assumed to be constant in the model. These parameters were estimated from general literature (Domenico and Schwarz 1997) and from groundwater flow studies in the coastal plain of the Netherlands (Oude Essink et al. 2010; De Louw et al. 2011; Zuurbier et al. 2014b) and Belgium (Lebbe 1999; Vandenbohede and Lebbe 2006).

First, an initialization model was constructed in which the development of a freshwater lens below the creek ridge was simulated, starting from an initially saline domain. The solute (salt) concentration was represented by ECw20. The maximum value of ECw20 corresponds to the maximum measured ECw20 (40.6 mS/cm at a depth of −18 m NAP), which is also the initial concentration of the initialization model. The density of the water with ECw20 = 40.6 mS/cm was calculated at 1,022.3 kg/m3, following the method of Post (2011). The minimum concentration of the groundwater is 0 and corresponds with a density of 1,000 kg/m3. The water density (ρ w) varies linearly with ECw20:

$$ {\rho}_{\mathrm{w}}=1,000+\beta\ \mathrm{E}{\mathrm{C}}_{\mathrm{w}20} $$
(3)

with β = (1,022.3 − 1,000) kg/m3/40.6 mS/cm.

The PCG package was used for solving the variable density groundwater flow equation. The hydraulic head closure criterion was 10−5 m and the flux closure criterion was 0.1 m3. The advection term of the solute transport equation was solved using the MOC solver, with a minimum number of nine particles per cell. A Courant number of 1.0 was used. The sink and source and the hydrodynamic dispersion terms of the solute transport equation were solved using the GCG solver, with a relative concentration closure criterion of 10−6.

Stress periods of 6 months were used, representing winter and summer conditions. The recharge rate in winter (1.5 mm/d) and the evapotranspiration rate in summer (−0.02 mm/d) were based on average data (period 1971–2000) from a nearby weather station—station Vlissingen, ~11 km from the study site (KNMI 2015). For the evapotranspiration, crop factors based on the crops on the fields of the farmers were taken from Hooghart and Lablans (1988), to correct for the reference Makkink evapotranspiration (Makkink 1957). At the vertical sides of the model, a head-dependent flux boundary condition was applied with the GHB package, using hydraulic heads and concentrations taken from a larger-scale model (Van Baaren et al. 2011). The conductance value was set to 1 m2/d. Testing showed that an order of magnitude higher or lower value for the GHB conductance did not significantly change the model results at the study site. This was in line with the expectations, as the boundaries were located relatively far from the CARD system. The bottom of the model is a no-flow boundary condition as it corresponds with the hydrogeological base. Ditches were simulated in the model using the RIVER package. Tile drainage was modelled using the DRAIN package. The conductance for the cells where the RIVER or DRAIN packages were used was based on the methods of De Lange (1999) and Pauw et al. (2014). More information on the input parameters of the DRAIN and RIVER conductance is given in the ESM.

Calibration of the initialization model and sensitivity analysis

The initialization model was ended at the moment the thickness of the freshwater lens below the creek ridge (D mix) reached a virtually stable value. The simulated end-thickness of the freshwater lens was compared with the thickness of the freshwater lens measured by the CPTs. Using the parameters that were described in the previous section, the thickness of the freshwater lens was overestimated in the model. Therefore, the model was calibrated to obtain a better match with the field observations. For the calibration, the groundwater recharge rate in winter was reduced, as this was considered the most uncertain parameter. Overestimation of the groundwater recharge is likely, as ponding and subsequent surface runoff frequently occurs at the field site due to the low permeability of the deposits close to the surface of the creek ridge, due to the fining-up sequence. As a result of ponding and surface runoff, the groundwater recharge is lower than the calculated value as described in the previous section (Healy 2010; Voortman 2010).

In addition to the model calibration, the sensitivity of other model parameters was investigated. However, due to the long runtime of the initialization model (approximately 1 week on a modern desktop PC), an extensive sensitivity analysis was not conducted. The effect of a higher evapotranspiration rate in summer was investigated by increasing it from 0.2 to 0.4 mm/d. This value is comparable with the anticipated increase of the evapotranspiration rate in summer to climate change (De Louw 2013; KNMI 2015). In addition, the hydraulic conductivity in the upper 15 m below the higher parts of the creek was multiplied by a factor of 2, for investigating the effect of the hydraulic conductivity in decreasing the thickness of the freshwater lens. Furthermore, a two times higher and two times lower value of the conductance of the RIVER and DRAIN cells was investigated, as these values were based on various assumptions (see ESM).

Scenario simulations to investigate the effect of the CARD system

The simulated hydraulic heads and concentrations at the end of the initialization simulation with the calibrated groundwater recharge in winter were used as initial condition for three scenario simulations where the effect of the CARD system on the thickness of the freshwater lens was simulated. In addition, a ‘reference’ simulation where the CARD system was not simulated was used for comparison with the three scenario simulations. In all four simulations, the yearly period May 1, 2013–May 1, 2014 was simulated using weekly stress periods. The recharge and evapotranspiration rates were based on data from a nearby weather station—station Vlissingen, ~ 11 km from the study site (KNMI 2015). The yearly period was repeated to cover a period of 10 years to study a long-term effect of the CARD system.

In scenario A the influence of the CARD system on the thickness of the freshwater lens was simulated. The RIVER package was used to represent infiltration and drainage by the drain tiles of the CARD system. The hydraulic head of each RIVER cell was based on the weekly average hydraulic heads that were measured at MLs 1, 2, and 4. For verification of this approach, the simulated artificial recharge was compared with the measured artificial recharge at the water collection point of the CARD system. Further information on the simulation of the CARD system is available as ESM.

The effect of a lower availability of fresh surface water for artificial recharge was investigated in scenario B by simulating a 0.1-m-lower water table during the winter period (i.e., November 1, 2013–May 1, 2014) compared to scenario A. In scenario C a 0.1-m-higher water table compared to scenario A was simulated to investigate the effect of a higher water table.

Results

Field measurements

Geology

In Fig. 5 and Table 2 the results of the CPTs at MLs 1–5 are shown. The R f in the CPTs at MLs 1, 2, and 5, at the higher parts of the creek ridge, is lower than 2 in the upper 1–2 m, indicating the fining-up sequence that is typically found below the higher parts of creek ridges. Below the fining-up sequence, the R f decreases, which indicates sandy tidal deposits. In the CPTs at MLs 3 (towards the lower part of the creek ridge) and 4 (at the lower part of the creek ridge), the R f is lower than 2 over a larger part of the upper section due to the presence of silt, clay and peat deposits, which form the semi-confining layer. Note that the semi-confining layer is also indicated by the GEOTOP model (Fig. 4).

Fig. 5
figure 5

ECbulk and R f measured by the cone penetration tests (CPTs) at measurement locations (MLs) 1–5 (Fig. 4)

Table 2 Estimated D mix from the CPTs and simulated D mix in the numerical model at MLs 1–5

The CPT at ML 1 indicates low-permeability deposits between −12 and −12.5 m NAP. This is also observed in the CPTs at MLs 2, 3, and 4, but at different depths between −8.0 and −15.0 m NAP. In the CPT at ML 5, low-permeability deposits at depth seem to be absent. These sections of low-permeability deposits are interpreted as bank deposits of the former creek system. This type of heterogeneity is not present in the GEOTOP model. In the ESM the R f is shown in combination with the vertical hydraulic conductivity at MLs 1–5.

ECbulk prior to the installation of the CARD system

In Fig. 6 the inversion results of the CVES transects A–A′, B–B′, and C–C′ are shown. The ECbulk values ≤ 2 mS/cm indicate the freshwater lens in the sandy deposits below the creek ridge. Further down, the ECbulk increases relatively fast, which indicates the transition towards the ambient saline groundwater. CVES transect C–C′ indicates a thicker freshwater lens compared to CVES transects A–A′ and B–B′. This is in line with the difference in D mix between the CPT at ML 5 and the CPTs at MLs 1 and 2 (Fig. 5; Table 2), and shows that the freshwater lens increases in thickness towards the south. CVES transect B–B′ indicates a 0.5–1.0 m lower D mix than the CPTs at MLs 1–4, which can be attributed to the lower vertical resolution in the CVES compared to the CPTs (Reynolds 1997).

Fig. 6
figure 6

Inversion results (ECbulk) of a CVES transect A–A′, b CVES transect B–B′, and c CVES transect C–C′

Monitoring of ECbulk and ECw20

In Fig. 7, the daily precipitation and the hydraulic heads relative to NAP at ML 1 are shown. From approximately mid-October 2013 until mid-April 2014, the water table could be maintained at a level between 0.2 and −0.2 m NAP. After mid-April 2014, the water table gradually dropped as the artificial recharge stopped. This was because the discharge of fresh surface water over the weir at the water collection point (Fig. 3) ceased.

Fig. 7
figure 7

a Daily precipitation (P, in mm/d) and b hydraulic head relative to NAP (h, in m) at ML 1. The surface elevation is 1.08 m NAP. The month indications on the horizontal axis indicate the first day of that month

The results of the ECw20 sampling at ML 1 are shown in Fig. 8. At September 20, 2013 (prior to the start of the artificial recharge), D mix (20.3 mS/cm) was located at about −12 m NAP. This agrees well with D mix derived from CPT at ML 1 (Fig. 5; Table 2). Between September 20, 2013 and March 28, 2014, the freshwater lens increased about 2 m in thickness (Fig. 8). The sampling results of May 1, 2014 indicate a slight decrease (about 0.05 m) of the freshwater lens between March 28, 2014 and May 1, 2014. At ML 4, the ECw20 at all depths were higher than 30 mS/cm and did not change significantly during the monitoring period, indicating that no significant freshening occurred at this location.

Fig. 8
figure 8

ECw20 (mS/cm) at six depths during five monitoring rounds at ML 1

In Fig. 9, the development of ECbulk, measured by the SMD at ML 2 during the monitoring period, is shown. The decrease of ECbulk indicates that freshening took place in between −8.0 and −12.0 m NAP. A remarkable feature of Fig. 9 is the relative constant ECbulk around z = −9.0 and around −10.5 m NAP. This is attributed to a high clay content in the deposits at these depths. Therefore, the permeability at this depth is expected to be low. The observed freshening can be explained by two processes. The first, which is considered the most likely processes, is that the freshening below −9.0 m NAP does not take place from above (i.e., through the clay-rich deposits), but by preferential flow through sections where the hydraulic conductivity is higher. This hypothesis is strengthened by the CPTs which indicated that the low-permeability sections at the higher parts of the creek ridge do not constitute a continuous unit. The second explanation is that the clay-rich section is still permeable enough to transmit fresh groundwater from above.

Fig. 9
figure 9

Development of the ECbulk in the period 1 May 2013–1 May 2014 relative to NAP, obtained from the SMD at ML 2. The month indications on the horizontal axis indicate the first day of that month. The erratic pattern of ECbulk between −12 and −15 m NAP is attributed to measurement errors and should therefore be ignored

Fig. 10
figure 10

ECbulk relative to depth at MLs 1, 5, and 6, obtained using the SLIMFLEX measurements

The SLIMFLEX measurement at ML 1 showed that D mix lowered by about 0.8 m in the period between December 10, 2013 and February 7, 2014. This agrees well with the sampling results for this ML (Fig. 8). The SLIMFLEX measurements at ML 5 show a lowering of D mix of about 0.5 m. At ML 6 (the reference ML) the SLIMFLEX measurements show no significant changes in ECbulk. This was expected, as the CARD system was not installed here.

Numerical simulations

Calibrated initialization model

The recharge rate in winter was reduced (calibrated) from 1.5 to 1.2 mm/d to obtain a reasonable match between the simulated D mix and the measured D mix by the CPTs, at MLs 1–5 (Table 2). At ML 1, the simulated D mix is about 1 m deeper than the measured D mix. This is attributed to the presence of low-permeability deposits at ML 1, indicated by the increase of the friction ratio R f in the CPT, and the absence of these low-permeability deposits in the numerical model. At MLs 2 and 3, the CPTs indicate more sections of low-permeability deposits than at ML 1. These sections are also not present in the numerical model; therefore, at these MLs, the simulated D mix is even deeper than the measured D mix, compared to ML 1. At ML 4, the difference between the simulated and measured D mix is comparable to the difference at ML 1. At this ML, the CPT indicated low-permeability deposits up to about −8.0 m NAP, i.e., the semi-confining layer. In the numerical model the semi-confining layer is also present. At ML 5 the simulated D mix is within the estimated uncertainty (0.5 m) of the measured D mix. At this ML, the CPT indicated no low-permeability deposits (Fig. 5). In the ESM, further information on the comparison between the CPTs and the numerical simulations is given.

In Fig. 11a, D mix is shown in plan view at the end of the calibrated initialization model. A total of 230 years of simulation time were needed to reach a steady-state position of D mix. D mix shows a clear relation with the elevation map (Fig. 4). At the creek ridge, the drainage level and the water table are higher than the surrounding lower lying areas; therefore, the freshwater lens is also thicker here. This is also indicated in Fig. 11b, where the simulated ECw20 at the end of the initialization model is shown in vertical cross section. In the lower-elevated areas with a low drainage level and low-permeability deposits at the surface, the lenses are thinner, presumably due to upward flow of saline groundwater (De Louw et al. 2011). The low water level (−1.2 m NAP) of the 10-m-wide ditch that dissects the creek ridge north of the CARD system has a negative impact on the thickness of the freshwater lens, as D mix below the ditch is found to be shallower than further away from the ditch.

Fig. 11
figure 11

a D mix in the initialization simulation in plan view, using the calibrated groundwater recharge. The black solid lines indicate ditches. The black dashed line indicates the location of the vertical cross section of the hydraulic conductivity distribution in the model (Fig. 4). The purple solid line indicates the cross section of the simulated ECw20, shown in b. The white dashed square shows the focus area of the results of the scenario models (Fig. 13 and ESM)

Sensitivity analysis

The sensitivity of the different parameters was investigated using the simulated ECw20 and D mix at MLs 1 and 5 (Fig. 12), as at these locations measurements of D mix were also available. The most important observation from Fig. 12 is that the variations in groundwater recharge and hydraulic conductivity have a larger influence on D mix than variations in the DRAIN and RIVER conductances. In the simulation with a higher groundwater recharge rate in winter, the freshwater lens at MLs 1 and 5 was about 1.5 m thicker compared to the simulation with the calibrated groundwater recharge. In the simulation with a higher evapotranspiration rate in summer, the freshwater lens at MLs 1 and 5 was about 1.5 m thinner than the freshwater lens in the simulation with the calibrated groundwater recharge. This suggests that climate change (i.e., a change in the groundwater recharge) can have a significant effect on the thickness of the freshwater lens.

Fig. 12
figure 12

Results of the calibration and sensitivity analysis of the initialization model. The simulated ECw20 as a function of depth (z) is shown for the simulation with the calibrated groundwater recharge (Cal) and the other simulations of the sensitivity analysis. Wet is the simulation where the groundwater recharge in winter is 1.5 mm/d, Dry is the simulation where the evapotranspiration is 0.4 mm/d, High K is the simulation where the hydraulic conductivity (both K h and K v) below the higher parts of the creek ridge is two times higher than in the other simulations, and Low cond and High cond are the simulations where the conductance of cells where the RIVER or DRAIN were used is a factor 2 lower and higher compared to the simulation with the calibrated groundwater recharge, respectively. The black dots indicate the measured D mix by the CPTs at MLs 1 and 5

Effect of the CARD system

In the reference simulation, where the CARD system was not operational, the yearly variation of D mix is 0.1 m at ML 1 (Fig. 13a). After 10 years of simulation, D mix is about 0.5 m lower than at the beginning of the reference simulation at this ML. This is also the maximum observed lowering of D mix over the whole freshwater lens below the creek ridge, after 10 years of simulation. This minor lowering of D mix can be explained by the slightly higher recharge rate that was applied in the reference simulation compared to the recharge rate of the initialization model. In scenario A, where the CARD system using the measured groundwater levels was simulated, the water table was on average 0.5 m higher than in the reference simulation during the winter period, which resulted in a much thicker freshwater lens (i.e., lowering of D mix) after 10 years (Fig. 13).

Fig. 13
figure 13

a Simulated ECw20 contours as a function of depth (z) and time at ML 1 in the reference simulation. b Simulated ECw20 contours as a function of depth and time at ML 1 in scenario A. c Difference in D mixD mix) between the reference simulation and scenario A, after 10 years of simulation time. The black solid lines indicate ditches. The green dots and the white numbers indicate the MLs

At ML 5, the measured 0.5 m lowering of D mix between December 10, 2013 and February 7, 2014 using SLIMFLEX (Fig. 10) agrees well with the simulated lowering of D mix (0.3 m) in the equivalent period during the first year (0.3 m). At ML 1, the difference between the measured and simulated lowering of D mix is larger. The sampling results (Fig. 8) indicate that between September 20, 2013 and March 28, 2014, D mix lowered about 2.0 m, whereas in scenario A, D mix lowered about 1.4 m. This difference is attributed to the initialization model, where the simulated D mix at ML 1 is about 1 m lower than the measured D mix (Table 2). The thicker the freshwater lens, the less the freshwater lens will increase in thickness upon a given increase of the water table, which can be deduced from the Badon Ghijben-Herzberg principle.

In Fig. 13b it is shown that at ML 1 the decrease of the thickness of the freshwater lens due to buoyancy during the summer months is much less than the increase of the thickness of the freshwater lens during the artificial recharge period. This effect also holds for the whole freshwater lens below the CARD system, which indicates a net increase of the thickness of the freshwater lens over a year. This freshening continues for at least 10 years (Fig. 13b). After 10 years, the freshwater lens has not reached its end-thickness, but the freshening rate has considerably decreased. In Fig. 13c it is shown that D mix has increased about 6 m in depth at ML 1 and about 8 m at ML 5 in scenario A after 10 years of simulation.

For brevity, the results of scenarios B and C (analogous to Fig. 13) are given in the ESM, and only the most important results of these scenarios are described here. As expected, in scenario B, in which the water table in the CARD system is 0.1 m lower than in scenario A, the increase of D mix in depth is less than in scenario A. After 10 years, D mix at MLs 1 and 5 has increased about 5 and 6 m compared to the reference simulation, respectively, and the freshwater lens has almost reached a new dynamic equilibrium. As expected, in scenario C, in which the water table in the CARD system is 0.1 m higher than in scenario A, the increase of D mix in depth compared to the reference simulation (8 m at ML 1 and 10 m at ML 5) after 10 years is larger than in scenario A. The differences in increase of the freshwater lens between the scenarios A, B, and C indicate the significant effect of raising the water table below the creek ridge in winter on the thickness of the freshwater lens.

In Fig. 14a–c, the simulated volume of artificial recharge (V AR) and volume of the increase of the freshwater lens below the CARD system (dV lens) are shown per year for the three scenarios. In scenario A, V AR is 25,052 m3 in the first year, which is about equal to the natural recharge that was simulated below the CARD system (27,444 m3) during the artificial recharge periods. Compared to the measured volume of artificial recharge (18,905 m3), the simulated V AR in the first year is overestimated. This overestimation can be explained by the presence of low-permeability deposits that were not present in the numerical model. The low-permeability deposits can influence the downward flow of the groundwater and, therefore, also the volume of artificial recharge.

Fig. 14
figure 14

ac The simulated volumes of the increase of the freshwater lens below the CARD system (dV lens; green) and artificial recharge (V AR; blue) per year, for three scenarios. df Cumulative values of dV lens and V AR for the three scenarios

The yearly V AR in scenario A decreases from 25,052 m3 in the first year to 16,702 m3 in the last year, which indicates the increase of the freshwater lens. In the first 4 years, the yearly dV lens is larger than the yearly V AR, which indicates the contribution of the natural groundwater recharge to the increase of the freshwater lens. In the last 6 years, the yearly dV lens is smaller than the yearly V AR. The value of the storage efficiency S eff, is defined as:

$$ {S}_{\mathbf{eff}}=\frac{\mathrm{d}{V}_{\mathrm{lens}}}{V_{\mathrm{AR}}}100 $$
(4)

In the last 6 years S eff decreases from 90 to 58 %.

As expected, the yearly V AR and dV lens are smaller in scenario B compared to scenario A. However, the yearly S eff is higher; in the first 7 years the yearly dV lens is larger than the yearly V AR and in the last 3 years the yearly S eff has decreased from 94 to 80 %. In scenario C the yearly V AR and dV lens is larger than in scenario A, but the yearly S eff is smaller. Only in the first 2 years is the dV lens larger than V AR. In the last 8 years the yearly S eff has decreased from 94 to 48 %.

In Fig. 14d–f, the cumulative V AR and dV lens are shown for the three scenarios. The total volume of the freshwater lens below the CARD system at the end of the initialization model is 373,088 m3. In scenario A the cumulative V AR (191,186 m3) and dV lens (186,945 m3) after 10 years are about equal. The total simulated natural recharge volume is 274,440 m3; hence, after 10 years, the cumulative dV lens is about 40 % of the sum of the total natural recharge and the cumulative V AR. As expected, after 10 years, the largest cumulative V AR and dV lens are present in scenario C, whereas the smallest V AR and dV lens are present in scenario B. The cumulative S eff in the three scenarios is higher than the yearly S eff. This is indicated by the difference in intersection of the lines of V AR and dV lens between Fig. 14a–c and Fig. 14d–f. In Fig. 14e, the lines of V AR and dV lens do not intersect. Hence, cumulatively, the total volume of the increase of the freshwater lens below the CARD system in scenario B is larger than the total volume of artificial recharge over a period of 10 years. Summarizing, the results of Figs. 13 and 14 indicate that the higher the water table is maintained using the CARD system, the larger the increase of the freshwater lens, the more water is needed for artificial recharge, and the lower the storage efficiency.

Discussion

Increase of extraction due to increase of the freshwater lens

The results of this study indicate the potential of the CARD system to increase a freshwater lens below a creek ridge on a long (10-year) term. Although the focus of the study was on increasing a freshwater lens using the CARD system, the implications for increased freshwater supply are discussed here. An analytical solution is used for this purpose, in view of the long simulation time of the numerical model that is described in this paper. As the area of the artificial recharge is large and the artificial recharge is applied in a relatively shallow aquifer, a horizontal well is considered for the extraction of the fresh groundwater. Compared to vertical wells, horizontal wells induce less drawdown per volume of abstracted groundwater (Yeh and Chang 2013), which in turn reduces saltwater upconing (Dagan and Bear 1968). The following analytical solution (Dagan and Bear 1968) is used for the quantification of saltwater upconing using a horizontal well:

$$ \zeta =\frac{Q_{\mathrm{d}}}{\delta \pi \sqrt{K_{\mathrm{h}}{K}_{\mathrm{v}}}}\underset{0}{\overset{\infty }{\int }}\frac{1}{\lambda}\frac{ \cosh \left[\lambda \left(a-d\right)\right]}{ \sinh \left(\lambda a\right)}\left[1- \exp\ \frac{-\lambda t}{\frac{n}{\delta {K}_{\mathrm{v}}} \coth \left(a\lambda \right)+\frac{n}{\delta {K}_{\mathrm{v}}} \coth \left(b\lambda \right)}\right] \cos \left(\lambda x\sqrt{\frac{K_{\mathrm{v}}}{K_{\mathrm{h}}}}\right)\mathrm{d}\lambda $$
(5)

in which Q d [L3/L/T] is the volumetric extraction rate of the horizontal well per unit length of the well, t [T] is the time since the start of pumping, x [L] is the horizontal coordinate (the well is located at x = 0), a [L] is the vertical distance between the well and the water table above the well at t = 0, d [L] is the vertical distance between the well and the depth of the sharp fresh–saline interface below the well at t = 0, b [L] is the vertical distance between the depth of the sharp fresh–saline interface and the hydrogeological base below the well at t = 0, ζ [L] is the vertical rise of the sharp interface relative to the depth of the sharp interface at t = 0, and λ is the variable of integration. Equation (5) is based on the sharp-interface approximation and is only valid up to certain percentage of d. For further details on this equation, the reader is referred to Dagan and Bear (1968) and Werner et al. (2009). Dagan and Bear (1968) recommended using Eq. (5) for ζ < d/3. For larger values of ζ, an undesired salinity can be expected in the well. As a conservative condition for avoiding an undesired salinity in the horizontal well, ζ < d/4 is used here.

As an approximation, D mix is taken here as the depth of the sharp fresh–saline interface. At ML 1, D mix was about −13 m NAP at the end of the initialization model (Fig. 11). Furthermore, the following parameters are chosen based on the initialization model: a = 13 m, b = 37 m, δ = 0.0223, d = 7 m (so the horizontal well is located at a depth of −6 m NAP), K h = 7.0 m/d, K v = 5.0 m/d, and n = 0.3. Below the well, at x = 0, the maximum pumping rate Q max can be calculated over a period (t) of 180 d in case ζ = d/4. Using these values, application of Eq. (5) yields a Q max of 0.33 m3/d, per meter length of the horizontal well.

The effect of the increase of the thickness of the freshwater lens (i.e., an increase of d and a decrease of b) on Q max is shown in Fig. 15. For scenario A, the 6–8 m increase of the thickness of the freshwater lens that was simulated after 10 years below the centre of the creek ridge (Fig. 13) allows for a 2–3 times higher Q max, as D mix in the initialization model was between −13 and −15 m NAP here. This indicates the potential of the CARD system to increase freshwater supply.

Fig. 15
figure 15

The increase of the maximum pumping rate (Q max) as a function of the increase of the thickness of the freshwater lens, inferred from Eq. (5)

Limitations of this study and recommendations for further research

Although the results of this study have indicated the potential of the CARD system to increase a freshwater lens below a creek ridge and to increase freshwater supply, the study has limitations. First of all, the effects of chemical, biological, and mechanical clogging of the CARD system were not investigated. Belcher and D’Itri (1995) described the ambiguous experiences regarding the influence of clogging on controlled drainage in their summary of a collection of field experiments on controlled drainage and sub-irrigation systems. Clogging will depend on the local conditions such as soil type, artificial recharge regime and the amount of suspended solids in the water. Further research is needed to determine the effect of clogging of the CARD system. Furthermore, although the CARD system is intended to be a relatively cheap method for artificial recharge, the economic aspects of the CARD system need to be investigated in detail.

The hydrological analysis that was presented also has limitations. An important reason for this is the long simulation time of the numerical models that were used. For this reason, the amount of simulations for the calibration, the sensitivity analyses and the scenarios were limited. Future work should focus on constructing a computationally efficient model which can be used to better estimate the uncertainty and sensitivity of the model parameters. Such a model can also be used to more extensively investigate the influence of the availability of surface water for artificial recharge on the thickness of the freshwater lens, and for a more thorough analysis of the increase of the maximum (sustainable) freshwater extraction from freshwater lenses. Although the aspects that control the thickness of the freshwater lens below a creek ridge were not investigated quantitatively to a large extent, some hydrogeological aspects are qualitatively discussed here.

It is expected that low-permeability layers, in particular their spatial coherence and permeability, play an important role in the thickness of the freshwater lens and, hence, the effect of the CARD system on the increase of the freshwater lens. Low-permeability layers can induce significant vertical hydraulic head gradients, such that the vertical pressure distribution in the freshwater lens is not hydrostatic. This has a negative effect on the thickness of the freshwater lens. Furthermore, low-permeability layers can retard or even prevent the increase of the freshwater lens. The field observations at ML 2 suggest an influence of thin low-permeability layers on the freshening process. For applying the CARD system elsewhere, it is important to realize that in hydrogeological models these thin layers may not always be present, as was the case in this study. Cone penetration tests offer a relatively cheap and effective technique with very high vertical resolution, for detecting (thin) low-permeability layers.

The increase of the thickness of the freshwater lens is furthermore influenced by the permeability of the sandy creek deposits and the permeability of the semi-confining layer. The lower the permeability of the sandy creek deposits and the confining layer, the less water (and energy) is needed to maintain an elevated water table to increase the thickness of the freshwater lens. A low permeability of the sandy creek deposits has a positive effect on the amount of water that is needed to maintain an elevated water table, but has a negative effect on the rate of the freshening. The rate and extent of the freshening furthermore depend on how much the water table can be raised during the year, as was indicated by Fig. 14. It is expected that an increase of the period where artificial recharge is applied increases the thickness of the freshwater lens, as this influences the time-averaged water table below the creek ridge, but this was not further investigated. Furthermore, for applying the CARD system elsewhere it is important that the unsaturated zone is thick enough to accommodate the rise of the water table.

Conclusions

The results of this case study indicate the potential of a controlled artificial recharge and drainage (CARD) system to increase a freshwater lens below a creek ridge. On the field site where the CARD system was tested, in the southwestern part of the Netherlands, the 13–15-m-thick freshwater lens is expected to increase 6–8 m within 10 years, as a result of the 0.5 m rise of the water table in winter in the CARD system, based on the results of numerical model simulations. A lower water table leads to a lower increase of the volume of the freshwater lens, but a higher storage efficiency (the volume of the increase of the freshwater lens relative to the volume of artificial recharge). A higher water table has the opposite effect. Based on an analytical solution, the increase of the freshwater lens after 10 years allows for a 2–3 times higher pumping rate using a horizontal well, which indicates the potential for the CARD system to increase freshwater supply.

It is expected that the most important factors that influence the effect of the CARD system on the increase of the freshwater lens are low-permeability layers, the availability of surface water for artificial recharge, the permeability of the sandy creek sediments, the permeability of the semi-confining layer, and the thickness of the unsaturated zone to accommodate the rise of the water table.