Characterisation of subglacial water using a constrained transdimensional Bayesian transient electromagnetic inversion

Subglacial water modulates glacier-bed friction and therefore is of fundamental importance when characterising the dynamics of ice masses. The state of subglacial pore water, whether liquid or frozen, is associated with differences in electrical resistivity that span several orders of magnitude; hence, liquid water can be inferred from electrical resistivity depth profiles. Such profiles can be obtained from inversions of transient (time-domain) electromagnetic (TEM) soundings, but these are often non-unique. Here, we adapt an existing Bayesian transdimensional algorithm (Multimodal Layered Transdimensional Inversion – MuLTI) to the inversion of TEM data using independent depth constraints to provide statistical properties and uncertainty analysis of the resistivity profile with depth. The method was applied to ground-based TEM data acquired on the terminus of the Norwegian glacier, Midtdalsbreen, with depth constraints provided by co-located ground-penetrating radar data. Our inversion shows that the glacier bed is directly underlain by material of resistivity 10² Ωm ± 1000 %, with thickness 5–40 m, in turn underlain by a highly conductive basement (10⁰ Ωm ± 15 %). High-resistivity material, 5×10⁴ Ωm ± 25 %, exists at the front of the glacier. All uncertainties are defined by the interquartile range of the posterior resistivity distribution. Combining these resistivity profiles with those from co-located seismic shear-wave velocity inversions to further reduce ambiguity in the hydrogeological interpretation of the subsurface, we propose a new 3-D interpretation in which the Midtdalsbreen subglacial material is partitioned into partially frozen sediment, frozen sediment/permafrost and weathered/fractured bedrock with saline water.

1 Introduction

Subglacial structure and material properties are one of several important controls on ice flow, both through composition and ice–material interactions. The potential for subglacial sediments to store and pressurise water is a key element in predicting the evolution of ice masses of all sizes, from small mountain glaciers to large polar ice sheets (Christoffersen et al., 2014; Siegert et al., 2018). Currently, the ability to develop accurate ice-flow models is limited by poor understanding of processes acting at the ice–bed interface and the composition of subglacial material. With increased knowledge of subglacial structure and sediment liquid water content, our ability to predict glacier retreat patterns would be greatly increased.

Non-invasive geophysical imaging methods are widely and successfully applied to characterise the internal properties of glacier ice and its immediate basal environment. Such methods (including reflection seismology and ground-penetrating radar) can underperform when characterising material properties beyond the first few metres of the glacier bed (Booth et al., 2012), yet subglacial aquifers, sediment accumulations and permafrost can extend to much greater depths (e.g. Mikucki et al., 2015; Hauck et al., 2001). Further still, inversions of isolated geophysical datasets are unconstrained and non-unique, with many models of the subsurface matching the observed dataset. Joint inversions using multiple independent datasets can constrain the model space, combining depth and resolution sensitivities from multiple datasets. Glaciological surveys often involve the acquisition of multiple geophysical datasets: given the typical absence of ground-truth data, imaging the target with several methods provides a more robust interpretation (e.g. Merz et al., 2016). However, these datasets are seldom combined numerically. In this paper, we provide a mechanism for the constrained inversion of transient electromagnetics (TEM), with depth constraints derived from ground-penetrating radar (GPR), to provide geophysical insight into the structure and water characteristics of the subglacial environment. This method is adapted from a transdimensional Bayesian framework, termed “MuLTI” (Multimodal Layered Transdimensional Inversion) and described in Killingbeck et al. (2018), originally applied to characterise subglacial sediment distribution from seismic surface wave data (Killingbeck et al., 2019). Here, we explore a similar concept for the TEM method.

Time-domain electromagnetic methods use electromagnetic fields to sound the subsurface structure. Here, we use the transient time-domain method (TEM) which indirectly probes the subsurface resistivity structure by measuring transient eddy currents induced by current transmitted through either a grounded wire, coincided loop or offset coil. The method has a depth sensitivity ranging from a few metres to kilometres, depending on the survey parameters used. Of particular relevance here is that electrical resistivity increases by several orders of magnitude when water in pores freezes (Hoekstra and McNeill, 1973), allowing resistivity methods to indicate the liquid water content of subsurface materials. TEM methods have been extensively applied for hydrogeophysical exploration to map groundwater resources (Auken et al., 2003), mapping permafrost on mountainous regions under debris-covered glaciers (Hauck et al., 2001), mapping arctic permafrost in Alaska (Minsley et al., 2012) and more recently mapping deep saline groundwater zones in Antarctica's Taylor Valley (Mikucki et al., 2015). These studies illustrate that characterising the resistivity of the subsurface offers a promising means of distinguishing material type and water content within the subglacial environment.

In common with most geophysical inversions, such resistivity profiles are non-unique: many profiles fit the data within error tolerance, and smoothing is usually employed to recover a single solution. Early inversion techniques for TEM data included non-linear least squares (Barnett, 1984) and an Occam-type regularisation method to obtain a smooth solution (Constable et al., 1987), but these were prone to being trapped in local minima with any large resistivity variations becoming smoothed. More recent inversion methods include laterally and spatially constrained algorithms to regularise the inversion and obtain solutions that agree with the expected geological variations (e.g. Christensen and Tølbøll, 2009; Vignoli et al., 2015; Auken et al., 2015). Yet, these methods do not provide detailed uncertainty analysis of the estimated model parameters and require a fixed number of layers in the model. The maximum depth of investigation (DOI) is generally estimated using methods such as half-space skin depth (Spies, 1989) or the Jacobian sensitivity matrix (Christiansen and Auken, 2012), though these do not consider the non-linear sensitivity of the DOI to conductivity structure. These limitations in uncertainty quantification, fixed model space and DOI estimation can be mitigated by transdimensional Bayesian sampling-based inverse methods. These produce an ensemble of models from which statistical properties of the model parameters, including model dimensions, can be inferred (Mosegaard and Tarantola, 1995; Blatter et al., 2018). The computed posterior probability density function (PDF) provides a robust measure of DOI, highlighting model uncertainty at each depth (Blatter et al., 2018). To further reduce the parameter space and improve vertical resolution, the inversion can be constrained with complementary depth information (e.g. from borehole records or other geophysical sources), which is most useful in the event that any internal layer represents a discontinuity in properties to which all techniques are sensitive.

In this paper, we derive the implementation of “MuLTI-TEM” (Multimodal Layered Transdimensional Inversion of Transient ElectroMagnetics) and its use for characterising the subglacial environment. After testing the method on a synthetic dataset, we analyse a TEM dataset acquired on the Norwegian glacier, Midtdalsbreen, an outlet of the Hardangerjøkulen ice cap, using complementary GPR data for constraining the ice thickness. Since the glacier bed represents a transition in both dielectric constant and electrical resistivity, the GPR depth constraint can be used directly in the TEM inversion. Recent results from Killingbeck et al. (2019) interpret the Midtdalsbreen subsurface from seismic shear-wave velocity (V_s) profiles as local accumulations of soft material (partially frozen glacial till) and permafrost overlying bedrock. Finally, we show a combined 3-D interpretation of previous shear-wave inversions (Killingbeck et al., 2019) and results output from MuLTI-TEM, and suggest the development of a joint resistivity–V_s depth-constrained inversion strategy.

2 Method

2.1 Transient electromagnetics

In TEM surveying, an electromagnetic field is generated by sending a periodic, modified square-wave current through a transmitter coil. When the current is on, a static electromagnetic field is established in the ground. The electromagnetic field is varied by terminating the current abruptly at the first quarter period, being reduced to zero for the second quarter period; the current is then reversed for the third quarter period, before being reduced to zero again for the final quarter period. This switch off induces eddy currents in the subsurface, initiating within the immediate vicinity of the transmitter, then spreading downwards. The eddy currents produce a secondary electromagnetic field which propagates back up through the subsurface, inducing a current in a receiver coil located at some distance from the transmitter. The receiver typically measures the induced secondary electromagnetic field in the transmitter off periods. The response of the subsurface is measured in terms of the decaying amplitude of the secondary electromagnetic field. This is recorded as a function of time, with later responses originating from greater depths. With regards to conductivity of the subsurface, the more conductive the subsurface, the larger the eddy currents and the larger the measured secondary electromagnetic field will be, implying a slower transient decay. By taking repeated measurements, a sounding curve, similar to DC resistivity soundings, is obtained (Geonics, 1994). The measured voltages versus time from the receiver coil are then used to constrain the resistivity profile with depth.

The maximum depth of investigation (DOI) (h), in metres, can be approximated by

$$\begin{array}{}\text{(1)}& h\approx \mathrm{0.55}{\left({\displaystyle \frac{IA\mathit{\rho }}{{\mathit{\eta }}_{\mathrm{v}}}}\right)}^{\mathrm{1}/\mathrm{5}},\end{array}$$

where I is the transmitter current, A is the loop area, η_v is the voltage noise level, and ρ is the average resistivity, in Ωm, of the underlying section (Spies, 1989). Equation (1) shows that the transmitter loop size is an important acquisition parameter controlling depth of investigation. Large loop soundings (e.g. $>\mathrm{40}\times \mathrm{40}$ m), where the receiver coil is located in the centre of the large transmitter coil, have been conducted on thick permafrost regions by Rozenberg et al. (1985) and Todd and Dallimore (1998). Our depth of target (0–80 m) allows for a smaller, more portable loop size to be used. In this study, we use an efficient 10×10 m loop, with the receiver 15 m offset from the centre of the transmitter loop (to stop electromagnetic interference with the receiver) shown and discussed further in Sect. 3.

2.2 MuLTI-TEM

MuLTI-TEM is a Bayesian inversion MATLAB code that determines the posterior distribution of resistivity as a function of depth. It is adapted from the MuLTI algorithm, developed for seismic surface wave inversions (see detailed information in Killingbeck et al., 2018). The data input, d, to MuLTI-TEM is the voltages (v) at each of the N time gates, measured as the mean recording in a stack window. The mean, through central limiting, is assumed to be normally distributed, with a variance σ², the variance of the measurements divided by the number of measurements in the stack window, so that the data and uncertainties can be written as

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \mathit{d}=[v_{\mathrm{1}},v_{\mathrm{2}}\mathrm{\dots },v_N],\\ & \mathit{\sigma }=[{\mathit{\sigma }}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{2}}\mathrm{\dots },{\mathit{\sigma }}_N].\end{array}\end{array}$$

The method used to find the posterior distribution of the resistivity profile, is outlined below:

$$\begin{array}{}\text{(3)}& p\left(\mathit{m}\mathrm{|}\mathit{d}\right)\propto p\left(\mathit{d}\mathrm{|}\mathit{m}\right)p\left(\mathit{m}\right),\end{array}$$

where p(m|d) is the posterior probability of the model (m) given d, p(m) is the prior information known about the model before the introduction of data, and p(d|m) is the likelihood. A Markov chain Monte Carlo (MCMC) methodology is used to sample the posterior distribution, traversing the space of admissible models with the statistics of the ensemble converging to the underlying posterior distribution, provided the chain of models is long enough.

We describe the 1-D variation of resistivity with depth as a piecewise constant function using Voronoi nuclei (see Killingbeck et al., 2018). Any available depth constraints separate the resistivity into different depth layers. In our case, in which constraints are drawn from high-resolution GPR data, we consider depth constraints to be exact since the accuracy of GPR depth estimation is at decimetre scale (Killingbeck et al., 2019) compared to the metre-scale resolution of the TEM. Within each layer, we define a single confined nucleus; aside from being confined to the given layer, this nucleus is otherwise unconstrained in depth. The number of confined layers, l, is equal to the number of layered depth constraints applied, with the addition of an assumed half space of constant resistivity extending to infinite depth. If no depth constraints are applied, then l=1, corresponding to a half space. We add in additional k nuclei in the model that are unconstrained in depth, termed “floating”. Our transdimensional framework allows the data to self-determine the required number of layers k (e.g. Bodin and Sambridge, 2009; Bodin et al., 2012; Livermore et al., 2018); thus, k is also an unknown that we determine.

The model vector, that describes the resistivity profile, is then

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\mathit{m}=[& {\mathrm{dp}}_{\mathrm{1}},{\mathrm{dp}}_{\mathrm{2}}\mathrm{\dots },{\mathrm{dp}}_k,\log \left(R_{\mathrm{1}}\right),\log \left(R_{\mathrm{2}}\right)\mathrm{\dots },\log \left(R_k\right),k,\\ & {\mathrm{dpc}}_{\mathrm{1}},{\mathrm{dpc}}_{\mathrm{2}}\mathrm{\dots },{\mathrm{dpc}}_l,\log \left({\mathrm{Rc}}_{\mathrm{1}}\right),\log \left({\mathrm{Rc}}_{\mathrm{2}}\right)\mathrm{\dots },\log \left({\mathrm{Rc}}_l\right)],\end{array}\end{array}$$

where dp_i is the floating nuclei depths, log (Rc_i) are the base-e log of the their respective resistivities, dpc_i is the confined nuclei depths, and log (Rc_i) is the base-e log of their respective resistivities. In our transdimensional framework, the number of floating nuclei (k) is a free parameter and self-determined in the algorithm.

For the choice of prior distribution in transdimensional calculations, it is worth noting that usually the geophysical properties of the cells (here the resistivity) and the cell depths are assumed independent, allowing a simple separated analytic form for the prior distribution (e.g. Bodin and Sambridge, 2009). This is followed in our simplest geometry with no GPR constraints, for which the prior distribution on the resistivity is depth independent and uniform with wide bounds on log (R) (e.g. R between 10⁰ and 10⁵ Ωm), to convey the fact that no prior information (beyond that which can be reasonably assumed for typical materials) is known about the subsurface. However, by interpreting any GPR-derived layers as different materials (Table 2) with much more narrowed ranges of resistivity, it is clear that a broad depth-independent prior distribution is no longer appropriate. Here, we allow the prior distribution of resistivity to depend on depth, by defining for each layer a different uniform distribution that reflects the tightened bounds from lithological information. This restricted prior distribution then significantly decreases the number of permissible models describing resistivity with depth, reducing model ambiguity from any given set of data. In terms of the model parameters, the prior of the resistivity for any nucleus is given by the specific layer that the nucleus is within (Killingbeck et al., 2018). Although a closed-form expression for the depth-dependent prior distribution cannot be easily formulated, in the algorithm, only the ratios of prior distributions are needed.

Lastly, the likelihood is defined by assuming that the measurements are normally distributed about values calculated from a forward model of TEM response (assuming a given resistivity profile) and the estimated standard deviation σ, where σ is dataset specific and defined for each different dataset. MuLTI-TEM uses the Leroi algorithm of the CSIRO and AMIRA project P223F (CSIRO and AMIRA, 2019) as a forward modeller to compare proposed subsurface models to the observed data. The Leroi algorithm is written in Fortran 95 and has a wide range of electromagnetic modelling capabilities; for more information, see Raiche (2008). We have used a simplified version of this algorithm and created a MEX file to call the code in MATLAB from within the MuLTI-TEM algorithm. See Table A1 and Fig. A1 in Appendix A for detailed TEM survey parameters input to MuLTI-TEM, defined from the Leroi forward modelling algorithm.

MuLTI-TEM numerically approximates the posterior distribution by creating an ensemble of models, traversing the model space and sampling the models with greater likelihood more often than models with a poor fit to the observed data. Provided the ensemble is sufficiently sampled, the numerically obtained posterior distribution will converge to the true posterior. This is achieved by constructing a Markov chain, each model in the chain being based on the previous model but randomly perturbed, the size of the perturbation being controlled by the user. We thin the Markov chain by using every 100th model when computing the distribution statistics, which suppresses any localised correlations of neighbouring models and speeds up convergence. MuLTI-TEM produces a variety of statistics of the resistivity ensemble including, but not limited to, the mean and mode (the most likely) solution and 95 % credible intervals as an estimate of its uncertainty, thus giving a profile with a quantified uncertainty.

3 Data acquisition

Data acquisition was performed on Midtdalsbreen, a NE-flowing outlet glacier of the Hardangerjøkulen ice cap in central-southern Norway (60.59^∘ N, 7.52^∘ E; Fig. 1a) in April–May 2018. Midtdalsbreen is surrounded by mountains of phyllite, crystalline granite and gneiss suggesting this as the underlying bedrock. The glacier is well suited for methodological development as it is logistically accessible, especially with multiple types of geophysical surveying equipment. More detailed information on previous glaciological and geophysical studies on Midtdalsbreen can be found in Andersen and Sollid (1971), Etzelmüller and Hagen (2005), Reinardy et al. (2013, 2019) and Willis et al. (2012).

Figure 1 (a) Location of Hardangerjøkulen ice cap in Norway and Google Earth image of the Midtdalsbreen survey area with nearest sources of TEM noise (town and railway line) highlighted. (b) Survey lines acquired during the 2018 field season with the seismic V_s results obtained at the top rock horizon displayed. The orange border around panel (b) identifies the same area as the orange box in panel (a). Note that panel (b) is rotated away from the north to enable optimal data comparison in later figures.

Figure 2 Survey configuration testing at the intersection of lines B and D. (a) 37 m × 37 m transmitter coil with receiver in the centre. (b) 10 m × 10 m transmitter coil with receiver 15 m offset. (c) 5 m × 5 m transmitter coil with receiver 12.5 m offset. (d) Raw data acquired at the intersection of B and D (237.5 Hz), from each survey configuration, plotted with background noise recorded with the transmitter turned off.

GPR, seismic and TEM surveys were performed around and over the glacier front (Fig. 1). All methods were acquired at each line highlighted, A–D, in the same field season. Lines B and C are located entirely on the glacier, whereas line A shows no glacier ice. Line D traverses through each of lines A, B and C, and extends beyond the glacier terminus. At the time of acquisition, the subsurface comprised snow (2–4 m thick) overlying a varying thickness (0–25 m) of glacier ice and a substrate of unknown subglacial material. This layered interpretation is based on the interpretation of the GPR dataset, which also suggests that the snow and ice layers show little variation in any of lines A, B and C.

Killingbeck et al. (2019) used MuLTI to jointly interpret the seismic and GPR data, defining regions of partially frozen sediment and hard bedrock based on subglacial shear-wave velocities (V_s<1000 m s⁻¹ for the former and >2000 m s⁻¹ for the latter. Figure 1b shows the seismic V_s results obtained at the glacier bed along all survey lines. This paper focuses initially on the joint inversion of the TEM and GPR data, and thereafter integrates the observations with the existing V_s distributions.

TEM data were acquired with a Geonics PROTEM47 system, consisting of a three-channel digital time-domain receiver, a TEM47 battery-powered transmitter and a 3-D multi-turn receiver coil. All survey parameter are listed in Table 1. For cross-glacier lines (A, B and C), the system was moved along the lines in 4 m intervals; for the longer down-glacier line (D), this was increased to 8 m. Multiple survey configurations were initially tested at the intersection of lines B and D to determine the optimal survey configuration for imaging the subglacial environment at Midtdalsbreen. These tests comprised the following (the maximum DOI of each test is estimated using Eq. 1, with voltage noise level 1 nV m⁻² (Fig. 2d), and the average underlying resistivity range was 0.1 to 1000 Ωm):

• a.

  37 m × 37 m square transmitter with receiver in centre, with estimated DOI 110 to 680 m (Fig. 2a);

• b.

  10 m × 10 m square transmitter with receiver 15 m away from centre of transmitter square, with estimated DOI 64 to 400 m (Fig. 2b); and

• c.

  5 m × 5 m square transmitter with receiver 12.5 m away from centre of transmitter square, with estimated DOI 45 to 300 m (Fig. 2c).

For accurate resistivity soundings, the raw signal should be above background noise levels (Fig. 2d). Background noise, measured with the transmitter coil turned off, is considered low at Midtdalsbreen since there are no large sources of electrical noise, e.g. power lines, buildings, roads, metal infrastructure, for ∼5 km (where the nearest town, Finse, is located; Fig. 1). The 10 m × 10 m square transmitter (Figs. 2b, 3) was chosen as the optimum survey configuration because

• i.

  it had a fast turn-off time (for imaging the shallow surface);

• ii.

  the raw signal (received voltage) recorded was sufficiently greater than the background noise (Fig. 2d);

• iii.

  it was easily deployed on the glacier and could be moved rapidly between the points of the survey lines (Fig. 3); and

• iv.

  the estimated DOI (ranging from 64 to 400 m depending on the underlying resistivity) is sufficient for imaging our target depth, subglacial sediments below ∼25 m thick ice.

Figure 3 (a) Survey configuration used for acquiring lines A, B, C and D on the glacier. (b) Image of the receiver unit on top of a rug to protect unit from snow and easily drag along the lines. (c) Image of transmitter unit sitting in bubble-wrap pocket used to protect unit and batteries from snow and cold.

Table 1 TEM survey parameters.

4 Application of MuLTI-TEM to a synthetic dataset

Synthetic TEM responses from a variety of models representing different possible glacial and subglacial structures of the Midtdalsbreen glacier were input into MuLTI-TEM for validation (models a–e in Fig. 4). Each model included layers of snow, ice and bedrock, with models (b)–(e) also including saturated subglacial sediment. Each layer was populated with representative resistivities from previous TEM studies (Mikucki et al., 2015). Certain models were designed to test particular aspects of the inversion: model (b) tested the maximum DOI using our specified survey design and synthetic, and model (c) tested whether the inversion can resolve a 5 m thin layer.

Figure 4 1-D synthetic block models created to simulate different subsurface scenarios expected at Midtdalsbreen.

Synthetic TEM responses were calculated from the 1-D block models using the Leroi forward modelling algorithm (Raiche, 2008), then normally distributed random noise, with a magnitude of 5 % of the signal at each time gate, was applied to all time gates, a similar noise model to Blatter et al. (2018); see Fig. 5. The simulated TEM survey configuration assumed a 10 m × 10 m square transmitter with a receiver 15 m away, consistent with the field acquisition.

Figure 5 Forward modelled responses for 1-D synthetic block models (a)–(e) with 5 % random noise applied. The lines within the circles represent the 5 % error bars.

The inversions were run using resistivity ranges shown in Table 2 using a maximum depth of 80 m, consistent with the estimated maximum DOI for this configuration. To highlight the benefit of additional depth constraints, MuLTI-TEM is run separately for an unconstrained case and a case in which the depths of snow-ice and ice-bed horizons are fixed. In total, 1 million iterations were sufficient for the posterior distribution to converge (a test of 2 million iterations produced the same posterior; some example test results are shown in Fig. A2). Additional figures in Fig. A3 show that the modal model from the depth-constrained inversion fits the data better than that from the non-constrained inversion, with, in general, a lower data misfit. Multiple chains were also tested with 1 million iterations using different initial conditions. For the constrained case, these produced similar posterior distributions with identical interpretation, indicating that only one chain was needed. For the unconstrained case, the posterior distributions differ slightly but are nevertheless qualitatively the same, suggesting that the unconstrained case is not yet converged (Fig. A4). More detailed inversion parameters used are documented in Table A2.

Table 2 Resistivity parameter boundaries used in MuLTI-TEM for the glacier feasibility study.

Posterior PDFs of the synthetic resistivity profiles produced from MuLTI-TEM are shown in Fig. 6. These are shown, within their 95 % credible interval, as coloured contours, where red indicates the most likely values. Consistent with many previous studies, both the unconstrained and constrained inversions indicate that the TEM method can resolve conductive structures much more accurately than resistive ones, highlighted by the much tighter PDF over the conductive sediment layer compared to the resistive ice layer. The unconstrained inversions (Fig. 6a) capture a similar structure to the true model, but they struggle to resolve true layer depths. The simple synthetic model (a) and thick resistive layered model (b) are relatively well resolved; however, the more complicated synthetic models with thin layers and large resistivity contrasts (c, d and e) are not. The depth-averaged resistivity errors within the subglacial layer (calculated from the difference of the modal and true solutions) of each synthetic non-constrained solution are (a) 275 Ωm, (b) 14 Ωm, (c) 30 Ωm, (d) 21 Ωm and (e) 2000 Ωm. The addition of depth constraints (Fig. 6b) improves the match throughout. The resistivity of the thin snow layers and conductive sediment layers is well-resolved in all synthetics, with depth-averaged resistivity errors within the subglacial layer reduced to (a) 8 Ωm, (b) 14 Ωm, (c) 23 Ωm, (d) 10 Ωm and (e) 25 Ωm (a factor of 34, 1, 1.3, 2.1 and 8 improvement on their unconstrained equivalents). Note that imaging beneath a conductive structure is difficult for the TEM method due to the attenuated signal (Blatter et al., 2018); however, our constrained inversion results show the bottom of the conductive layers (in b–e) to be much better resolved, i.e. to within <10 m.

Figure 6 Posterior distributions of resistivity determined from MuLTI-TEM inversion: (i) without depth constraints and (ii) with depth constraints. The models correspond to those shown in Fig. 4a–e, highlighted by the black line. The colour scale represents the probability density distribution of resistivity within the 95 % credible interval.

The TEM method is generally more sensitive to conductance (the product of conductivity and thickness) rather than the layer conductivity or thickness alone (Geonics, 1994). Therefore, modelling is challenged by a non-unique problem, for example, with thinner, more conductive layers producing a similar TEM signal to thicker, less conductive layers. The addition of depth constraints greatly reduces this non-uniqueness, enabling more accurate solutions to be obtained at all depths. However, an example where this TEM inversion struggles is when a thin conductive layer exists above a resistive basement (Fig. A6). In this example (Fig. A6a), a thin layer 1 m thick and resistivity 1 Ωm has conductance of 1 S, equivalent to a 10 m thick semi-conductive layer of 10 Ωm resistivity. Inverting this synthetic example with depth constraints in MuLTI-TEM shows that such a thin conductive layer above a resistive basement cannot be fully resolved even by the constrained inversion, appearing as a much thicker, less conductive layer.

This feasibility study highlights the significant added value of depth constraints when a complex resistivity structure is expected. It demonstrates that MuLTI-TEM is promising for the potential distributions of resistivity beneath Midtdalsbreen.

5 Application of MuLTI-TEM to the Midtdalsbreen dataset

5.1 1-D resistivity profiles

Using the data collected at Midtdalsbreen, we produced 1-D resistivity profiles for the soundings acquired at the midpoint of lines A, B and C (Fig. 7). Anomalous data points either induced by residual current still in the transmitter (at early time gates) or below the background noise level (at later time gates) were removed, shown as “X” symbols in Fig. 7b. Inversions were run with depth constraints taken from the GPR dataset for the snow-ice and glacier-bed interfaces (red and blue horizons, respectively, in Fig. 7). The same inversion parameters were used as the synthetic study (Table A2) but with the maximum depth extended to 160 m to test the limit of DOI. The data variance (σ) was kept at 5 % of the signal at each time gate as this was a good representation of the data variance of each data point calculated from the three stack recordings acquired in the field.

Figure 7 Results of the 1-D soundings acquired at the midpoint of lines A, B and C inverted using MuLTI-TEM. (a) Resistivity posterior probability distributions, with the maximum depth of investigation (DOI) plotted as the dotted black line. The solid blue and red lines highlight the snow–ice and ice–material depths. (b) Comparison of the observed data and 200 randomly chosen forward models from the model ensemble. The black “X” symbols show anomalous data points removed. (c) Posterior distribution of number of nuclei.

Posterior resistivity distributions are shown in Fig. 7a; comparisons of data fit with 200 randomly chosen forward models from the model ensemble are shown in Fig. 7b, and posterior distributions of the number of nuclei are shown in Fig. 7c. The estimated uncertainty for the most likely solution is calculated as one-half of the interquartile range at each depth, used due to non-normal PDFs. As is clear from Fig. 7a, conductive layers identified within the subglacial material are well resolved, with a tight PDF and low uncertainty. However, resistive layers (e.g. ice) have a wide PDF with large uncertainty estimates. For example, the uncertainty of the resistive ice layer is typically estimated as $\sim \pm {\mathrm{10}}^{\mathrm{4}}$ Ωm, whereas that of the conductive subglacial material is $\sim \pm {\mathrm{10}}^{\mathrm{2}}$ Ωm. The maximum DOI can also be identified in these distributions, expressed where the posterior distribution extends across the prior resistivity boundaries applied in a given layer (Table 2). This is 90 m for line A, 87 m for line B and 76 m for line C.

The 1-D inversions show a ∼10 m thick layer, of 10² Ωm resistivity, directly under the ice, underlain in turn by conductive material of 10⁰–10¹ Ωm for lines B and C. In contrast, line A (off the ice margin; see Fig. 1b) shows a ∼70 m thick resistive layer, 10⁴–10⁵ Ωm, immediately beneath the snow layer, underlain by either a very conductive thin layer (∼1 Ωm) or a conductive material that extends to greater depth (these cannot be distinguished as it is close to the limit of DOI). We ensured that our methodology can discriminate the thick conductive layers observed in lines B and C, and thin conductive layers underlain by resistive material approaching the DOI, using forward modelling. When compared (Fig. A5) to the data at the midpoint of line C, the thicker conductive model has the closest resemblance to the observed data. This analysis suggests a spatially varying pattern of subglacial resistivities from the front of the glacier up to line C.

5.2 2-D resistivity profiles

MuLTI-TEM is used to invert multiple independent 1-D soundings acquired along lines A, B, C (4 m intervals) and D (8 m intervals). Again, anomalous data points either induced by residual current in the transmitter (at early time gates) or below the background noise level (at later time gates) were removed, shown as the grey regions in Fig. 8 (left column). The raw signal acquired is generally above the background noise level for all time gates, except some anomalous points in the centre of line A, corresponding with anomalous points at the NE end of line D, highlighted in the left column of Fig. 8.

Figure 8 2-D inversion outputs for lines A–D from multiple 1-D MuLTI-TEM inversions. Left column: received voltages input to MuLTI-TEM; central column: most likely 2-D resistivity profiles; right column: estimated uncertainty (half the interquartile range of the posterior distribution). Snow and ice horizons are plotted in blue and red, respectively.

When using a central-loop TEM survey configuration, the 1-D response is simply located at the centre of the transmitter, where the receiver is located. However, with an offset transmitter–receiver survey configuration used in this study, the location of the 1-D sounding is a subject of debate. Some place the location below the receiver and others midway between the transmitter and receiver (Hoekstra and Blohm, 1990). The entire section between the transmitter and receiver is expected to influence the measurements, especially at late times, as the current is diffuse. In what follows, we assume the 1-D location of each sounding to be at the centre point between the transmitter and receiver, although we note subsurface conditions near the transmitter may have a slightly larger influence on the received voltage measure at early times, when the current loop radius is approximately the same as the transmitter loop radius and not overlapping the receiver offset position.

Inversions were run with depth constraints supplied from snow and ice horizons picked from the GPR data and using the same parameters as the synthetic study. We verified convergence of the solutions by running another Markov chain and increasing the chain length to 1.5 million iterations: all tests reproduced the same posterior distribution. Consistent with initial observations in the 1-D analysis, the 2-D resistivity profiles (Fig. 8, central column) highlight a wide range of subglacial resistivity values, from 10⁰ to 10⁵ Ωm, both in along- and cross-glacier profiles; however, there are some key consistent observations between these profiles, including

• a ∼10² Ωm layer directly below the ice, for ice thicknesses <20 m (mainly observed in line B), which varies in thickness,

• a high-resistivity layer, 10⁴–10⁵ Ωm, in line A which matches line D at their intersection, and

• a lowermost layer of highly conductive material, ∼ 10⁰–10¹ Ωm, that generally extends to the DOI.

The estimated uncertainties for the mode resistivity solutions are displayed in Fig. 8 (rightmost column). This reiterates previous observations made from the synthetic study and 1-D analysis that conductive layers identified within the subglacial material are well resolved by the TEM method with low uncertainty; however, TEM methods struggle to fully resolve the more resistive layers which therefore have a larger uncertainty. We note that marginal uncertainties along a 2-D line can also be displayed as a 3-D probability cube, with axes representing resistivity, line distance and depth, and the colour bar representing the probability (e.g. Ray et al., 2014). This aids visualisation of the Bayesian solution and its uncertainty, particularly useful when characterising an anomalous target from a constant background resistivity, e.g. a subglacial aquifer or lake underlain by bedrock.

6 Interpretation and discussion

6.1 Joint interpretation of MuLTI-TEM with MuLTI seismic results

The variability of resistivity (10⁰–10⁵ Ωm) in our TEM profiles suggests a complex subsurface structure, in which subglacial water may be in liquid and frozen states. However, the nature of the matrix – whether sediment or bedrock – cannot be determined from resistivity alone. Resistivity is related to the resistivity of the pore fluids divided by the fractional porosity. A commonly used approximation is given by Archie's law, which states that

$$\begin{array}{}\text{(5)}& R=a{R}_W/{\mathit{\varnothing }}^m,\end{array}$$

where R is the bulk resistivity of a saturated porous medium, ∅ is the porosity, R_W is the pore fluid resistivity, and m and a are empirical quantities determined by the geometry of the pores (Archie, 1952). However, it is difficult to distinguish material type, such as sediment or bedrock, from resistivity alone. By contrast, seismic shear-wave methods are sensitive to the shear modulus, or stiffness, of a material but are insensitive to water content. However, a combined interpretation of resistivity and V_s profiles can be used to define a mutually consistent system to characterise the material properties and water content of the subglacial environment. We do this using complementary seismic data acquired alongside our TEM acquisitions in 2018.

Table 3 V_s and resistivity ranges for subglacial material lithologies, used in analysis of both seismic and TEM. Material types have been defined from King et al. (1988), Mikucki et al. (2015), and Killingbeck et al. (2019).

Figure 9 Joint interpretation of V_s and resistivity profiles for lines A–D. Left column: modal V_s solution. Central column: modal resistivity solution. An example of areas with high V_s and high R is shown in line A, low V_s and low R is shown in line B, and high V_s with low R is shown in line D. Right column: estimated subglacial material when applying V_s and resistivity conditions of Table 3. Note the sparse, disconnected areas identified as bedrock in line A are regarded as a misallocation.

The initial interpretation of the seismic data was presented in Killingbeck et al. (2019). This study excluded certain phase velocities on the grounds that they were too high (Fig. A7), but this merits re-evaluation when compared with the co-located observation of high resistivity. Therefore, in this integrated interpretation, the high-phase velocities are included, thus providing broader bandwidth dispersion curves.

Observing trends of V_s and resistivity in our profiles, three clear patterns emerge within the subglacial material:

• i.

  zones of low V_s and low resistivity,

• ii.

  zones of high V_s and high resistivity, and

• iii.

  zones of high V_s with low resistivity (Fig. 9; leftmost and centre columns).

These patterns have been used to define three different material types (Table 3). From previous electromagnetic and seismic studies (e.g. King et al., 1988; Schneider et al., 2013; Wu et al., 2017), liquid water in the pores of unconsolidated material has a low resistivity and low V_s (we define this as partially frozen sediment), whereas frozen water in pores is very resistive, with a high V_s (defined as frozen sediment/permafrost). We assume the bedrock is comprised of phyllite, crystalline granite and gneiss with a high V_s and high resistivity. However, the resistivity profiles show the bedrock to have a very low resistivity, suggesting it could be highly fractured and weathered with saline water in the fractures. The presence of saline water can decrease the electrical resistivity by as much as 9 orders of magnitude (Olhoeft, 1981).

The resistivity and V_s profiles are linearly interpolated such that they have mutually consistent sample intervals (1 m) and depth extents (40 m). The TEM and seismic profiles were acquired in the same field season along the same lines; the acquisition parameters were chosen so the two methods could be directly comparable. However, the location of each 1-D TEM sounding and seismic common midpoint gather is offset by 2 m; therefore, we linearly interpolated and resampled both the resistivity and V_s solutions, originally sampled every 4 m (lines A, B and C) and 8 m (line D), to every metre, thus making them directly comparable. Killingbeck et al. (2019) consider abrupt lateral variations in the V_s profiles to be noise and interpret the broader variation in lateral and vertical character as the representative velocity structure. Therefore, since the lateral resolution of the V_s profiles (estimated from the length of the geophone spread; 30–40 m) is larger than that of the resistivity profiles (estimated from the horizontal spacing between 1-D soundings; 4–8 m), we apply lateral smoothing to the V_s profiles during the joint analysis. After these steps, we obtain a smooth joint interpretation of the predicted subglacial material shown in Fig. 9 (rightmost column) and Fig. 10.

Figure 10 (a) 3-D cross section of lines A–D, showing the subglacial material estimated from applying V_s and resistivity conditions stated in Table 3. (b) Depth slice through 3-D cross at the top rock horizon. (c) Depth slice through 3-D cross section at 38 m.

The joint interpretation shows zones of mainly sediment and thick permafrost (>40 m) at the front of the glacier, along line A, matching observations at the corresponding intersection point with line D. Note the sparse, disconnected areas identified as bedrock in line A are regarded as a misallocation, potentially due to abrupt lateral variations in V_s, which we regard as noise. Unfrozen sediment occurs directly below the ice at line B, with a varying thickness of 30 m at the eastern end to 5–10 m towards the western end. In contrast, a mixture of frozen and unfrozen sediment is observed directly below the ice at line C, typically ∼10 m thick. Underlying the ice and sediment layers in all lines is the conductive bedrock, with its structure shown clearly in the cross-glacier line (D). This also highlights a thin zone of frozen sediment directly under the glacier tongue, matching observations from the GPR data shown in Reinardy et al. (2019), suggesting there is a frozen tongue.

6.2 Discussion

MuLTI-TEM combines a probabilistic approach with external depth constraints to mitigate ambiguous, non-unique solutions found in conventional TEM inversions. It provides a robust quantitative uncertainty analysis of any chosen model at all depth levels, also providing an accurate estimate of DOI using the posterior distribution. The addition of depth constraints improves the characterisation of material, particularly beneath conductive layers and enables a faster convergence of the solution, as demonstrated in the synthetic study. We note other methods could be used to enhance the efficiency of the transdimensional inversion, potentially providing better convergence rates, such as proposing the birth parameters from the prior (instead of a Gaussian distribution), e.g. Dosso et al. (2014). Further still, having access to the full posterior distribution enables subsets of the posterior model probabilities to be selected, testing various hypothesis about the model structure (Ray and Key, 2012).

Nevertheless, the success of MuLTI-TEM depends fundamentally on the input data quality and its suitability for the specific target imaged. With TEM methods, it is often not possible to separately determine the conductivity and thickness; only the conductance (product of thickness and conductivity) can be determined. Therefore, thorough synthetic modelling should be undertaken before acquiring data in the field to determine if the survey design and time range of measurements are sufficient and suitable to detect specific targets.

Although this paper focuses on the specific TEM survey design used in this study, a ground-based 10×10 m transmitter with receiver 15 m away, MuLTI-TEM can be used with most TEM datasets. The Leroi forward modelling code can be used in frequency or time-domain mode to model most TEM transmitter–receiver combinations, ground-based or airborne, and needs only to be adapted to the user-specific TEM configuration (Raiche, 2008). Equally, depth constraints are provided as numerical inputs and can therefore be supplied from any external source, e.g. borehole measurements or any complementary inference from an external geophysical data source. In our Midtdalsbreen case study, the uncertainty in the depth constraints applied is negligible (decimetre accuracy from GPR data) compared to the observed data uncertainty (metre accuracy from TEM), motivating us to fix the internal interface depths. However, there remains a finite resolution in GPR data; hence, we are considering a modification to the MuLTI-TEM code to make it compatible with uncertain interface depths. This would also benefit depth constraints supplied from more uncertain data sources, thus making MuLTI-TEM more broadly applicable.

This paper has presented how a joint analysis of three geophysical datasets can increase our understanding of the material in the subsurface and provide a more detailed interpretation. Using GPR information as a depth constraint, we have combined insight from TEM and seismic shear-wave methods to provide a detailed characterisation of the material beneath the margins of Midtdalsbreen. Critically, TEM data reveal hydrological properties to which the seismic analysis was insensitive, whereas the seismic data indicate the varying stiffness of the subglacial material. Future extensions of this interpretative strategy could include petrophysical relationships to obtain and/or guide interpretations of the volumetric proportions of water, ice and air in the subsurface (e.g. Hauck et al., 2008). A further promising extension would be a modification to calculate the joint distribution of resistivity and V_s (rather than only the marginal distributions discussed in this paper), which could lead to a more accurate understanding of the subsurface structure (utilising the structural similarities between resistivity and seismic velocity; e.g. Wisén and Christiansen, 2005). Such a combined approach would also result in more detailed analysis of the Midtdalsbreen margin, including a probabilistic facies classification, leading to a framework by which aquifer properties, such as porosity, water content and pore fluid conductivity/salinity, beneath large ice masses could be quantified. This could have a direct impact on basal parameters used as input to ice-flow models for a better prediction of ice motion over time and hence future sea level rise.

7 Conclusions

The material properties of the subglacial environment, in particular their water content and saturation, can be characterised by inferences of their resistivity, obtained from TEM measurements. However, conventional TEM inversions provide solutions that are non-unique with no quantification of uncertainty estimates in depth and resistivity. This paper has presented an inversion algorithm, MuLTI-TEM, used to overcome such problems. Our method uses a transdimensional Bayesian inversion approach adapted from the MuLTI algorithm (Killingbeck et al., 2018), which incorporates independent depth constraints to limit the solution space reducing ambiguity. Synthetic testing of multiple different scenarios representing a small glacier underlain by sediment showed the addition of depth constraints greatly improves numerical convergence. This results in constrained solutions having a large improvement in the depth-average uncertainty of the output model, an average factor of 15 improvement on their unconstrained equivalents, with little computational power needed to obtain these results.

A joint interpretation, using V_s and resistivity boundaries, of the MuLTI-TEM results with MuLTI seismic surface wave results, presented in Killingbeck et al. (2019), considers three subglacial material classifications: sediment (V_s<1600 m s⁻¹, 50 Ωm < R<500 Ωm), permafrost (V_s>1600 m s⁻¹, R>500 Ωm) and weathered/fractured bedrock with saline water in the fractures (V_s>1900 m s⁻¹, R<50 Ωm). Their spatial extent, within Midtdalsbreen's subglacial environment, shows a mixture of sediment and permafrost directly below the ice and in the moraine at the front of the glacier, underlain by bedrock.

MuLTI-TEM is highly versatile, being compatible with most TEM survey designs, ground-based or airborne, as the Leroi forward modelling code can model most transmitter–receiver combinations, along with the depth constraints being provided from any external source. This study presents novel methodologies, through MuLTI-TEM and MuLTI, by which other glacier and ice-sheet subglacial material can be explored, highlighting the importance of acquiring multiple geophysical datasets for accurately characterising the subglacial environment.

Appendix A

Table A1 TEM survey parameters input to MuLTI-TEM, defined from the Leroi forward modelling algorithm.

Figure A1 Schematic image of the transmitter waveform, highlighting the waveform parameters defined in MuLTI-TEM described in Table A1.

Table A2 Inversion parameters used in MuLTI-TEM for the synthetic feasibility study and 1-D and 2-D real data inversions, explained further in Killingbeck et al. (2018). Burn-in number is the number of iterations discounted at the start of the chain to remove any dependencies of the initial conditions. Sigma resistivity, change, move and birth are user-specified parameters that determine the magnitude of the four different perturbations that can be applied (change resistivity, move nucleus, give birth to a new nucleus and remove a nucleus).

Figure A2 Posterior distributions of resistivity determined from MuLTI-TEM inversion for the synthetic models (b), (c) and (e) without depth constraints applied and with (a) 1 million iterations and (b) 2 million iterations.

Figure A3 Synthetic inversion results for all models in Fig. 4a–e, showing comparison of data fit and posterior distribution of number of nuclei plots.

Figure A4 Posterior distribution of resistivity for synthetic model (d) using three different chains for the non-constrained and constrained inversions with MuLTI-TEM, highlighting the independence of the constrained solutions on chain index, whereas the non-constrained distributions show a weak variation across chains.

Figure A5 (a) Example synthetic models created with a conductive thin layer on top of a resistive basement (green) and a conductive basement (blue), compared to the original synthetic model (d). (b) Calculated responses, using the Leroi forward modelling code, of the synthetics compared to the observed data at the midpoint of line C.

Figure A6 MuLTI-TEM inverted result for the (a) conductive thin layer synthetic model and the (b) conductive basement synthetic model, shown in Fig. A5.

Figure A7 (a) Example of re-picked dispersion curves along line A at 20, 40 and 70 m. The noisy high-phase velocities observed >55 Hz (originally not picked) are now thought to be real, matching the high-resistivity observations along line A. Red circles highlight the extra high frequencies picked. (b) Most likely resistivity profile with locations of dispersion curves marked by the black lines.