Hostname: page-component-76fb5796d-wq484 Total loading time: 0 Render date: 2024-04-29T01:41:41.068Z Has data issue: false hasContentIssue false
  • English
  • Français

Basal hydrofractures near sticky patches

Published online by Cambridge University Press:  21 September 2022

Hanwen Zhang ,
Timothy Davis Open the ORCID record for Timothy Davis [Opens in a new window] ,
Richard F. Katz ,
Laura A. Stevens Open the ORCID record for Laura A. Stevens [Opens in a new window]  and
Dave A. May
Hanwen Zhang*
Affiliation:
Department of Earth Sciences, University of Oxford, Oxford, UK
Timothy Davis
Affiliation:
Department of Earth Sciences, University of Oxford, Oxford, UK
Richard F. Katz
Affiliation:
Department of Earth Sciences, University of Oxford, Oxford, UK
Laura A. Stevens
Affiliation:
Department of Earth Sciences, University of Oxford, Oxford, UK
Dave A. May
Affiliation:
Scripps Institution of Oceanography, University of California San Diego, La Jolla, CA, USA
*
Author for correspondence: Hanwen Zhang, E-mail: hanwen.zhang@earth.ox.ac.uk
Rights & Permissions [Opens in a new window]

Abstract

Basal crevasses are macroscopic structural discontinuities at the base of ice sheets and glaciers that arise by fracture. Motivated by observations and by the mechanics of elastic fracture, we hypothesise that spatial variations in basal stress (in the presence of basal water pressure) can promote and localise basal crevassing. We quantify this process in the theoretical context of linear elastic fracture mechanics. We develop a model evaluating the effect of shear-stress variation on the growth of basal crevasses. Our results indicate that sticky patches promote the propagation of basal crevasses, increase their length of propagation into the ice and, under some conditions, give them curved trajectories that incline upstream. A detailed exploration of the parameter space is conducted to gain a better understanding of the conditions under which sticky-patch-induced basal crevassing is expected beneath ice sheets and glaciers.

Keywords

Calvingcrevassesglacier modellingice-sheet modelling
Type
Article
Information
Journal of Glaciology , Volume 69 , Issue 275 , June 2023 , pp. 475 - 486
Check for updates
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
Copyright © The Author(s), 2022. Published by Cambridge University Press

1. Introduction

Crevasses are deep cracks in glaciers, ice sheets and ice shelves. Based on their position in the ice, they can be classified into surface crevasses and basal crevasses. On ice shelves, the propagation of surface and basal crevasses is a precursor to calving and rifting, and thus affects the stability of the ice shelf (Bassis and Ma, Reference Bassis and Ma2015; Lai and others, Reference Lai2020; Lipovsky, Reference Lipovsky2020). In grounded ice, basal crevasses may have important influences on subglacial hydrology and glacier dynamics (Walter and others, Reference Walter, Dalban Canassy, Husen and Clinton2013). Crevasses are also a potential explanation for seismic activity detected in glaciers and ice sheets (Hudson and others, Reference Hudson2020). Clustered, deep icequakes with non-double-couple sources have been detected in alpine glaciers (Walter and others, Reference Walter, Dalban Canassy, Husen and Clinton2013; Helmstetter and others, Reference Helmstetter, Moreau, Nicolas, Comon and Gay2015). Because these icequakes have a ‘substantial, if not dominant, isotropic component in the moment tensor’, they cannot be explained by shearing in stick-slip motion along the ice–bed interface. Instead, Walter and others (Reference Walter, Dalban Canassy, Husen and Clinton2013) attributed these icequakes to opening-mode fractures deep within the ice. Basal crevasses have been studied less than surface crevasses because observations of basal crevasses are few, and complicated, heterogeneous conditions characterise the ice–bed interface. The propagation of basal crevasses can be very sensitive to basal conditions, including subglacial hydrological conditions and bedrock rigidity (Jimenez and Duddu, Reference Jimenez and Duddu2018).

Various observations document the existence of basal crevasses. At Bench Glacier, Alaska, crevasses with connections to the bed were detected by both drilling experiments and radar imaging (Harper and others, Reference Harper, Bradford, Humphrey and Meierbachtol2010). Such crevasses, filled by high-pressure water, serve as important components of the subglacial hydrology system. In West Antarctica, Wearing and Kingslake (Reference Wearing and Kingslake2019) detected relic basal crevasses in grounded regions of the Henry Ice Rise in the Ronne Ice Shelf using ice-penetrating radar. To consistently explain the formation of these buried relic crevasses and the ice rise, they proposed a conceptual model where the floating ice shelf re-grounds on high points of the bedrock, leading to upstream thickening and downstream crevassing. Here, the ice–bed contacts are areas where the basal shear stress is higher than the nominally stress-free ice–sea water interface.

The complexity of basal crevasses arises from the heterogeneous condition of the ice–bed interface. Basal conditions of grounded ice sheets such as subglacial water pressure, basal shear stress and temperature vary spatially and temporarily. For some ice streams of the Antarctic Ice Sheet, rib-like patterns of high basal shear stress are inferred on the basis of inversions of surface data (Sergienko and Hindmarsh, Reference Sergienko and Hindmarsh2013). Other studies have argued that spatial and temporal variation of ice flow and surface velocity are a consequence of basal ‘sticky spots’ (referred to here as sticky patches) that have higher basal shear stress than their surroundings (Stokes and others, Reference Stokes, Clark, Lian and Tulaczyk2007). The existence of sticky patches has significant influence on the nonuniformity of ice-flow dynamics (Wolovick and others, Reference Wolovick, Creyts, Buck and Bell2014).

Here we argue that for grounded regions of glaciers and ice sheets, variation of basal-stress conditions can promote basal crevassing, while high basal water pressure remains the prerequisite. We quantify and explore the fracturing process induced by a combination of water pressure and sticky patches. This arrangement is different from the process of pure hydrofracturing, as considered by Smith (Reference Smith1976) and Van der Veen (Reference Van der Veen1998). At long timescales, glacial ice deforms as a viscous fluid, with a viscosity modelled by Glen's law (Glen, Reference Glen1955), which states that the viscosity has a power-law relation to the shear stress. The fracturing process, however, typically occurs on a shorter timescale over which ice behaves elastically. Therefore, linear elastic fracture mechanics (LEFM) theory is widely used to explore the existence and growth of cracks in ice. Van der Veen (Reference Van der Veen1998) applied LEFM to develop elastic models of surface and basal crevassing, with the assumption that basal crevasses are mode-I (i.e. opening) cracks propagating vertically in ice. More recently, Jimenez and Duddu (Reference Jimenez and Duddu2018) reassessed a key component of previous LEFM models of crevassing, finding that the basal boundary conditions influence the crevasses. For grounded ice sheets, Jimenez and Duddu (Reference Jimenez and Duddu2018) recommended a LEFM weight function from Tada and others (Reference Tada, Paris and Irwin2000) to study opening-mode fractures analytically.

The application of LEFM is not limited to opening-mode fractures, however. In several studies, LEFM is applied to calculating mixed-mode stress intensity factors (SIFs) in rift propagation on ice shelves in two dimensions (Hulbe and others, Reference Hulbe, LeDOUX and Cruikshank2010) and three dimensions (Lipovsky, Reference Lipovsky2020). Here, we use LEFM theory in the context of two-dimensional elasticity to explore the mixed-mode (in-plane opening mode and shearing modes) basal crevasses arising from sticky patches.

The paper is organised as follows. In the model section we introduce the mathematical theory and its numerical implementation. The results section illustrates the essential mechanics and shows the dependence of crack propagation on physical parameters. We first consider crack propagation under the mode-I fracture assumption followed by a consideration of curved fracture trajectories produced by mixed-mode fracturing. We find that the propagation of basal crevasses is controlled by basal water pressure, the size of the sticky patch and the stress variation between the sticky patch and the neighbouring bed. We then evaluate how these parameters affect the trajectories of basal crevasses. The discussion section explores the thermal consequences of basal crevasses, potential applications of the model to real ice sheets and limitations of the model.

2. The model

We model a grounded ice-sheet sliding due to gravity, influenced by a sticky patch. Under the assumption of plane strain (i.e. there is no cross-stream elastic strain), we model the ice sheet with a two-dimensional, infinite elastic strip sliding down a slope. Figure 1 shows a schematic diagram of the mathematical model, which considers a finite length of the strip that extends across a sticky patch. The total system is decomposed into two different sub-problems, based on the two components of the gravity vector in the tilted coordinate system. The first sub-problem addresses the sliding state, where steady sliding is driven by down-slope component of gravity g x and resisted by uniform basal stress (Fig. 1a). The second sub-problem is the sticky state, where the sticky patch provides additional basal drag, leading to tensile stress that partially offsets the static compression due to g z (Fig. 1b).

Fig. 1. Schematic diagram of the computational domain: a 2L × H strip extracted from the ice sheet on a slope with angle α. (a) The first sub-problem where ice is sliding at a constant rate under a uniform basal drag τ0. (b) The second sub-problem where there are a stick patch and a hydrofracture downstream. The sticky patch is represented by an interval [ − W, W] on the x-axis with excess basal stress Δτ. A water-filled crack is initiated on the downstream end of the sticky patch. The crack can be a vertical line with length Z C (as illustrated in the figure), or a curve, depending on the criterion used to determine its propagation.

In the first sub-problem, the ice strip with height H and length 2L (L ≫ H) is sliding at a constant speed. A uniform basal shear stress τ0, representing background basal drag, is imposed on (− L, L). The uniform basal stress τ0 exactly balances the down-slope component of gravity g x, leading to steady sliding with internal deformation (to maintain torque balance). This is consistent with the fact that the ice-sheet motion can be decomposed into basal sliding and internal deformation (Van der Veen, Reference Van der Veen2013).

In the sticky state of panel (b), we focus on the influence of the sticky patch by treating it in superposition to the steady, sliding state. To model the stickiness of the sticky patch, an excess shear stress Δτ is imposed on x ∈ [ − W, W] at z = 0. Since the down-slope component of gravity is already balanced by τ0 in panel (a), this excess traction needs to be balanced. We add an uniform increment of traction WΔτ/L in the opposite direction on x ∈ [ − L, L] at z = 0 – along the whole bottom boundary of the finite ice strip. This extra increment of uniform traction vanishes for W/L → 0; for W/L ≪ 1 it has little effect on basal stress estimates around the sticky patch.

In panel (b), the excess traction on − W to W at z = 0 creates tension on the downstream side of the sticky patch. If that tension, combined with basal water pressure, is large enough to overcome the ice overburden pressure plus the fracture toughness of ice, basal crevassing should occur. In natural ice sheets, the basal drag on the sticky patch Δτ can be much larger than the background basal drag τ0 (Sergienko and Hindmarsh, Reference Sergienko and Hindmarsh2013). Therefore, in regions close to the sticky patch, the elastic deformation caused by the uniform basal drag is negligible compared with the effect of the patch. From this point forward, we will only focus on the model in panel (b), which accounts for the effect of the sticky patch.

Based on LEFM theory, the SIFs are calculated and used to predict the maximum penetration of the crack (Broek, Reference Broek1982). SIFs are parameters that describe the stress distribution and magnitude close to the crack tip. A more detailed discussion of SIFs can be found in the following subsections.

2.1 2D elastic model of basal hydrofracture

Focusing on sub-problem (b) shown in Figure 1(b), the equation governing the conservation of momentum in ice is

(1)$${\boldsymbol \nabla}\cdot{\boldsymbol \sigma}-\rho_{\rm i} g_{z}\hat{{\boldsymbol z}} = {\bf 0},\; $$

where σ is the Cauchy stress tensor (tensile stress is positive), ρ i is the density of ice, g z = g cos α is the z-component of gravity and $\hat {{\boldsymbol z}}$ is the unit vector in positive-z direction. With the assumption that slope α ≪ 1, we take g z ≈ g in the following calculations. Both side boundaries of the ice strip are loaded by a depth-dependent compression ρ ig(z − H) that represents the overburden stress due to the weight of the overlying ice,

(2)$${\boldsymbol \sigma}\cdot{\boldsymbol n} = \rho_{{\rm i}}g \left(z-H\right){\boldsymbol n}\quad \rm{\,for}\ x = \pm L,\; $$

where n is the outward-pointing unit normal vector of the domain. The top boundary is assumed to be traction free, neglecting the atmospheric pressure and other traction on the surface,

(3)$${\boldsymbol \sigma}\cdot{\boldsymbol n} = {\bf 0} \quad \rm{\,for}\ z = H.$$

On the ice–bed interface (z = 0), we impose a zero-displacement boundary condition in the z-direction and the excess traction Δτ in x-direction,

(4)$${\boldsymbol u}\cdot{\boldsymbol n} = 0,\; $$
(5)$$\eqalign{{\boldsymbol t}\cdot{\boldsymbol \sigma}\cdot{\boldsymbol n} & = \tau\,\left(x,\; 0\right) = \left\{\matrix{ \displaystyle -{W\over L}\Delta\tau & \left\vert x\right\vert \geq W,\; \cr \displaystyle \left(1-{W\over L}\right)\Delta\tau & \left\vert x\right\vert < W,\; }\right.}$$

where Δτ represents the stress variation caused by a sticky patch, u is the displacement and ${\boldsymbol t} = -\hat {{\boldsymbol x}}$ is the unit tangent vector to the boundary. On the bottom, t is in the negative x-direction. The extra term − WΔτ/L is for force balance, as discussed above. Physically, we balance the excess shear stress on the sticky patch using a uniform stress on a much larger area. In our computation, the ratio W/L is set to be 0.1, which is small enough to make the effect of the extra term unimportant.

Crack walls are loaded by water pressure that is, in general, dependent on the subglacial hydrology. In this study, we assume static basal water pressure and represent subglacial hydrology in terms of the flotation fraction, which is the ratio of basal water pressure to basal overburden pressure,

(6)$$f = {\,p_{{\rm w}}\over \rho_{{\rm i}}g H} = {\rho_{{\rm w}}g H_{{\rm w}}\over \rho_{{\rm i}}g H},\; $$

where H w is the piezometric head measured in borehole experiments as p w = ρ wgH w (Harper and others, Reference Harper, Bradford, Humphrey and Meierbachtol2010). The fracture is assumed to initiate at x = W, where tensile stress is maximum, and reach a height Z C that is to be determined. Using Eqn (6), the boundary conditions on the crack walls (x = W, 0 ≤ z ≤ Z C) can be written in terms of f and z,

(7)$${\boldsymbol \sigma}\cdot{\boldsymbol n} = \left\{\matrix{ \displaystyle \rho_{{\rm i}}g H\left({\rho_{{\rm w}}\over \rho_{{\rm i}}}{z\over H} - f\right){\boldsymbol n} & \rm{at}\ {\it z} < {\it H}_{\rm w},\; \cr {\bf 0}\hfill & \rm{at}\ {\it z} \ge {\it H}_{\rm w}. }\right.$$

The transition at z = H w occurs because water rises to that height in the crevasse.

In formulating the plane–strain elastic constitutive relation, we assume that strain occurs in response to deviations from the overburden stress, as suggested by Cathles (Reference Cathles2015). This approach was used by Lipovsky (Reference Lipovsky2020) to model rift propagation in floating ice. A perturbation stress tensor T (sometimes referred to as resistive stress in glaciology), is introduced as the total Cauchy stress tensor σ minus the ice overburden stress,

(8)$${\boldsymbol T} = {\boldsymbol \sigma} + p_{{\rm i}}{\boldsymbol I},\; $$

in which p i = ρ ig(H − z) accounts for the ice overburden pressure and − p iI is the ice overburden stress. The elastic constitutive law linearly relates the strain to the perturbation stress as

(9)$${\boldsymbol T} = {E\nu\over \left(1 + \nu\right)\,\left(1-2\nu\right)} \rm{tr}\,\left({\boldsymbol \epsilon}\right)\,{\boldsymbol I} + {{\it E}\over 1 + \nu}{\boldsymbol \epsilon},\; $$

where

(10)$${\boldsymbol \epsilon}( \boldsymbol u) = {1\over 2}\left(\nabla \boldsymbol u + \nabla \boldsymbol u^{T} \right).$$

To define the constitutive law, two parameters are needed to account for ice properties. Here we use the Young's modulus E = 10 GPa, and a Poisson's ratio ν = 0.33 (Van der Veen, Reference Van der Veen1998). Further details associated with the problem description in terms of the perturbation stress are provided in Appendix A.

2.2 Stress intensity factors

In LEFM theory, SIFs K I, K II are used to describe the additional stress near the tip of a crack that is due to the presence of the crack (Broek, Reference Broek1982),

(11)$$\sigma_{i j}^{{I}}( r,\; \theta) = {K_{{I}}\over \sqrt{2\pi r}}f_{i j}^{{I}}\left(\theta\right) + {O}\left(r^{{1\over 2}}\right),\; $$
(12)$$\sigma_{i j}^{{II}}( r,\; \theta) = {K_{{II}}\over \sqrt{2\pi r}}f_{i j}^{{II}}\left(\theta\right) + {O}\left(r^{{1\over 2}}\right),\; $$

where I, II represent two in-plane modes of cracks (opening and shearing, respectively), r and θ are coordinates in a local plane-polar coordinate system with origin at the crack tip, and $f_{ij}^{{I}}$, $f_{ij}^{{II}}$ are tensor-valued functions describing the angular dependence of the excess stress. The O(r 1/2) terms are related to far-field stress and are omitted when r is small.

The predicted extent of the crack that is in equilibrium with the total stress field, including its direction (if not constrained to be vertical), is determined by the criterion of crack propagation. The broadly accepted G-criterion states that cracks grow in the direction along which the maximum potential energy is released (Broek, Reference Broek1982). In the context of this study, the potential energy is the strain energy stored in the elastic deformation. The elastic energy released when a crack extends is measured by the release rate G, defined as the energy released per unit of crack extension. The crack extends when G is greater than or equal to G C, a threshold value for crack growth. Equilibrium is attained when this condition is no longer met and the crack ceases to propagate.

G is related to the SIFs by

(13)$$G = G_{{I}} + G_{{II}} = {1-\nu^{2}\over E}\left(K_{{I}}^{2} + K_{{II}}^{2}\right),\; $$

where G I and G II denote the energy release rate for mode-I and mode-II components, with G C given by

(14)$$G_{{\rm C}} = {1-\nu^{2}\over E}K_{{I}, {\rm C}}^{2},\; $$

where K I,C is the fracture toughness, a material property measured in experiments. The G-criterion (G > G C) can be expressed in terms of the SIFs as $K = \sqrt {K_{{I}}^{2} + K_{{II}}^{2}}> K_{{I}, {\rm C}}$. Therefore, extension of a crack occurs when the total SIF exceeds fracture toughness K I,C. We assume that K I,C = 100 kPa m1/2 (Rist and others, Reference Rist1996), which is an estimate widely used in ice-fracture problems.

In the results section, we first assume a vertical, pure mode-I crack and employ a simple criterion for its propagation. In particular, we neglect the mode-II component by assuming that K II = 0. Thus the criterion for crack extension reduces to K I > K I,C; if satisfied, the crack can grow vertically in length. Later we reconsider mixed-mode crack growth with K II ≠ 0 and assess the difference associated with the different criteria.

2.3 Non-dimensionalisation

In the results section we focus on solutions of non-dimensional problems; we denote the non-dimensional version of a variable using the same symbol appended with a prime ′. We choose to scale length (and displacement) by H and stress by ρ igH. Hence, the conservation of momentum (Eqns (1)–(7)) and stress intensity factors are non-dimensionalised by the following scales:

(15)$$\boldsymbol x^{\prime} = {\boldsymbol x\over H},\; \quad \boldsymbol u^{\prime} = {\boldsymbol u\over H},\; \quad {\boldsymbol \sigma}^{\prime} = {{\boldsymbol \sigma}\over \rho_{{\rm i}}g H},\; \quad \Delta\tau^{\prime} = {\Delta\tau\over \rho_{\rm i} g H},\; \quad K_{{I}, {II}, {\rm C}}^{\prime} = {K_{{I}, {II}, {\rm C}}\over \rho_{\rm i} g H^{3/2}}.$$

Beside the flotation fraction representing basal water pressure, there are two important, non-dimensional parameters that control fracturing, W′ = W/H and Δτ′, which represent the size and excess shear stress of the sticky patch, respectively. The non-dimensional SIFs and non-dimensional fracture toughness are used in the fracture criterion. Note that although the dimensional fracture toughness is a material property measured by experiments, the non-dimensional fracture toughness is scaled by ρ igH 3/2, giving it a dependence on the ice-sheet thickness.

2.4 Numerical implementation

For the vertical-line-crack problem, the governing equations are solved using an open-source finite element (FE) software library, FEniCS (Logg and Wells, Reference Logg and Wells2010; Logg and others, Reference Logg, Mardal and Wells2012; Langtangen and Logg, Reference Langtangen and Logg2017), with meshes generated by Gmsh (Geuzaine and Remacle, Reference Geuzaine and Remacle2009). The basal crevasse is represented by a straight, triangular notch in the mesh, perpendicular to the bottom boundary. The mesh is locally refined near the tip of the notch and the bottom boundary. The mesh is composed of triangular elements, with element size varying from 5 × 10−3 times the crack length (near the tip) to 0.2 ice thicknesses (away from the tip). For an 1000 m-thick ice sheet with a 500 m-high basal crevasse, the element sizes are 2.5 m near the tip and 200 m away from the tip. SIFs are calculated by the displacement correlation method (DCM) (Chan and others, Reference Chan, Tuba and Wilson1970; Banks-Sills and Sherman, Reference Banks-Sills and Sherman1986) with Richardson extrapolation (Guinea and others, Reference Guinea, Planas and Elices2000). The code is benchmarked by comparison with a weight function from Jimenez and Duddu (Reference Jimenez and Duddu2018) (see Appendix B).

The FEniCS code relies on meshing the interior of the elastic domain, and requires local mesh refinement at the fracture tip to resolve the stress singularity there. To perform simulations of fracture growth where the path is not specified a priori, this method would require repeated remeshing. Dynamic remeshing is feasible with specialised software but is unreliable and incurs a significant computational cost. To avoid such issues and to provide an additional means to verify the numerical results, we use a separate code for consideration of curved fracture trajectories. This code employs the displacement discontinuity method (DDM) (Crouch and Starfield, Reference Crouch and Starfield1983), a scheme that is based on the boundary element method (BEM).

In the DDM approach, only boundaries subject to given conditions must be meshed. The interior of the domain is assumed to be a uniform, linear elastic material that deforms according to Green's functions forced by the facets of the boundary mesh. This method avoids reliance on meshing software. To extend the fracture at a given step, a new straight-line segment is connected to the former tip of the fracture. The orientation of this new segment is defined by the optimal direction of fracture growth. The new segment must have an appropriately small length, the choice of which defines the resolution of the model. At each growth step, we test for multiple orientations of growth and choose that with the highest strain-energy release rate, as defined by the G-criterion (Dahm, Reference Dahm2000). The BEM code used here is developed by Davis (Reference Davis2017). As with the FE models, in the DDM models we subject the base of the glacier domain to the condition of no vertical displacement.

3. Results

Results are presented in two parts. The first is for vertical fractures and uses the FE model; the second is for mixed-mode, curving fractures and uses the DDM model.

3.1 Mode-I fracture growth

We consider a vertical, mode-I fracture as part of a reference case where the sticky patch has a width two times the ice thickness (W′ = 1). The non-dimensional excess shear stress is Δτ′ = 0.3, which means the dimensional excess shear stress is 0.3ρ igH. The reference flotation fraction is f = 0.7, which we think is a reasonable estimate (Engelhardt and others, Reference Engelhardt, Humphrey, Kamb and Fahnestock1990; Harper and others, Reference Harper, Bradford, Humphrey and Meierbachtol2010; Andrews and others, Reference Andrews2014). Cases with larger f (0.7 ≤  f ≤ 1) are also included in the following results. Figure 2 shows the three components of the perturbation stress that arise due to the sticky-patch excess stress and an existing basal crevasse. Among these, we are most interested in $T_{xx}^{\prime }$ and $T_{xz}^{\prime }$, which are related to the fracture propagation. As shown in panel (a), there is horizontal tension concentrated on the downstream (right) end of the sticky patch and compression concentrated on the upstream (left) end. In panel (b), the vertical normal stress $T_{zz}^{\prime }$ is also concentrated near the ends of the sticky patch and the crack tip. This is a consequence of the no-vertical-displacement condition Eqn (4) imposed on the bottom boundary. Panel (c) shows the pattern of shear stress, which is localised to a region around the sticky patch. In this case, the stress intensity is K I/K I,C = 8.75. Thus the crack is unstable and will continue growing vertically.

Fig. 2. Perturbation stress in a cracked ice sheet with W′ = 1, Δτ′ = 0.3, and $Z_{\rm C}^{\prime } = 0.3$. Panels show (a) $T_{xx}^{\prime }$, (b) $T_{zz}^{\prime }$ and (c) $T_{xz}^{\prime }$. Tension is defined as having a positive sign. In the zoom-in box a white rectangle is added to highlight the crack, which is actually much narrower, making it difficult to see otherwise.

We are interested in the maximum stable crack length permitted by the fracture criterion, and how this depends on ice-sheet thickness and sticky-patch width. To investigate, we extend our calculation to cases with varying crack length $Z_{{\rm C}}^{\prime }$ and assume that the criterion for crack growth is $K_{{I}}^{\prime }> K_{{I}, {\rm C}}^{\prime }$ (Van der Veen, Reference Van der Veen1998). Figure 3 shows the normalised SIF K I/K I,C as a function of $Z_{{\rm C}}^{\prime }$ near the sticky patch for the reference case (W′ = 1, Δτ′ = 0.3) and another case with smaller Δτ′ = 0.2. Comparison of the two cases shows that increasing Δτ′ leads to larger K I. In each case, with the crack length $Z_{{\rm C}}^{\prime }$ increasing, $K_{{I}}^{\prime }$ increases to its maximum due to tension near the sticky patch, then begins to drop and eventually reaches negative values due to overburden pressure. The maximum crack length is the value of $Z_{{\rm C}}^{\prime }$ such that $K_{{I}}^{\prime } = K_{{I}, {\rm C}}^{\prime }$ (i.e. the intersection of the $K_{{I}}^{\prime }$ curve and vertical $K_{{I}, {\rm C}}^{\prime }$ line in Fig. 3). For the crack with this value of $Z_{{\rm C}}^{\prime }$ to be in a stable equilibrium, the crack should be resistive to perturbations. If we add a positive perturbation $\Delta Z_{{\rm C}}^{\prime }> 0$ to the crack length $Z_{{\rm C}}^{\prime }$, the perturbed crack length should return to the stable state, which means $K_{{I}}^{\prime }( Z_{{\rm C}}^{\prime } + \Delta Z_{{\rm C}}^{\prime }) < K_{{I}, {\rm C}}^{\prime }$. Thus, the condition for a stable crack is ${{\rm d} K_{{I}}}^{\prime }/{{\rm d} Z_{{\rm C}}^{\prime }}< 0$. Note there is K I < K I,C at small crack lengths. To propagate a basal crevasse, an initial flaw or crack of adequate length and orientation is required (Van der Veen, Reference Van der Veen1998). Here we simply assume a pre-existing vertical line crack with K I > K I,C.

Fig. 3. Normalised stress intensity factor KI/KI,C versus dimensionless crack length with W′ = 1. Two cases with different Δτ′ are considered: Δτ′ = 0.2 (red) and Δτ′ = 0.3 (blue). Vertical axis represents the dimensionless crack length. Horizontal axis represents the dimensionless SIFs normalised by the fracture toughness, $K_{{I}, {\rm C}}^{\prime } = K_{{I}, {\rm C}}/( \rho _{\rm i} g H^{3/2})$, where H is chosen to be 100 m. The black dash-dotted line represents the fracture toughness.

Figure 4 shows the maximum stable crack length $Z_{{\rm C, max}}^{\prime }$ as a function of dimensionless problem parameters Δτ′ and W′ for an ice thickness of 100 m. In panel (a), flotation fraction is f = 0.7. The maximum crack length increases monotonically with W′ or Δτ′. For sticky patches with W′ = 1, a minimum Δτ′ ≈ 0.2 is required to have a crack with length $Z_{{\rm C}}^{\prime }\sim 0.2$. Panel (b) shows another case where the ice is closer to flotation (i.e. f = 0.9). Here, fracture can be propagated by smaller values of W′ and Δτ′, compared with panel (a). The sticky patch creates a localised horizontal tension that promotes the growth of the vertical line crack, as shown by the reference case in Figure 2. Evidently, this effect is sensitive to the non-dimensional parameters W′ and Δτ′.

Fig. 4. Maximum crack length $Z_{{\rm C, max}}^{\prime }$ as a function of W′ and Δτ′. Note that in the criterion, the non-dimensional fracture toughness depends on the ice thickness H. Here we set $H = 100\ \rm {m}$. For larger thicknesses, $K_{{I}, {\rm C}}^{\prime }$ gets smaller, leading to longer cracks. Two cases with different f are considered: (a) f = 0.7, (b) f = 0.9. When the crack length $Z_{{\rm C}}^{\prime } = 1$, the crack dissects the full ice-sheet thickness.

3.2 Mixed-mode fracture growth

A limitation of the models above is that they assume basal crevasses are mode-I cracks that only propagate vertically and, to simplify the calculation, neglect the mode-II component. This approach is consistent with previous LEFM models of basal crevasses that are pure hydrofractures (Van der Veen, Reference Van der Veen1998). However, with excess shear stress arising from the sticky patch, basal crevasses are expected to have both mode-I and mode-II contributions. Relaxing our assumption that the crack propagates vertically, we now consider the curved, quasi-static path of fractures using the BEM implementation. Note that for curved fracture paths, we use the G-criterion for fracture growth: the fracture stops when either G < G C or when there is closure of fracture walls just behind the fracture's tip (K I ≤ 0).

For an ice thickness H = 1000 m, we impose an initial vertical line crack with length 30 m and propagate the crack in its curved trajectories determined by the above criterion. Figure 5 shows the mixed-mode fracture paths with different values of dimensionless excess shear stress Δτ′ and flotation fraction f. The three panels represent three sizes of the sticky patch: W′ = 0.1 for panel (a), W′ = 1.0 for panel (b) and W′ = 10.0 for panel (c).

Fig. 5. Fracture paths calculated under different excess basal shear stress Δτ′ and flotation fraction f, with three values of W′: (a) W′ = 0.1; (b) W′ = 1.0; (c) W′ = 10.0. An ice thickness H = 103 m is used to scale the fracture toughness. The background vector field indicates the local direction of the maximum compressive stress, with the vector length scaled by the deviatoric stress $\left \vert {T_{1}^{\prime }-T_{2}^{\prime }}\right \vert$, where $T_{1}^{\prime }$ and $T_{2}^{\prime }$ are the local principal stresses calculated from T′. In panels (a) and (b), the fracture paths (colours) are calculated by BEM, with the principal stress trajectory represented by the dashed curve. Colours represent four different combinations of Δτ′ and f as in the legend. In panel (c), since the principal stress trajectory is nearly a vertical line, we no longer conduct the BEM simulation and just assume purely vertical fracture paths. Thus, the maximum crack length is determined as we have done in Figure 3. Colours represent three different combinations of Δτ′ and f as in the legend.

In panel (a) we consider a sticky patch whose width is one-fifth of the ice thickness (W′ = 0.1). There is a zoom-in plot of the fracture paths in the black box. The four curves are quasi-static paths with four combinations of flotation fraction f and dimensionless excess shear stress. The paths deviate from the vertical path assumed in models above. In particular, they incline upstream under the influence of the shear stress. Comparing the four curves in panel (a), we find that Δτ′ and f have little effect on the direction of the paths.

In panel (b), the length of the panel is kept the same as panel (a). In addition, we keep the same combinations of f and Δτ′. With W′ = 1 (sticky patch width twice the ice thickness; an order of magnitude larger than in panel (a)), the fractures align closer to the vertical and extend further into the ice sheet.

Panel (c) shows the fracture paths when W′ = 10. In this case the sticky patch is 20 times wider than the ice thickness. The magnitude of the excess shear stress Δτ′ required for fracture propagation is reduced to about 0.035. Meanwhile, since $T_{xx}^{\prime }\gg T_{zz}^{\prime }$, the principal stress trajectories are nearly vertical lines, as indicated by the background stress field. Thus, fracture paths can be approximated by vertical lines, as we did previously. The crack length becomes very sensitive to the magnitude of Δτ′, since a small perturbation to Δτ′ (from 0.033 to 0.036) causes a large variation of the crack length.

It is possible to predict the trajectory of the curved cracks without solving the BEM model, using only stresses computed in an uncracked domain. This may be computationally convenient in combination with Stokes-flow models of ice sheets (Krug and others, Reference Krug, Weiss, Gagliardini and Durand2014; Yu and others, Reference Yu, Rignot, Morlighem and Seroussi2017). As in the models discussed above, the excess shear stress due to the sticky patch is imposed as a boundary condition on the uncracked domain. We then calculate the perturbation stress in the uncracked domain and plot the principal stress trajectories as the dashed curves in Figure 5. We use the perturbation stress due to the excess shear only, because the overburden stress and water pressure are isotropic and do not contribute to the orientation of the principal stresses. Comparison between these trajectories and the fracture paths predicted by the BEM shows that the former are accurate approximations of fracture paths under the conditions considered here.

If there is no deviatoric stress from sources other than the sticky patch, the magnitude of Δτ′ does not contribute to the direction of the fracture path. Instead, the direction of trajectories depends only on the ratio W′. This is consistent with the results predicted by the BEM and indicates that for real basal crevasses affected by sticky patches, the direction of the fracture is predominantly controlled by the relative size of the sticky patch. Since the area of sticky patches can be of the order of 1–100 km2 (Stokes and others, Reference Stokes, Clark, Lian and Tulaczyk2007), we believe the panels in Figure 5 are representative cases that reveal the influence of patch size. If sticky patches are usually W′ ≫ 1 in real ice sheet, then crevasses are mostly vertical and BEM would not be necessary for future studies.

4. Discussion

We have developed and analysed a model of basal crevassing associated with sticky patches at the bed of an elastic glacier or ice sheet. Our model, based on LEFM theory, evaluates the role of shear-stress variations and makes predictions of crack lengths and trajectories. As shown above, the growth of such basal cracks depends on the flotation fraction f, the non-dimensional size of the sticky patch W′ and the non-dimensional stress variation Δτ′.

4.1 Initiation of basal crevasses

When modelling both vertical and curved basal crevasses in the results section, we assumed pre-existing small cracks with K I > K I,C. They are a prerequisite for further fracture propagation. By using a simple Griffith fracture model (Tada and others, Reference Tada, Paris and Irwin2000)

(16)$$K_{{I}} = K_{{I}, {\rm C}} = \sigma_{xx}\sqrt{\pi Z_{{\rm C}}},\; $$

we estimate the minimum crack length Z C required for growth, given a local effective stress σxx across the fracture. In the BEM simulations we assumed a 30 m-long initial crack, corresponding to a 10 kPa local stress; this combination guarantees K I > K I,C. If tensile stress reaches 326 kPa, basal fractures might initiate on the boundaries of grains 3 cm in diameter. This large value of stress seems unlikely to be achieved under normal circumstances.

4.2 Potential applications to real ice sheets

For real ice sheets, basal shear stress patterns are difficult to observe directly. In several cases, they have been estimated by inversion of surface data under specific assumptions. Sergienko and Hindmarsh (Reference Sergienko and Hindmarsh2013) showed that for a region of the Antarctic ice sheet with large ice thickness (H ≈ 103 m), the basal shear stress has rib-like patterns of variation. The width of such sticky patches varies from one ice thickness (W′ = 0.5) to 10 ice thicknesses (W′ = 5). The excess shear stress (200–300 kPa) estimated for these locations is small compared with the overburden pressure (Δτ′ ~ 0.03). In this case, the effect of the sticky patch depends on its size and local water pressure, as discussed in the Results section. For small patches (W′ ~ 1), such excess shear stress only affects the direction of fracture propagation, as shown in Figure 5. Basal crevassing occurs when the flotation fraction f reaches 1, which is consistent with the conclusions of Van der Veen (Reference Van der Veen1998). For larger patches with W′ ~ 10, basal crevasses would initiate on the downstream end at a lower water pressure f ~ 0.7.

For alpine glaciers with an ice thickness of order 100 m, basal shear stress was found to be in the range of 0–200 kPa (Brædstrup and others, Reference Brædstrup, Egholm, Ugelvig and Pedersen2016). In the non-dimensional parameter space, Δτ′ is of order 0.1, larger than in Antarctic settings. Moreover, W′ is of the order of 1–10, which is similar to that of sticky patches in Antarctic ice sheets. Thus, we predict that sticky patches play a more important role in determining the length and direction of basal crevasses in alpine glaciers. For further investigations of the basal hydrofractures around sticky patches, we need to go beyond the inversion results, understand the specific causes of the sticky patches and include the essential physics in our model.

In natural glaciers and ice sheets, there are many factors that might create sticky patches, including bedrock bumps, till-free areas, well-drained tills and basal freeze-on. All of these would lead to localised, high basal friction (Stokes and others, Reference Stokes, Clark, Lian and Tulaczyk2007). An important friction phenomenon is the stick-slip motion of ice, detected in Whillans Ice Stream (WIS) in West Antarctica. Wiens and others (Reference Wiens, Anandakrishnan, Winberry and King2008) investigated the WIS stick-slip motion and related it to a sticky patch on the bed. Furthermore, Sergienko and others (Reference Sergienko, MacAyeal and Bindschadler2009) argued that WIS can be considered as a typical stick-slip system controlled by basal friction, where the sticky spot nucleates the stick-slip cycle. Therefore, for sticky spots associated with stick-slip ice motion, we can estimate the stress variation from parameters in relevant friction experiments (McCarthy and others, Reference McCarthy, Savage and Nettles2017; Lipovsky and others, Reference Lipovsky2019; Zoet and others, Reference Zoet2020), rather than from inversion from surface data based on viscous rheology. If we interpret Δτ in terms of friction coefficient variation Δμ, we find that

(17)$$\Delta\tau^{\prime} = \Delta\mu N^{\prime},\; $$

where N′ = 1 − f is the dimensionless effective normal stress. The basal stress variation is controlled by both friction coefficient variation Δμ and the flotation fraction f. The magnitude of Δμ during the stick-slip motion of ice can be measured experimentally. By keeping a constant normal stress of N = 500 kPa between the ice sample and the bedrock asperity, Zoet and others (Reference Zoet2020) found that the friction-coefficient decrease during stick-slip motion is between 0.1 and 0.4. We assume that the measurements of Δμ also apply to the sticky patch discussed in our model.

Well-drained till could serve as a sticky patch. For a well-drained till surrounded by a water saturated layer (Stokes and others, Reference Stokes, Clark, Lian and Tulaczyk2007), the water pressure on the till would be smaller than the surroundings. For such a till, we assume that the excess friction coefficient Δμ is 0.4 and flotation fraction f is 0.7. Then, according to (17), the excess basal shear stress is Δτ′ ~ 0.1, which would lead to a tensile-stress concentration on the downstream end of the patch. Meanwhile, in the surroundings, the local subglacial water pressure is expected to be higher (f ≥ 0.7). For the case considered above, when f reaches 0.9, basal crevassing is likely to occur.

4.3 Thermal implications of basal crevasses

The temperature structure of ice has important effects on ice dynamics. Using high-vertical-resolution sensing, Law and others (Reference Law2021) reported spatial heterogeneity of englacial ice temperature and deformation. They found a basal temperate-ice zone with thickness that varies from 5 to 73 m at two locations separated by only 9 km at Store Glacier, an outlet glacier of the Greenland Ice Sheet. Their study indicates spatially varying basal thermal conditions over distances of a few ice thicknesses. Injection of water-filled basal crevasses can locally modify the thermal profile of the ice sheet (Luckman and others, Reference Luckman2012) and is a potential explanation of the spatially-varying temperate ice layer. The thermal structure of basal crevassing, which is similar to that of dykes in rock (Daniels and others, Reference Daniels, Bastow, Keir, Sparks and Menand2014), has been modelled in several studies (Jarvis and Clarke, Reference Jarvis and Clarke1974; McDowell and others, Reference McDowell, Humphrey, Harper and Meierbachtol2021). Their approach recognises that water-filled basal crevasses in sub-temperate ice propagate on a short timescale, followed by rapid refreezing of water inside the crack. McDowell and others (Reference McDowell, Humphrey, Harper and Meierbachtol2021) modelled the refreezing of a basal crevasse as an instantaneous heat source in a one-dimensional heat-conduction system, using the analytical solution of Carslaw and Jaeger (Reference Carslaw and Jaeger1959). We use the same analytical solutions for a two-dimensional thermal structure of a basal crevasse. Details of this calculation are provided in Appendix C.

To estimate the heat released by refreezing, we return to the dimensional problem and assume a static, water-filled vertical crack with crack length Z C = 50 m in a 100 m-thick, subtemperate ice sheet. The crack width w is assumed to be 10 cm uniformly along the crack. The water inside the crack instantaneously refreezes at t = 0, releasing an amount of heat q i per unit length per unit depth into the page (i.e. in the direction normal to the crack). The heat is estimated as

(18)$$q_{{\rm i}} = \rho_{{\rm i}} w L = 3\times 10^{7}\ \rm{J}\, \rm{m}^{-2},\; $$

where ρ iw is the mass of water per unit area of the fracture and L is the latent heat of solidification. Mathematically, the refreezing process is assumed to be an instantaneous source at t = 0. The temperature rise caused by refreezing, ΔT, is held at 0 K on the surface boundary with the atmosphere, the basal boundary and at the limit of x → ±∞. The temperature rise will asymptotically decay to zero after a long cooling process.

Figure 6 shows the perturbation in temperature due to the refreezing of water in a single basal crevasse over 20 years. The surrounding ice undergoes a rapid warming at t = 0, followed by a long cooling period until the temperature drops back to the background state ΔT = 0 (McDowell and others, Reference McDowell, Humphrey, Harper and Meierbachtol2021). Panel (a) shows that after 5 years of diffusion, the temperature perturbation due to refreezing is localised around the crevasse and decreases sharply within ± 30 m. Panels (b) and (c) show the temperature perturbation after 10, 15 and 20 years. After 20 years, ΔT drops back below 0.05 K and the system is again close to the background state.

Fig. 6. Perturbation of a single basal crevasse on background temperature field. t = 0 is the time of crevasse opening and refreezing. (a) t = 5 years. (b) t = 10 years. (c) t = 15 years. (d) t = 20 years.

If a sticky patch is fixed on the bedrock beneath a sliding ice sheet, the patch could generate a series of basal crevasses, as shown in Figure 9 (Appendix C). Thus, we next consider a case in which the spacing between these crevasses is equal to the half-width of the sticky patch W = 100 m. For these equally spaced basal crevasses, the temperature field is a linear superposition of the effect of each crevasse. In the mathematical model, the basal crevasse is initiated on the downstream end of the sticky patch and subsequently advected by sliding. Details are provided in Appendix C.

The thermal effect of a series of basal crevasses is shown in Figure 7. The downstream perturbation will gradually smooth out after decades of cooling. The stable pattern of temperature rise depends on the volume of water inside the crack. The BEM simulations show that the thickness of basal crevasses is of order 0.1 m, which is much smaller than the upper limit of basal crevasse thickness observed by Harper and others (Reference Harper, Bradford, Humphrey and Meierbachtol2010). Here we simply assume that the crack width is 0.1 m. The initial, localised temperature rise will decay rapidly as it smooths out. Refreezing in basal crevasses is a possible factor influencing the temperature profile of basal ice, leading to localised heating around the relic basal crevasses, followed by cooling back toward a steady state (McDowell and others, Reference McDowell, Humphrey, Harper and Meierbachtol2021).

Fig. 7. Thermal effect of a series of basal crevasses that is 2–3 km downstream from the sticky patch. After a long time of diffusion, the localised temperature field smooths out.

Alternatively, if the basal fracture is embedded in temperate ice and hence doesn't freeze rapidly, it becomes a persistent mechanical perturbation to the ice sheet. A series of such basal crevasses will be carried downstream to the grounding line. If an ice shelf extends seaward from the grounding line, the basal crevasses will weaken the shelf to lateral shear stresses associated with buttressing and vertical shear associated with calving (Bassis and Ma, Reference Bassis and Ma2015). The recent numerical model developed by Berg and Bassis (Reference Berg and Bassis2022) suggests that the advection of crevasses could increase the calving rate and promote glacier retreat.

4.4 Limitations of the model

The model presented has three significant limitations. First, it is based on equilibrium equations of elasticity. The ice sheet is assumed to be an isotropic, elastic body without any internal viscous deformation that is commonly computed by a Stokes model. Therefore, our model only accounts for the physics of fracturing on a short timescale in which ice is dominated by elastic rheology. This approach might miss important interactions between the sticky patch and viscous deformation that would modify the stress field. Basal crevasses can be deformed by the internal deformation caused by viscous flow. For an ice sheet with a thickness H = 1000 m, the velocity difference between the surface and the bed can be up to several kilometres per year, which can significantly deform the basal crevasses whose length is of the same order as the ice thickness. Thus, a viscoelastic rheology is required to study basal crevasses on a longer timescale. Moreover, it is important to include the spatial and temporal variation of subglacial hydrology (Harper and others, Reference Harper, Humphrey, Pfeffer, Fudge and O'Neel2005; Hewitt, Reference Hewitt2013), which is simplified to a static water pressure in the current model.

The second main limitation is that the model does not account for the three-dimensional effects of sticky patches. In a real ice sheet, some sticky patches may have a round shape instead of a long, rib-like pattern. In that case, it is more appropriate to study them in a 3-D domain, which includes both vertical and lateral extension of cracks, rather than as a plane-strain problem in the x–z plane.

The third limitation relates to the fracture toughness. In applying the LEFM approach, we implicitly assume a brittle medium. However, a brittle-to-ductile transition occurs at a critical grain size (Schulson and others, Reference Schulson, Lim and Lee1984). Regions and layers of coarse-grained ice exist in ice sheets (Gow and others, Reference Gow1997). Thus, in some glaciers and ice sheets there may exist ductile ice that eludes brittle fracture. This heterogeneity would modify the predicted crack lengths. Furthermore, Sinha (Reference Sinha1978) argued that the Young's modulus of ice E depends on temperature, even on a short loading timescale. Further investigation is therefore required to determine whether the thermal profile is important for ice elasticity and fracture.

A more complete understanding of how basal crevasses grow and interact with viscous flow will require a three-dimensional, viscoelastic model including variations in fracture toughness.

5. Conclusions

Besides basal water pressure (Van der Veen, Reference Van der Veen1998), stress variations of sticky patches at the ice–bed interface can promote the propagation of basal crevasses. In previous studies, basal crevasses were assumed to be mode-I hydrofractures under pure horizontal tension (Van der Veen, Reference Van der Veen1998; Krug and others, Reference Krug, Weiss, Gagliardini and Durand2014; Jimenez and Duddu, Reference Jimenez and Duddu2018). Assuming water pressure smaller than the flotation condition, we examined the effect of sticky patches on basal crevassing in a grounded glacier or ice sheet. We found that sticky patches can provide stress required for crack extension. Alongside the flotation fraction, such sticky-patch-assisted crevassing depends on two non-dimensional parameters: (1) the ratio W′ of the sticky-patch half-width to the ice thickness, and (2) the ratio Δτ′ of excess shear stress to the basal ice overburden pressure.

With a sufficient variation of basal shear stress, the direction of basal fracture is controlled by the relative size of the sticky patch. When the width of the sticky patch is much larger than the ice thickness, basal crevasses grow nearly vertically and are essentially mode-I fractures. When the width of the sticky patch is smaller than the ice thickness, however, curved basal crevasses grow with trajectories inclined upstream. In this case, principal stresses can be used to approximate crack trajectories.

For real glaciers or ice sheets with complicated geometries, our model can be combined with the basal stress pattern to investigate how stress variation promotes basal crevassing. Compared to stress estimates from inversion calculations (Sergienko and Hindmarsh, Reference Sergienko and Hindmarsh2013; Brædstrup and others, Reference Brædstrup, Egholm, Ugelvig and Pedersen2016), qualitative analysis shows that for a flotation fraction of about 0.9, the basal-stress pattern in alpine glaciers and ice sheets may play an important role in determining the growth of basal crevasses. To better understand the basal crevasses controlled by sticky patches, future research could incorporate viscous ice flow, spatially resolved subglacial hydrology and detailed fracturing properties of ice.

Acknowledgements

We thank B. Lipovsky and C.-Y. Lai for their guidance on the application of LEFM to ice, R. Greve and A. Rempel for their insightful editorial comments, and two anonymous reviewers for helpful suggestions. The general idea that we explore arose during a conversation with P. Christoffersen and R. Law. This research received funding from the European Research Council under Horizon 2020 research and innovation programme grant agreement number 772255.

Code availability

We use the stable version of FEniCS (Logg and Wells, Reference Logg and Wells2010; Logg and others, Reference Logg, Mardal and Wells2012; Langtangen and Logg, Reference Langtangen and Logg2017) (https://quay.io/repository/fenicsproject/stable) and the open-source boundary element code CutAndDisplace (https://github.com/Timmmdavis/CutAndDisplace.git) developed by Davis (Reference Davis2017). The code for the runs in this study is available at the online repository (https://doi.org/10.5281/zenodo.6364838).

Appendix A. 2-D elasticity in terms of perturbation stress

The domain is a notched ice strip with length 2L as shown in Figure 1. We set L ≫ W to make sure the crack is far from the boundaries to avoid any edge effects. Substituting T = σ + p iI into the governing equation and boundary conditions, we can express the stress equilibrium equation in terms of the perturbation stress as

(A.1)$${\boldsymbol \nabla}\cdot{\boldsymbol T} = \bf 0.$$

For the top boundary, a traction-free boundary condition is imposed as

(A.2)$${\boldsymbol T}\cdot{\boldsymbol n} = {\bf 0} \quad \rm{\,for}\ {\it z} = {\it H},\; $$

where n is the outward-pointing unit normal vector of the domain. On both sides, we impose traction-free boundary condition for the perturbation stress,

(A.3)$${\boldsymbol T}\cdot{\boldsymbol n} = {\bf 0} \quad\rm{\,for}\ {\it x} = \pm {\it L}.$$

On the crack walls (x = W, 0 ≤ z ≤ Z C), compression caused by static water pressure is imposed as

(A.4)$${\boldsymbol T}\cdot{\boldsymbol n} = \left\{\matrix{ \displaystyle\rho_{{\rm i}}g H\left[\left({\rho_{{\rm w}}\over \rho_{{\rm i}}}-1\right){z\over H} + 1 - f\right]{\boldsymbol n} & \rm{\,for}\ {\it z} < {\it H}_{\rm w},\; \cr \displaystyle \rho_{{\rm i}}g H\left(-{z\over H} + 1\right){\boldsymbol n} & \rm{\,for}\ {\it z}\ge {\it H}_{\rm w},\; }\right.$$

where ρ ifH/ρ w = H w is the hydraulic head. The bottom boundary condition remains the same as (5) and (4), and since the overburden stress doesn't contribute to the shear stress and the elastic deformation,

(A.5)$$\eqalign{{\boldsymbol t}\cdot{\boldsymbol T}\cdot{\boldsymbol n} & = \tau\,\left(x,\; 0\right) = \left\{\matrix{ \displaystyle-{W\over L}\Delta\tau & \left\vert x\right\vert \geq W,\; \cr \displaystyle\left(1-{W\over L}\right)\Delta\tau & \left\vert x\right\vert < W,\; }\right.}$$
(A.6)$${\boldsymbol u}\cdot{\boldsymbol n} = 0.$$

Here ${\boldsymbol t} = -\hat {{\boldsymbol x}}$ is the unit tangent vector to the bottom boundary. Note that in Eqn (A.5) an extra term − WΔτ/L is added on the bottom boundary in order to maintain the total force balance of the ice strip. Because L ≫ W, we can neglect the near-field effect of this force-balance term.

Different from standard elasticity, we set up a constitutive relation between the perturbation stress tensor T and the strain tensor ε as in Eqn (9).

Appendix B. Benchmark

The weight function method has proved to be a useful tool to calculate stress intensity factors in ice sheets under certain basal boundary conditions (Tada and others, Reference Tada, Paris and Irwin2000; Jimenez and Duddu, Reference Jimenez and Duddu2018). In essence, this involves calculating stress intensity factors by integrating the stress along a hypothesised crack in an uncracked domain with an appropriate weight function. The advantage of the weight function method is that the computation is performed in the uncracked domain with no discontinuity and singularity caused by cracks. For certain simple domain and crack geometry and boundary conditions, the weight function method can give SIFs with high accuracy.

Jimenez and Duddu (Reference Jimenez and Duddu2018) has provided the appropriate weight function from Tada and others (Reference Tada, Paris and Irwin2000) to be used for basal cracks in a grounded ice sheet,

(B.1)$$K_{{I}} = \int_{0}^{Z_{C}}\sigma_{xx}\, \left(z \right)\,G_1\left(\lambda,\; \gamma\right)\,\rm{d}{\it z},\; $$
(B.2)$$G_1\left(\lambda,\; \gamma\right) = {2\over \sqrt{2H}}\sqrt{{\tan\left({\pi\lambda \over 2}\right)\over 1-{\cos{\pi\lambda \over 2}}/{\cos{\pi\lambda\gamma \over 2}}}} \,\left\{1 + 0.297\sqrt{1-\gamma^{2}}\,\left[1-\cos\left({\pi\over 2}\lambda\right)\right]\right\},\; $$

where λ = Z C/H is the non-dimensional crack length and γ = z/Z C is the z coordinate normalised by the crack length. This weight function can be used to predict the mode-I stress intensity factor K I in grounded ice. A similar weight function in Tada and others (Reference Tada, Paris and Irwin2000) can also be applied to calculating mode-II stress intensity factor,

(B.3)$$K_{{II}} = \int_{0}^{Z_{{\rm C}}}\sigma_{xz} \left(z \right)\,G_2\left(\lambda,\; \gamma\right)\rm{d}{\it z},\; $$
(B.4)$$G_2\left(\lambda,\; \gamma\right) = {2\over \sqrt{2H}}\sqrt{{\tan\left({\pi\lambda \over 2}\right)\over 1-{\cos{\pi\lambda \over 2}}/{\cos{\pi\lambda\gamma \over 2}}}} \,\left\{1 + 0.297\sqrt{1-\gamma^{2}}\,\left[1-\cos\left({\pi\over 2}\lambda\right)\right]\right\}{\sin{{\pi \lambda\gamma \over 2}}\over \sin{{\pi \lambda \over 2}}}.$$

To verify the implementation of the DCM method we consider a problem where there is a mixed-mode (mode-I and mode-II) water-filled crack. SIFs are calculated by both the DCM method and the weight function method in Figure 8. In the DCM calculation, the element size is set to be 0.005 times the crack length near the crack tip and 0.2 ice thicknesses away from the tip.

Fig. 8. Stress intensity factors K I and K II when Δτ′ = 0.3 and f = 0.7. The red lines show SIFs calculated by the weight function method (‘WF’ in the legend means ‘weight function’). The blue dots show SIFs calculated by DCM.

Fig. 9. Schematic of a series of basal crevasses produced on a sticky patch.

For K I, these two methods agree over most crack lengths. However, there is large deviation for K II at large crack length.

Appendix C. Thermal structure of basal crevasses

C.1 Thermal structure of a single basal crevasse

Refreezing of water inside basal crevasses is a factor potentially affecting the thermal profile of ice. Using analytical solutions of Carslaw and Jaeger (Reference Carslaw and Jaeger1959), this section shows how a single basal crevasse or a series of basal crevasses affects the thermal structure of ice.

In order to simplify the computation, we neglect advection and other englacial heat sources and focus on heat diffusion after refreezing (Luckman and others, Reference Luckman2012). The background temperature profile is assumed to have no effect on the heat conduction caused by refreezing. That means that the simple analytical model just accounts for the effect of refreezing and doesn't make any prediction of the net temperature profile. Let ΔT be the temperature perturbation induced by a basal crevasse. The ice sheet is simplified to an isotropic infinite strip 0 < z < H. On the boundaries we assume zero temperature perturbation (ΔT = 0) at z = 0, z = H and x = ±∞. The governing equation of thermal conduction in terms of ΔT is

(C.1)$${\partial^{2} \Delta T\over \partial x^{2}} + {\partial^{2} \Delta T\over \partial z^{2}} = {1\over \kappa}{\partial \Delta T\over \partial t}.$$

The boundary conditions are

(C.2)$$\Delta T\left(\pm \infty,\; z\right) = 0,\; $$
(C.3)$$\Delta T\left(x,\; 0\right) = \Delta T\left(x,\; H\right) = 0.$$

Table 1 shows the values of parameters used in the calculation. At time t = 0, an amount of heat q i per unit length per unit depth into the page is released by a line heat source from (0, 0) to (0, Z c). In order to get an estimated value of that released heat, we need to estimate the volume of water that refreezes at t = 0. Based on the BEM simulations, the width of the crack w is approximately 0.1 m and hence q i is estimated as

(C.4)$$q_{{\rm i}} = \rho_{{\rm i}} L w = 3\times 10^{7}\, \rm{J}\, \rm{m}^{-2},\; $$

where L is the latent heat of melting.

Table 1. Constants used in calculation of temperature around basal crevasses

Thus we have the full thermal problem in an infinite strip with an initial line heat source representing refreezing in real basal crevasses.

The solutions are given by Carslaw and Jaeger (Reference Carslaw and Jaeger1959) as

(C.5)$$\Delta T = {q_{\rm i}\over \pi\rho_{\rm i} C_{\rm p}\sqrt{\kappa \pi t}}\exp{\left(-{\left(x-x^{\prime}\right)^{2}\over 4\kappa t}\right)} \sum_{{n = 1}}^{ + \infty}{1\over n}\left(1-\cos\left({{n\pi Z_{{\rm C}}\over H}}\right)\right)\sin\left({{n\pi z\over H}}\right)\exp\left({{-\kappa n^{2} \pi^{2} t\over H^{2}}}\right).$$

C.2 Thermal structure of a series of equally spaced basal crevasses

A stable sticky patch can cause a local stress variation and potentially generate a series of basal crevasses. The thermal effect of a single basal crevasse, which is shown above, can be advected downstream and superposed with the effect of other basal crevasses, resulting in temperature variation on a longer timescale and larger area. Mathematically, based on ΔT (x, z, t) that we have obtained above, the net effect of a series of basal crevasses can be obtained by a simple linear superposition of the ΔT of each crevasse.

Assuming these crevasses are generated at a fixed spacing w s = W = 100 m in an 100 m-thick ice sheet, starting from t = 0, we fix the coordinate system on the sticky patch and move the ice sheet at v s = 100 m × y−1. Thus the crevasses are also generated at a fixed time interval Δt = W/v s = 1 y. The net temperature perturbation ΔT net is

(C.6)$$\Delta T_{{\rm net}} = \sum_{n = 0}^{\infty} \Delta T\left(x-x_{n}^{\prime},\; z,\; t-t_{n}^{\prime}\right),\; $$

where n is the number of the crack, $x_{n}^{\prime } = n w_{s}$ is the position of the crack n, $t_{n}^{\prime } = -n\Delta t$ is the time when the crack n formed and refroze.

References

Andrews, LC and 7 others (2014) Direct observations of evolving subglacial drainage beneath the Greenland Ice Sheet. Nature 514(7520), 80–83. doi: 10.1038/nature13796CrossRefGoogle ScholarPubMed
Banks-Sills, L and Sherman, D (1986) Comparison of methods for calculating stress intensity factors with quarter-point elements. International Journal of Fracture 32(2), 127–140. doi: 10.1007/BF00019788CrossRefGoogle Scholar
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.003CrossRefGoogle Scholar
Berg, B and Bassis, J (2022) Crevasse advection increases glacier calving. Journal of Glaciology 1–10. doi: 10.1017/jog.2022.10CrossRefGoogle Scholar
Brædstrup, C, Egholm, D, Ugelvig, S and Pedersen, V (2016) Basal shear stress under alpine glaciers: insights from experiments using the iSOSIA and Elmer/Ice models. Earth Surface Dynamics 4(1), 159–174. doi: 10.5194/esurf-4-159-2016CrossRefGoogle Scholar
Broek, D (1982) Elementary Engineering Fracture Mechanics, 3rd ed. Netherlands: Martinus Nijhoff Publishers.Google Scholar
Carslaw, HS and Jaeger, JC (1959) Conduction of Heat in Solids. UK: Clarendon P Oxford.Google Scholar
Cathles, LM (2015) Viscosity of the Earth's mantle. Princeton: Princeton University Press. doi: 10.1515/9781400867684.Google Scholar
Chan, S, Tuba, I and Wilson, W (1970) On the finite element method in linear fracture mechanics. Engineering Fracture Mechanics 2(1), 1–17. doi: 10.1016/0013-7944(70)90026-3CrossRefGoogle Scholar
Crouch, SL and Starfield, A (1983) Boundary Element Methods in Solid Mechanics. London: Goerge Allen and Unwin Publishers. doi: 10.1115/1.3167129).CrossRefGoogle Scholar
Dahm, T (2000) Numerical simulations of the propagation path and the arrest of fluid-filled fractures in the Earth. Geophysical Journal International 141(3), 623–638. doi: 10.1046/j.1365-246x.2000.00102.xCrossRefGoogle Scholar
Daniels, KA, Bastow, I, Keir, D, Sparks, R and Menand, T (2014) Thermal models of Dyke intrusion during development of continent–ocean transition. Earth and Planetary Science Letters 385, 145–153. doi: 10.1016/j.epsl.2013.09.018CrossRefGoogle Scholar
Davis, T (2017) A new open source boundary element code and its application to geological deformation: exploring stress concentrations around voids and the effects of corrugation orientation on fault slip in 3D (Master's thesis). Aberdeen University.Google Scholar
Engelhardt, H, Humphrey, N, Kamb, B and Fahnestock, M (1990) Physical conditions at the base of a fast moving Antarctic ice stream. Science 248(4951), 57–59. doi: 10.1126/science.248.4951.57CrossRefGoogle ScholarPubMed
Geuzaine, C and Remacle, JF (2009) Gmsh: a 3-D finite element mesh generator with built-in pre-and post-processing facilities. International Journal for Numerical Methods in Engineering 79(11), 1309–1331. doi: 10.1002/nme.2579CrossRefGoogle Scholar
Glen, JW (1955) The creep of polycrystalline ice. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 228(1175), 519–538. doi: 10.1098/rspa.1955.0066Google Scholar
Gow, A and 6 others (1997) Physical and structural properties of the Greenland Ice Sheet Project 2 ice core: a review. Journal of Geophysical Research: Oceans 102(C12), 26559–26575. doi: 10.1029/97JC00165CrossRefGoogle Scholar
Guinea, GV, Planas, J and Elices, M (2000) KI evaluation by the displacement extrapolation technique. Engineering Fracture Mechanics 66(3), 243–255. doi: 10.1016/S0013-7944(00)00016-3CrossRefGoogle Scholar
Harper, JT, Bradford, JH, Humphrey, NF and Meierbachtol, TW (2010) Vertical extension of the subglacial drainage system into basal crevasses. Nature 467(7315), 579–582. doi: 10.1038/nature09398CrossRefGoogle ScholarPubMed
Harper, JT, Humphrey, NF, Pfeffer, WT, Fudge, T and O'Neel, S (2005) Evolution of subglacial water pressure along a glacier's length. Annals of Glaciology 40, 31–36. doi: 10.3189/172756405781813573CrossRefGoogle Scholar
Helmstetter, A, Moreau, L, Nicolas, B, Comon, P and Gay, M (2015) Intermediate-depth icequakes and harmonic tremor in an Alpine glacier (Glacier d'Argentière, France): evidence for hydraulic fracturing? Journal of Geophysical Research: Earth Surface 120(3), 402–416. doi: 10.1002/2014JF003289CrossRefGoogle Scholar
Hewitt, I (2013) Seasonal changes in ice sheet motion due to melt water lubrication. Earth and Planetary Science Letters 371–372, 16–25. doi: 10.1016/j.epsl.2013.04.022CrossRefGoogle Scholar
Hudson, TS and 5 others (2020) Breaking the ice: identifying hydraulically forced crevassing. Geophysical Research Letters 47(21), e2020GL090597. doi: 10.1029/2020GL090597CrossRefGoogle Scholar
Hulbe, CL, LeDOUX, C and Cruikshank, K (2010) Propagation of long fractures in the Ronne Ice Shelf, Antarctica, investigated using a numerical model of fracture propagation. Journal of Glaciology 56(197), 459–472. doi: 10.3189/002214310792447743CrossRefGoogle Scholar
Jarvis, GT and Clarke, GK (1974) Thermal effects of crevassing on Steele glacier, Yukon Territory, Canada. Journal of Glaciology 13(68), 243–254. doi: 10.3189/S0022143000023054CrossRefGoogle Scholar
Jimenez, S and Duddu, R (2018) On the evaluation of the stress intensity factor in calving models using linear elastic fracture mechanics. Journal of Glaciology 64(247), 759–770. doi: 10.1017/jog.2018.64CrossRefGoogle 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-2014CrossRefGoogle Scholar
Lai, CY and 7 others (2020) Vulnerability of Antarctica's ice shelves to meltwater-driven fracture. Nature 584(7822), 574–578. doi: 10.1038/s41586-020-2627-8CrossRefGoogle ScholarPubMed
Langtangen, HP and Logg, A (2017) Solving PDEs in Python: The FEniCS Tutorial I. Berlin: Springer International Publishing. doi: 10.1007/978-3-319-52462-7.Google Scholar
Law, R and 10 others (2021) Thermodynamics of a fast-moving Greenlandic outlet glacier revealed by fiber-optic distributed temperature sensing. Science Advances 7(20), eabe7136. doi: 10.1126/sciadv.abe7136CrossRefGoogle ScholarPubMed
Lipovsky, BP and 6 others (2019) Glacier sliding, seismicity and sediment entrainment. Annals of Glaciology 60(79), 182–192. doi: 10.1017/aog.2019.24CrossRefGoogle Scholar
Lipovsky, BP (2020) Ice shelf rift propagation: stability, three-dimensional effects, and the role of marginal weakening. The Cryosphere 14(5), 1673–1683. doi: 10.5194/tc-14-1673-2020CrossRefGoogle Scholar
Logg, A, Mardal, KA and Wells, G (2012) Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Vol. 84. Berlin, Heidelberg: Springer Science & Business Media.CrossRefGoogle Scholar
Logg, A and Wells, GN (2010) DOLFIN: automated finite element computing. ACM Transactions on Mathematical Software (TOMS) 37(2), 1–28. doi: 10.1145/1731022.1731030CrossRefGoogle Scholar
Luckman, A and 5 others (2012) Basal crevasses in Larsen C Ice Shelf and implications for their global abundance. The Cryosphere 6(1), 113–123. doi: 10.5194/tc-6-113-2012CrossRefGoogle Scholar
McCarthy, C, Savage, H and Nettles, M (2017) Temperature dependence of ice-on-rock friction at realistic glacier conditions. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 375(2086), 20150348. doi: 10.1098/rsta.2015.0348CrossRefGoogle ScholarPubMed
McDowell, IE, Humphrey, NF, Harper, JT and Meierbachtol, TW (2021) The cooling signature of basal crevasses in a hard-bedded region of the Greenland Ice Sheet. The Cryosphere 15(2), 897–907. doi: 10.5194/tc-15-897-2021CrossRefGoogle Scholar
Rist, M and 5 others (1996) Experimental fracture and mechanical properties of Antarctic ice: preliminary results. Annals of Glaciology 23, 284–292. doi: 10.3189/S0260305500013550CrossRefGoogle Scholar
Schulson, E, Lim, P and Lee, R (1984) A brittle to ductile transition in ice under tension. Philosophical Magazine A 49(3), 353–363. doi: 10.1080/01418618408233279CrossRefGoogle Scholar
Sergienko, OV and Hindmarsh, RC (2013) Regular patterns in frictional resistance of ice-stream beds seen by surface data inversion. Science 342(6162), 1086–1089. doi: 10.1126/science.1243903CrossRefGoogle ScholarPubMed
Sergienko, OV, MacAyeal, DR and Bindschadler, RA (2009) Stick–slip behavior of ice streams: modeling investigations. Annals of Glaciology 50(52), 87–94. doi: 10.3189/172756409789624274CrossRefGoogle Scholar
Sinha, NK (1978) Short-term rheology of polycrystalline ice. Journal of Glaciology 21(85), 457–474. doi: 10.3189/S002214300003361XCrossRefGoogle Scholar
Smith, R (1976) The application of fracture mechanics to the problem of crevasse penetration. Journal of Glaciology 17(76), 223–228. doi: 10.3189/S0022143000013563CrossRefGoogle Scholar
Stokes, CR, Clark, CD, Lian, OB and Tulaczyk, S (2007) Ice stream sticky spots: a review of their identification and influence beneath contemporary and palaeo-ice streams. Earth-Science Reviews 81(3-4), 217–249. doi: 10.1016/j.earscirev.2007.01.002CrossRefGoogle Scholar
Tada, H, Paris, P and Irwin, G (2000) The Analysis of Cracks Handbook, Vol. 2. New York: ASME Press.Google Scholar
Van der Veen, C (1998) Fracture mechanics approach to penetration of bottom crevasses on glaciers. Cold Regions Science and Technology 27(3), 213–223. doi: 10.1016/S0165-232X(98)00006-8CrossRefGoogle Scholar
Van der Veen, CJ (2013) Fundamentals of Glacier Dynamics. Boca Raton, FL: CRC Press, Taylor and Francis Group.CrossRefGoogle Scholar
Walter, F, Dalban Canassy, P, Husen, S and Clinton, JF (2013) Deep icequakes: what happens at the base of alpine glaciers? Journal of Geophysical Research: Earth Surface 118(3), 1720–1728. doi: 10.1002/jgrf.20124CrossRefGoogle Scholar
Wearing, MG and Kingslake, J (2019) Holocene formation of Henry Ice Rise, West Antarctica, inferred from ice-penetrating radar. Journal of Geophysical Research: Earth Surface 124(8), 2224–2240. doi: 10.1029/2018JF004988CrossRefGoogle Scholar
Wiens, DA, Anandakrishnan, S, Winberry, JP and King, MA (2008) Simultaneous teleseismic and geodetic observations of the stick–slip motion of an Antarctic ice stream. Nature 453(7196), 770–774. doi: 10.1038/nature06990CrossRefGoogle ScholarPubMed
Wolovick, MJ, Creyts, TT, Buck, WR and Bell, RE (2014) Traveling slippery patches produce thickness-scale folds in ice sheets. Geophysical Research Letters 41(24), 8895–8901. doi: 10.1002/2014GL062248CrossRefGoogle Scholar
Yu, H, Rignot, E, Morlighem, M and Seroussi, H (2017) Iceberg calving of Thwaites Glacier, West Antarctica: full-Stokes modeling combined with linear elastic fracture mechanics. The Cryosphere 11(3), 1283–1296. doi: 10.5194/tc-11-1283-2017CrossRefGoogle Scholar
Zoet, L and 6 others (2020) Application of constitutive friction laws to glacier seismicity. Geophysical Research Letters 47(21), e2020GL088964. doi: 10.1029/2020GL088964Google Scholar
Figure 0

Fig. 1. Schematic diagram of the computational domain: a 2L × H strip extracted from the ice sheet on a slope with angle α. (a) The first sub-problem where ice is sliding at a constant rate under a uniform basal drag τ0. (b) The second sub-problem where there are a stick patch and a hydrofracture downstream. The sticky patch is represented by an interval [ − W, W] on the x-axis with excess basal stress Δτ. A water-filled crack is initiated on the downstream end of the sticky patch. The crack can be a vertical line with length ZC (as illustrated in the figure), or a curve, depending on the criterion used to determine its propagation.

Figure 1

Fig. 2. Perturbation stress in a cracked ice sheet with W′ = 1, Δτ′ = 0.3, and $Z_{\rm C}^{\prime } = 0.3$. Panels show (a) $T_{xx}^{\prime }$, (b) $T_{zz}^{\prime }$ and (c) $T_{xz}^{\prime }$. Tension is defined as having a positive sign. In the zoom-in box a white rectangle is added to highlight the crack, which is actually much narrower, making it difficult to see otherwise.

Figure 2

Fig. 3. Normalised stress intensity factor KI/KI,C versus dimensionless crack length with W′ = 1. Two cases with different Δτ′ are considered: Δτ′ = 0.2 (red) and Δτ′ = 0.3 (blue). Vertical axis represents the dimensionless crack length. Horizontal axis represents the dimensionless SIFs normalised by the fracture toughness, $K_{{I}, {\rm C}}^{\prime } = K_{{I}, {\rm C}}/( \rho _{\rm i} g H^{3/2})$, where H is chosen to be 100 m. The black dash-dotted line represents the fracture toughness.

Figure 3

Fig. 4. Maximum crack length $Z_{{\rm C, max}}^{\prime }$ as a function of W′ and Δτ′. Note that in the criterion, the non-dimensional fracture toughness depends on the ice thickness H. Here we set $H = 100\ \rm {m}$. For larger thicknesses, $K_{{I}, {\rm C}}^{\prime }$ gets smaller, leading to longer cracks. Two cases with different f are considered: (a) f = 0.7, (b) f = 0.9. When the crack length $Z_{{\rm C}}^{\prime } = 1$, the crack dissects the full ice-sheet thickness.

Figure 4

Fig. 5. Fracture paths calculated under different excess basal shear stress Δτ′ and flotation fraction f, with three values of W′: (a) W′ = 0.1; (b) W′ = 1.0; (c) W′ = 10.0. An ice thickness H = 103 m is used to scale the fracture toughness. The background vector field indicates the local direction of the maximum compressive stress, with the vector length scaled by the deviatoric stress $\left \vert {T_{1}^{\prime }-T_{2}^{\prime }}\right \vert$, where $T_{1}^{\prime }$ and $T_{2}^{\prime }$ are the local principal stresses calculated from T′. In panels (a) and (b), the fracture paths (colours) are calculated by BEM, with the principal stress trajectory represented by the dashed curve. Colours represent four different combinations of Δτ′ and f as in the legend. In panel (c), since the principal stress trajectory is nearly a vertical line, we no longer conduct the BEM simulation and just assume purely vertical fracture paths. Thus, the maximum crack length is determined as we have done in Figure 3. Colours represent three different combinations of Δτ′ and f as in the legend.

Figure 5

Fig. 6. Perturbation of a single basal crevasse on background temperature field. t = 0 is the time of crevasse opening and refreezing. (a) t = 5 years. (b) t = 10 years. (c) t = 15 years. (d) t = 20 years.

Figure 6

Fig. 7. Thermal effect of a series of basal crevasses that is 2–3 km downstream from the sticky patch. After a long time of diffusion, the localised temperature field smooths out.

Figure 7

Fig. 8. Stress intensity factors KI and KII when Δτ′ = 0.3 and f = 0.7. The red lines show SIFs calculated by the weight function method (‘WF’ in the legend means ‘weight function’). The blue dots show SIFs calculated by DCM.

Figure 8

Fig. 9. Schematic of a series of basal crevasses produced on a sticky patch.

Figure 9

Table 1. Constants used in calculation of temperature around basal crevasses