Characterisation of subglacial water using a constrained transdimensional Bayesian Time Domain Electromagnetic Inversion

Subglacial water influences 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 10 from electrical resistivity depth profiles. Such profiles can be obtained from inversions of time domain electromagnetics (TEM) soundings, but these are often non-unique. Here, we adapt an existing Bayesian transdimensional algorithm (‘MuLTI’) to the inversion of TEM data constrained by 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 15 shows that the glacier bed is directly underlain by material of resistivity 10 Ωm ±100 %, with thickness 5-40 m, in turn underlain by a highly conductive basement (10 Ωm ±15 %). High resistivity material, 5x10 Ω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 co-located seismic shear-wave velocity inversions to further reduce ambiguity in the hydro-geological interpretation of the subsurface, we propose a new 3D interpretation of the Midtdalsbreen subglacial material partitioned into 20 partially frozen sediment, frozen sediment/permafrost and weathered/fractured bedrock with saline water.


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 25 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 is 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 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 andHauck et al., 2001). Further still, inversions of isolated geophysical datasets are unconstrained and 5 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 time-domain 10 electromagnetics (TEM), with and depth constraints derived from ground penetrating radar (GPR), to are combined in a constrained inversion 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' 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. 15 TEM methods use electromagnetic fields to investigate subsurface resistivity structure by measuring eddy currents induced by current transmitted through either a grounded-wire, coincided loop or offset transmitter coil. The method has a depth sensitivity ranging from a few meters 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 20 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 . 25 In TEM, the resistivity-depth profile is not directly measured but rather inferred from the measured eddy currents. In common with most geophysical inversions, this profile is non-unique, implying that many profiles fit the data within error tolerance, and averaging 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 regularization 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. 30 More recent inversion methods include laterally-and spatially-constrained algorithms to regularize 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 andTarantola, 1995, Blatter et al., 2018). The computed posterior probability density function (pdf) 5 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), provided that the constrained depth is consistent across datasets.
In this paper, we derive the implementation of 'MuLTI-TEM' (Multimodal Layered Transdimensional Inversion of Time Domain Electromagnetics) and its use for characterising the subglacial environment. After testing the method on a 10 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 (Vs) profiles, as local accumulations of soft material (partially frozen glacial till) and permafrost overlying bedrock. 15 Finally, we show a combined 3D interpretation of previous shear wave inversions (Killingbeck et al., 2019) and results output from MuLTI-TEM, and suggest the development of a joint resistivity-Vs depth-constrained inversion strategy.

Time Domain Electromagnetics
In TEM surveying, an electromagnetic field is generated by sending a periodic, modified square-wave, current through a 20 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 timedependence induces eddy currents in the subsurface, initiating within the immediate vicinity of the transmitter then spreading downwards. The eddy currents induce a secondary electromagnetic field which propagates back up through the subsurface, 25 inducing a current in a receiver coil located at some distance from the transmitter. The receiver 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. By taking repeated measurements a sounding curve, 30 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) (ℎ), in meters, can be estimated by: where is the transmitter current, is the loop area, is the voltage noise level and is the average resistivity, in Ωm, of the underlying section (Spies, 1989). Equation (1) shows the transmitter loop size is an important acquisition parameter controlling depth of investigation. Large loop soundings (e.g., > 40x40 m), where the receiver coil is located in the centre of 5 the large transmitter coil, have been conducted on thick permafrost regions by Rozenberg and others (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 x 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 section 3.

MuLTI-TEM 10
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 ('Multimodal Layered Transdimensional Inversion), developed for seismic surface wave inversions (see detailed in Killingbeck et al. (2018)). The data input, d, to MuLTI-TEM are the measured voltages ( ) at each timegate ( ), together with an estimate of their uncertainty ( ) derived from the variance of each data point calculated from the stack recordings: 15 The method used to find the posterior distribution of the resistivity profile is outlined below, where (  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 25 which constraints are drawn from high resolution GPR data, we consider depth constraints to be exact since the accuracy of GPR depth estimation is centimetre scale (Killingbeck et al., 2019), compared to the meter scale resolution of the TEM~100times smaller than the thinnest resolvable layer in 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 30 to infinite depth. If no depth constraints are applied then l = 1, corresponding to a half space. We add in an 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, thus k is also an unknown that we determine.
The model vector, that describes the resistivity profile, is then , log ( 1 ), log ( 2 ) … . log ( ), , 1 , 2 … . , log ( 1 ), log ( 2 ) … . log ( )], The prior distribution on the model parameters depend on which layer each nuclei is within. If no depth constraints of the subsurface interfaces are available (that is, there is only a single layer of the half space), the prior distributions on the 10 resistivity is uniform with wide bounds on log(R) (e.g., R between 10 0 -10 5 Ωm), as no prior information (beyond that which can be reasonably assumed for typical materials) is known about the subsurface. However, when depth constraints also provide lithological information, the range of R can be tightened, thus significantly reducing the model parameter space.
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 15 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 Fig. A1 in the appendix for the detailed Leroi input file on which our 20 simplified code is based.
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 25 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.  TEM data were acquired with a Geonics PROTEM 47 system, consisting of a 3 channel digital time-domain receiver, a TEM-47 battery powered transmitter and a 3D multi-turn 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 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 5 (the maximum DOI of each test is estimated using Eq. (1), with voltage noise level 1 nV/m 2 (Fig. 2d)  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 x 10 m square transmitter (Fig. 2b, Fig. 3) was chosen as the optimum survey configuration because: i) it had a fast turn off time (for imaging the shallow surface), 5 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 m to 400 m depending on the underlying resistivity) is sufficient for imaging our target depth, subglacial sediments below ~25 m thick ice. 10

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 5 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 10 (Raiche, 2008), then normally distributed random noise, with a magnitude of 5% of the signal at each timegate, was applied to all time gates, a similar noise model to Blatter et al. (2018)   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 separately run for an unconstrained case, and a case in which the depths of snow-ice and ice-bed horizons are fixed. One 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 figure A2). Additional figures in A3 show that the modal model from the 5 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 one 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 A1. 10 Table 2. Resistivity parameter boundaries used in MuLTI-TEM for the glacier feasibility study.

Material
Resistivity boundaries (Ωm) improvement on their unconstrained equivalents). Note, imaging beneath a conductive structure is difficult for the TEM 25 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 meters.
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 5 (Fig.A6). In this example (Fig. A6 left panel), 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 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 10 is expected. It demonstrates MuLTI-TEM is promising for the potential distributions of resistivity beneath Midtdalsbreen.

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 timegates) or below the background noise level (at later timegates) were removed, shown as "X" in Fig. 7(ii). Inversions were run with 5 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 A1), 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 timegate as this was a good representation of the data variance of each data point calculated from the 3 stack recordings acquired in the field.
Posterior resistivity distributions are shown in Fig. 7(i), comparison of data fit with 200 randomly chosen forward 10 models from the model ensemble are shown in Fig. 7(ii) and posterior distribution of the number of nuclei are shown in Fig.   7(iii). 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. 7(i), 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 ~ ± 10 4 Ωm, whereas that of the 15 conductive subglacial material is ~ ± 10 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 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 Line B and C. In contrast, Line A (off the ice margin, see Fig. 1b) shows a ~70 m-thick 20 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 25 spatially varying pattern of subglacial resistivities from the front of the glacier up to Line C.

2D Resistivity Profiles
MuLTI-TEM 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 timegates) or below the background noise level (at later timegates) were removed, shown as the grey regions in Fig. 8 (left column). The raw signal 10 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 Fig. 8 left column. When using a central loop 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 5 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 10 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 along-and cross-glacier 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 15 in thickness,  a high resistive layer, 10 4 -10 5 Ωm, in line A which matches line D at their intersection, and  an 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 20 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. 6 Interpretation and Discussion

Joint interpretation of MuLTI-TEM 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 matrixwhether sediment or bedrockcannot be determined from resistivity alone. Resistivity is related to the resistivity of the pore fluids divided by the fractional porosity. 5 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, 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 10 modulus, or stiffness, of a material but are insensitive to water content. However, a combined interpretation of resistivity and Vs 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 merit re-evaluation when compared with the co-15 located observation of high resistivity. Therefore, in this integrated interpretation, the high phase velocities are include, thus providing broader bandwidth dispersion curves.
Observing trends of Vs and resistivity in our profiles, three clear patterns emerge within the subglacial material: i) zones of low Vs and low resistivity, ii) high Vs and high resistivity, and 20 iii) high Vs with low resistivity ( Fig. 9; leftmost and centre columns).
These patterns have been used to define 3 different material types (Table 3). From previous electromagnetic and seismic studies (e.g., King et al., 1988;Schneider et al., 2013;Wu at al., 2017), liquid water in the pores of unconsolidated material has a low resistivity and low Vs (we define this as partially frozen sediment), whereas frozen water in pores is very resistive with a high Vs (defined as frozen sediment/permafrost). We assume the bedrock comprises of phyllite, crystalline granite and gneiss with 25 a high Vs 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 magnitudes (Olhoeft, 1981).

Material Vs range (m/s) Resistivity range (Ωm)
Partial frozen sediment/till Vs < 1600 50 < R < 500 Frozen sediment/permafrost Vs > 1900 R > 500 Weathered/fractured bedrock with saline water Vs > 1900 R < 50 The resistivity and Vs profiles are linearly interpolated such that they have mutually consistent sample intervals (1 m) and depth extents (40 m). The TEM and seismic were acquired in the same field season along the same lines, the acquisition 5 parameters were chosen so the two methods could be directly comparable. However, the location of each 1D TEM sounding and seismic common midpoint gather are offset by 2 m, therefore, we linearly interpolated and resampled both the resistivity and Vs solutions, originally sampled every 4 m (Line A, B and C) and 8 m (Line D), to every meter thus making them directly comparable. Killingbeck et al. (2019) consider abrupt lateral variations in the Vs 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 10 Vs 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 Vs 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, 15 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 Vs 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 20 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.

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, with 5 a better fit to the observed data (lower misfit), as demonstrated in the synthetic study.
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 determine separately 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 is 10 sufficient and suitable to detect specific targets.
Although this paper focuses on the specific TEM survey design used in this study, ground based 10 x 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 15 inputs and can therefore be supplied from any external source e.g., borehole measurements or any complementary inference from an external geophysical data source. 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 20 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 would include a fully-coupled joint inversion strategy for resistivity and Vs, which could lead to a more accurate understanding of the subsurface structure (utilizing the structural similarities between resistivity and seismic velocity (e.g., Wisén and Christiansen, 2005). Alternatively, petrophysical relationships could be derived to obtain and/or guide 25 interpretations of the volumetric proportions of water, ice and air in the subsurface (e.g., Hauck et al., 2008). Such a combined approach would lead to more detailed analysis of the Midtdalsbreen margin, 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 inputted to ice-flow models for a better prediction of ice motion over time, and hence future sea level rise. 30

Conclusions
The material properties of the subglacial environment, in particular their water content and saturation, can be characterised by resistivity observations, 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 the inversion algorithm 'MuLTI-TEM', used to overcome such problems. Our method uses a transdimensional Bayesian inversion approach 5 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 and a reduction in misfit. 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. 10 A joint interpretation, using Vs 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 (Vs < 1600 m/s, 50 Ωm < R < 500 Ωm), permafrost (Vs > 1600 m/s, R > 500 Ωm) and weathered/fractured bedrock with saline water in the fractures (Vs > 1900 m/s, R < 50 Ωm). Their spatial extent, within the 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. 15 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 combination, 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-sheets subglacial material can be explored, highlighting the importance of acquiring multiple geophysical datasets for accurately characterising the subglacial environment. 20

Code Availability
MuLTI-TEM can be found at: https://github.com/eespr/MuLTI-TEM, DOI 10.5281/zenodo.3351505. Figure A1. Leroi input file specifically created to match the survey parameters used to acquire TEM data at Midtdalsbreen. The simplified Leroi forward modelling code used in MuLTI-TEM is created based on these parameters. Table A1. Inversion parameters used in MuLTI-TEM for the synthetic feasibility study and 1D and 2D 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).

Inversion Parameter Value
Number of Layers (non-constrained) 1 Number of Layers (constrained) 3 Weighting (data variance, σ) 5% of the signal at each timegate  and support while using the TEM equipment. SK prepared the manuscript with contributions from all co-authors.