Hostname: page-component-89b8bd64d-n8gtw Total loading time: 0 Render date: 2026-05-08T09:09:47.244Z Has data issue: false hasContentIssue false
  • English
  • Français

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

Published online by Cambridge University Press:  26 August 2020

Lizz Ultee Open the ORCID record for Lizz Ultee [Opens in a new window] ,
Colin Meyer  and
Brent Minchew Open the ORCID record for Brent Minchew [Opens in a new window]
Lizz Ultee*
Affiliation:
Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA
Colin Meyer
Affiliation:
Thayer School of Engineering, Dartmouth College, Hanover, NH, USA
Brent Minchew
Affiliation:
Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA
*
Author for correspondence: Lizz Ultee, E-mail: ehultee@mit.edu
Rights & Permissions [Opens in a new window]

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.

Keywords

Crevassesglacier mechanicsglacier modeling

Information

Type
Article
Information
Journal of Glaciology , Volume 66 , Issue 260 , December 2020 , pp. 1024 - 1033
Check for updates
Creative Commons
Creative Common License - CCCreative Common License - BYCreative Common License - NCCreative Common License - SA
This is an Open Access article, distributed under the terms of the Creative Commons Attribution-NonCommercial-ShareAlike licence (http://creativecommons.org/licenses/by-nc-sa/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the same Creative Commons licence is included and the original work is properly cited. The written permission of Cambridge University Press must be obtained for commercial re-use.
Copyright
Copyright © The Author(s), 2020. Published by Cambridge University Press

1. Introduction

Iceberg calving is a key source of uncertainty in future projections of global sea-level rise (Church and others, Reference Church, Stocker, Qin, Plattner, Tignor, Allen, Boschung, Nauels, Xia, Bex and Midgley2013). 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, Reference Pralong, Funk and Lüthi2003; Duddu and Waisman, Reference Duddu and Waisman2013; Åström and others, Reference Åström2013; Krug and others, Reference Krug, Weiss, Gagliardini and Durand2014; Benn and others, Reference Benn2017). 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, Reference Currier and Schulson1982; Lee and Schulson, Reference Lee and Schulson1988; Xian and others, Reference Xian, Chu, Scavuzzo and Srivatsan1989; Druez and others, Reference Druez, McComber and Tremblay1989) and field observations (Vaughan, Reference Vaughan1993), 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, Reference Bassis and Ma2015; Ma and others, Reference Ma, Tripathy and Bassis2017; Borstad and others, Reference Borstad, McGrath and Pope2017).

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, Reference Petrovic2003) for laboratory and lake ice; Druez and others (Reference Druez, McComber and Tremblay1989) 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, Reference Vaughan1993). Part of the discrepancy may be the grain-size dependence of the tensile strength of ice (Currier and Schulson, Reference Currier and Schulson1982). 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, Reference Currier and Schulson1982) but only at multi-meter scale in existing field and remote-sensing observations (Colgan and others, Reference Colgan2016). Another source of discrepancy is that stress is difficult to observe in situ (Pfeffer and others, Reference Pfeffer, Humphrey, Amadei, Harper and Wegmann2000; Colgan and others, Reference Colgan2016), so estimates derived from glacier-scale observations, such as Vaughan (Reference Vaughan1993), 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, Reference Reeh, Christensen, Mayer and Olesen2003; Borstad and others, Reference Borstad, McGrath and Pope2017). Further, Dempsey and others (Reference Dempsey, Adamson and Mulmule1999) 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 (Reference Åström2013, Reference Åström2014) and Benn and others (Reference Benn2017) 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, Reference Benn2017). 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, Reference Currier and Schulson1982), it is unclear that the stress threshold σc is the same physical quantity measurable in the laboratory or field.

By contrast, Krug and others (Reference Krug, Weiss, Gagliardini and Durand2014) use a tensile stress threshold for fracture of 0.20 MPa in their continuum damage mechanics framework. Pralong and others (Reference Pralong, Funk and Lüthi2003), Pralong and Funk (Reference Pralong and Funk2005) and Duddu and Waisman (Reference Duddu and Waisman2013) use a similar threshold, with tensile strength ranging from 0.20 to 0.50 MPa, based on the macroscale condition derived by Vaughan (Reference Vaughan1993). This entire range is below even the lowest laboratory estimate of ice tensile strength (Petrovic, Reference Petrovic2003). 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 (Reference Banwell, Willis, Macdonald, Goodsell and MacAyeal2019) 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 (Reference Banwell, MacAyeal and Sergienko2013) used the elastic-plate analysis derived by Sergienko (Reference Sergienko2005) to compute the fracture spacing associated with supraglacial lake drainage on Antarctic Peninsula ice shelves. Evatt and Fowler (Reference Evatt and Fowler2007) computed fracture spacing associated with cauldron collapse on the basis of viscous beam theory derived in Evatt and others (Reference Evatt, Fowler, Clark and Hulton2006). Even extraterrestrial ice subsidence has supported studies of ice stress and fracture: Walker and Schmidt (Reference Walker and Schmidt2015) 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, Reference Walker, Bassis and Liemohn2012a).

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, Reference Gudmundsson2011; Goldberg and others, Reference Goldberg, Schoof and Sergienko2014; MacAyeal and others, Reference MacAyeal, Sergienko and Banwell2015; Robel and others, Reference Robel, Tsai, Minchew and Simons2017), 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, Reference Porter2018) 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, Reference Guđmundsson, Högnadóttir and Rossi2016) and finally is drained by glacial outburst floods (jökulhlaups) lasting hours to days (Björnsson, Reference Björnsson1992). 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, Reference Einarsson2016).

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, Reference Amante and Eakins2009), 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, Reference Wessel, Smith, Scharroo, Luis and Wobbe2013).

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 m3 s−1 (Jóhannesson and others, Reference Jóhannesson2016). 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, Reference Porter2018) 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, Reference Porter2018); 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, Reference Porter2018) 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, Reference Guđmundsson, Magnússon, Högnadóttir, Pálsson and Rossi2018). 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−1, followed by a prolonged period of slower (0.5 m h−1) 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, Reference Reeh, Christensen, Mayer and Olesen2003; Gudmundsson, Reference Gudmundsson2011; Rosier and others, Reference Rosier, Gudmundsson and Green2014). 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, Reference Ugural2018). 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, Reference Ugural2018):

(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, Reference Vaughan1995; 4–10 GPa in Rist and others, Reference Rist, Sammonds, Oerter and Doake2002).

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 σ1 > 0 is compression and σ1 < 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, Reference Porter2018) 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 (Reference Guđmundsson, Högnadóttir and Rossi2016) 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, Reference Porter2018). 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, Reference Gudmundsson2011; Goldberg and others, Reference Goldberg, Schoof and Sergienko2014; MacAyeal and others, Reference MacAyeal, Sergienko and Banwell2015; Banwell and MacAyeal, Reference Banwell and MacAyeal2015; Robel and others, Reference Robel, Tsai, Minchew and Simons2017); in particular, Gudmundsson (Reference Gudmundsson2011) shows that Maxwell viscoelasticity is an appropriate simplification from Burgers viscoelasticity as implemented by Reeh and others (Reference Reeh, Christensen, Mayer and Olesen2003).

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 (Reference Howell, Kozyreff and Ockendon2009), 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, Reference Goldberg, Schoof and Sergienko2014; Robel and others, Reference Robel, Tsai, Minchew and Simons2017). 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, Reference Nye1970; Iken, Reference Iken1981; Fowler, Reference Fowler1986) and is applied in the Burgers viscoelastic model of Reeh and others (Reference Reeh, Christensen, Mayer and Olesen2003). 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, Reference Jull and McKenzie1996; Segall, Reference Segall2010).

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, Reference Mathews and Walker1964).

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, Reference Landau and Lifshitz1959):

(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, Reference Landau and Lifshitz1959; Howell and others, Reference Howell, Kozyreff and Ockendon2009). 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−1, far slower than the deformation observed within the cauldron during its collapse (> 1 mh−1; 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, Reference Landau and Lifshitz1959).

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, Reference Howell, Kozyreff and Ockendon2009), which could be remedied with a small ($\sim 15\percnt$) correction for thick-plate deformation (Wang and others, Reference Wang, Nazmul and Wang2004).

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, Reference Porter2018) 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 (Reference Evatt and Fowler2007) 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, Reference Howell, Kozyreff and Ockendon2009), 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, Reference Rice and Liebowitz1968; Weertman, Reference Weertman1973). The minimum length a 0 of pre-existing cracks that could concentrate stress depends on the stress and fracture toughness (see e.g. Liu and Miller, Reference Liu and Miller1979):

(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 0 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, Reference Broberg1999). 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 0 ≈ 5 mm (Budd and Jacka, Reference Budd and Jacka1989), and a fracture toughness for glacier ice of K 0 = 150 kPa m1/2 (Rist and others, Reference Rist1999), 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 (Reference Guđmundsson, Högnadóttir and Rossi2016) and Porter and others (Reference Porter2018) 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, Reference Glen1955), 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, Reference Goldsby and Kohlstedt2001) indicate that deviatoric stresses on the order of 10 MPa produce strain rates on the order of 10−3 s−1, so the center of a cauldron of radius 1000 m could subside at rates exceeding 1 m s−1 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−1, which is insufficient for the sustained rapid subsidence (1−3 m h−1) 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 (Reference Petrovic2003). 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, Reference Guđmundsson, Sigmundsson, Björnsson and Högnadóttir2004; Evatt and Fowler, Reference Evatt and Fowler2007), 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 (Reference Vaughan1993) and as used in continuum damage mechanics studies (Pralong and others, Reference Pralong, Funk and Lüthi2003; Duddu and Waisman, Reference Duddu and Waisman2013; Krug and others, Reference Krug, Weiss, Gagliardini and Durand2014), 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 (Reference Evatt and Fowler2007) 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 (Reference Åström2013), 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, Reference Krug, Weiss, Gagliardini and Durand2014; Jouvet and others, Reference Jouvet, Picasso, Rappaz, Huss and Funk2011) 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, Reference Mobasher, Duddu, Bassis and Waisman2016), tidal variability of ice stream flow over days to weeks (Walker and others, Reference Walker, Christianson, Parizek, Anandakrishnan and Alley2012b; Rosier and others, Reference Rosier, Gudmundsson and Green2014, Reference Rosier, Gudmundsson and Green2015; Robel and others, Reference Robel, Tsai, Minchew and Simons2017) and ice stream and ice shelf motion over months to years (Reeh and others, Reference Reeh, Christensen, Mayer and Olesen2003; Gudmundsson, Reference Gudmundsson2011; Goldberg and others, Reference Goldberg, Schoof and Sergienko2014; MacAyeal and others, Reference MacAyeal, Sergienko and Banwell2015; Banwell and MacAyeal, Reference Banwell and MacAyeal2015).

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.

Acknowledgements

The data presented in this work can be downloaded from the ArcticDEM data repository (http://data.pgc.umn.edu/elev/dem/setsm/ArcticDEM/mosaic/v3.0/2m/) and upon request from the Icelandic Meteorological Office. Code for simple viscoelastic cauldron model and crevasse analysis can be inspected and downloaded from a public GitHub repository (http://github.com/ehultee/VE-cauldrons). GPS observations during the collapse of the cauldron were made within the framework of FutureVolc, funded by the European Union's Seventh Programme under grant agreement No. 308377. We thank Benedikt G. Ófeigsson and Vilhjálmur Kjartansson for their contributions to field work and data processing. We thank Tómas Jóhannesson (Icelandic Meteorological Office) and four anonymous reviewers for their comments on the manuscript. The authors have declared that no conflict of interest exists.

References

Amante, C and Eakins, BW (2009) ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. Technical Report NESDIS NGDC-24, NOAA Technical Memorandum, National Geophysical Data Center, NOAA (doi 10.7289/V5C8276M), accessed November 2012.Google Scholar
Åström, JA and 6 others (2013) A particle based simulation model for glacier dynamics. The Cryosphere 7(5), 1591–1602. doi 10.5194/tc-7-1591-2013.CrossRefGoogle Scholar
Åström, JA and 10 others (2014) Termini of calving glaciers as self-organized critical systems. Nature Geoscience 7, 874–878. doi 10.1038/NGEO2290.CrossRefGoogle Scholar
Banwell, AF and MacAyeal, DR (2015) Ice-shelf fracture due to viscoelastic flexure stress induced by fill/drain cycles of supraglacial lakes. Antarctic Science 27(6), 587–597. doi 10.1017/S0954102015000292.CrossRefGoogle Scholar
Banwell, AF, MacAyeal, DR and Sergienko, OV (2013) Breakup of the Larsen B Ice Shelf triggered by chain reaction drainage of supraglacial lakes. Geophysical Research Letters 40(22), 5872–5876. doi 10.1002/2013GL057694.CrossRefGoogle Scholar
Banwell, AF, Willis, IC, Macdonald, GJ, Goodsell, B and MacAyeal, DR (2019) Direct measurements of ice-shelf flexure caused by surface meltwater ponding and drainage. Nature Communications 10(1), 730. doi 10.1038/s41467-019-08522-5.CrossRefGoogle ScholarPubMed
Bassis, JN and Ma, Y (2015) Evolution of basal crevasses links ice shelf stability to ocean forcing. Earth and Planetary Science Letters 409, 203–211. doi 10.1016/j.epsl.2014.11.003.CrossRefGoogle Scholar
Benn, DI and 7 others (2017) Melt-under-cutting and buoyancy-driven calving from tidewater glaciers: new insights from discrete element and continuum model simulations. Journal of Glaciology 63(240), 691–702. doi 10.1017/jog.2017.41.CrossRefGoogle Scholar
Björnsson, H (1992) Jökulhlaups in Iceland: prediction, characteristics and simulation. Annals of Glaciology 16, 95–106. doi 10.3189/1992AoG16-1-95-106.CrossRefGoogle Scholar
Borstad, C, McGrath, D and Pope, A (2017) Fracture propagation and stability of ice shelves governed by ice shelf heterogeneity. Geophysical Research Letters 44(9), 4186–4194. doi 10.1002/2017GL072648.CrossRefGoogle Scholar
Broberg, KB (1999) Cracks and Fracture. Cambridge, UK: Cambridge University Press, ISBN 978-0-12-134130-5Google Scholar
Budd, WF and Jacka, TH (1989) A review of ice rheology for ice sheet modelling. Cold Regions Science and Technology 16(2), 107–144. doi 10.1016/0165-232X(89)90014-1.CrossRefGoogle Scholar
Church, JA and 12 others (2013) Sea level change. In Stocker, TF, Qin, D, Plattner, GK, Tignor, M, Allen, SK, Boschung, J, Nauels, A, Xia, Y, Bex, V and Midgley, PM (eds), Climate Change 2013: The Physical Science Basis. Cambridge, United Kingdom and New York, NY, USA: Cambridge University Press, pp. 1137– 1216.Google Scholar
Colgan, W and 6 others (2016) Glacier crevasses: observations, models and mass balance implications. Reviews of Geophysics 54(1), 119–161. doi 10.1002/2015RG000504.CrossRefGoogle Scholar
Currier, JH and Schulson, EM (1982) The tensile strength of ice as a function of grain size. Acta Metallurgica 30(8), 1511–1514. doi 10.1016/0001-6160(82)90171-7.CrossRefGoogle Scholar
Dempsey, JP, Adamson, RM and Mulmule, SV (1999) Scale effects on the in-situ tensile strength and fracture of ice. Part II: first-year sea ice at Resolute, N.W.T. International Journal of Fracture 95(1), 347. doi 10.1023/A:1018650303385.CrossRefGoogle Scholar
Druez, J, McComber, P and Tremblay, C (1989) Experimental results on the tensile strength of atmospheric ice. Transactions of the Canadian Society for Mechanical Engineering 13(3), 59–64. doi 10.1139/tcsme-1989-0010.CrossRefGoogle Scholar
Duddu, R and Waisman, H (2013) A nonlocal continuum damage mechanics approach to simulation of creep fracture in ice sheets. Computational Mechanics 51(6), 961–974. doi 10.1007/s00466-012-0778-7.CrossRefGoogle Scholar
Einarsson, B and 5 others (2016) A spectrum of jökulhlaup dynamics revealed by GPS measurements of glacier surface motion. Annals of Glaciology 57(72), 47–61. doi 10.1017/aog.2016.8.CrossRefGoogle Scholar
Evatt, GW and Fowler, AC (2007) Cauldron subsidence and subglacial floods. Annals of Glaciology 45(1), 163–168. doi 10.3189/172756407782282561.Google Scholar
Evatt, GW, Fowler, AC, Clark, CD and Hulton, NRJ (2006) Subglacial floods beneath ice sheets. Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 364(1844), 1769–1794. doi 10.1098/rsta.2006.1798.Google ScholarPubMed
Fowler, AC (1986) A sliding law for glaciers of constant viscosity in the presence of subglacial cavitation. Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 407(1832), 147–170. doi 10.1098/rspa.1986.0090.Google Scholar
Glen, JW (1955) The creep of polycrystalline ice. Proceedings of the Royal Society of London A 228(1175), 519–538. doi 10.1098/rspa.1955.0066.Google Scholar
Goldberg, DN, Schoof, C and Sergienko, OV (2014) Stick-slip motion of an Antarctic ice stream: the effects of viscoelasticity. Journal of Geophysical Research: Earth Surface 119(7), 1564–1580. doi 10.1002/2014JF003132.Google Scholar
Goldsby, DL and Kohlstedt, DL (2001) Superplastic deformation of ice: experimental observations. Journal of Geophysical Research: Solid Earth 106(B6), 11017–11030. doi 10.1029/2000JB900336.CrossRefGoogle Scholar
Gudmundsson, GH (2011) Ice-stream response to ocean tides and the form of the basal sliding law. The Cryosphere 5(1), 259–270. doi 10.5194/tc-5-259-2011.CrossRefGoogle Scholar
Guđmundsson, MT, Högnadóttir, T and Rossi, C (2016) Report D5.7 - Ice evolution time series of elevation changes during caldera collapse and geothermal activity in Iceland. Technical report, FutureVolc, available online at http://futurevolc.hi.is/sites/futurevolc.hi.is/files/Pdf/Deliverables/f%v_d_5_7_to_submit_low_1.pdf.Google Scholar
Guđmundsson, M, Magnússon, E, Högnadóttir, T, Pálsson, F and Rossi, C (2018) Hættumat vegna jökulhlaupa í Skaftá: Skaftárkatlar, saga og þróun 1938–2018 [Hazard assessment for jökulhlaups in Skaftá: The Skaftá Cauldrons, history and development]. Technical Report Rep. RH-16-2019, Institute of Earth Sciences, University of Iceland, available online at https://www.vedur.is/media/vedurstofan-utgafa-2018/VI_2018_017_rs.pdf.Google Scholar
Guđmundsson, MT, Sigmundsson, F, Björnsson, H and Högnadóttir, T (2004) The 1996 eruption at Gjálp, Vatnajökull ice cap, Iceland: efficiency of heat transfer, ice deformation and subglacial water pressure. Bulletin of Volcanology 66(1), 46–65. doi 10.1007/s00445-003-0295-9.CrossRefGoogle Scholar
Howell, P, Kozyreff, G and Ockendon, J (2009) Applied Solid Mechanics. New York, USA: Cambridge University Press.Google Scholar
Iken, A (1981) The effect of the subglacial water pressure on the sliding velocity of a glacier in an idealized numerical model. Journal of Glaciology 27(97), 407–421. doi 10.3189/S0022143000011448.CrossRefGoogle Scholar
Jóhannesson, T and 6 others (2016) A record jökulhlaup from the subglacial lake beneath the eastern Skaftá cauldron, Vatnajökull ice cap, Iceland. Sixth Mars Polar Science Conference, pp. 1–2, available online at https://www.hou.usra.edu/meetings/marspolar2016/pdf/6109.pdf.Google Scholar
Jouvet, G, Picasso, M, Rappaz, J, Huss, M and Funk, M (2011) Modelling and numerical simulation of the dynamics of glaciers including local damage effects. Mathematical Modelling of Natural Phenomena 6, 263–280. doi 10.1051/mmnp/20116510.CrossRefGoogle Scholar
Jull, M and McKenzie, D (1996) The effect of deglaciation on mantle melting beneath Iceland. Journal of Geophysical Research: Solid Earth 101(B10), 21815–21828. doi 10.1029/96JB01308.CrossRefGoogle Scholar
Krug, J, Weiss, J, Gagliardini, O and Durand, G (2014) Combining damage and fracture mechanics to model calving. The Cryosphere 8(6), 2101–2117. doi 10.5194/tc-8-2101-2014.CrossRefGoogle Scholar
Landau, L and Lifshitz, EM (1959, 1986) Theory of Elasticity, volume 7 of Course of Theoretical Physics, 2nd Edn. New York, USA: Butterworth-Heinemann.Google Scholar
Lee, RW and Schulson, EM (1988) The strength and ductility of ice under tension. Journal of Offshore Mechanics and Arctic Engineering 110(2), 187–191, ISSN 0892-7219 doi: doi 10.1115/1.3257049.CrossRefGoogle Scholar
Liu, HW and Miller, KJ (1979) Fracture toughness of fresh-water ice. Journal of Glaciology 22(86), 135–143. doi 10.3189/S0022143000014118.CrossRefGoogle Scholar
Ma, Y, Tripathy, CS and Bassis, JN (2017) Bounds on the calving cliff height of marine terminating glaciers. Geophysical Research Letters 44(3), 1369–1375. doi 10.1002/2016GL071560.CrossRefGoogle Scholar
MacAyeal, DR, Sergienko, OV and Banwell, AF (2015) A model of viscoelastic ice-shelf flexure. Journal of Glaciology 61(228), 635–645. doi 10.3189/2015JoG14J169.CrossRefGoogle Scholar
Mathews, J and Walker, RL (1964) Mathematical Methods of Physics. New York, NY and Amsterdam, Netherlands: W. A. Benjamin Inc, ISBN 978-0805370027.Google Scholar
Mobasher, ME, Duddu, R, Bassis, JN and Waisman, H (2016) Modeling hydraulic fracture of glaciers using continuum damage mechanics. Journal of Glaciology 62(234), 794–804. doi 10.1017/jog.2016.68.CrossRefGoogle Scholar
Nye, JF (1970) Glacier sliding without cavitation in a linear viscous approximation. Proceeding of the Royal Society of London. Series A, Mathematical and Physical Sciences 315(1522), 381–403. doi 10.1098/rspa.1970.0050.Google Scholar
Petrovic, JJ (2003) Mechanical properties of ice and snow. Journal of Materials Science 38(1), 1–6, ISSN 1573-4803. doi: doi 10.1023/A:1021134128038.CrossRefGoogle Scholar
Pfeffer, WT, Humphrey, NF, Amadei, B, Harper, J and Wegmann, J (2000) In situ stress tensor measured in an Alaskan glacier. Annals of Glaciology 31, 229–235. doi 10.3189/172756400781820354.CrossRefGoogle Scholar
Porter, C and 29 others (2018) ArcticDEM v1. Harvard Dataverse (doi 10.7910/DVN/OHHUKH).Google Scholar
Pralong, A and Funk, M (2005) Dynamic damage model of crevasse opening and application to glacier calving. Journal of Geophysical Research: Solid Earth 110, article number B01309. doi: 10.1029/2004JB003104.CrossRefGoogle Scholar
Pralong, A, Funk, M and Lüthi, MP (2003) A description of crevasse formation using continuum damage mechanics. Annals of Glaciology 37, 77–82. doi 10.3189/172756403781816077.CrossRefGoogle Scholar
Reeh, N, Christensen, EL, Mayer, C and Olesen, OB (2003) Tidal bending of glaciers: a linear viscoelastic approach. Annals of Glaciology 37(1), 83–89. doi 10.3189/172756403781815663.CrossRefGoogle Scholar
Rice, JR (1968) Mathematical analysis in the mechanics of fracture. In Liebowitz, H ed. Fracture: An Advanced Treatise, Volume 2, Chapter 3. New York, USA: Academic Press, pp. 191–311.Google Scholar
Rist, MA and 6 others (1999) Experimental and theoretical fracture mechanics applied to Antarctic ice fracture and surface crevassing. Journal of Geophysical Research 104(B2), 2973–2987. doi 10.1029/1998JB900026.CrossRefGoogle Scholar
Rist, MA, Sammonds, PR, Oerter, H and Doake, CSM (2002) Fracture of Antarctic shelf ice. Journal of Geophysical Research: Solid Earth 107(B1), ECV 2–1–ECV 2–13. doi 10.1029/2000JB000058.CrossRefGoogle Scholar
Robel, AA, Tsai, VC, Minchew, BM and Simons, M (2017) Tidal modulation of ice shelf buttressing stresses. Annals of Glaciology 58(74), 12–20. doi 10.1017/aog.2017.22.CrossRefGoogle Scholar
Rosier, SHR, Gudmundsson, GH and Green, JAM (2014) Insights into ice stream dynamics through modelling their response to tidal forcing. The Cryosphere 8(5), 1763–1775. doi 10.5194/tc-8-1763-2014.Google Scholar
Rosier, SHR, Gudmundsson, GH and Green, JAM (2015) Temporal variations in the flow of a large Antarctic ice-stream controlled by tidally induced changes in the subglacial water system. The Cryosphere 9(4), 1649–1661. doi 10.5194/tc-9-1649-2015.CrossRefGoogle Scholar
Segall, P (2010) Earthquake and volcano deformation. Princeton, NJ, USA: Princeton University Press.CrossRefGoogle Scholar
Sergienko, OV (2005) Surface melting of ice shelves and icebergs. Ph.D. thesis, University of Chicago, available online at https://search.proquest.com/docview/305419795.Google Scholar
Ugural, AC (2018) Plates and Shells: Theory and Analysis, 4th Edn. CRC Press: Applied & Computational Mechanics, ISBN 9781138032453.Google Scholar
Vaughan, DG (1993) Relating the occurrence of crevasses to surface strain rates. Journal of Glaciology 39(132), 255–266. doi 10.3189/S0022143000015926.CrossRefGoogle Scholar
Vaughan, DG (1995) Tidal flexure at ice shelf margins. Journal of Geophysical Research: Solid Earth 100(B4), 6213–6224. doi 10.1029/94JB02467.CrossRefGoogle Scholar
Walker, C, Bassis, J and Liemohn, M (2012 a) On the application of simple rift basin models to the south polar region of Enceladus. Journal of Geophysical Research: Planets 117, E07003. doi 10.1029/2012JE004084.CrossRefGoogle Scholar
Walker, RT, Christianson, K, Parizek, BR, Anandakrishnan, S and Alley, RB (2012 b) A viscoelastic flowline model applied to tidal forcing of Bindschadler Ice Stream, West Antarctica. Earth and Planetary Science Letters 319-320, 128–132. doi 10.1016/j.epsl.2011.12.019.CrossRefGoogle Scholar
Walker, CC and Schmidt, BE (2015) Ice collapse over trapped water bodies on Enceladus and Europa. Geophysical Research Letters 42(3), 712–719. doi 10.1002/2014GL062405.Google Scholar
Wang, CM, Nazmul, IM and Wang, Q (2004) Exact bending solutions of axisymmetric reissner plates in terms of classical thin plate solutions. Advances in Structural Engineering 7(2), 129–145. doi 10.1260/1369433041211075.CrossRefGoogle Scholar
Weertman, J (1973) Can a water-filled crevasse reach the bottom surface of a glacier? Proceedings of Cambridge Symposium on the Hydrology of Glaciers, 95, pp. 139–145, International Association of Hydrological Sciences, Cambridge, UK.Google Scholar
Wessel, P, Smith, WHF, Scharroo, R, Luis, JF and Wobbe, F (2013) Generic mapping tools: improved version released. Eos, Transactions American Geophysical Union 94(45), 409–410. doi 10.1002/2013EO450001.CrossRefGoogle Scholar
Xian, X, Chu, ML, Scavuzzo, RJ and Srivatsan, TS (1989) An experimental evaluation of the tensile strength of impact ice. Journal of Materials Science Letters 8(10), 1205–1208. doi 10.1007/BF01730071.CrossRefGoogle Scholar
Figure 0

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

Figure 1

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).

Figure 2

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).

Figure 3

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.

Figure 4

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.

Figure 5

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.

Figure 6

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.

Figure 7

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.