How to model crevasse initiation? Lessons from the artificial drainage of a water-filled cavity on the Tête Rousse Glacier (Mont Blanc range, France)

Crevasses play a crucial role in glacier-related hazards by facilitating water intrusion into the ice body and potentially triggering the collapse of large ice masses. However, the stress conditions governing their initiation and propagation remain uncertain. In particular, there is ongoing debate regarding the most relevant stress invariants to define fracture initiation (the failure criterion) and the corresponding failure strength, i.e. the stress threshold beyond which crevasses form. Laboratory estimates are hampered by the difficulty of reproducing natural glacier conditions, while in situ studies encounter uncertainties when converting strain or strain rate into stress estimates. This study investigates crevasse initiation processes by analyzing the artificial drainage of a water-filled cavity on Tête Rousse Glacier in 2010. Using the finite element code Elmer/Ice, we simulate the drainage and subsequent cavity refilling over three consecutive years. Given the well-constrained cavity geometry and water levels, stress fields are inferred directly from the force balance, removing the need to convert deformation data into stress estimates. Simulated stress distributions are compared with a pattern of circular crevasses mapped around the cavity after the first drainage event. Our results suggest that crevasses started forming in the delayed non-linear viscous regime governed by the Glen-Nye flow law, rather than as a result of the immediate elastic response to drainage. Additionally, by evaluating four failure criteria commonly used in glaciology, we show that the maximum principal stress criterion, with a stress threshold of 100 to 130 kPa, provides the best match to the observed crevasse field.

1 Introduction

Crevasses are essential components of the cryo-hydrologic system that have significant implications for glacier-related hazards due to their ability to: (i) channel water and the thermal energy it carries from the glacier surface to its bed (Irvine-Fynn et al., 2011; Chudley et al., 2021), with implication on basal friction and, consequently, glacier stability (Faillettaz et al., 2011); (ii) retain water, thereby increasing the risk of outburst floods (Bondesan and Francese, 2023); and (iii) serve as pre-existing fractures that may suddenly propagate through the ice, potentially triggering the collapse of large ice masses (Faillettaz et al., 2015; Chmiel et al., 2023). In addition, crevasses are involved in two processes that control the rate at which ice sheets lose mass into the ocean: iceberg calving (Colgan et al., 2016) and accelerated flow due to the softening of shear margins (e.g., Borstad et al., 2016; Surawy-Stepney et al., 2023; Sun and Gudmundsson, 2023). A proper representation of crevasses in numerical models is therefore a prerequisite for reliable forecasts of glacier instability on the one hand, and of the contribution of ice sheets to future sea level rise on the other hand.

Crevasses are fractures in the ice body that occur when internal stresses exceed a critical threshold marking the onset of failure (Colgan et al., 2016). Although stress fields are routinely computed in ice flow models, establishing a robust link between these stress fields and the occurrence of crevasses is not straightforward. This difficulty arises from the fact that, while failure is relatively well-defined under uniaxial stress, no universal failure theory exists for triaxial stress states. Collapsing a complex three-dimensional stress field into a single scalar equivalent stress, which can be compared to a stress threshold (or failure strength) assumed to be an intrinsic and easily measurable material property, requires the definition of a failure criterion. To ensure independence from the choice of coordinate system, such criteria are typically formulated using combinations of the stress tensor invariants and its eigenvalues. Several failure criteria have been proposed in the literature to explain and model crevasse formation, but no clear consensus has emerged (e.g. Vaughan, 1993; Choi et al., 2018; Mercenier et al., 2019; Chudley et al., 2021; Grinsted et al., 2024; Wells-Moran et al., 2025).

Furthermore, failure strength estimates from laboratory experiments show considerable discrepancies when compared to those derived from in situ observations. Laboratory estimates for tensile strength range from 700 to 3100 kPa for laboratory and lake ice (Petrovic, 2003), while field-based estimates typically fall within the range of 90 to 320 kPa (Vaughan, 1993; Grinsted et al., 2024; Wells-Moran et al., 2025), with the exception of Ultee et al. (2020), who reported a tensile strength on the order of 1000 kPa. Laboratory estimates are likely to be affected by the lack of representativeness of samples (e.g., sample size, crystal fabric, grain size distribution, impurity content, and density) and of the applied stresses, which may not accurately reflect natural conditions (Vaughan, 1993; Petrovic, 2003). In contrast, studies based on in situ or remote sensing observations rely on converting strain rate measurements into stress estimates (Vaughan, 1993; Grinsted et al., 2024; Wells-Moran et al., 2025). This process introduces significant uncertainties due to: (i) the need for assumptions regarding ice rheology; (ii) spatial and temporal averaging of velocities used to produce the strain rate field; (iii) the limitation of observations to the surface; and (iv) the fact that, when advection is significant, the local stress field at a crevasse location may not reflect the conditions under which it was formed.

Because they constitute full-scale experiments characterized by a sudden evolution of stresses which causes the formation of circular crevasses, ice subsidence events have been used to advance knowledge of ice fracture processes (e.g., Evatt and Fowler, 2007; Ultee et al., 2020). For instance, Ultee et al. (2020) analyzed the 2015 eastern Skaftá cauldron collapse (Vatnajökull ice cap, Iceland) following the sudden drainage of a large subglacial lake to estimate the tensile strength of ice. However, due to missing data regarding the exact geometry of the cavity, such studies typically rely on surface deformation observations to infer stresses, or on simplified models that approximate the poorly constrained three-dimensional (3D) geometry of the problem as, e.g., a 2D viscous beam (Evatt and Fowler, 2007) or an idealized circular viscoelastic plate (Ultee et al., 2020). Furthermore, the transient nature of the drainage is generally disregarded as data on the temporal evolution of water level in the cavity is usually unavailable or poorly constrained.

In this study, we use the carefully monitored artificial drainage of a water-filled cavity on Tête Rousse Glacier (Mont Blanc range, France) to constrain the failure criterion and stress threshold that best reproduce the field of circular crevasses mapped at the surface during the summer following the pumping operation. To this end, we performed numerical simulations using the finite-element code Elmer/Ice. Because our model is continuous, it does not explicitly resolve individual crevasses. Instead, what we refer to as crevasse initiation corresponds to the local transition from an intact ice state to a state where fractures emerge at the macroscopic scale. From the perspective of Continuum Damage Mechanics (CDM), the framework commonly used in continuum ice-flow models to represent the averaged effect of crevasses on glacier flow, this essentially means constraining the source term of the damage-evolution equation that governs damage production. Since the evolution of the water level and the cavity volume (and, to a lesser extent, its geometry) are well-constrained, water pressure against the cavity wall can be prescribed as a time-dependent boundary condition, and stress fields are inferred directly from the force balance. In other words, our numerical experiments are entirely force-driven. The resulting deformations serve as independent control variables to assess the nature of the mechanical response, which stands in sharp contrast to studies that infer stress fields from observed displacements or velocities.

The paper is organized as follows: in Sect. 2, we present the study site and the available data.; in Sect. 3, we introduce the model and describe the numerical experiments; results are presented in Sect. 4 and discussed in Sect. 5 in light of previous work.

2 Study Site and Data

2.1 Study Site and historical background

Tête Rousse Glacier is a small glacier of the Mont Blanc massif (French Alps, 45°55^′ N, 6°57^′ E). It is approximately 200 m wide and 500 m long in the east–west direction, covering a total surface area estimated at 0.08 km² as of 2007 (Fig. 1a). Its surface topography ranges from 3110 m above sea level (a.s.l.) to 3260 m a.s.l. This glacier would likely have remained largely unnoticed among the more prominent glaciers of the surroundings, were it not for the catastrophic event of July 1892, when the collapse of a water pocket drained ∼200 000 m³ of water and ice, claiming 175 lives in the village of Saint-Gervais-Le Fayet. Between 2007 and 2010, extensive geophysical surveys revealed the presence of a new subglacial cavity containing an estimated water volume of 53 500 m³. In two boreholes drilled above the cavity, artesian outflows rising 10–20 cm above the glacier surface were observed, indicating that the water pressure exceeded the ice overburden pressure. In light of these findings, public authorities ordered the drainage of the cavity, which was carried out between 25 August and 10 October 2010 using down-hole pumps. A total of 47 700 m³ of water was extracted from the glacier, leaving a few thousand cubic meters that could not be pumped out. In August 2011, circular crevasses were detected at the glacier surface around the cavity (Fig. 1b). The cavity naturally refilled during summer 2011 and needed to be drained again in September of the same year. The same thing happened in summer 2012. Meanwhile, the cavity's geometry evolved towards a reduced volume, due to ice creep during periods of low water levels. Several detachments of ice blocks from the cavity roof were also observed in the field. Volume loss was further accelerated by the partial collapse of the cavity roof on 14 August 2012 (Fig. 1c), which likely caused the refreezing of significant volumes of water in the pores between the ice block debris. Consequently, the cavity's volume decreased from an estimated 53 500 m³ in 2010 to 12 750 m³ in 2013. This volume reduction, combined with the control of water pressure through numerous holes drilled from the surface, led to the decision to cease further drainage operations after 2013. Since the water-filled cavity was discovered in 2010, the Tête Rousse Glacier has been meticulously monitored, leading to the collection of extensive data and the publication of numerous studies (Vincent et al., 2010; Legchenko et al., 2011; Gagliardini et al., 2011; Vincent et al., 2012; Gilbert et al., 2012; Legchenko et al., 2014; Vincent et al., 2015; Garambois et al., 2016).

Figure 1 (a) Elevation contours of Tête Rousse Glacier as of 2011 overlaid on an aerial orthophoto (IGN) from 30 September 2023; the thick black line outlines the modeled domain; the cyan line marks the cavity outline; red lines indicate circular crevasses mapped in summer 2011, whereas magenta lines correspond to crevasses unrelated to the cavity drainage (not shown in subsequent figures); black dots and crosses represent the stake networks used to survey displacements in 2010 and 2011, respectively; blue diamonds denote stakes used for surface mass balance measurements; triangles indicate meteorological stations: the Automatic Weather Station and rain gauge (orange triangles) and the Surface Energy Balance station (red triangle); green stars mark the borehole locations used for the 2010 temperature measurements; the blue pictogram on the northern moraine indicates the approximate position from which photo (b) was taken; the black arrow indicates the general flow direction but note that, away from the cavity, the surface velocity are generally below 1 m a⁻¹. (b) Photograph of the glacier taken in summer 2011, showing circular crevasses (highlighted by yellow arrows) surrounding the cavity roof (approximately outlined by the cyan circle). (c) Photograph taken at the end of summer 2012, after the collapse of a small part of the cavity roof.

2.2 Data

Below, we provide a brief overview of the data used in this study. We refer the reader to the literature introduced above for detailed information on data acquisition.

2.2.1 Surface, bedrock and subglacial cavity topographies

The surface Digital Elevation Model (DEM) of the glacier was obtained from laser-scan measurements conducted over the entire glacier surface on 10 August 2011 (Vincent et al., 2015). The bedrock DEM was derived from several Ground Penetrating Radar (GPR) surveys (Garambois et al., 2016). The cavity topography was reconstructed as follows: (1) an initial estimate of the cavity geometry was derived from sonar measurements performed in September 2010 and September 2011 (Vincent et al., 2015), supplemented by GPR surveys conducted in May 2010, May 2011 and August 2011 (Garambois et al., 2016); (2) manual adjustments were applied to achieve an initial volume of approximately 51 000 m³, which falls between the 53 500 m³ of water estimated from surface nuclear magnetic resonance (SNMR) measurements performed in September 2009 and June 2010 and the 47 700 m³ pumped out of the cavity in 2010 (Legchenko et al., 2014). As a result, the initial position of the cavity within the glacier and its volume are well constrained, but its exact shape remains more uncertain.

2.2.2 Surface mass balance

Between 2010 and 2013, annual and seasonal mass balances were measured at eight stakes on Tête Rousse Glacier (blue diamonds, Fig. 1a). These measurements were used to calibrate a degree-day model which, combined with meteorological data acquired both in situ (triangles, Fig. 1a) and at the Chamonix station, enabled the reconstruction of daily Surface Mass Balance (SMB) across the glacier. More details on these measurements in Vincent et al. (2015).

2.2.3 Ice temperature

Seven boreholes were drilled using hot water drilling in July 2010 along a central longitudinal section of the glacier (Fig. 1a). In each borehole, temperature was measured throughout the entire ice thickness using thermistor chains with an accuracy of ±0.1 °C (Gilbert et al., 2012). A reference temperature field was constructed from these measurements, assuming no variations of temperature in the transverse (south-north) direction. Specifically, a 2D vertical temperature field was generated by projecting temperatures measured at all thermistors onto a common west-east longitudinal cross-section. Linear interpolation was applied between measurement points, while simulation results from Gilbert et al. (2012) were used to extrapolate temperatures in regions beyond the interpolation grid. The resulting 2D vertical temperature field is shown in Fig. S4 in the Supplement.

2.2.4 Water-level measurements

The evolution of the mean water level within the cavity over time was reconstructed from piezometer measurements described in Vincent et al. (2015). Short-term fluctuations were smoothed, and data gaps were filled using linear interpolations (see Fig. 2a in Vincent et al., 2015). The resulting water level curve, that we use as a forcing for our prognostic simulations, is shown in Fig. 2.

Figure 2 Reconstructed water level inside the cavity from the beginning (20 August 2010) to the end (5 September 2013) of the prognostic simulations. Grey shaded areas correspond to pumping operations. Colored dots correspond to days for which equivalent stress profiles are plotted in Fig. 6. The red diamond corresponds to the day for which maps of equivalent stress and pressure are given in Figs. 5 and 9, respectively.

2.2.5 Topographic measurements

The displacements in the vicinity of the cavity roof in response to water-level fluctuations was surveyed using a network of 27 stakes in 2010 (black dots in Fig. 1a) and 29 stakes in 2011 (black crosses in Fig. 1a). Several stakes were lost over the season. Stake positions were measured using a total station with an accuracy of ±0.005 m (Vincent et al., 2015).

2.2.6 Crevasse positions

Circular crevasses were observed around the cavity roof in August 2011 after the snow had melted away (yellow arrows Fig. 1b). These crevasses were mapped using a differential GPS on the 30 and 31 August 2011 (red lines, Fig. 1a). Two of these crevasses, one downstream and one upstream of the cavity, were instrumented with stakes placed opposite each other to monitor their opening and closing. Measurements taken between 9 September and 21 October 2011 confirmed that their behavior was governed by water level fluctuations within the cavity. Specifically, when the cavity was pressurized (9–23 September 2011), the upstream crevasse closed at a rate of up to 5 mm d⁻¹, whereas after drainage (starting 23 September 2011), it reopened at a rate of up to 3 mm d⁻¹. A few additional crevasses were mapped much further upstream of the glacier (magenta lines, Fig. 1a), but they are unlikely to be related to cavity drainage. Although the crevasses were discovered late in the season due to prolonged snow cover, we are confident that they were mapped within less than 1 m of their original formation sites, given the very low surface horizontal velocities ($<\mathrm{1}\mathrm{m}{\mathrm{a}}^{-\mathrm{1}}$).

3 Methods

In this Section, we give an overview of the model, introduce the various failure criteria investigated, and describe the experimental setup. The full set of model equations is presented in Appendix A, while variables, parameters, and constants used in our model are summarized in Table 1.

Table 1 List of variables, parameters and constants used in our model.

3.1 Model description

In this study, our primary goal is to establish a correlation between the occurrence of circular crevasses and the distribution and magnitude of modeled stresses. We are thus interested in modeling the evolution of the stress and displacement fields in response to fluctuations of water pressure within the cavity. This requires solving the momentum balance equation together with a constitutive law, which relates stresses to strain and/or strain rates. Any load applied to polycrystalline ice will result in an instantaneous elastic deformation accompanied by a time-dependent creep (viscous) deformation. Over short timescales, glacier ice is thus usually viewed as a compressible isotropic linear-elastic solid. Beyond the so-called Maxwell timescale (see Sect. 5.1), viscous deformations largely dominate over elasticity and ice is best described as an incompressible non-linear viscous fluid. Some authors have adopted viscoelastic models to capture both regimes within a unified framework (e.g., Christmann et al., 2019; Ultee et al., 2020; Hageman et al., 2024), whereas others have modeled the two regimes separately (e.g., Mosbeux et al., 2020). In this paper, we show that the latter approach is sufficient for our purposes. However, as discussed in Sect. 5, we also performed a few additional viscoelastic simulations (Zwinger et al., 2020), whose results are presented in the Supplement, to further support this conclusion. In all cases, ice is assumed to be isotropic.

The linear elastic regime of isotropic ice is characterized by a linear increase of strain with applied stress according to Hooke's law. The key parameters of Hooke's law are Young's modulus E and Poisson's coefficient ν. Experimental studies investigating the propagation of sound waves in ice have reported Young's modulus values of about 9 GPa (Cuffey and Paterson, 2010). Authors who have used simple elastic beam or plate models to reproduce the observed flexure of ice shelves/ice sheet margins in response to ocean tides have inferred lower values, e.g., E=0.88 ± 0.35 GPa for Vaughan (1995), E=0.8–3.5 GPa for Schmeltz et al. (2002). In their study, Ultee et al. (2020) tested the range E=0.1–10 GPa, with E=1 GPa being reported as the representative value. In contrast, the standard value for Poisson's ratio, ν=0.3, is more widely accepted among authors (Cuffey and Paterson, 2010). It is important to note that for any ν<0.5, ice is implicitly treated as a compressible material. In this case, strict mass conservation would require updating the density to account for volume changes. However, we assume these deformations to be sufficiently small that density can be considered constant. We solve the resulting system of equations for the three components of the displacement field. The domain geometry evolves instantaneously in response to changes in external forces and is directly inferred from the displacement field.

The viscous behavior of ice is described by the Glen-Nye flow law, a Norton-Hoff-type relationship that, assuming ice incompressibility, relates deviatoric stresses to strain rates through a stress-dependent viscosity, making it inherently non-linear. Viscosity is highly sensitive to temperature. Furthermore, Continuum Damage Mechanics (CDM) can be used to account for the degradation of ice mechanical properties at the mesoscale caused by the nucleation of cracks that are too small to be explicitly resolved in the continuum model. This is achieved by prescribing a dependence of viscosity on a physical state variable known as damage (e.g., Krug et al., 2014). Damage initiates in regions where the prescribed failure criterion is met, is advected with the ice flow, and feeds back on ice viscosity (see Appendix A). In this work, we use the damage model implemented within Elmer/Ice by Krug et al. (2014). Unlike the elastic response, where deformation is instantaneous, the viscous regime is characterized by a transient evolution of the geometry toward a new steady state following a change in external forces. The strain rates associated with these time-delayed viscous deformations continuously feed back on ice viscosity, and thus on the stress field, through the non-linear rheology of ice. We incorporate the geometry evolution into the model by solving a free surface equation on the upper surface, which accounts for flow divergence and surface mass balance. To capture the dynamic evolution of the subglacial cavity in response to variations in water pressure, we also treat the lower surface of the glacier as a free surface (neglecting accretion/melting) with a lower bound given by the bedrock elevation. Concretely, the lower surface of the glacier is either in contact with the bedrock or forms the walls of the cavity. When the water level in the cavity rises, previously grounded ice may begin to float. Conversely, if the water level drops, floating ice can re-ground. This contact problem is analogous to the dynamics of a marine ice sheet's grounding line, the boundary where grounded ice begins to float and forms an ice shelf. To accurately simulate the evolution of the cavity geometry in transient simulations, we employ dedicated libraries developed within Elmer/Ice that have been specifically designed to address this type of problem (Durand et al., 2009).

The computational domain is bounded by three surfaces: the glacier upper surface, the lateral boundary, and the glacier lower surface. We treat the upper surface as a stress-free boundary. At the lateral boundary, we enforce a non-penetration condition in the normal direction, while we apply a free-slip condition in the tangential direction. The lower surface of the glacier is subject to two distinct boundary conditions, depending on whether the ice is in contact with the bedrock or forms part of the cavity walls. Within the cavity, we use the measured water level to impose hydrostatic water pressure while we apply a no-displacement (resp. no-velocity) condition elsewhere. Given the very low observed surface velocities and the presence of cold ice in the lower part of the glacier, we anticipate little to no sliding at the ice/bedrock interface (Vincent et al., 2015). For transient simulations within the viscous framework, solving the contact problem using Elmer/Ice libraries requires prescribing a sliding law (Durand et al., 2009). As a consequence, we implement a linear Weertman law with a very strong friction coefficient. The resulting sliding velocities are systematically much lower than 1 m a⁻¹.

3.2 Failure criteria

We evaluate four of the most frequently used failure criteria in glaciology-related studies, which are introduced below. These criteria differ in how they define the equivalent stress, that is, a scalar quantity that serves as a surrogate for the complex three-dimensional stress state and that can be compared to a scalar failure threshold. Within a CDM framework, the difference between the equivalent stress and this threshold determines whether damage is initiated locally, through the formulation of the source term (Eq. A11) in the damage evolution law. We denote the three Cauchy principal stresses as σ_I, σ_II, and σ_III, adopting the convention ${\mathit{\sigma }}_{\text{III}}<{\mathit{\sigma }}_{\text{II}}<{\mathit{\sigma }}_{\mathrm{I}}$. Stresses are defined as positive in tension and negative in compression. For the remainder of the paper, it is important to keep in mind that, because the upper surface is a free surface, one of the three principal stresses should theoretically vanish at the surface, although small deviations from zero may occur due to numerical discretization. Consequently, the maximum principal stress is systematically either positive or zero (within numerical residuals). An area where the maximum principal stress is zero therefore corresponds to a compressive regime, as it implies that both σ_II and σ_III are negative.

3.2.1 Coulomb criterion

The Coulomb criterion, initially formulated to describe shear failure in rocks, has since been adapted to ice (MacAyeal et al., 1986; Vaughan, 1993; Weiss and Schulson, 2009; Wells-Moran et al., 2025). The Coulomb equivalent stress is expressed as:

$$\begin{array}{}\text{(1)}& {\mathit{\sigma }}_{\text{eq,C}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{I}}-{\mathit{\sigma }}_{\text{III}}}{\mathrm{2}}}+\mathit{\mu }{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{I}}+{\mathit{\sigma }}_{\text{III}}}{\mathrm{2}}},\end{array}$$

where μ denotes the internal friction coefficient, set here to μ=0.1, following MacAyeal et al. (1986) and Vaughan (1993). The first term on the right-hand side represents the Tresca criterion, widely applied in engineering, which is symmetric under tension and compression and does not depend on pressure. The second term accounts for frictional resistance along fault planes, increasing compressive strength relative to tensile strength, and introducing a pressure dependence to the failure criterion.

3.2.2 von Mises criterion

The von Mises (vM) criterion is an energy-based criterion frequently used to model calving (e.g., Morlighem et al., 2016; Choi et al., 2018; Mercenier et al., 2019) and to explain or model crevasse opening (e.g., Vaughan, 1993; Albrecht and Levermann, 2012, 2014; Grinsted et al., 2024; Wells-Moran et al., 2025). Initially proposed within plasticity theory to predict the yielding of ductile materials, this criterion states that yielding begins when the elastic energy of distortion (i.e., the deviatoric component of deformation energy) reaches a critical value. Expressed solely in terms of deviatoric stresses, the von Mises criterion remains independent of isotropic pressure and is symmetric in tension and compression. The von Mises equivalent stress reads:

$$\begin{array}{}\text{(2)}& {\mathit{\sigma }}_{\text{eq,vM}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[({\mathit{\sigma }}_{\mathrm{I}}-{\mathit{\sigma }}_{\text{II}}{)}^{\mathrm{2}}+({\mathit{\sigma }}_{\mathrm{I}}-{\mathit{\sigma }}_{\text{III}}{)}^{\mathrm{2}}+({\mathit{\sigma }}_{\text{II}}-{\mathit{\sigma }}_{\text{III}}{)}^{\mathrm{2}}\right]},\end{array}$$

where deviatoric principal stresses have been rewritten in terms of Cauchy principal stresses.

3.2.3 Maximum principal stress criterion

The Maximum Principal Stress (MPS) criterion, also known as Rankine criterion, defines the equivalent stress solely based on the largest eigenvalue of the Cauchy stress tensor, σ_I. Notably, σ_I can be either tensile or compressive, except at the upper free surface, where it is positive or zero due to the traction-free boundary condition. Most studies applying this criterion to glacier and ice sheet damage consider only tensile stresses, even at depth, because ice is stronger in compression than in tension (Petrovic, 2003). In such a case, the equivalent stress reads:

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_{\text{eq,MPS}}=max(\mathrm{0},{\mathit{\sigma }}_{\mathrm{I}}).\end{array}$$

This criterion is pressure-dependent as the Cauchy principal stresses include the effect of pressure.

3.2.4 Hayhurst criterion

The criterion was initially proposed by Hayhurst (1972) to describe the creep rupture of metallic alloys at a fraction of their yield stress when tested at temperatures between 40 % and 60 % of their melting point. It is arguably the most widely used criterion in glaciology-related applications of Continuous Damage Mechanics (e.g., Pralong and Funk, 2005; Jouvet et al., 2011; Duddu and Waisman, 2012, 2013; Mobasher et al., 2016; Jiménez et al., 2017; Mercenier et al., 2018, 2019; Huth et al., 2021, 2023; Ranganathan et al., 2025). The Hayhurst equivalent stress is given by:

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_{\text{eq,H}}=\mathit{\alpha }{\mathit{\sigma }}_{\mathrm{I}}+\mathit{\beta }{\mathit{\sigma }}_{\text{eq,vM}}+(\mathrm{1}-\mathit{\alpha }-\mathit{\beta })\text{tr}\left(\mathit{\sigma }\right),\end{array}$$

indicating that the Hayhurst criterion is a linear combination of the maximum principal stress criterion, the von Mises criterion, and the isotropic pressure $p=-\text{tr}\left(\mathit{\sigma }\right)/\mathrm{3}$. The parameters α and β control the relative contributions of these terms and must satisfy:

$$\begin{array}{}\text{(5)}& \mathrm{0}\le \mathit{\alpha },\mathit{\beta },(\mathrm{1}-\mathit{\alpha }-\mathit{\beta })\le \mathrm{1}.\end{array}$$

Here, we follow Pralong and Funk (2005) and take α=0.21 and β=0.63.

3.3 Experimental setup

Below, we present the numerical setup and describe the two types of experiments, diagnostic and prognostic, on which this study relies.

3.3.1 Numerical setup

We solve all equations of the model using the open-source finite-element code Elmer (https://github.com/ElmerCSC/elmerfem, last access: 17 March 2026) and its glaciological extension, Elmer/Ice (Gagliardini et al., 2013). We construct two initial computational domains from the dataset introduced in Sect. 2.2.1: one including the subglacial cavity and one without any cavity. In both cases, we extract the upper glacier surface from the 2011 surface DEM. For the domain without a cavity, we extract the lower glacier surface directly from the 2014 bedrock DEM. For the domain including the cavity, we use the reconstructed cavity topography to initialize the bottom free surface in the region where the ice is not in contact with the bed. For both computational domains, we generate the finite-element mesh by vertically extruding a 2D horizontal footprint comprising 5842 linear triangles over 16 vertical layers, from the lower to the upper glacier surface. To ensure comparable meshes, we refine the horizontal footprint around the projected cavity contour even in the cavity-free case, leading to a typical node spacing of 2 m within the expected cavity area, increasing to 12 m at the glacier's lateral boundaries. As a result, the final mesh consists of 53 727 nodes (93 472 elements) for both computational domains. To optimize CPU time, we divide the mesh into 8 partitions for parallel computing.

3.3.2 Diagnostic experiments

We conduct a first set of diagnostic experiments to assess the signature of an empty cavity on the surface stress field. This approach serves as a means of quantifying the response of the stress field to the instantaneous drainage of a cavity initially in pressure equilibrium, which we consider to be an unrealistic worst case scenario in terms of stress shock. We run all simulations in pairs: each simulation is first conducted on the cavity-free domain and then repeated on the domain containing a cavity. In the latter case, we assume the cavity to be empty, meaning that water pressure is set to zero throughout the cavity. We consider both the elastic and viscous frameworks. For the elastic framework, we test two end-member values of Young's modulus that bound the typically accepted range: E=1 GPa and E=9 GPa. We argue that using lower Young's modulus values to fit observations in a purely elastic model may overlook the viscous component of the observed deformation (see Sect. 5.1). For the viscous framework, we consider first the usual non-linear Glen-Nye flow law, i.e. n=3 in Eq. (A5). We conduct a first pair of simulations using the reference temperature field introduced in Sect. 2.2.3. To assess how temperature influences the results, we repeat the pair of diagnostic simulations twice with uniform ice temperatures of T=0 °C and $T=-\mathrm{2}\mathrm{°}\mathrm{C}$, respectively. In addition, we consider the case of a linear viscous law, i.e. n=1 in Eq. (A5), instead of the non-linear Glen-Nye law. Since the parameters of the Arrhenius law (Eq. A8) were derived for n=3 (Cuffey and Paterson, 2010), adopting n=1 requires an adjustment of the fluidity. We test two uniform fluidity values: A=0.4 ${\mathrm{MPa}}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$ and A=0.8 ${\mathrm{MPa}}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$. We choose these values so that the computed velocities remain of the same order as those obtained using the non-linear Glen-Nye flow law. We emphasize that these linear viscous simulations are not intended to reproduce realistic glacier behavior, as linear viscosity is an inadequate representation of ice rheology. Rather, their purpose is to illustrate the role of the non-linear, rate-thinning nature of ice in stress concentration that may lead to crevasse formation.

We carry out a separate set of diagnostic experiments within the elastic framework to quantify the expected purely elastic deformation associated with the drainage of the water in the cavity between two specific dates. We run pairs of simulations on the computational domain with cavity, applying water pressure corresponding to the reconstructed water level at considered pair of dates. We then subtract the displacement field computed for the later date from that computed for the earlier one. While these simulations neglect the viscous evolution of the geometry over time, the resulting displacement field offers an estimate of the order of magnitude of the total elastic deformations expected between the two dates.

3.3.3 Prognostic experiments

To simulate the evolution of the stress field and deformations over time in response to cavity water pressure changes, we perform transient simulations within the viscous framework using the non-linear Glen-Nye law with n=3 only, as linear viscosity is not expected to realistically represent ice behavior. These simulations run continuously from 20 August 2010 to 5 September 2013, covering the three drainage sequences of 2010 (25 August to 10 October), 2011 (26 September to 9 October), and 2012 (23 September to 9 October), along with the refilling phases in between (Fig. 2). The simulations start six days before the first pumping sequence to allow for a short relaxation period, reducing deformation anomalies linked to the initial model configuration. These anomalies notably arise from the upper surface being representative of August 2011, while the initial water level and cavity volume correspond to that of August 2010, along with uncertainties in the exact cavity shape. Since the cavity was out of pressure equilibrium when discovered, as evidenced by artesian outflow observed in two boreholes, a longer spin-up period would have led to cavity expansion and downstream migration, resulting in an initial geometry that would deviate too much from observations. All pumping sequences are simulated with a one-day timestep, while the periods in between are simulated with a five-day timestep. As discussed in Sect. 5.1, these timesteps are longer than the Maxwell time expected near the cavity, and are therefore consistent with our choice, for these transient experiments, to focus on the purely viscous mechanical response without attempting to resolve the short-lived elastic response. We apply water pressure in the cavity, i.e. at nodes on the basal boundary which are not in contact with the bedrock. We derive this pressure from the reconstructed water level shown in Fig. 2, following Eq. (A16). We apply the reconstructed daily SMB introduced in Sect. 2.2.2 across the glacier surface.

We perform a reference simulation without including damage to identify the stress conditions leading to crevasse initiation in intact ice, which, from a CDM perspective, corresponds to the binary switch from undamaged ice (D=0) to the onset of damage (D>0). We run this simulation using the reference temperature field introduced above. We then repeat it replacing the spatially-distributed temperature field with uniform temperatures of T=0 °C and $T=-\mathrm{2}\mathrm{°}\mathrm{C}$. Following this, we run a set of simulations based on the reference simulation with the damage model enabled, considering each of the four failure criterion introduced above alternatively. These simulations include the feedback of damage on viscosity according to Eqs. (A7)–(A9). Since the amount of damage created depends critically on the values assigned to the damage enhancement factor B and stress threshold σ_th according to Eqs. (A10) and (A11), we test eight combinations of these parameters for each failure criterion. We generate these combinations using a Latin Hypercube Sampling (LHS), with σ_th spanning the range 50–200 kPa and B spanning the range 0.5–2.0 (Krug et al., 2014). The resulting values of B and σ_th are summarized in Table 2.

Table 2 Combinations of parameters B and σ_th tested.

4 Results

4.1 Signature of the cavity on diagnostic surface stresses

Figure 3 shows the signature of an empty cavity (i.e., the difference between the solutions obtained in the empty-cavity and no-cavity configurations) in terms of the diagnostic maximum principal stress computed at the surface of the glacier for the three constitutive laws considered. Results are presented for E=1 GPa, $A=\mathrm{0.4}{\mathrm{MPa}}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$, and the reference temperature for the linear elastic, linear viscous, and non-linear Glen-Nye laws, respectively. As shown in Figs. S5–S7 in the Supplement, the computed diagnostic stress field is independent of the law parameters within each of the three constitutive laws; only the resulting displacements/velocities are affected. Specifically, Young's modulus for the linear elastic law and fluidity for the linear viscous law act as scaling factors for displacements and velocities, respectively. For example, elastic displacements computed with E=1 GPa are exactly 9 times those obtained with E=9 GPa. Similarly, velocities computed with $A=\mathrm{0.8}{\mathrm{MPa}}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$ are twice those obtained with $A=\mathrm{0.4}{\mathrm{MPa}}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$. The same applies to the non-linear Glen-Nye law, where different temperature fields lead to different velocity fields, while the stress field remains unaffected. As discussed in Sect. 5, the fact that the computed stress field is independent of the constitutive law parameters stems from the fact that our experiments are force-driven rather than displacement or velocity-driven.

Figure 3 Maximum principal stress anomaly modeled at the surface due to the presence of an empty cavity and for the three constitutive laws considered: compressible linear elastic with ν=0.3 (left), incompressible linear viscous (middle), and incompressible non-linear Glen-Nye law (right). The cavity contour is shown in cyan. Circular crevasses mapped in summer 2011 are marked in green.

Regardless of the constitutive law, the presence of an empty cavity induces an anomaly in the surface stress field, characterized by negative values on the cavity roof and a surrounding crown of elevated tensile stresses. The negative anomaly is below 5 kPa in the elastic case, below 30 kPa in the linear viscous case, and below 60 kPa when using the Glen-Nye law. We recall that, at the surface, the maximum principal stress is either positive or zero. Areas of negative anomaly therefore correspond to zones where σ_I is slightly positive in the no-cavity geometry but reaches zero (up to numerical residuals) in the empty-cavity geometry. This is highlighted in Figs. S5–S7 which show σ_I for the empty-cavity configuration only. We also recall that σ_I=0 implies that both σ_II and σ_III are negative, indicating that the cavity roof becomes strongly compressive at the surface in the empty-cavity geometry.

The crown of elevated tensile stress surrounding the cavity is continuous when stresses are computed using the Glen-Nye law and shows a good agreement with the mapped circular crevasses. In contrast, it appears discontinuous over narrow regions in the linear viscous case and over broader regions in the elastic case. In the latter, some circular crevasses are located in areas without any apparent stress anomaly. The magnitude of the anomalies is similar for the elastic and linear viscous cases, with peak stress anomalies reaching approximately 40 kPa for the linear elastic response and 50 kPa for the linear viscous response. In contrast, the Glen-Nye law results in significantly larger stress anomalies that extend over a much wider area. In fact, except for the very upstream and downstream regions, the entire glacier surface is affected by the presence of an empty cavity. In this case, the highest recorded maximum principal stress anomaly reaches approximately 150 kPa.

4.2 Modeled surface displacements

Figure 4 compares the three components of total displacement measured at stakes between 14 September and 6 October 2010 (i.e., during the 2010 pumping sequence) with those simulated using the elastic framework with E=1 GPa and the Glen-Nye law with the spatially-distributed T field. Table 3 summarizes the minimum and maximum displacements measured and simulated above the cavity over the three periods. While both constitutive laws yield similar displacement patterns indicative of subsidence of the cavity roof, the amplitude of the elastic displacements is systematically two orders of magnitude too low: a few millimeters whereas measurements indicate displacements of the order of 10 cm over the period. As reported above, simulations with E=9 GPa yield displacements that are nine times smaller and therefore remain systematically below 1 mm. In contrast, in line with Gagliardini et al. (2011), the Glen-Nye law produces displacements that are very close to that observed in the x and y directions, although vertical displacements are slightly overestimated by the model. The peak of vertical displacement in the model is also slightly downstream compared to that measured at the stakes, which may indicate inaccuracies in the reconstructed cavity shape. Given the temperature field in the vicinity of the cavity, replacing the spatially-distributed T field by a uniform temperature of T=0 °C (resp. $T=-\mathrm{2}\mathrm{°}\mathrm{C}$) produces slightly larger (resp. slightly smaller) displacements (Table 3).

Figure 4 Total displacements in the x (top row), y (middle row) and z (bottom row) directions, obtained from (left) diagnostic simulations based on the linear elastic framework with E=1 GPa and ν=0.3, (middle) transient simulation using the Glen-Nye law with the spatially distributed T field, and (right) linearly interpolated from total station measurements at stakes between 14 September 2010 and 6 October 2010. Note that displacement units are in millimeters for the elastic case and in centimeters for the other cases, and that the colorscale is not the same across the horizontal and vertical directions. In the right column, black dots indicate positions of stakes. The cavity contour is shown in black. The general flow direction away from the cavity is indicated by the black arrow at the top of the figure.

Table 3 Minimum and maximum displacements above the cavity in the x (top lines), y (middle lines), and z (bottom lines) directions obtained over three periods from diagnostic simulations based on the linear elastic framework, transient simulation using the Glen-Nye law, and total station measurements at stakes. All displacements are given in centimeters.

Two similar figures, covering the periods from 9 to 28 September 2011 (when the cavity was filled with water and under pressure) and from 28 September to 21 October 2011 (spanning the 2011 pumping sequence), are available in the Supplement (Figs. S8 and S9). Overall, the same observations hold as for the 2010 pumping sequence described above. We note, however, that the overestimation of the rate of subsidence of the cavity simulated by the viscous model is higher for the 2011 pumping sequence than for the 2010 pumping sequence. This could be attributed to uncertainties regarding the cavity shape that increase over time, particularly with the detachment of ice blocks from the roof. This detachment, observed in the field but not accounted for in the model, would have lightened the roof of the cavity, thereby reducing the gravitational forces compared to those estimated in the model.

We also note that the model based on the Glen-Nye law is less effective at capturing the cavity roof uplift observed in September 2011 than it is at representing the subsidence during the pumping sequences. This uplift episode corresponds to a period when the cavity was filled with pressurized water exerting an upward force on the roof, as visible in the stake measurements (Figs. 7 and S8). The measured displacement field, however, is significantly noisier than during subsidence periods, with a few stakes showing subsidence while neighboring stakes exhibit uplift. The elastic model, on the other hand, does capture the uplift but, as before, underestimates the displacement by one to two orders of magnitude.

To summarize, these results show that the deformation observed at the glacier surface in response to pressure variations in the cavity is of viscous type, with the elastic component being largely negligible over the time scales of interest. As a consequence, the results presented in the following of this section rely solely on the Glen-Nye law. The relevance of performing viscoelastic simulations to capture the transition from the purely elastic to the purely viscous regime is discussed in Sect. 5.1.

4.3 Failure criterion and Stress Threshold

In Sect. 4.1, we analyzed the diagnostic surface stress field under the assumption of instantaneous cavity drainage. This represents an unrealistic worst-case scenario because: (1) in reality, cavity drainage is not instantaneous (for example, the 2010 pumping operations lasted 45 d); and (2) the resulting stress field is transient as stresses are redistributed over time by viscous deformation of the domain. Here, we present results from transient simulations performed with the Glen-Nye law, which capture the time-dependent evolution of stresses in response to water pressure fluctuations in the cavity.

Figure 5 shows the equivalent stress at the surface computed at the end of the 2010 drainage with the four failure criteria introduced in Sect. 3.2. Both the MPS and Hayhurst criteria produce a crown of elevated σ_eq surrounding the cavity roof. In contrast, the area directly above the cavity roof is characterized by a circular region of strictly zero σ_eq with the MPS criterion, while the Hayhurst criterion yields reduced but nonzero values in this zone. Conversely, the vM criterion produces very high σ_eq values in this region, with the highest values occurring in crevasse-free areas. Similarly, although the Coulomb criterion results in overall lower σ_eq values compared with the vM criterion, the same observation holds: crevasses occur in regions of relatively low σ_eq, while the highest values of σ_eq are found over the crevasse-free cavity roof. This behavior contrasts sharply with the Hayhurst and MPS criteria, for which crevasses are systematically located in regions of maximum σ_eq.

Figure 5 Equivalent stress at the surface at the end of the 2010 pumping operation, computed using the four failure criteria presented in Sect. 3.2. The cavity contour is shown in cyan. Circular crevasses mapped in summer 2011 are marked in red. Transects AA', BB', and CC', along which equivalent stresses are plotted in Fig. 6, are indicated in white.

Figure 6 illustrates the evolution of surface equivalent stress along three transects throughout the 2010 pumping sequence for the four considered failure criteria. A corresponding figure for the period between the 2010 and 2011 pumping operations, during which the cavity refilled, is available in the Supplement (Fig. S10). When discovered the cavity was under pressure as demonstrated by the strong positive maximum principal stress (∼100 kPa) prevailing over the cavity roof at the very beginning of the simulation. Immediately after the drainage begins, the maximum principal stress shifts to zero, indicating a compressive stress state (i.e., σ_II and σ_III are negative). In crevassed regions, σ_I is high (σ_I≳90 kPa), and the other nonzero principal stress is of second order in the calculation of both σ_eq,vM and tr(σ). As a consequence, in these regions (and only there), Eq. (4) yields σ_eq,H≈σ_I (first and third rows of Fig. 6). This explains the similarity between the patterns obtained with the MPS and Hayhurst criteria in the crevassed areas and highlights that σ_I alone is the relevant quantity for predicting crevasse formation in this case. It is also noteworthy that the equivalent stress profiles obtained at the surface with the Coulomb and von Mises criteria are similar in shape but differ in magnitude, with the von Mises criterion producing much higher equivalent stresses.

Figure 6 Evolution of surface equivalent stress along transects AA' (left column), BB' (middle column), and CC' (right column) throughout the 2010 pumping sequence, for the four failure criteria: MPS (first row, Eq. 3), von Mises (second row, Eq. 2), Hayhurst (third row, Eq. 4), and Coulomb (last row, Eq. 1). Transects are reported in Fig. 5 (white lines). Vertical lines indicate crevasses, while the grey shaded area represents the cavity roof. For each day when a profile is plotted, the corresponding water level is marked by a colored dot in Fig. 2.

Crevasses are systematically found in areas corresponding to peaks of σ_I (i.e., the crown of high tensile stresses described above). Only the upstream end of the CC' transect and, to a lesser extent, the downstream end of BB' transect exhibit a slight shift between crevasse positions and stress peaks. As for vertical displacement peaks, this may reflect uncertainties in the exact shape of the cavity. In addition, the value of equivalent stress at which crevasses tend to occur is remarkably robust: crevasses are found in places where σ_I is between 100 and 130 kPa. Only the crevasse at the upstream end of the AA' profile seems to occur at higher stress (σ_I≈165 kPa) if we focus on the stress field computed at the end of the pumping operation. However, this crevasse could have started to initiate earlier in the pumping operation, when peak stresses were lower (around mid September 2010, Fig. 6), although this observation was not reported by those operating in the field. In contrast, as inferred above from Fig. 5, no clear correlation is observed between the presence of crevasses and singularities in equivalent stress when computed using the von Mises or Coulomb criterion.

4.4 Damage feedback on viscosity

The results presented above suggest that a maximum principal stress criterion, combined with a stress threshold σ_th of 100–130 kPa, is the most suitable approach to represent the initiation of crevasses in our model. From a CDM perspective, this provides crucial qualitative information on where damage may initiate, but it does not offer quantitative insights into how fast damage should grow. In the framework of CDM, the damage variable D ranges between 0 (undamaged) and 1 (fully damaged), and directly affects ice viscosity through Eq. (A9). The rate of damage production is proportional to the gap between the local value of the equivalent stress and the stress threshold (Eqs. A10 and A11). The proportionality constant, known as damage enhancement factor B, is an unknown parameter that directly affects the degree of ice softening expected locally. The prognostic simulations with damage introduced in Sect. 3.3 aim to determine whether B can be constrained using displacement data from the cavity roof monitored between August 2010 and September 2013. Figure 7 compares the mean vertical velocity measured at stakes above the cavity roof with simulations using the Glen-Nye law, both without and with damage feedback on viscosity. Simulations including damage consider all of the four failure criteria. Each panel corresponds to one of the eight combinations of parameters (B,σ_th) tested (Table 2). We emphasize that, while blue dots represent mean vertical velocities derived from measurements performed at stakes above the cavity, error bars do not represent measurement precision but the range between the fastest and slowest mean vertical velocities recorded across these stakes during each measurement period (Fig. 7). Consistent with the results reported in Sect. 4.2 and regardless of whether damage is accounted for, all prognostic simulations capture remarkably well both the acceleration of cavity roof subsidence during the 2010 pumping and the subsequent slowdown associated with the cavity refill in late spring/summer 2011. As already mentioned above, the slow uplift of the cavity roof observed in late summer 2011 is less well captured, while the rate of subsidence associated with the 2011 pumping operation is exaggerated in the model, likely due to errors in the cavity roof thickness following the detachment of ice blocks. The second striking result is that, independently of the parameter pair (B,σ_th), the damage model based on a MPS criterion, and to a lesser extend the one based on the Hayhurst criterion, produces the same behavior of the cavity roof as the model without damage. This occurs because, with these two criteria, damage develops in the crown surrounding the cavity roof, but it does not have time to advect to the roof itself due to the very low surface velocity. To illustrate this point, maps of damage obtained at the end of the 2010–2011 refill sequence for the four criteria, using parameter combination 7, are provided in the Supplement (Fig. S11). As a result, the viscosity of the cavity roof remains unaffected, making it impossible to constrain the value assigned to B based on this experiment.

Figure 7 Mean vertical velocity above the cavity roof over time, simulated without damage for different temperature conditions: spatially distributed T field (black solid line), uniform T=0 °C (black dotted line), and uniform $T=-\mathrm{2}\mathrm{°}\mathrm{C}$ (black dashed line). Simulations including damage are shown for the four failure criteria: MPS (blue), Hayhurst (orange), von Mises (green), and Coulomb (pink), all using the spatially distributed T field. Each of the eight panels corresponds to a specific combination of parameters (B,σ_th) as listed in Table 2. Blue dots indicate mean vertical velocities derived from stake measurements, calculated between consecutive observation dates. Time intervals range from 4 to 8 d during the three drainage sequences, and from 9 d (summer 2011, before pumping) to 41 d (summer 2012, before pumping) otherwise. Error bars do not represent measurement precision but the range between the fastest and slowest mean vertical velocities recorded across these stakes during each measurement period. Grey shaded areas correspond to pumping operations.

Although both criteria should be discarded based on previous results, it is still noteworthy that, because the von Mises and Coulomb criteria produce damage directly on the cavity roof (Fig. S11), the deformation of the latter is strongly affected by the inclusion of damage when using either of these criteria, provided that the stress threshold is not too high and the parameter B is sufficiently large. More specifically, since σ_eq,C never exceeds 150 kPa, all combinations with σ_th higher than this value (combinations 4, 6 and 8) result in no damage under the Coulomb criterion. In contrast, because σ_eq,vM exceeds 200 kPa over most of the cavity roof, all simulations relying on the von Mises criterion produce damage in this region, which then deforms more readily, particularly when the gap between σ_eq and σ_th is large (e.g. combinations 2, 3, 7) and/or when B is high (e.g. combinations 1, 5, 7, 8). These results illustrate that adopting an inappropriate failure criterion in a CDM framework can lead to large errors in the modeled velocities, which can be of the same order as, or even larger than, errors resulting from uncertainties in the ice temperature field.

5 Discussion

5.1 Elastic or Viscous?

The force balance determines the total stress within a given geometry. For a given set of body and contact forces, only two factors can influence the stress distribution in this geometry: (i) the assumption regarding compressibility, and (ii) the potential feedback of stresses through strains or strain rates on the rheological parameters (that is, the potential non-linearity of the constitutive law). This explains why the diagnostic stress solution obtained on the empty-cavity geometry is independent of Young's modulus in the elastic framework and of fluidity in the linear viscous framework (Figs. S5 and S6). In contrast, the resulting strains and strain rates are highly sensitive to these parameters and can thus be used as control variables for comparison with observations. This stands in sharp contrast to studies that infer stress states in crevassed areas from observed surface strain or strain rates, where the inferred stresses are highly sensitive to assumptions about ice stiffness or fluidity (Vaughan, 1993; Ultee et al., 2020; Grinsted et al., 2024; Wells-Moran et al., 2025). Of course, in transient viscous simulations, feedback mechanisms occur: differences in velocity responses for different fluidities lead to differences in geometry evolution, which in turn affect subsequent stress fields, and so on.

The fact that the diagnostic stress responses to the instantaneous drainage of the cavity obtained using a linear elastic, a linear viscous, and a non-linear Glen-Nye rheology differ indicates that at least one of the aforementioned conditions is not satisfied. Specifically, the discrepancy between the elastic and linear viscous solutions arises from conflicting assumptions regarding ice compressibility: ice is considered compressible in the elastic regime (ν=0.3), and incompressible in the viscous regime. An additional simulation using an elastic law with ν=0.5 (i.e. the incompressible limit) confirms this point clearly (Fig. S3 in the Supplement). Likewise, the non-linear nature of the Glen-Nye law explains why it produces higher and more localized stresses around the cavity roof than those computed with the two other rheologies. This occurs because the effective viscosity decreases in regions of high deviatoric stress, making these areas less capable of sustaining load. As a result, stress is redistributed to surrounding regions where the ice remains stiffer. This raises the question of the nature of the mechanical response that should be considered when modeling the initiation of crevasses in virgin ice: elastic or viscous (Glen-Nye)?

It is important to keep in mind that, while elastic deformations occur instantaneously, viscous deformations are delayed in time, and so is their feedback on ice viscosity. It follows that the stress field computed using the Glen-Nye flow law represents the actual stress state in the domain only after sufficient time has elapsed for viscous deformations to feed back on viscosity. This characteristic timescale corresponds to the Maxwell time, and can be estimated as:

$$\begin{array}{}\text{(6)}& {\mathit{\tau }}_{\mathrm{m}}={\displaystyle \frac{\mathrm{2}\mathit{\eta }(\mathrm{1}+\mathit{\nu })}{E}}.\end{array}$$

Figure 8 shows the distribution of the Maxwell time at the surface and along a vertical profile of the Tête Rousse glacier, assuming an empty cavity, the Glen-Nye flow law, and a Young's modulus of E=1 GPa. The same figure for a Young's modulus of E=9 GPa is provided in the Supplement (Fig. S12). As pointed out in earlier work (e.g., Podolskiy et al., 2019; Hageman et al., 2024), the non-linear nature of the Glen-Nye law results in a highly heterogeneous Maxwell time, which is particularly short in regions of high deviatoric stresses. In our case, the Maxwell time spans four orders of magnitude: while it reaches several months in low-stress regions, it is on the order of ten hours to one day around the cavity. In other words, beyond approximately one day, the stress solution obtained with the Glen-Nye law is representative of the actual stress state prevailing around the cavity. This is confirmed by additional non-linear viscoelastic simulations presented in the Supplement. If crevasses had opened within this timescale after drainage, they would likely have been noticed and reported by field operators because of safety concerns. Yet no such observations were made: crevasses were only detected the following summer after snowmelt, approximately ten months after the drainage operation (although they likely formed earlier but remained hidden beneath snow). Together with the results presented in Sect. 4.1, showing that only the Glen-Nye law produces stress patterns and magnitudes consistent with the observed crevasse field, this supports the idea that crevasse onset occurred in the pure viscous regime.

Figure 8 Maxwell time at the glacier surface (left) and across a vertical profile (right) evaluated for an empty cavity according to Eq. (6) with the Glen's law and for E=1 GPa. The cavity contour is shown in cyan. Circular crevasses mapped in summer 2011 are marked in red. The transect along which the vertical profile is extracted is reported by the dashed black line in left panel.

5.2 What stress threshold?

Results presented in Sect. 4.3 suggest a failure strength on the order of 100–130 kPa. This value is inferred from a specific configuration for which we are fortunate to have a unique and well-documented dataset. We do not claim that this estimate is universal or directly applicable to all settings. First, the stress threshold used in models must be consistent with the adopted failure criterion. Second, for a given criterion, the failure strength may vary slightly depending on factors such as grain size, deformation rate, and temperature (Schulson and Duval, 2009). Nonetheless, our estimate is broadly consistent with failure strengths inferred from other field-based studies in Greenland (Grinsted et al., 2024) and Antarctica (Vaughan, 1993; Wells-Moran et al., 2025), where ice flow is generally faster and ice temperatures generally lower. These studies are discussed in Sect. 5.3. A notable exception is the work of Ultee et al. (2020), who inferred a tensile strength of σ_th∼1 MPa from observations of the 2015 eastern Skaftá cauldron collapse (Vatnajökull ice cap, Iceland). We propose a possible explanation for this discrepancy below.

In their work, Ultee et al. (2020) approximate the cauldron as a circular plate and solve the force balance for this plate assuming a linear viscoelastic rheology. They perform a series of simulations exploring a range of parameters and compare the total deformation (i.e. instantaneous elastic plus delayed viscous) modeled over four days at the ice surface with the observed one. Because the considered plate is radially symetric, whereas the actual ice thickness was not uniform across the cavity roof, the sampled parameter space includes geometrical variables, most notably the ice thickness above the cauldron, which directly affects the load applied to the plate. We argue that adopting a linear rather than a non-linear (Glen-Nye) viscosity leads to an underestimation of the viscous contribution to the observed surface deformation. Indeed, the strong feedback of strain rate on effective viscosity is expected to markedly enhance ice flow over the cauldron roof, an effect that is not captured by a linear viscoelastic model. Consequently, to reproduce the observed total deformation, the model must compensate for the underestimated viscous deformation by overestimating the elastic contribution. Concretely, this overestimated elastic deformation likely arises from the chosen parameter set. In particular, if the Young's modulus and cauldron radius are assumed to be known, larger elastic displacements can be produced by increasing the applied load, for example by adjusting the poorly constrained value of the cauldron's ice thickness, which in turn leads to overestimated stresses.

5.3 What failure criterion?

Results presented in Sect. 4.4 highlight the strong feedback of damage on ice viscosity, which makes the choice of an appropriate failure criterion critical in a CDM framework. In our study, the MPS criterion proves to be the most effective predictor of crevasse initiation. The Hayhurst criterion also performs well, primarily because, in crevassed areas, the Hayhurst equivalent stress nearly reduces to the maximum principal stress, which thus remains the sole relevant predictor. In contrast, the Coulomb and von Mises criteria do not correlate well with mapped surface crevasses. While our result is consistent with the widely accepted notion that Mode I fracture (opening) is the dominant mechanism for crevasse formation in most settings (Colgan et al., 2016), it challenges findings based on glacier-scale in situ and remote sensing observations (Vaughan, 1993; Grinsted et al., 2024; Wells-Moran et al., 2025). These studies apply a similar approach, converting observed surface strain rates into deviatoric stresses using the Glen-Nye law, and constructing failure envelopes that enclose stress states associated to non-crevassed areas while excluding stress states associated to crevassed areas. Applying this method across 17 polar and alpine locations, Vaughan (1993) found that the Coulomb and von Mises criteria, combined with a failure strength between 90 to 320 kPa, best matched observations. Comparable results were reported by Wells-Moran et al. (2025) for Antarctic Ice Shelves, who refined the failure strength range to 202–263 kPa (for n=3). Similarly, Grinsted et al. (2024) applied the method to the Greenland Ice Sheet and found that a von Mises criterion with a failure strength of 265 ± 73 kPa provides an adequate fit to the data.

As discussed by Wells-Moran et al. (2025), we suggest that these findings are influenced by the fact that they rely on observations from settings where isotropic pressure is negligible compared to deviatoric stresses. As a result, the role of isotropic pressure in preventing crevasse opening cannot be used to rule out pressure-independent criteria. For example, as Wells-Moran et al. (2025) point out for Antarctic ice shelves, a von Mises criterion would likely not fit basal crevasse data well if such data were available, since overburden pressure is expected to act against crevasse opening. In our case study, although we focus on surface stresses, the subsidence of the cavity roof induces a stress pattern in which the surface isotropic pressure deviates significantly from zero (Fig. 9). It is indeed important to keep in mind that isotropic pressure, defined as the mean of the Cauchy normal stresses, generally differs from the overburden pressure ρ_ig(z_s−z). Specifically, the surface pressure over the cavity roof shifts from slightly tensile (i.e., negative) at the beginning of the simulation (Fig. 9b), when the cavity is pressurized, to largely compressive (i.e., positive) by the end of the 2010 pumping operation, when the cavity is empty (Fig. 9a,b). In this region, no crevasses are observed, yet the von Mises criterion still predicts high equivalent stresses due to its pressure independence. A similar behavior is observed for the Coulomb criterion, which is mostly (but not fully) pressure-independent given the low value of μ=0.1 adopted here. This finding contradicts the applicability of pressure-independent criteria and instead supports the concept of a failure strength that increases with pressure, consistent with the experimental results of Nadreau and Michel (1986) for confining pressure values typically found in glaciers.

Figure 9 (a) Same as Fig. 5, but showing pressure instead of equivalent stress. (b) Same as Fig. 6 but showing pressure instead of equivalent stress. Note that we adopt the convention whereby the pressure p is defined as positive in compression and negative in tension.

Interestingly, Grinsted et al. (2024) found that a Schmidt-Ishlinsky criterion outperforms the von Mises criterion in describing glacier ice failure. The Schmidt-Ishlinsky criterion is essentially a MPS criterion, except that it considers deviatoric stresses rather than Cauchy stresses. However, since isotropic pressure is arguably small in their case, the MPS and Schmidt-Ishlinsky criteria would be expected to yield similar results. The stress threshold associated with this latter criterion in their study, 158 ± 44 kPa, is consistent with our estimate of 100–130 kPa.

5.4 From crevasse initiation at the surface to full-thickness fracturation

As in most studies aiming to assess a failure criterion and stress threshold for crevasse initiation (Vaughan, 1993; Ultee et al., 2020; Grinsted et al., 2024; Wells-Moran et al., 2025), we focused on the stress fields at the surface, where crevasses are commonly conceptualized to initiate, at least in ablation areas (Colgan et al., 2016). We emphasize that our study addresses crevasse initiation specifically, not fracture propagation, which involves different mechanisms. In particular, our finding that crevasses initiate in the delayed viscous regime rather than as a result of the immediate elastic response following a perturbation does not extend to fracture propagation. Indeed, once a crevasse reaches a critical depth, stress concentrations at the crack tip become the dominant control of the propagation, and elastic processes take over as the primary mechanism, although viscous deformation may still contribute (Hageman et al., 2024). Because fractures are not explicitly represented in the computational domain, continuum-based models cannot directly capture these stress concentrations. Even if fractures were resolved explicitly, the stress distribution near crevasse tips would be highly sensitive to mesh resolution. The inability of continuum models to capture stress singularities caused by the presence of crevasses raises questions about the reliability of the modeled stress field in crevassed regions. In our case, however, no crevasses were reported during the 2010 drainage operations, suggesting that the computed stress fields during this sequence are unlikely to be significantly affected by this limitation.

Studies that have attempted to model calving directly within a viscous framework through CDM, thereby omitting the stress singularities at crack tips, have found that a von Mises criterion produces the most realistic behavior in terms of calving front evolution (Choi et al., 2018; Mercenier et al., 2019). We suggest that this is because only a pressure-independent criterion can allow a water-free fracture to propagate through the full thickness of the ice body when stress concentrations are overlooked. Although this may be an acceptable approach to implement an “averaged” calving rate in large-scale ice sheet models, we argue that it bypasses the fundamental physics of ice fracture and is therefore not universally applicable for forecasting specific ice mass collapse events.

Linear Elastic Fracture Mechanics (LEFM) provides a more physically grounded, albeit still partly empirical, alternative for representing fracture propagation while retaining a continuum framework (van der Veen, 1998). It relies on computing the stress intensity factor, a physical quantity that characterizes the magnitude of the stress concentration at a crevasse tip, based on the far-field stress state, the crevasse geometry, and its initial depth. This stress intensity factor is then compared to a material property, the fracture toughness, which determines whether the fracture will propagate further. Because LEFM requires prescribing an initial crack depth as a starting condition, several authors have combined viscous flow models to describe crack nucleation with LEFM to simulate subsequent fracture propagation (e.g., Krug et al., 2014; Yu et al., 2017; Zarrinderakht et al., 2022). In particular, the hybrid LEFM-CDM strategy used by Krug et al. (2014) for iceberg calving applications is consistent with our findings and represents a meaningful improvement over using CDM alone. However, an even more physically rigorous alternative would be to adopt an energetic formulation rooted in Griffith's theory of brittle fracture, and providing a unified framework for modeling crack initiation, propagation, merging, and branching (Francfort and Marigo, 1998). The phase-field method is precisely such an approach and is therefore particularly promising for glacier fracture modeling. At present, though, its development remains in an early stage, and existing applications are limited to idealized synthetic configurations (e.g., Sun et al., 2021; Clayton et al., 2022; Sondershaus et al., 2023).

6 Conclusions

In this study, we use the rich dataset from the careful monitoring of the artificial drainage of a water-filled cavity on Tête Rousse Glacier to investigate crevasse initiation processes. Specifically, we use numerical experiments to determine the failure criterion and its associated stress threshold that best match the observed circular crevasses that appeared following the 2010 pumping operation. From a Continuum Damage Mechanics perspective, this is equivalent to constrain the source term of the damage evolution equation, i.e. the stress state and magnitude beyond which an undamaged area (D=0) begins to accumulate damage (D>0). In contrast to previous studies that inferred stresses from observed deformations, our good knowledge of the cavity geometry and of the evolution of water level enables a direct computation of stresses from the force balance.

The magnitude and spatial pattern of stresses predicted in the purely elastic regime are insufficient to explain crevasse initiation. Only after the delayed viscous deformation has fed back on viscosity do stress levels compatible with crevasse formation (those computed with the non-linear Glen-Nye law) emerge. This feedback, driven by the non-linear rate-thinning nature of ice, causes stress redistribution and results in higher and more concentrated tensile stresses around the cavity in the delayed viscous regime than in the immediate elastic response, a result that may appear counter-intuitive if one assumes that the initial elastic response always produces the largest stresses. Yet this interpretation is supported by Maxwell-time estimates around the cavity, which are on the order of a day or less: had crevasses opened within this timescale, they would likely have been reported. The fact that none were reported therefore supports the conclusion that crevasses start forming in the delayed viscous regime rather than as a result of the immediate elastic response.

Adopting the usual non-linear Glen-Nye law with n=3, we show that a Maximum Principal Stress criterion, combined to a stress threshold of approximately 100 to 130 kPa, works well for modeling crevasse initiation. The specific stress pattern associated with such a subsidence event combined to the absence of crevasses on the cavity roof allows us to rule out pressure-independent criteria, such as the von Mises criterion, which observation-based studies applied to settings with small isotropic pressure cannot do. Regarding fracture propagation through the full ice thickness, our results support approaches in which crevasse nucleation is described by a viscous-flow model supplemented with a continuum damage mechanics scheme, while subsequent propagation is handled using linear elastic fracture mechanics. Although recent applications of the phase field method to idealized glaciological problems show promise, the method remains in an early stage and requires significant further development.

Our study focuses on the stress fields at the surface, where crevasses are expected to initiate in pure ice settings. In accumulation areas, however, crevasses are expected to form within the firn layer rather than in ice. Firn is characterized by its compressibility and by a viscosity that decreases rapidly and non-linearly with decreasing density (Gagliardini and Meyssonnier, 1997). As a result, unlike in ice, peaks of tensile stress in firn are not located at the surface but rather a few meters below, where the stiffer firn can sustain greater stress than the highly deformable surface layer, and where the overburden pressure is still moderate. We anticipate that the Maximum Principal Stress criterion remains valid in firn, but that the stress threshold must be parameterized as a function of density. This parameterization is currently underway.

Appendix A: Model equations

Here, we present the mathematical model. The computational domain is three-dimensional, with the coordinate system defined as follows: the x axis corresponds to the west-east direction, the y axis to the south-north direction, and the z axis to the vertical, pointing upwards. The variable t represents time.

A1 Field equations and constitutive laws

The distribution of stresses in a body subjected to body forces and external loads is governed by the momentum balance. Neglecting inertial effects, the latter reads

$$\begin{array}{}\text{(A1)}& \text{div}\mathit{\sigma }+{\mathit{\rho }}_{\mathrm{i}}\mathit{g}=\mathrm{0},\end{array}$$

where σ is the Cauchy stress tensor, ρ_i the ice density, and g the gravitational acceleration vector. The Cauchy stress tensor can be decomposed into its deviatoric S and isotropic pI parts as $\mathit{\sigma }=\mathbf{S}-p\mathbf{I}$, with p the isotropic pressure and I the identity matrix. In this work, we adopt the standard convention in fluid mechanics whereby the pressure p is defined as positive in compression and negative in tension. Equation (A1) alone is not sufficient to determine σ, as it consists of three scalar equations for the six independent components of the stress tensor. To close the system, a constitutive law must be prescribed. Below, we detail the linear elastic and viscous constitutive laws that have been used in this study.

A1.1 Linear elastic framework

Assuming deformations are small, the linear elastic behavior of ice is described by Hooke's law, which relates the Cauchy stress tensor σ to the linearized strain tensor ϵ as follows:

$$\begin{array}{}\text{(A2)}& \mathit{\sigma }=\mathit{\lambda }\text{tr}\left(\mathit{ϵ}\right)\mathbf{I}+\mathrm{2}\mathit{\mu }\mathit{ϵ},\end{array}$$

where λ and μ are the first and second Lamé parameters given by:

$$\begin{array}{}\text{(A3)}& \mathit{\lambda }={\displaystyle \frac{E\mathit{\nu }}{(\mathrm{1}+\mathit{\nu })(\mathrm{1}-\mathrm{2}\mathit{\nu })}},\mathit{\mu }={\displaystyle \frac{E}{\mathrm{2}(\mathrm{1}+\mathit{\nu })}},\end{array}$$

where E is Young’s modulus and ν Poisson’s ratio. The components of the linearized strain tensor are expressed as a function of the components of the displacement vector u as follows:

$$\begin{array}{}\text{(A4)}& {\mathit{ϵ}}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right).\end{array}$$

The system of Eqs. (A1)–(A2)–(A4) is now closed and can be solved for the three components of the displacement field. The evolution of the geometry follows from the computed displacement field.

A1.2 Viscous framework

The viscous behavior of ice is described by the Glen-Nye flow law (Glen, 1955), which expresses the deviatoric stress tensor S as a function of the strain rate tensor $\dot{\mathit{ϵ}}$:

$$\begin{array}{}\text{(A5)}& \mathbf{S}=\mathrm{2}\mathit{\eta }\dot{\mathit{ϵ}},\end{array}$$

where the strain rate tensor components are defined as a function of the components of the velocity vector v as:

$$\begin{array}{}\text{(A6)}& {\dot{\mathit{ϵ}}}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial v_i}{\partial x_j}}+{\displaystyle \frac{\partial v_j}{\partial x_i}}\right).\end{array}$$

The effective viscosity η is expressed as:

$$\begin{array}{}\text{(A7)}& \mathit{\eta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(E_{\mathrm{D}}A{)}^{-\mathrm{1}/n}{\dot{\mathit{ϵ}}}_{\mathrm{e}}^{(\mathrm{1}-n)/n},\end{array}$$

where ${\dot{\mathit{ϵ}}}_{\mathrm{e}}=\sqrt{\frac{\mathrm{1}}{\mathrm{2}}{\dot{\mathit{\epsilon }}}_{ij}{\dot{\mathit{\epsilon }}}_{ij}}$ is the second invariant of the strain rate, n the flow law exponent, E_D an enhancement factor, and A the fluidity parameter. The temperature dependence of ice fluidity follows an Arrhenius-type relationship:

$$\begin{array}{}\text{(A8)}& A=A_{\mathrm{0}}\exp \left(-{\displaystyle \frac{Q}{RT}}\right),\end{array}$$

with A₀ a reference fluidity, Q the activation energy, R the gas constant, and T the temperature (in K). In addition, for simulations that are run with the damage module activated, the feedback of damage on viscosity is accounted for via the enhancement factor through the following relationship:

$$\begin{array}{}\text{(A9)}& E_{\mathrm{D}}={\displaystyle \frac{\mathrm{1}}{(\mathrm{1}-D{)}^n}},\end{array}$$

where D is the damage variable. Following Krug et al. (2014), damage evolution is governed by the advection equation:

$$\begin{array}{}\text{(A10)}& {\displaystyle \frac{\partial D}{\partial t}}+\mathit{v}\mathrm{\nabla }D=\left\{\begin{array}{ll}B\mathit{\chi }& \text{if }\mathit{\chi }>\mathrm{0},\\ \mathrm{0}& \text{otherwise}.\end{array}\right.\end{array}$$

Here, B is a damage enhancement factor, and χ is the damage criterion, defined as:

$$\begin{array}{}\text{(A11)}& \mathit{\chi }={\displaystyle \frac{{\mathit{\sigma }}_{\text{eq}}}{\mathrm{1}-D}}-{\mathit{\sigma }}_{\text{th}},\end{array}$$

where σ_th is the stress threshold and σ_eq the equivalent stress. As thoroughly explained in Sect. 3.2, the expression of σ_eq as a function of the stress tensor invariants depends on the selected failure/damage criterion, to which the stress threshold must also be adapted. The formulation in Eq. (A11) implies that ice already experiencing damage is more susceptible to further damage compared to intact ice. This can be interpreted as an empirical approach to account for stress concentration at crevasse tips, which are not explicitly modeled. However, we stress that this implementation alone is insufficient to propagate damage at depth when a maximum principal stress criterion is adopted. A healing criterion could easily be incorporated into this framework (e.g., Pralong and Funk, 2005), but it is not considered here as very little information is available in the literature on the dynamics of healing in ice (Colgan et al., 2016).

The system of Eqs. (A1)–(A5)–(A6) is not closed. This is because, unlike Hooke's law, the Glen-Nye flow law accounts only for the deviatoric part of the Cauchy stress tensor, leaving the isotropic pressure p as an additional unknown. To fully close the system, we incorporate the mass conservation equation. In the viscous framework, ice is assumed to be incompressible, which simplifies the mass conservation equation to div v=0. With this additional constraint, the system is fully determined and can be solved for the three components of the velocity field and the isotropic pressure. The velocity field governs mass transfer, driving the delayed evolution of the geometry toward a new steady state in response to changes in external forcings. Concretely, the computational domain is vertically limited by two free surfaces, namely the upper ice/atmosphere interface $z=z_{\mathrm{s}}(x,y,t)$, and the lower interface $z=z_{\mathrm{b}}(x,y,t)$ that can either be grounded on the bedrock or forms walls of the cavity. An advection equation governs the evolution of these free surfaces over time, which has the general form:

$$\begin{array}{}\text{(A12)}& {\displaystyle \frac{\partial z_j}{\partial t}}+v_x{\displaystyle \frac{\partial z_j}{\partial x}}+v_y{\displaystyle \frac{\partial z_j}{\partial y}}-v_z=a_j,\end{array}$$

where v_x, v_y, v_z are the velocity components at the considered interface, $z_j(x,y,t)$ is the elevation of the considered free surface, and $a_j(x,y,t)$ is the accumulation/ablation rate. Here, the accumulation a_s over the upper interface is prescribed based on field measurements, while melting/accretion a_b in the cavity is neglected. For the lower interface, Eq. (A12) is completed by the non-penetration condition:

$$\begin{array}{}\text{(A13)}& z_{\text{bed}}(x,y,t)\le z_{\mathrm{b}}(x,y,t)<z_{\mathrm{s}}(x,y,t),\end{array}$$

where $z_{\text{bed}}(x,y,t)$ is the elevation of the bedrock. In the vicinity of the cavity, ice cannot penetrate the bedrock but can move away when water pressure becomes sufficiently high. Mathematically, this corresponds to solving a contact problem. At a given point on the lower interface, $\mathit{x}=(x,y,z_{\mathrm{b}})$, ice is considered grounded if it is in contact with the bedrock and the stress exerted by the ice exceeds the water pressure, which is formulated as:

$$\begin{array}{}\text{(A14)}& z_{\mathrm{b}}(x,y,t)=z_{\text{bed}}(x,y,t)\text{ and }{\mathit{\sigma }}_{\text{nn}}{|}_{z_{\mathrm{b}}}>p_{\mathrm{w}}(z_{\mathrm{b}},t),\end{array}$$

where the normal stress at the lower interface is given by ${\mathit{\sigma }}_{\text{nn}}{|}_{z_{\mathrm{b}}}=\mathit{n}\cdot (\mathit{\sigma }\cdot \mathit{n}){|}_{z_{\mathrm{b}}}$, with n denoting the unit normal vector to the interface pointing outward, and p_w representing the water pressure. Conversely, the ice is considered floating if it is either above the bedrock or in contact with it but subjected to a water pressure greater than or equal to ${\mathit{\sigma }}_{\text{nn}}{|}_{z_{\mathrm{b}}}$, which is formulated as:

$$\begin{array}{}\text{(A15)}& z_{\mathrm{b}}(x,y,t)>z_{\text{bed}}(x,y,t)\text{ or }{\mathit{\sigma }}_{\text{nn}}{|}_{z_{\mathrm{b}}}\le p_{\mathrm{w}}(z_{\mathrm{b}},t).\end{array}$$

In Eqs. (A14) and (A15), water pressure is directly deduced from the measured water level z_w (given in meters above sea level) as follows:

$$\begin{array}{}\text{(A16)}& p_{\mathrm{w}}(z_{\mathrm{b}},t)={\mathit{\rho }}_{\mathrm{w}}g\left(z_{\mathrm{w}}\right(t)-z_{\mathrm{b}}(x,y,t\left)\right),\end{array}$$

where ρ_w is the water density.

A2 Boundary conditions

The computational domain is delimited by three boundaries: the upper interface, the lateral boundary, and the lower interface.

A2.1 Upper interface

Neglecting the atmospheric pressure acting on the surface, the upper boundary satisfies a stress-free condition:

$$\begin{array}{}\text{(A17)}& (\mathit{\sigma }\cdot \mathit{n}){|}_{z_{\mathrm{s}}}=\mathrm{0}.\end{array}$$

A2.2 Lateral boundary

On the lateral boundary a no-flux/no-displacement Dirichlet condition applies in the normal direction. This condition is expressed as:

$$\begin{array}{}\text{(A18)}& \mathit{u}\cdot \mathit{n}=\mathrm{0},\end{array}$$

for the elastic framework, and

$$\begin{array}{}\text{(A19)}& \mathit{v}\cdot \mathit{n}=\mathrm{0},\end{array}$$

for the viscous framework. In the tangential direction, a free-slip/free-displacement condition is prescribed. This Neumann condition is expressed as:

$$\begin{array}{}\text{(A20)}& {\mathit{\tau }}_{\mathrm{l}}=\mathrm{0},\end{array}$$

where ${\mathit{\tau }}_{\mathrm{l}}=(\mathit{\sigma }\cdot \mathit{n}){|}_l-{\mathit{\sigma }}_{\text{nn}}{|}_l\mathit{n}$ is the shear stress on the lateral boundary. Note that, because the lateral boundary is very shallow (a few meters at most), the influence of this boundary condition on the simulated displacement and velocity fields is limited to the immediate vicinity of the boundary.

A2.3 Lower interface

Two types of boundary conditions apply at the lower interface, depending on whether the ice is in contact with the bedrock or forms part of the cavity walls. Very little to no sliding is expected at the ice/bedrock interface (Vincent et al., 2015). Consequently, for the elastic framework the no-displacement Dirichlet condition $\mathit{u}(x,y,z_{\mathrm{b}})=\mathrm{0}$ is imposed wherever ice is grounded. For the viscous framework, the numerical implementation of the contact problem in Elmer/Ice requires prescribing a friction law (Durand et al., 2009). We use a linear Weertman law, which relates the basal shear stress ${\mathit{\tau }}_{\mathrm{b}}=(\mathit{\sigma }\cdot \mathit{n}){|}_{z_{\mathrm{b}}}-{\mathit{\sigma }}_{\text{nn}}{|}_{z_{\mathrm{b}}}\mathit{n}$ to the sliding velocity ${\mathit{v}}_{\mathrm{b}}=\mathit{v}-(\mathit{v}\cdot \mathit{n})\mathit{n}$ as:

$$\begin{array}{}\text{(A21)}& {\mathit{\tau }}_{\mathrm{b}}=C_{\mathrm{W}}{\mathit{v}}_{\mathrm{b}},\end{array}$$

where the friction coefficient is set to a high value, $C_{\mathrm{W}}=\mathrm{0.1}\mathrm{MPa}\mathrm{a}{\mathrm{m}}^{-\mathrm{1}}$. As a result, the simulated sliding velocities remain systematically much below 1 m a⁻¹. Additionally, a non-penetration Dirichlet condition is imposed in the normal direction, expressed as:

$$\begin{array}{}\text{(A22)}& \mathit{v}\cdot \mathit{n}=\mathrm{0}.\end{array}$$

For cavity walls, regardless of the chosen framework, water pressure is applied without shear stress, resulting in the following Neumann condition:

$$\begin{array}{}\text{(A23)}& (\mathit{\sigma }\cdot \mathit{n}){|}_{z_{\mathrm{b}}}=-p_{\mathrm{w}}(z_{\mathrm{b}},t)\mathit{n}.\end{array}$$