The Cryosphere A new model for estimating subsurface ice content based on combined electrical and seismic data sets

Detailed knowledge of the material properties and internal structures of frozen ground is one of the prerequisites in many permafrost studies. In the absence of direct evidence, such as in-situ borehole measurements, geophysical methods are an increasingly interesting option for obtaining subsurface information on various spatial and temporal scales. The indirect nature of geophysical soundings requires a relation between the measured variables (e.g. electrical resistivity, seismic velocity) and the actual subsurface constituents (rock, water, air, ice). In this work, we present a model which provides estimates of the volumetric fractions of these four constituents from tomographic electrical and seismic images. The model is tested using geophysical data sets from two rock glaciers in the Swiss Alps, where ground truth information in form of borehole data is available. First results confirm the applicability of the so-called 4-phase model, which allows to quantify the contributions of ice-, waterand air within permafrost areas as well as detecting solid bedrock. Apart from a similarly thick active layer with enhanced air content for both rock glaciers, the two case studies revealed a heterogeneous distribution of ice and unfrozen water within Muragl rock glacier, where bedrock was detected at depths of 20–25m, but a comparatively homogeneous ice body with only minor heterogeneities within Murtèl rock glacier. Correspondence to: C. Hauck (christian.hauck@unifr.ch)


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 an 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 enhanced occurrence of slope instabilities, where 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 the modelling of the future evolution of the ground thermal regime in permafrost areas and geotechnical assessment of the hazard potential due to degrading permafrost.Especially temporally and spatially varying ice and liquid water contents in the subsurface will alter the thermal, hydraulic and geotechnical characteristics of the permafrost material and have to be known for initialisation and validation of permafrost models.However, due to the commonly difficult access for permafrost locations in Published by Copernicus Publications on behalf of the European Geosciences Union.
C. Hauck et al.: A new model for estimating subsurface ice content 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 is 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 measured geophysical variables (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 only few petrophysical relationships available for full 4-phase systems, e.g. the approach by Wang and Schmugge (1980) for the dielectrical properties.Simple relations have been developed for the electrical properties, e.g. the well-known relation between the electrical resistivity of the material, the pore-water resistivity, the porosity and the saturation known as Archie's Law (Archie, 1942) and for the elastic properties such as the slownessaveraging 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 (for unconsolidated sediments see e.g.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 com-bine 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 the 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 , liquid 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).The porosity 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 moderate temperatures below zero.In the presence of electrolytic conduction in the pore water, ice can be seen as an insulator similar to the air.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 space can cause large increases in seismic velocity compared to the velocity when the interstitial water is unfrozen, especially in high-porosity sediments and debris (Hilbich, 2010).Since ice is much stiffer than water, the wave velocity is tightly coupled to the iceto-water ratio (Timur, 1968;Carcione and Seriani, 1998).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 properties ρ w , v r , v w , v a , v i as well as the (material dependent) free parameters 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 volumetric ice content f i we obtain: Similarly, equations for the volumetric air content f a and the volumetric water content f w can be derived: In addition to the dependence on the material parameters 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: 1. a porosity dependent model, where a porosity model as one of the four unknowns is prescribed, and 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 parameters 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 properties ρ w , v r , v w , v a , v i 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)  (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 (Table 1).

Fig. 2.
As in Fig. 1, but for a porosity of 5 %. 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).
To simplify our approach and to reduce the number of degrees of freedom of the problem, we used fixed values for the Archie parameter (a = 1, m = 2, n = 2) throughout the whole model grid.King et al., 1988;Hauck and Kneisel, 2008).

4PM calculations in this study
For a porosity of 0.5 (Arenson and Springman, 2005 found porosities between 0.4-0.8 from boreholes in Muragl rock glacier) this leads to a solution space between around 600-4500 m s −1 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 300 m s −1 (Fig. 1a).For the calculation of the water content, an inherent problem of using Archie's Law becomes apparent in Fig. 1c: the calculated 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 −1 (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 −1 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 ( 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 having rather narrow and distinct ranges of possible velocity/resistivity data pairs (very low resistivities and medium velocities for water; very high resistivities and very low velocities for air; cf. also Figs. 1 and 2).
To analyse the implications of these ambiguities, the full solution space of Eqs. ( 6)-( 8) was explored using the Matlab symbolic math toolbox.The solution space comprises the 4PM result (ice, water, air and rock content) for each individual resistivity/velocity pair for a given set of parameters (resistivity of pore water, rock velocity, Archie parameter etc).As the general model is strongly underdetermined, more than one combination of the four phases may result in the same observed resistivity/velocity pair.The calculation of all theoretically possible solutions for the four phases for a given resistivity/velocity data pair of the measured tomograms are then used to analyse the uniqueness of the final 4PM result and to determine how well the four phases are constrained by the data.
Evaluations for several examples of resistivity/velocity data show that water and air contents are usually well constrained due to their marked differences in resistivity and velocity values (shown for two examples in Fig. 3).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 fractions, but not the individual contributions.
To evaluate this for the 2-D tomograms of the geophysical surveys, all possible solutions are calculated for each measured 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).For cases where minimum and maximum values of all possible solutions for the 4 phases con- verge, the general model approach allows the determination of all 4 phases (including the rock content, i.e. the porosity).By this, regions with higher porosity can be differentiated from regions with lower porosity (see below).

Sensitivity to prescribed parameters
As mentioned above, the material properties ρ w , v r , v w , v a , v i and the free parameters in Archie's Law a, m, n have to be prescribed for both model approaches.Whereas wellestablished laboratory results for v w , v a , v i exist (Table 1) and a is set to 1 in most previous studies, values for m, n, ρ w , and v r are material/site dependent and could also vary strongly over the model domain.Figure 4 shows the sensitivity of the 4PM ice content result to uncertainties in these four parameters as well as porosity for the two velocity/resistivity pairs of Fig. 3. Hereby, each parameter was varied individually while the others were held constant.
It is seen that the sensitivity of the ice content estimates is comparatively low for the two Archie exponents m and n, but larger for the pore water resistivity, especially for conductive www.the-cryosphere.net/5/453/2011/The Cryosphere, 5, 453-468, 2011 pore waters.In permafrost studies, this parameter is often comparatively easy to obtain, as small water flows from rock glaciers or moraines are often present in summer which can be sampled and measured.Large uncertainties are present regarding the correct relation between prescribed porosity and P-wave velocity of the rock material, similar to the trade-off between rock and ice content in Fig. 3.For comparatively high measured bulk velocities (black crosses), v r is not sensitive but the ice content depends strongly on the prescribed porosity (the higher the porosity, the higher the ice content).For low measured bulk velocities (red circles) the ice content depends on v r but is insensitive to uncertainties in porosity.This again highlights the importance of consistent estimates for porosity and velocity of the rock material (see also below).

Geophysical inversion
The two approaches outlined above are based on successfully determined data sets of electrical resistivity and seismic P-wave velocity on a 2-D grid.Even assuming that the petrophysical relationships (Eqs. 2 and 5) are justified and that the free parameters can be estimated, the resistivity and velocity models determined by inverted surface measurements are still non-unique and have likely been affected by the inver-sion process.The resolution of tomographic surveys plays a major role in the retrieval of subsurface characteristics.Resistivity values of massive ground ice and of model regions with strong resistivity contrasts (e.g. at the lower boundary of the active layer) can often be plagued by inversion artefacts (Hilbich et al., 2009;Hauck et al., 2003).As the inversion algorithm and the choice of inversion parameters determine how well the inverted model will reproduce the real distribution, reliability studies addressing the influence of inversion parameters on the resulting resistivity and velocity distribution have to be usually performed (Rings and Hauck, 2009;Hilbich et al., 2009;Kneisel et al., 2008).Day-Lewis et al. (2005) refer to this loss of information caused by the inversion process, lack of sufficient prior information and survey geometry as "correlation loss" and computed it with respect to its effect on inferred petrophysical parameters.A similar analysis of the correlation loss of both, resistivity and seismic data, petrophysical relationships in dependence of both inversion algorithms, and its effect on all four phases would be clearly beyond the scope of this paper.However, the inherent dependence of the 4PM results on the quality and ambiguity of the inversion models has to be taken into account.

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 Muragl rock glacier and a shorter 145 m longitudinal profile at Murtèl rock glacier (Fig. 5, see also Maurer and Hauck, 2007).Seismic velocities were obtained using a spread of 120 geophones and 53 shot points at Muragl rock glacier and a spread of 120 geophones and 44 shot points at Murtèl rock glacier 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 Muragl rock glacier two 70 m deep boreholes (BH1/99 and BH2/99) are present along the survey line.One further borehole (BH4/99) is located around 50 m upslope (Fig. 5b).At Murtèl rock glacier two boreholes are present some 50 m upslope of the survey line (Fig. 5d).Ground temperature is measured in the boreholes using a thermistor string attached to a data logger and is available every 6 h, 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).Hereby, the ice content was determined by melting the sample, analysing the volume of the solid particles and calculating the volume of the air and ice by assuming that the specific gravity of ice at a temperature close to the melting point is 0.917 Mg m −3 .The volume of the unfrozen water was neglected.Due to uneven sample surfaces, Arenson and Springman (2005) cautioned that the calculated volumetric ice contents can be under-and the volumetric air contents overestimated for samples with high contents of solid particles.

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 Muragl rock glacier 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) Muragl rock glacier shows much more spatial heterogeneity than Murtèl rock glacier.Nevertheless, we believe that a comparison between the results for both rock glaciers clearly shows the advantages and limitations of the 4PM approach for applications in mountain permafrost terrain.

Muragl rock glacier
Figure 6 shows the results of the electrical resistivity and seismic velocity inversion models for Muragl rock glacier.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 240 m 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 7 shows the calculated ice-, water-and air contents for a constant porosity model of 50 % over the same model domain as in Fig. 6.To obtain the distributions shown in Fig. 7, 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 20-56 % within the uppermost 15 m of borehole BH 4/99.As temperatures in BH 1/99 (B1) do not show permafrost conditions, we assumed a slightly higher content of solids in the left-hand part of the profile 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. 7 delineate those model blocks, where no physically plausible solution of Eqs. ( 6)-( 8) could be found with the prescribed set of parameters, i.e. volumetric fractions were negative, 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. 6b, 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. 7 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. 8 where the results of Fig. 7 are shown in relation to the porosity (e.g. % ice content per available pore space).
As can be seen from Figs. 7c and 8c the calculated air content within Muragl rock glacier 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. 7a) 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. 8a).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 20-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.7b and 8b) is close to zero near the predicted ice core and greater than 30 % below the uppermost 15 m and outside of the ice core.Finally, the non-zero ice contents to the left of B1 are an artefact of the constant porosity model of 50 %, as this region is outside of the rock glacier morphology with firm bedrock present near the surface (cf.Fig. 5b).

Verification and interpretation
Borehole temperatures at the end of summer (Fig. 9, taken from the PERMOS network) show values above the freezing point below 5 m 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 16 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 for BH 1/99 (B1) at depths of around 10 m (Figs.7a and 8a).model correctly recognises this region as a low ice content region, but the absolute ice content values are around 20 % higher 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 (active layer) and maximal ice contents of 30-40 % between 5 and 20 m depth.Finally, the water table is accurately predicted for BH 1/99 (B1) but could not be estimated from the 4PM for BH 2/99 (B2).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. 7.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 10 shows the calculated volumetric fractions for the general model as (Fig. 10a-d) minimum values and (Fig. 10e-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. 7  and 8).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. 10d and 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 11 shows the calculated ice-, water-and air contents for Murtèl rock glacier.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) to the non-permafrost area below (hori-zontal 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 Muragl rock glacier, the 4PM results show a uniform ice body within the rock glacier body below an active layer of 3-5 m (Figs.11a and 12a).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.11c and 12c) 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 microtopographic effect is present.In front of the rock glacier tongue, the dominant fraction of the available pore space is made up by water (Figs.11b and 12b).This is in good agreement with surface observations of a small water flow originating at the base of the rock glacier tongue.Figure 13 shows the results for the general model.As for Muragl rock glacier, the minimum and maximum values give similar results for the air content as well as for the water content.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 (between horizontal distance 130-150 m in Fig. 13d and e).Finally, maximum ice contents (∼100 % saturation, Fig. 12a) are slightly higher for Murtèl rock glacier than for Muragl rock glacier (∼90 %, Fig. 8a), in accordance with Arenson and Springman (2005).

Verification and interpretation
For the interpretation, the analysis of the volumetric fraction values relative to the available pore space is most suitable (Fig. 12).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 Murtèl rock glacier 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. 12, 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 Muragl rock glacier, 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 crossprofile data sets (Hilbich et al., 2009).

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 Muragl rock glacier).For interpreting the geoelectric and seismic tomograms, we employed the two 4PM approaches (prescribed porosity and general 4PM) to compute the spatial water-, iceand air content distributions, showing a spatially heterogeneous ice distribution for Muragl rock glacier and a rather homogenous ice body underneath the active layer for Murtèl rock glacier.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 air-and 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 (cf.Fig. 3): the reason can be found in the similarities of the two solid phases, ice and rock, concerning their  1) and their electrical resistivities do not enter Archie's Law used in the present formulation of the 4PM, a differentiation between ice and rock remains difficult without a priori information regarding the porosity.However, the higher the seismic velocity of the rock material, the larger the contrast between the velocities of ice and rock, and the higher the ability of the general 4PM approach to distinguish between the two solid phases.
A good example of the above 4PM characteristic is the detection of firm bedrock below the two 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 ice-and 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 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 slownessaveraged 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 are the key results from model performance tests and two applications to rock glaciers: -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 non-conducting 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.

Fig. 1 .
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. 3 .
Fig. 3. Relation between the estimates of the four phases for the general model approach without prescribed porosity.All possible solutions are shown for two pairs of P-wave velocity and electrical resistivity.(black) 3500 m s −1 and 100 K m representing high ice content conditions and (red) 1000 m s −1 and 10 K m representing low ice content conditions.

Fig. 4 .
Fig. 4. Sensitivity of ice content estimates to uncertainties in the free Archie parameter and material properties, which have to be prescribed.The sensitivity is shown for the two data pairs of Fig. 3.Each parameter is varied individually while reference values are used for the others.

Fig. 5 .
Fig. 5. Topographical map of the study areas (source: swisstopo) and photographs with borehole locations and survey lines for (a-b) Muragl rock glacier and (c-d) Murtèl rock glacier.Photos by C. Hilbich.

Fig. 7 .
Fig. 7. 4-phase model results for Muragl rock glacier 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 black horizontal lines.(a) Ice content, (b) water content and (c) air content.

Fig. 10 .
Fig. 10.General 4-phase model results for Muragl rock glacier.The locations of the two boreholes are indicated by black vertical linesthe depth of the water table by the black horizontal lines.(Left) Solutions with minimum values and (right) solutions with maximum values.Note, that the range of the colour scales is different.

Fig. 11 .
Fig. 11.4-phase model results for Murtèl rock glacier calculated with a porosity model of 50 % for the rock glacier and 30 % for the permafrost-free area in front of the tongue.

Fig. 13 .
Fig. 13.General 4-phase model results for Murtèl rock glacier.(Left) Minimum values and (right) maximum values.Note, that the range of the colour scales is different.

Table 1 .
Model parameters used for 4PM calculations for the two rock glaciers Murtèl and Muragl in the Eastern Swiss Alps (after Hauck et al.:A new model for estimating subsurface ice content geophysical characteristics.Because the seismic velocities of ice and rock are comparatively similar in general (Table www.the-cryosphere.net/5/453/2011/TheCryosphere, 5, 453-468, 2011C.