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

*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 Liemohn2012*a*).

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

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^{3} 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.

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

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

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

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

### 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

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

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:

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

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

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

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

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,

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

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

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

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

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.

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

(see also Landau and Lifshitz, Reference Landau and Lifshitz1959).

From Eqn (20), we can define the slope φ and in-plane displacement *u* _{r} as

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

respectively, and tensile stress toward the cauldron center is then

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

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.

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

## 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

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

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 m^{1/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 Alley2012*b*; 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.