Resolution capacity of geophysical monitoring regarding permafrost degradation induced by hydrological processes

Geophysical methods are often used to characterise and monitor the subsurface composition of permafrost. The resolution capacity of standard methods, i.e. Electrical Resistivity Tomography and Refraction Seismic Tomography, depends hereby not only on static parameters such as measurement geometry, but also on the temporal variability in the contrast of the 10 geophysical variables (electrical resistivity and P-wave velocity). Our study analyses the resolution capacity of Electrical Resistivity Tomography and Refraction Seismic Tomography for typical processes in the context of permafrost degradation using synthetic and field data sets of mountain permafrost terrain. In addition, we tested especially the resolution capacity of a petrophysically-based quantitative combination of both methods, the so-called 4-phase model, and by this analysed the expected changes in water and ice content upon permafrost thaw. The results from the synthetic data experiments suggest a 15 higher sensitivity regarding an increase in water content compared to a decrease in ice content, and a potentially larger uncertainty originating from the individual geophysical methods than from the combined evaluation with the 4-phase model. In the latter, a loss of ground ice can be detected quite reliably, whereas artefacts occur in the case of increased horizontal or vertical water flow. Analysis of field data from a well-investigated rock glacier in the Swiss Alps successfully visualised the seasonal ice loss in summer, and the complex spatially variable ice-, waterand air content changes in an interannual 20 comparison.


Introduction
Geophysical methods are frequently applied in investigations of mountain permafrost (Hauck and Kneisel, 2008).In particular electrical resistivity tomography (ERT) and refraction seismic tomography (RST) are well suited to differentiating between frozen and unfrozen conditions (for reviews see Scott et al., 1990;Kneisel et al., 2008;Hauck, 2013).In monitoring applications, ERT and RST allow phase changes to be resolved between liquid water and ice as well as saturated and unsaturated conditions over a wide range of temporal and spatial resolutions (e.g.Hilbich et al., 2008Hilbich et al., , 2011;;Maurer and Hauck, 2007;Krautblatter et al., 2010;Hilbich, 2010;Kneisel et al., 2014;Supper et al., 2014;Dräbing, 2016).
Among the various types of landforms in which mountain permafrost occurs, rock glaciers occupy a specific place, partly due to their role as permafrost indicators, and can be easily detected in the landscape due to their kinematic behaviour, with movement rates typically ranging from a few centimetres up to several metres per year.The increase in creep velocities observed for many rock glaciers in the European Alps for the last 20-30 years, and in particular for the last decade (e.g.Kellerer-Pirklbauer and Kaufmann, 2012;Delaloye et al., 2008;PERMOS, 2016), makes the investigation of the role of liquid water for the dynamics of rock glaciers increasingly important (e.g.Haeberli et al., 2006;Jansen and Hergarten, 2006;Krainer and Mostler, 2002).Thermal and hydro-mechanical effects of water infiltration may significantly control the seasonal and short-time pattern of creep velocities, e.g. by the reduction of shear resistance by infiltrated meltwater (e.g.Krainer and He, 2006;Perruchoud and Delaloye, 2007;Ikeda et al., 2008).

B. Mewes et al.: Resolution capacity of combined geophysical monitoring
Whereas water infiltration at the surface of rock glaciers could be estimated by means of indirect interpretation of ground thermal regimes (e.g.Humlum, 1997;Sawada et al., 2003;Staub and Delaloye, 2016), in-situ measurements are rare, in particular because a standard sensor deployment is often difficult due to the coarse-grained surface layer (Buchli et al., 2013).To our knowledge, direct measurements of water infiltration into the permafrost body are nonexistent so far.Indirect methods such as geophysical techniques (e.g.ERT and RST, but also nuclear magnetic resonance (NMR) and self-potential (SP); see Hauck, 2013) could in principle help to estimate the liquid water content within the permafrost body and to identify preferential flow paths (Scapozza et al., 2008;Scherler et al., 2010;Lehmann-Horn et al., 2011;Voytek et al., 2016).They can also efficiently be used to monitor permafrost changes over time (e.g.Hilbich et al., 2011).However, the resolution capacity of geophysical data for localized anomalies, such as preferential flow paths of meltwater through permafrost and talik formation, depends on various factors such as depth of the anomaly, resistivity and velocity contrasts and measurement geometry.A reliable interpretation of small-scale geophysical anomalies (or their temporal changes) remains, therefore, difficult (Hilbich et al., 2009).Current challenges for the application of geophysical methods to rock glaciers include (a) the resolution potential for small-scale anomalies and processes (ideally down to the depth of the shear horizon at 15-20 m or more), and (b) the delineation of the complex internal structure and the quantification of the overall ice/water content in a rock glacier (Maurer and Hauck, 2007;Monnier et al., 2011;Monnier and Kinnard, 2013;Hausmann et al., 2007;Merz et al., 2015).
The aim of this study is to analyse the qualitative and quantitative resolution potential of repeated geophysical measurements (ERT and RST) for the monitoring of hydrological and glaciological processes with a focus on snowmelt infiltration and ice melt.Here we focus on standard geometries that are commonly used in mountain permafrost monitoring applications such as in the context of permafrost degradation within rock glaciers.We set up four synthetic models with simplified geometries simulating idealized water flow paths and ice loss scenarios.After inverting the synthetic data sets to yield tomograms of specific resistivity and seismic P-wave velocity, we combine the results from both methods to run a petrophysical model, which estimates the quantitative distribution of ice and water contents within the rock glacier (Hauck et al., 2011).With this so-called 4-phase model (4PM) and the simplified scenarios of changes in ice and water content, we investigate the performance of the 4PM regarding its spatial resolution capacity, which is partly a function of the individual resolution capacities of ERT and RST tomograms.Because resolution capacity depends first of all on acquisition geometry (sensor spacing), the synthetic experiments are conducted for different spacings including the rather coarse standard spacing of 4-5 m used on many rock glaciers worldwide to achieve large penetration depths.
Finally, the results will be analysed in comparison with real field data from the Becs de Bosson rock glacier, a wellinvestigated rock glacier in the Swiss Alps (e.g.Tenthorey, 1993).Seasonal and interannual ERT and RST monitoring data with coarse spacings during a 2-year period from 2012 to 2013 are used to model the volumetric fractions of ice, water and air with the 4PM to illustrate that the resolution potential regarding hydrological processes in the context of ice content changes can be quite small for standard coarsescale monitoring set-ups.

Field site and available data
The Becs de Bosson rock glacier (see Fig. 1) is the biggest active rock glacier in the Val de Réchy, western Swiss Alps.It covers an area of 246 000 m 2 in NW exposition over an altitudinal range between 2620 and 2850 m a.s.l.The rock glacier mainly consists of calc-schist debris, whereas the landform extends over a basement of quartzite and gneiss locally superimposed by evaporitic rocks (dolomite, gypsum) (Tenthorey, 1993).The measured mean annual air temperature on the rock glacier is −1 • C (1997-2014), which is approximately 1-1.5 • C warmer than in the decades prior to the 1980s (according to available long-term time series in the Swiss Alps, http://www.meteoswiss.ch).Most of the precipitation falls as snow from October to May, constituting a snowpack that frequently exceeds 2 m before the beginning of the melt season in early May.In early July, the snow has usually completely disappeared on most of the rock glacier, while some snow patches may persist in the rooting zone until end of summer.
Large parts of the rock glacier rooting zone contain no permafrost due to the prior extension of a glacier during the Little Ice Age (Mari et al., 2013;Tenthorey, 1993;Delaloye, 2004).The rock glacier tongue is still frozen and affected by permafrost creep.The active layer is estimated to be about 4 m thick along the frozen section and consists of small debris mixed with big boulders (Tenthorey, 1993;Delaloye, 2004).Creep rates are very heterogeneous but gradually decrease downwards from about 2.4 m a −1 (in 2012-2013) in the central part to 0.3 m a −1 near the rock glacier tongue.After the average creep velocity had already doubled between the 1980s and the early 2000s (Kääb, 2005;Perruchoud and Delaloye, 2007), GNSS surveys revealed a continuous velocity increase in recent years for selected points along the ERT/RST profile (Fig. 1).Seasonal creep velocity analysis revealed an acceleration during the time of snowmelt (Perruchoud and Delaloye, 2007) and a drop in creep velocity in July followed by an acceleration phase in late summer.
Repeated ERT and RST surveys are conducted along a longitudinal profile reaching from the rooting zone (no permafrost) to across the rock glacier tongue (see Fig. 1) with the aim to investigate the change of its spatial ground ice distribution and the influence that ice and liquid water content may have on its kinematic behaviour.For this study, coincident ERT and RST data from 21 July 2012, 25 July 2013 and 27 September 2013 were used.Note that the ERT profile is longer here (315 m) than the RST profile (189 m).As the 4PM (see Sect. 4) approach needs ERT and RST data acquired along the same horizontal grid as input for the calculation of ice, water and air contents, we focused our analysis on the part of the profile where both data sets are available, i.e. the length of the seismic profile.
ERT data were collected with a GeoTom instrument (Geolog, Augsburg) along 71 permanently installed electrodes with 4.5 m spacing using the Wenner and Wenner Schlumberger configuration.Galvanic coupling was good with contact resistances < 50 m for every electrode.However, data misfit (absolute error in the inversion) was still between 6.9 and 14 % after four iterations, illustrating the difficult measurement conditions on a rock glacier (contact resistances) as well as the large and spatially variable resistivity contrasts within the rock glacier.They are nevertheless comparable to similar surveys on rock glaciers presented in previous studies (e.g, Marescot et al., 2003;Hilbich et al., 2009;Maurer and Hauck, 2007).The RST data were measured with a Geode24 instrument (Geometrics) in roll-along mode using 43 receiver locations with 4.5 m spacing and 28 shot locations with ∼ 9 m spacing and including offset shot locations at ±4 m located between every second geophone pair as well as four off-end locations at both ends of the profile.Generally, about 15-20 individual shots were stacked at each shot point location to ensure good data quality for first-arrival picking.Shot points and geophone locations were marked and measured by dGPS to reduce the geometric error by the repeated surveys.The overlapping part of the ERT and RST profiles used for the 4PM calculation yields an overall length of 160 m.

Geophysical monitoring and forward-inverse modelling cycles
Geophysical techniques such as refraction seismic tomography and electrical resistivity tomography are suitable field methods for gaining a better understanding of (invisible) subsurface processes from the surface of the landform (e.g.infiltration of meltwater or degradation of ground ice) (Hilbich et al., 2009;Merz et al., 2015).Their non-invasive nature also allows repeated measurements of the same profile (Hilbich, 2010;Hilbich et al., 2011;Kneisel et al., 2014), which can considerably improve the identification of hydrological processes in alpine permafrost systems (Langston et al., 2011;Wright et al., 2009), non-permafrost catchments (Kobierska et al., 2015) and in 1-D infiltration studies (Scherler et al., 2010).Rock glaciers, however, belong to the most challenging landforms for geophysical monitoring.Firstly, their rough surface characteristics consisting of very coarsegrained material and boulders, their small-scale topography and their high ice contents make good-quality measurements and data inversions difficult to achieve (Hilbich et al., 2009;Hauck and Vonder Mühll, 2003).Secondly, the movement of the landforms disturbs permanent installations and significantly alters the measurement geometry with time (Wilkinson et al., 2010;Loke et al., 2013).2009) demonstrated the general applicability of repeated ERT surveys for monitoring permafrost evolution on rock glaciers and highlighted the use of synthetic modelling to aid the sometimes ambiguous interpretation of the results.A systematic use of synthetic modelling in socalled forward-inverse cycles was presented by Fortier et al. (2008) and Kneisel et al. (2008) to ameliorate the process of data acquisition, inversion and interpretation for geophysical investigations in alpine or arctic settings.
Forward modelling of a synthetic model allows us to generate synthetic measurement data sets (in our case apparent electrical resistivity for ERT and seismic travel times for RST), which can then be inverted to investigate the response of the complete measurement-inversion procedure to a given subsurface structure and measurement geometry.By this, the potential resolution of structures and processes of interest as well as the unwanted occurrence of inversion artefacts can be analysed (e.g.Hilbich et al., 2009).
In this study, we use forward-inverse cycles to compare the inversion results of synthetic data sets to the inversion results of real field data as well as to analyse the potential resolution of three geophysical approaches regarding processes which are relevant for permafrost degradation: (1) ERT alone, (2) RST alone and (3) ERT and RST combined in the 4-phase model.While forward-inverse cycles for permafrost applications have already been introduced in previous studies (Fortier et al., 2008;Kneisel et al., 2008;Hilbich et al., 2009), a forward-inverse based appraisal of a combination of methods (as in the 4PM) has never been attempted before.It also allows the performance of the 4PM to be explicitly analysed with respect to a given synthetic model and measurement geometry.

Work flow
Figure 2 illustrates the main steps in the work flow applied in this study.The measured field data, apparent resistivities and travel times are inverted in Step 1 to obtain tomograms with the spatial distribution of specific resistivities and P-wave velocities.Inspired by the inverted field data, a synthetic reference model was set up, which was chosen to be (a) consistent with the measured data of Step 1, as well as (b) simplified regarding geometry and small-scale structure for use as a starting point for various scenarios.Realistic values for ρ s and v p were then assigned to all layers (see next section).These values were taken from the literature (see Table 1) and adapted in order to fit the inverted field data as best as possible.
In Step 2, synthetic data sets of ρ a and seismic travel times were generated from the reference model through forwardmodelling, applying the same measurement geometry as for the field data, and were subsequently inverted using the same inversion parameters as in Step 1. Through comparison of the inverted tomograms of synthetic reference model and field data and repetition of the subcycle consisting of forward, inverse and appraisal problem, the synthetic model can itera-tively be adapted until it produces tomograms of ρ s and v p , which correspond to those of the field data within acceptable bounds (step 3).
Three modifications of the baseline synthetic model were then developed (see next section) to represent typical mountain permafrost processes in the context of degradation.In Step 4, the ERT and RST tomograms of the reference model and the three modified scenarios served as input to the 4PM and the individual fractions of ice, water and air were calculated from (a) the geophysical values in the synthetic models, (b) the inverted synthetic models and (c) the inverted field data.Finally, the results of the 4PM for all three approaches were compared and analysed regarding their resolution capacity for the different processes represented in the three scenarios (step 5).

Synthetic models and data inversion
The synthetic models were created using the example of the Becs de Bosson rock glacier.Figure 3a shows the inverted tomograms of ERT (top) and RST (bottom) measurements from 25 July 2013, with black lines representing the individual interpretation of the tomograms.All field and synthetic data were inverted with the RES2DINV (Loke, 2015) and ReflexW (Sandmeier, 2015) software using a robust L1 Gauss-Newton inversion for the geoelectric data (and standard inversion parameters as described, e.g. in Hilbich et al., 2009), and the simultaneous iterative reconstruction technique for the seismic data (inversion parameters chosen according to Hilbich, 2010).Certain layers are visible in both tomograms, but the two methods also deliver complementary information.Layer 1 is clearly detected by both methods and is interpreted as the active layer consisting of coarse blocks.However, this layer is considerably thinner in the ERT (about 3-6 m) than in the RST tomogram (5-10 m).In the upslope part of the profile the active layer is characterized by lower resistivities and partly also lower velocities, which is interpreted as finer-grained material with higher moisture content (Layer 2).Layer 3 is an ice-rich layer with roughly 10 m thickness, which is very prominent in the ERT data (very high resistivities of > 100 k m), but also visible in the RST data (∼ 3500 m s −1 ).Towards the upslope part, both resistivity and velocity gradually decrease and indicate decreasing ice content.Zones 4 and 5 are therefore interpreted as icepoor and ice-free zones with similar thickness, zone 4 can, however, rarely be identified in the RST data.Whereas both methods give clear indications that the bedrock is present at approximately 15 m depth, only the seismic method is able to delineate the details of the bedrock topography.
A simplified structure of these layers and their corresponding geophysical properties is then translated into the baseline synthetic model (SYN1) with the aim to reproduce the observed subsurface composition for a shortened part of the profile.Note that an exact correspondence of the synthetic model to the field case of Becs de Bosson rock glacier is Table 1.Geophysical properties used in synthetic reference models.The values have been taken from literature (Barsch, 1996;Hauck and Kneisel, 2008;Hilbich et al., 2009 and references therein) and adapted to match the inverted field data in an optimal way.

Material ρ w v p
Unfrozen sediment (moist -dry) 3500-10 000 m 330-650 m s −1 Ice-rich permafrost 200 000-2 000 000 m 3500 m s −1 Ice-poor permafrost 200 000 m 3300 m s −1 Pore water (in sediment) 100-500 m 1500 m s −1 Bedrock 3000 m 6000 m s −1 not crucial for the appraisal of the 4PM approach.It merely serves as illustrating example.Figure 3b (top panel) shows the resulting model: an active layer, an intermediate layer with three lateral zones of different ice contents (ice-rich, icepoor and no ice), and the bedrock at the bottom of the profile.The baseline sensor spacing for all scenarios is 4.5 m, as this is the set-up at the Becs de Bosson rock glacier and it represents the commonly used sensor geometries in other high alpine permafrost investigations (between 4 and 5 m spacing; see Hauck and Kneisel, 2008).In a sensitivity experiment (Sect.5) this spacing is decreased in two steps down to 1125 m in order to investigate the influence of sensor spacing on the inversion results and finally on the 4PM.
Based on this baseline model three different scenarios are set up: -SYN2 is a saturated 0.5 m thick layer on top of the ice core.
-SYN3 is two 0.5 m wide vertical preferential flow paths, of which one crosses the ice core.All flow paths end up in a horizontal saturated layer at 10 m depth, beneath the ice core.
-SYN4 is a degradation of the ice core by 1 m thickness from above with a corresponding thickening of the active layer.
Whereas SYN4 represents the classical permafrost degradation scenario (see e.g.Hilbich et al., 2009), scenarios SYN2 and SYN3 depict two hydrological situations that were hypothesized for rock glaciers within the seasonal cycle of thawing and freezing (e.g.The geophysical properties in the synthetic models are assigned in the forward-modelling module of ReflexW for the seismic data (Sandmeier, 2015) and the data generator RES2DMOD for the ERT data (Loke, 2002).As the geoelectric forward-modelling software RES2DMOD does not allow topography, no topography was included in both synthetic models.Five percent of Gaussian noise was added to the synthetic measurement data to simulate real (i.e.non-ideal) field data.The resulting synthetic measurement data are then inverted using RES2DINV (Loke, 2015) and ReflexW (Sandmeier, 2015) software as described above for the field data.No further a priori information was prescribed in the inversion to avoid the introduction of (wanted or unwanted) additional biases into the 4PM evaluation.

4-phase model (4PM)
The 4PM is a petrophysical model that relates specific electrical resistivities and seismic P-wave velocities using two well-known petrophysical mixing rules (Archie, 1942;Timur, 1968) to obtain volumetric fractions of air, water, ice and rock (Hauck et al., 2011).As electrical resistivities are well suited to differentiating between ice and water, and seismic velocities are sensitive to transitions between air and ice or rock, the combination of these methods is in principle suitable for estimating the lateral heterogeneity and quantifying the volumetric contents of ice-rich bodies in a rock glacier.The model has already been successfully applied to alpine permafrost bodies (Hauck et al., 2011;Schneider et al., 2013;Pellet et al., 2016) and uses the following equations to determine the volumetric ice (f i ), water (f w ) and air content (f a ) for a given porosity model (x, z) ( = 1 − f r ; f r being the rock content): (1) where a(= 1 in many applications), m (cementation exponent) and n (saturation exponent) are empirically determined parameters (Archie, 1942), ρ w is the resistivity of the pore water, v r , v w , v a , v i are the theoretical P-wave velocities of the four components, and ρ(x, z) and v(x, z) are the inverted resistivity and P-wave velocity distributions, respectively.The free parameters in Eqs. ( 1)-( 3) are the so-called Archie parameters n, m and ρ w , the seismic velocities of the four materials and the porosity distribution (x, z), the latter being estimated as direct measurements rarely exist.
A comparatively simple porosity model was chosen with a 10 m thick 40 % layer representing the blocky material and the bedrock layer below (10 %; see top panel in Fig. 6).A bedrock porosity of 10 % is probably overestimated compared to the real conditions of Becs de Bosson rock glacier (given as gneiss, quartzite; see Sect. 2).However, as the underlying petrophysical equations of the 4PM are badly defined for very low porosities, 10 % was chosen for numerical reasons.Apart from this minimum constraint, sensitivity analyses of the 4PM dependence on porosity show that the absolute ice content values follow an almost linear relationship with available pore space.Consequently, the calculated saturation values do not change much for small variations in porosity.Further details of 4PM sensitivities are given in Hauck et al. (2011) and Pellet et al. (2016).Note also that we refrained from adding porosity heterogeneities in the coarse blocky layer so as not to bias the suggested ice content differences in the synthetic models of the baseline and scenarios.
Values for the seismic velocities v w , v a , v r and v i as well as n and m are taken from the literature (Schön, 2015;Knödel et al., 2013;Zimmerman and King, 1986) (see Table 2).The highly sensitive parameter ρ w (Hauck et al., 2011) was estimated using experience from various electrical conductivity measurements in alpine permafrost regions (unpublished).4PM plausibility check is achieved by verifying that physically consistent solutions are available for a specific data pair of ρ s and v p and the parameter specified (e.g.sum of f i , f a and f w is ≤ , and all volumetric contents are ≤ 0).The resulting ice, water and air content distributions can then be analysed as absolute values or saturation values relative to the available pore space, the latter being less dependent on uncertain porosity models.For the field and synthetic examples of this study, volumetric contents are displayed as changes to the baseline model SYN1 to illustrate the different processes of the scenarios (SYN2, SYN3, SYN4).
In this paper, we conducted three different steps to evaluate the performance of the 4PM: a.We run the 4PM using the pure synthetic models (without inversion) as a demonstration of the proof-ofconcept of the 4PM.
b.We use the inverted synthetic models as input for the 4PM to analyse the general performance of the 4PM.
c. From (a) and (b) we identify primary error sources of the model and attribute them to either the input data or the 4PM itself.Artefacts become apparent by comparing both input data types the resolution capacity of the different geophysical methods as well as the occurrence of inversion.The resolution capacity of temporal changes in the subsurface conditions is then examined by a time-lapse analysis, which compares 4PM results Table 2. Overview of the values attributed to the free parameters in the three equations of the 4PM in this study.
Parameter Value of the three scenarios to the baseline model.Finally, the procedure is repeated with reduced sensor geometries to identify the importance of measurement geometry on inversion and 4PM for the given scenarios.4c, d).An appearance of a water saturated layer on top of the ice core would, therefore, result in a decrease in resistivity (red horizontal bar in SYN2), but an increase in P-wave velocity (blue bar in SYN2).Ice loss (SYN4) can be detected by a contemporaneous decrease in resistivity and velocity, and the preferential flow pattern would be visible by resistivity decrease and velocity increase in the upper part, but resistivity and velocity decrease in the deeper part of the profile (SYN3).Velocity and resistivity values vary up to 50 % compared to the baseline model, with generally higher magnitudes for the RST data.All three scenarios show a unique pattern of changes in ρ s and v p and should, therefore, in general be distinguishable by the combination of ERT and RST.By construction, the defined zones in the synthetic models have sharp edges and show large contrasts to their neighbouring zones with different geophysical characteristics.

Inverted synthetic models
Figure 5a and b show the inverted models of the synthetic data presented in Fig. 4 as well as their relative changes with respect to the baseline model (Fig. 5c, d).Due to the standardly used regularization constraints of the inversion (see Hauck and Vonder Mühll, 2003;Hilbich et al., 2009)    models represent a smoothed version of their counterparts in Fig. 4, but the simulated processes are not always visible by looking at the ERT and RST tomograms alone.They first become visible in the analysis of the significant changes in ρ s and v p , seen in Fig. 5c and d.The spatial resolution of most features is weak due to the coarse electrode/geophone spacing of 4.5 m in comparison to the size of the anomalies.However, the unique pattern of resistivity and velocity changes identified in Fig. 4 (marked by the blue, green and pink rect-angles, respectively) is also clearly visible in the inverted tomograms.The thickness of the layers with water saturation (SYN2 and SYN3) and ice loss (SYN4) is strongly overestimated.Moreover, the depth at which ice loss is detected in SYN4 is too shallow, especially in the RST tomogram.The vertical pattern of the preferential flow paths is clearly visible in the velocity changes and to a lesser amount in the resistivity changes, but horizontal flow paths below the ice body are poorly resolved.From Fig. 5c and d it is also visible that in- version artefacts are produced with both methods for almost all scenarios.The thaw depth of 5m is delineated reasonably well in the ERT results (Fig. 5a) and slightly overestimated in the RST data (Fig. 5b)

Synthetic models
Figure 6 shows the calculated 4PM results for the three theoretical scenarios of Fig. 4a and b.The simulated water, ice and air content patterns are correctly reproduced by the 4PM, as should be expected, because the 4PM uses the resistivity and velocity values of Table 1, which were also used to transfer the synthetic models to their respective ERT and RST counterparts.However, to achieve this good match two specific constraints were necessary: a.Low ρ s and v p values at the surface often cause erroneous values in the 4PM.In mountain permafrost studies, these low values at the surface correspond to the thaw depth (ice content = 0), which can therefore be estimated by the so-called 3-phase model (3PM; see Pellet et al., 2016).In such cases, the ice content is automat-ically set to 0 % to better constrain the calculation of water and air contents.
b.A well-chosen porosity model including the bedrock topography is crucial for modelling realistic ice contents (see Sect. 3.4 and the uppermost panel in Fig. 6a).A combined analysis of ρ s and v p and the 3PM-derived porosity distribution (assuming an ice content = 0) allows the bedrock topography to be delineated: in cases of high resistivities and velocities 3500 m s −1 at greater depth a low porosity (corresponding to bedrock) can be assumed.
Following from the good match in Fig. 6, the 4PM-calculated changes in ice, water and air contents between synthetic reference model and scenarios also show the expected perfect match with respect to those shown in Fig. 4c, d (not shown).

Inverted synthetic models
When the 4PM is not applied to synthetic models but to the inverted data, i.e. after one complete forward-inverse cycle, the individual uncertainties with respect to measurement geometry and inversion (see Fig. 6) and with respect to the formulation of the 4PM become apparent (Fig. 7).Most simulated processes are reproduced but with weaker spatial resowww.the-cryosphere.net/11/2957/2017/The Cryosphere, 11, 2957-2974, 2017 lution and smaller amplitudes.The best results are obtained for SYN4, which unambiguously shows the ice loss, which is compensated by an increase in air content.However, the size of this anomaly is strongly exaggerated, and its depth is too shallow.In addition, the overall magnitude of changes is lower (up to 25 % in SYN2 and SYN4 as opposed to the correct value of 40 %).Although the depth of the saturated layer in the SYN2 scenario is correctly delineated, the strong increase in water content, caused also by an overestimated thickness of the saturated layer in the resistivity inversion (Fig. 6), is partly compensated by an erroneous decrease in ice content.In general this means that detected ice losses in the 4PM could also originate from inversion artefacts (mainly in the ERT data) caused by strong resistivity contrasts, as in the SYN2 scenario (saturated flow above the ice body).In addition, the 5 % Gaussian noise produced a similar anomaly pattern on the right-hand boundary of the tomogram, which is also a pure inversion artefact.
The vertical preferential flow paths (SYN3) can roughly be delineated and the resulting pattern is similar (but smoothed) to the synthetic models.The horizontal saturated layer below the ice body, however, produces a similar artefact to that for the SYN2 experiment, with a too thick and too strong decrease in ice content that compensates for a similar increase in water content.

Geophysical monitoring
The synthetic modelling results of the previous sections have indicated that the three permafrost and hydrology-related processes can be detected on rock glaciers with a combination of ERT and RST, although the size and location of the anomalies may be incorrectly delineated, especially when coarse sensor spacings are used (4.5 m in the present case).In this section, the approach is now tested on real field data from the Becs de Bosson rock glacier.Figure 8 shows the ERT and RST inversion results for three measurement dates at the beginning and end of summer 2012 and 2013.The consistency of all tomograms is clearly seen with a strong contrast between the active layer (lower resistivity and velocity values) and ice body (higher resistivity and velocity values).Lateral variability exists both in the ERT as well as RST data, corresponding to the lateral variability of the ice core (see detailed interpretation in Sect.3.3).Temporal changes are supposed to be caused solely by water infiltration and seasonal ice melt processes.Note that the thaw depth is significantly smaller in the ERT results than in the RST results (for July 2012 ∼ 2-3 m in the ERT and ∼ 4-5 m in the RST data).
Figure 10 shows the percentage resistivity and velocity changes for a 1-year period (beginning of summer measurements, upper panel) and over a summer season in 2013 (lower panel).It is seen that the percentage P-wave velocity changes are much larger than corresponding resistivity changes, highlighting again the general usefulness of RST monitoring in permafrost applications (Hilbich, 2010;   -7).Dräbing, 2016).Resistivity generally decreased between July 2012 and July 2013 with a predominant increase in Pwave velocities, with the largest changes in the downslope part near the rock glacier tongue (Fig. 9a).This change pattern corresponds to the SYN2 experiment in Figs. 4 and 5, indicative of an increase in water content in the active layer.Increasing water content in an unfrozen blocky layer would lead to a resistivity decrease and a P-wave velocity increase, as the P-wave velocity of water is higher than that of air (see Tables 1 and 2).However, the seismic results (Fig. 9b upper panel) also show regions with decreasing P-wave velocity (which partly correspond to zones with unchanging or slightly increasing resistivity), indicating that there are also regions with decreasing water content present, especially in the upslope part of the rock glacier.
The general increase in saturation between July 2012 and July 2013 can be explained by the temporal difference in complete snowmelt during these 2 years.According to ground surface temperature (GST) measurements  have led to wetter conditions in July 2013 at the time of the geophysical measurements compared to 2012.
For the seasonal evolution during the year 2013 (Fig. 9, lower panel), the generally decreasing resistivity and velocity patterns correspond to the SYN4 experiment (Figs. 4 and 5 right panel).These results visualize the thawing of the active layer between the two measurement dates at the end of July and end of September, which corresponds roughly to the snow-free period.Borehole temperatures from the nearby PERMOS station Lapires BH98/01 at a similar altitude indicate a propagation of the thawing front from about 1.0 to 4.9 m between 25 July and 27 September 2013 (PERMOS, 2016).The contemporaneous decrease of resistivity and velocity during that period is therefore indicative of seasonal ice loss in the active layer and corresponds to the process described in the synthetic model experiment SYN4, albeit at a shallower depth.Similarly to the conclusions drawn from the synthetic experiments, resistivity changes occur at shallower depths as velocity changes but can be assumed to be caused by the same process (i.e.active layer thawing).No vertical structures are observed in Fig. 9; hence there are no indications of vertical preferential flow through the ice core.

4-phase model
Figure 10 shows the 4PM results for the two periods shown in Fig. 9. Now, the processes, but also the uncertainties described above, can be analysed in more detail.The seasonal ice loss due to the advance of the thawing front during summer is clearly seen in Fig. 10b with a calculated volumetric ice content loss of 10-15 %.This ice content loss is compensated by an increase in air content (drying of the subsurface), especially in the central and upper parts of the profile.Water content changes are small, indicating a predominant run-off of the meltwater, except for a small increase on top of the presumed ice layer in the lower part of the profile.Lateral variability is high, which might be due to the large heterogeneity of the rock glacier or due to the small-scale noise in the respective measurements, which can be amplified during inversion (e.g.Hauck et al., 2003).As the resistivity contrast in the downslope part of the field data is much smaller than in the synthetic data (comparing Fig. 8 left panel and Fig. 5a), the inversion artefacts are assumed to be weaker or non-existent; i.e. the spatial dimensions of the changes in ice, water and air contents are assumed to be more realistic.Note that the overall changes are the same order of magnitude as the synthetic experiments within 10-15 %.
The interannual changes between 2012 and 2013 (Fig. 10a) produce a more complex pattern in the 4PM results.Ice content changes in the uppermost ∼ 5 m (corresponding to the active layer) are predominantly negative and positive below.Compensation is calculated mainly through increased (upper part) and decreased (lower part) air contents, but also through a slight water content increase at the top of the ice body.The latter could be indicative of a change in saturation conditions in the active layer as hypothesized in the previous section.
It has to be noted though that the compensation between ice and air with a residual change in water content may have two causes, one being process-based and the other purely nu-merical: (a) due to the coarse-grained material and the slope of the rock glacier, seasonally or interannually melting ice would result in quick infiltration and run-off on top of the ice layer, which causes a short-term decrease of resistivity due to the presence of meltwater.After the water left the system, resistivities may increase due to the increased air voids or return to their former values if consolidation and settlement of the surface material takes place.(b) The 4PM model formulation and the strong dependence of resistivity on water content (as opposed to ice, rock and air content, which are all considered as electrical isolators) lead to a tendency of compensating effects between air and ice contents with independent determination of water content (see the structure of Eqs. ( 1)-( 3) and a thorough sensitivity analysis shown in Hauck et al., 2011).

Discussion
The geophysical monitoring and 4PM study using synthetic data presented above may lead the way to an operational geophysical monitoring of ice degradation and hydrologic processes within rock glaciers (or other coarse-grained permafrost material), although several limitations have to be faced.In order to distinguish between limitations caused by the inversion process, those caused by the 4PM model formulation and those caused by the geometry of the geophysical surveys, the results of inverted synthetic data and not inverted ("true") synthetic data are discussed separately.In addition, the simulations were repeated with smaller sensor spacings (see below and Supplement).
4PM modelling of the true synthetic data show almost perfect results for the three different processes (Fig. 6) However, it has to be noted that a realistic porosity model is necessary to obtain realistic ice contents from the 4PM (see Hauck et al., 2011;Pellet et al., 2016).If no additional data are available (e.g.boreholes) it is crucial to approximate the bedrock topography of the rock glacier using the geophysical data or a 3-phase model simulation (3PM).Based on the assumption that no ice is present, the 3PM determines the porosity distribution, and bedrock can be identified from low-porosity zones with intermediate resistivities (e.g. a few 1000 m) and high P-wave velocities ( 3500 m s −1 ) (Pellet et al., 2016).From this information, the bedrock topography can be delineated and included in the porosity model of the 4PM (Fig. 6).It has to be repeated, however, that the models presented in Fig. 6 describe the expectable ideal case, in which some often unknown input parameters for the 4PM were prescribed (e.g. the P-wave velocity of the rock material and the pore water resistivity) by using the same values as for the construction of the synthetic models.This is not always the case in real applications.
Vertical structures of increasing water content with respect to the baseline model are resolved relatively well.The horizontal layer beneath the ice core in SYN3 is difficult to dis-tinguish as 4PM-calculated water content changes are very small.This was expected, as changes in water content are predominantly seen in the ERT results, but an ice core acts as an electrical isolator, strongly decreasing the sensitivity of the method underneath.In addition, detection of a shallow water layer beneath an ice core also poses problems to RST, as seismic velocities would decrease with depth in this case and are therefore difficult to detect by refraction seismics (see Hilbich, 2010).
Largest errors are introduced by horizontally saturated layers (such as in SYN2), which are overestimated in thickness due to inversion artefacts.This leads to an erroneous decrease in ice content and in false ice loss detection in this scenario.In contrast, the true ice loss scenario in SYN4 is correctly reproduced, with a slight underestimation of the decreased ice content and an overestimation of its vertical extent.It also predicts correctly the resulting increase in air content in the part of the tomogram concerned.In spite of the limitations described above, most processes can still be detected in the 4PM results, although interpretation remains difficult.
The spatial resolution of the 4PM calculations for the inverted synthetic data is much lower than for the true noninverted synthetic data, and, therefore, the magnitude of changes is underestimated in the former experiments (see Fig. 7).This is, on the one hand, caused by the smoothing effect of the inversion process.In the absence of a priori information, such as borehole data, standard inversion algorithms like the robust L1 Gauss-Newton used in this study have to apply smoothing constraints as regularization due to the underdetermined nature of the inversion problem (e.g.Loke and Barker, 1995).On the other hand, the low resolution and underestimated magnitude of change are due to the coarse sensor geometry, which was intentionally chosen to be identical with the field data to simulate real field conditions.The sensor spacing of 4.5 m for the monitoring at Becs de Bosson is comparable to many other field sites in permafrost geophysics, where the profile length (ideally covering the entire landform) and the investigation depth (reaching the base of the landform) have similarly high priority to spatial resolution.
To analyse the specific effects of acquisition geometry, the sensor spacing of the synthetic set-up was decreased in two steps from 4.5 m down to 1.25 m.All results of this sensitivity analysis (corresponding to Figs. 5 and 7) are shown in the Supplement.In the following, we illustrate the sensitivity using the case of SYN4 (Figs. 11 and 12).Surprisingly, the overall impact of the increase in resolution is low, although the decrease in sensor spacing leads to moderately improved RST tomograms (Fig. 11b) and slightly fewer inversion artefacts in the ERT tomogram (Fig. 11a).However, neither the decrease of spacing from 4.5 to 2.25 m nor the further decrease to 1.125 m had any significant impact on the 4PM results, showing that a smaller geophone/electrode spacing does not automatically lead to a higher resolution capacity (Fig. 12).This is also due to the spatial dimension of the  anomalies in the synthetic scenarios as well as the fixed inversion parameter settings.In a real field study, these settings (regularization parameter, mesh size etc.) would have to be changed according to the different sensor geometries.A thorough analysis of all sensitivities regarding data acquisition geometry, ERT and RST inversion parameter and 4PM parameter is, however, beyond the scope of the present paper.It should, however, be emphasized that the resolution capacity is expected to increase significantly with smaller spacings in case of shallower anomalies, e.g. for field cases with smaller active layer depths or for landforms with less ice content (resulting in smaller resistivity contrasts).
Assuming that the small magnitudes of change observed in our field data (Fig. 9) correspond to the poorly resolved changes in our synthetic experiments which originate from much stronger magnitudes, our results strongly suggest that observed processes are expected to be (i) much more pronounced in reality but with (ii) a smaller spatial dimension than illustrated by the geophysical inversions and the 4PM.Further, a systematic difference in the resolution of the thaw depth was observed, with a significantly smaller thickness in the ERT data compared to the RST data (see Figs. 5 and 8).The different inversion characteristics (smoothing and sensitivity, e.g.Marescot et al., 2003;Maurer and Hauck, 2007;Hilbich et al., 2009;Hilbich, 2010) of ERT and RST lead to an overall overestimation of the spatial extent of both, conductive and resistive anomalies in the inversion results.
In the current version of 4PM and for a completely prescribed porosity model, water content is only determined through Eq. ( 1) and not corrected by the RST data.Pellet et al. (2016) have shown that this can be improved through coupled 3PM and 4PM simulations, as well as near-surface calibration of the 4PM by measured soil moisture values.Future work should also explore the possibility of using alternative formulations of the resistivity mixing rule that takes into account the resistivity of the rock material (Python, 2015).Finally, using alternative ERT inversion algorithms may also improve the delineation of sharply limited horizontal features.
The analysis of the results from the inverted synthetic data illustrates the difficulties for the interpretation of 4PM calculations of real field data.Inversion artefacts which are transferred to the 4PM may lead to false interpretations of the subsurface composition and the underlying processes (like degradation of ice cores or the existence of saturated layers).In combination with higher measurement errors in permafrost environments (up to 20 % instead of the 5 % used in the synthetic experiments; see e.g.Rosset et al., 2013) and only little a priori information about the real subsurface structures, an interpretation of 4PM results can be much more complex and difficult.
However, the results of this study also demonstrate that an interpretation of geophysical tomograms can be significantly facilitated through 4PM-based geophysical monitoring compared to an interpretation of single tomograms be-cause small-scale anomalies and processes are hardly resolvable from visual inspection of a tomogram alone (see Fig. 5a  and b).

Conclusion
We presented a detailed study about the resolution capacity of standard geophysical monitoring methods for the detection and analysis of ice and water content changes in the context of mountain permafrost evolution.The study is based on synthetic and field data and analyses the individual performance of electrical resistivity tomography (ERT) and refraction seismic tomography (RST) inversions as well as their combined analysis within the so-called 4-phase model (4PM).Synthetic modelling and field results are shown for ground ice and water saturation changes illustrated with a real field case from a rock glacier in the western Swiss Alps.In the synthetic experiments, three scenarios were analysed: (i) the development of a saturated layer on top of the ice core, (ii) preferential flow within the active layer and permafrost and (iii) ground ice loss, i.e. permafrost degradation.Key results of this study include the following: -Our approach, including complete forward-inverse cycles of synthetic data for the 4PM, proved to be useful for estimating its resolution capacity regarding typical permafrost evolution scenarios.Limitations arising from the individual inversion and the sensor spacing as well as from the coupled approach with the 4PM can be distinguished.
-In general, monitoring adds information over singular surveys, as not only can physical processes be assessed, but also the negative impact of singular, and geometryinduced inversion artefacts tend to be diminished.This is especially the case when using the combined analysis within the 4PM.
-The 4PM is able to reproduce the true (that is, noninverted) synthetic model scenarios almost perfectlyinaccuracies in the final 4PM results of the inverted scenarios stem rather from the individual inversion themselves.
-Inversion artefacts from the individual surveys are transported into 4PM analysis and lead to potential misinterpretation, such as wrongly detected ice loss in the synthetic case of the development of a saturated water layer on top of the ice core.
-For all geophysical monitoring with the 4PM a suitable porosity model is of great importance.Because porosity stays constant with time, small errors in the porosity model will have only small impacts in a 4PM monitoring context.-Due to the smoothing effect of most inversion algorithms, magnitudes of changes are often underestimated and vertical extents are often overestimated within the 4PM.Using different inversion algorithms with less smoothing should improve the magnitudes of change, but can be more susceptible for inversion artefacts.
-When applying the findings of the synthetic model experiments to field data from Becs de Bosson rock glacier, heterogeneous saturation changes were found in the thaw layer between July 2012 and July 2013, and the seasonal propagation of the thaw layer (seasonal ice melt) between July and September 2013 could be well delineated.
Competing interests.The authors declare that they have no conflict of interest.
Special issue statement.This article is part of the special issue "The evolution of permafrost in mountain regions".It is not associated with a conference.

Figure 1 .
Figure 1.Map and photograph of the Becs de Bosson rock glacier.The red line on map and photograph marks the position of the full geophysical profile, while the green line on the photo shows the joint ERT/RST profile, which is used for the 4PM analysis of this study.Map source: Swiss Federal Office of Topography swisstopo.

B.
Mewes et al.: Resolution capacity of combined geophysical monitoring Hilbich et al. (

Figure 2 .
Figure 2. Overview of the work flow applied in this study.Here, ρ a denotes the apparent resistivity, ρ s the specific resistivity and v p the P-wave velocity.

Figure 3 .
Figure 3. Creation of synthetic models based on inverted field data of the Becs de Bosson rock glacier: (a) inverted ERT and RST tomograms of 25 July 2013 with superposition of interpreted zones, (b) simplified synthetic baseline model (SYN1) based on the interpretation in (a) and different scenarios of hydrological processes (SYN2 to SYN4).

Figure
Figure4aand b show the ERT and RST synthetic models for the baseline and the scenario experiments.By subtracting the baseline model SYN1 from all scenarios the resistivity and velocity changes for the various scenarios are determined (Fig.4c, d).An appearance of a water saturated layer on top of the ice core would, therefore, result in a decrease in resistivity (red horizontal bar in SYN2), but an increase in P-wave velocity (blue bar in SYN2).Ice loss (SYN4) can be detected by a contemporaneous decrease in resistivity and velocity, and the preferential flow pattern would be visible by resistivity decrease and velocity increase in the upper part, but resistivity and velocity decrease in the deeper part of the profile (SYN3).Velocity and resistivity values vary up to 50 % compared to the baseline model, with generally higher magnitudes for the RST data.All three scenarios show a unique pattern of changes in ρ s and v p and should, therefore, in general be distinguishable by the combination of ERT and RST.By construction, the defined zones in the synthetic models have sharp edges and show large contrasts to their neighbouring zones with different geophysical characteristics.

Figure 4 .
Figure 4. Synthetic model scenarios for (a) ERT and (b) RST, and percentage changes between (c) ρ s and (d) v p of the three scenarios with respect to the baseline model.

Figure 5 .
Figure 5. Inverted synthetic model sections of the synthetic model scenarios of Fig. 4 for (a) ERT and (b) RST, and percentage changes between (c) ρ s and (d) v p of the three scenarios with respect to the baseline model.The positions of the introduced changes between the different scenarios are indicated by black lines.

Figure 6 .
Figure 6.4PM results for the synthetic model scenarios (see Fig. 4a, b).The blue, green and pink rectangles correspond to the same rectangles in Figs. 4, 5 and 7.The prescribed porosity model is shown in the uppermost panel.

Figure 7 .
Figure 7. 4PM calculated changes in volumetric fractions for the three scenarios relative to the baseline model for the inverted synthetic data.The positions of the introduced changes between the different scenarios are indicated by black lines.

Figure 8 .
Figure 8. Inverted geophysical monitoring data from the Becs de Bosson rock glacier for three measurement dates in 2012 and 2013.

Figure 9 .
Figure 9. Temporal changes in (a) resistivity and (b) P-wave velocity between the measurements shown in Fig. 8. Upper panel: interannual changes between 21 July 2012 and 25 July 2013.Lower panel: seasonal changes between 25 July and 27 September 2013.The coloured rectangles correspond to the resistivity and velocity change patterns in the synthetic experiments (see Figs. 5-7).
(PER-MOS, 2016), snowmelt was completed about 8-15 days later in 2013 than in 2012 (16 June 2013 vs. 8 June 2012 at a snow-poor location close to the lower end of the profile; 24 July 2013 vs. 8 July 2012 at a snow-rich location close to the upper part of the profile).The later snowmelt in 2013 may www.the-cryosphere.net/11/2957/2017/The Cryosphere, 11, 2957-2974, 2017

Figure 10 .
Figure 10.Changes in volumetric fractions for the two time intervals of inverted monitoring data of the Becs de Bosson rock glacier: (a) 1 year (July 2012 to July 2013) and (b) 2 months (July 2013 to September 2013).Note the reduced range of the colour scale.

Figure 11 .
Figure 11.Comparison of (a) resistivity and (b) velocity changes of inverted tomograms for scenario SYN4 for three different sensor spacings.Top: 4.5 m, middle: 2.25 m, bottom: 1.125 m.The position of the introduced change is indicated by black lines.

Figure 12 .
Figure 12.Comparison of 4PM-calculated changes in volumetric fractions for the inverted synthetic data of scenario SYN4 for three different sensor spacings: (a) 4.5 m, (b) 2.25 m, and (c) 1.125 m.The position of the introduced change is indicated by black lines.
Mewes et al.: Resolution capacity of combined geophysical monitoring -In general, vertical structures can be resolved well, whereas the vertical dimensions of horizontal saturated areas are overestimated, especially by ERT.Ice loss is easily detectable with the 4PM, showing that the model is a valuable tool for geophysical monitoring of rock glaciers in the context of permafrost degradation.