1. Introduction
Icefalls pose a considerable risk to humans, settlements and infrastructures. Their destructive power is usually greater in winter, as they can drag snow and ice in their wake. In the Alps one of the most tragic icefall events occurred in 1965 in Switzerland, when a major part of Allalingletscher broke off and killed 88 employees of the Mattmark dam construction site (Reference Röthlisberger and KasserRöthlisberger, 1981; Reference Raymond, Wegmann and FunkRaymond and others, 2003). Following this event, interest in the instability of hanging glaciers grew within the alpine glaciological community. In 1973, the first successful icefall prediction was performed by Reference FlotronFlotron (1977) and Reference Röthlisberger and KasserRöthlisberger (1981) at the Weisshorn hanging glacier, which regularly poses a threat to the village of Randa, in Valais, Switzerland. Due to climatic variations, some hanging glaciers undergo rapid changes, leading either to isolated catastrophic events or to new situations with no historical precedent. It is difficult to perform direct measurements on these steep, heavily crevassed and avalancheendangered glaciers. Measurements are often sparse and fragmentary, and therefore difficult to interpret. They have always been performed after clear signs of destabilization, so the conditions prevailing before an unstable state are unknown. The factors responsible for the destabilization of large ice masses are the strength of the ice and the stresses in the zone of fracture. However, the physics of the ice fracture and the feedback mechanisms between crevassing, ice deformation and load distribution are complex and mostly unknown. Lack of theory and sparse measurements make an accurate stability assessment difficult.
To cope with these difficulties, we developed a numerical model describing the progressive maturation of a mass towards a gravitydriven instability, which combines basal sliding and cracking (described by Reference Faillettaz, Sornette and FunkFaillettaz and others, 2010). Our primary hypothesis was that gravitydriven ruptures in natural heterogeneous material are characterized by a common triggering mechanism resulting from the competition between frictional sliding and tension cracking.
The present paper is devoted to the application of this general numerical model to a particular gravitydriven instability: the breakingoff of hanging glaciers. The gigantic breakingoff of the Altels hanging glacier is an interesting case for several reasons:
It is the largest breakoff recorded in the Alps.
It was well documented by Reference ForelForel (1895), Reference HeimHeim (1895) and Reference Du PasquierDu Pasquier (1896), within the limits of their knowledge – their observations relate to the rupture of the glacier and include photographs of the glacier before and after the collapse.
The causes of this glacier instability are not entirely understood, despite the reanalysis of Reference Röthlisberger and KasserRöthlisberger (1981) using data collected by Reference ForelForel (1895), Reference HeimHeim (1895) and Reference Du PasquierDu Pasquier (1896). However, it is suspected that hot summers prior to the event triggered the breakoff.
It is the only welldocumented breakoff of a cold hanging glacier where progressive warming of the ice/bed interface towards temperate conditions has likely driven the phenomenon.
Any breakoff of this size is of interest in the context of global climate change.
Our study should improve the understanding of the processes leading to this catastrophic phenomenon. Moreover, it should give new insights into the probable causes of this particular gravitydriven instability and reveal precursors of this event. Section 2 describes the study site and gives a qualitative description of the rupture. Section 3 describes how the model is implemented, and the main qualitative and quantitative results are presented in section 4.
2. Study Site: Altels
2.1. Generality and description of the event
The Altels (Berner Oberland, Switzerland) summit is 3629 m a.s.l. and has a pyramidal shape. The northwestern flank is 1500 m high with a 35–40° slope (Fig. 1). It consists of relatively smooth malm limestone (see Fig. 3). In the mid19th century, this face was largely covered with an unbalanced ramp glacier (i.e. the snow accumulation is mostly compensated by breakoff; Reference Pralong and FunkPralong and Funk, 2006), located between 3629 and 3000 m a.s.l. In the early morning of 11 September 1895, a large part of the glacier broke off and tumbled down. This icefall lasted ∼1 m in and the accompanying thunder could be heard in Kandersteg (∼10 km away). Many people thought it was an earthquake. This catastrophic breakoff was carefully described and reported by Reference ForelForel (1895), Reference HeimHeim (1895) and Reference Du PasquierDu Pasquier (1896) and, later, by Reference Röthlisberger and KasserRöthlisberger (1981).
The volume of the breakoff was estimated as 4 × 10^{6} m^{3}, which is the largest known glacial fall in the Alps. The resulting ice avalanche ravaged the high mountain pasture situated underneath (the Spittelmatte) and killed six people and 170 cattle. Due to its high velocity (430 kmh^{−1} (Reference HeimHeim, 1895; Reference Röthlisberger and KasserRöthlisberger, 1981)), the avalanche piled up to 300 m on the opposite slope, towards the Üschinengrat (Fig. 2). An area of ∼1 km^{2} of the pasture was buried under an ice layer 3–5 m thick. As Reference ForelForel (1895) reported, a similar event had occurred at the same place in 1782, killing four people and hundreds of domestic animals (Reference Raymond, Wegmann and FunkRaymond and others, 2003). Reference ForelForel (1895) pointed out that the 1895 summer was warmer than usual. The event is unlikely to occur again, because Altelsgletscher has almost disappeared (a tongue remains on its left side, which will melt away in the near future; Fig. 3a).
2.2. Rupture and probable causes
The crown crack was a huge regular paraboliclike arch with a width of ∼580 m (Figs 3b and 4b). The ice thickness was ∼40 m at the crown crack and 20 m at the glacier terminus (Fig. 3b). The final rupture took place and propagated along the bedrock.
Reference ForelForel (1895) analyzed the causes of the rupture. He suggested the extremely hot previous summers may have been linked to this event. This could explain why the glacier was no longer mostly frozen at the bed and lost adhesion, due to a thin film of water between the bedrock and the ice.
Reference Röthlisberger and KasserRöthlisberger (1981) reanalyzed the documented field observations of Reference ForelForel (1895), Reference HeimHeim (1895) and Reference Du PasquierDu Pasquier (1896) to infer the thermal conditions of Altelsgletscher before its breakoff (Fig. 5). On the basis of Figures 3b and 4b, Reference Röthlisberger and KasserRöthlisberger (1981) concluded that: (1) The central part of the glacier was temperate, as the bedrock was completely icefree in this area after the breakoff. This is an indication that the glacier could slip off in this area. Therefore the glacier was not frozen to the bedrock there. (2) Above the bergschrund, the glacier was frozen to its bedrock, as well as at the margins near the terminus, where some remaining ice (likely still frozen to the bedrock) can be recognized (white spots). He argued that the glacier geometry is likely to explain such thermal conditions at the bedrock. In the central part, a large snowcovered zone can be recognized. In this moderately steep area, surface meltwater could percolate into the cold glacier and refreeze. The release of latent energy leads to an increase in englacial temperatures. If this process is strong enough, it could explain why the glacier bed became temperate in this part. The margins (above the bergschrund and near the terminus) are much steeper and mostly snowfree, so meltwater can hardly penetrate into the glacier and there is no warming by release of latent energy.
Reference Röthlisberger and KasserRöthlisberger (1981) analyzed the retaining forces of the hanging glacier and the possible causes of the breakoff. He found:
Adhesion to the bedrock: Given the dimensions of the unstable ice slab compared to its thickness, the basal properties should be crucial to the stability. The existence of temperate conditions at the bedrock/glacier interface is a destabilizing process that favors sliding of the glacier on its bed. When the glacier slides over bumps or bedrock asperities, cavities form above the bedrock, leading to a loss of adhesion between the glacier and the bedrock. Moreover, surface melting and high water pressure at the interface favor and accelerate the enlargement of the cavities. This could lead to destabilization of the glacier in a very short time (from 1 day to 1 week). Precursors of such a destabilization process are hardly detectable, as the presence of such cavities cannot be easily determined from the surface.
Southern support of the side glacier: The side glacier could also be important to the global stability of the hanging glacier. In Figure 4a it is possible to distinguish opened crevasses on the side glacier prolongated into Altelsgletscher 10 months before the rupture. These crevasses could have formed because of sliding of the side glacier, which acts as a support for the whole glacier. Such a sliding phenomenon was observed in 1895 and again in 1927/28, when the sliding speed was 25–30 m a^{−1} (Reference Röthlisberger and KasserRöthlisberger, 1981).
Support of the frozen margins: It is possible that the lateral support of the frozen margins played a considerable part in the stability of the glacier, especially because of the significant surface melting of the glacier 2 years before the rupture (from a photograph by Reference Du PasquierDu Pasquier, 1896). Reference HeimHeim (1895) and Reference Du PasquierDu Pasquier (1896) analyzed air temperature data from different locations in Switzerland and found that the 1895 summer was exceptionally warm. With the prevailing climate conditions, sufficient meltwater could have been produced to weaken the basal frozen support.
Traction at the crown crack (cohesion): The crown crack was already open, to a limited depth, 10 months before the rupture (Fig. 4a), indicating large tensile stresses.
2.3. Temperature and precipitation before 1895
As Reference ForelForel (1895), Reference HeimHeim (1895) and Reference Röthlisberger and KasserRöthlisberger (1981) pointed out, very warm summers (generating much more meltwater than usual) are suspected to have favored the instability, so we investigated the evolution of climatic conditions at Altelsgletscher before its 1895 breakoff. A database with homogenized continuous daily time series of temperature and precipitation since 1865 is available at the Federal Office of Meteorology and Climatology MeteoSwiss. The time series are corrected for systematic biases (e.g. due to the relocation of weather stations or changing measurement techniques (Reference Begert, Schlegel and KirchhoferBegert and others, 2005)). A total of 12 MeteoSwiss stations are available; we used the two closest stations, one located in Sion (35 km southwest) and the other in Bern (60 km north).
Air temperature is known to be relatively well correlated over large distances (Reference Begert, Schlegel and KirchhoferBegert and others, 2005) and can therefore be extrapolated with confidence. We evaluated the temperature at 3000 m a.s.l. by applying a temperature gradient of −6°C km^{−1} per 1000 m on the daily mean values from the two meteorological stations. The extrapolation of solid precipitation is more problematic. As a first approximation, we computed the mean precipitation from the two stations and considered it to be solid when the extrapolated temperature at 3000 m was <0°C.
Daily snow and ice melt rates can be assumed proportional to the positive degreeday (PDD) (Reference HockHock, 2003). Because the albedo of ice is lower than that of snow, the melting rate will be higher after the annual snow cover has disappeared. Total solid precipitation and the sum of the PDDs are indications of the annual mass balance of the glacier, i.e. low annual solid precipitation rates will lead to an earlier disappearance of the snow cover and therefore to higher melt rates, because of a longer time with low albedo.
Time series of PDDs since 1870 (Fig. 6) indicate higher values in the 1870s and an increasing trend over the 5 years before 1895. Moreover, solid precipitation decreased during this period, indicating less snow accumulation. Five years before the breakoff, the glacier experienced smaller annual mass balances (decrease of solid precipitation combined with a larger PDD). These results are compatible with the descriptions of Reference HeimHeim (1895) and Reference Röthlisberger and KasserRöthlisberger (1981).
3. Numerical Modeling
The aim of the present work is to reanalyze this event by applying a new numerical model designed for describing natural gravitydriven instabilities (Reference Faillettaz, Sornette and FunkFaillettaz and others, 2010). This model allows us to test the different hypotheses proposed to explain the breakoff of Altelsgletscher.
3.1. Model description
We use a model describing the progressive maturation of a mass towards a gravitydriven instability, which combines basal sliding and cracking. Our hypothesis is that gravitydriven ruptures in natural heterogeneous materials are characterized by a common triggering mechanism, resulting from the competition between frictional sliding and tension cracking. Heterogeneity of material properties and dynamical interaction of damage and cracks along the sliding layer seem to have a significant influence on the global behavior.
This numerical model is based on the discretization of the natural medium in terms of blocks and springs forming a twodimensional network sliding on an inclined plane. Each block, which can slide, is connected to its four neighbors by springs that can fail, depending on the history of displacement and damage. We develop physically realistic models describing the frictional sliding of the blocks on the supporting surface, and the tensile failure of the springs between blocks acts as a proxy for crack opening. Frictional sliding is modeled with a state and velocityweakening friction law with a threshold. This means that solid friction is not used as a parameter but as a process evolving with the concentration of deformation and properties of sliding interfaces. Crack formation is modeled with a timedependent cumulativedamage law with thermal activation, including stress corrosion. In order to reproduce cracking and dynamical effects, all equations of motion (including inertia) for each block are solved simultaneously.
This model improves the multiblock model of Reference Andersen, Sornette and LeungAndersen and others (1997) and Reference Leung and AndersenLeung and Andersen (1997) in two ways. First, we use the state and velocityweakening friction law instead of a constant (or just state or velocityweakening) solid friction coefficient. Second, rather than a static threshold for spring failures, we model the progressive damage accumulation via stress corrosion and other thermally activated processes aided by stress. Both improvements make the numerical simulations significantly longer, but have the advantage of embodying, rather well, the known empirical phenomenology of sliding and damage processes. Adding the state and velocitydependent friction law and timedependent damage processes allows us to model the interplay between sliding and cracking between blocks and the overall selforganizing of the system of blocks (Reference Faillettaz, Sornette and FunkFaillettaz and others, 2010).
The geometry of the system of blocks interacting via springs and with a basal surface is shown in Figure 7. In summary, the model includes the following characteristics:

Frictional sliding on the ground.

Heterogeneity of basal properties.

Possible tension rupture by accumulation of damage.

Dynamical interactions between cracks or damage along the sliding layer.

Geometry and boundary conditions.

Interplay between frictional sliding and cracking.
The steps describing how the development of the instability is modeled are shown in Figure 8. As explained above, two phases have to be distinguished:

A quasistatic (quiescent) phase corresponding to the nucleation of block sliding and bond rupture.

A dynamical (active) phase corresponding to the sliding phase of the blocks and the failure of the bonds.
3.2. Geometric parameters
We first have to consider the geometric input parameters for modeling the Altelsgletscher. The glacier is discretized into a system of regular cubic blocks. The application of this model to the time evolution of an unstable glacier implies that snow accumulation at the surface of the glacier is neglected, corresponding to a timescale of 1–2 years. This assumption seems to be justified for Altelsgletscher, as its northwestern face has a steep slope (35–40°) and is subjected to strong winds, drifting snow away. In addition, typical snow and ice accumulation rates are low compared with the overall thickness (∼30 m) of the glacier.
In order to obtain a realistic description of the damage and fragmentation process that may develop in the ice mass, we need a sufficient number of blocks. As a compromise between reasonable sampling and numerical speed, we use a model composed of 70 × 70 blocks for this example. It is possible to evaluate the size of the unstable part from analysis of Figure 4 and from direct observations. Accordingly, the glacier surface area was chosen to be ∼4 km^{2}, with a mean ice depth of 30 m. As we consider a model composed of 70 × 70 blocks, each block corresponds to a discrete mesh 30 m thick, 30 m long and 30 m wide. The weight of each block is ∼24.75 × 10^{6} kg (density 917 kgm^{−3}).
The slope of the bedrock, 𝜙, ranges from 35° (lower part) to 40° (upper part) (Fig. 9). To account for the curvature of the bedrock slope, we used a digital elevation model (DEM) supporting the blocks described above (Fig. 10).
We now describe the two key processes in the model, the friction and damage laws, applied to blocks and bonds, respectively.
3.3. Friction law between the discrete blocks and the basal surface
Figure 3 and the discussion in section 2.2 suggest that the central part of the glacier was temperate at the time of its destabilization. The rupture developed and propagated within the ice (on the margin) as well as above the bedrock (central part). These two zones need to be considered here, depending on where the modeled discrete ice blocks lie (either on ice or on the bedrock). The friction law between the discrete blocks and the basal surface should also capture both cases, ice/ice friction and ice/bedrock friction. Guided by the rate and state friction character of rock (Reference DieterichDieterich, 1994) and ice (Reference Fortt and SchulsonFortt and Schulson, 2009; Reference Lishman, Sammonds and FelthamLishman and others, 2011), a rate and statedependent friction law seems to be adequate to describe both ice/ice friction and ice/bedrock friction. The version of the rate/statevariable constitutive law, currently most accepted as being in reasonable agreement with experimental data on solid friction, is known as the Dieterich–Ruina law (Reference DieterichDieterich, 1994):
where the state variable, θ, is usually interpreted as proportional to the surface of contact between asperities of the two surfaces. The constant μ _{0} is the friction coefficient for a sliding velocity and a state variable θ _{0}. A and B are coefficients that depend on material properties.
The time evolution of the state variable, θ, is described by
where D _{c} is a characteristic slip distance, usually interpreted as the typical size of asperities. The friction law, Equation (1), with Equation (2), accounts for the fundamental properties of a broad range of surfaces in contact, namely that they strengthen (age) logarithmically when ageing at rest, and tend to weaken (rejuvenate) when sliding.
The primary parameter that determines stability, A  B, is a material property. Physically, A B indicates the sensitivity of the friction coefficient to velocity change: a negative value indicates velocityweakening leading to an unstable slip; a positive value indicates velocitystrengthening leading to a stable slip. Reference Faillettaz, Sornette and FunkFaillettaz and others (2010) showed that the case A = B is of special interest because it retrieves the main qualitative features of the two classes, and also because, empirically, A is very close to B. Assuming A = B is therefore reasonable and ensures more robust results. The critical time, t _{f}, signaling the transition from subcritical sliding to dynamical inertial sliding is given in this case (Reference Faillettaz, Sornette and FunkFaillettaz and others, 2010) by
for μ > μ _{0}. (Details of how this equation is obtained are given in the Appendix.) Note that t _{f} → ∞ for μ ≤ μ _{0}, where t _{f} is the time when the block starts sliding, μ _{0} is a constant friction coefficient, A is a constant parameter depending on material properties and θ _{0} is the state parameter at steady state. The parameter μ is evaluated for each block with the definition of the solid friction, μ = T/N, where T is the tangential force determined from the positions of its connected neighbors and N is the normal component of the weight of the block. The state parameter, θ _{0} in Equation (3), is given by (see Appendix)
where is generally interpreted as the initial low velocity of a sliding mass, before it starts to accelerate towards its dynamical instability and where D _{c} can be interpreted as a characteristic slip distance over which different asperities come in contact.
Three parameters have to be correctly evaluated to model the frictional processes within the glacier: μ _{0}, A and θ _{0}. An increasing friction coefficient with increasing sliding velocity (velocitystrengthening) at low velocities has been observed for polished iceonice (Reference Kennedy, Schulson and JonesKennedy and others, 2000; Reference Montagnat and SchulsonMontagnat and Schulson, 2003), iceongranite (Reference Barnes, Tabor and WalkerBarnes and others, 1971) and iceonice along a Coulombic fault (Reference Fortt and SchulsonFortt and Schulson, 2009). A decreasing friction coefficient with increasing sliding velocity (velocityweakening) at high velocities has also been observed (Reference Barnes, Tabor and WalkerBarnes and others, 1971; Reference Kennedy, Schulson and JonesKennedy and others, 2000). The transition between these two regimes occurs at a velocity of 10^{−5} m s^{−1} (Reference Fortt and SchulsonFortt and Schulson, 2009). These two different behaviors are generally explained by two physical mechanisms, depending on the sliding velocity regimes:
The first mechanism is the waterlubrication mechanism (produced by frictional heat at the sliding surface) working at sliding velocities above ∼0.01 m s^{−1}. The waterlubrication mechanism is characterized by the low viscous resistance of a water film produced by frictional heat at the sliding interface (Reference Barnes, Tabor and WalkerBarnes and others, 1971; Reference Kennedy, Schulson and JonesKennedy and others, 2000; Reference Montagnat and SchulsonMontagnat and Schulson, 2003; Reference Maeno and ArakawaMaeno and Arakawa, 2004).
The second mechanism is the adhesion and plastic deformation of ice at the friction interface, which is present at velocities below ∼0.01 m s^{−1} (Reference Kennedy, Schulson and JonesKennedy and others, 2000; Reference Montagnat and SchulsonMontagnat and Schulson, 2003; Reference Maeno and ArakawaMaeno and Arakawa, 2004).
By analogy with rock physics, a rate and statedependent friction model for fresh and sea ice was proposed recently (Reference Fortt and SchulsonFortt and Schulson, 2009; Reference Lishman, Sammonds and FelthamLishman and others, 2011). Such a model enables us to explain the change in the friction coefficient when varying sliding velocity along Coulombic shear faults or slip histories. The experiments performed by Reference Fortt and SchulsonFortt and Schulson (2009) at −10°C on sliding along Coulombic shear faults in ice suggest that the ice/ice friction coefficient is velocitydependent and varies between 0.6 and 1.4, depending on the applied sliding velocity along the fault. As these last results are the closest to our case, we used a mean value of in our calculation.
Once the block slides, the dynamics is controlled by a kinetic friction coefficient, which is in general smaller than the static coefficient, . At low temperatures, Reference Kennedy, Schulson and JonesKennedy and others (2000) and Reference Fortt and SchulsonFortt and Schulson (2009) found a relatively constant friction coefficient around 0.6 for sliding velocities less than 10^{−4} m s^{−1}, and rapidly decreasing values down to 0.1 at velocities more than 10^{−4} m s^{−1}. As velocities more than 10^{−4} m s^{−1} are not expected in our model describing the nucleation phase of the catastrophic rupture, we assume a kinetic friction coefficient, μ _{d} = 0.6.
Next, we turn to the evaluation of coefficient A (and B, since we have made the simplifying assumption that B = A) in the rate and statedependent friction law (Equation (1)). Laboratory experiments suggest that A is smaller than μ _{0} by typically one and sometimes up to two orders of magnitude for rock (Reference ScholzScholz, 1998, Reference Scholz2002; Reference Ohmura and KawamuraOhmura and Kawamura, 2007), ranging from 0.01 to 0.2. A recent experiment on saline ice suggests that A could be significantly higher for ice (A ≈ 0.3; Reference Lishman, Sammonds and FelthamLishman and others, 2011). As we do not have access to strong experimental constraints and since the friction law should be valid for both ice/ice and ice/bedrock friction, we choose A ≈ 0.1, corresponding to onetenth of the static friction coefficient, which is the upper limit given by Reference Ohmura and KawamuraOhmura and Kawamura (2007) for rock. This choice seems reasonable for both cases.
The last parameter to be determined is θ _{0} in Equation (4). In the case of a glacier, the sliding velocity is typically of the order of centimeters per day. Therefore a roughly correct estimate is . D _{c} can be interpreted as a characteristic slip distance over which different asperities come in contact. It is difficult to decide this value. Recent seismological literature reports that D _{c} lies in the range of tens of centimeters to meters for earthquakes (Reference Mikumo, Olsen, Fukuyama and YagiMikumo and others, 2003; Reference Zhang, Iwata, Irikura, Sekiguchi and BouchonZhang and others, 2003). We arbitrarily choose D _{c} ≃ 1 m. Inserting this value in Equation (4), we obtain θ _{0} = 100 days.
To account for the heterogeneity and roughness of the sliding surface, the state variable θ _{i} is reset to a new random value once the dynamical sliding stops. This random value should not be chosen too low, so a block which just stopped sliding does not switch immediately into a new dynamical phase (i.e. t _{f} = 0). Thus, we assign θ _{i} = νθ _{0}, with ν uniformly distributed between 0.5 and 1.5.
Finally, the possible water lubrication of the bedrock leading to a progressive decrease of the ice/rock frictional resistance is modeled using a progressive decrease of the friction coefficient, μ _{0}, at the interface between the simulated ice block and the basal surface. In this way, we account for the progressive onset of the temperate central zone by changing only one parameter.
3.4. Creep law
As discussed in section 3.1, bonds are modeled as linear springs in parallel with an Eyring dashpot. The springs transmit the forces associated with the relative displacements of the blocks. The spring stiffness has also to be evaluated in order to reflect the elastic property of the bulk ice mass. In continuous elasticity, Hooke’s law of elasticity relates stress, σ, and strain, ϵ, via Young’s modulus, E: σ = Eϵ. This leads to the expression: σ _{bulk} = E _{ice}(δL/L). This stress is applied to a side surface of the block, S = HL (where H is the height and L the length of the surface), leading to an equivalent force in the bulk, F _{bulk} = E _{ice}(δL/L)S. A linear spring is subjected to forces given by F _{bond} = K _{bond} δL. In order to find an equivalent behavior, these two forces have to be of the same order, leading to a spring stiffness given by K _{bond} = E _{ice} H. Usually values for E _{ice} are reported to be 9 GPa (Reference Petrenko and WhitworthPetrenko and Whitworth, 1999; Reference PetrovichPetrovich, 2003). However, there is a disagreement of an order of magnitude between measurements of E in the laboratory (9 GPa) and field observations (∼1 GPa), as argued by Reference VaughanVaughan (1995).
Depending on the applied stress, ice has either linear viscous or nonlinear viscous behavior. In glaciers, ice creep is usually described by a nonlinear viscous rheology, called Glen’s flow law (Reference HutterHutter, 1983, and references therein). This law relates strain rate and stress in the secondary creep regime under steadystate conditions. It is thus not possible to use this law to describe tertiary creep and the rupture of a bond. However, Reference HutterHutter (1983) states that the socalled Prandtl–Eyring flow model shows behavior compatible with Glen’s flow law at low stresses. Nonlinear viscous behavior is introduced in our model with an Eyring dashpot in parallel with a linear spring of stiffness, K _{bond} (which is analogous to the Prandtl–Eyring flow model), following Reference Nechad, Helmstetter, El Guerjouma and SornetteNechad and others (2005). Its deformation, e, is governed by the Eyring dashpot dynamics:
where the stress in the dashpot, s _{dashpot}, is given by
Here s is the total stress applied to the bond and P(e) is the damage accumulated within the bond during its history leading to a cumulative deformation, e. P(e) can be equivalently interpreted as the fraction of representative elements within the bond which have broken, so the applied stress, s, is supported by the fraction, 1 − P(e), of unbroken elements. Following Reference Nechad, Helmstetter, El Guerjouma and SornetteNechad and others (2005), we postulate the following dependence of the damage, P(e), on the deformation, e:
where e _{01} , e _{02} and ξ are constitutive properties of the bond material.
Finally, by combining the previous equations, Reference Faillettaz, Sornette and FunkFaillettaz and others (2010) found a creep model that computes the critical time (i.e. failure of the bond) as a function of the stress experienced by the bond, s, given by:
where
and
Creep properties are defined by the parameters K, β, e _{01}, e _{02} and ξ, that we need to fix for our simulations.
We need to find the most appropriate parameters to describe the creep behavior of ice. The ice of natural glaciers has a complex polycrystalline structure composed of crystals of different sizes. From Equation (7), a fraction, 1 – (e _{01}/e _{02}) ^{ξ} , of all the representative elements present undergo abruptfailure immediately after the stress is applied. Since ice does not show significant damage immediately after being loaded, a reasonable assumption is e _{01} = e _{02}. The relative heterogeneity of the material is introduced with the parameter ξ. Despite its complex structure, ice is a fairly homogeneous material compared to a fiber matrix composite. The more homogeneous a material, the greater ξ. In the following, we set ξ = 10 which means that ice is not very heterogeneous; e.g. Reference Nechad, Helmstetter, El Guerjouma and SornetteNechad and others (2005) used ξ ≈ 2 for a fiber matrix composite. Moreover, as e _{01} = e _{02}, it is ξ that influences parameter s ^{⋆}, i.e. the critical stress above which damage starts. The change in s ^{⋆} values is small compared with that of ξ (e.g. ξ = 2 =>. s ^{⋆} = 7.5 × 10^{−4} E, ξ = 5 => s^{*} = 2.4 × 10^{−4} E, ξ = 10 =>. s ^{⋆} = 1.4 × 10^{−4} E).
The other parameters describing the deformation of the Eyring dashpot under an applied stress are β for the nonlinear term and K for the linear term. It is difficult to infer such parameters for ice. Reference Nechad, Helmstetter, El Guerjouma and SornetteNechad and others (2005) used β = 50 × 10^{−9} Pa^{−1} and K = 10^{5} s^{−1} for a fiber matrix composite. As ice is significantly weaker, we arbitrarily choose β = 10^{−7} Pa^{−1} and K = 10^{−3} s^{−1}. With a tensile strength of ice of 1 MPa, we obtain, from Equation (5), de/dt = 10^{−4} s^{−1}, which is coherent with the behavior of polycrystalline ice (Reference Schulson and DuvalSchulson and Duval, 2009).
4. Numerical Results
The aim of the numerical simulations is to test the different causes of the rupture summarized by Reference Röthlisberger and KasserRöthlisberger (1981) and described in section 2.2. In particular, we intend to provide answers to the following:

What is the influence of glacier geometry on the instability?

Can such a rupture occur without changes of the basal properties?

How rapidly does the instability develop?

Can we expect some precursors?
4.1. Description of the simulation
Blocks are distributed in a regular mesh on the inclined plane, so bonds are initially stressfree. At each timestep we evaluate, for each block, the local bed slope from the DEM of the Altels area (Fig. 10). This aims to mimic the real topography of Altelsgletscher. To determine the causes of the glacier collapse, we test the different contributions described by Reference Röthlisberger and KasserRöthlisberger (1981) (section 2.2).
Table 1 summarizes the parameters used in our simulations.
4.2. Qualitative results
Is the 1895 Altels collapse solely due to glacier thickening?
To answer this question, we tested the case of a constant friction coefficient in the presence of a progressive increase of block weight. This means that mass is added at a constant rate on each block. In all the simulations performed, results show that an instability starts from the upper part of the glacier, in contradiction to the observations (Fig. 11). This could be explained by the bedrock topography; the slope is indeed steeper in the upper part, inducing an initial sliding of the upper blocks. Then this instability propagates downwards and the whole glacier collapses.
The progressive thickening is therefore not the cause of the 1895 Altels breakoff event.
Is the 1895 Altels breakoff due to uniform warming conditions at the bedrock?
At such altitudes the glacier is expected to be cold, i.e. frozen onto its bedrock. But we saw above that it experienced successive extremely hot summers, which could have initiated the breakoff. Warming conditions could lead to a lubrication at the bedrock, due to penetration of meltwater and the consequent increase in basal water pressure. This can cause decoupling of the glacier from its bedrock and a decrease in friction between ice and bedrock. There are two ways to model uniform warming conditions in our model: (1) uniformly decreasing the friction coefficient under each block and (2) decreasing the tangential component of the weight of each block. We tested both approaches, and found very similar qualitative results. Again, as expected, the whole glacier collapses, starting from its upper part (Fig. 11), for the same topographical reasons as explained above.
Differential evolution of the basal properties likely affects the stability of this glacier.
Is the 1895 Altels collapse due to a local decrease in the friction coefficient?
We now show that the most likely scenario to reproduce the geometry of the observed rupture is to progressively modify the basal properties in a restricted zone, corresponding to the likely temperate area at the bedrock (Reference Röthlisberger and KasserRöthlisberger, 1981).
In the following, the friction coefficient is decreased at different rates, δμ/δt, on three different surface areas (Fig. 12), corresponding to the assumed temperate bedrock zone. This simulates water penetrating the glacier and lubricating the ice/bedrock interface.
Figure 13 shows the evolution of the set of blocks in the regime where the instability develops. Different phases can be distinguished during this simulation:
In the following, the value of μ _{0} is set larger than arctan(max(𝜙)), where 𝜙 is the slope evaluated from the DEM of the glacier bed. In this way, the glacier is assumed to be stable when the simulation starts. At each numerical timestep, the friction coefficient, μ _{0}, is decreased by δμ/δt.

1. Initially, blocks situated in the zone where the friction coefficient is decreased start sliding where the bedrock slope is largest. Sliding of these blocks leads to a change in stress experienced by the bonds. This internal bond deformation propagates upstream (Fig. 13a).

2. Then the glacier starts accommodating its new stress state, resulting in a quiet phase. The number of sliding blocks progressively increases, leading to stronger interactions and to synchronization of the sliding blocks (Fig. 13d and e).

3. After a certain time (depending on δμ/δt) the lattice starts fracturing. A large crack appears perpendicular to the main slope around the middle of the lattice (Fig. 13b and c). This corresponds to the opening of the crevasse just below the bergschrund (Fig. 4b).

4. The final instability develops. The blocks located below the upper crevasse start to accelerate and a fracture of side bonds propagates in the bedrock slope direction, forming an unstable slab, which finally slides off (Fig. 13f).
Initiation of the instability
For each of the three process zones (i.e. the area where the friction coefficient is decreased) and for different rates of decrease of the friction coefficient (RDFCs), δ(μ)/δ(t), we show the number of sliding blocks at each timestep in Figure 14. It appears that, in all cases, the instability evolves rapidly, typically within 1–2 days. Two different regimes can be distinguished: a first quiescent one, in which isolated blocks slide, and a second active one, in which blocks start to slide in clusters, leading to the final collapse.
The number of sliding blocks after the onset of the instability depends on the RDFC. The greater the RDFC, the larger the number of sliding blocks. When the RDFC is small, the glacier has time to adapt to changes in basal conditions. In this case, the size of the initial unstable cluster is strictly given by the size of the process zone.
Damage evolution within the glacier
In order to measure the damage evolution within the glacier, we count the number of surviving bonds during each simulation. Figure 15 shows the number of surviving bonds within the lattice at each timestep for different process zones and different RDFCs. A rapid increase in the damage a few days after the initiation of the instability can be observed for all simulations. Moreover, this increase does not seem to depend on the RDFC. Results show that, once the behavior enters the active regime, bonds start to fail, leading to the opening of the crown crevasse. That the crevasse opens very rapidly, typically within a few days, was observed by Reference HeimHeim (1895).
Energy analysis
In order to detect precursors to the rupture, we investigated the evolution of the energy stored in the bonds, the kinetic energy (corresponding to the flow velocity) and the radiated energy (i.e. energy released by the rupture of bonds) during the destabilization process (Fig. 16).
Six phases can be distinguished. (1) The glacier remains in a stable phase. (2) Initiation of the instability. (3) The number of sliding blocks drastically increases. (4) The number of sliding blocks reaches a maximum and the energy stored in the bonds increases. Note that the increase of kinetic and radiated energies starts to be visible. (5) The energy stored in the bonds drops because the crown crevasse opens. The number of sliding blocks remains unchanged and the instability is now established. Both the kinetic and radiated energies increase. (6) After an increase of the energy stored in the bonds, clusters of bonds fail.
These results show that no precursor to the instability can be inferred from the time evolution for the energies considered earlier than a few days before the breakoff.
4.3. Quantitative results
The results obtained with the spring–block model are summarized in the following.
Occurrence of the instability as a function of δ(μ)/δ(t) To assess whether the RDFC, δ(μ)/δ(t)), influences the final time of rupture, we performed independent runs with different rates and evaluated their respective rupture times. The results show that the rupture occurs earlier for greater RDFC, which is not surprising (Fig. 17). However, the time of rupture does not depend linearly on RDFC but follows a power law with an exponent, b, of −0.82. The exponent was estimated using the maximumlikelihood fitting method with goodnessoffit tests based on the Smirnov test (e.g. Reference Clauset, Shalizi and NewmanClauset and others, 2009). This means that for small RDFC the glacier has time to adapt to these changes and the final instability arises later than in the case of large RDFC.
The influence of the surface area of the process zone, S, on the time of rupture, t _{i}, was tested and we found an inverse effect compared with the RDFC. Specifically, a glacier for which a large RDFC acts on a small process zone becomes unstable after the same time as a glacier subjected to a small RDFC applied to a large process zone. This can be summarized by plotting the reduced variable t _{i} S ^{0.78} as a function of (δμ/δt)^{−0.82}, as shown in Figure 18.
The determination of these two parameters is not yet possible, especially for the area of the process zone. The collapse of three power laws shown in Figure 18 can be rewritten as t _{i} ∼ [S(δμ/δt)] ^{−ν} , where ν ≈ 0.8. This law expresses a combined ‘size’ effect (through the term S) and a ratedependence effect (through the term δμ/δt). Similar behavior is found in most heterogeneous mechanical systems going to failure (Reference CollinsCollins, 1993; Reference CarpinteriCarpinteri, 1996). Interestingly, the combined dependence of the failure time, t _{i}, on the unstable area, S, and on the rate, δμ/δt, is through their product, suggesting that the driving mechanism for the failure time is the total shear force applied to the unstable area.
5. Conclusions
The Altelsgletcher fall of 1895 is the largest known icefall in the Alps. The mechanisms leading to this event are not fully understood. With a new model, developed by Reference Faillettaz, Sornette and FunkFaillettaz and others (2010), of the progressive maturation of a heterogeneous mass towards a gravitydriven instability, characterized by the competition between frictional sliding and tension cracking, we have contributed to a better understanding of this event. We used an array of sliding blocks on an inclined (and curved) basal topography, which interact via elastic–brittle springs. A realistic state and ratedependent friction law was used for the block/bed interaction. We modeled the material properties of the mass and its progressive damage eventually leading to failure, by means of a laboratorybased stress corrosion law governing rupture of the springs.
Our simulations show that the only way to reproduce the particular arch shape of the crown crevasse was to reduce the basal friction coefficient in a limited area. Such a breakoff arises because of the onset of a weak zone at the interface between the glacier and its bedrock, probably due to infiltration of meltwater trapped within the glacier. Climatic observations indicate that the air temperature increased during summers in the 3 years before the event, supporting this assumption.
A twostep behavior can be seen from our simulations: (1) a quiescent phase without visible changes with a duration depending on the RDFC, followed by (2) an active regime with a rapid increase of basal motion within a few days before the breakoff. As a consequence of the increased basal motion, a crown crevasse opens (as was observed) and the final rupture occurs. This means the destabilization process of a hanging glacier due to progressive warming of the ice/bed interface towards a temperate regime is expected to occur without any visible sign, until a few days prior to the collapse.
The area of the process zone and the RDFC have an equivalent relative influence on the time of onset of the instability. A small process zone area with a large RDFC will lead to the same behavior as a large process zone area with a small RDFC. From a practical point of view, knowledge of both parameters is needed to predict the onset of such an instability. Unfortunately, the a priori determination of these parameters is far from possible, particularly the area of the process zone.
Reference Faillettaz, Pralong, Funk and DeichmannFaillettaz and others (2008, Reference Faillettaz, Sornette and Funk2011) showed that seismic measurements could help to predict the approaching mechanical instability of cold hanging glaciers with some seismic precursors (e.g. by using changes in the statistical behavior of icequakes). Unfortunately, we could not infer any seismic precursors prior to the instability of Altelsgletscher from our modeling results.
In a more general context, climate change may affect the stability of cold hanging glaciers. Moreover, as the rupture process takes some time to develop and external precursors are only visible a few days prior to the breakoff, some cold hanging glaciers could already be in the unstable phase where the instability is developing. Early warning of such events is still far from being possible.
Acknowledgements
Comments by two anonymous reviewers and the editor, T.H. Jacka, substantially contributed to the clarity of the manuscript. This work was partly funded by the European Union seventh framework program’s (EUFP7) project Assessing Climate impacts on the Quantity and quality of WAter (ACQWA; grant 212250).
Appendix
Initiation of sliding for a single block
In section 3.3 we described the subcritical sliding process of a given block interacting via state and velocitydependent solid friction with its inclined basal surface. When the subcritical sliding velocity, dδ/dt, diverges (we refer to the time when this occurs as the ‘critical time’, t _{f}, for the frictional sliding instability), this indicates a change towards a dynamical sliding regime where inertia (the block mass and its acceleration) has to be taken into account.
In this appendix we explicitly calculate how the critical time, t _{f}, is obtained and define its dependence on the parameters and boundary conditions. Let us call the total shear (or normal) force exerted on a given block, where is the force exerted by a neighboring spring bond, and N _{weight} and T _{weight} are the normal and tangential forces due to the weight of the block. We then have
where 𝜙 is the angle of the basal surface supporting the block. Therefore,
As mentioned in section 3.3, A − B is usually very small for a natural material: A− B ≈ ±0.02. For the sake of simplicity, we assume A = B. As discussed in section 3.3, this choice is not restrictive as it recovers the two important regimes (Reference Helmstetter, Sornette, Grasso, Andersen, Gluzman and PisarenkoHelmstetter and others, 2004). This leads to
whose solution is
Combining Equations (2) and (A4), we obtain
Integrating
and, using Equation (A4), we obtain
This expression exhibits the usual regimes: a finite time singularity is obtained for . In this case, Equation (A7) can be rewritten as
With
We can simplify Equation (A8) using the condition that, for μ = tan 𝜙 = μ _{0}, we should have t _{f} → ∞ But, when μ = μ _{0} then exp [(μ – μ _{0}) /A] = 1 and thus, for the condition t _{f} → ∞ to hold, we need
The final expression for the critical time, t _{f}, signaling the transition from subcritical sliding to dynamical inertial sliding is, for μ > μ _{0},
while t _{f} → ∞ for μ ≤ μ _{0}. Note that the dependence on has disappeared due to the relation of Equation (A10).
To summarize, a given configuration of blocks and spring tensions determines the value of and N ≡ N _{weight}, and therefore of μ through Equation (A1). Knowing μ and given the other material parameters, θ _{0}, μ _{0} and A, we determine the time, t _{f}, for the transition to the dynamical regime for that block via Equation (A11).
General algorithm
The simulation of the frictional process for each block proceeds as follows:

1. A given configuration of blocks and spring tensions determines the value of and N ≡ N _{weight} for each block and therefore their solid friction coefficient, μ, corresponds to the ratio T/N.

2. Knowing μ for a given block together with other material parameters (θ _{0}, μ _{0} and A) for that block, the time for the transition to the dynamical sliding regime, t _{f}, is calculated using Equation (3). The value of t _{f} is the waiting time until the next block starts to slide.

3. When the block undergoes a transition into the dynamical sliding regime at time t _{f}, its subsequent dynamics should obey Newton’s law.

4. The dynamical slide of the block goes on as long as the velocity of the block remains positive. When its velocity reaches zero, we assume that the block is no longer sliding. To account for the heterogeneity and roughness of the sliding surface, we assume that the state variable, θ _{0}, is reset to a new random value after dynamical sliding stops. This random value is taken to reflect the characteristics of new asperities constituting the fresh surface of contact.

5. After a dynamical slide, the forces exerted by the springs that connect the block to its neighbors are updated, as is the new gravitational force (if the basal surface has a curvature), the new value of μ is obtained, the time counter for frictional creep is reset to zero, and a new process of slow frictional creep develops over the new waiting time, t _{f}, that is, in general, different from the previous one.
In summary, simulation of the damage process leading to bond rupture between blocks proceeds as follows:

1. Given an initial configuration of all the blocks within the network, the elastic forces exerted by all bonds can be calculated from their extension/compression.

2. For each bond i subjected to an initial stress s _{0}(i), we calculate the corresponding critical time, t _{c,0}(i) at which it would rupture if neither of the two blocks connected to it moved in the meantime. For those bonds where s _{0}(i) ≤ s ^{⋆}, defined in Equation (8), t _{c,0}(i) is infinite.

3. Some bonds will eventually fail, modifying the force balance on their blocks and accelerating the transition to the sliding regime, after which the stresses in the bonds connected to the same blocks are modified.