https://orcid.org/0000-0001-7407-9741
SUMMARY
Quantitative estimation of pore fractions filled with liquid water, ice and air is crucial for a process-based understanding of permafrost and its hazard potential upon climate-induced degradation. Geophysical methods offer opportunities to image distributions of permafrost constituents in a non-invasive manner. We present a method to jointly estimate the volumetric fractions of liquid water, ice, air and the rock matrix from seismic refraction and electrical resistivity data. Existing approaches rely on conventional inversions of both data sets and a suitable a priori estimate of the porosity distribution to transform velocity and resistivity models into estimates for the four-phase system, often leading to non-physical results. Based on two synthetic experiments and a field data set from an Alpine permafrost site (Schilthorn, Bernese Alps and Switzerland), it is demonstrated that the developed petrophysical joint inversion provides physically plausible solutions, even in the absence of prior porosity estimates. An assessment of the model covariance matrix for the coupled inverse problem reveals remaining petrophysical ambiguities, in particular between ice and rock matrix. Incorporation of petrophysical a priori information is demonstrated by penalizing ice occurrence within the first two meters of the subsurface where the measured borehole temperatures are positive. Joint inversion of the field data set reveals a shallow air-rich layer with high porosity on top of a lower-porosity subsurface with laterally varying ice and liquid water contents. Non-physical values (e.g. negative saturations) do not occur and estimated ice saturations of 0–50 per cent as well as liquid water saturations of 15–75 per cent are in agreement with the relatively warm borehole temperatures between −0.5 and 3 ° C. The presented method helps to improve quantification of water, ice and air from geophysical observations.
1 INTRODUCTION
Climate-induced degradation of permafrost can release substantial amounts of soil organic carbon into the atmosphere (e.g. Schuur et al. 2015) and increase the probability of slope failures in Alpine regions (e.g. Huggel et al. 2012). Understanding hydrological processes in permafrost systems is crucial to parametrize numerical models that simulate the evolution and potential carbon feedback of terrestrial permafrost as well as to assess the hazard potential of permafrost degradation on a physical basis.
While borehole information is expensive and limited to discrete locations, geophysical imaging offers opportunities to derive quantitative and non-invasive insights on permafrost characteristics at high spatial and temporal resolution. Electrical resistivity and acoustic velocity of a medium are sensitive to the phase change of water between its liquid, frozen, and gaseous states. Electrical resistivity tomography (ERT) and refraction seismic tomography (RST) are thus widely used in cryospheric geophysical applications (Hauck & Kneisel 2008).
Hilbich et al. (2008) analysed a 7-yr-long ERT and borehole temperature monitoring data set at the Schilthorn, Swiss Alps and characterized both short-term (e.g. seasonal active layer dynamics) and long-term effects (e.g. ground ice degradation as a consequence of the extraordinary hot European summer in 2003). Dafflon et al. (2016) combined ERT with frequency-domain electromagnetic induction data, core analysis, and digital surface models to estimate the spatial distribution of shallow permafrost in an Alaskan tundra environment. Oldenborger & LeBlanc (2018) imaged changes in unfrozen water content in accordance with temperature measurements alongside airport infrastructure.
While ERT is highly sensitive to unfrozen water content, quantitative estimates of ice from electrical data alone are impeded since ice and air are practically electrical insulators and may not be distinguished from a resistive rock matrix. Fortunately, air and ice usually differ by an order of magnitude in their acoustic velocities, in addition to velocity changes of up to 2000 m s–1 between frozen and unfrozen water (Hilbich 2010), making RST a valuable additional method in cryospheric geophysics (e.g. Harris & Cook 1986; Krautblatter & Draebing 2014; Steiner et al. 2019). For example, Dou & Ajo-Franklin (2014) performed seismic measurements in the arctic permafrost region and identified a low velocity zone, which is likely related to partially thawn saline permafrost and in agreement with ERT measurements by Hubbard et al. (2013). Merz et al. (2016) used ERT and RST measurements together with other geophysical as well as geotechnical data for a multidisciplinary characterization of an Alpine rock glacier. The authors combined the insights obtained from the different methods in a qualitative interpretation and concluded that joint inversions of different geophysical data sets offer opportunities for the quantification of permafrost composition in future studies.
To take advantage of the complementary sensitivities of ERT and RST for the quantification of permafrost constituents, Hauck et al. (2011) presented a petrophysical four-phase model (4PM) incorporating estimates of liquid water, ice and air contents for given separate inversions of electrical and seismic refraction data sets and a pre-existing porosity estimate. Mewes et al. (2017) assessed the resolution capacity of the 4PM, which depends on the individual resolution capacities of ERT and RST, with regard to the estimation of liquid water and ice contents for typical processes in the context of permafrost degradation. The authors emphasized that artifacts in the individual inversions can impair the physical plausibility of the estimated constituents and potentially lead to misinterpretation. They further found that smoothness regularization applied to the individual inversions can underestimate the magnitude of changes. Pellet et al. (2016) improved the 4PM through soil moisture calibration. Instead of prescribing the porosity, the authors developed an estimate of its distribution based on a three-phase model in the unfrozen part of the subsurface. They found that one remaining challenge in the application of the 4PM is the possibility of simulated ice occurrence even in regions where it is highly unlikely, highlighting the need for means to impose physical constraints during inversion. To this end, we present an approach that uses apparent resistivities and seismic traveltimes simultaneously in a petrophysically coupled joint inversion to improve the quantification of permafrost constituents from geophysical observations.
2 THEORY AND METHOD
2.1 Acoustic and electrical properties of partially frozen ground
Schematic on the estimation of water, ice and air from ERT and RST data. (a) Conventional inversion of both data sets with subsequent petrophysical transformation. (b) Petrophysical joint inversion honoring both data sets and petrophysical relations during parameter estimation.
2.2 Petrophysical joint inversion for permafrost constituents
The first term quantifies the misfit between observed data |$\boldsymbol {d}$| and the model response |$\boldsymbol {\mathcal {F}}(\boldsymbol {m})$| incorporating the reciprocals of the individual data errors on the diagonal of data weighting matrix |$\boldsymbol {W}_\text{d}$|. Note that the model response |$\boldsymbol {\mathcal {F}}(\boldsymbol {m})$| contains petrophysical transformation according to eqs (2) and (3) followed by independent solutions of the RST and ERT forward problems. The second term represents a smoothness regularization applied to the model vector |$\boldsymbol {m}$|, where |$\boldsymbol {W}_\text{m}$| is a block matrix holding four first-order roughness operators on its diagonal to promote smoothness in the distribution of each constituent of the four-phase system. The third term is an additional regularization term to fulfill the volume conservation constraint in eq. (1). Here, |$\boldsymbol {W}^\text{sum}_\text{p}$| is a block matrix of four adjacent identity matrices acting on the untransformed petrophysical parameter vector |$\boldsymbol {p}$| to penalize solutions for which the sum of the four volumetric fractions deviates from unity. The fourth term represents a damping regularization and allows to incorporate a priori information on the petrophysical target parameters by penalizing deviations from a given reference model |$\boldsymbol {p}_\text{0}$|. Here, |$\boldsymbol {W}_\text{p}$| is a square matrix with either zeros or ones along its diagonal depending on which model parameters are sought to stay close to the reference model |$\boldsymbol {p}_\text{0}$|. The fourth term is optional. For example, we discuss joint inversion results with and without a prescribed porosity distribution throughout this paper. The former uses γ = 0, whereas the latter uses γ = β and |$\boldsymbol {W}_\text{p} = \text{diag}([\boldsymbol {0},\boldsymbol {0},\boldsymbol {0},\boldsymbol {I}])$| to penalize solutions for which the rock content distribution deviates from its prior estimate. The dimensionless factors α and β scale the influence of the regularization terms. β is chosen large enough to prohibit non-physical solutions, while α is chosen to fit the data within error bounds.
3 SYNTHETIC EXAMPLES
3.1 Synthetic model and data
To evaluate the performance of the joint inversion approach, a three-layer model is considered (Fig. 2a). Parameters are defined in terms of water, ice, air and rock contents and subsequently transformed into velocity and resistivity distributions for the generation of synthetic data. The model represents a typical layered scenario encountered in Alpine periglacial environments comprising a 5 -m-thick, unfrozen active layer (i.e. the seasonal thaw layer), a 10 -m-thick partially thawn layer with laterally changing ice and liquid water contents, and a frozen bedrock underneath. The porosity is decreasing from 40 to 20 per cent layer-wise with depth. Velocity and electrical resistivity are calculated according to eqs (2) and (3) with the parameters listed in Table 1 and generally increase with depth due to the decreases in liquid water and air contents, combined with increases in ice and rock contents.
(a) True model, (b) conventional inversion results and (c) joint inversion results with a priori knowledge of the porosity distribution. All models are cut off below the lowermost ray path and in regions where the respective volumetric fraction is negative. Sensors are marked as black semicircles. Note that the electrical resistivity (second row) is displayed on a logarithmic colour scale. If not otherwise indicated in the lower left of each panel, the distributions shown in (b) and (c) directly result from the respective inversions. The annotation Transformed means that the respective quantity is obtained through transformation of the actual inversion results using petrophysical equations.
| Archie parameters | Constituent velocities | ||||
|---|---|---|---|---|---|
| Parameter | Value | Unit | Parameter | Value | Unit |
| ρw | 150 | Ωm | vw | 1500 | m s–1 |
| n | 2 | – | vi | 3500 | m s–1 |
| m | 1.3 | – | va | 330 | m s–1 |
| vr | 5500 | m s–1 | |||
| Archie parameters | Constituent velocities | ||||
|---|---|---|---|---|---|
| Parameter | Value | Unit | Parameter | Value | Unit |
| ρw | 150 | Ωm | vw | 1500 | m s–1 |
| n | 2 | – | vi | 3500 | m s–1 |
| m | 1.3 | – | va | 330 | m s–1 |
| vr | 5500 | m s–1 | |||
| Archie parameters | Constituent velocities | ||||
|---|---|---|---|---|---|
| Parameter | Value | Unit | Parameter | Value | Unit |
| ρw | 150 | Ωm | vw | 1500 | m s–1 |
| n | 2 | – | vi | 3500 | m s–1 |
| m | 1.3 | – | va | 330 | m s–1 |
| vr | 5500 | m s–1 | |||
| Archie parameters | Constituent velocities | ||||
|---|---|---|---|---|---|
| Parameter | Value | Unit | Parameter | Value | Unit |
| ρw | 150 | Ωm | vw | 1500 | m s–1 |
| n | 2 | – | vi | 3500 | m s–1 |
| m | 1.3 | – | va | 330 | m s–1 |
| vr | 5500 | m s–1 | |||
While the acoustic velocities of water, ice, and air are agreed upon in the literature, vr is strongly dependent on the type of rock (or soil). The chosen value of 5500 m s–1 may represent a metamorphic rock such as a pyritic paragneiss (e.g. Draebing & Krautblatter 2012). Common literature values are used for the Archie parameters (m and n), whereas the pore water resistivity ρw was chosen on the basis of laboratory measurements in an Alpine permafrost context (e.g. Hauck & Kneisel 2008). All parameters in Table 1 are chosen to be constant throughout the model domain and not estimated during inversion. We refer to Hauck et al. (2011) for a sensitivity analysis of vr, ρw, m and n in the context of ice estimation.
Synthetic traveltime data are generated assuming 53 geophones spaced by 2.5 m and collocated shot positions resulting in 2756 shot-receiver pairs. Ray tracing utilizes the shortest path method with three additional nodes on each edge of a triangular model cell to increase the ray path accuracy following the approach outlined by Giroux & Larouche (2013). Traveltime noise was added as additive Gaussian white noise (AGWN) with a standard deviation of 0.5 ms.
Electrode positions coincide with geophone locations. A dipole–dipole data set with dipole lengths of one, two and four unit electrode spacings is simulated neglecting absolute geometric factors above 5000 m resulting in 1414 apparent resistivities. Forward modelling employs quadratic shape functions, an enlarged forward modeling domain to avoid boundary effects, and is detailed by Rücker et al. (2006). A normally distributed relative error of 5 per cent was added to the simulated apparent resistivities. All conventional and joint inversion results presented in the following use four times larger smoothing in the horizontal direction to promote the expected layered structure and describe the synthetic data sets within their respective error bounds.
3.2 Inversion with correct porosity estimate
Conventional inversion results of traveltimes and apparent resistivities (Fig. 2b) are obtained using a Gauss–Newton scheme with smoothness regularization as detailed in Rücker et al. (2017). Velocity and electrical resistivity tomograms are transformed into distributions of water, ice and air using the 4PM and assuming that the porosity structure is known.
The different volumetric fractions are qualitatively and quantitatively in good agreement with the true model. Solely in the conductive and low velocity top layer, where the individual inversions exhibit small-scale variability close to the sensors, the 4PM produces small regions of negative ice content down to –6 per cent and overestimates the water content by up to 9 per cent.
Fig. 2(c) shows the results obtained with the developed petrophysical joint inversion approach. The quantitative agreement to the true model (Fig. 2a) is slightly better in comparison to the conventional inversion and the layer boundaries are more pronounced. Moreover, non-physical values do not occur close to the sensors.
3.3 Inversion with incorrect porosity estimate
We evaluate the inversion performance without detailed knowledge of the porosity distribution, as common for field applications (Fig. 3). Note that the true model has not changed, that is Figs 2(a) and 3(a) show identical distributions. For the conventional inversion, a homogeneous rock matrix content of 70 per cent (i.e. ϕ = 0.3) is assumed (lowermost panel in Fig. 3b). This estimate is correct for the middle layer, but under- and overestimates the porosities of the top and bottom layer by 10 per cent, respectively, leading to non-physical ice content estimations in the top layer (reaching –14 per cent). In turn, ice contents in the bottom layer are strongly overestimated. This overestimation stems from the assumption of homogeneous porosity, as the high velocity below 15 m depth can no longer be explained by a porosity decrease and is compensated by additional ice. Furthermore, the air content in this region is non-physical (i.e. slightly below zero).
As Fig. 2 but without a priori knowledge of the porosity distribution.
In the corresponding joint inversions, homogeneous rock content was only used as the starting model for fr and allowed to vary by ±15 per cent during the inversion. Estimates of water, ice and air contents are significantly improved through the joint inversion approach and do not exhibit any non-physical values (Fig. 3c). Moreover, the inversion indicates a porosity decrease with depth revealing that the measured data cannot be explained with the homogeneous starting model of fr. The high ice content in the centre of the model cannot be reconstructed and is compensated by an increase in rock content.
3.4 Model parameter interdependency
To quantify the interdependency of water, ice, air and rock content, we consider the model covariance matrix of the coupled inverse problem. The model covariance elucidates how errors in the data are propagated into errors of the estimated model parameters. Diagonal elements represent variances of the corresponding parameters, while off-diagonal values indicate correlations between pairs of model parameters (e.g. Aster et al. 2012). For better illustration, the parameters are grouped into five discrete blocks as outlined in Fig. 3(a).
To highlight the need for petrophysical constraints, we compare the covariances of the coupled inverse problem without (Fig. 4a) and with volume conservation constraints (Fig. 4b). In the unconstrained covariance matrix, the air content has the lowest diagonal values, partly due to generally small air contents in the model. In addition, the velocity of air differs significantly from the other constituent velocities (Table 1) leading to its generally good discriminability. For the remaining fractions, strong off-diagonal elements appear particularly in regions where the amount of unfrozen water is relatively high (parameters i, ii and iv). In comparison, the covariances are considerably lower in regions where the water content is low and no air exists (parameters iii and v). Covariances within one parameter group and between parameter groups are drastically reduced after applying volume conservation constraints in eq. (1) (Fig. 4b). Significant variances remain in the upper unfrozen part of the model (parameter i), where sensitivities are generally higher, indicating that errors in the data are more strongly reflected in errors in the model. In addition, strong variances are visible between ice and rock contents for the other parameters, that is, at greater depth.
Model covariance matrix for the petrophysical joint inversion. Red and black indicate strong positive and negative parameter covariance, white indicates no correlation. The Roman numbering of the model parameters is depicted in Fig. 3(a).
4 FIELD DATA EXAMPLE
To demonstrate the applicability of the developed joint inversion approach to a real permafrost scenario, we consider a data set acquired at the Schilthorn site located in the Bernese Alps in Switzerland. The lithology of the Schilthorn massif mainly consists of ferruginous sandstone schists (Imhof et al. 2000). The bedrock weathering produced a fine-grained surface debris layer with a thickness up to 5 m (Hilbich et al. 2008). As part of the PACE project, a drilling campaign and first RST and ERT measurements were performed in 1998 to assess the presence of permafrost (Vonder Mühll et al. 2000). Borehole temperature measurements revealed warm permafrost, that is, near subzero temperature, where the liquid water content can be relatively high. The Schilthorn site became a reference monitoring site of PERMOS (Swiss Permafrost Monitoring Network, PERMOS 2016) including numerous measurements of temperature in three boreholes, ground surface temperature, soil moisture, ERT, RST, snow thickness, wind speed and direction, as well as solar radiation.
An ERT monitoring profile was initiated in 1999, automatized in 2009, extended in 2012 and is still maintained and functional (for data acquisition and preprocessing details see Hilbich et al. 2008, 2011; Mollaret et al. 2018). The ERT acquisition system is composed of 49 fixed electrodes spaced by 2 m connected to a GeoTom device (GEOLOG, Germany) used with a Wenner configuration. We only consider the first half of this permanent profile, where collocated geophone positions exist to ensure equal lateral coverage of both methods. The seismic signal is generated by a sledge hammer striking a steel plate at 25 shot positions located in between each geophone and recorded through a Geode device (Geometrics, USA). First breaks were manually picked and seismic traveltime errors were estimated to range between 0.1 and 0.5 ms (Hilbich 2010). We use 0.3 ms as an error estimate for the RST data and 3 per cent relative error for the ERT data. Similarly as for the inversion of synthetic data, all conventional and joint inversion results use four times larger smoothing in the horizontal direction to promote layered subsurface structures.
We consider a data set measured on 19 August 2014 presented by Pellet et al. (2016). On that day, the thaw layer depth in the two boreholes available along the profile was observed around 2.1 m depth. In warm recent years, maximum thaw depths at the end of summer reached 9–10 m indicative of considerable permafrost degradation (PERMOS 2019). Fig. 5(a) shows the conventional inversion results with a subsequent application of the 4PM using the parameters from Pellet et al. (2016) as listed in Table 2 and a constant porosity of 53 per cent based on in situ measurements performed by (Scherler 2006). To allow comparability with previous Schilthorn studies (e.g. Hilbich et al. 2008; Pellet et al. 2016), water, ice and air contents are displayed as saturations, that is divided by porosity.
Tomograms of the Schilthorn field data sets obtained through (a) conventional inversion, (b) joint inversion and (c) joint inversion with borehole constraints. If not otherwise indicated in the lower left of each panel, the shown distributions directly result from the respective inversions. All models are cut off below the lowermost ray path and in regions where the respective saturation is negative or exceeds one. Sensors are marked as black circles. Note that the electrical resistivity (second row) is displayed on a logarithmic colour scale. Boreholes and associated thaw depths are marked as vertical and horizontal lines, respectively. The black boxes superimposed on the ice saturation in (c) mark the regions where the ice content has been constrained to zero.
| Archie parameters | Constituent velocities | ||||
|---|---|---|---|---|---|
| Parameter | Value | Unit | Parameter | Value | Unit |
| ρw | 60 | Ωm | vw | 1500 | m s–1 |
| n | 2.4 | – | vi | 3500 | m s–1 |
| m | 1.4 | – | va | 300 | m s–1 |
| vr | 6000 | m s–1 | |||
| Archie parameters | Constituent velocities | ||||
|---|---|---|---|---|---|
| Parameter | Value | Unit | Parameter | Value | Unit |
| ρw | 60 | Ωm | vw | 1500 | m s–1 |
| n | 2.4 | – | vi | 3500 | m s–1 |
| m | 1.4 | – | va | 300 | m s–1 |
| vr | 6000 | m s–1 | |||
| Archie parameters | Constituent velocities | ||||
|---|---|---|---|---|---|
| Parameter | Value | Unit | Parameter | Value | Unit |
| ρw | 60 | Ωm | vw | 1500 | m s–1 |
| n | 2.4 | – | vi | 3500 | m s–1 |
| m | 1.4 | – | va | 300 | m s–1 |
| vr | 6000 | m s–1 | |||
| Archie parameters | Constituent velocities | ||||
|---|---|---|---|---|---|
| Parameter | Value | Unit | Parameter | Value | Unit |
| ρw | 60 | Ωm | vw | 1500 | m s–1 |
| n | 2.4 | – | vi | 3500 | m s–1 |
| m | 1.4 | – | va | 300 | m s–1 |
| vr | 6000 | m s–1 | |||
With a prescribed and constant porosity, the conventional approach can only satisfy the low velocities near the surface with non-physical values for ice and air saturations reaching –142 and +207 per cent, respectively (Fig. 5a). In contrast, non-physical values do not occur in the petrophysical joint inversion results, as lower porosities near the surface (Fig. 5b) are allowed during parameter estimation.
Since the joint inversion is formulated in terms of the petrophysical target parameters, it becomes possible to include prior information, for example, constraints on water content based on soil moisture measurements or temperature-dependent constraints on ice occurrence. Here, we demonstrate the latter by constraining ice contents to zero, where the measured borehole temperatures (shown in Figs 6a and 7a) are positive. This is achieved by creating a box with a radius of 2 m around the corresponding depth interval of the borehole (as indicated for ice saturation in Fig. 5c) and adding cells in this box to the damping constraint in eq. (6). To allow for sharp transitions in the vertical direction across the known thaw depth, we disable smoothing regularization across the lower boundary of this box (by setting the corresponding entries of |$\boldsymbol {W}_\text{m}$| to zero). A similar approach has been used by Wunderlich et al. (2018) to constrain ERT inversion results to direct push electrical conductivity logs.
(a) Temperature measured on 19 August 2014 in borehole SCH_5198 and (b–e) corresponding four-phase constituents derived from conventional inversion and petrophysical joint inversion without and with borehole constraints. The horizontal and vertical grey lines mark the thaw depth and zero line, respectively. The borehole location is depicted in the lower left panel of Fig. 5.
(a) Temperature measured on 19 August 2014 in borehole SCH_5000 and (b–e) corresponding four-phase constituents derived from conventional inversion and petrophysical joint inversion without and with borehole constraints. The horizontal and vertical grey lines mark the thaw depth and zero line, respectively. The borehole location is depicted in the lower left panel of Fig. 5.
While generally comparable to the unconstrained version (Fig. 5b), no ice appears in the upper 2 m around the boreholes in the joint inversion result with borehole constraints (Fig. 5c). The decrease in ice also affects deeper parts of the model and is compensated by an increase in rock content (i.e. a decrease in porosity), which in turn leads to higher water saturations. The air saturation in all three approaches is high in the thawed part of the model and quickly approaches zero at larger depths.
To highlight quantitative differences between the different inversion approaches, Figs 6 and 7 show the measured borehole temperatures with depth together with the inversion results shown in Fig. 5 extracted at the borehole locations. The effect of logarithmic barriers can be seen in Fig. 6(c) for instance, where in contrast to the conventional approach, joint inversion results do not cross the zero line.
While the conventional inversion and joint inversion without borehole constraints show smooth transitions of ice saturation in the upper few meters for example, the added borehole constraints result in sharp transitions, that is, a better delineation of the thawn layer (e.g. Fig. 7c). Again, it can be seen how the decrease in ice (e.g. Fig. 7c) is causing an increase in rock content (e.g. Fig. 7e), which in turn increases water saturation (e.g. Fig. 7b).
5 DISCUSSION
Transformation of conventionally inverted velocity and resistivity tomograms into estimates of water, ice and air can lead to non-physical results in the presence of data errors and incorrect porosity estimates. For the Schilthorn case, this was also found by Pellet et al. (2016), who then derived porosity in the unfrozen part from a three-phase model calibrated to shallow soil moisture measurements and assumed a gradient porosity in the deeper part. With the assumption of porosity decrease with depth and soil moisture measurements, Pellet et al. (2016) obtained tomograms comparable to our results from the petrophysical joint inversion, which used no prior information on porosity.
Petrophysical joint inversion combines the information of RST and ERT measurements and leads to quantitatively improved images by honoring petrophysical relations and volume conservation during parameter estimation. Inversion of synthetic data without a porosity estimate (Fig. 3c) and analysis of the model covariances has revealed a strong ambiguity between ice and rock contents. This ambiguity also became apparent in the application to field data and is in agreement with the findings by Hauck et al. (2011). The authors analytically explored the range of possible values for ice, water, air, and rock contents for a given pair of resistivity and velocity values. By comparing the spread of solutions, it was found that air and water content can be discriminated quite well even if porosity is unknown, while there is a strong ambiguity between ice and rock contents.
The discriminability between ice and rock matrix could be improved by incorporation of additional freeze–thaw sensitive data sets such as complex electrical resistivity measurements. Results of first laboratory (Wu et al. 2013; Kemna et al. 2014) and field studies (Grimm & Stillman 2015; Mudler et al. 2019) hold promise for ice quantification, as the spectral response of frozen ground is strongly affected by the electrical polarization characteristics of ice. Ambiguity could further be reduced in a monitoring context, where the porosity can be assumed to be constant within the observed period (Hauck et al. 2017). Inclusion of multiple timesteps into a time-lapse joint inversion therefore represents a promising extension to this work.
The approach presented is not restricted to the empirical model by Archie (1942) and the time-averaging equation by Timur (1968) and does not attempt to address their general shortcomings. For example, a general and often overlooked problem is that the empirical factors (e.g. the saturation exponent) are commonly assumed to be spatially constant in field applications. Future extensions of our work could incorporate advanced petrophysical formulations more representative for the studied field sites.
For the case of saline permafrost for example, Wu et al. (2017) presented a modified formulation, which takes structural soil changes as well as temperature-dependent salinity changes (i.e. electrical conductivity changes) of unfrozen water during freeze-thaw transitions into account. Based on a field survey on a rock glacier and associated laboratory measurements, Duvillard et al. (2018) emphasized that surface conduction can be significant and has to be taken into account when estimating liquid water content in environments where electrolytic conduction is not dominating the bulk electrical conductivity. For the case of low-porosity hard rocks, Draebing & Krautblatter (2012) presented a modified Timur equation that accounts for changes in matrix velocity due to ice pressure. Dou et al. (2017) presented a two-end-member mixing approach accounting for the coexistence of frame-strengthening and pore-filling ice to describe the P-wave velocity in saturated, unconsolidated saline permafrost, where conventional slowness averaging would be inadequate.
6 CONCLUSIONS
We have developed a petrophysical joint inversion approach that uses seismic traveltimes and apparent resistivities to image distributions of water, ice, air and rock content. Since petrophysical relations and volume conservation are honored during parameter estimation, our approach produces physically meaningful results, even in the absence of correct porosity estimates, and thereby outperforms post-inversion transformation of conventional tomograms. A significant advantage is that the inversion constraints can be formulated in terms of the petrophysical target parameters facilitating the flexible use of a priori information and direct incorporation of non-geophysical data, for example, constraints on water content inferred from soil moisture measurements, into the inversion.
An application to a field data set from the Schilthorn, Swiss Alps, has revealed physically plausible tomograms in agreement with previous studies without relying on suitable a priori porosity estimates. We conclude that our method contributes to improved quantification of water, ice and air from geophysical observations and will therefore be of direct use for researchers and practitioners in cryogeophysical and hydrogeophysical applications.
Yet joint inversion alone is not able to overcome the inherent petrophysical ambiguities between ice and rock matrix, which remain to be addressed in future studies through additional (non-)geophysical observations, advanced petrophysical formulations and monitoring applications. To facilitate adoption and further development of the method, we made the algorithm available under a permissive BSD license. It is based on the open-source modeling and inversion library pyGIMLi (Rücker et al. 2017). The implementation as well as scripts to reproduce results and figures of this paper can be found at https://github.com/florian-wagner/four-phase-inversion.
ACKNOWLEDGEMENTS
The first author gratefully acknowledges funding received by the Dr. Erich Ritter foundation in cooperation with the Water Science Alliance, which allowed preparation of this paper during a research visit at Lawrence Berkeley National Laboratory inspired by many fruitful discussions with Baptiste Dafflon, Sebastian Uhlemann and other colleagues in the Earth and Environmental Sciences Area. The constructive comments by editor Jörg Renner, Adam Booth and three anonymous reviewers helped to improve earlier versions of this manuscript. We thank Carsten Rücker for his dedication to the development of pyGIMLi.
