A new model for quantifying subsurface ice content based on geophysical data sets

Introduction Conclusions References


Introduction
Permafrost underlies most polar and many mountainous regions of the Earth.Its ground thermal regime is primarily determined by the regional atmospheric and local (micro-)climatic conditions, site-specific surface characteristics, the topography (especially in high-mountain areas) and the heat flux from the earth's interior.In the context of climate change there is increased concern that temperature changes in the subsurface will be caused via changes in the ground-surface temperature, which in turn is determined by the energy balance at the surface.In many permafrost regions a climate induced warming of the subsurface and a corresponding thickening of the seasonally thawing surface layer (active layer) has already been confirmed by measurements (e.g., Frauenfeld et al., 2004;Marchenko et al., 2007;Brown and Romanovsky, 2008).In mountain regions with increasing subsurface temperatures, strong active layer thickening has been observed in specific years (Hilbich et al., 2008;Harris et al., 2009), which may lead to an increased occurrence of slope instabilities, when ground ice close to the melting point exists.
Based on this observational evidence, it is now recognized that a detailed knowledge of the internal structure of permafrost is required for modelling of the future evolution of the ground thermal regime in permafrost areas and geotechnical assessment of the hazard potential due to degrading permafrost.However, due to the commonly difficult access for permafrost locations in mountain areas and the corresponding logistical and financial difficulties in obtaining high-quality data sets, reliable estimates of the subsurface material compositions are still very scarce.
Partly or permanently frozen subsurface material can be composed of four different phases: two solid phases (rock/soil matrix and ice), liquid (unfrozen pore water) and gaseous phases (air-filled pore space and cavities).The physical properties, such as thermal conductivity or heat capacity, of the different phases are markedly different.
The physical properties of permafrost material are therefore determined by the relative contributions of the different phases to the bulk material.For a reliable thermal modelling in seasonal and permanently frozen ground the exact composition of the various phases must be known.However, except for laboratory analysis of field probes obtained near the surface or within boreholes (e.g., Arenson and Springman, 2005), the composition of the subsurface material can only be inferred through indirect geophysical investigations (for a review see Hauck and Kneisel, 2008 and references herein).
Due to the complexity of the subsurface, a combination of geophysical methods (e.g.electrical resistivity tomography and refraction seismic tomography) is favoured in most geophysical studies to minimise ambiguities in the interpretation of the results.A possible option to reduce the interpretational ambiguities is to impose relationships between the unknown model parameters (electrical resistivity, permittivity, seismic velocity, density) and the actual subsurface materials (rock, water, air, ice).Even though several model approaches exist to determine the bulk electric and elastic properties for 3-phase mediums there are no petrophysical relationships available for full 4-phase systems.Simple relations have been developed for the electrical properties, e.g. the well-known relation between the electrical resistivity of the probe, the pore-water resistivity, the porosity and the saturation known as Archie's Law (Archie, 1942) and for the elastic properties such as the time-average equation for seismic P-wave velocities by Wyllie et al. (1958) and its extension to the frozen phase by Timur (1968).These relationships were originally only validated for a restricted range of materials (e.g.unconsolidated sediments, Zimmerman and King, 1986).In later studies theoretical concepts for simple pore geometries were developed including both electric and elastic properties of the material (e.g., Sheng, 1990), thereby justifying the empirical relations mentioned above for a wider parameter range.
In permafrost studies, quantitative combinations of electric and seismic data sets were introduced by McGinnis et al. (1973), who used the resistivity information for calculating the increase in seismic P-wave velocity due to the frozen layer.Oberholzer et al. (2003) used an index, based on the ratio of resistivity and seismic P-wave velocity, to differentiate between frozen and unfrozen morainic material in the Swiss Alps.Hauck and Wagner (2003) applied fuzzy logic to combine resistivity and seismic velocity data sets for identifying regions with ground ice occurrences.Hausmann et al. (2007) used gravimetry, seismics and GPR to quantify the ice content in an alpine rock glacier.However, to our knowledge no physically-based relation between electric and elastic parameters commonly measured in surface geophysical surveys and the four phases present in permafrost material exists for practical applications in permafrost studies up to now.
In this study we will present a first approach for a so-called 4-phase model (4PM), which is amenable for explicit calculation of the four phases based on 2-D tomographic electrical and seismic measurements.For these two methods, numerous data sets from mountain permafrost sites in Europe are available.Two different approaches with increasing complexity will be presented.In a first step a porosity dependent model will be introduced, which estimates volumetric ice-, water-and air content relative to the available pore space.For this model, unique solutions exist for each pair of resistivity and seismic P-wave velocity data obtained by geophysical measurements.In a second step the 4PM is used to compute all possible solutions of the 4 phases without prescribing porosity and based on the same geophysical data sets.Here, in general no unique solution exists, but the general approach can be used to identify those model regions, where physically plausible solutions exist and where all 4 phases are constrained by the data (e.g. because one of the volumetric contents approaches zero).The 4-phase model is applied to two rock glaciers in the Eastern Swiss Alps, where corresponding borehole information for validation is available.

Theory and 4-phase model approach
In our 4-phase model (4PM) we assume that each cell of the 2-D model domain consists of the sum of the volumetric fractions for rock f r , water f w , ice f i and air f a .For each model cell the equation: must be fulfilled.In order to determine the specific fractions the model is further based on the two geophysical relationships mentioned above (Archie's law and an extended Timur equation) which will be introduced in the following.Under the condition that the clay fraction within the subsurface material is negligible, Archie's second law relates the resistivity ρ (in Ωm) of a 3-phase medium (rock matrix, liquid, pore space filled with air) to the resistivity of the pore water ρ w , the porosity Φ and the saturation, that is the fraction of the pore space occupied by liquid water S w : where a, m and n are empirically determined parameters (Archie, 1942) and the saturation can be expressed in terms of the volumetric fractions introduced in Eq. (1) using: In our model we assume that Eq. (2) still holds for partly frozen material of medium below-zero temperatures, because in the presence of electrolytic conduction in the pore water, ice can be seen as an insulator similar to the air and the rock matrix.In the mountain permafrost sites of the European Alps subsurface temperatures seldom fall below −5 • C, where unfrozen water can still be present.
The presence of ice in the pore spaces can cause large increases in seismic velocity compared to the velocity when the interstitial water is unfrozen, especially in highporosity sediments and debris (Hilbich, 2010).Since ice is much stiffer than water, the wave velocity is tightly coupled to the ice-to-water ratio.Rock is much stiffer than either ice or water; therefore the wave velocity is also a decreasing function of the porosity Φ (Zimmerman and King, 1986).Finally, due to the low velocity of P-waves in air, seismic velocity of the bulk material depends also on the saturation, that is, on the air-filled pore space.To include this in our model we chose a time-averaged approach, i.e. an extension of Timur's equation to 4 phases, which states that the reciprocal of the P-wave velocity (the so-called slowness) of a mixture, 1/v, is equal to the sum of the slownesses of the respective components, each weighted by its volumetric fraction: Similar to Eq. (2), Eq. ( 5) can only be seen as an approximation to the real conditions in the subsurface.Several other approaches exist to model the wave velocity of a 3phase medium (rock, ice, water), but all models assume that the composite density is the volume-weighted average of the densities of the constituents (for an overview and comparison see Carcione and Seriani, 1998).We assume that Eq. ( 5) is a fair approximation for the heterogeneous conditions in Alpine permafrost terrain -it certainly does not hold for sediments with a considerable amount of clay particles.Then, the dependence of the P-wave velocity on water content becomes non-trivial (Santamarina et al., 2005;Fratta et al., 2005).Equations ( 1), ( 2) and ( 5) form a set of equations for the four unknown volumetric fractions in dependence of the material constants ρ w , v r , v w , v a , v i as well as the (material dependent) free constants in Archie's Law a, m and n.The 2-D fields of the electrical resistivity ρ and the seismic P-wave velocity v are determined by tomographic measurements in the field.Combining Eqs. ( 1)-( 5) and solving for the ice content f i we obtain: Similarly, equations for the air content f a and the water content f w can be derived: In addition to the dependence on the material constants and the obtained data sets, Eqs. ( 6)-( 8) depend also on the rock content f r (=1−Φ).Due to the three governing equations and the four unknown variables, there is in general no unique solution to the system of equations.For special cases, such as saturated ground conditions (f a =0) or unfrozen conditions (f i =0) the number of unknowns can be reduced.For the general case of four phases, two different approaches are investigated within this study: 2. a general approach, where for a given data pair of electrical resistivity and P-wave velocity all possible 4PM solutions are calculated to determine the range of possible ice/water contents and to find the inherent model ambiguities.

Porosity dependent model
A first approach to solve Eqs. ( 6)-( 8) is to prescribe one of the four unknowns, e.g. the porosity (in addition to the constants in Archie's Law) homogeneously over the model domain and calculate the remaining volumetric fractions relative to the available pore space (Hauck et al., 2008).This approach is justified, if the porosity is known from a borehole log and can be assumed to be distributed homogeneously within the model domain.The material constants can be taken from literature or can be estimated in the laboratory using field samples (e.g., Sch ön, 2004).
Figures 1 and 2 show the three volumetric fractions calculated from the 4PM as a function of electrical resistivity and seismic P-wave velocity for a set of material properties corresponding to the debris-ice mixture often found at mountain permafrost sites (Table 1) and for porosities of 0.5 (Fig. 1) and 0.05 (Fig. 2).The values in Table 1 are taken from a series of seismic and resistivity field experiments during the EU-PACE project (Permafrost and Climate in Europe) and from literature (cf.Hauck and Kneisel, 2008;Kneisel et al., 2008).Values in Figs. 1 and 2 are only given for all resistivity and P-wave velocity data pairs which lead to a physically consistent solution of the set of Eqs. ( 6)-( 8).
For a porosity of 0.5 (characteristic for e.g.rock glaciers) this leads to a solution space between around 600-4500 m/s and above 1000 Ωm, depending on the prescribed values given in Table 1.A quick evaluation of the model performance can be performed by analysing the extreme cases of high air, high water and high ice contents.High air contents within the pore space are only found for very low velocities and very high resistivities, due to the electrically isolating characteristics of air and due to the seismic velocity in air of 330 m/s (Fig. 1a).For the calculation of the water content, an inherent problem of using Archie's Law becomes apparent in Fig. 1c water content does not depend on the seismic velocity even though its dependence is prescribed through Eq. ( 5).The reason can be found in Eq. ( 2), which relates the electrical resistivity only to water content and porosity (ice and air being electrically isolating).In prescribing the porosity, the water content depends only on the observed resistivity and the material properties.However, as this strong dependence of the electrical resistivity on the pore water content is one of the reasons for the good applicability of electrical methods in permafrost studies, this simplification is considered to be justified.In contrast to the air and water content, high ice contents can be found within a larger region of the functional space (Fig. 1b).Pure ice exhibits very high resistivities up to several MΩm and seismic velocities in a range between 3500 and 4000 m/s (e.g., Hauck and Kneisel, 2008).
For a porosity of 0.05 (Fig. 2), as for bedrock with low weathering grade, the solution space changes strongly to higher velocities and resistivities.Now, physically plausible solutions exist only for P-wave velocities between 3200-5800 m/s and resistivities above a few 10 kΩm, reflecting the typical ranges of bedrock material (see e.g., Sch ön, 2004;Hauck and Kneisel, 2008).The qualitative dependencies of air, ice and water content on high/low velocities and resistivities are the same as in Fig. 1, with high air contents for comparatively low velocities and high water contents for comparatively low resistivities.

General model
Theoretically, an infinite number of porosity values will lead to a solution of the equation system (Eqs.6-8), if suitable values for the Archie parameters m, n and a are used.In practice, porosity values leading to a solution of Eqs. ( 6)-( 8) are often constrained to a narrow range of possible values, as negative values for the air and ice content would occur especially for very high (or low) porosity values in the presence of high (or low) measured P-wave velocities (Eqs.6 and 7).In addition, due to the similar P-wave velocities of ice and rock and the fact that they both act as electrical insulators in Eq. ( 2), ice and rock will be more difficult to discern than water and air, the latter 795 having rather narrow and distinct ranges of possible velocity/resistivity data pairs (cf. Figs. 1 and 2).
To analyse this, the full solution space of Eqs. ( 6)-( 8) was explored using the Matlab symbolic math toolbox.This solution space can then be used to extract all theoretically possible solutions for the four phases for a given resistivity/velocity data pair of the measured data sets.Evaluations for several examples show that water and air contents are usually well constrained due to their marked differences in resistivity and velocity values.On the other hand, there is strong ambiguity between ice and rock content: the lower the ice content the higher the rock content.For these cases the sum of the ice and rock contents remains nearly constant meaning that the model can only determine the sum of ice and rock volumes, but not the individual contributions.
To evaluate this for the 2-D tomograms of the geophysical surveys, all possible solutions are calculated for each data pair.We then determine minimum and maximum values of each of the 4 phases and for all model blocks.Next, regions can be identified where a narrow range of solutions exists for some or all of the phases indicating reliable model results.If minimum and maximum values of e.g.ice content are close to 0 and 100%, respectively, no conclusion can be drawn regarding the ice content without further data (such as a porosity estimate).On the other hand this general model approach allows the determination of all 4 phases (including the rock content, i.e. the porosity) for cases where minimum and maximum values of all possible solutions for the 4 phases are similar.By this, regions with higher porosity can be differentiated from regions with lower porosity.

Site description and data sets
Both approaches were applied to geophysical data sets from two well studied rock glaciers in the Upper Engadine, Swiss Alps, namely the rock glaciers Muragl and Murt èl described in more detail in Maurer and Hauck (2007).At both field sites borehole data and additional geophysical data are available for model validation (Musil et al., 2002;  Maurer et al., 2003;Arenson and Springman, 2005;Hilbich et al., 2009).The tomographic geoelectric and refraction seismic data sets for the case studies presented here were obtained during several field campaigns and include an almost 300 m cross profile near the tongue of rock glacier Muragl and a shorter 145 m longitudinal profile at rock glacier Murt èl (Fig. 3, see also Maurer and Hauck, 2007).
Seismic velocities were obtained using a spread of 120 geophones and 53 shot points at rock glacier Muragl and a spread of 120 geophones and 44 shot points at rock glacier Murt èl with charges of 200-400 g explosives as source.Electrical resistivities were obtained using multi-electrode instruments with 36 (Muragl) and 30 (Murt èl) electrodes.As we focus on the development of the 4-phase model in this paper, we will not repeat specifics regarding measurement accuracy, inversion parameters and validation of resistivity and P-wave velocity results.We refer the reader to Maurer and Hauck (2007) for the details on data acquisition and inversion of both seismic and electrical data sets.
At rock glacier Muragl two 70 m deep boreholes (BH1/99 and BH2/99) are present along the survey line.Two further boreholes are located around 50 m upslope (Fig. 3a).
At rock glacier Murt èl two boreholes are present some 50 meters upslope of the survey line (Fig. 3b).Ground temperature data are available from all boreholes, laboratory analysis of core samples regarding ice-, air-and rock content are only available from boreholes BH 4/99 at Muragl and BH 1/00 at Murt èl (Arenson and Springman, 2005).

Results
Results for both rock glaciers and both model approaches are shown using the inversion results described in Maurer and Hauck (2007).We will discuss the results for rock glacier Muragl in more detail than Murt èl, because (1) the differences between porosity dependent and general model results are similar for both rock glaciers and ( 2) clearly shows the advantages and limitations of the 4PM approach for applications in mountain permafrost terrain.

Rock Glacier Muragl
Figure 4 shows the results of the electrical resistivity and seismic velocity inversion models for rock glacier Muragl.Even though some areas with anomalously high or low resistivity and velocity values are readily identified (e.g. the low velocity region near borehole BH1/99 (B1) or the high resistivity region between 200 and 250 horizontal distance), the overall distribution is non-trivial and resistivity and velocity models show markedly different patterns.A delineation of regions with enhanced ice-, water-and air content can not be determined directly from the two inversion models, but requires additional information such as borehole data or additional geophysical data sets (cf.Maurer and Hauck, 2007).

Porosity dependent model
Figure 5 shows the calculated ice-, water-and air contents for a constant porosity model of 50% over the same model domain as in Fig. 4. To obtain the distributions shown in Fig. 5, Eqs. ( 6)-( 8) were computed with the material constants shown in Table 1.The homogeneous porosity model of 50% was chosen to simulate the large air-, water-or ice-filled voids between the boulders of the rock glacier.Arenson and Springman (2005) found volumetric solid contents between 40-60% within the uppermost 15 m of borehole BH 4/99 making 50% a reasonable assumption as a mean value.Clearly, a constant porosity model likely oversimplifies the highly heterogeneous subsurface conditions within a rock glacier.The white areas in Fig. 5 delineate those model blocks, where no physically plausible solution of Eqs. ( 6)-(8) could be found, i.e. volumetric fractions were negative or did not add up to 1, violating the necessary conditions of Eq. ( 1).These regions are predominantly found below 10 m depth and at the boundaries of the model, suggesting either bedrock occurrences, violating the assumption of 50% porosity, or inversion artefacts due to the diminished sensitivity near the model boundaries.Comparing these regions to the seismic velocity results in Fig. 4a, it is seen that they are often coinciding with anomalously high P-wave velocities.As was shown in Figs. 1 and 2, the solution space for a specific porosity is bounded for high velocities, indicating the higher the porosity the lower the maximum velocity still giving a physically plausible solution.
Consequently, the white areas in Fig. 5 indicate the presence of rock occurrences with considerable smaller porosities than the assumed 50% at larger depths.
Even though the constant porosity imposes some interpretational constraints, the results from this porosity dependent model may be interpreted as indication to what extent the available pore space is filled with the respective phases, independent of the absolute value and the distribution of the porosity within the rock glacier.This is emphasised in Fig. 6 where the results of Fig. 5 are shown in relation to the porosity (e.g. % ice content per available pore space).
As can be seen from Figs. 5c and 6c the calculated air content within rock glacier Muragl is small except for the uppermost 5-10 m, where absolute air contents up to 20% are found.Data from borehole BH 4/99 indicate volumetric air contents between 3-10% between 5 and 15 m depth (Arenson and Springman, 2005) which is similar to the results obtained with the model.As high air contents in boreholes are difficult to measure due to the instability of the core caused by the large voids, the modelled air content values may even be more realistic than the data from the borehole.Maximum ice content (Fig. 5a) in the 4PM is found to the right of borehole BH 2/99 (B2) with values up to 45%, which is close to full ice saturation for a model porosity of 50% (Fig. 6a).Further ice occurrences are predicted between the two boreholes, but are limited to the uppermost 10-20 m.Again, this agrees well with data from the borehole BH 4/99, where ice contents between 50-70% were found between 5 and 10 m depth (with a corresponding rock content of around 30-40%, which is less than the prescribed value of 50%) and 40% at 15 m depth (where a rock content of 50% was found).The modelled water content (Figs.5b and 6b) is close to zero near the predicted ice core and greater than 30% below the uppermost 15 m and outside of the ice core.

Verification and interpretation
Borehole temperatures (Fig. 7) presented in Maurer et al. (2003) show values above the freezing point throughout the entire length of borehole BH 1/99 (B1) indicating an ice content of zero in this part of the profile.The water table was found at a depth of around 18 m.In contrast, temperatures in borehole BH 2/99 (B2) showed positive values in the uppermost 5 m (active layer), isothermal conditions near the freezing point down to 20 m depth and unfrozen conditions below.The water table was found around 23 m depth which is close to the lower model domain of the 4PM.
The estimations from the 4PM are in remarkably good agreement with the borehole results except for the low, but non-zero ice contents predicted in the 4PM results near BH 1/99 (B1) (Figs.5a and 6a).The model correctly recognises this region as a low ice content region, but the absolute ice content values appear larger than the borehole results would indicate.At the location of BH 2/99 (B2) the 4PM accurately predicts low ice content values in the uppermost 5 m (unfrozen active layer) and maximal ice contents of 30-40% between 5 and 20 m depth.Below 20 m depth the ice content decreases to values smaller than 15% and the water content values exceed 30%, indicating the water table.Similarly, the water table is accurately predicted for BH 1/99 (B1).
The 4PM results confirm the interpretation from the original geophysical surveys in combination with the borehole validation data, i.e. the location of borehole BH 1/99 within a zone of already degraded permafrost and the location of borehole BH 2/99 at the margin of the remaining ice core shown by the high ice contents in Fig. 5.This highly heterogeneous distribution of ground ice within the rock glacier is further confirmed by results from additional crosshole GPR tomography measurements between the boreholes (Musil et al., 2006).

General model
Figure 8 shows the calculated volumetric fractions for the general model as (a)-(d) minimum values and (e)-(h) maximum values of the solution space.Large differences between minimum and maximum values indicate that model results are poorly constrained by the data, and small differences indicate that the volumetric fractions are well constrained even when porosity is not prescribed.
As a first result, areas with no solution (white regions) are now strongly diminished due to the freely varying porosity.This confirms the above hypothesis, that the decreasing porosity at larger depths (including the occurrence of firm bedrock) is responsible for the inability of the porosity dependent model to find solutions are larger depths (Figs. 5 and 6).Secondly, minimum and maximum values for the air and water content are very similar (note that the colour scale for these two fractions is only between 0-20%) indicating low ambiguities for these two phases regarding the correct prescription of porosity.
Comparing the minimum and maximum values for ice and rock content reveals a trade-off between these two phases.For both phases values between 0-100% are possible indicating the inability of the model to discern between rock and ice over almost the entire model domain.An exception can be found for certain regions at larger depths, where minimal rock contents are between 60-100% and maximal ice contents are between 0-40% (see Fig. 8d,e).In these regions bedrock occurrences are highly probable.Apart from these regions, no results concerning ice and rock occurrences can be determined from the general model due to the ambiguity between ice and rock.

Porosity dependent model and general model
Figure 9 shows the calculated ice-, water-and air contents for rock glacier Murt èl.
In this case, the profile line was oriented longitudinal to the flow direction of the rock glacier, extending over the tongue of the rock glacier (horizontal distance 140-150 m) 801 to the non-permafrost area below (horizontal distance >150 m).Because of this, the porosity model was divided into a rock glacier part (porosity 50%) and a non rock glacier part (>145 m horizontal distance: porosity 30%) for the material outside the rock glacier.The original geoelectric and seismic tomograms are discussed in detail in Maurer and Hauck (2007).
In contrast to rock glacier Muragl, the 4PM results show a uniform ice body within the rock glacier body below an active layer of 3-5 m (Figs.9a and 10a).The lower boundary of the ice body was not found due to the limited penetration depth of the geophysical surveys.Therefore, the homogeneous porosity model seems to be more consistent with the data and regions without physically plausible solutions are restricted to an anomaly within the rock glacier (90 m horizontal distance) and the region in front of the rock glacier.
The calculated air content (Figs.9c and 10c) is only significant within the active layer, reaching maximal values of 20%.Maximal values can be found under the so-called ridges, the characteristic small hills of the rock glacier morphology.Correspondingly, the water content in the active layer is maximal in the troughs, where ice occurrences can be found as well.This feature is in good agreement with theory and direct observations, with increased summer melting under the ridges where enhanced radiation due to the micro-topographic effect is present.In front of the rock glacier tongue, the dominant fraction of the available pore space is made up by water (Figs.9b and 10b).This is in good agreement with surface observations of a small water flow originating at the base of the rock glacier tongue.
As for rock glacier Muragl, the minimum and maximum values of the general model are similar for the ice and water contents, respectively (Fig. 11).The range between minimum and maximum values for the regions with high water and air anomalies is around 15-20%, indicating only a small ambiguity regarding the unknown porosity.On the other hand, rock and ice exhibit the full possible range between 0-100% making a distinction between ice and rock occurrences impossible for most part of the profile.Notable exceptions can be found for the above mentioned anomaly near horizontal distance 90 m, where high minimal rock contents indicate the presence of a large boulder within the rock glacier.Furthermore, high minimal rock contents below the tongue of the rock glacier indicate the presence of bedrock at 15-20 m depth (around horizontal distance 150 m in Fig. 11d,e).

Verification and interpretation
For the interpretation, the analysis of the volumetric fraction values relative to the available pore space is most suitable (Fig. 10).The unsaturated conditions (significant air content) in the active layer are clearly resolved, as is the presence of ground water near the front of the rock glacier (horizontal distance 160), where a spring can be found at the surface.From the relative ice content it is seen that the singular rock occurrence within the rock glacier near 90 m is frozen.
The stratigraphic characteristics of all boreholes available on rock glacier Murt èl are similar, in BH 2/87 (located 50 m upslope from the upper end of the survey line) the following structure was encountered at successively greater depths: 0-3 m -large boulders and air-filled voids, 3-15 m -an ice-rich layer, 15-30 m -a mixture of ice and sand, 30-52 m -a mixture of sands and boulders and bedrock below (Vonder M ühll and Holub, 1992;Arenson et al., 2002).As can be seen from Fig. 10, this in good agreement with the calculated ice and air contents from the 4PM.Results from ice, air and solid particle content measurements within the two other boreholes (additional 10 m away from the survey line) indicate supersaturated ice conditions (up to 90% volumetric ice content) between 5 and 20-25 m depth and air contents of generally less than 10% (Arenson and Springman, 2005).
In comparison to rock glacier Muragl, this rock glacier is still exhibiting a solid supersaturated ice body within most of its morphology -a result which is confirmed by additional electrical resistivity tomography monitoring and cross-profile data sets (Hilbich et al., 2009).803

Discussion
The field data sets from the two rock glaciers showed highly variable distributions of electrical resistivities and seismic velocities including large horizontal and vertical gradients, as well as anomalies with limited extent (especially for rock glacier Muragl).
For interpreting the geolectric and seismic tomograms, we employed the two 4PM approaches (prescribed porosity and general 4PM approach) to compute the spatial water-, ice-and air content distributions, showing a spatially heterogeneous ice distribution for rock glacier Muragl and rather homogenous ice body underneath the active layer for rock glacier Murt èl.Both results are consistent with borehole data.This good overall performance of the 4PM is basically due to the strong contrasts in resistivity between frozen and unfrozen material (delineating the ice body from the unfrozen parts of the rock glacier) and the strong contrast in P-wave velocity between airand ice-filled block layers (delineating the active layer and the non-permafrost material from the permafrost).On the other hand, the general model approach showed clearly the still inherent ambiguity between rock and ice occurrences based on electric and seismic data sets: the reason can be found in the similarities of the two solid phases, ice and rock, concerning their geophysical characteristics.Because the seismic velocities of ice and rock are comparatively similar (Table 1) and their electrical resistivities do not enter the set of equations used in the 4PM a differentiation between ice and rock remains difficult without a priori information.
A notable exception is the detection of firm bedrock below the rock glaciers and the occurrence of very large boulders within the ice core.Due to the markedly different porosity of firm bedrock compared to the on average 50% porosity of the block-ice mixture, the possible solutions of the 4PM are constrained to high rock contents, wherever high seismic velocities were observed, as otherwise no physically plausible solution can be obtained (cf.Fig. 2).Similar unambiguous results of the general model approach were found for water and air contents, which are well constrained by the 4PM, even if the porosity is not prescribed.
Interpretation can be difficult where the model fails to differentiate between ice and rock and no a priori information about porosity is present.An improved formulation of the 4PM, e.g. by using an electrical relationship that includes the resistivity of the bedrock, may overcome this problem.Similarly, a possible bias in the predicted iceand unfrozen water content at larger depths is present for the constant porosity model.
Because the seismic tomogram is based on the assumption of an increasing P-wave velocity with depth, low seismic velocities (corresponding to high air contents) are rarely found at larger depths in the inversion models.From a geological point of view this is often not unrealistic, as the compaction of the ground material usually increases with depth.Similarly, the degree of weathering of bedrock, with a corresponding increased porosity, decreases with depth.In such cases a porosity model with decreasing porosity values with depth has to be applied.However, rock glaciers can sometimes exhibit large air voids at greater depths, where large blocks are still present, but the ice content is diminished.In these cases, the air content can be underestimated by the 4PM.

Conclusions
A new model for quantifying subsurface ice, water and air content based on geophysical data sets has been presented.This so-called 4-phase model ( 4PM) is based on a combination of Archie's Law and an extension of Timur's time-averaged mixing rule for seismic velocities.Tomographic electrical resistivity and P-wave velocity data sets serve as input for the 2-D model for ice-, water-, air-and rock content.Two approaches were applied: 1. a porosity dependent model, where a simple porosity model is prescribed and the remaining 3 phases are calculated relative to the available pore space and 2. a general model approach, where rock content may vary over the model domain to identify ambiguities between the prediction of the 4 phases and to delineate the previously unknown distribution of porosity (i.e.rock content).
Here -Both, the porosity dependent and the general 4PM approach result in similar distributions of water-and air content within the rock glacier bodies.Reliable ice contents can only be determined, if the porosity is correctly prescribed, as the differentiation between rock and ice occurrences in the general 4PM is not always constrained by the geophysical data.However, only the general 4PM approach enables the determination of depth and geometry of the bedrock layer and major rock occurrences within the rock glacier, wherever large porosity contrasts are present.
-The main current limitation of the 4PM is therefore the difficulty in discerning ice and rock due to their similar seismic velocities and the simplified formulation of Archie's Law, where the electrical resistivity of neither rock nor ice is included.An improved 4PM formulation should aim for a more sophisticated electrical mixing rule taking into account the different electrical characteristics of the three nonconducting phases: air, ice and rock.
-For the presented rock glacier case studies the lateral and vertical extent of the ice body as well as the air inclusions between the coarse boulders of the unfrozen top layer and the water table could be reasonable well delineated.The respective calculated volumetric contents coincide well with data from nearby boreholes.
-The calculated vertical and horizontal variability of the volumetric contents is large in the case of Muragl and small for Murt èl rock glacier, both being in good agreement with borehole results and complementary geophysical data sets.
-In both cases a significant air content is only found near the surface and especially underneath the ridges of the rock glaciers.The water content is low except near the troughs and in the unfrozen/degraded areas of the rock glaciers.
Finally, it should be noted, that the 4PM can be seen as visualisation tool for complementary geophysical data sets with respect to the four phases commonly found in permafrost regions.It does not change or improve the original tomographic images.The latter could be achieved by combining the 4PM with a joint inversion approach of seismic and electric data sets.Nevertheless, as the delineation and quantification of ground ice, its spatial variability, as well as the detection of isolated air inclusions, is an important prerequisite for thermal modelling of the future permafrost evolution and for the stability analysis and hazard assessment of frozen mountain slopes, we believe that the 4PM may prove to become a useful tool in permafrost studies.
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

789
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | . The porosity 791 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | 1. a porosity dependent model, where a porosity model as one of the four unknowns is prescribed, and 793 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | : the calculated Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | glacier Muragl shows much more spatial heterogeneity than rock glacier Murt èl.Nevertheless we believe that a comparison between the results for both rock glaciers 797 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

799
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper |

Discussion
Paper | Discussion Paper | Discussion Paper | Discussion Paper | are the key results from model performance tests and two applications to rock 805 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | glaciers:

Fig. 1 .Fig. 2 .
Fig. 1.Volumetric fractions of (a) air, (b) ice and(c) water within the pore spaces as a function of electrical resistivity and seismic P-wave velocity for a porosity of 50% and a set of material properties corresponding to the debris-ice-rock mixture often found at mountain permafrost sites (Table1).

Fig. 5 .Fig. 6 .
Fig. 5. 4-phase model results for rock glacier Muragl calculated with a constant porosity model of 50%.The locations of the two boreholes are indicated by black vertical lines -the depth of the water table by the white horizontal lines.(a) Ice content, (b) water content and (c) air content.