Basal hydrofractures near sticky patches

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.

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, 2015; Lai and others, 2020; Lipovsky, 2020). In grounded ice, basal crevasses may have important influences on subglacial hydrology and glacier dynamics (Walter and others, 2013). Crevasses are also a potential explanation for seismic activity detected in glaciers and ice sheets (Hudson and others, 2020). Clustered, deep icequakes with non-double-couple sources have been detected in alpine glaciers (Walter and others, 2013; Helmstetter and others, 2015). 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 (2013) 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, 2018).

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, 2010). Such crevasses, filled by high-pressure water, serve as important components of the subglacial hydrology system. In West Antarctica, Wearing and Kingslake (2019) 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, 2013). 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, 2007). The existence of sticky patches has significant influence on the nonuniformity of ice-flow dynamics (Wolovick and others, 2014).

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 (1976) and Van der Veen (1998). At long timescales, glacial ice deforms as a viscous fluid, with a viscosity modelled by Glen's law (Glen, 1955), 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 (1998) 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 (2018) 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 (2018) recommended a LEFM weight function from Tada and others (2000) 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, 2010) and three dimensions (Lipovsky, 2020). 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 τ₀. (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 τ₀, representing background basal drag, is imposed on (− L, L). The uniform basal stress τ₀ 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, 2013).

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 τ₀ 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 τ₀ (Sergienko and Hindmarsh, 2013). 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, 1982). 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, 2010). 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 (2015). This approach was used by Lipovsky (2020) 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, 1998). 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, 1982),

(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, 1982). 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 m^1/2 (Rist and others, 1996), 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, 2010; Logg and others, 2012; Langtangen and Logg, 2017), with meshes generated by Gmsh (Geuzaine and Remacle, 2009). 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⁻³ 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, 1970; Banks-Sills and Sherman, 1986) with Richardson extrapolation (Guinea and others, 2000). The code is benchmarked by comparison with a weight function from Jimenez and Duddu (2018) (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, 1983), 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, 2000). The BEM code used here is developed by Davis (2017). 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, 1990; Harper and others, 2010; Andrews and others, 2014). 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, 1998). 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, 1998). Here we simply assume a pre-existing vertical line crack with K _I > K _I,C.

Fig. 3. Normalised stress intensity factor K_I/K_I,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, 1998). 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 = 10³ 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, 2014; Yu and others, 2017). 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 km² (Stokes and others, 2007), 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, 2000)

(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 (2013) showed that for a region of the Antarctic ice sheet with large ice thickness (H ≈ 10³ 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 (1998). 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, 2016). 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, 2007). An important friction phenomenon is the stick-slip motion of ice, detected in Whillans Ice Stream (WIS) in West Antarctica. Wiens and others (2008) investigated the WIS stick-slip motion and related it to a sticky patch on the bed. Furthermore, Sergienko and others (2009) 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, 2017; Lipovsky and others, 2019; Zoet and others, 2020), 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 (2020) 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, 2007), 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 (2021) 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, 2012) 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, 2014), has been modelled in several studies (Jarvis and Clarke, 1974; McDowell and others, 2021). 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 (2021) 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 (1959). 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, 2021). 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 (2010). 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, 2021).

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, 2015). The recent numerical model developed by Berg and Bassis (2022) 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, 2005; Hewitt, 2013), 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, 1984). Regions and layers of coarse-grained ice exist in ice sheets (Gow and others, 1997). 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 (1978) 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, 1998), 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, 1998; Krug and others, 2014; Jimenez and Duddu, 2018). 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, 2013; Brædstrup and others, 2016), 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.