A practical quantification of error sources in regional-scale airborne groundwater salinity mapping

Our approach is summarized in six steps, and explained in order in sections 2.2 to 2.7: (1) synthetic data generation by transforming two existing 3D models (lithological and chloride) into an ECb volume; (2) simulate AEM acquisition using the ECb model (forward modelling) and subsequent inversion; (3) simulate lithological data collection from the 3D lithological model; (4) 3D interpolation of inversion results and lithological data; and (5) conversion to ECw and chloride by transforming interpolated results into ECw, and finally into 3D chloride models; (6) comparison with the synthetic chloride model and error analysis: finally, the estimated chloride estimates are compared against the synthetic model. Steps 2–6 are performed for different flightline spacings and different densities of ground-based information (geological drillings), where results are presented in section 3. The process is illustrated in figure S1 (stacks.iop.org/ERL/15/074002/mmedia). Hereafter, these steps are explained in more detail.

The synthetic data is based on two existing high-resolution 3D datasets from southwestern Netherlands (figure 1). The latitudinally oriented Zeeuws-Vlaanderen strip (Zeelandic Flanders in English) covers an area of ∼60 × 15 km and lies within the southern edge of the North Sea basin. The shallow subsurface (up to ∼50 m) consists of gently northward dipping Neogene and Quaternary sediments (Stafleu et al 2011). Overall, lithologies comprise younger fine sands, clays and peats and deeper, gently dipping sands and silts (Vos 2015). Along the North Sea coast sandy coastal barriers are marked by localized topographic highs. Holocene sea-level transgressions caused extensive groundwater salinization (Delsman et al 2014) however subsequent land reclamation has allowed the formation of shallow freshwater lenses (Berendsen 2005). Given these events are recent in geological terms, distinct zones of fresh and saline groundwater with a relatively thin brackish zone are present throughout (Delsman et al 2018).

Zoom In Zoom Out Reset image size

The chloride model used here (figure 1(b)) is the result of a provincial scale FEM survey that was flown in 2014/2015 in the Netherlands as a part of a research project called FRESHEM Zeeland (Delsman et al 2018). The resulting 3D groundwater salinity model is publicly available and has a resolution of 50 m horizontal and 0.5 m vertical, covering an area of ∼2000 km². Specific details of survey parameters are available in Delsman et al (2018), and summarised in table 2. AEM surveys decrease in resolution at depth, therefore it should be noted that the 0.5 m resolution throughout the model is the result of a secondary interpolation step. In the study area the survey was flown at mostly 300 m line spacing, in areas of interest spacing was decreased to 100 m resulting in ∼3000 line km of HEM data. Flightlines are shown in figure 1(a). Flightline spacing is selected according to a compromise between target sizes and budget. Tie-lines are flown perpendicularly to these for data processing at a wider spacing.

A high-resolution 3D lithological voxel model (100 × 100 × 0.5 m resolution, up to 50 m depth from ground surface) of the area is available from the GeoTOP geological model as part of the Dutch national database of geological data (dinoloket.nl—accessed November 2016), and provided by the Geological Survey of the Netherlands (GSN-TNO). As well as a stratigraphic model, GSN-TNO's GeoTOP model contains the probability of lithoclasses occurring at any voxel and is therefore a probabilistic heterogeneous style model (lithoclasses and formation factors outlined in table 1 and shown in figure 1(c)). A more detailed description of the model is available in (Stafleu et al 2011). Apparent FF values were assigned to these lithologies and are based on nearby field measurements (De Louw et al 2011) (table 1). We recognize that the EC of saturated sediment is complex, such as surface conductivity caused by clay particles (Revil and Glover 1998), however for this work, we neglect the effect of the clay to surface conductivity since the exact values are heavily depended on the area of investigation (Revil et al 2017). Previous research examines this effect (Delsman et al 2018), but for this regional analysis we did not consider it an important factor.

To forward model (and subsequently invert) an AEM survey, a synthetic ECb model (i.e. the combined effect of both groundwater and lithology) is required. Using a combination of the two existing 3D models, a realistic and heterogeneous-style ECb model was created, the result of which is shown in figure 1(d). This process was undertaken in two steps: (1) the 3D chloride model was transformed to ECw using an empirical linear relationship from field measurements (De Louw et al 2011) and corrected back to a reference groundwater temperature of 11 °C from an initial 25 °C; (2) ECw was transformed into ECb by applying FF values based on the geological model. The resulting regional-scale 3D ECb synthetic model has a cell resolution of 50 m horizontal and 0.5 m vertical, covering an extent of ∼15 × 60 km down to 50 m depth from ground surface.

To simulate a typical survey across the 3D model, a distribution of measurement locations was needed. These were taken from the same survey that produced the 3D salinity model (Delsman et al 2018) in order to replicate typical heliborne survey characteristics (figure 1(a)). In total ∼800 000 measurements locations over 3000 line km were used from 361 flightlines, these were typically flown at between 100–300 m flightline spacing depending on the area. To facilitate processing times, the measurements were decimated to one every ∼40 m from every ∼4 m (totalling ∼100 000 measurements)—this was found to be a good balance between respecting the FEM system's ∼100 m footprint (Yin et al 2014) and hardware practicalities. Given that the signal measured by the AEM system is inherently 3D, and the forward modelling code used in this is paper is 1D (Auken et al 2005)—a pseudo-3D sampling technique was implemented. To achieve this, 100 m diameter circles were created around each measurement location to approximate the footprint of the FEM system (Yin et al 2014). Within each circle ECb was averaged, this was repeated at 0.5 m depth-slice intervals down to 50 m depth (i.e. at vertical model resolution), resulting in 100 000 synthetically generated AEM acquisition points on flight lines, each with 100 ECb values with depth (every 0.5 m until 50 m). Figure 2 illustrates the pseudo-3D sampling technique used, as well as an illustration of the resulting AEM coverage based on different flightline spacing.

Zoom In Zoom Out Reset image size

The resulting data were forward modelled using AarhusINV (Auken et al 2005). Source and receiver spacing and frequencies were selected based on the RESOLVE HEM system; orientation was selected based on typical and previously used hydrogeophysical applications e.g. (Steuer et al 2009, Gunnink et al 2012, Delsman et al 2018). Finally, 5% random gaussian noise was added to the forward response to approximate HEM noise levels (Farquharson et al 2003, Green and Lane 2013).

The 'acquired' data, were inverted using AarhusINV (Auken et al 2005). Based on a commonly used inversion method (e.g. Auken et al 2008, Chongo et al 2015, Delsman et al 2018) and corresponding input parameters that were found suitable for hydrogeophysical studies (King et al 2018) (table 1), a laterally constrained inversion (LCI) was used (Auken and Christiansen 2004, Siemon et al 2009). The LCI method is a pseudo-2D method, whereby neighbouring model properties are constrained to produce laterally coherent inversion results along separate flightlines. The minimum-structure inversion style used here inverts for ECb values only, where layer thicknesses and depths remain constant, and has been found to accurately reproduce salinity distributions when compared with available ground data (King et al 2018, Delsman et al 2018).

The inversion starting model consisted of 20 layers, the thickness of the top layer was set to 0.68 m and subsequent layers increasing in thickness logarithmically until 50 m. Below 50 m the final layer is assumed to extend infinitely. A starting model conductivity of 1 S m⁻¹ was selected based on common knowledge that low-lying coastal zones are (electrically) conductive environments.

Simulating the acquisition of geological information, or indeed realistic quantities of existing data that may be available during hydrogeophysical mapping is not trivial. For example, the GeoTOP model used to construct the synthetic geological model in this paper is the result of expertly interpreted boreholes and geophysical data, resulting in a large-scale 3D probabilistic model with different properties (Stafleu et al 2011, Gunnink et al 2012). For this reason, a practical approach based on the current distribution of publicly available borehole data was used. Here we sampled synthetic borehole locations and depths from existing borehole data in the area available on Dinoloket, the Dutch Geological Survey (TNO) data portal. It was found that within the synthetic model area, cumulative vertical borehole depths totalled ∼42 000 m (or ∼60 vertical metres km⁻²), with a depth distribution highlighted in figure 3.

Zoom In Zoom Out Reset image size

Typically, most available field data was found to be very shallow (<5 m), with very few holes deeper than ∼40 m. To replicate this characteristic, Inverse Transform Sampling (ITS) was used which iteratively generated a set of borehole depth values that matched this distribution. This was undertaken for a set of prescribed total vertical borehole depths in six total length classes as follows: 100 m (0.13 m km⁻²), 1000 m (1.36 m km⁻²), 5000 m (6.80 m km⁻²), 10 000 m (13.61 m km⁻²), 20 000 m (27.23 m km⁻²) and 42 000 m (57.19 m km⁻²)—based on practical criteria up to the actual amount available. Five ITS realisations were used for each total length class to account for sampling bias and to quantify uncertainty due to placement location. Finally, lithological classifications were sampled from the synthetic voxel model by each randomly selected borehole at 0.5 m vertical resolution till borehole depth. In effect drilling into the model and undertaking a full lithological interpretation. This resulted in 6 × 5 = 30 synthetic borehole datasets containing FF information. Here it is assumed that the process of transforming lithological information into apparent FF is free of instrument or human error.

In order to assess the effect of different flightline spacing on mapping results, the ECb inversions were subsequently split into three separate databases: (1) 300 m spacing, 100 m in specific areas (∼3000 line km), (2) 600 m flightline spacing (∼1500 line km), and (3) 1200 m spacing (∼750 line km). The 100 m line-spaced data was kept together with the 300 m, as this reflects a realistic data distribution and is the actual distribution of the original survey. This grouping allows an analysis of the effects of a higher AEM data density in specific areas. Finally, each database was interpolated into 3D volumes of ECb using the method of King et al (2018) and is based on 2D Kriging. We recognise that there are many other sophisticated methods to interpolate data, however this approach was found to accurately preserve the characteristics of inversion results between flightlines and offered a straightforward and recognised approach. The model base was assigned according to the DOI of the airborne system which varied between 5 to 60 m between conductive and resistive respectively, this was calculated from the inversion output at each 1D model location (Vest Christiansen and Auken 2012). Minimum-curvature gridding was then used to interpolate DOI estimates, data below these were removed from further analysis for all further models. A statistically robust, simple and repeatable 3D interpolation method was utilised to interpolate the FF values at the borehole data throughout the synthetic model. For this an established technique called Sequential Indicator Simulation (SIS) was used (Gomez-Hernandez and Strivastave 1990), based on Indicator formulation. Separate indicator semivariogram models were fitted for all six lithoclasses from each borehole iteration (a separate semivariogram model with different horizontal and vertical ranges for each of the 30 borehole datasets). For all 30 lithological spatial distributions, 10 SIS realisations were simulated. The median (or p50) of the 10 simulated lithologies at each model cell was used for further analysis, resulting in six FF value classes. The interpolation resulted in 30 3D lithological spatial distributions for which appropriate FF values were assigned (for each of the six classes of specified total drilling length with five randomised locations for each).

Each of the 30 interpolated FF spatial distributions were multiplied by the three interpolated ECb inversion results to obtain 90 3D ECw spatial distributions (6 × 3 = 18 combinations of total borehole and flightline lengths with 5 random borehole sets per combination), which were then converted back to 3D chloride (as per figure S1).

With practical hydrogeophysical mapping considerations in mind, the resulting 90 3D ECw spatial distribution models were compared against the synthetic model in three ways. All evaluations were undertaken over the same domain, where values beneath the calculated DOI from inversion results were removed. First, mean absolute error (MAE) of groundwater salinity estimates as mg l⁻¹ chloride were calculated as per the following:

$$\begin{equation*}\mathrm{MAE}\, = \,\frac{1}{n}\mathop \sum \limits^{i = 1}_n \left| {{x_i} - x} \right|\end{equation*}$$

where x_i is the prediction and x is the true value. Second, groundwater salinity volumes were calculated according to the chloride classifications set out by Stuyfzand (1986) and include the following: (1) 0–1000 mg l⁻¹, or 'fresh', (2) 1000–3000 mg l⁻¹, or 'brackish', and (3) >3000 mg l⁻¹, or 'saline'. Finally, based on this classification, interface depths for the fresh-brackish-saline regions were extracted as 500, 1500 and 3000 mg l⁻¹ chloride respectively. Saline inversion was present in some areas (saline above fresh groundwater), therefore multiple interfaces were frequently encountered at the same horizontal location. As a result, only the shallowest, or first encountered interfaces were extracted as iso-surfaces for further analysis.

The mean absolute error (MAE) of mapped chloride is presented in figure 4, showing average error for all data configurations tested. To better isolate error per data type, in figure S2 we assumed a perfect inversion model and only adjusted the amount of lithological information.

Zoom In Zoom Out Reset image size

Figure 4 illustrates that increasing borehole density to a value of more than ∼10 m km⁻² results in relatively little improvement of mapped chloride. Geological data reduces mapping error by ∼100 mg l⁻¹ chloride between the highest and lowest data densities for all flightline configurations and at maximum 150 mg l⁻¹ chloride in case of perfect geological knowledge. Reducing flightline spacing consistently reduces error by ∼70 mg l⁻¹ chloride per 100 m reduction. The reduction in salinity estimation error by adding geological information is significant if the AEM inversion is without error (figure S2), the combined error is dominated by the error in AEM inversion. Thus, the most effective way of reducing the error in chloride concentration mapping is by decreasing flightline spacing.

Calculations were undertaken based on the classifications outlined in section 2.6, where a porosity factor of 0.3 was applied based on a sandy aquifer (Kirsch 2009), as the most commonly occurring sediment within the study area is fine sand. The result is a total groundwater volume of ∼5.71 × 10⁶ km³, with 3.13 × 10⁶ km³ (54%) fresh, 2.2 × 10⁶ km³ (39%) saline and 3.82 × 10⁵ km³ (6.7%) brackish groundwater volumes within the reference model. The AEM data acquisition and inversion process alone (assuming perfect geological knowledge) results in a substantial thickening of the brackish zone, adding ∼6.1 × 10⁵ km³ (10%) to brackish zone volume estimates (figure S4), underestimating the fresh and saline volume by ∼3.27 × 10⁵ km³ (4%) and ∼2.83 × 10⁵ km³ (6%) respectively. Figure 5 presents averaged volume differences against the reference model for each geological data density class with standard error (over the six borehole datasets per total borehole length class) and for a flightline distance of 300 m. Total error for each estimate has been subtracted by those of the reference model—thus positive values represent an overestimation.

Zoom In Zoom Out Reset image size

Results for 600 m and 1200 m flightline distances share a similar magnitude of error. Figure S5 shows that if no AEM acquisition and inversion error is included, the reduction of errors in estimated volumes are notable. With less geological data a trend of overestimating the brackish zone at the expense of saline volumes is seen. Quantitatively, the addition of geological data reduces volume estimate error between the 0.13 m km⁻² and 57.19 m km⁻² borehole classes by ∼8.0 × 10³ km³ (0.3%), ∼4.0 × 10⁴ km³ (10.5%), and ∼5.0 × 10⁴ km³ (2.3%) for the fresh, brackish and saline classes respectively. It is clear however in figure S4, that that the inversion results in a significant overestimation of the brackish zone, therefore the addition of geological information reduces this error, but only marginally relative to the error introduced by the inversion.

The analysis of errors in volume estimates show that the largest error comes from the AEM acquisition and inversion method. Contrary to errors in salinity concentration estimates, over regional scales this error is not very sensitive to flightline spacing due to data averaging. Laterally less extensive features however, such as a small (1–2 km wide) fresh-water lens are likely to be more sensitive to this. Also, the reduction of the volume error from additional geological data is insignificant. Thus, if one is interested in estimating volumes of fresh, saline and brackish groundwater by airborne AEM, large flightline distances and little geological information is the most efficient option with the current methods of AEM acquisition and inversion.

When resolving the 500, 1500 and 3000 mg l⁻¹ chloride interfaces, possible ECb values from inversion results may range between 0.06–0.13 S m⁻¹, 0.11–0.25 S m⁻¹ and 0.2–0.4 S m⁻¹ respectively, in case the FF is unknown and between the values illustrated in figure S6. The result is that there exists an interface depth uncertainty range of a specified vertical thickness as a result of the inversion method used. As volume estimates show (figure S4), there is a significant smoothing effect from the inversion process—thus it follows that this will affect the thickness of the FF uncertainty range. This effect is demonstrated in cross-section, where the uncertainty limits are presented in figure 6 for the reference model and for the 300 m flightline spacing inversion model.

Zoom In Zoom Out Reset image size

Figure 6 highlights how the vertical uncertainty range is narrow (∼1.7 m MAE) in the synthetic ECb model, with substantial thickening (∼2.9 m MAE) in the inverted ECb model. It is also apparent that the inversion has caused the interface to become shallower overall, which is the result of an overestimation of the brackish zone by the inversion. To understand magnitude of the contribution of flightline distance and borehole density on the error in interface estimation, the MAE of interface positions against the reference model for the 1500 mg l⁻¹ chloride are shown in figure 7. The 500 and 3000 mg l⁻¹ chloride plots are supplied as supplementary information (figures S7 and S8).

Zoom In Zoom Out Reset image size

Without the inversion, error margins are small (<0.5 m) for all interfaces and the addition of geology adds little improvement (between 0.05–0.1 m). The inversion and increasing flightline spacing however adds a more significant error to the mapping results: ∼3, 1.5 and 1 m for the 500, 1500 and 3000 mg l⁻¹ chloride interfaces respectively. For every 300 m reduction in flightline spacing, approximately 0.1, 0.25 and 0.4 m of accuracy are added for 500, 1500 and 3000 mg l⁻¹ chloride interfaces respectively—indicating that decreasing flightline spacing is useful for the more accurate delineation of deeper interfaces. However, the addition of geological information does little to help improve error and may even worsen it (for the 500 and 1500 mg l⁻¹ interfaces). This contra-intuitive effect may be caused by the relatively large number of shallower boreholes, which resulted in an oversampling of clay lithologies that are predominantly close to the surface. As clay is has a lower FF (2.5), this results in more brackish groundwater to be unjustifiably classified as fresh and enhances the shallow bias in the 500 and 1000 mg l⁻¹ interface depths. Naturally, this oversampling is more apparent when less lithological data are available.