Large-scale, probabilistic salinity mapping using airborne electromagnetics for groundwater management in Zeeland, the Netherlands

The global population is concentrated in coastal lowlands, owing to fertile soils and favourable transport connections (Aerts et al 2009, Nicholls and Small 2002). Sustainable use and protection of the available fresh water, however, present clear challenges in these coastal areas (Ferguson and Gleeson 2012, Michael et al 2017). The Dutch Province of Zeeland is a prime example. Groundwater resources in Zeeland have been largely salinised during marine transgressions, and freshwater is scarce. The available fresh groundwater is generally limited to elevated areas such as the dune areas and creek ridges. As a result of land reclamation, peat mining and subsidence (Pierik et al 2017), 25% of the area currently lies below mean sea level. The area is both an agricultural and tourism hotspot, putting stress on the available fresh groundwater resources. Successful management of these scarce resources requires a clear and detailed picture of their spatially variable occurrence.

From their first successful application in the late 1970s (Sengpiel and Meiser 1981), airborne electromagnetic induction (EM) surveys have become a popular method to map saltwater intrusion and delineate fresh groundwater occurrences (Paine and Minty 2005, Siemon 2006). Airborne EM is especially suited to these kind of surveys, as the electrical conductivity of the subsurface is determined by properties of (1) the porous medium, such as porosity and the electrical conductivity of the grains and (2) the pore water, such as its volumetric content and its salinity (Kirsch 2006). The airborne EM systems currently in use have been shown to yield good results when compared to ground-based geophysical and borehole measurements (Chongo et al 2015, De Louw et al 2011, Fitterman and Deszcz-Pan 1998, Pedersen et al 2017, Rasmussen et al 2013 Steuer et al 2009, Sulzbacher et al 2012).

While airborne EM provides a detailed image of the conductivity of the subsurface, this is seldom sufficient for water management purposes. Rather, the interest of water users and managers is focused on groundwater salinity. In the Netherlands, chloride concentration is historically used as a proxy for water salinity, as chloride is the major conservative anion found in saline groundwater bodies in the Dutch subsurface, being mixtures of seawater and freshwater. Only few studies have attempted to quantify groundwater salinity from airborne EM surveys. Paine (2003) combined airborne EM data with borehole measurements to estimate chloride mass in a saline-water plume. Faneca Sànchez et al (2012) use Archie's formation factors (Archie 1942) to translate airborne EM data to groundwater salinity as input for a three-dimensional groundwater modelling study.

The uncertainties associated with the inversion of airborne EM measurements have been widely noted, as the inversion problem is mathematically ill-posed and prone to non-uniqueness (King et al 2018, Sengpiel and Siemon 2000, Siemon et al 2009b, Viezzoli et al 2009). Researchers have worked to regularise the inversion through introduction of a-priori data, or by applying lateral (i.e. along a flight line) (Auken et al 2005, Siemon et al 2009a) or spatial (i.e. along and across flight lines) (Steuer 2008, Viezzoli et al 2009) constraints. Recent efforts also focused on joint-inversion of geophysical and hydrological data (Christensen et al 2017, Herckenrath et al 2013). While these approaches have successfully constrained and improved the inversion of airborne EM measurements, non-uniqueness will persist and analysis of parameter uncertainty remains important (Minsley 2011). Stochastic approaches to the inversion of hydrogeophysical data, allowing a more thorough inclusion and analysis of uncertainty, are therefore becoming more widespread (e.g. Chen et al 2008, Laloy et al 2012). In ground based geophysics, Jardani et al (2013) present a joint stochastic inversion of direct-current resistivity, self-potential data and well salinity measurements. Minsley (2011) presents a Bayesian approach to inversion of airborne EM data, but, perhaps due to the required computational load, the approach has not been widely adopted. A Monte-Carlo based uncertainty propagation analysis could provide a relatively straightforward and computationally manageable way to characterise uncertainty. Moreover, such an approach not only allows characterisation of the uncertainty associated with the inversion of airborne EM data, but also the subsequent translation to groundwater salinity.

In this paper, we present the approach and results of the FRESHEM Zeeland study: a 3D probabilistic mapping of groundwater salinity based on a large-scale (1800 km², 9640 flight line kilometres) airborne EM survey. The novelty of this work lies in the detailed translation of bulk ground conductivity to groundwater salinity, using available high-resolution lithological information. Moreover, our approach includes a Monte-Carlo based probabilistic uncertainty propagation analysis to quantify the uncertainty associated with the different steps in our approach. A locally varying anisotropic kriging approach was applied to preserve local scale fresh groundwater features in the 3D result. Survey results provide an unprecedented view of groundwater salinity in the Dutch province of Zeeland, including previously unknown features of the Zeeland groundwater salinity distribution such as fresh groundwater extending out into the sea, and significant fresh groundwater resources below salinised groundwater. FRESHEM results are used by local stakeholders to support the management of fresh groundwater.

2.1. General methodology and uncertainty assessment

Our methodology comprises the following steps (figure 1): (1) airborne EM survey, including data processing, (2) 1D inversion of airborne EM data to bulk conductivity, (3) separation of bulk conductivity in contributions of groundwater salinity and host material, (4) transformation of groundwater conductivity to groundwater chloride concentration, (5) 3D interpolation of chloride concentration. Steps 3 to 5 are embedded within an uncertainty assessment framework.

Zoom In Zoom Out Reset image size

For the calculation from bulk conductivity resulting from the 1D inversion (result step 2) to the groundwater chloride concentration (result step 4) a Monte Carlo approach was used. In this approach the calculation steps were performed 600 times for each data point, using parameter-values drawn from their distribution (see figure 1 for list of stochastic parameters). This resulted in 600 values of chloride concentration for each data point, schematised in a probability density function (pdf) for use in the 3D interpolation (step 5).

2.2. Airborne EM survey and data processing (step 1)

The survey was performed using the frequency-domain, helicopter-borne electromagnetic (HEM) system (RESOLVE), operated by the airborne geophysics group of the German Federal Institute for Geosciences and Natural Resources (BGR). The system consists of five horizontal-coplanar (HCP, 8 m coil separation) coil systems and one vertical-coaxial (VCX, 9 m coil separation) coil system; measuring frequencies range from 380 Hz to 130 kHz. The survey covered the entire province of Zeeland, the Netherlands (total area 1800 km², 9640 flight line km) and was conducted in three two-week periods in 2014 and 2015, in addition to a one-week period in 2009. We consider the time between these surveys allowable, given the relative stationarity of the groundwater salinity distribution (De Louw et al 2011). Flight line separation was generally 300 m, some high-interest areas were surveyed with flight line separations of 100 or 200 m (figure 2). Average sensor altitude was about 50 m above ground level. The amount of infrastructure present necessitated frequent deviations from the planned flight paths; power lines, railway tracks and major roads were crossed at higher elevation. Raw measurement data was carefully processed to minimise influences not related to subsurface conductivity, such as thermal drift (Siemon 2009) and cultural noise (Siemon et al 2011). Data of inadequate quality, mostly caused by infrastructure (power lines, radar stations, major roads), were not used in further calculation of the chloride concentration.

Zoom In Zoom Out Reset image size

Table 1. Applied inversion approaches, used parameters, software and references.

Inversion	Description	Software	Reference
SSI	Single site Marquardt-Levenberg inversion, 19 layers with fixed thicknesses (except cover layer) on top of homogeneous half-space. Starting model from apparent conductivity/centroid depth sounding curve, vertical smoothing achieved by regularisation.	BGR proprietary	(Sengpiel and Siemon 2000, Siemon 2012, 2006)
LCIsmooth	Laterally constrained inversion, 20 layer smooth model (L2 norm inversion scheme). The regularization penalizes vertical changes in conductivity, resulting in a vertical smooth conductivity model. This inversion scheme is considered a 1.5D inversion scheme. Forward modelling is still 1D, but the inversion is 2D. This way, nearby inversions cells are processed simultaneously and lateral weights during the inversion are added, in such a way that the model consists of a 2D section.	Aarhus Workbench	(Auken et al 2005, 2009)
LCIsharp	Laterally constrained inversion, 20 layer sharp model. Different to LCIsmooth as the regularization penalizes the number of vertical conductivity transitions of a certain size, resulting in conductivity models with relatively sharp vertical conductivity transitions.	Aarhus Workbench	(Auken et al 2005, 2009)

2.3. Inversion of airborne EM data (step 2)

2.4. Derivation of groundwater conductivity (step 3)

In deltaic environments, sedimentary composition is not clear-cut, and different sediments show significant surface conductivities (electrical conductivity via the electrical double layer surrounding the mineral grains) (Revil et al 2017). We calculated groundwater conductivity following Waxman and Smits (1968):

$$\begin{equation} \sigma _w = F(\sigma _b - \sigma _s ), \end{equation} \tag{ 1 }$$

with σ_w being groundwater conductivity (mS cm⁻¹), F formation factor, σ_b bulk conductivity (mS cm⁻¹) and σ_s surface conductivity (mS cm⁻¹). The relation between σ_b and σ_w becomes non-linear at low groundwater conductivities, and applying equation (1) therefore underestimates σ_w for low values of σ_w (Revil et al 2017). However, as further differentiation within the freshwater category was not an objective of this study, we applied the linear equation (1) over the entire conductivity range.

Lithology for the entire study area was derived from GeoTOP (Stafleu et al 2011), a 3D voxel (100 × 100 × 0.5 m) model of stratigraphy and lithology of Zeeland to a depth of about 50 m below mean sea level (MSL), developed by the Geological Survey of the Netherlands. GeoTOP is a stochastic model based on interpolation of 23000 stratigraphically interpreted corings in the study area (about 13 corings per km²) (Stafleu et al 2011). For every voxel the most likely stratigraphy (following Weerts et al 2005) is reported, as well as the probability of the voxel being best represented by each of the following lithological classes: peat, clay, clayey sand or sandy clay, fine sand, medium sand, coarse sand, gravel or shells. Information on lithology below 50 m below MSL was obtained from the layer-based hydrogeological model REGIS-II (Van der Meulen et al 2013). REGIS-II is a nationwide model, indicating depth and thickness of sandy aquifers and clayey aquitards on a 100 × 100 m grid. Note that we accept these high-resolution lithological models—with their uncertainty—as truth for this study; no attempt was made to use bulk conductivity results to improve the lithological information (see for instance Gunnink et al (2012)).

We characterised parameters F and σ_s for each lithological class in GeoTOP using laboratory analysis of 71 samples obtained during a sampling campaign in May 2016 in Zeeland (Revil et al 2017). Each of these samples was equilibrated with six sodium chloride solutions with known σ_w (ranging from 0.31–220 mS cm⁻¹ (at 25 °C)). The σ_b for each solution was measured with a high-precision impedance analyser using the four-terminal method of Zimmermann et al (2008). Equation (1) was then fitted to the data to obtain sample-specific values for F and σ_s. To derive values for F and σ_s per lithological class, we then classified the lithology of each sample, and fitted equation (1) to the combined samples of a given lithological class. The standard deviation of the regression was used in the Monte Carlo procedure as uncertainty. While differences in F and σ_s are expected given the sometimes distinct differences in depositional environment, data was insufficient to further discriminate between lithologies of different stratigraphical units.

2.5. Transformation to groundwater chloride concentration (step 4)

We transformed groundwater conductivity to chloride concentration using the linear relation established by (De Louw et al 2011), based on 79 Zeeland groundwater samples:

$$\begin{equation} Cl_w = \alpha \sigma _w - \beta , \end{equation} \tag{ 2 }$$

with Cl_w chloride concentration (in mg Cl L⁻¹), σ_w groundwater conductivity referenced to 25 °C (mS cm⁻¹), α and β regression parameters with values of 360 (σ = 6) and 450 (σ = 190) respectively (for the Monte Carlo approach a normal distribution of these parameters is assumed). Equation (2) is valid from a σ_w of about 2 mS cm⁻¹, above which the linearly related Cl⁻ and Na⁺ ions dominate the electrical conductivity. We applied this relation for the entire σ_w range after linearly correcting σ_w for temperature, assuming a constant groundwater temperature of 11 °C (annual average). We finally classified chloride concentrations in ten discrete classes, to account for underestimations due to the use of equations (1) and (2) over the entire σ_w range, to facilitate interpolation (see section 2.6) and to avoid overly precise use of the final data.

2.6. Interpolation to 3D model (step 5)

Resulting data of chloride probabilities at depth intervals at HEM measurement points were interpolated using ISATIS© geostatistical software (Bleines et al 2004) to obtain a 3D voxel model (50 × 50 × 0.5 m) of groundwater chloride concentration. The interpolation approach sought to best deal with (1) the bimodality of the chloride data, with high frequencies of both fresh and saline concentrations, and only limited occurrence of brackish concentrations, (2) non-random spatial distribution of values, located along flight lines, (3) the occurrence of small, elongated fresh groundwater bodies below sandy creek ridges (De Louw et al 2011, Pauw et al 2015), and (4) the interpolated quantity being probabilities rather than discrete values.

Chloride probabilities were first resampled to an average probability density function of chloride concentration per voxel. To minimise 'smoothing' of the bimodal chloride distribution, the chloride pdf was discretized to a cumulative density function at ten chosen threshold values (being the ten chloride class boundaries of section 2.5), enabling interpolation to the 3D model by indicator kriging (Goovaerts 1997). Fresh groundwater bodies underneath narrow elongated creek ridges give rise to clear local anisotropy in the spatial pattern of groundwater salinity. Along these creek ridges, interpolation would result in high uncertainties if this local anisotropy is not accounted for. We therefore used a locally varying anisotropy (LVA) approach (te Stroet and Snepvangers 2005) and estimated the direction of anisotropy and the corresponding spatial correlation from the salinity data along the flight lines (figure 3). The angle of the anisotropy ellipse was determined from the largest distance from each 'fresh' voxel (chloride concentration below a threshold of 1000 mg L⁻¹) to a 'brackish' voxel (above a threshold of 3000 mg L⁻¹). Length of the short axis was determined by the smallest distance to a 'brackish' voxel, length of the long axis was fixed at 1000 m. Anisotropy parameters were then interpolated using an inverse distance weighted technique to obtain an LVA field for the full extent of the 3D model (figure 3). The LVA field defines both the anisotropy of the neighbourhood selection as well as the variogram anisotropy used in the final kriging interpolation.

Zoom In Zoom Out Reset image size

2.7. Uncertainty contributions and validation

We assessed the contribution of each stochastic variable (inversion scheme, lithological class, formation factor, surface conductivity, and conductivity to chloride translation) to the total uncertainty present in the flight line results (classified chloride concentration) before interpolation. We defined uncertainty as the mean range between the p10 (10th) and p90 (90th percentile) values, expressed in number of chloride class difference. In a representative subarea (see symbol U in figure 2), we sequentially performed steps 2–4 of the outlined procedure (figure 1), while for each repetition only one of the stochastic variables was allowed to vary. Separate (marginal) uncertainties were then compared to the overall uncertainty. We calculated the contribution of the 3D interpolation (step 5) to the uncertainty of the end result by comparing total uncertainty in this 3D result versus the uncertainty in the result before interpolation.

We validated the end result using (1) groundwater chloride analyses for a total of 1294 monitoring well screens, and (2) depth of the fresh-saline interface for a total of 178 interpreted geophysical measurements (borehole logging and electrical cone penetration tests) (locations in figure 4). For (1), we compared classified measured chloride concentrations to the chloride class of the 3D result averaged over the well screen depth, using the classes 'fresh' (< 1500 mg Cl L⁻¹), 'brackish' (1500–10000 mg Cl L⁻¹) and 'saline' (> 10000 mg Cl L⁻¹). For (2), we compared the top and the bottom of the visually interpreted fresh-saline transition zone in the geophysical measurement to the top and the bottom of the visually interpreted fresh-saline interface in the 3D result. All available historical measurements were included in the validation, given the relative stationarity of the groundwater salinity distribution, except at selected locations where significant changes in the groundwater system were expected due to known historical events (e.g. discontinuation of groundwater extraction, inundation of polders). Here, measurements prior to the event were discarded.

Zoom In Zoom Out Reset image size

The constructed voxel model of chloride concentration probability covers an area of about 1800 km², with vertical extent ranging from 30–160 m below ground surface (BGS). Figure 5(a) shows a slice of the voxel model, for p50 (median) chloride probability. The result shows clear and plausible patterns of groundwater salinity, with deep freshwater lenses in the dune area and at the slightly elevated creek ridge areas (De Louw et al 2011), and shallow occurrence of saline groundwater in low-lying areas. The sharp spatial delineation of fresh groundwater features is apparent from a derived vertical 1500 mg Cl L⁻¹ contour map (figure 5(b)). Our investigation shows the limited occurrence of fresh groundwater in Zeeland, as 75% of the area (70%–80% p25–p75 range) lacks fresh groundwater below 15 m BGS. Even though investigation depth in airborne surveys is generally limited when shallow saltwater is encountered (Viezzoli et al 2010), our survey indicated the presence of fresh groundwater below saline water in a number of areas across Zeeland (figure 5(c)). This result was corroborated with limited available deep bore logs. Widespread occurrences of fresh groundwater below saline groundwater in the south of Zeeland (indicated by I in figure 5(c)) coincide with the presence of mapped impermeable clay layers (based on REGIS (van der Meulen et al 2013)), which hindered the infiltration of saline water during marine transgressions. Also notable is the occurrence of fresh groundwater extending outward into the sea, driven by the head gradient from the nearby dune area (II in figures 5(a) and (c)). We found uncertainty in the 3D chloride distribution to be most pronounced in no data areas (power lines, major roads, build-up areas), and the lateral edges of freshwater lenses, where chloride gradients are steep (figure 5(d)).

Zoom In Zoom Out Reset image size

Figure 6 shows intermediate results of p50 chloride concentration, bulk conductivity and lithology along the flight line indicated in figure 5(b). The profile shows a representative cross-section of Zeeland that extends from the dune area, with fresh groundwater resources present to a depth of around 40 m MSL, to the low-lying hinterland, with shallow fresh groundwater resources limited to narrow creek ridges (De Louw et al 2011). The delineation of the transition zone between fresh and saline groundwater is sharper in the p50 chloride concentration than in the bulk conductivity, owing partly to the effect of lithology in the bulk conductivity profile. The presence of clay layers, e.g. below the dune area, is absent from the chloride concentration profile, as their influence on the bulk conductivity is successfully accounted for in the FRESHEM approach. Note that, conversely, bulk conductivity variations within regions of known equal salinity could be used to delineate lithological variations and, for instance, validate the lithological prediction of the GeoTOP model. This was, however, outside the scope of this research. Available groundwater chloride analyses and ground-based geophysical results are also plotted in figure 6, FRESHEM results show good visible correspondence with these measurements.

Zoom In Zoom Out Reset image size

Validation results for chloride analyses show that median results corresponded for 78% with classified chloride analyses (93% for p10–p90 range). Differences between measures and modelled chloride values were most pronounced for the intermediate, 'brackish', concentration class. This finding is consistent with the brackish zone occupying only a small interval between fresh and saline groundwater: small errors or a shift in the vertical location of the interface will therefore promptly cause errors in validation. Furthermore, most differences occurred at shallow depths (< 5 m BGS), corresponding to the larger spatial variation in groundwater salinity at shallow depths. No spatial patterns in error distribution were observed. On average, the depth of the fresh-saline interface in the results was predicted about half a meter deeper compared to local geophysical measurements (mean error (ME): 0.56 m, root mean squared error (RMSE): 1.5 m).

The contributions of the different stochastic variables to the overall uncertainty was analysed by comparing the mean difference in chloride class between the p10 and p90 (10th and 90th percentile) for each stochastic variable separately. Results are listed in table 2. Note that the values in table 2 for the stochastic variables separately do not add up to the uncertainty when all variables are combined due to interdependencies between variables. Largest uncertainty is caused by the three inversion methods, then surface conductivity and the lithological class. The added uncertainty from the 3D interpolation amounts to about 38% (calculated as (1.70–1.06)/1.70) of the total uncertainty in the 3D voxel model.

All results have been made publicly available through a user-friendly website of the Province of Zeeland (available at https://kaarten.zeeland.nl/map/freshem), and the results have been taken into use by relevant stakeholders. The 3D groundwater chloride distribution has so far been used: (1) to guide farmers in their use of fresh groundwater resources, (2) as a starting concentration in the groundwater model of Zeeland (van Baaren et al 2016), (3) as a new basis for groundwater extraction zoning, and (4) as input for feasibility maps of different agricultural aquifer storage and recovery measures. Interestingly, stakeholders were well able to handle the uncertainty in the results, after providing guidance and deriving tangible uncertainty measures.

In this paper we present the construction of a high-resolution 3D voxel model of groundwater chloride concentration from a large-scale airborne EM survey in the Province of Zeeland, the Netherlands. We applied a forward uncertainty analysis and propagated uncertainties associated with the different steps in our methodology. We validated our result and found good correspondence with available chloride concentration measurements and ground-based geophysical measurements. The results help stakeholders to sustainably manage the scarce available fresh groundwater resources in the region.

This study highlights the importance of applying uncertainty analysis in airborne EM surveys, given the demonstrated considerable uncertainty in our results. Uncertainty mainly originated from the inversion and the 3D interpolation. Uncertainty was largest at transition zones between fresh and saline groundwater, both laterally and vertically, while generally low in fresh or saline zones. This indicates the ability of airborne EM surveys to effectively map the occurrence of fresh groundwater resources in estuarine deltaic areas like Zeeland, while the delineation of the exact transition between fresh and saline groundwater is more challenging. This finding was corroborated by validation of the results, where deviations between the FRESHEM result and available chloride measurements occurred mainly in intermediate salinities. Note, however, that results of our uncertainty analysis are conditional on our prior assumptions, i.e. our chosen probabilistic parameters and the probabilities we assigned to each of them. Uncertainties in the raw airborne EM measurements were for instance not included, nor will the three included inversion methods span the entire plausible range of inversion results. Moving from uncertainty propagation towards conditioning uncertainties on available ground measurements, perhaps building on the work of Minsley (2011), should improve the quantification of uncertainty.

The significant surface conductivities measured in the considered sediments (Revil et al 2017) necessitates accounting for the conductivity of the host material in this and similar delta settings. The availability of a high-resolution, probabilistic lithological model allowed us to translate bulk ground resistivities into groundwater salinities, while accounting for the conductance of the host material. This approach proved able to correct bulk conductivity for lithological differences, and to remove most lithological 'artefacts' from the chloride concentration result. Still, the translation is far from straightforward, and requires a pedotransfer function relating conductive properties to a lithological class. We used the dataset reported by Revil et al (2017) as a basis for such a pedotransfer function. Differences in conductive properties within the lithological classes were, however, high, reflecting their relatively broad definition (Weerts et al 2005). A need therefore exists for improving pedotransfer functions for the conductive properties of sediments; the ability to relate conductive properties to more easily measured quantities such as the cation-exchange capacity may provide a way forward (Revil et al 2017). Moreover, applying such relationships in more data-scarce regions is still challenging and requires future work, e.g. on approaches that aim to infer lithological properties from the airborne EM signal (Gunnink et al 2012).

The availability of a large-scale 3D chloride concentration voxel model increased the understanding of the groundwater chloride distribution in Zeeland, and of the processes governing this distribution. Fresh groundwater resources underlying saline groundwater are, for example, more widespread than previously known; their presence under aquitards points to the role of these aquitards in preventing the infiltration of saline water during past marine transgressions. The translation of bulk conductivity, the common end-result of airborne EM surveys (e.g. Kirkegaard et al 2011, Pedersen et al 2017), to groundwater chloride concentration, the quantity stakeholders are interested in, likely improved uptake of the results by stakeholders. Results now guide the use of scarce fresh groundwater resources in the area, as well as the implementation of recently developed managed aquifer recharge approaches (Pauw et al 2015, Zuurbier et al 2015), and overall increase awareness on freshwater scarcity in the region.

This work was funded by Provincie Zeeland, Waterschap Scheldestromen, Evides Waterbedrijf, Rijkswaterstaat Zee and Delta, Vereniging van Zeeuwse Gemeenten, Vlaams Nederlandse Scheldecommissie, Zuidelijke Land en Tuinbouw Organisatie and the Deltafonds. This paper benefited from the helpful suggestions of two anonymous reviewers.