Abstract

In coastal-zone fields with a high groundwater level and sufficient rainfall, freshwater lenses are formed on top of saline or brackish groundwater. The fresh and the saline water meet at shallow depth, where a transition zone is found. This study investigates the mixing zone that is characterized by this salinity change, as well as by cation exchange processes, and which is forced by seepage and by rainfall which varies as a function of time. The processes are first investigated for a one-dimensional (1D) stream tube perpendicular to the interface concerning salt and major cation composition changes. The complex sequence of changes is explained with basic cation exchange theory. It is also possible to show that the sequence of changes is maintained when a two-dimensional field is considered where the upward saline seepage flows to drains. This illustrates that for cation exchange, the horizontal component (dominant for flow of water) has a small impact on the chemical changes in the vertical direction. The flow’s horizontal orientation, parallel to the interface, leads to changes in concentration that are insignificant compared with those that are found perpendicular to the interface, and are accounted for in the 1D flow tube. Near the drains, differences with the 1D considerations are visible, especially in the longer term, exceeding 100 years. The simulations are compared with field data from the Netherlands which reveal similar patterns.

Introduction

Deltas worldwide are areas in which fresh and saline water meet and mix. Because of the generally dense population in coastal areas, the availability and quality of freshwater is important for drinking water and agricultural production. This availability and the threats caused by salinization have received worldwide attention over the last few years, for example Alvarez and Dapeña (2015) describe the main causes of salinization in coastal zones. The consequences of an increasing sea level and climate change for salinization of coastal zones has been investigated over the past years by Giambastiani et al. (2013), Langevin and Zygnerski (2013), Colombani et al. (2014) and Feseker (2007) and De Louw et al. (2011, 2013) for deltas in Italy, Florida (USA), Argentina, Germany and The Netherlands respectively. Scott et al. (2014) investigated the development of shallow brackish water in Bangladesh. Lenses, formed by infiltrating rainwater ‘floating’ on saline groundwater, are often a major source of freshwater. The freshwater lenses, sometimes called Badon Ghijben-Herzberg (BGH) lenses (Drabbe and Badon Ghijben 1889; Herzberg 1901) after the first investigators of these freshwater bodies, are formed by rainwater in areas that are somewhat elevated with respect to mean sea level (MSL), and have an annual precipitation surplus. These lenses are found on (barrier) islands (e.g. Underwood et al. 1992), and in dune areas (e.g. Stuyfzand 1993; Vandenbohede et al. 2008; Röper et al. 2012). In low-lying areas, with a soil surface close to mean sea level (MSL), thin lenses are found, on top of upward seeping saline groundwater (Maas 2007; Eeman et al. 2011; De Louw et al. 2011). Also in the absence of saline groundwater, lenses can be formed on deeper, upwelling groundwater and differ from that groundwater by its chemical signature, which can be ecologically relevant (Cirkel et al. 2014).

These thin lenses, referred to as fresh or rainwater (RW) lenses, are studied because of their vulnerability during expected future droughts. Sensitivity (Eeman et al. 2011), occurrence and spatial variability (De Louw et al. 2011) and response to recharge variations, including delay and damping, have been investigated in detail (Eeman et al. 2012). Geological and topographical variations have been shown to have a major impact on the thickness of RW lenses (De Louw et al. 2011).

De Louw et al. (2013) studied the dynamics and related mixing in two agricultural fields in Zeeland, The Netherlands, and are studied in this paper as well. They define the extent of a RW lens as the entire groundwater body from the base of the mixing zone (Bmix, Fig. 1, below which no mixing with rainwater occurs) to the water table. This definition is adopted in the current study, and implies that RW lenses generally do not contain truly freshwater. They range from brackish (mostly 0–2 g/l Cl, but in some cases up to 8 g/l Cl) at the groundwater level (GWL), to saline (>10 g/l Cl) at a depth of 2–4 m. De Louw et al. (2013) found that variations in the position of the mixing zone, and mixing zone salinities, are small and vary on a seasonal timescale, which is attributed to the slow transient oscillatory flow regime in the deeper part of the lens. The flow and mixing processes are much faster near the water table, which fluctuates on a daily basis in response to recharge and evapotranspiration. Similarly, variations of drain-water salinity are very dynamic as they respond to individual rain events.

Conceptual visualization of a rainwater (RW) lens on top of upward-seeping saline groundwater. The RW-lens extends from the base of the mixing zone (Bmix) to the water table. At the depth of the center of the mixing zone (Dmix) the salinity is 50 % of the seepage water salinity (ECs) (taken from De Louw et al. 2013)

The ion exchange associated with this dynamic mixing process affects the hydrochemistry that determines the water quality in the root zone as well as the drain water composition. European regulations, such as the EU Water Framework Directive, and the vulnerability of deltaic systems in view of climate change demand a better understanding of chemical processes in RW lenses (Bates et al. 2008).

Although the authors are not aware of any published studies on hydrochemical processes in thin RW lenses in areas with saline seepage, hydrochemical processes have been studied in thick BGH lenses because of their influence on drinking-water quality. Stuyfzand (1993) investigated the hydrochemistry of dune and adjacent areas, where mixing of saline and freshwater is an ongoing, although slow process. He found that hydrochemical processes in mixing zones are in a transient state, even where groundwater flow is relatively stable, due to the slow kinetics of dissolution and precipitation processes. Marconi et al. (2011) and Mollema et al. (2013) showed for small coastal wetlands and forests in the Po delta and the Ravenna area, respectively, both in Italy, that the hydrochemical conditions are different. Both used the Base Exchange index (BEX index) of Stuyfzand (2008) to characterize the water types found, and concluded that the combination of land use, climate and sea level rise is causing salinization in the aquifers in this area. Mollema et al. (2013) point out the influence of seepage of hypersaline water on coastal freshwater systems caused by artificial drainage of agricultural land.

Hydrochemical conditions in RW lenses differ from those in a BGH lens, where mixing occurs at a much larger depth, often with anoxic conditions and relatively steady situations. RW lenses, on the contrary, are much more affected by oxic rainwater that mixes with anoxic seepage water. They are also a much smaller size than BGH lenses in horizontal (10’s of meters instead of 100’s of meters) and vertical (less than a few meters instead of 10–100’s of meters) planes, causing travel times of less than 5 years (De Louw et al. 2013) instead of decades to centuries (Stuyfzand 1993; Vaeret et al. 2011; Oude Essink 1996).

Foppen and Griffioen (1995) stressed that the role of drains in hydrochemical processes needs elucidation. A clear demonstration of the importance of RW-lenses on the local hydrogeochemistry of regions with upward groundwater seepage is given by Cirkel et al. (2014). They show the resilience to eutrophication by fresh groundwater, saturated by calcite and rich in sulphate, due to the precipitation of pyrite and calcite invoked by processes in the organic matter rich top soil. In another study, Cirkel et al. (2015) simulated the movement and shape of oscillating solute fronts with a significant change in concentration, taking into account exchange of Na+ and Ca2+. They found that the transition zone thickness is unaffected by the ion-exchange process, but the concentrations of individual cations change profoundly within the zone, influenced by the phase of oscillation. Hence, such patterns may also be expected in the case of fresh/salt gradients. Many other studies focus on exchange processes (e.g. Bolt and Bruggenwert 1976); however, these studies concern steady flow situations.

The objective of this study is to investigate the dominant hydrochemical process (ion exchange) as affected by mixing and fluctuations of rainwater infiltration and upward saline seepage. Such situations abound in practice in coastal regions (with fresh/salt gradients), but are similar with those in brook valleys (with different chemical signatures of shallow and deeper groundwater) and even with drinking-water-well extraction/injection cycles. For lenses, the gradients as affected by mixing, fluctuating flow directions and exchange, are ecologically and agronomically important. Because a preliminary analysis of field data revealed spatiotemporal complexity, this study focused on a simplified one-dimensional (1D) vertical situation, to begin with: the development of the fluctuating transition zone in the solution and the soil adsorption complex where seepage and infiltration water meet. The main ions, Na+, Ca2+, Mg2+ and K+ are considered, which is an extension of the work by Cirkel et al. (2015). This extension is not trivial, because consideration of more cations leads to much more combinations of initial and boundary conditions that affect the mixing zone. Using the insights obtained with this first model assessment, a more complicated two-dimensional (2D) flow pattern was analyzed, using parameters typical for a site in Zeeland in the Netherlands. Theory based on geochemistry and numerical transport modeling was in this way applied to the chemically complex, cyclic situation of a 2D tile-drained field.

Study sites

The two field sites have been described by De Louw et al. 2013. Both sites are situated on the island of Schouwen-Duivenland in the south-western part of the Netherlands (Fig. 2a). The mean annual precipitation and potential (Makkink 1957) evapotranspiration for the period 1980–2010 are 815 and 613 mm, respectively (data from the Royal Netherlands Meteorological Institute, KNMI). Groundwater is predominantly saline, and originates from Holocene transgressions during the periods 7,500–5,000 BP and AD 350–1000 (Vos and Zeiler 2008; Post and Kooi 2003). Around AD 1000 the island formed, and the first polders and drainage of the land started in the 12th century. Intensification and deepening of the drainage system developed in the 20th century. Tile drainage, as used for both fields, started in the 1970s. Upward seepage of saline groundwater occurs on the lower part of site A and on site B, where the land surface lies below mean sea level. RW-lens dynamics and mixing behavior on the sites (A and B) were monitored from March 2009 to December 2010 (De Louw et al. 2013, their Fig. 2b); for detailed information on the dynamics in space and time of these lenses, reference to their paper is made.Hydrochemical data were collected during two field campaigns in March and May 2009 and soil samples were taken in June.

a Location of study area. b Cluster of rhizons and piezometers installed at different depths. c Set up of monitoring network of site A, the cross section with water table, vertical flow direction and Bmix according to De Louw et al. (2011). d Setup of monitoring network of site B (top view) (Source: De Louw et al. 2013)

Site A

Site A is a 300-m-long agricultural field located at the transition from a former salt marsh, presently at an elevation of −0.7 m MSL, to an inverted creek paleo-channel, presently at −0.2 m MSL (Fig. 2c). Only data from the low-lying clayey zone are presented in this study. The field is bordered by ditches with a surface-water level which is maintained at an elevation of −2.0 m MSL, and is tile-drained at a depth of 1.0 m below soil surface with a drain distance of 10 m. Fine-grained to medium coarse-grained sand is found between 2.5 and 30 m depth (part of a regionally extensive aquifer), confined by 2.5 m of heterogeneous low-permeability sediments comprised of peat, clay and sandy clay at the former marsh (Fig. 2c). The measured hydraulic freshwater heads in the former marsh area were permanently higher in the sandy aquifer at 4 m depth than at 1.5 m depth in the aquitard, resulting in a continuous upward groundwater flow (De Louw et al. 2013). The RW lenses below the low-lying clayey zone are thin, between 1.5 and 3 m.

Site B

Site B is located in a reclaimed salt marsh at an elevation of −2.0 m MSL and is bordered by a ditch with a surface-water level maintained at −3.1 m MSL. Drainage tiles with a length of ∼80 m are located at a depth of 0.70 m with a horizontal separation of 10 m (Fig. 2d). The upper 2.5 m of the subsoil consist of silty clay and a 10–15-cm thick peat layer at about 0.5 m depth. Shrinkage cracks in the clayey top soil were observed during the drier part of the year (May–October). Between 3.0 and 4.5 m depth, a sequence of clayey fine sand and sandy clay is found. These low-permeability sediments form a semi-confining unit, below which the same regionally extensive sandy aquifer as at site A is found. Permanent upward seepage occurs everywhere below the agricultural field, and the RW-lenses have a thickness between 1.5 and 2 m (De Louw et al. 2013).

Hydrochemical data collection and laboratory measurements

At sites A and B, 14 and 11 sampling locations were selected, respectively, in March and May 2009. These were chosen based on the expected differences between location near and in between drains and ditches and illustrated as MP (measurement point) in Fig. 2c,d and at different depths between 1 and 4 m (Fig. 2b). In the field, O2, temperature, pH, electrical conductivity (EC) and alkalinity were measured directly. Cation samples were collected and analyzed in the laboratory. Concentration measurements were done using ion chromatography (SO42−) and ICP-AES (Na+, K+, Ca2+, Mg2+). In this study, only a limited number of the measurements have been used, from the dominating cations Na+, K+, Ca2+, Mg2+ and the pH and SO42−, for the sampling sites and depths that can be related to the drainage profile described in the previous section. For site A, results of MP 1, 4 and 6 were used and for site B for MP 2, 3, 4 and 5.

The data in Table 1 reveal a significant difference in salinity between the locations. Dilution due to mixing with rain or river water as the sediment formed causes site A to have concentrations of approximately two thirds of seawater salinity, whereas at site B it is only half. This mixing with other water sources can also explain the relative abundance of magnesium in the soil at site B.

Table 1

Composition of water use for input, based on field measurements and equilibrium calculations

Species

Unsaturated zone water

Seepage water, site A

Seepage water, site B

(mmol/kgw)

(mmol/kgw)

(mmol/kgw)

Na

3.4

319.9

270.5

K

0.26

6.9

4.7

Ca

3.4

18.1

15.2

Mg

1.8

21.1

30.6

Cl

7.9

375.4

344.5

SO4

0

10.3

7.5

CO2

1.6

2.3

0.96

pH

7.0

6.7

7.1

pe

14.8

14.6

14.7

In June 2009, soil samples were taken at both sites at depths of 55, 95, 135, 165 and 195 cm at MP 3 and MP 4 at site A (Fig. 2c) and MP 3 and MP 5 at site B (Fig. 2d). Cores were analyzed in the lab according to Tan (1996). The concentration of exchangeable cations (Na+, K+, Ca2+, Mg2+) and effective cation exchange capacity (CEC, cmol+/kg) were determined by ion displacement using an adapted version of the unbuffered compulsive exchange procedure, as proposed by Gillman (1979). ICP-AES was used to measure concentrations of the exchangeable ions. Results were corrected for the presence of excess salts in the examined cores.

Theory and methods

This section presents the used flow and transport models, describes the input for the hydrochemical model, and introduces the theoretical concept used to interpret model results.

Flow and transport models

To understand and visualize the flow patterns and exchange processes in the studied fields, the data gathered in the field were used as the basis to construct two hydrochemical flow and transport models. The PHT3D software package (Prommer and Post 2010) is used to combine the possibilities of flow and transport from MODFLOW and MT3D (Zheng and Wang 1999), with the calculation of the chemical processes by PHREEQC2. For a 1D column, exchange and how this is influenced by alternating flow direction and input of recharge from the top and seepage from the bottom, is visualized. With a 2D model the study site B is represented, where flow towards drains and constant upward seepage of saline water creates different mixing patterns and associated cation exchange processes (Eeman et al. 2011) compared to the 1D column. Although there is a difference in the density of the unsaturated zone water and the seepage water, this has been neglected in the models, since De Louw et al. (2011) and Post et al. (2007) showed that the studied field systems are controlled by head gradients and geology, buoyancy effects are not significant.

1D column basic model

A 1D vertical-column MODFLOW model was made, consisting of 41 square cells with sides of 0.1 m (Fig. 3), resembling the flow-tube concept used by PHREEQC2. All side boundaries were no-flow and no-solute flux boundaries. Inflow and outflow, due to alternating flow directions, were limited to the top (recharge) and bottom (seepage) boundaries. Inflow from the top results in an equal amount of outflow at the bottom and vice versa. The column remains completely filled with water at all times, and effective porosity was set to 35 %. Apart from porosity, hydrogeological parameters do not influence this artificial model, since water is forced by the flux conditions at top and bottom boundaries. This column does not represent a particular site, and its aim is to illustrate the effect of fluctuations on a fresh/saline front in a basic way. The flow was calculated using stress periods of 1 day and transport and chemistry time steps were 1 h.

For an initially saline column, filled with seepage water, three different alternating flow patterns were applied to investigate both the influence of the average flux and variations around the average. The patterns, shown in Fig. 4 for a period of 2 years, are repeated to calculate the processes in the column for 300 years. This long period is needed since ion exchange progresses very slowly compared to flow and transport. The water compositions of seepage and recharge water are based on observations from study site B and A respectively and elaborated in the next section.

Flux patterns used in the column simulations. Positive values relate to downward flow

The first simulation (red line in Fig. 4), applied an average zero flux and a step-wise change of flow direction of 1 mm/day downward (recharge, implying outflow of saline water at the bottom) and 1 mm/day upward (seepage, leading to outflow of relatively freshwater at the top). Alternation of flow direction occurred every 3 months. The mixing in this simulation resembles the process in the center of the mixing zone, in De Louw et al. 2013. The second simulation (green line in Fig. 4) was similar, only the average flux is 1 mm/day downward, comparable to the mixing in the top of the mixing zone, in De Louw et al. 2013. Again a step-wise 3 month change of flow direction has been used of 1.5 mm/day downward and 0.5 mm/day upward. The third simulation used the daily precipitation measured at field site B, combined with Makkink evaporation data from the nearest KNMI weather station, also on a daily basis (WiIheminadorp, 17 km away). The average flux of these data is 1 mm/d downward, equal to the second step-model. The amount of evaporation in this case equals the amount of seepage water coming from the bottom. The blue line in Fig. 4 shows the variation around this average.

Two-dimensional field model of profile between drains

The 2Dmodel of site B is linked to the 1D model, and an extension from the model proposed by De Louw et al. (2013). Hydraulic properties for this model are illustrated by Fig. 5. Half of a cross-section between drains was modelled, 5 m wide and 4.1 m thick, with a mesh of 41 layers of 0.1 m thick and 21 columns of 0.25 m wide. Hydraulic properties were based on auger hole descriptions and verified using REGIS II, the database for Dutch hydrogeological information. Total porosity at the site was measured at 49 %, effective porosity was estimated at 25 %, in line with findings of Stephens et al. (1998). The anisotropy factor, horizontal conductivity over vertical conductivity (kh/kv), was set to 5 for the entire model (Vernes and Van Doorn 2005). Vertical upward seepage into the model is calculated based on field measurements, at a value of 0.29 mm/day. Recharge and evaporation used are the same as described for the third simulation of the 1D model: field measurements combined with KNMI data from Wilhelminadorp. Recharge was applied at the highest active cell in MODFLOW, using the option that makes cells inactive as the head drops below their bottom. Cells are re-activated again as the head in the underlying cell rises above the cell’s bottom (Harbough et al. 2000). A specific yield (S) v depth relation was applied, illustrated in the inset of Fig. 5. The drain depth was at 0.7 m and the side boundaries were no-flow and no-solute flux boundaries. All inflow, i.e. recharge and vertical upward seepage, were applied at the top and bottom model boundaries, respectively; water leaves the system by tile-drainage and evapotranspiration.

a The schematic representation of the model of a RW-lens between two tile-drains at site B. b Specific-yield–depth relation. Source: De Louw et al. 2013

Hydrochemistry models

The water quality in the freshwater lenses at site B is evaluated by combining field data and lab measurements (see section ‘Study sites’) with the flow and transport models (see section ‘Flow and transport models’), using the hydrochemical package PHREEQC (Parkhurst and Appelo 1999).

The data from site B are the most valuable, due to the setup of the measurements above and between drains, the larger number of data at these MPs, and the calibration of the transport model based on drain data (De Louw et al. 2013). Although the 1D simulation is a generic one rather than a site-specific representation, the hydrochemical data from site B were adopted for the seepage water, which is also true for the 2D field model of site B.

Initially, all models are filled with saline seepage water. The initial composition of the cation exchange complex in the models was calculated with PHREEQC2, using the composition of the seepage water based on averages of the different MPs on both sites, additionally, representative CEC values for the upper 2 m of the profiles were also used. The small amount of available soil data, two locations per site, did not justify more detail. The CEC was on average much larger for site A (0.5 eq/l) than for site B (0.13 eq/l), which is explained by the substantial amount of peat in the subsoil of site A (a layer of approximately 0.7 m thick, see Fig. 2c). The upward seeping water entering the model at the bottom equals the seepage water quality that initially fills the model.

The water entering the model at the top is similar to rainwater, but enriched in Ca2+ while passing the unsaturated zone, according to the unsaturated water sample from site A. The reason that a sample from site A here was used is that no suitable sample of unsaturated zone water was available at site B. The unsaturated zone water in the calculations was assumed to be equilibrated with CO2 and calcite, to determine the composition of the input recharge water in the model. CO2 pressure was determined using the average partial pressure of the five samples that were taken on site A, according to Henry’s law (NIST 2015). Results for the three water types used as input for the models are in Table 1.

Significant chemical precipitation is very unlikely to occur in view of the very small (<1) saturation indices (SI) of the unsaturated zone and seepage water at both sites for calcite and dolomite. For all other minerals, the SI was negative; however, to rule out errors due to the impact of chemical precipitation, this was tested in preliminary simulations. There was no significant impact on model outcomes and therefore chemical precipitation has not been accounted for. If dissolution or precipitation occurs, this is likely to be a non-equilibrium process (particularly for dissolution), which will affect the concentrations of the cations; however, it modifies but does not in essence change the processes. Another important aspect is that the effect of DOC/SOM-mediated redox processes was not accounted for, which does not imply that this can be considered as an insignificant process. However, DOC-facilitated transport may affect the solute balances, but not the impact of cation exchange as a predominant mechanism in the root zone.

Analysis of the simulations

To describe exchange processes, the Gaines-Thomas equation is used, which is available in PHREEQC2, and for Na+ and Ca2+ is given by

where X stands for surface site. To interpret the outcomes of the hydrogeochemical models, the systems is made dimensionless as was proposed by Cirkel et al. (2015), using C for the total concentration of cations in solution, f for the fraction of dissolved species i with concentration ci and N for the fraction of the cation exchange complex occupied by adsorbed species i

To investigate the modelled and observed changes in soil water and exchange complex due to varying lens thickness, the 1D oscillation is considered, as recently investigated by Cirkel et al. (2015), comparing relative normalized amounts of Ca2+ in solution (f) and on the complex (N). As explained in the transport section, initial conditions are saline. The Ca2+ enriched rainwater that reaches the groundwater (Table 1) is often found in coastal zones, due to the marine history of the sediment (Appelo and Postma 2005). As this Ca2+ enriched water enters the flow domain, two changes occur: (1) the salinity decreases and (2) the fraction \( {f}_{{\mathrm{Ca}}^{2+}} \) in solution increases. Only the results for the three main cations—Ca2+, Mg2+, and Na+—are described, and after preliminary calculations it was decided to ignore K+ in view of its tenfold smaller concentrations (Table 1) and accordingly small impact; however, K+ is included in the total concentration C.

Results

Transport in a 1D column with time varying boundaries

To understand the results obtained for the column simulations, the concept is first presented and illustrated in Fig. 6. Due to the decrease of salinity, the preference of the cation exchange complex for Ca2+ increases (Bolt and Bruggenwert 1976; their Figure 7.5) and some of the dissolved Ca2+ adsorbs, whereas some Na+ desorbs. The increase of \( {f}_{{\mathrm{Ca}}^{2+}} \) also causes Ca2+ to adsorb. The net effect is therefore that it may depend sensitively on the quantity of incoming Ca2+ and on the quantity that is exchanged with Na+ due to the two changes in the solution (different cationic ratios and different salinity), whether the net effect is an increase or a decrease of \( {f}_{{\mathrm{Ca}}^{2+}} \); however, in either case, the fraction \( {N}_{{\mathrm{Ca}}^{2+}} \) of Ca2+ on the complex increases. Next, as the fresh/salt interface oscillates, the original saline solution displaces the freshwater that infiltrated, causing some re-adsorption of Na+ but also removing in this process some of the desorbed Na+ from the soil at the top of the flow domain. If a full oscillation is considered, some Na+ will be removed and some Ca2+ will have entered the flow domain, and particularly at shallow depths, the salinity will decrease somewhat. In Fig. 6b, which schematically depicts the trajectory of Ca2+ changes, the impact of changing salinity and cationic composition is illustrated.

a A vertical profile with two observation depths and b the corresponding relation between N and f as both waters start mixing due to an alternating flow direction

As freshwater enters at the top, the changes are largest there (red line); however, due to dispersion, and particularly if the net downward flux is positive, these changes will propagate downwards (green line). Gradually, the fraction of Ca2+ in the solution and at the cation exchange complex will increase and to which level this occurs depends on the salinity profile that will develop between the fresh incoming and saline initially present pore water.

Situation 1: alternating flow, no net flux

Figure 7a illustrates the development of dimensionless adsorbed fractions (N) as a function of dissolved (f) fractions of divalent Ca2+, Mg2+ and Na+, for the simulated regular alternation (every 3 months) of flow direction (1 mm/day), without net displacement of water. The initial saline situation is indicated by the black circle. Different colors correspond with different depths. Figure 8a shows the relation between the dissolved fractions f and the total concentration in the water, C, for the same simulation.

Fractions f of dissolved Ca2+ (left panel), Mg2+ (middle panel) and Na+ (right panel) versus the total amount of dissolved salts at different depths, for a simulations of 300 years. a Step-wise change of flow direction every 3 months (average 0 mm/day); b Step-wise change of flow direction every 3 months (average 1 mm/day)—the sharp angle is reached after approximately 10 years, depending on depth. c Actual net precipitation simulation (average 1 mm/day)—the sharp angle is reached after approximately 5 years, depending on depth

The complex depth-dependent pattern of Fig. 7a can be interpreted on the basis of sequential processes. For the top of the domain (zero depth), the major short term impact of incoming water is that the total salt concentration decreases, increasing the preference of the adsorption complex for divalent Ca2+ and Mg2+ at the expense of monovalent Na+. Both divalent cations therefore adsorb and their f-values decrease. Because the distribution ratio (amount of cations at the complex divided by their amount in solution, is large (\( \frac{\mathrm{CEC}}{C} \)= 16.5 in this model), the change in N-values is small compared with the change of f -values. Na+ on the other hand, initially shows an increase in f, as it is exchanged due to the preferentially adsorbing divalent cations. As the infiltrating freshwater reaches only a limited depth, and due to dispersion, the decrease of salinity is smaller for larger depths, which is recognizable in the smaller short-term impact at larger depths. After a relatively small number of oscillations, the fluctuation of the total salt concentration stabilizes, though at a different level for each depth.

In the longer term, another process becomes important and that is the cation exchange, caused by the different f -values of initially present and incoming freshwater. As Table 1 and Fig. 7 reveal, Ca2+ is far more abundant in the incoming freshwater than Na, and also more abundant than Mg. Accordingly, the f -value of Ca2+ starts to grow, and Ca2+ adsorbs, predominantly at the expense of Na, which is also the case for Mg2+, which likewise exchanges with Na+ at the complex, but to a lesser degree than Ca2+. In the long term, when Na+ has been by and large expelled from the exchange complex, Ca2+ also starts to outcompete Mg2+, and the N -value of Mg2+ declines again. As is apparent from Fig. 7 for Mg, the maximum N -value attained increases with increasing depth, which is partly the result of passing Mg2+ that is desorbed higher in the soil profile, and partly due to the complex interplay of changing concentrations of Ca2+, Mg2, Na+, and total salinity.

Figure 8a is in agreement with the aforementioned interpretation: in the short term (related to the half cycle of 3 months), f -values decrease and increase sequentially, due to varying salinity, where the salinity variation is depth dependent. As divalent cations in soil build up their mass, which is possible due to the nonlinearity of the cation exchange, the f -values of Ca2+ and Mg2+ increase and that of Na+ decreases; however, this occurrence is strongly retarded in time, due to the large distribution ratio. Figure 9a more clearly shows that the composition of the complex N shifts from Na+ towards Ca2+, while Mg2+ increases only slightly, while the changes at larger depth are almost negligible.

Situations 2 and 3: alternating flow with net downward flux

Figure 7b,c again illustrates the development of divalent Ca2+ and Mg2+. The net water displacement is 1 mm/day based on the regular pattern of alternating flow every 3 months (situation 2) or on alternation triggered by meteorological data (situation 3), as described in section ‘Theory and methods’.

Although the processes for situation 1 (Figs. 7a, 8a, 9a) are similar to situation 2 and 3 (Figs. 7b,c, 8b,c, 9b,c), the differences are striking: the saw-tooth pattern is far less distinct in the case of a net water displacement and the changes occur over a different range of f and N values. This illustrates the importance of net water displacement, which clearly dominates over the changing flow direction that occurs on the timescale of months.

For Fig. 7b,c, the major short-term impact is still determined by the decrease in total concentration due to the incoming water. Na+ initially increases in the solution (f) since the changing isotherm causes Na+ to be pushed from the complex by Ca+ and Mg+. The large decrease in f of Ca+ and Mg+ compared to a much smaller decrease of N for Na+, due to the large distribution ratio, is visible as well. However, in this case this effect is significant to a depth of 1.5 m, since the infiltration reaches progressively deeper due to the net water displacement in the downward direction.

In the longer term, the oscillations in f are absent, because the salinity at the top of the domain never reaches the initial situation again, a direct consequence of the net downward flux. Compared to the zero net displacement simulation, Fig. 7b,c shows a rather smooth development of the exchange process between monovalent Na+ and divalent Mg2+ and Ca2+. The initial differences in the N –f trajectory for the considered depths disappear and the patterns converge to the same one, with a delay related to depth. The interplay of concentrations is slightly different for the regular oscillation (Fig. 7b) and the natural precipitation simulations (Fig. 7c). Changes in N occur faster at larger depths for the precipitation model, and the maximum concentration of Mg2+ increases sharper with depth when compared to the regular model. This can be explained by the shorter periods of upward flow, which cause less reversal of the exchange process compared to the tri-monthly alternating flow direction, causing also a difference in the balance between the three main cations.

Figure 8b,c confirms the interpretation of Fig. 7. Total concentration C drops quickly for both simulations due to the incoming freshwater and does not visibly increase again. However, relative concentrations of Ca2+ and Mg2+ develop much slower, continuously increasing for Ca2+ and decreasing again after a maximum for Mg2+. The larger velocity of exchange with depth for the natural precipitation simulation is shown by the left panels of Fig. 8b,c, where the lines showing the deeper positions are more developed in the bottom (precipitation) than in the middle (regular) panel. Just as for Fig. 7, Fig. 8 shows a much smoother development of N, f and C for more realistic flow patterns (given the recharge–seepage situation in the field) compared to the artificial zero-displacement simulation. This difference is not so explicit in Fig. 9b,c. The temporal dominance of Mg2+ with increasing depth is illustrated clearly for the complex. Due to the exchange of Ca2+ higher in the column, less is available at larger depth, causing more Mg2+ to exchange with Na+.

Movement of the front

In Fig. 10, the movement of the total concentration front is shown for the three simulations. Figure 10a shows solid lines for the front position in situation 1 after 3 months of downward flow, and dashed lines for the position of the front after 3 months of upward flow. Sets of solid/dashed lines are shown after each 5-year interval. After 50 years, the lines approach their equilibrium position with a broad transition zone that covers the entire considered domain. This is in line with the smooth development in Figs. 7a and 8a. Figure 10b,c, situation 2 and 3 respectively, show that the net downward flux causes a fast downward progression of the front. In both figures, solid lines are drawn for each year, showing that for the regular oscillation simulation (middle) the domain (until 4 m depth) is almost completely filled with freshwater after 10 years. For natural precipitation, this case is already even faster after 5 years, which can also be seen in Figs. 7b,c and 8b,c, where the sharp angles show that decrease in C ends as the fronts have passed.

a Front movement for the simulation with 0 mm average. The solid lines, approaching the lower limiting straight line, show the development of the front developing each 5 years, for 3 months of upward and downward flux. The dotted lines, approaching the upper limit of the black diagonal bar show the upper limit front that develops after ca. 50 years; the lower end of this bar is the limit after 3 months of downward flow. b Front movement for the 1 mm downward average step simulation. The solid lines show the annual progress of the front. c Front movement for the 1 mm downward average precipitation/evaporation simulation. The solid lines again show the annual front progress

For all cases, Fig. 10 confirms that the short-term impact of the oscillating front both with and without net downward flux can be attributed to the movement of this front, whereas the long-term impact occurs after the front has reached its steady-state position in situation 1, or has passed in situation 2 and 3. Figure 10b,c shows that natural precipitation (bottom) causes faster flushing than the regular model with average downward flux of 1 mm/day (middle), which is in line with the interpretation of Figs. 7 and 8. The middle panel also shows that the bottom will not completely freshen in the regular case with a 3 month alternation. A zone of approximately 15 cm at the lower boundary is refilled with seepage water each cycle (0.5 mm/day, 35 % porosity, 91 days) and is therefore a consequence of putting the boundary at −4 m.

Transport and exchange in a drained field model

Whereas the 1D simulations give an impression of the complexity of the geochemical interactions, they also neglect a crucial aspect: the dimensionality of flow in real situations. Strictly vertical flow may be in the middle between two drains (at the hydrological divide, HD); however, closer to the drains the water movement may be dominated by horizontal flow (Eeman et al. 2011). This section, investigates whether this affects the transport.

Exchange process

Figures 11, 12 and 13 illustrate the results for the 2D field cross section. The results for the HD and halfway between the HD and the drain are very similar, and therefore only the midfield position is shown. The explanation for these similar results is found in the shape of the transition zone between fresh and salt water, which is almost horizontal, apart from the 2 m nearest to the drain. There is a significant difference between the horizontal flow components at the HD, where it is almost zero, and halfway between the drain and the HD, where the horizontal flow component dominates. The reason that the horizontal flow component has a limited impact on the composition of the transition zone compared with the 1D cases is that this incoming water has more or less the same composition as the water that it pushes away sideward. For the vertical flow component, this is not the case, as higher and lower water is fresher and more saline, respectively. Near the drain, the gradient of the transition zone increases, and flow lines are forced to converge (see also Eeman et al. 2011), leading to increased mixing of water types, less freshening and a slower progress of the exchange between Na+, Ca2+, and Mg2+.

Site B: Development of the ratio of dissolved vs adsorbed Ca2+ for different depths over time. a Shows the situation in the middle of the field at the HD and b shows the situation close to the drain (at 70 cm depth). Points are drawn each 10 days, to include variation within years. The clay layer starting at 170 cm depth forms a boundary, below which little freshwater is found, shown by the very small changes at 185 cm depth after 200 years

Front positions a in the middle of the field and b at the drain, for the site B model. Lines were drawn each week for 200 years. Lines that seem cut off are caused by drying out of the model in a dry period

For Ca2+, the exchange over time and depth are visualized from 85 to 185 cm in Fig. 11. The drain is at 70 cm depth and the clay layer with low permeability (1 mm/day) starts at 170 cm depth (Fig. 5). Above the drain, the profile is mostly water-unsaturated, which leads to inactive cells in the model, for which the chemical balance is not calculated. The top of the relatively impermeable clay layer is a hydrogeological barrier. Little freshwater enters this layer, as shown by Fig. 14a, which is in agreement with the lower limit of the mixing zone found at approximately 2 m depth, according to De Louw et al. (2013), and therefore is also valid for Mg2+ and Na+.

a Field measurements of total concentration (sum of concentrations of Na, Ca2+, Mg2+ and K+). Locations of MPs are in Fig. 2. This figure can be compared to the Cl-fronts illustrated by De Louw et al. (2013) in their Figure 4, but only cations instead of anions are given for this case. b Comparison of field observations and the 2D simulation for site B. The lines in different colors show the model outcomes for the same date as the field measurements with the same color

When Figs. 12a and 7c are compared, the similarity is striking. Both the short- and long-term processes described for the 1D domain are very recognizable in the case of a 2D field situation with drains and upward seepage instead of a 1D domain in which the saline and freshwater are artificially pumped up and down. Even the velocity of the exchange for Ca2+ is rather similar when the region just below the top of the saturated zone (85 cm depth) is compared to the 25 cm depth line in the 1D domain. At larger depth, the exchange is slower than for the 1D domain, which can be explained by the decreasing vertical displacement with depth due to the constant upward flow of the seepage water. The situation deeper is more comparable Fig. 7a, since net vertical displacement becomes very small. The water in the middle of the transition zone moves mainly horizontally towards the ditch, which is in line with the findings of De Louw et al. (2013) concerning the dynamics of the mixing zone.

For Mg2+, the temporary increase of adsorption on the complex at larger depth is not as profound for the field situation. A slight temporary increase of Mg2+ is still visible, but the relative availability of Ca2+ compared to Mg2+ at larger depth seems larger in a field situation than for a 1D column. This corresponds with the difference in flow pattern between 1D and 2D: whereas the seepage water, containing Mg2+, is forced to flow upward in the column, it follows a nearly horizontal flow direction in the mixing zone for a field situation (De Louw et al. 2013). Less Mg2+ is therefore supplied after the freshwater lens is established, compared to the 1D situation with alternating flow directions.

Movement of the front

For the two sides of the simulated domain (the HD and near the drain) for site B, the front ranges are illustrated in Fig. 13. There is a limited bandwidth in which the total concentration (C) front varies dependent on weather and season. Near the drain, its depth of 0.7 m below the soil surface is determining for the depth of the mixing zone, whereas in the middle of the field, the clay layer between 1.7 and 3 m depth provides a lower limit for the front (Fig. 5).

Model and field data

Figure 14a shows the field data for total concentration, which in Fig. 14b are used to illustrate the match between the 2D model and the field data. The comparison is reasonable: the mixing zone is at the same depth for the field and the data, and for both the model and field, a shallower mixing zone and higher total concentration towards the phreatic level are found near the drain (black and red) compared to the middle between drains (green and blue). This shows that the chemical model captures the main features occurring in the field, since flow and transport were calibrated by De Louw et al. (2013) on drain data that were not presented in this paper. For site A, not enough data were available to calibrate the 2D model; therefore, only the field data and drain depth are shown at the same scale.

Influence of hydrogeological and chemical differences on the exchange process

Since the composition of the saline end members for two sites, B and A, as well as CEC values for both soils are available, a qualitative assessment can be made of the influence of hydrogeology and chemical composition. It is assumed that the calcium-enriched precipitation is equal for both sites, based on measurements at site A and the similarity of the amount of shells found in both top soils. Based on the pore-water composition and CEC (see section ‘Hydrochemistry models’), the equilibrium situation was calculated for the complex and the unsaturated zone water. Although this equilibrium situation may not be reached, it provides a direction of change for the complex composition for both sites. The equilibria for the saline start, differing per site, and “end” situation on the complex are illustrated for the different ions in Fig. 15a.

a Adsorbed fractions (N) of equivalents of Na+, Mg2+, K+ and Ca2+ on the complex, in equilibrium with seepage water (blue and orange) and unsaturated zone water (green). b The current fractions on the complex for the field samples at two depths for both sites (averaged over two tested sub-samples). NaX and CaX2 are both clearly in between the seepage and unsaturated zone equilibrium. Site A at 1.6 m depth is close to the unsaturated zone status. MgX2 shows higher levels, which is in line with preferential exchange of Na+ for Mg2+, followed by exchange of Mg2+ for Ca

Magnesium is far more abundant in the seepage water of site B than at site A (see section ‘Hydrochemistry models’), illustrated by Fig. 15b. Magnesium has increased compared to the initial situation, as explained in section ‘Transport and exchange in a drained field model’. However, the initial difference between sites remains up to this point. Table 2 compares the simulated ranges for shallow depth, coarsely schematized, of N. This is illustrated in Fig. 12 with the measured fractions of Na+, Ca2+, and Mg2+ at site B, where the sample at 1.6 m depth is from MP3 (see Fig. 2d) between drains, and the sample at 3 m depth is taken near a drain at MP5.

Table 2

Comparison of approximate simulated N value at 115-cm depth (upper three rows) and observed fractions (bottom two rows) on the complex of site B. Differences in depth are caused by limited data availability

CaX2

MgX2

NaX

Mid

Drain

Mid

Drain

Mid

Drain

Model start

0.2

0.2

0.3

0.25

0.5

0.5

Model after 100 years

0.5

0.4

0.4

0.35

0.25

0.3

Model after 200 years

0.8

0.6

0.2

0.25

0.0

0.1

Depth of observation (m)

1.6

3.0

1.6

3.0

1.6

3.0

Observed

0.55

0.48

0.35

0.5

0.09

0.03

Table 2 provides two clear suggestions concerning the relation between the situation at site B and the model. Firstly, although the depth of the simulated values presented in Table 2 (1.15 m) is shallow compared to the observation depths (1.6 and 3 m) the N values are comparable to the expected situation after approximately 120 years. This indicates that the field situation has developed over a longer time than the simulated time, which is in line with the work of Vos and Zeiler (2008) as cited in section ‘Study sites’. A very precise comparison in not possible since water management (dikes, polder levels) as well as hydrogeological conditions (sea level rise, land subsidence, e.g. Post et al. 2007) varied over time.

Secondly, the drain impact on water flow cannot already have influenced the composition of the exchange complex of soil. After all, these drains were installed only decades ago, which is little compared to the timescale of significant exchange variations according to the simulations. For the same reason, it may be understandable that even at the depth of 3 m, the composition of the exchange complex has changed towards a fresher situation then would be expected based on the simulations. Apparently, freshwater has been able to reach a much larger depth somewhere during the geological formation of the land.

Discussion and conclusions

In this paper, the cation exchange and transport in thin rainwater lenses on top of upward seepage of saline water was investigated in the case of an alternating flow pattern caused by recharge, evaporation, saline seepage and horizontal flow towards tile drains. Because a preliminary analysis of field data revealed some complexity, an abstracted 1D analog was first studied. After interpreting these results, it was possible to understand the cation exchange occurring in a more realistic 2D case that better represented field conditions.

The 1D models show the exchange that develops over time, starting from an initially saline situation with alternating inflow of calcium-enriched rainwater from the top and saline water from the bottom:

In the short term, a large change in solution composition is caused by the initial salt shock: freshwater enters the saline-water-filled domain.

In the longer term, gradual increases of calcium on the complex and in the solution are caused by the larger affinity of the complex for calcium in freshwater and the gradual supply of calcium by the freshwater.

Magnesium initially increases on the complex at the cost of sodium; however, as calcium in the solution gradually increases, it outcompetes magnesium on the complex.

Differences between the oscillation patterns of the 1D models are illustrated:

Oscillations of flow direction in combination with a zero net flux leads to a strongly oscillating pattern and, ultimately (i.e. after centuries) to a dispersed front with a stable alternation.

Oscillations with a 1 mm/day net downward flux behave similarly for an artificial step model and a natural net precipitation pattern. Both lead to a rather fast passage of the fresh/saline front up to just above the bottom of the column. Exchange progresses quickly compared to the zero net displacement situation, but decreases sharply in velocity with depth.

The field model with drains spaced at distances of 10 m, and constant upward seepage of saline water, can be interpreted using the 1D models:

Up to a depth below surface of ca. 1.2 m, there is a sequence and velocity of exchange comparable with the 1D model with natural net recharge

The limited downward movement of the front, inhibited by the constant upward seepage of saline water, causes the process at depths > ca. 1.2 m to resemble the 1D model with zero net flux.

The horizontal flow component is shown to have little influence on the exchange of saline and fresh ions, which can be understood from the fact that only vertical flow causes changes in pore water composition.

Comparison of field data with the model for site B shows that the main features of the field are represented:

In line with the findings of De Louw et al. (2011, 2013), the mixing zone position depends largely on the presence of aquitards combined with the position of the drain.

Comparing observations and model results, the exchange process seems to have progressed towards freshwater equilibrium for rather longer; even at larger depth (3 m) the exchange has progressed considerably.

After functioning for 30–40 years, the drains have only changed the composition of the pore water and not yet that of the soil complex.

Analysis of main cation concentrations in the soil of both study sites reveals:

Both for site A and B, observed complex composition confirms the expected increase of calcium, decrease of sodium, and intermediate increase of magnesium.

The initial difference in magnesium between the sites remains intact up to now, illustrating a long-term effect of the composition of the seepage water on the complex.

Both sites show that there is no equilibrium; however, the complex has developed a lot from the initially saline conditions towards a freshwater dominated complex

History for both study sites is known only approximately. After initial lens formation since AD 1300, intensified deeper drainage at the beginning of the 20th century was followed by tile drainage between 10 and 40 years ago. Given these big changes over the last century, and the ongoing changes in polder systems, both in water management (e.g. drainage level) and external changes (e.g. increasing seepage caused by relative sea level rise and climate change), a real equilibrium is unlikely. In this paper, cation exchange has been emphasized; however, in the case of large amounts of organic matter (e.g. peat layers), redox and chemical precipitation-dissolution reactions may complicate the hydrochemistry of the subsoil.

Notes

Acknowledgements

We thank the students that helped in the field and the laboratory gathering all necessary data, especially Bram de Vries, Marina Dudley-Southern and Julia Klaassen. We also thank Vincent Post for his help during the field campaign and his advice on the setup of the chemical model. This project was partly funded by the Dutch Technology Foundation STW and the Ministry I & M through the Program (14299) WaterNexus.

References

Alvarez MDP, Dapeña CE (2015) The role of evapotranspiration in the groundwater hydrochemistry of an arid coastal wetland (Península Valdés, Argentina). Sci Total Environ 506–507:299–307CrossRefGoogle Scholar

Eeman S, Leijnse A, Raats PAC, Van der Zee SEATM (2011) Analysis of the thickness of a fresh water lens and of the transition zone between this lens and upwelling saline water. Adv Water Resour 34:291–302CrossRefGoogle Scholar

Copyright information

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.