Tensile strength of glacial ice deduced from observations of the 2015 eastern Skaftá cauldron collapse, Vatnajökull ice cap, Iceland

Abstract

The representation of iceberg calving in numerical models is a key source of uncertainty in century-scale sea-level rise projections. Parameters central to model representations of calving, including the tensile strength of glacier ice, remain poorly constrained. Grain-size and sample-size dependence make it difficult to reconcile laboratory and in situ estimates of ice tensile strength. Further, assumptions of various numerical models obscure comparison of the ‘strength’ parameter with a physically observable value. Here, we address the problem of fracture during calving using an analogous natural laboratory: a viscoelastic analysis of observed surface deformation and associated stresses in the 2015 collapse of eastern Skaftá cauldron, Vatnajökull ice cap, Iceland. We find that the ice within the cauldron could have experienced instantaneous elastic stress on the order of several MPa. We observe surface crevasses at the cauldron edges and center, but find that large areas of ice remain intact despite high stress. Our findings suggest a tensile strength of glacier ice on the order of 1 MPa, consistent with laboratory estimates but exceeding previous glacier-specific estimates.

1. Introduction

Iceberg calving is a key source of uncertainty in future projections of global sea-level rise (Church and others, 2013). Recent modeling efforts, including the use of continuum damage mechanics and discrete element models, have made considerable progress toward a dynamically consistent representation of calving. However, many models still rely on poorly constrained parameters. One parameter appearing in several studies is the strength of glacier ice (e.g. Pralong and others, 2003; Duddu and Waisman, 2013; Åström and others, 2013; Krug and others, 2014; Benn and others, 2017). Here, we refer to ice strength as a bulk property equaling the maximum stress intact ice can sustain before fractures appear at macroscale, and we focus specifically on the tensile stress regime. Although there are some constraints on the tensile strength of glacier ice from laboratory experiments (Currier and Schulson, 1982; Lee and Schulson, 1988; Xian and others, 1989; Druez and others, 1989) and field observations (Vaughan, 1993), the ranges derived from the different methods barely overlap. Further, models disagree about what part of the observed range is relevant to glacier modeling and even what style of deformation – viscous, elastic or a combination – determines the stresses associated with fracture (e.g. Bassis and Ma, 2015; Ma and others, 2017; Borstad and others, 2017).

A major challenge in constraining the strength of glacier ice is that estimates derived from laboratory experiments do not agree with the limited estimates from in situ observations. Laboratory estimates range from 0.7 to 3.1 MPa (review by Petrovic, 2003) for laboratory and lake ice; Druez and others (1989) found tensile strength as high as 5 MPa for laboratory-grown glaze ice. In situ observations of glaciers, meanwhile, suggest a tensile strength ranging from 0.09 to 0.32 MPa for glacier ice (Vaughan, 1993). Part of the discrepancy may be the grain-size dependence of the tensile strength of ice (Currier and Schulson, 1982). Laboratory experiments conducted with finer-grained ice than is typical in glaciers may overestimate the tensile strength. In addition, failure can be detected at millimeter scale in the laboratory (Currier and Schulson, 1982) but only at multi-meter scale in existing field and remote-sensing observations (Colgan and others, 2016). Another source of discrepancy is that stress is difficult to observe in situ (Pfeffer and others, 2000; Colgan and others, 2016), so estimates derived from glacier-scale observations, such as Vaughan (1993), must assume a rheology to convert observable strain rates into stresses. Assuming a viscous, Glen's law-type rheology will produce quantitatively different stresses, and therefore different estimates of ice strength, than would an elastic or viscoelastic rheology (Reeh and others, 2003; Borstad and others, 2017). Further, Dempsey and others (1999) found that the measured tensile strength of sea ice scaled with sample size, implying that tests across a wide range of sample scales are needed to relate fracture mechanics observed in the laboratory with field-scale processes.

A second challenge in determining a consistent tensile strength of glacier ice is that various numerical models select different ice strengths according to the demands of their fracture- or damage-modeling framework. For example, Åström and others (2013, 2014) and Benn and others (2017) implement an elastic stress threshold to break ‘beams’ connecting elements of ice. The stress threshold σ_c ranges from 0.1 to 1 MPa, within the range of both laboratory and field estimates. Yet the numerical model used in those studies, the Helsinki Discrete Element Model, uses a Young's modulus an order of magnitude too low for glacier ice to offset the effect of model elements several orders of magnitude larger than glacier ice grains (Benn and others, 2017). Because elastic stress is dependent on Young's modulus (see Eqns (1–3) below), and the strength of glacier ice is grain-size dependent (Currier and Schulson, 1982), it is unclear that the stress threshold σ_c is the same physical quantity measurable in the laboratory or field.

By contrast, Krug and others (2014) use a tensile stress threshold for fracture of 0.20 MPa in their continuum damage mechanics framework. Pralong and others (2003), Pralong and Funk (2005) and Duddu and Waisman (2013) use a similar threshold, with tensile strength ranging from 0.20 to 0.50 MPa, based on the macroscale condition derived by Vaughan (1993). This entire range is below even the lowest laboratory estimate of ice tensile strength (Petrovic, 2003). The authors write that the stresses their model produces are too small to reach any higher damage threshold. However, the stresses they model with a Glen's law rheology are active over timescales much longer than the viscoelastic relaxation time of ice (Maxwell time t _r ≈ 8–12 h, see Table 1), and thus do not capture instantaneous, elastic fluctuations of stress that could be more important in exceeding the stress threshold for crevasse propagation. Indeed, Banwell and others (2019) estimated maximum elastic stress due to the flexure of an Antarctic ice shelf to reach 0.5 MPa without any associated fracturing.

Table 1. Material parameters and settings used in the analysis of Sections 3.1 and 3.2

*min/max corresponding to ranges of η, ν, E.

Ice subsidence events offer a novel avenue to examining ice strength at geophysical scale. A number of theoretical results have taken advantage of the defined loading at ice cauldrons and supraglacial lake sites. For example, Banwell and others (2013) used the elastic-plate analysis derived by Sergienko (2005) to compute the fracture spacing associated with supraglacial lake drainage on Antarctic Peninsula ice shelves. Evatt and Fowler (2007) computed fracture spacing associated with cauldron collapse on the basis of viscous beam theory derived in Evatt and others (2006). Even extraterrestrial ice subsidence has supported studies of ice stress and fracture: Walker and Schmidt (2015) computed the stress and surface morphology associated with collapse of ice over trapped water pockets on icy satellites, using a model of ice shell flexure with an elastic fracturing layer (Walker and others, 2012a).

Here, we build on previous studies of subsidence events to estimate the surface stress associated with a single ice cauldron collapse event. We apply a new viscoelastic analysis to interpret detailed observational data of the 2015 collapse of eastern Skaftá cauldron, Vatnajökull ice cap, Iceland. Our Maxwell viscoelastic rheology accounts for both short-term elastic and longer-term viscous deformation (for viscoelastic treatment of glacier ice see Gudmundsson, 2011; Goldberg and others, 2014; MacAyeal and others, 2015; Robel and others, 2017), but we focus on the elastic stress active over the short timescale of crevasse nucleation and propagation. Our analysis is well constrained by observations, including high-resolution, time-dependent digital elevation models of the cauldron surface to constrain the magnitude of collapse (Porter and others, 2018) and an in situ GPS record that constrains the timescale of collapse. We are thus able to estimate the maximum stress that the cauldron collapse could have produced, and we compare the maximum stress field with the observed crevasse locations to constrain glacier tensile strength.

2. Physical setting of eastern Skaftá cauldron

The Eastern and Western Skaftá ice cauldrons are two spots of elevated geothermal heat flux located in the southwest of Vatnajökull ice cap, Iceland (Fig. 1). Locally warm conditions at the base of 400 m thick ice lead to enhanced subglacial melting and eventual flotation, creating a 300 m thick ‘internal ice shelf’ confined on all sides by grounded, temperate ice. Water builds up at the base of the cauldrons for 2–5 years (Guđmundsson and others, 2016) and finally is drained by glacial outburst floods (jökulhlaups) lasting hours to days (Björnsson, 1992). Sudden drainage of the cauldrons leads to high strains and stresses in the cauldron ice, producing rings of fractures. In quiescent years, the cauldrons can be identified by persistent, kilometer-wide depressions on the glacier surface (Einarsson and others, 2016).

Fig. 1. Location of eastern Skaftá cauldron (blue triangle) on the Vatnajökull ice cap and in Iceland (inset). In the inset map of Iceland, dark stripes indicate volcanic regions, white patches indicate glaciers and ice caps, and the red rectangle indicates the region of interest. Surface topography is from ETOPO1 (Amante and Eakins, 2009), with contours at intervals of 250 m a.s.l (light) and 1000 m a.s.l. (heavy contours). Light blue sinuous lines are rivers, white is ice cover, and brown and green are non-ice surface. The figure was made using the Generic Mapping Toolbox (GMT; Wessel and others, 2013).

Between 29 September and 3 October 2015, subglacial water that had accumulated over 5 years drained from the eastern Skaftá cauldron in a jökulhlaup of record proportion: peak discharge downstream in the river Skaftá exceeded 3000 m³ s⁻¹ (Jóhannesson and others, 2016). The cauldron collapse created a roughly circular surface depression, approximately 110 m deep at its center and 2.7 km in diameter. Rings of fractures are visible at the ice surface in optical imagery and in the 2 m resolution Arctic Digital Elevation Model (‘ArcticDEM’; Porter and others, 2018) snapshot acquired on 10 October 2015, just 7 days after the end of the drainage event. Figures 2 and 3 show satellite, aerial and cross-sectional views of the cauldron and its crevasses.

Fig. 2. Comparison of the surface of eastern Skaftá cauldron in (a) October 2012 and (b) October 2015 based on hillshade view of ArcticDEM surface elevation data (Porter and others, 2018); and (c) an oblique aerial view (photo by Ragnar Axelsson, used with permission) of the cauldron following its collapse. In panel (b), a red triangle indicates the location of a GPS station maintained by the Icelandic Meteorological Office, a blue line indicates the transect shown in Fig. 3, and an area with no data appears white. Horizontal scale shown in panel (a) is maintained in panel (b).

Fig. 3. (a) Surface elevation on 15 October 2012 (labeled ‘Pre-collapse’) and 10 October 2015 (labeled ‘Post-collapse’) from ArcticDEM (Porter and others, 2018) along transect P − P′ shown in Figure 2. Red triangle indicates GPS station location, with a representative vertical path during subsidence indicated by dotted line. (b) GPS record of net subsidence from initial elevation of 1660 m; (c) GPS vertical displacement rate during the 2015 collapse. Horizontal axis labels on lower panels indicate 2015 calendar date.

A GPS station, placed near the center of the cauldron by the Icelandic Meteorological Office, recorded vertical displacement during the collapse with temporal resolution of 5 s (Guđmundsson and others, 2018). Figure 3 shows the net subsidence and subsidence rate recorded between 27 September and 10 October 2015. We note two distinct phases of collapse: an initial rapid collapse, with peak subsidence rates of 3 m h⁻¹, followed by a prolonged period of slower (0.5 m h⁻¹) settling.

3. Maximum stresses in intact and fractured ice

We will use two complementary methods to compute the maximum instantaneous stress associated with the cauldron collapse. Motivated by the temporal pattern of subsidence described above, we treat glacier ice as a viscoelastic material (e.g. Reeh and others, 2003; Gudmundsson, 2011; Rosier and others, 2014). The maximum instantaneous stress in a viscoelastic cauldron collapse is elastic in nature and sustained during an initial period shorter than the Maxwell time t _r. As the system approaches time t _r, viscous dissipative effects become dominant. In both regimes, stress is highest at the surface and base of the ice, a distance of half the ice thickness h from a central neutral plane of deformation (e.g. Ugural, 2018). In Section 3.1, we ignore time-dependence to estimate the largest-magnitude linear elastic stress that could have arisen from the observed surface curvature. The GPS record (Fig. 3) suggests, however, that there was a non-negligible period of slow subsidence lasting several days, consistent with viscous settling. In Section 3.2, we account for these viscous effects using a Maxwell linear viscoelastic rheology applied to an idealized circular cauldron.

3.1. Purely elastic stress estimate

We first investigate the maximum stress possible under a linear elastic collapse. Stress concentration around any pre-existing fractures would tend to limit stress near the surface of a collapsing ice plate; here, we examine the large-magnitude stresses that could be produced if the entire surface of the eastern Skaftá cauldron were intact prior to collapse. Assuming initially intact ice also facilitates an estimate of surface stress directly from the observed surface slope.

In the purely elastic regime, the normal stresses σ_xx, σ_yy and shear stress τ_xy at the surface of a sagging two-dimensional ice plate are related to the surface curvature (Ugural, 2018):

(1)$$\sigma_{xx} = -{Eh\over 2 \lpar 1-\nu^2\rpar } \left(\kappa_{xx} + \nu \kappa_{yy}\right)\comma\; $$

(2)$$\sigma_{yy} = -{Eh\over 2 \lpar 1-\nu^2\rpar } \left(\kappa_{yy} + \nu \kappa_{xx}\right)\comma\; $$

(3)$$\tau_{xy} = -{Eh\over 2\lpar 1 + \nu\rpar } \kappa_{xy}\comma\; $$

where h is the ice thickness, ν is Poisson's ratio, E is Young's modulus and curvature κ _ij is the second derivative with respect to spatial coordinates i, j of surface elevation S. Representative values of h, ν, E are given in Table 1. We note that our analysis uses parameter values calibrated to reflect the magnitude and speed of observed cauldron deformation, which should not be interpreted to constrain their true material value. For example, our effective Young's modulus E = 1.0 GPa produces appropriate subsidence but is at the low end of estimated material Young's modulus of ice (e.g. 0.8 GPa in Vaughan, 1995; 4–10 GPa in Rist and others, 2002).

From the stress components of Eqns (1–3), we compute the maximum principal stress

(4)$$\sigma_1 = {\sigma_{xx} + \sigma_{yy}\over 2} + \sqrt{\left({\sigma_{xx} - \sigma_{yy}\over 2} \right)^2 + \tau_{xy} ^2}\comma\; $$

with the sign convention that σ₁ > 0 is compression and σ₁ < 0 is tension.

We used Eqns (1–4) to deduce the largest-magnitude elastic stress that the 2015 eastern Skaftá cauldron collapse could have generated. First, we approximated an ‘intact’ ice surface from the 10 October 2015 ArcticDEM data. We applied a median filter with a 10 m × 10 m window to the 2 m ArcticDEM surface (Porter and others, 2018) and used a high-pass filter with 1 m threshold to identify and mask crevasses. The resulting mask is shown in Figure 4a. We then fit a two-dimensional, 5th-order B-spline to the filtered and masked surface elevation (Fig. 4b) using built-in functions of SciPy v1.2.1 and calculated the surface curvatures $\kappa _{ij} = \partial _{ij}^2 S$. Next, we calculated the stress components with Eqns (1–3) and maximum principal stress with Eqn (4) (Fig. 4c). Finally, we examined the maximum principal stress in the areas we had identified as crevassed (or ‘fractured’) and un-crevassed (or ‘intact’).

Fig. 4. (a) Mask distinguishing intact ice (dark gray, 3 964 764 of 4 443 505 pixels or 89% of the surface) from unmasked fractured ice (402 396 pixels or 9.1% of the surface); (b) smooth interpolated post-collapse surface elevation; and (c) corresponding maximum principal stress field for eastern Skaftá cauldron. All images include hillshading from ArcticDEM to reveal surface crevasses, and hatching indicates no-data areas in the ArcticDEM observations (1.7% of the surface). Ticks on the outside of each panel appear at 500 m intervals.

3.2. Distribution of elastic stresses in intact versus fractured areas

Figure 4 shows the post-collapse surface and pattern of maximum principal stress we computed. There are alternating, roughly concentric areas of high tensile ( − ) and high compressive (+) stress. The edge of the cauldron generally shows high tensile stress, though stress is lower and ice remains intact in two regions where the radial symmetry is distorted (Figs 4a, c). There is another area of high stresses, both compressive and tensile, near the center of the cauldron, where a bump creates steeper surface curvature. Guđmundsson and others (2016) suggest that the bump is an area of thicker ice rather than a bedrock protrusion. We observe crevasses in the high-stress areas at both the edge and the center of the cauldron. In between, the ice surface appears intact.

Figure 5 summarizes the distribution of stresses in areas coinciding with intact (gray) or fractured (white) ice. We sampled the maximum principal stress field shown in Figure 5 for all points in the domain on the 2 m × 2 m ArcticDEM grid. In areas of intact ice, which accounts for 89% of the ice surface, the distribution of surface maximum principal stress peaks near 0 MPa. However, the intact ice area also includes locations of higher stress; more than 20% of the intact ice sampled is found where maximum principal stress exceeds 1 MPa in tension. By contrast, the surface maximum principal stress distribution for crevassed areas (9.1% of the ice surface) is flatter, with greatest frequency around 5 MPa tension. The stress distributions for intact and crevassed areas are distinct: Less than 3% of the fractured sample had maximum principal stress 0 − 500 kPa in tension, and tensile stresses up to 2.5 MPa are more common in intact than in crevassed areas.

Fig. 5. Normalized histogram of maximum elastic stresses within the cauldron, at locations identified as intact (dark gray) or fractured (white) from the ArcticDEM surface observations (Porter and others, 2018). Red shading denotes the tensile regime and blue shading the compressive regime. Vertical dashed line indicates 1 MPa tensile stress, which we suggest as the tensile strength of glacier ice in Section 4.

3.3. Stress associated with viscoelastic collapse

Because the 2015 eastern Skaftá cauldron collapse took place over several days, viscous deformation likely played a role in producing the observed post-collapse surface. As a result, deducing stress from final observed surface curvature as in Section 3.1 will tend to overestimate the maximum principal stresses that could have been active during the elastic phase of collapse. We now introduce a Maxwell viscoelastic rheology to account for both viscous and elastic effects. Several previous authors have applied Maxwell viscoelasticity to glacial ice under transient loading (Gudmundsson, 2011; Goldberg and others, 2014; MacAyeal and others, 2015; Banwell and MacAyeal, 2015; Robel and others, 2017); in particular, Gudmundsson (2011) shows that Maxwell viscoelasticity is an appropriate simplification from Burgers viscoelasticity as implemented by Reeh and others (2003).

A Maxwell viscoelastic material combines viscous and elastic elements in series. At short timescales (t < t _r), the deformational response to forcing is elastic, while at longer timescales (t > t _r) viscous deformation dominates. Following Howell and others (2009), the Maxwell constitutive relation for linear viscoelasticity is

(5)$${\eta\over \mu} {\partial \sigma_{ij}\over \partial t} + \sigma_{ij} = {\eta\over \mu} {\partial\over \partial t} \left(\lambda \delta_{ij} \varepsilon_{kk} + 2\mu \varepsilon_{ij} \right)\comma\; $$

with η the dynamic viscosity, μ the shear modulus, σ_ij the Cauchy stress tensor, ɛ _ij the strain tensor, λ the first Lamé parameter and δ_ij the Kronecker delta. The elastic moduli of Eqn (5),

(6)$$\mu = {E\over 2 \lpar 1 + \nu\rpar }\comma\; \qquad \lambda = {\nu E\over \lpar 1 + \nu\rpar \lpar 1-2\nu\rpar } = {2\mu \nu\over 1-2\nu}\comma\; $$

are defined in terms of Young's modulus E and Poisson's ratio ν. The ratio of dynamic viscosity η to shear modulus μ defines the characteristic Maxwell relaxation time $t_{ {\rm {\rm r}}}$ over which stress decays in the material subject to constant strain loading:

(7)$$t_{{\rm {\rm r}}} = {\eta\over \mu} = {2 \eta \lpar 1 + \nu\rpar \over E}.$$

Previous authors studying timescales between the Maxwell time and the long-timescale viscous limit have allowed non-linear, Glen's law creep in the viscosity η (Goldberg and others, 2014; Robel and others, 2017). Here, by contrast, we study the cauldron system close to the Maxwell time t _r, such that the response to forcing is predominantly elastic with viscous deformation becoming apparent only later. For this reason, we take viscosity η to be linear. Linear viscosity also appears as a simplifying assumption in notable previous theoretical studies of glacial flow (Nye, 1970; Iken, 1981; Fowler, 1986) and is applied in the Burgers viscoelastic model of Reeh and others (2003). Lastly, choosing a linear viscoelastic constitutive relation simplifies our analysis by exploiting the Laplace transform correspondence of linear elastic and viscoelastic constitutive relations (Jull and McKenzie, 1996; Segall, 2010).

The Laplace transform ${\cal L}$ of a function g(t) and its inverse ${\cal L}^{-1}$ are given by

(8)$${\cal L}\lcub g\rcub = \overline{g}\lpar s\rpar = \int_0^{\infty} g\lpar t\rpar \, {\rm e}^{-st}\, {\rm d}t \comma\; $$

(9)$${\cal L}^{-1}\lcub \overline{g}\rcub = g\lpar t\rpar = {1\over 2\pi i} \int_{c-i\infty}^{c + i\infty} \overline{g}\lpar s\rpar \, {\rm e}^{sx}\, {\rm d} s\comma\; $$

respectively, where $t \in {\opf R}$, the Laplace variable $s \in {\opf C}$, and Re(s) = c (see e.g. Mathews and Walker, 1964).

Assuming negligible initial stress in the cauldron ice, the Laplace transform of Eqn (5) is

(10)$$\overline{\sigma}_{ij} = 2 {\eta s\over 1 + t_r s} \overline{\varepsilon}_{ij} + {\lambda + {2 \over 3} \mu + \lambda t_r s\over 1 + t_r s} \overline{\varepsilon}_{kk}\delta_{ij}\comma\; $$

where bars denote transformed variables. Defining the transformed Lamé parameters

(11)$$\overline{\mu} = \mu \left({t_{\rm r} s\over 1 + t_{\rm r} s}\right)\comma\; \quad \overline{\lambda} = \lambda + {2 \mu\over 3 \lpar 1 + t_{\rm r} s\rpar }\comma\; $$

assuming isotropy and neglecting shear, Eqn (10) takes the form of the constitutive relation for an isotropic, linear elastic material in transformed space (Landau and Lifshitz, 1959):

(12)$$\overline{\sigma}_{ij} = 2 \bar{\mu} \overline{\varepsilon}_{ij} + \bar{\lambda}\overline{\varepsilon}_{kk} \delta_{ij}.$$

This elastic-viscoelastic correspondence allows us to derive a linear viscoelastic model by taking the inverse Laplace transform of an elastic deformation equation. We therefore proceed with analyzing a linear elastic collapse in Laplace space.

For simplicity, we approximate the cauldron as a circular plate. A radially symmetric plate deforms according to the plate-bending equation,

(13)$$\bar{D} \nabla^4 \bar{w} = -\bar{\,f}\comma\; $$

where $\bar {D}$ is the transformed bending modulus, $\bar {w}$ is the transformed deflection from initial position, and $\,\,\bar {\!\!f}$ is the transformed loading (Landau and Lifshitz, 1959; Howell and others, 2009). Here, the plate is sagging downward under its own weight, and the loading in physical space is simply f = ρ _i g h, with ρ _i the density of glacier ice, g the acceleration due to gravity and h the ice thickness. The transformed loading is $\bar {f} = f/s$, with s the Laplace variable.

The Laplace-transformed elastic bending modulus in Eqn (13) is

(14)$$\bar{D} = {\bar{E} h^3\over 12 \lpar 1- \bar{\nu}^2\rpar }\comma\; $$

with bars again denoting transformed variables, and the transformed Young's modulus and Poisson's ratio

(15)$$\bar{E} = 2\bar{\mu} + {\bar{\mu}\bar{\lambda}\over \bar{\mu} + \bar{\lambda}}\comma\; \quad \bar{\nu} = {\nu\over s}.$$

Note that $\bar {D}$ depends on the Laplace variable s via the transformed parameters ($\bar {\lambda }\comma\; \bar {\mu }\comma\; \bar {\nu }$) so our eventual viscoelastic solution will have a time-dependent bending modulus D(t).

We choose a coordinate system with origin at the center of the collapsing cauldron, radial dimension r increasing outward and vertical dimension z increasing upward from the central neutral plane. In this coordinate system, the axisymmetric differential operator $\nabla ^4$ is

(16)$$\nabla^4 = {\partial^4\over \partial r^4} + {2\over r}{\partial^3\over \partial r^3} - {1\over r^2}{\partial^2\over \partial r^2} + {1\over r^3}{\partial\over \partial r}\comma\; $$

and we have the general solution to Eqn (13):

(17)$$\bar{w} = -{\bar{\,f}\over 64 \bar{D}} r^4 + c_1 \ln{r} + c_2 r^2 + {c_3\over 2} r^2 \ln{r} + c_4\comma\; $$

with c _i constants.

Around the edges of the cauldron (r = R), collapsing ice meets intact ice. Background viscous flow of the intact ice is on the order of 0.2 m d⁻¹, far slower than the deformation observed within the cauldron during its collapse (> 1 mh⁻¹; see Fig. 3), and we do not include it in this analysis. We therefore apply clamped boundary conditions, i.e.

(18)$$\left. \bar{w} \right\vert_{r = R} \equiv 0\comma\; $$

(19)$$\left. {\partial \bar{w}\over \partial r} \right\vert_{r = R} \equiv 0.$$

Subject to these conditions and the requirement that w(r = 0) be finite, we find the particular solution

(20)$$\bar{w} = -{\bar{\,f}\over 64\bar{D}} \left(r^2 - R^2 \right)^2\comma\; $$

(see also Landau and Lifshitz, 1959).

From Eqn (20), we can define the slope φ and in-plane displacement u _r as

(21)$$\varphi\lpar r\rpar = {\partial \bar{w}\over \partial r} = {\bar{\,f} r\over 16\bar{D}} \left(r^2 - R^2 \right)\comma\; $$

(22)$$u_r \lpar r\rpar = -\zeta \varphi\lpar r\rpar \comma\; $$

where ζ denotes vertical distance from the neutral plane of the plate, such that at the ice surface ζ = h/2. The in-plane radial and hoop strains are

(23)$$\varepsilon_{rr} = {\partial u_r\over \partial r} = -{\bar{\,f}\zeta\over 16 \bar{D}} \left(R^2 - 3r^2\right)\comma\; $$

(24)$$\varepsilon_{\theta \theta} = {u_r\over r} = - {\bar{\,f} \zeta\over 16 \bar{D}} \left(r^2 - R^2 \right)\comma\; $$

respectively, and tensile stress toward the cauldron center is then

(25)$$\sigma_{rr} = {\bar{E}\over 1- \bar{\nu}^2} \left(\varepsilon_{rr} + \nu \varepsilon_{\theta \theta}\right).$$

At last, taking the inverse Laplace transform, we find the expressions for viscoelastic deflection and stress:

(26)$$w\lpar r\comma\; t\rpar = -{\,f R^4\over 64 \, D\lpar t\rpar } \left({r^2\over R^2} -1 \right)^2\comma\; $$

(27)$$\sigma_{rr} \lpar r\comma\; t\rpar = {E f \zeta\over 16 \left(1- \nu^2\right)D\lpar t\rpar } \left(\lpar \nu + 1\rpar R^2 - \lpar \nu + 3\rpar r^2 \right).$$

In our simple viscoelastic cauldron collapse model, we use the symbolic mathematics package SymPy v0.7.6 to find $D\lpar t\rpar = {\cal L}^{-1}\lcub \bar {D}\lpar s\rpar \rcub$. Representative values for all parameters are summarized in Table 1. The thickness of the ice plate (h/R ≈ 0.2) suggests that some non-linear geometric effects are present in the natural system but not captured in our simple model (Howell and others, 2009), which could be remedied with a small ($\sim 15\percnt$) correction for thick-plate deformation (Wang and others, 2004).

Figure 6 shows three example transects taken across the cauldron, and the deformed surfaces we compute for each using Eqn (26). We have selected representative transects that cross the deepest point of the cauldron, start and end in intact ice around the cauldron edge, and include some visible crevasses in the 2015 surface. All three transects show an initial elastic drop (solid purple curves) accounting for part of the observed deformation, with subsequent viscous profiles (dashed curves) settling closer to the observed final configuration. The final viscoelastic profiles at time t = 4 d reasonably approximate the observed surface on transects I and II. That is, Eqn (26) produces realistic deformation and we can expect the corresponding stress computed from Eqn (27) will be realistic. On transect III, which has a shorter effective cauldron radius R, the viscoelastic profiles we compute underestimate the true subsidence. To mitigate over- and under-estimates of deformation and stress from slight asymmetry of the cauldron, we calculate stress using an idealized effective cauldron radius as described below.

Fig. 6. (a) Difference in eastern Skaftá cauldron surface elevation post-collapse versus pre-collapse (i.e. the difference of Figs 2a, b). (b) Three transects with observed surface elevations from 2012 (dotted black lines) to 2015 (solid black lines), and surfaces of idealized elastic (solid purple) and viscoelastic (dashed) collapse. Viscoelastic profiles shown are at 2 and 4 days after onset of collapse. All transects share horizontal and vertical scale, with 5:1 vertical exaggeration.

3.4. Radial stress and crevasse locations

We define the deepest point in the ArcticDEM observations as the cauldron center, and the area of minimal difference between 2012 and 2015 observations as the cauldron edge. We take the mean radial distance from edge to center as the cauldron ‘radius’ for the axisymmetric approximation. We calculate peak instantaneous stress along an idealized cauldron radius, at the moment the lake level drops and the cauldron loses support, according to Eqn (27). We then sample the observed ice surface elevation along 100 evenly-spaced radii from the center to the cauldron edge, and we identify crevasses using a one-dimensional analog to the crevasse-detection algorithm in Section 3.1.

Figure 7 summarizes the radial location of crevasses and corresponding peak radial stress along all 100 sampled cauldron radii. Again, as in Section 3.1, we find crevasses clustered in areas of high peak stress, with large areas of intact ice in between. Surface stress is lowest where the stress regime transitions from compressive (positive y-values in Fig. 7) to tensile (negative y-values). We note an area of intact ice on all radii between approximately 700 < r < 1050 m, where the peak radial surface stress we compute ranges from 15 MPa compressive to 8 MPa tensile. The peak radial stresses we compute using this approach span a slightly larger range than the surface maximum principal stresses computed in Section 3.1, but they are of comparable magnitude.

Fig. 7. Peak surface radial stress σ_rr (black curve) as a function of radial coordinate r (Eqn (27)). Vertical lines show locations of observed crevasses, with line color indicating stress regime. Positive stress values and blue colors indicate compression; negative stress values and red colors indicate tension. A gray overlay indicates a region of intact ice (no crevasses observed) at effective radii 700 < r < 1050 m.

4. Discussion

The cauldron collapse we study here resulted in pronounced deformation of the ice surface and large, deep surface crevasses (Figs. 2, 3). High-resolution observations of the ice surface from before and after the collapse allow us to estimate the surface elastic stress field throughout the cauldron (Section 3.1) and identify the distribution of maximum principal stress in crevassed and intact areas (Figs. 4, 5). With our linear viscoelastic model (Section 3.2), we are able to produce realistic post-collapse ice surface profiles (Fig. 6). The same model gives us an independent estimate of the largest-magnitude surface stress during collapse (Fig. 7), which agrees with our first-order elastic estimate in order of magnitude and radial pattern.

The crevasses we identify in the ArcticDEM surface (Porter and others, 2018) are several meters wide and tens of meters deep. Most crevasses are in concentric arcs around the outer rim of the cauldron, though there is another set of crevasses apparent near the cauldron center. The remaining cauldron surface appears generally intact. Indeed, our median filtering algorithm (Section 3.1) identifies 89% of the surface as intact ice, 9% as fractured and the remaining 2% areas of missing data. We find that crevasses coincide with areas where maximum stress was higher (Fig. 5).

The maximum instantaneous stresses we estimate in Sections 3.1 and 3.2 are large, of order 1–10 MPa. For comparison, the viscous beam stress calculated by Evatt and Fowler (2007) for rapid cauldron collapse is of order 1 MPa. Our methods in both sections do tend to overestimate stress, though for different reasons. In Section 3.1, we ignore time dependence and assume an instantaneous, fully elastic collapse. Because eventual viscous deformation would tend to reduce stress (Howell and others, 2009), Eqn (4) will tend to overestimate the elastic stress that contributed to the observed surface deformation. Nevertheless, we can constrain the magnitude of overestimate. For example, assuming cylindrical symmetry for simplicity, we can approximate the stress components of Eqns (1–2) as

(28)$$\sigma_{\rm e} \sim {Eh\over 2 \lpar 1-\nu^2\rpar } {1 \over R} {w\over R}\comma\; $$

with $\partial ^2_r w \sim \lpar {1}/{R}\rpar \lpar {w}/{R}\rpar$ replacing the curvature κ _ij and all other terms as before. We can then use representative parameter values from Table 1 in Eqn (28) to deduce that elastic subsidence of only 10 m would be sufficient to induce elastic stress up to 1 MPa. To induce instantaneous stress of up to 10 MPa requires predominantly-elastic subsidence on the order of 100 m. The GPS record (Fig. 3) confirms that subsidence of 10 m or more took place during the initial elastic phase of collapse (30 September), but that total subsidence approaching 100 m was not reached until the viscous phase (3 October). Thus, we conclude that instantaneous stress of order 10 MPa is an overestimate due to our simple method and unlikely in reality, but that instantaneous stress of 1 MPa and higher likely did arise during the elastic collapse.

In Section 3.2, we account for viscous effects in Eqn (27) but overestimate stress by assuming initially intact ice. Existing near-surface crevasses would reduce the effective thickness of the ice plate (h in Eqn (14)) and prevent stress transmission at the ice surface, where stress in an intact plate would be highest. Furthermore, pre-existing fractures would tend to concentrate stresses and thereby relieve stress in surrounding intact ice (Rice, 1968; Weertman, 1973). The minimum length a ₀ of pre-existing cracks that could concentrate stress depends on the stress and fracture toughness (see e.g. Liu and Miller, 1979):

(29)$$a_0 = {1\over \pi} \left({K_0 \sqrt{1-\nu^2}\over \Upsilon \sigma}\right)^2\comma\; $$

where σ is the applied stress, K ₀ is the fracture toughness, $\sqrt {1-\nu ^2}$ is a correction for plane strain stress conditions and $\Upsilon$ is a geometric factor. Here we use $\Upsilon = 1.12$ for an edge crack in semi-infinite geometry (Broberg, 1999). We can use Eqn (29) to compute the maximum stress σ that could be sustained before fracture. In initially intact ice, ice grain boundaries themselves could serve as initial ‘cracks’ concentrating stresses. With a typical grain size of glacier ice as the initial crack length, a ₀ ≈ 5 mm (Budd and Jacka, 1989), and a fracture toughness for glacier ice of K ₀ = 150 kPa m^1/2 (Rist and others, 1999), we find that applied stress σ ≥ 1 MPa would concentrate along grain boundaries to nucleate fractures. Pre-existing fractures of ~ 5 m in length could have concentrated stresses as low as 30 kPa, reducing even further the largest-magnitude stress that could have been active at the surface. Although observations in Guđmundsson and others (2016) and Porter and others (2018) do not indicate pre-existing surface fractures, we cannot rule out the presence of initial flaws at depth or smaller than the 2–4 m spatial resolution of those datasets. Based on initial flaw size analysis with Eqn (29), we conclude that peak instantaneous stress of order 10 MPa is an unrealistic overestimate, but that stress of order 1 MPa could have occurred and produced the surface crevasses observed after the 2015 event.

An alternative interpretation of the eastern Skaftá cauldron collapse might describe the observed subsidence as an entirely viscous phenomenon, with no elastic component. According to the usual power-law viscous rheology invoked for glacier ice (Glen, 1955), strain rate $\dot {\varepsilon }$ increases with the third power of deviatoric stress τ, which suggests that the cauldron ice could subside rapidly under the stress induced by loss of water pressure below. We would expect lower maximum stress in this case due to continual viscous-regime deformation. For example, experimental results (summarized in Goldsby and Kohlstedt, 2001) indicate that deviatoric stresses on the order of 10 MPa produce strain rates on the order of 10⁻³ s⁻¹, so the center of a cauldron of radius 1000 m could subside at rates exceeding 1 m s⁻¹ and rapidly relieve stress. However, the same relationship implies that stress on the order of 1 MPa would produce cauldron-center subsidence of only 0.3 m h⁻¹, which is insufficient for the sustained rapid subsidence (1−3 m h⁻¹) indicated in the GPS record (Fig. 3). Thus, Glen's law viscosity also suggests peak instantaneous stresses at the cauldron surface of 1 MPa and greater.

Even in light of the overestimates inherent in our methods, our analyses show that the collapse of eastern Skaftá cauldron and similar events could produce peak instantaneous stress equaling or exceeding the 0.7−3.1 MPa tensile strength of ice estimated in laboratory experiments summarized by Petrovic (2003). Despite such high stress, 89% of the cauldron surface area appears intact. Indicators of ice yield observed in other settings, such as the complete ring fractures and nearly vertical inner walls visible on the 1996 Gjálp eruption cauldrons (Guđmundsson and others, 2004; Evatt and Fowler, 2007), are not visible in the 2015 eastern Skaftá surface observations. Furthermore, we find intact ice with higher frequency than fractured ice at peak tensile stress of 1 MPa and even higher (Fig. 5). If the threshold stress for fracture were as low as 0.2 − 0.5 MPa, as suggested by the glacier-specific estimates of Vaughan (1993) and as used in continuum damage mechanics studies (Pralong and others, 2003; Duddu and Waisman, 2013; Krug and others, 2014), we would expect more surface crevasses throughout the cauldron than the eastern Skaftá observations indicate. Both our analysis and the viscous analysis of Evatt and Fowler (2007) support a tensile strength of glacier ice that is consistent with laboratory values rather than previous glacier-specific estimates.

Although we deduce tensile strength from the elastic stress active over the short timescale of crevasse nucleation and propagation, our results are generalizable to glacier flow modeling. A glacier tensile strength on the order of 1 MPa is easily implementable in short-term, process-scale modeling such as Åström and others (2013), which already uses elastic bonds between ice ‘grains’ and tests fracture thresholds up to 1 MPa. Stress fluctuations in longer-term, viscous ice-sheet modeling (Krug and others, 2014; Jouvet and others, 2011) are generally too small to reach a fracture threshold of 1 MPa. Yet future refinements in continuum representation of fracture and damage could implement a fracture threshold informed by the observations and analysis we present here. In particular, a Maxwell viscoelastic rheology accounts for both short-term elastic effects and longer-term viscous deformation. That is, localized, short-term (t < t _r) increases in stress can generate an elastic response and propagate crevasses even as background flow remains viscous in response to the global stress field. We do not suggest that all ice-sheet models be redeveloped to incorporate viscoelasticity. However, previous authors have successfully applied a viscoelastic rheology to model ice damage evolution over hours to days (Mobasher and others, 2016), tidal variability of ice stream flow over days to weeks (Walker and others, 2012b; Rosier and others, 2014, 2015; Robel and others, 2017) and ice stream and ice shelf motion over months to years (Reeh and others, 2003; Gudmundsson, 2011; Goldberg and others, 2014; MacAyeal and others, 2015; Banwell and MacAyeal, 2015).

5. Conclusions

We have applied two complementary methods to constrain the tensile strength of glacier ice from remote sensing and in situ observations. Our analysis suggests that the 2015 collapse of eastern Skaftá cauldron, Iceland, induced tensile stress on the order of 1 MPa over much of the cauldron surface. That stress, together with pre-existing flaws, produced a set of crevasses around the cauldron rim and center but left much of the cauldron ice apparently intact. Our findings support an estimate of ice tensile strength on the order of 1 MPa, broadly consistent with laboratory estimates but exceeding previous estimates from in situ observations of glaciers. As numerical model development advances toward more physically consistent representation of glacier flow and fracture, we suggest that model parameters be brought in line with natural-scale observations such as those presented here.