Failure strength of glacier ice inferred from Greenland crevasses

Ice fractures when subject to stress that exceeds the material failure strength. Previous studies have found that a von Mises failure criterion, which places a bound on the second invariant of the deviatoric stress tensor, is consistent with empirical data. Other studies have suggested that a scaling effect exists, such that larger sample specimens have a substantially lower failure strength, implying that estimating material strength from laboratory-scale experiments may be insufficient for glacier-scale modeling. In this paper, we analyze the stress conditions in crevasse onset regions to better understand the failure criterion and strength relevant for large-scale modeling. The local deviatoric stress is inferred using surface velocities and reanalysis temperatures, and crevasse onset regions are extracted from a remotely sensed crevasse density map. We project the stress state onto the failure plane spanned by Haigh–Westergaard coordinates, showing how failure depends on mode of stress. We find that existing crevasse data are consistent with a Schmidt–Ishlinsky failure criterion that places a bound on the absolute value of the maximal principal deviatoric stress, estimated to be 158±44 kPa. Although the traditional von Mises failure criterion also provides an adequate fit to the data with a von Mises strength of 265±73 kPa, it depends only on stress magnitude and is indifferent to the specific stress state, unlike Schmidt–Ishlinsky failure which has a larger shear failure strength compared to tensile strength. Implications for large-scale ice flow and fracture modeling are discussed.

1 Introduction

Understanding the mechanics of ice fracture is important for predicting the stability of ice sheets and glaciers, since fractures lead to crevassing, calving, and moulins (Colgan et al., 2016): calving is a major component of the mass budget of ice sheets and a potential source of rapid unstable retreat; crevassed and damaged ice is expected to flow more readily (Borstad et al., 2013; MacAyeal et al., 1986; Albrecht and Levermann, 2014; Sun et al., 2017); moulins couple basal conditions to surface conditions (Das et al., 2008) by channeling warmer surface meltwater through the englacial system and potentially accelerating rates of ice loss (cryohydrologic warming; e.g., Phillips et al., 2010; Solgaard et al., 2022). Yet in spite of the broad literature on glacier ice, important material properties that constrain ice fracturing remain understudied.

In classical failure theory, fractures form when an appropriate measure of the internal stress magnitude exceeds a critical value (material failure strength). Fracturing leads to a rapid elastic response with a redistribution of the internal stress, tending to concentrate near material defects and crack tips. Fracture propagation, and therefore ultimately the size of crevasses, is controlled by how the local stress field evolves in the presence of a crack. This kind of propagation has previously been modeled using linear elastic fracture mechanics that relies on material parameters such as the fracture toughness (van der Veen, 1998; Petrovic, 2003; Albrecht and Levermann, 2014).

There is, however, considerable discrepancy between failure strength estimates derived from lab experiments (Hawkes and Mellor, 1972; Haynes, 1978; Petrovic, 2003) and field data (Kehle, 1964; Ambach, 1968; Vaughan, 1993; Petrovic, 2003). Lab experiments found the tensile strength of ice to be 1.7–3.2 MPa for temperatures above −37 °C, whereas estimates inferred from field measurements are an order of magnitude smaller, 90–320 kPa (Vaughan, 1993). While such difference could be partly due to material differences like crystal fabric, grain size distribution, or impurity content, it is more likely the result of sample size affecting failure strength: larger samples have substantially reduced tensile strength compared to smaller samples (Dempsey, 1996; Dempsey et al., 1999), which has been attributed to how the probability of weakest links scale with volume (Petrovic, 2003). Specifically, a factor of 10 increase in the characteristic length scale has been found to reduce the apparent tensile strength by a factor of 3 (Dempsey et al., 1999).

Several failure criteria have been proposed in the literature. Vaughan (1993) used failure maps to justify that both Coulomb and von Mises failure criteria conform well with empirical crevasse data. Lab experiments have meanwhile found that the compressive strength is an order of magnitude greater than the tensile strength (Haynes, 1978; Petrovic, 2003), which is in apparent contradiction with von Mises failure theory. Furthermore, the von Mises criterion is independent of pressure, but Nadreau and Michel (1986) found that failure stresses increase with hydrostatic pressure.

In summary, there is a discrepancy between the failure strength inferred from laboratory and field observations, and multiple distinct failure criteria have been proposed. In this paper, we use a large dataset of Greenland crevasses to estimate a macro-scale failure criterion for naturally occurring glacier ice, relevant for large-scale modeling.

Figure 1 The spatial distribution of crevasses in Greenland from Chudley et al. (2021) is shown in semitransparent lavender blue. The crevasse onset regions along the upstream edge of crevasse fields are shown in magenta. The hatched region in the interior of the ice sheet has been excluded from the analysis.

2 Data

Chudley et al. (2021) derived a complete Greenland Ice Sheet crevasse map by processing high-resolution elevation data from the ArcticDEM v3 mosaic (Porter et al., 2022), which in turn is based on remote sensing data collected over the period from 2007–2015. In this paper, we use the 200 m resolution crevasse density dataset (Fig. 1), which contains the area fraction of 2×2 m pixels that has been classified as crevassed (Chudley, 2022). This optically derived dataset is insensitive to snow-filled crevasses and it therefore likely underestimates the crevasse extent at higher elevations (Chudley et al., 2021).

Strain rates are derived from ice velocity measurements from the Greenland Ice Sheet Velocity Mosaic (Joughin et al., 2016, 2018), which is a multi-mission velocity average spanning the years from 1995 to 2015 (Fig. 2). Surface air temperatures at 2.5 km resolution are taken from the CARRA reanalysis data (Schyberg et al., 2021) averaged over the period from 1991–2020 (Fig. 2).

We compare the calculated von Mises stress to BedMachine v5 ice thickness (H) (Morlighem, 2022), ArcticDEM ice sheet elevation (z) (Porter et al., 2022), along-flow acceleration since 1985 ($\dot{v}$) (Grinsted et al., 2022; Grinsted, 2022), and the seasonal amplitude in ice velocities ($V_{\mathrm{peak}}/V_{\mathrm{winter}}$) derived from PROMICE Sentinel-1 ice velocities (Solgaard and Kusk, 2021; Solgaard et al., 2021).

Figure 2 Average surface air temperatures from CARRA (a), ice velocities from MEAsUREs (b), and near-surface von Mises stress calculated assuming solid ice rheology (c).

3 Methods

Our aim is to relate crevassing to the local stress environment where crevasses initiate. Evidence suggests that most crevasses initiate at a depth of 15–30 m (Colgan et al., 2016), which in a Greenland context can be considered near surface. We will therefore assume that the horizontal velocities at initiation depth are equal to surface velocities. We furthermore assume that vertical shear components are negligible and write the strain rate tensor as

$$\begin{array}{}\text{(1)}& \dot{\mathit{ϵ}}=\left(\begin{array}{ccc}{\dot{\mathit{ϵ}}}_{xx}& {\dot{\mathit{ϵ}}}_{xy}& \mathrm{0}\\ {\dot{\mathit{ϵ}}}_{xy}& {\dot{\mathit{ϵ}}}_{yy}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\dot{\mathit{ϵ}}}_{zz}\end{array}\right),\end{array}$$

where the horizontal components (${\dot{\mathit{ϵ}}}_{xx}$, ${\dot{\mathit{ϵ}}}_{xy}$, ${\dot{\mathit{ϵ}}}_{yy}$) are calculated from the observed surface velocities (Joughin et al., 2016). Crevasses tend to form at lower elevations where temperatures are warmer and firn densities increase more quickly with depth. Theoretical arguments suggest that crevasses initiate below the firn as greater stresses are possible in ice (van der Veen, 1998). Indeed, this is supported by observations of crevasses initiating from refrozen melt layers (Scott et al., 2010; Christoffersen et al., 2018). We therefore assume that fractures initiate in the solid ice and hence that incompressibility applies, ${\dot{\mathit{ϵ}}}_{zz}=-{\dot{\mathit{ϵ}}}_{xx}-{\dot{\mathit{ϵ}}}_{yy}$. Finally, we disregard the effect of pressure on the failure envelope (Nadreau and Michel, 1986) following Vaughan (1993) since only surface crevasses are considered here; that is, crevasse formation subject to relatively small ice pressure.

Glen's flow law for solid ice relates the strain rate tensor to the deviatoric stress tensor (τ_ij) as

$$\begin{array}{}\text{(2)}& {\dot{\mathit{ϵ}}}_{ij}=A\left(T\right){\mathit{\tau }}_{\mathrm{e}}^{n-\mathrm{1}}{\mathit{\tau }}_{ij},\end{array}$$

where A is a temperature-dependent rate factor, n is the flow exponent, ${\mathit{\tau }}_{\mathrm{e}}=\sqrt{I_{\mathrm{2}}}=\sqrt{{\mathit{\tau }}_{ij}{\mathit{\tau }}_{ij}/\mathrm{2}}$ is the effective deviatoric stress, and I₂ is the second invariant of the deviatoric stress tensor. We calculate τ_e from the effective shear strain rate ${\dot{\mathit{ϵ}}}_{\mathrm{e}}$ using ${\mathit{\tau }}_{\mathrm{e}}=[{\dot{\mathit{ϵ}}}_{\mathrm{e}}/A(T){]}^{\mathrm{1}/n}$. Here, we take the canonical value for the flow exponent, n=3, and set the rate factor following Cuffey and Paterson (2010). For simplicity, we assume that the temperature at the crevasse initiation depth can be approximated by the CARRA annual average surface air temperature. With these assumptions, we can calculate the deviatoric stress tensor and the corresponding von Mises stress

$$\begin{array}{}\text{(3)}& {\mathit{\tau }}_{\mathrm{vM}}=\sqrt{\mathrm{3}{I}_{\mathrm{2}}}=\sqrt{{\displaystyle \frac{({\mathit{\tau }}_{\mathrm{1}}-{\mathit{\tau }}_{\mathrm{2}}{)}^{\mathrm{2}}+({\mathit{\tau }}_{\mathrm{2}}-{\mathit{\tau }}_{\mathrm{3}}{)}^{\mathrm{2}}+({\mathit{\tau }}_{\mathrm{3}}-{\mathit{\tau }}_{\mathrm{1}}{)}^{\mathrm{2}}}{\mathrm{2}}}},\end{array}$$

where τ₁, τ₂ and τ₃ are the principal deviatoric stresses.

3.1 The π plane

The stress conditions where failure occurs can be represented by a three-dimensional threshold surface, spanned by the three principal stress components. For pressure-insensitive materials, the threshold surface must be reflection symmetric along the hydrostatic axis (${\mathit{\sigma }}_{\mathrm{1}}={\mathit{\sigma }}_{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{3}}$). This symmetry implies that the threshold surface can, without ambiguity, conveniently be represented on a two-dimensional subspace (plane) by a threshold curve (or “failure envelope”). A popular choice is the so-called π plane view that corresponds to a plane perpendicular to the hydrostatic axis.

In the following, we will determine the failure criteria of ice by plotting principal deviatoric stresses in the π plane and search for a match with hypothesized polygon profiles representing different failure criteria (Kolupaev, 2018). The π plane is the plane spanned by the Haigh–Westergaard coordinates ξ₂ and ξ₃, related to the principal deviatoric stresses through the transformation

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathit{\xi }}_{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}({\mathit{\tau }}_{\mathrm{1}}-{\mathit{\tau }}_{\mathrm{3}})\text{(5)}& {\displaystyle }{\mathit{\xi }}_{\mathrm{3}}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{6}}}}(-{\mathit{\tau }}_{\mathrm{1}}+\mathrm{2}{\mathit{\tau }}_{\mathrm{2}}-{\mathit{\tau }}_{\mathrm{3}}).\end{array}$$

In this plane, the von Mises failure criterion appears as a circle with the radius $r=\sqrt{\mathrm{2}/\mathrm{3}}{\mathit{\tau }}_{\mathrm{vM}}$ since

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{\xi }}_{\mathrm{2}}^{\mathrm{2}}+{\mathit{\xi }}_{\mathrm{3}}^{\mathrm{2}}={\displaystyle \frac{({\mathit{\tau }}_{\mathrm{1}}-{\mathit{\tau }}_{\mathrm{2}}{)}^{\mathrm{2}}+({\mathit{\tau }}_{\mathrm{2}}-{\mathit{\tau }}_{\mathrm{3}}{)}^{\mathrm{2}}+({\mathit{\tau }}_{\mathrm{3}}-{\mathit{\tau }}_{\mathrm{1}}{)}^{\mathrm{2}}}{\mathrm{3}}}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\mathit{\tau }}_{\mathrm{vM}}^{\mathrm{2}},\end{array}$$

where τ_vM is the material failure strength (critical stress).

An alternative visualization (Kolupaev, 2018, Chap. 3.4), common in the glaciological literature, is to focus on the plane where σ₃=0 (i.e., plot σ₁ vs. σ₂) (e.g., Vaughan, 1993). In this plane, the von Mises envelope appears as an ellipse tilted 45°. However, Glen's flow law is concerned with τ_ij, not σ_ij, so it is common to plot τ₁ vs. τ₂ instead, calculated from horizontal velocities (e.g., Chudley et al., 2021). In this case, results must be considered carefully to avoid misinterpretations: the approach implicitly assumes that τ₃=0, which is not valid for horizontally diverging flows. If τ₃=0 is not fulfilled we cannot expect points to fall on a von Mises ellipse even if generated by a von Mises failure process. We therefore advocate for the π plane visualization as it is independent of pressure (i.e., plots for σ or τ gives the same result) and importantly takes all principal stress components into consideration.

3.2 Regions of interest

Crevasses are transported with the flow, and thus not all crevasses are a product of the local stress conditions. As we are interested in the failure strength, we create a mask to select locations where crevasses have recently opened. These are defined as being located on the uphill edge of large-scale crevasse patches. We therefore apply a smoothed Sobel edge detection filter to a binary representation of the crevasse density map and select the upstream edge by requiring that the along-flow crevasse density gradient to increase. We further require the crevasse density to be smaller than 5 %. We label these as “crevassing onset” regions.

The high-elevation ice sheet interior has been excluded from the analysis using a manually traced polygon (Grinsted, 2024), as we gauge that the false positive rate in the crevasse product is greatest there (see Fig. 1). Regions where the ice thickness is below 200 m are also excluded due to concerns of whether the velocity product has sufficient spatial resolution to resolve the local strain rate.

Strain rates and stresses are calculated from long-term average ice velocities. This may not accurately reflect the conditions under which crevasses were formed since trends or seasonal variability in ice flow are unaccounted for (Grinsted et al., 2022; Solgaard et al., 2022) (elaborated on below). We therefore separately examine onset regions with “steady flow”, defined as fulfilling $V_{\mathrm{peak}}/V_{\mathrm{winter}}<\mathrm{2}$ and $\dot{v}<\mathrm{2}$ m a⁻².

4 Results

The spatial distribution of crevasses and the automatic detection of crevasse onset regions are shown in Fig. 1. Maps over the calculated von Mises stress, the ice velocity, and the surface temperature are shown in Fig. 2. Calculated von Mises stresses are found to be greater in onset regions compared to crevassed regions in general and over the ice sheet as a whole (Fig. 3). Onset regions with steady flow are characterized by even greater von Mises stresses, ${\mathit{\tau }}_{\mathrm{vM}}=\mathrm{265}\pm \mathrm{73}$ kPa (Fig. 3). The regional differences in the von Mises stress distribution is shown in the Appendix (Fig. A1). The distribution of principal stresses in onset regions with steady flow are shown in Fig. 4. A similar plot for the distribution of strain rates is shown in Fig. B1. In Fig. A2, we show how the von Mises stress in onset regions relates to other field variables.

Figure 3 The von Mises stress distribution over different subsets of the ice sheet. The counts per bin are all normalized to have the same peak height. Horizontal bars show the 5 %, 50 %, and 95 % percentiles of the distribution.

Figure 4 Empirical failure map showing a π plane density map of stresses in crevasse onset regions with steady flow. The grey color scale is linear in the counts per bin. The empirical median in 5° windows is shown in black. The median von Mises stress is shown as a pink circle. The Schmidt–Ishlinsky median maximum absolute deviatoric stress is shown as a green hexagon. Tensile, compressive, and shear directions have been labeled with T, C, and S, respectively.

5 Discussion

5.1 Uncertainties

Stresses are calculated assuming that ice temperatures, at the depth of crevasse initiation, can be approximated by surface air temperature as simulated by the CARRA regional reanalysis (Fig. 2). It is, however, not clear whether this model is accurate, and there will likely be regional biases in the modeled temperature. Further, ice temperatures at depth may deviate from surface temperatures due to advection and redistribution of energy via melt and refreezing (Løkkegaard et al., 2023; Harrington et al., 2015). We gauge the sensitivity of our stress estimates to temperature by calculating ${\left(A\left(T_{\mathrm{1}}\right)/A\left(T_{\mathrm{2}}\right)\right)}^{\mathrm{1}/n}$ (see Eq. 2), suggesting that a 1 °C change in temperature may result in inferred deviatoric stresses changing by 4 %–8 %. Note that canonical parameters for Glen's flow law were here assumed that disregard the effect of impurities or crystal fabric on ice rheology.

We find that crevasse onset regions are characterized by larger stresses than in crevasse fields in general (Fig. 3). This is not surprising as many crevasses have been transported with the flow after being formed upstream. It is therefore clear that the von Mises stress distribution in crevasse fields does not, in general, accurately reflect the failure stress. We moreover find substantial regional differences in the onset von Mises stress distributions (Fig. A1). Ideally, the estimated failure stress should be independent of region, which suggests that we are not accurately estimating the in situ stress everywhere. We look for clues that might explain regional discrepancies by examining the how the onset region stress co-varies with a range of local conditions (Fig. A2), finding a strong anti-correlation between the von Mises stress and the seasonal amplitude in ice velocities (rank correlation of −0.59). We interpret this as our method systematically underestimating the failure stress in regions with strong seasonal variations in near-surface strain rates, and we therefore argue that onset regions with steady flow more accurately represent the failure stress. Excluding regions with non-steady flow results in a narrower stress distribution that has an improved separation from the stress distribution over ice sheet as a whole (Fig. 3). We therefore estimate the failure stress to be ${\mathit{\tau }}_{\mathrm{vM}}=\mathrm{265}\pm \mathrm{73}$ kPa from onset regions with steady flow. This separation is not as clear in Chudley et al. (2021), which is based on a regional analysis of the same crevasse dataset but does not exclude areas with non-steady flow.

We note that not every location where estimated stresses exceed the failure stress appear to be crevassed (Fig. 2). This may be due to limitations in the crevasse dataset which is insensitive to snow-filled crevasses or because exceeding the failure criterion is a necessary but not sufficient condition for crevasse formation – conditions must also favor crack propagation. Furthermore, in some interior locations, cracks may initiate in the firn pack, thus requiring a failure criterion for firn and a compressible rheology (e.g., Gagliardini and Meyssonnier, 1997; Petrovic, 2003).

5.2 Failure criterion

The empirical failure map in Fig. 4 shows that in steady-onset regions, crevasses predominantly open when subject to large shear or tension and much less so due to compression. This does not necessarily reflect the relative failure strength of the different modes but may simply reflect that the greatest surface stresses occur in shear zones. We therefore estimate the empirical failure envelope as the median of the data in 5° windows in the π plane (solid black line in Fig. 4). The empirical envelope shows that ice in tensile and compressive deformation is nearly equally strong, with the tensile failure strength only being 4 % greater than the compressive strength (see Fig. 4). The empirical envelope compares reasonably well with the von Mises failure criterion, which appears as a circle in the π plane. However, the data clearly have a hexagonal pattern, which we compared against a library of common failure criteria (Kolupaev, 2018, Chap. 3). The hexagonal shape, with corners in the shear directions, strongly suggests that glacier ice follows a Schmidt–Ishlinsky failure criterion (Burzynski, 1928; Schmidt, 1932; Yu, 1983; Kolupaev, 2018).

In the Schmidt–Ishlinsky hypothesis, it is the largest principal deviatoric stress, rather than the magnitude of the second invariant, that defines the failure criterion:

$$\begin{array}{}\text{(7)}& {\mathit{\tau }}_{\mathrm{SI}}=max(\mathrm{|}{\mathit{\tau }}_{\mathrm{1}}\mathrm{|},\mathrm{|}{\mathit{\tau }}_{\mathrm{2}}\mathrm{|},\mathrm{|}{\mathit{\tau }}_{\mathrm{3}}\mathrm{|}),\end{array}$$

where the material failure strength (critical stress) is derived from the inscribed radius r of the hexagon profile in the π plane using ${\mathit{\tau }}_{\mathrm{SI}}=\sqrt{\mathrm{2}/\mathrm{3}}r=\mathrm{158}\pm \mathrm{44}$ kPa. The exscribed radius, corresponding to the shear failure strength, is given by $R=\sqrt{\mathrm{4}/\mathrm{3}}r$, i.e., the Schmidt–Ishlinsky criterion implies that ice is 15 % stronger in shear relative to tensile stresses. Other polycrystalline materials – such as mild steel, copper, nickel alloy, titanium, and stainless steel – have been found to have Schmidt–Ishlinsky-like behavior (see Kolupaev, 2018, Table 2.1).

Alternative criteria have also been proposed in the literature, such as the Coulomb criterion and the maximum strain energy dissipation criterion (MacAyeal et al., 1986; Vaughan, 1993), but are inconsistent with our data as they imply that ice is weakest for shear. We quantify the misfit between previously proposed theoretical failure criteria and the empirical envelope by calculating the root-mean-square log deviation (logRMS). Glen's flow law allows for only deviatoric stress to be estimated, so we are prevented from considering (calculating misfits of) pressure-dependent criteria. The pressure invariant version of the Coulomb criterion is known as the Tresca criterion. We find that the Schmidt–Ishlinsky fits the data best with a logRMS error of 0.02. For comparison, the von Mises and Tresca criteria have logRMS errors of 0.04 and 0.08, respectively. Further, measurement noise acts as to blur the octahedral failure profile, thus softening the corners of the hexagonal profile. We therefore argue that the Schmidt–Ishlinsky criterion is an even better failure model for glacier ice than it appears from the empirical envelope.

In the maximum strain energy dissipation criterion (Vaughan, 1993), it is the rate of deformational work ($P={\dot{\mathit{ϵ}}}_{ij}{\mathit{\tau }}_{ji}/\mathrm{2}$), which is hypothesized to be limited, rather than the stress state that can be withstood by the material. Vaughan (1993) uses Glens flow law (Eq. 2) to calculate τ corresponding to a $\dot{\mathit{ϵ}}$. Unfortunately, estimating the τ failure envelope corresponding to a given threshold value for P becomes ambiguous as there will be a separate failure curve for every temperature. We therefore cannot evaluate the fit using the same logRMS metric we used for the other criteria. However, we calculate the standard deviation of log (P) to be 1.2 in crevasse onset regions with steady flow. This large spread indicates that there is not a single threshold value for P, and we therefore have little confidence in relevance of the maximum strain energy dissipation criterion. For comparison, the log-standard deviations of the estimated threshold values for the von Mises and the Schmidt–Ishlinsky criteria, both have log-standard deviations of 0.3, and thus they are much better fits to the data.

We find that ice thickness correlates positively with the von Mises stress (Fig. A2), contrary to expectations based on the established volume scaling effect. We speculate that the failure strength might approach a limiting value for large ice sample sizes as predicted by, e.g., the multi-fractal scaling law (Carpinteri et al., 1995; Dempsey et al., 1999). This would imply that our reported critical stress values can be used in large-scale ice sheet models without adjustment for the scale effect.

5.3 Modeling crevasse fields

The von Mises and Schmidt–Ishlinsky failure envelopes do not deviate strongly from each other, but, for a given stress magnitude, the failure criterion might be fulfilled in one case and not the other depending on the stress state: the von Mises criterion is indifferent to the stress state, whereas the Schmidt–Ishlinsky criterion favors failure by tension or compression for a given stress magnitude; that is, glacier ice appears stronger when subject to shear as opposed to tensile stresses (see Fig. 4).

The choice of failure criterion may therefore impact where crevasses are formed in large-scale models. A popular approach is to model the evolution of the crevasse density field (or damage phase field) as a function of local conditions and couple it back into ice viscosity (e.g., Albrecht and Levermann, 2012, 2014; Borstad et al., 2016; Sun et al., 2017; Kachuck et al., 2022). In such models, crevasse density evolution is represented as a production–healing process, where production depends on a scalar measure of the local stress state, which, if exceeding some critical value, is taken to grow proportionally to the local principal spreading rate (strain rate). Our results are therefore directly applicable to such modeling efforts: we provide both the failure criterion (Eq. 7) and critical value ${\mathit{\tau }}_{\mathrm{SI}}=\mathrm{158}\pm \mathrm{44}$ kPa. We leave it for future work to test whether modeled crevasse density fields using the Schmidt–Ishlinsky criterion in fact provide the best fit with observations.

5.4 Third stress tensor invariant

The Schmidt–Ishlinsky criterion implies that the failure of ice depends not only on the second, but also the third invariant of the deviatoric stress, $I_{\mathrm{3}}={\mathit{\tau }}_{ij}{\mathit{\tau }}_{jk}{\mathit{\tau }}_{ki}/\mathrm{3}$, since Eq. (7) can equivalently be written as (Yu, 1983)

$$\begin{array}{}\text{(8)}& {\displaystyle }\left(I_{\mathrm{3}}+I_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{SI}}-{\mathit{\tau }}_{\mathrm{SI}}^{\mathrm{2}}\right)\left(I_{\mathrm{3}}-I_{\mathrm{2}}{\mathit{\tau }}_{\mathrm{SI}}+{\mathit{\tau }}_{\mathrm{SI}}^{\mathrm{2}}\right)=\mathrm{0}.\end{array}$$

The plastic behavior and failure of materials are often linked, which suggests that the third invariant should be included in the flow law for ice; a dependency that has previously been proposed (Glen, 1958; Van Der Veen and Whillans, 1990; Morland and Staroszczyk, 2019; Baker, 1987). Indeed, Steinemann (1954) found that ice deforms slower in shear than expected from uniaxial deformation tests (Glen, 1958) and suggested that deformation data are not consistent with Glen's canonical co-axial flow law with a dependence on only I₂ but should also depend on I₃. Although this seems to be consistent with the Schmidt–Ishlinsky failure pattern found here (Fig. 4), constraining the stress response functions in a flow law depending on both I₂ and I₃ has proven difficult from available data (Morland and Staroszczyk, 2019; Staroszczyk and Morland, 2022). On that note, Staroszczyk and Morland (2022) recently showed that a quadratic (non-coaxial) flow law depending only on I₂ can also be constructed to improve the fit with deformation experiments. Nonetheless, we believe it is worth considering whether the failure criterion might also inform on the form of the flow law in addition to correlating with deformation data.

6 Conclusions

We automatically identified crevasse onset regions from a dataset of Greenland crevasses (Fig. 1), where we argue local stress conditions must have exceeded the failure stress. We inferred the local stress conditions using Glen's flow law combined with observed ice velocities and average surface air temperatures (Fig. 2), disregarding regions with seasonally variable ice flow. We estimated the failure strength of glacier ice to be ${\mathit{\tau }}_{\mathrm{vM}}=\mathrm{265}\pm \mathrm{73}$ kPa in crevasse onset regions with steady flow (Fig. 3). This is compatible with the 90–320 kPa tensile strength estimated by Vaughan (1993). The corresponding π`plane failure map (Fig. 4) suggests, however, that the mechanical failure of glacier ice is best modeled using the Schmidt–Ishlinsky failure criterion, where failure occurs once the maximum absolute deviatoric stress exceeds 158±44 kPa. While we argue that the a Schmidt–Ishlinsky criterion is the better model for ice failure, we note that the deviations to the von Mises criterion are quite small (Fig. 4). It may therefore still be a good approximation to model ice failure with a von Mises criterion on average. This, however, disregards the effect that the mode of deformation has on the failure strength, where the Schmidt–Ishlinsky criterion implies that ice is 15 % stronger in shear relative to tensile stresses.

This study is, to our knowledge, the first to propose a Schmidt–Ishlinsky failure criterion for glacier ice. More work is needed to validate this hypothesis using, e.g., forward modeling.

Appendix A: Regional von Mises stress

The failure strength of ice should largely be independent of location. We therefore plot the onset von Mises stress distribution by region (Zwally et al., 2012). We find that there is substantial variability between regions (Fig. A1). To illuminate the source of the scatter we examine how the onset von Mises stress varies with other fields in a pair plot (Fig. A2). We transform variables that have long-tailed distributions with either a logarithm or hyperbolic tangent function prior to visualization. We find that the von Mises stress anti-correlates with seasonal velocity amplitudes. This suggests that the von Mises stress, calculated from long-term average velocities, is biased low in some regions as it does not capture the seasonal peak stress which may drive crevasse formation.

Figure A1 The calculated von Mises stress for opening crevasses varies substantially between different basins of the ice sheet. The counts per bin are all normalized to have the same peak height. Horizontal bars show the 5 %–50 %–95 % percentiles of the distribution.

Figure A2 Pair plot showing the relationships between various quantities extracted over regions where crevasses are opening. The variables plotted are the von Mises stress (τ_vM), the seasonal peak velocity relative to the winter velocity ($V_{\mathrm{peak}}/V_{\mathrm{winter}}$), the mean velocity from MEAsUREs (v), the long-term acceleration in the along-flow direction (A), the surface temperature (T_surf), and the ice thickness (Thick.). Some variables with long-tailed distributions have been transformed using logarithms or hyperbolic tangent.

Appendix B: Critical strain rates

In this paper, we focus on the stress conditions that result in failure of glacier ice. We do not observe the stress state, but infer it from observed strain rates and our understanding of the mechanical properties of ice. In this section, we present the raw strain rate data. To investigate how the critical strain rate depends on the mode of deformation we define strain rate π plane coordinates in an equivalent way to Eq. (4) but using the three invariants of ${\dot{\mathit{ϵ}}}_{ij}$ rather than τ_ij. Figure B1 shows the observed strain rates in crevasse onset regions with steady flow. Vaughan (1993) found that there is a strong impact on temperature on the critical strain rate, which can largely be accounted for by the rate factor in Glen's flow law (Eq. 2). We therefore calculate empirical failure envelopes in three different temperature intervals (Fig. B1) and find that the critical strain rate varies by an order of magnitude between them. This is reflected in a long-tailed strain rate distribution. All three empirical failure envelopes show a similar flower-like pattern where the critical strain rate is greater for shear than for tensile and compressive deformation.

Figure B1 Empirical failure map showing a π plane density map of strain rates in crevasse onset regions with steady flow. The π plane coordinates are calculated from strain rates rather than deviatoric stresses. The color scale is linear in the counts per bin. The empirical median in 10° windows for different temperature intervals is shown in thick lines. Tensile, compressive, and shear directions have been labeled with T, C, and S, respectively.