Failure strength of glacier ice inferred from Greenland crevasses

large-scale


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., 2017b); moulins couple basal conditions to surface conditions (Das et al., 2008) by channeling warmer surface melt water 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 under-studied.
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 1 https://doi.org/10.5194/egusphere-2023-1957Preprint.Discussion started: 19 September 2023 c Author(s) 2023.CC BY 4.0 License. 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 modelled 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 MPa to 3.2 MPa for temperatures above −37 • C, whereas estimates inferred from field measurements are an order of magnitude smaller, 90 kPa to 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 10 increase in the characteristic length scale has been found to reduce the apparent tensile strength by a factor three (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 to von Mises failure theory.Furthermore, both von Mises and Coulomb failure criteria are 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 data set of Greenland crevasses to estimate a macroscale failure criterion for naturally occurring glacier ice, relevant for large-scale modelling.

Data
Chudley et al. ( 2021) derived a Greenland Ice Sheet wide 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 200m resolution crevasse density data set (fig. 1) which contains the area fraction of 2x2m pixels that has been classified as crevassed (Chudley, 2022).Strain rates are derived from ice velocity measurements from the Greenland Ice Sheet Velocity Mosaic (Joughin et al., 2016(Joughin et al., , 2018) ) which is a multi-mission velocity average spanning the years from 1995 to 2015 (fig.2).Surface air temperatures at 2.5km resolution are taken from the CARRA reanalysis data (Schyberg et al., 2021) averaged over the period from 1991-2020 (fig.2).
We compare our results to BedMachine v5 ice thickness (H) (Morlighem, 2022), ArcticDEM ice sheet elevation (z) (Porter et al., 2022), along flow acceleration since 1985 ( v) (Grinsted et al., 2022;Grinsted, 2022), and the seasonal amplitude in ice velocities (V peak /V winter ) derived from PROMICE Sentinel-1 ice velocities (Solgaard and Kusk, 2021;Solgaard et al., 2021).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-30m (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 neglible and write the strain rate tensor as where the horizontal components ( εxx , εxy , ε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, εzz = − εxx − ε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, crevasses subject to relatively small pressure.
Glen's flow law for solid ice relates the strain-rate tensor to the deviatoric stress tensor (τ ij ) as where A is a temperature-dependent rate factor, n is the flow exponent, τ e = √ I 2 = τ ij τ ij /2 is the effective deviatoric stress, and I 2 is the second invariant of the deviatoric stress tensor.We calculate τ e from the effective shear strain rate εe using τ e = A(T ) εn e .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 where τ 1 , τ 2 and τ 3 are the principal deviatoric stresses.

π plane
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 ξ 2 and ξ 3 , related to the principle deviatoric stresses through the transformation In this plane, the von Mises failure criterion appears as a circle with radius r = 2/3τ vM since where τ vM is the material failure strength (critical stress).

Regions of interest
Crevasses are transported with the flow, and so 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 as we gauge that the false positive rate in the crevasse product is greatest there.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 peak /V winter < 2 and v < 2 m a −2 .

Results
The

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 modelled 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., 2022).We gauge the sensitivity of our stress estimates to temperature by calculating (A(T 1 )/A(T 2 )) 1/n (see eqn. 2), suggesting that a 1 • C change in temperature may result in inferred deviatoric stresses changing by 4 % to 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, and so this 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 therefore argue that onset regions with steady flow more accurately represent the failure stress, estimated to be τ vM = (265 ± 73) kPa (fig.3).
We note that not every location where estimated stresses exceed the failure stress appear to be crevassed (fig.2).Exceeding the failure criterion is a necessary but not sufficient condition for crevasse formation, conditions must also favour crack propagation.And 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).

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 compression.The data compares reasonably well with the von Mises failure criterion which appears as a circle in the π-plane.However, the data clearly has a hexagonal pattern, suggesting 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: where the material failure strength (critical stress) is derived from the inscribed radius r of the hexagon profile in the πplane using τ SI = 2/3 r = (158 ± 44) kPa.The exscribed radius, corresponding to the shear failure strength is given by R = 4/3 r, i.e. the Schmidt-Ishlinsky criterion implies that ice is 15% stronger in shear relative to tensile stresses.
The estimated threshold values for the von Mises and the Schmidt-Ishlinsky criteria have nearly identical log-standard deviations (0.3) and thus both are good fits to the data.However, 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 a better failure model for glacier ice than it appears from the profile.
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 (not shown).
We find that ice thickness correlates 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.

Modelling 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 andLevermann, 2012, 2014;Borstad et al., 2016;Sun et al., 2017a;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 modelling efforts: we provide both the failure criterion (7) and critical value τ SI = (158 ± 44) kPa.We leave it for future work to test whether modelled crevasse density fields using the Schmidt-Ishlinsky criterion in fact provide the best fit with observations.

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 3 = τ ij τ jk τ ki /3, since eqn.(7) can equivalently be written as (Yu, 1983) The plastic behaviour 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;Veen and Whillans, 1990 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 2 , but should depend on I 3 , too.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 2 and I 3 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 2 can also be constructed to improve the fit with deformation experiments.Nonetheless, we believe it is worth considering whether the failure criterion might, too, inform on the form of the flow law in addition to correlating with deformation data.

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 tensile strength of glacier ice to be τ vM = (265 ± 73) kPa in crevasse onset regions with steady flow (fig.
3).This is compatible with the 90 kPa to 320 kPa estimated by Vaughan (1993).The corresponding π-plane failure map (fig. 4) suggests, however, that the mechanical failure of glacier ice is best modelled 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 modelling.

Figure 1 .
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.

Figure 2 .
Figure 2. Average surface air temperatures from CARRA (left), Ice velocities from MEAsUREs (middle), and near surface von Mises stress calculated assuming solid ice rheology (right).
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, τ vM = (265 ± 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 figure 4. In fig.A2, we show how the von Mises stress in onset regions relates to other field variables.https://doi.org/10.5194/egusphere-2023-1957Preprint.Discussion started: 19 September 2023 c Author(s) 2023.CC BY 4.0 License.

Figure 3 .
Figure 3.The von Mises stress distribution over different subsets of the ice sheet.Bars show the 5%, 50% and 95% percentiles of the distribution.

Figure 4 .
Figure 4. Empirical failure map showing a π-plane density map of stresses in crevasse onset regions with steady flow.The empirical median in 5 • windows is shown in black; The median von Mises stress is shown as a pink circle; And the Schmidt-Ishlinsky median maximum absolute deviatoric stress is shown as a green hexagon.Tensile, Compressive and Shear directions have been labelled with T, C, and S.

Figure A1 .
Figure A1.The calculated von Mises stress for opening crevasses varies substantially between different basins of the ice sheet.Bars show the 5%-50%-95% percentiles of the distribution.

Figure A2 .
Figure A2.Pair plot showing the relationships between various quantities extracted over regions where crevasses are opening.