the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Characterisation of subglacial water using a constrained transdimensional Bayesian transient electromagnetic inversion
Adam D. Booth
Philip W. Livermore
C. Richard Bates
Landis J. West
Subglacial water modulates glacierbed 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 (timedomain) electromagnetic (TEM) soundings, but these are often nonunique. 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 groundbased TEM data acquired on the terminus of the Norwegian glacier, Midtdalsbreen, with depth constraints provided by colocated groundpenetrating radar data. Our inversion shows that the glacier bed is directly underlain by material of resistivity 10^{2} Ωm ± 1000 %, with thickness 5–40 m, in turn underlain by a highly conductive basement (10^{0} Ωm ± 15 %). Highresistivity material, 5×10^{4} Ω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 colocated seismic shearwave velocity inversions to further reduce ambiguity in the hydrogeological interpretation of the subsurface, we propose a new 3D interpretation in which the Midtdalsbreen subglacial material is partitioned into partially frozen sediment, frozen sediment/permafrost and weathered/fractured bedrock with saline water.
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 iceflow 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.
Noninvasive 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 groundpenetrating 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 nonunique, 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 groundtruth 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 groundpenetrating 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.
Timedomain electromagnetic methods use electromagnetic fields to sound the subsurface structure. Here, we use the transient timedomain 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 debriscovered 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 nonunique: 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 nonlinear least squares (Barnett, 1984) and an Occamtype 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 halfspace skin depth (Spies, 1989) or the Jacobian sensitivity matrix (Christiansen and Auken, 2012), though these do not consider the nonlinear sensitivity of the DOI to conductivity structure. These limitations in uncertainty quantification, fixed model space and DOI estimation can be mitigated by transdimensional Bayesian samplingbased 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 “MuLTITEM” (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 shearwave velocity (V_{s}) profiles as local accumulations of soft material (partially frozen glacial till) and permafrost overlying bedrock. Finally, we show a combined 3D interpretation of previous shearwave inversions (Killingbeck et al., 2019) and results output from MuLTITEM, and suggest the development of a joint resistivity–V_{s} depthconstrained inversion strategy.
2.1 Transient electromagnetics
In TEM surveying, an electromagnetic field is generated by sending a periodic, modified squarewave 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
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 MuLTITEM
MuLTITEM 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 MuLTITEM 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 σ^{2}, 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
The method used to find the posterior distribution of the resistivity profile, is outlined below:
where p(md) 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(dm) 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 1D 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 highresolution 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 metrescale 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 selfdetermine 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
where dp_{i} is the floating nuclei depths, log (Rc_{i}) are the basee log of the their respective resistivities, dpc_{i} is the confined nuclei depths, and log (Rc_{i}) is the basee log of their respective resistivities. In our transdimensional framework, the number of floating nuclei (k) is a free parameter and selfdetermined 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^{0} and 10^{5} Ω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 GPRderived layers as different materials (Table 2) with much more narrowed ranges of resistivity, it is clear that a broad depthindependent 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 closedform expression for the depthdependent 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. MuLTITEM 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 MuLTITEM algorithm. See Table A1 and Fig. A1 in Appendix A for detailed TEM survey parameters input to MuLTITEM, defined from the Leroi forward modelling algorithm.
MuLTITEM 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. MuLTITEM 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.
Data acquisition was performed on Midtdalsbreen, a NEflowing outlet glacier of the Hardangerjøkulen ice cap in centralsouthern 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).
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 shearwave velocities (V_{s}<1000 m s^{−1} for the former and >2000 m s^{−1} 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 threechannel digital timedomain receiver, a TEM47 batterypowered transmitter and a 3D multiturn receiver coil. All survey parameter are listed in Table 1. For crossglacier lines (A, B and C), the system was moved along the lines in 4 m intervals; for the longer downglacier 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^{−2} (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 turnoff 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.
Synthetic TEM responses from a variety of models representing different possible glacial and subglacial structures of the Midtdalsbreen glacier were input into MuLTITEM 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.
Synthetic TEM responses were calculated from the 1D 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.
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, MuLTITEM is run separately for an unconstrained case and a case in which the depths of snowice and icebed 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 depthconstrained inversion fits the data better than that from the nonconstrained 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.
Posterior PDFs of the synthetic resistivity profiles produced from MuLTITEM 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 depthaveraged resistivity errors within the subglacial layer (calculated from the difference of the modal and true solutions) of each synthetic nonconstrained 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 wellresolved in all synthetics, with depthaveraged 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.
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 nonunique 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 nonuniqueness, 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 semiconductive layer of 10 Ωm resistivity. Inverting this synthetic example with depth constraints in MuLTITEM 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 MuLTITEM is promising for the potential distributions of resistivity beneath Midtdalsbreen.
5.1 1D resistivity profiles
Using the data collected at Midtdalsbreen, we produced 1D 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 snowice and glacierbed 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.
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 onehalf of the interquartile range at each depth, used due to nonnormal 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 1D inversions show a ∼10 m thick layer, of 10^{2} Ωm resistivity, directly under the ice, underlain in turn by conductive material of 10^{0}–10^{1} Ω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^{4}–10^{5} Ω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 2D resistivity profiles
MuLTITEM is used to invert multiple independent 1D 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.
When using a centralloop TEM survey configuration, the 1D 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 1D 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 1D 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 1D analysis, the 2D resistivity profiles (Fig. 8, central column) highlight a wide range of subglacial resistivity values, from 10^{0} to 10^{5} Ωm, both in along and crossglacier profiles; however, there are some key consistent observations between these profiles, including

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

a highresistivity layer, 10^{4}–10^{5} Ωm, in line A which matches line D at their intersection, and

a lowermost layer of highly conductive material, ∼ 10^{0}–10^{1} Ω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 1D 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 2D line can also be displayed as a 3D 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.1 Joint interpretation of MuLTITEM with MuLTI seismic results
The variability of resistivity (10^{0}–10^{5} Ω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
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 shearwave 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.
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 reevaluation when compared with the colocated observation of high resistivity. Therefore, in this integrated interpretation, the highphase 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 1D 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 1D 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.
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 crossglacier 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
MuLTITEM combines a probabilistic approach with external depth constraints to mitigate ambiguous, nonunique 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 MuLTITEM 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 groundbased 10×10 m transmitter with receiver 15 m away, MuLTITEM can be used with most TEM datasets. The Leroi forward modelling code can be used in frequency or timedomain mode to model most TEM transmitter–receiver combinations, groundbased or airborne, and needs only to be adapted to the userspecific 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 MuLTITEM code to make it compatible with uncertain interface depths. This would also benefit depth constraints supplied from more uncertain data sources, thus making MuLTITEM 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 shearwave 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 iceflow models for a better prediction of ice motion over time and hence future sea level rise.
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 nonunique with no quantification of uncertainty estimates in depth and resistivity. This paper has presented an inversion algorithm, MuLTITEM, 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 depthaverage 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 MuLTITEM results with MuLTI seismic surface wave results, presented in Killingbeck et al. (2019), considers three subglacial material classifications: sediment (V_{s}<1600 m s^{−1}, 50 Ωm < R<500 Ωm), permafrost (V_{s}>1600 m s^{−1}, R>500 Ωm) and weathered/fractured bedrock with saline water in the fractures (V_{s}>1900 m s^{−1}, 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.
MuLTITEM is highly versatile, being compatible with most TEM survey designs, groundbased 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 MuLTITEM and MuLTI, by which other glacier and icesheet subglacial material can be explored, highlighting the importance of acquiring multiple geophysical datasets for accurately characterising the subglacial environment.
MuLTITEM can be found at https://doi.org/10.5281/zenodo.3471638 (Killingbeck and Livermore, 2019).
SFK, ADB, PWL and LJW designed the project. SFK and ADB acquired the data. SFK and PWL developed MuLTI TEM. CRB provided advice and support while using the TEM equipment. SFK prepared the manuscript with contributions from all coauthors.
The authors declare that they have no conflict of interest.
The timedomain electromagnetic equipment was supplied by NERC Geophysical Equipment Facility, loan 1090. Fieldwork was greatly assisted by Emma Pearce, James Killingbeck and Kjell Magne Tangen. Alan Hobbs and all NERC GEF staff are thanked for their support and advice throughout the project.
This research has been supported by the UK NERC SPHERES DTP (grant no. NE/L002574/1). Fieldwork was funded by the research project “Snow Accumulation Patterns on Hardangerjøkulen Ice Cap (SNAP)”, itself funded by the European Union's Horizon 2020 project INTERACT, under grant agreement no. 730938.
This paper was edited by Ulrike Werban and reviewed by Anandaroop Ray and one anonymous referee.
Andersen, J. L. and Sollid, J. L.: Glacial chronology and glacial geomorphology in the marginal zones of the glaciers Midtdalsbreen and Nigardsbreen, south Norway, Norsk Geogr. Tidsskr., 25, 1–38, 1971.
Archie, G. E.: Classification of carbonate reservoir rocks and petrophysical considerations, AAPG Bull., 36, 278–298, 1952.
Auken, E., Jørgensen, F., and Sørensen, K. I.: Largescale TEM investigation for groundwater, Explor. Geophys., 34, 188–194, 2003.
Auken, E., Christiansen, A. V., Kirkegaard, C., Fiandaca, G., Schamper, C., Behroozmand, A. A., Binley, A., Nielsen, E., Effersø, F., Christensen, N. B., and Sørensen, K.: An overview of a highly versatile forward and stable inverse algorithm for airborne, groundbased and borehole electromagnetic and electric data, Explor. Geophys., 46, 223–235, 2015.
Barnett, C. T.: Simple inversion of timedomain electromagnetic data, Geophysics, 49, 925–933, 1984.
Blatter, D., Key, K., Ray, A., Foley, N., Tulaczyk, S., and Auken, E.: Transdimensional Bayesian inversion of airborne transient EM data from Taylor Glacier, Antarctica, Geophys. J. Int., 214, 1919–1936, 2018.
Bodin, T. and Sambridge, M.: Seismic tomography with the reversible jump algorithm, Geophys. J. Int., 178, 1411–1436, 2009.
Bodin, T., Sambridge, M., Tkalčić, H., Arroucau, P., Gallagher, K., and Rawlinson, N.: Transdimensional inversion of receiver functions and surface wave dispersion, J. Geophys. Res.Sol. Ea., 117, B02301, https://doi.org/10.1029/2011JB008560, 2012.
Booth, A. D., Clark, R. A., Kulessa, B., Murray, T., Carter, J., Doyle, S., and Hubbard, A.: Thinlayer effects in glaciological seismic amplitudeversusangle (AVA) analysis: implications for characterising a subglacial till unit, Russell Glacier, West Greenland, The Cryosphere, 6, 909–922, https://doi.org/10.5194/tc69092012, 2012.
Christensen, N. B. and Tølbøll, R. J.: A lateral model parameter correlation procedure for onedimensional inverse modelling, Geophys. Prospect., 57, 919–929, 2009.
Christiansen, A. V. and Auken, E.: A global measure for depth of investigation, Geophysics, 77, WB171–WB177, 2012.
Christoffersen, P., Bougamont, M., Carter, S. P., Fricker, H. A., and Tulaczyk, S.: Significant groundwater contribution to Antarctic ice streams hydrologic budget, Geophys. Res. Lett., 41, 2003–2010, 2014.
Constable, S. C., Parker, R. L., and Constable, C. G.: Occam's inversion: A practical algorithm for generating smooth models from electromagnetic sounding data, Geophysics, 52, 289–300, 1987.
CSIRO and AMIRA: P223Suite Open source exploration geophysical electromagnetic modelling tools, available at: http://p223suite.sourceforge.net/, last access: 1 February 2019.
Dosso, S. E., Dettmer, J., Steininger, G., and Holland, C. W.: Efficient transdimensional Bayesian inversion for geoacoustic profile estimation, Inverse Probl., 30, 114018, https://doi.org/10.1088/02665611/30/11/114018, 2014.
Etzelmüller, B. and Hagen, J. O.: Glacier–permafrost interaction in Arctic and alpine mountain environments with examples from southern Norway and Svalbard, Geol. Soc. Spec. Publ., 242, 11–27, 2005.
Geonics: PROTEM47D operating manual, Geonics Ltd., Mississauga, Ontario, Canada, 1994.
Hauck, C., Guglielmin, M., Isaksen, K., and Vonder Mühll, D.: Applicability of frequencydomain and timedomain electromagnetic methods for mountain permafrost studies, Permafrost Periglac., 12, 39–52, 2001.
Hauck, C., Bach, M., and Hilbich, C.: A 4phase model to quantify subsurface ice and water content in permafrost regions based on geophysical data sets, in: Proceedings Ninth International Conference on Permafrost, Fairbanks, Vol. 1, edited by: Kane, D. L. and Hinkel, K. M., Institute of Northern Engineering, University of Alaska Fairbanks, 675–680, 2008.
Hoekstra, P. and Blohm, M. W.: Case histories of timedomain electromagnetic soundings in environmental geophysics, Geotechnical and Environmental Geophysics, 2, 1–15, 1990.
Hoekstra, P. and McNeill, D.: Electromagnetic probing of permafrost, in: Proceedings of Second International Conference on Permafrost, 517–526, Yakutsk, Russia, 1973.
Killingbeck, S. and Livermore, P.: MuLTITEM, Zenodo, https://doi.org/10.5281/zenodo.3471638, 2019.
Killingbeck, S. F., Livermore, P. W., Booth, A. D., and West, L. J.: Multimodal Layered Transdimensional Inversion of Seismic Dispersion Curves with Depth Constraints, Geochem. Geophy. Geosy., 19, 4957–4971, 2018.
Killingbeck, S. F., Booth, A. D., Livermore, P. W., West, L. J., Reinardy, B. T. I., and Nesje, A.: Subglacial sediment distribution from constrained seismic inversion, using MuLTI software: Examples from Midtdalsbreen, Norway, Ann. Glaciol., 60, 206–219, 2019.
King, M. S., Zimmerman, R. W., and Corwin, R. F.: Seismic and electrical properties of unconsolidated permafrost, Geophys. Prospect., 36, 349–364, 1988.
Livermore, P. W., Fournier, A., Gallet, Y., and Bodin, T.: Transdimensional inference of archeomagnetic intensity change, Geophys. J. Int., 215, 2008–2034, 2018.
Merz, K., Maurer, H., Rabenstein, L., Buchli, T., Springman, S. M., and Zweifel, M.: Multidisciplinary geophysical investigations over an alpine rock glacier, Geophysics, 81, WA1–WA11, 2016.
Mikucki, J. A., Auken, E., Tulaczyk, S., Virginia, R. A., Schamper, C., Sørenson, K. I., Doran, P. T., Dugan, H., and Foley, N.: Deep groundwater and potential subsurface habitats beneath an Antarctic dry valley, Nat. Commun., 6, 6381, https://doi.org/10.1038/ncomms7831, 2015.
Minsley, B. J., Abraham, J. D., Smith, B. D., Cannia, J. C., Voss, C. I., Jorgenson, M. T., Walvoord, M. A., Wylie, B. K., Anderson, L., Ball, L. B., DeszczPan, M., Wellman, B. K., and Ager, T. A.: Airborne electromagnetic imaging of discontinuous permafrost, Geophys. Res. Lett., 39, L02503, https://doi.org/10.1029/2011GL050079, 2012.
Mosegaard, K. and Tarantola, A.: Monte Carlo sampling of solutions to inverse problems, J. Geophys. Res.Sol. Ea., 100, 12431–12447, 1995.
Olhoeft, G. R.: Electrical properties of granite with implications for the lower crust, J. Geophys. Res.Sol. Ea., 86, 931–936, 1981.
Raiche, A.: Leroi version 8.0, Fortran 95 program, developed for AMIRA project P223F, available at: http://p223suite.sourceforge.net/descriptions.html (last access: 12 February 2019), 2008.
Ray, A. and Key, K.: Bayesian inversion of marine CSEM data with a transdimensional self parametrizing algorithm, Geophys. J. Int., 191, 1135–1151, 2012.
Ray, A., Key, K., Bodin, T., Myer, D., and Constable, S.: Bayesian inversion of marine CSEM data from the Scarborough gas field using a transdimensional 2D parametrization, Geophys. J. Int., 199, 1847–1860, 2014.
Reinardy, B. T. I., Leighton, I., and Marx, P. J.: Glacier thermal regime linked to processes of annual moraine formation at Midtdalsbreen, southern Norway, Boreas, 42, 896–911, 2013.
Reinardy, B. T. I., Booth, A. D., Hughes, A. L. C., Boston, C. M., Åkesson, H., Bakke, J., Nesje, A., Giesen, R. H., and Pearce, D. M.: Pervasive cold ice within a temperate glacier – implications for glacier thermal regimes, sediment transport and foreland geomorphology, The Cryosphere, 13, 827–843, https://doi.org/10.5194/tc138272019, 2019.
Rozenberg, G., Henderson, J. D., Sartorelli, A. N., and Judge, A. S.: Some aspects of transient electromagnetic soundings for permafrost delineation, in: Workshop on Permafrost Geophysics, Special Report 855, 74–90, Cold Reg. Res. Lab, Hanover, NH, 1985.
Schneider, S., Daengeli, S., Hauck, C., and Hoelzle, M.: A spatial and temporal analysis of different periglacial materials by using geoelectrical, seismic and borehole temperature data at Murtèl–Corvatsch, Upper Engadin, Swiss Alps, Geogr. Helv., 68, 265–280, https://doi.org/10.5194/gh682652013, 2013.
Siegert, M. J., Kulessa, B., Bougamont, M., Christoffersen, P., Key, K., Andersen, K. R., Booth, A. D., and Smith, A. M.: Antarctic subglacial groundwater: a concept paper on its measurement and potential influence on ice flow, Geol. Soc. Spec. Publ., 197–213, 2018.
Spies, B.: Depth of investigation in electromagnetic sounding methods, Geophysics, 54, 872–888, 1989.
Todd, B. J. and Dallimore, S. R.: Electromagnetic and geological transect across permafrost terrain, Mackenzie River delta, Canada, Geophysics, 63, 1914–1924, 1998.
Vignoli, G., Fiandaca, G., Christiansen, A. V., Kirkegaard, C., and Auken, E.: Sharp spatially constrained inversion with applications to transient electromagnetic data, Geophys. Prospect., 63, 243–255, 2015.
Willis, I. C., Fitzsimmons, C. D., Melvold, K., Andreassen, L. M., and Giesen, R. H.: Structure, morphology and water flux of a subglacial drainage system, Midtdalsbreen, Norway, Hydrol. Process., 26, 3810–3829, 2012.
Wisén, R. and Christiansen, A. V.: Laterally and mutually constrained inversion of surface wave seismic data and resistivity data, J. Environ. Eng. Geoph., 10, 251–262, 2005.
Wu, Y., Nakagawa, S., Kneafsey, T. J., Dafflon, B., and Hubbard, S.: Electrical and seismic response of saline permafrost soil during freezeThaw transition, J. Appl. Geophys., 146, 16–26, 2017.
 Abstract
 Introduction
 Method
 Data acquisition
 Application of MuLTITEM to a synthetic dataset
 Application of MuLTITEM to the Midtdalsbreen dataset
 Interpretation and discussion
 Conclusions
 Appendix A
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Method
 Data acquisition
 Application of MuLTITEM to a synthetic dataset
 Application of MuLTITEM to the Midtdalsbreen dataset
 Interpretation and discussion
 Conclusions
 Appendix A
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References