Ice viscosity governs hydraulic fracture that causes rapid drainage of supraglacial lakes

Full-thickness crevasses can transport water from the glacier surface to the bedrock where high water pressures can open kilometre-long cracks along the basal interface, which can accelerate glacier flow. We present a first computational modelling study that describes time-dependent fracture propagation in an idealised glacier causing rapid supraglacial lake drainage. A novel two-scale numerical method is developed to capture the elastic and viscoelastic deformations of ice along with crevasse propagation. The fluid-conserving thermo–hydro–mechanical model incorporates turbulent fluid flow and accounts for melting and refreezing in fractures. Applying this model to observational data from a 2008 rapid-lake-drainage event indicates that viscous deformation exerts a much stronger control on hydrofracture propagation compared to thermal effects. This finding contradicts the conventional assumption that elastic deformation is adequate to describe fracture propagation in glaciers over short timescales (minutes to several hours) and instead demonstrates that viscous deformation must be considered to reproduce observations of lake drainage rates and local ice surface elevation changes. As supraglacial lakes continue expanding inland and as Greenland Ice Sheet temperatures become warmer than −8 °C, our results suggest rapid lake drainage events are likely to occur without refreezing, which has implications for the rate of sea level rise.

1 Introduction

Mass loss through melting and iceberg calving from the Greenland Ice Sheet (GrIS) will continue increasing throughout the century in response to atmospheric warming (Bamber et al., 2019; Bevis et al., 2019; Pattyn, 2018). Under the highest baseline emissions scenario, RCP8.5, mass loss from the GrIS could raise global sea level by up to 90 mm ± 50 mm by 2100 (Goelzer et al., 2020). While this mass loss is dominated by increased runoff derived from surficial melting (Hofer et al., 2020), the influence of meltwater increase on ice dynamics remains unresolved and unaccounted for in ice sheet models. To better parameterise the coupling between melting and ice dynamics in ice sheet models, it is essential to develop and utilise advanced physics-based models to better understand the processes at the smaller (glacier) scale.

When crevasses reach the ice sheet's base, they can transfer liquid water from the surface to the base of ice sheets, where it can influence ice dynamics by modulating pressures within the subglacial drainage system. Meltwater produced on the surface of glaciers and ice sheets often collects in depressions, forming supraglacial lakes that can drain rapidly when intersected by a crevasse (Das et al., 2008; Doyle et al., 2013; Chudley et al., 2019; Selmes et al., 2011; Smith et al., 2015; Christoffersen et al., 2018). When filled with water, crevasses can continue to propagate deeper in the ice sheet until reaching the bed (van der Veen, 2007; Weertman, 1973), thereby creating direct connections between the supraglacial and subglacial drainage systems. It has been observed that only a limited amount of water is required to create this connection (Krawczynski et al., 2009; Selmes et al., 2011). However, when these hydraulically driven fractures connect to supraglacial lakes, they have the ability to transport massive amounts of water to the ice sheet's base (Lai et al., 2021). These large volumes of water quickly overwhelm and pressurise the subglacial drainage system to initiate horizontal basal fracture propagation and uplift, changes in basal friction, and acceleration of the overlying ice (Das et al., 2008; Stevens et al., 2015; Doyle et al., 2013; Liang et al., 2012; Chudley et al., 2019; Shannon et al., 2013), which can extend tens of kilometres downglacier (Andrews et al., 2018; Hoffman et al., 2011; Mejía et al., 2021).

Crevasses are usually modelled as opening (mode I) fractures penetrating through the ice thickness under the action of tensile stress (van der Veen, 2007). Most existing models for estimating the depth of water-filled crevasses rely on the assumption that the flow of water associated with fracture propagation is such that hydrostatic conditions are valid and the thermal process is slow enough that melting and freezing can be neglected. Based on this assumption, analytical linear-elastic-fracture-mechanics-based (LEFM-based) models were used to estimate the depth of water-filled crevasses, namely, the Nye zero stress (Jezek, 1984; Benn et al., 2007; Nick et al., 2010), dislocation-based LEFM (Weertman, 1971), and stress-intensity-factor-based LEFM models (Smith, 1976; Van Der Veen, 1998). Recently, stress-intensity-factor-based LEFM models for estimating water-filled crevasse depth have been proposed using finite-element-based (Jimenez and Duddu, 2018) and boundary-element-based methods (Zarrinderakht et al., 2022). Alternatively, continuum damage mechanics and phase field models for water-filled crevasse propagation have been more recently developed (Duddu et al., 2020; Clayton et al., 2022). In the abovementioned studies, except for Zarrinderakht et al. (2022), the water column height within the crevasse rather than the volume of water was prescribed to drive crevasse propagation.

The mechanical response of glacier ice is usually described by the isotropic linear-elastic model over short timescales (microseconds to a few minutes) for simplicity, but over longer timescales (hours to months), ice response is better described by the nonlinear viscoelastic model (i.e. dependent on the strain rate, temperature, and time) known as Glen's law (Glen, 1955). A few studies in the literature incorporated viscous deformation of ice and investigated longer-term processes such as the slow propagation of crevasses over weeks to months (Poinar et al., 2017) or moulin formation post-hydrofracture (Andrews et al., 2022). Other studies described hydraulically driven (vertical) crevasse propagation that subsequently leads to horizontal fracture propagation at the ice–bedrock interface and ice sheet uplift, using approximate analytic solutions (Tsai and Rice, 2010, 2012) or using large-scale mass and momentum balances (Rice et al., 2015; Andrews et al., 2022). Within numerical simulations, reduced-order models have been used to include the interactions between basal water and basal friction (Pimentel and Flowers, 2011). Some studies utilised the finite-element method to accurately simulate the deformations of the ice sheet including creep (Hewitt et al., 2012; De Fleurian et al., 2014; Crawford et al., 2021). However, there exists no comprehensive numerical model capable of describing both the vertical and horizontal fracturing and uplifting due to the turbulent flow of pressurised meltwater within the fracture, along with melting and refreezing.

In this work, we present a comprehensive numerical model for hydraulic fracturing that takes into account elastic-viscoplastic deformation and thermal processes driven by turbulent water flow (see Methods for full details). We use this model to investigate how the choice of a viscoelastic or linear-elastic ice sheet rheology influences fracture propagation, thermal processes (refreezing and frictional heating within crevasses), and ice sheet uplift. We also examine the requirements for fluid supply from supraglacial lakes and the timescales involved during the fracturing process by comparing model outputs with observations from rapid lake drainage on the Greenland Ice Sheet. Ultimately, our results highlight the important role of creep deformation in facilitating rapid supraglacial lake drainage events, thereby supporting the incorporation of a viscoelastic ice rheology even on short (hours) timescales, whereas the effects of melting were less relevant at these timescales.

Figure 1 Schematic diagrams describing the physical phenomena associated with lake drainage driven by hydraulic fracture. (a) In the 3D real case scenario, the vertical crevasse below the lake propagates and eventually reaches the bedrock upon which fluid flow is sustained by horizontal fracture propagation and basal uplift. (b) Assuming uniformity in the out-of-plane direction, the problem can be idealised assuming 2D plane-strain conditions to reduce the computational burden.

2 Models and methods

To investigate the hydraulic fracturing process induced by supraglacial lakes, the 3D geometry from Fig. 1 is simplified by assuming crevasse opening is uniform over the full out-of-plane width (see Fig. 1). This allows the domain to be approximated as a 2D domain, shown in Fig. 2. Here, we note that this assumption requires the out-of-plane length of the crevasse to be large, such that basal uplift is produced by the unidirectional fluid flow away from the vertical crevasse. An alternative assumption to reduce the problem from three dimensions to two dimensions would be to consider the crevasse as a vertical cylindrical conduit to allow for radial subglacial water flow. Using an axisymmetric representation would be appropriate for crevasses with much shorter surface lengths (i.e. the out-of-plane direction from Fig. 1) compared to the length of the horizontal/radial basal crack, such that the vertical crevasse can be considered a conduit propagating downwards. In contrast, the plane-strain representation used here is suitable for when the surface crevasse spans longer distances, such that the crack is considered a plane propagating downwards. While our 2D model for horizontal basal crack propagation and basal uplift is valid for the axisymmetric assumption, it would not be appropriate to assume that the vertical crevasse is radially symmetric because it would no longer represent a planar fracture undergoing opening but rather the creation of a cylindrical conduit. This is an important distinction because of the associated processes governing fracture propagation. For a typical planar or vertical crevasse formed under plane strain, mechanical stresses drive crevasse propagation. In contrast, under axisymmetric conditions, the cylindrical conduit would propagate downwards due to fracture, melting, and erosion processes. While it is possible to model it as a cylindrical moulin evolution (Trunz et al., 2022), it is difficult to describe the fracture mechanics using the existing framework. Of course a 3D model able to capture both these phenomena would be ideal, but it would be computationally too expensive, so we utilise the 2D plane-strain approximation, focusing on the propagation of fractures driven by stresses and only considering the lateral melting of the fracture faces to ultimately align with the observed morphology of crevasses.

Figure 2 Schematic sketch showing the macro-scale domain description, boundary conditions, and the in-plane x–y coordinates assuming 2D plane-strain conditions. The ice and rock domain is denoted by Ω, and the fracture interface is denoted by Γ_d. The inset image shows the micro-scale domain and the corresponding ξ–η coordinates used for describing turbulent fluid flow and thermal conduction.

It is further assumed that crevasse length is much larger than its width, defined as the spacing between crevasse walls measured at the ice surface. This assumption allows the separation of two distinct spatial scales: the glacier scale on which deformations occur and the fracture scale that determines the fluid flow and thermal effects (melting or refreezing, frictional heating, and conduction). The 2D domain consists of a 6 km × 980 m portion of an ice sheet on a 200 m thick rock layer. The thickness of the rock layer is chosen to be large enough that the stresses around the basal crevasse are not spuriously altered by the zero-vertical-displacement boundary condition applied to the bottom surface of the rock layer. An initial 30 m deep crevasse connected to a supraglacial lake is assumed to be present at the ice surface. Given the density difference between the water and surrounding ice, a 30 m crevasse depth is sufficient to overcome the tensile strength of the ice and initialise the hydraulically driven crevasse propagation.

The domain geometry is a 2D idealisation of a lake drainage event recorded in Das et al. (2008), utilising the same ice thickness and ice material properties (see Table 1). To enforce a no-slip basal boundary condition during crevasse propagation, we assume the ice–bedrock interface is frozen. Although the ice–bedrock interface is not frozen in the North Lake area of Greenland modelled here, we implement the frozen-bed assumption to isolate hydraulically induced basal crevasse growth. A full discussion of this rationale is provided in Sect. 2.1.2. A further simplification made is that no defects exist within the ice, and the rock layer is perfectly impermeable, preventing any fluid leak-off from the base of the ice sheet. Imposing no leak-off from the water-filled crevasse deviates from realistic conditions, where water may slowly seep into the bedrock layer at the base. This leak-off will (slightly) reduce the water pressure within the crevasse while the surrounding bedrock is saturated, slightly slowing down the propagation of the horizontal crack. While we acknowledge that these effects might play a significant role under certain circumstances, accurately capturing the full glacial drainage system (e.g. subglacial drainage channels) is beyond the scope of this study. In the discussion section, we comment on the potential implications of this simplification.

A finite-element formulation is used to accurately capture the fracturing and subsequent uplift process. This formulation includes the linear-elastic and nonlinear viscoelastic deformations of the ice.¹ This allows for the inclusion of both shorter-term (seconds to minutes) elastic deformations and longer-term (hours to months) viscous relaxation. As the nonlinear viscous response is dependent on the magnitude of the deviatoric stresses, its effect becomes relevant in areas with rapidly changing stress states (i.e. crack tips) even on shorter timescales. We explicitly include the conservation of mass within the crack to track the fluid flow as it moves through the crevasse. In addition, thermal processes (melting, frictional heating, and thermal diffusion) are included by means of a new small-scale formulation. In the remainder of this section, these different components will be discussed in more detail.

2.1 Momentum balance and constitutive models

The ice sheet through which hydrofracture occurs is modelled as the domain Ω in Fig. 2. This domain uses the plane-strain assumption to simplify the 3D geometry from Fig. 1 to the 2D geometry in Fig. 2. The domain is composed of a layer of ice resting on a layer of deformable rock, which are described through their displacements u and the history-dependent viscous strains ε_v. Within this domain a discontinuity is present, Γ_d, which represents both the vertical crevasse and the horizontal basal cracks (in both directions away from the vertical crevasse). The current condition of this fluid-filled fracture is described through its pressure p, and its history is included through the time-since-fracture t₀ and the thickness of wall melting/refreezing h_melt (positive for melting, negative for freezing).

Hydraulic fracture occurs rapidly over a short timescale. As such, the ice will exhibit both elastic deformations due to sudden fracture propagation, and some viscous relaxation with time. To describe these two mechanisms, ice is considered to be a deformable solid material (Duddu and Waisman, 2012), in contrast to the commonly used non-Newtonian viscous liquid for long-term ice flow simulations (Larour et al., 2012; Lipscomb et al., 2019). Within this description, the momentum balance is solved to obtain the displacements u throughout the domain:

$$\begin{array}{}\text{(1)}& {\mathit{\rho }}_{\mathit{\pi }}\ddot{\mathit{u}}-{\mathbf{L}}^{\text{T}}\mathit{\sigma }={\mathit{\rho }}_{\mathit{\pi }}\mathit{g},\end{array}$$

where the density ρ_π corresponds to either the bedrock density ρ_r or the ice density ρ_i; $\mathit{g}=[\mathrm{0}-\mathrm{9.81}\mathrm{m}{\mathrm{s}}^{-\mathrm{2}}{]}^{\mathrm{T}}$, the gravitational acceleration vector; and L is the matrix mapping the displacement to linearised strain vector $\mathit{\epsilon }=[{\mathit{\epsilon }}_{xx}{\mathit{\epsilon }}_{yy}{\mathit{\epsilon }}_{zz}{\mathit{\epsilon }}_{xy}{]}^{\mathrm{T}}=\mathbf{L}\mathit{u}$. The rock is assumed to be a linear-elastic solid, whereas ice is assumed to be a Norton–Hoff type elastic-viscoplastic material. Thus, the stress in ice is defined by the linear-elastic strain, which is obtained by offsetting the total strain with the viscous strain history. The constitutive law for ice and rock is defined by

$$\begin{array}{}\text{(2)}& \mathit{\sigma }={\mathbf{D}}_{\mathit{\pi }}{\mathit{\epsilon }}_{\text{e}}={\mathbf{D}}_{\mathit{\pi }}(\mathit{\epsilon }-{\mathit{\epsilon }}_{\mathrm{v}})={\mathbf{D}}_{\mathit{\pi }}\mathbf{L}\mathit{u}-{\mathbf{D}}_{\mathit{\pi }}{\mathit{\epsilon }}_{\mathrm{v}},\end{array}$$

where the isotropic linear-elastic stiffness matrix D_π is dependent on the material parameters of ice and rock. The viscous strains ε_v are obtained by integrating the viscous strain rate obtained from Glen's law (Glen, 1955) as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\dot{\mathit{\epsilon }}}_{\mathrm{v}}=& A{\left({\mathit{\sigma }}_{\text{dev}}^{\mathrm{T}}{\mathit{\sigma }}_{\text{dev}}\right)}^{\frac{n-\mathrm{1}}{\mathrm{2}}}{\mathit{\sigma }}_{\text{dev}}\\ =& A{\left(\left({\mathit{u}}^{\mathrm{T}}{\mathbf{L}}^{\mathrm{T}}-{\mathit{\epsilon }}_{\mathrm{v}}^{\mathrm{T}}\right){\mathbf{D}}^{\mathrm{T}}{\mathbf{P}}^{\mathrm{T}}\mathbf{PD}(\mathbf{L}\mathit{u}-{\mathit{\epsilon }}_{\mathrm{v}})\right)}^{\frac{n-\mathrm{1}}{\mathrm{2}}}\\ & \mathbf{PD}(\mathbf{L}\mathit{u}-{\mathit{\epsilon }}_{\mathrm{v}}),\end{array}\end{array}$$

where the projection matrix P is used to obtain the deviatoric part of the stress vector σ_dev=Pσ. The creep coefficient A is set to zero for the rock layer to solely allow for viscous deformation within the ice, and its temperature dependence is described using the Arrhenius law in the material properties section below. We assume temperature changes caused by the water-filled crevasse to be localised close to the crevasse, so the temperature field and the creep coefficient A are time-invariant within the bulk of the ice sheet.

Throughout this work, ice is considered a viscoelastic solid, which encompasses the special case of a linear-elastic solid. Within the viscoelastic model, Eqs. (2) and (3) are used to capture both the immediate elastic response and the slower viscous strains occurring over time. In contrast, the linear-elastic model neglects the viscous strains, setting ${\dot{\mathit{\epsilon }}}_{\mathrm{v}}$ to zero throughout the simulation.

2.1.1 Viscous and elastic timescales

The inclusion of viscous strains introduces a timescale over which the mechanical behaviour of ice changes from a compressible linear-elastic solid to incompressible viscoplastic solid. For materials described by a linear viscous model, this timescale is referred to as the relaxation time (also referred to as the Maxwell timescale), given by (Jellinek and Brill, 1956)

$$\begin{array}{}\text{(4)}& {\mathit{\tau }}_{\mathrm{m}}={\displaystyle \frac{\mathit{\eta }}{\frac{E}{\mathrm{2}(\mathrm{1}+\mathit{\nu })}}}={\displaystyle \frac{\mathrm{2}(\mathrm{1}+\mathit{\nu })}{EA{\mathit{\sigma }}_{\text{dev}}^{n-\mathrm{1}}}},\end{array}$$

where the effective viscosity for a Glen's law viscoelastic material is given by ${\mathit{\eta }}^{-\mathrm{1}}=A{\mathit{\sigma }}_{\text{dev}}^{n-\mathrm{1}}$ (with A being temperature-dependent via Eq. 21).

For the material properties used in this study (see Sect. 2.4), the range of this timescale is shown in Fig. 3. The immediate deformation response of ice is solely determined by the linear-elastic strain; however, over a time comparable to the Maxwell timescale, the viscous strains become dominant in driving ice deformation. The rate of increase in creep strain greatly varies depending on the local deviatoric stress. Away from the crack tip, the deviatoric stress is small (< 0.05 MPa), so the Maxwell timescale would be large (on the order of hours to days), implying that elastic strain dominates the deformation response. However, at the crack tip, the deviatoric stress is larger (≈ 0.21 MPa), so the Maxwell timescale would be small (on the order of minutes), implying that viscoelastic strains can influence the deformation response during hydraulic fracture. This is evident from Fig. 4, where the range of deviatoric stress values during a simulation is used to evaluate the range of Maxwell timescales existing within the glacier domain at two different ice temperatures.

Figure 3 Range of Maxwell timescales following from Eq. (4) and the material parameters from Table 1, assuming an effective deviatoric stress ranging from near zero to the tensile strength of the ice. The shaded area indicates the range of the timescale due to variations in temperature (and thus creep coefficient) with depth. A logarithmic scale is used for the vertical axis, ranging from 6 min to 100 h.

Figure 4 Range of Maxwell timescales following from Eq. (4) and the material parameters from Table 1, evaluated after 2 h of hydraulic fracture propagation. Deformations are magnified by 1000, and a logarithmic colour scale is used.

2.1.2 Cohesive fracture model

The force balance at the fracture interface is described by the traction $\mathit{\tau }=\mathit{\sigma }\cdot \mathit{n}$, which is decomposed into cohesive and pressure components as

$$\begin{array}{}\text{(5)}& \mathit{\tau }={\mathit{\tau }}_{\text{CZM}}+p\mathit{n},\end{array}$$

where p is the pressure of the water within the crevasse acting normal to the fracture faces. It is assumed that the pressurised fracture propagates solely in mode I (tension-driven fracture). When modelling sharp cracks within a linear-elastic material, the stresses at the crack tip become singular. Propagation criteria in linear-elastic fracture mechanics are then required to take this singularity into account, and capturing the stress intensity around the crack requires a very fine mesh. However, experimental observations in elasto-plastic materials reveal the existence of a finite fracture process zone or cohesive zone, wherein the stress is bounded by a characteristic material strength. To alleviate these physical and numerical issues with modelling sharp cracks within the finite-element method, a cohesive zone model (CZM) is used to regularise the stresses at the crack tip region. The CZM is a type of damage mechanics model, in which the traction response of the interface is degraded based on the crack separation or opening width due to damage evolution, and has been widely used to study fracture/delamination of composite materials (Ghosh et al., 2019) and hydraulic fracture in rocks (Hageman and de Borst, 2021; Hageman et al., 2019). Here we use the exponential traction–separation law to define the cohesive traction for mode I cracking as

$$\begin{array}{}\text{(6)}& {\mathit{\tau }}_{\text{CZM}}=-f_{\mathrm{t}}\mathit{n}\exp \left(-\textcolor{red}{}\mathit{u}\textcolor{red}{}\cdot \mathit{n}{\displaystyle \frac{f_{\mathrm{t}}}{G_{\mathrm{c}}}}\right).\end{array}$$

This cohesive traction depends on the fracture release energy G_c and the tensile strength of the ice f_t, unlike linear-elastic fracture mechanics that ignores tensile strength. Before cracking, the ice ahead of the crack tip is assumed to be fully intact and undamaged. The crack propagates once the stress within the ice, normal to the prescribed crack direction, exceeds the tensile strength: σ_yy>f_t for the horizontal crack and σ_xx>f_t for propagation of the vertical crevasse. As the weight of the ice induces compressive (negative) stresses within the ice sheet, the total change in stress needed to propagate the crack is therefore given by $(\mathit{\sigma }-{\mathit{\sigma }}_{\mathrm{0}})\cdot \mathit{n}>f_{\mathrm{t}}+{\mathit{\rho }}_{\mathrm{i}}g(H-y)$. For the horizontal crack, this implies that even through a frozen rock–ice interface is assumed with a tensile strength, the majority of the stresses that need to be overcome are a result of the weight of the ice and not the tensile strength. This frozen bed also imposes a no-slip boundary condition between the ice and the bed, requiring displacements in the ice to be matched by the basal rock it is contacting. Upon fracture, this constraint is released, and the ice and the rock can move independently in both normal direction (causing crack opening) and tangential direction (causing slip between the ice and rock layers). It should be noted that assuming the bed to be frozen has implications for the downward crevasse propagation and crevasse opening width after the vertical crevasse reaches the base, which is only driven by the elastic and viscous deformations. In contrast, were frictional sliding allowed at the glacier bed even before the onset of horizontal cracking and uplift, an additional opening width would be created due to the two sides of the ice sheet sliding apart. The effect of the glacier-sliding-induced crevasse widening would be significant if the basal friction were weak.

We use an extrinsic-type CZM implementation, in which the interface elements are inserted dynamically into the finite-element mesh as the crack propagates. As a result, the cohesive zone model is only applied post-cracking, while the not-yet-fractured interface behaves as an intact material. This is unlike other conventional intrinsic CZM implementations, wherein the interface elements are inserted a priori ahead of the crack tip along the potential crevasse path, which may result in additional and nonphysical displacements around the crack tip (Boone and Ingraffea, 1990). The exponential traction-separation law used here defines a length scale $\mathrm{\ell }\approx E{G}_{\mathrm{c}}/f_{\mathrm{t}}^{\mathrm{2}}$ (ℓ≈ 2.3 m for our case). This length scale gives an indication of the fracture process zone ahead of the traction-free crack, where the normal traction varies nonlinearly depending on the crack separation or opening width. Taking the length scale ℓ close to zero approximates brittle fracture but requires extremely small interface elements to accurately capture the traction both ahead of and behind the crack tip. In contrast, taking larger values of ℓ leads to a deviation from fully brittle fracture, instead emulating ductile effects near the crack tip, but allows larger elements to be used while still adequately capturing the stresses and propagation behaviour. Thus, through the choice of an appropriate length scale, our implementation allows for reasonably sized elements to be used, facilitating the simulation of cracks on the glacier scale in a thermodynamically consistent and computationally tractable manner.

2.2 Thermo–hydro–mechanical flow model within the fracture

Within the fracture, fluid flow is considered to be driven by gravity and a small amount of inlet pressure imposed at the crevasse mouth (i.e. the top of the domain). This pressure is required to allow the crevasse to open initially, after which the main driving force is the density difference between the glacier ice and water pressurising the fracture. The fluid flow induces heating of the fracture walls due to friction, while the surrounding ice leaches heat from the fluid flow. This behaviour is included through a two-scale scheme: resolving the fluid mass conservation within the fracture at the macro-scale (implemented through a standard finite-element scheme) while resolving the thermal processes as micro-scale effects (numerically resolved on a per-integration-point basis).

2.2.1 Pressure-driven flow model

The fluid (meltwater) flow is assumed to be solely contained within the fracture, and both the intact ice and rock are taken as impermeable and nonporous. Thus, we neglect the loss of fluid through fracture walls due to any flow through porous media. The fracture local coordinate frame $(\mathit{\xi },\mathit{\eta })=\mathbf{R}(x,y)$ is defined, using the rotation matrix R to align the coordinate frame with the fracture direction s and fracture normal n. A common assumption for flows within fractures is a laminar flow profile, resulting in a fluid flux given by the cubic law (Witherspoon et al., 1980). While this is accurate for most cases in hydraulic fracturing of rock where typical crack apertures are on the order of millimetres, crevasses have been observed to attain openings on the order of metres. To more accurately model the flow through such large crevasses, it is assumed that the combination of fracture aperture, wall roughness, and driving pressure are sufficiently large to cause the flow within the crack to exhibit a turbulent flow profile. The total fluid flux flowing through a fracture with opening h is therefore given through a Gauckler–Manning–Strickler-type flow law (Gauckler, 1867; Strickler, 1981), resulting in a volume flux q produced by turbulent flow as (Tsai and Rice, 2010, 2012)

$$\begin{array}{}\text{(7)}& -h\left({\displaystyle \frac{\partial p}{\partial \mathit{\xi }}}-{\mathit{\rho }}_{\mathrm{w}}\mathit{g}\cdot \mathit{s}\right)={\displaystyle \frac{f}{\mathrm{4}}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{h^{\mathrm{2}}}}\left|q\right|q,\end{array}$$

where ρ_w is the water density and g⋅s is the gravity component acting on the fluid. By comparing the fluid flux within the crevasse during the simulations, it has been verified that assuming a turbulent flow is a good approximation. For instance, the fluid flux near the inlet stabilises around 1000 ${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$ (see results section and Fig. 12), providing a Reynolds number for this flow as $\mathit{Re}={\mathit{\rho }}_{\mathrm{w}}\left|q\right|/{\mathit{\mu }}_{\mathrm{w}}\approx$ 3 × 10⁵, well above the typical transition point from laminar to turbulent flow (Re>𝒪(10³)). We note that while this fluid flux is indeed turbulent in the majority of the crevasse, it is unlikely to be close to the crack tips. For simplicity, we use this turbulent relation throughout the crevasse; however, slightly more realistic results could potentially be obtained by dynamically switching between laminar and turbulent flow models based on the current fluid flux. Although we do not present any results here for other flow models, simulations using a laminar flow model have also been performed (with the fluid flux scaling with $q=h^{\mathrm{3}}/\mathrm{12}\mathit{\mu }\partial p/\partial \mathit{\xi }$). In these simulations it was observed that the higher water flow rate resulted in faster crevasse propagation, indicating that the choice of flow model can directly impact the timescale over which the propagation process occurs. As altering the friction flow parameters can similarly alter the timescale imposed by the water flow, the presented simulation will use a single set of friction flow parameters taken from Tsai and Rice (2010). This ensures that there is no artificial fitting of our results to observations by tuning the flow parameters; instead our approach relies on including physically relevant processes to better capture reality. Even though the values used for these parameters are realistic, we do acknowledge that better fits of observed results could be attained by fine-tuning these friction parameters.

The total aperture of the fracture is obtained based on the displacement jump across the discontinuity and melting layer thickness as

$$\begin{array}{}\text{(8)}& h=h_{\text{melt}}+\mathit{n}\cdot \textcolor{red}{}\mathit{u}\textcolor{red}{}\end{array}$$

and the friction factor f related to the (constant) wall roughness k_wall and the reference friction factor f₀ as (Strickler, 1981)

$$\begin{array}{}\text{(9)}& f=f_{\mathrm{0}}{\left({\displaystyle \frac{k_{\text{wall}}}{h}}\right)}^{\frac{\mathrm{1}}{\mathrm{3}}}.\end{array}$$

Together, Eqs. (7) and (9) allow for the fluid transport within the crack to be described explicitly as

$$\begin{array}{}\text{(10)}& \begin{array}{rl}q=& -\mathrm{2}{\mathit{\rho }}_{\mathrm{w}}^{-\frac{\mathrm{1}}{\mathrm{2}}}{k}_{\text{wall}}^{-\frac{\mathrm{1}}{\mathrm{6}}}{f}_{\mathrm{0}}^{-\frac{\mathrm{1}}{\mathrm{2}}}{h}^{\frac{\mathrm{5}}{\mathrm{3}}}{\left|{\displaystyle \frac{\partial p}{\partial \mathit{\xi }}}-{\mathit{\rho }}_{\mathrm{w}}\mathit{g}\cdot \mathit{s}\right|}^{-\frac{\mathrm{1}}{\mathrm{2}}}\\ & \left({\displaystyle \frac{\partial p}{\partial \mathit{\xi }}}-{\mathit{\rho }}_{\mathrm{w}}\mathit{g}\cdot \mathit{s}\right).\end{array}\end{array}$$

To conserve the meltwater within the fracture, this flux must satisfy the mass balance equation (Réthoré et al., 2006; Carrier and Granet, 2012; de Borst, 2017; Hageman and de Borst, 2021):

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial q}{\partial \mathit{\xi }}}+\dot{h}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}}{\dot{h}}_{\text{melt}}+{\displaystyle \frac{h}{K_{\mathrm{w}}}}\dot{p}=\mathrm{0},\end{array}$$

where the four terms account for the changes in fluid flux, additional volume created through deformations and melting, fluid produced by the melting process, and the fluid compressibility, respectively. The compressibility term uses the bulk modulus of water K_w to allow it to be slightly compressible. While this term is near-negligible in practice (especially with the properties used in this paper), it provides a damping-like term within the numerical solution scheme and helps to stabilise the simulations. Through Eq. (11), the pressure within the crevasse is temporarily decreased when the crevasse opens, limiting the rate of crevasse opening and pressurisation based on the available fluid flow. This process introduces a timescale to the hydraulic fracturing process, which is often absent when a set meltwater height is used (Clayton et al., 2022; Sun et al., 2021), wherein the pressure is directly prescribed based on the local water height relative to the crack tip.

To complete the model, a free inflow condition is imposed at the inlet of the fracture (i.e. at the crevasse mouth) through a penalty-like approach:

$$\begin{array}{}\text{(12)}& q_{\text{inlet}}=k_{\mathrm{p}}(p_{\text{ext}}-p),\end{array}$$

where p_ext is the constant pressure at the inlet, taken as equal to the hydrostatic pressure of a 10 m deep lake, and k_p is a penalty parameter chosen large enough to accurately enforce this inflow constraint. We elect to enforce this pressure boundary condition in a weak sense through a penalty approach such that the fluid influx at the inlet can be easily recorded. It should be noted, however, that using Eq. (12) with a high value for k_p is equivalent to directly substituting the pressure at the inlet for the boundary condition, and it has been verified that the pressure at the inlet is approximately equal to the applied p_ext. For post-processing purposes, this flux is additionally integrated over time to obtain the total volume of water flowing into the crevasse as

$$\begin{array}{}\text{(13)}& Q_{\text{total}}=\underset{t}{\int }{q}_{\text{inlet}}\mathrm{d}t.\end{array}$$

As a 2D domain is considered, this produces the total volume of fluid that has entered the crevasse per unit of out-of-plane width (given in m³ m⁻¹). For the section where results are compared to observations from the literature, this inflow per unit width is converted to total inflow by assuming a prescribed out-of-plane width of the crevasse.

Equation (11) is solved to capture the macro-scale behaviour of the fluid contained within the fracture, whereas Eq. (10), defining the fluid flux, together with the thermal effects described in the next section, are solved to capture the micro-scale behaviour on a point-by-point basis. This novel and efficient two-scale approach is needed to account for the dynamic feedback between fluid flux and wall melting, which precludes a closed-form expression for the fluid flux.

2.2.2 Thermal model

The thermal processes within the fracture are described by assuming the heat transfer through advection within the crack to be negligible due to the near-freezing temperature of the lake water and limited opening width. As a result, the water within the fracture is required to be in thermal equilibrium with the surrounding ice, with any imbalance resulting in either additional ice melting or water freezing within the crack. This melting rate is described by

$$\begin{array}{}\text{(14)}& {\mathit{\rho }}_{\mathrm{i}}\mathcal{L}{\dot{h}}_{\text{melt}}-j_{\text{ice}}-j_{\text{flow}}=\mathrm{0},\end{array}$$

using the latent heat of fusion ℒ. Two contributions to the heat flux are considered in the above equation: (1) j_ice represents the heat being released by the ice into the water, which is always negative due to the surrounding ice by definition being below the freezing point, and (2) j_flow is the heat flux produced through turbulent flow, which is always positive. If the turbulent heat flux is larger than the ice absorption, $j_{\text{flow}}>-j_{\text{ice}}$, the ice will melt, as will be the case for a warm ice sheet that is at a temperature close to the ice melting point with a large fluid flux flowing through the crevasse. The opposite happens for a relatively cold ice sheet with little fluid flow, $j_{\text{flow}}<-j_{\text{ice}}$, where the fracture walls will start freezing. We emphasise that the thermal process is self-reinforcing: once the walls start melting a larger fluid flux is enabled, resulting in a larger amount of heat production, which leads to further melting of the walls. Conversely, when the walls start to freeze, the fracture will be more restrictive to fluid flow, causing the heat production to be more limited.

The flow-produced heat flux is given by (Andrews et al., 2022)

$$\begin{array}{}\text{(15)}& j_{\text{flow}}=-q\left({\displaystyle \frac{\partial p}{\partial \mathit{\xi }}}-{\mathit{\rho }}_{\mathrm{w}}\mathit{g}\cdot \mathit{s}\right),\end{array}$$

which scales with $(\frac{\partial p}{\partial \mathit{\xi }}-{\mathit{\rho }}_{\mathrm{w}}\mathit{g}\cdot \mathit{s}{)}^{\mathrm{3}/\mathrm{2}}$ and $h^{\mathrm{5}/\mathrm{3}}$, considering the definition of fluid flux q in Eq. (10). As such, it is expected that this heat flux becomes significant for large opening heights, which are more easily achieved for thicker ice sheets due to the increased overpressure sustained by the lake water. Furthermore, the gravity contribution to the driving force will cause increased melting for vertical cracks. Horizontal parts of the crack, in contrast, will exhibit reduced frictional heating, especially near the crack tips where the opening height and fluid flux are limited.

For the heat absorbed by the surrounding ice, it is assumed that once the fracture surfaces are created, heat conducts away normal to the fracture, described through the 1D heat conduction equation:

$$\begin{array}{}\text{(16)}& {\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{p}}\dot{T}-k{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial {\mathit{\eta }}^{\mathrm{2}}}}=\mathrm{0},\end{array}$$

using the heat capacity c_p and the thermal conductivity k, and using the surface-normal coordinate η. Because the ice–water boundary is consistently at T_w = 0 °C, an analytic solution for this equation is given by (White, 2006)

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{T(\mathit{\eta },t)-T_{\mathrm{\infty }}}{T_{\mathrm{w}}-T_{\mathrm{\infty }}}}=\text{erfc}\left({\displaystyle \frac{\mathit{\eta }}{\mathrm{2}\sqrt{\frac{k}{{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{p}}}(t-t_{\mathrm{0}})}}}\right),\end{array}$$

using erfc to indicate the complementary error function, T_∞ the temperature of the ice and t₀ the time from which heat conduction occurs, equal to the time of fracture. By taking the derivative with respect to η, and the temperature in degrees Celsius (such that T_w = 0 °C), the thermal flux at any point with distance η to the wall is given as

$$\begin{array}{}\text{(18)}& \begin{array}{rl}j\left(\mathit{\eta }\right)=& -k{\displaystyle \frac{\partial T}{\partial \mathit{\eta }}}=k\cdot {\displaystyle \frac{\mathrm{2}{T}_{\mathrm{\infty }}}{\mathrm{2}\sqrt{k/{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{p}}}\sqrt{t-t_{\mathrm{\infty }}}}}\\ & \cdot \text{exp}{\left(-{\displaystyle \frac{\mathit{\eta }}{\mathrm{2}\sqrt{k/{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{p}}}\sqrt{t-t_{\mathrm{\infty }}}}}\right)}^{\mathrm{2}},\end{array}\end{array}$$

where the factor π within this equation is a result of the derivative of complementary error function erfc from Eq. (17) and is not related to any axisymmetry assumptions: $\partial \text{erfc}\left(x\right)/\partial x=-\mathrm{2}\text{exp}\left(x^{\mathrm{2}}\right)/\sqrt{\mathit{\pi }}$. Taking the gradient at the wall, η=0, and accounting for the fact that each fracture contains two surfaces from which heat is absorbed, then results in the heat being absorbed into the ice (White, 2006):

$$\begin{array}{}\text{(19)}& j_{\text{ice}}=-\mathrm{2}k{\displaystyle \frac{\partial T}{\partial \mathit{\eta }}}={\displaystyle \frac{\mathrm{2}{k}^{\frac{\mathrm{1}}{\mathrm{2}}}{T}_{\mathrm{\infty }}{\mathit{\rho }}_{\mathrm{i}}^{\frac{\mathrm{1}}{\mathrm{2}}}{c}_{\mathrm{p}}^{\frac{\mathrm{1}}{\mathrm{2}}}}{{\mathit{\pi }}^{\frac{\mathrm{1}}{\mathrm{2}}}(t-t_{\mathrm{0}}{)}^{\frac{\mathrm{1}}{\mathrm{2}}}}}.\end{array}$$

As a result, the temperature changes throughout the ice sheet do not need to be resolved, and it is instead sufficient to keep track of the time since the local fracture was created. This time-since-fracture is then used to determine the heat flux at the wall during the simulation, greatly reducing the computational cost, and can be used to estimate the temperature close to the fracture. It should be noted, however, that by assuming the heat conduction to be localised near and normal to the crack, the fracture always needs to propagate through ice at its initial temperature. As a result, once the fracture propagates, the new crack surfaces start to release heat into the surrounding ice, and the surfaces start to freeze instantly. This is in contrast to assuming heat to be conducted in all directions away from the crack, in which case the ice ahead of the crack tip will be partially heated, depending on the crack propagation speed. As such, the model described here is solely accurate if crack propagation exceeds heat conduction, which, considering the low thermal conductivity of ice and small temperature differences, is a fairly reasonable criterion.

As a result of using analytic expressions for 1D heat conduction into the ice instead of separately simulating the temperature throughout the ice sheet, we explicitly separate the thermal problem on a micro-scale from the mechanical problem, which is considered the macro-scale. This assumes the thermal effects are localised near the crevasse, such that the overall mechanical material parameters (e.g. tensile strength and creep coefficient) are not altered by local changes in temperature and that the overall geometry is not changed by the melting process in the fractures. To establish the validity of this separation of length scales between heat conduction and mechanical fracture propagation, the length scale influenced by the heating due to the crevasse is estimated as ${\mathrm{\ell }}_{\text{thermal}}=\sqrt{tk/{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{p}}}$, which for the 2 h duration of our simulations is ℓ_thermal≈ 8 cm. As this length scale is orders of magnitude smaller compared to the crevasse length (and well below the element length used), using this separation of scales is a good approximation. This furthermore shows the strength of the described two-scale model: if we were to explicitly simulate the thermal processes throughout the complete ice sheet and wanted to include the crevasse, our spatial discretisation would need to be fine enough to capture this thermal length scale and would thus require centimetre-sized elements around the crevasse. Instead, by formulating the conduction through the analytic expression from Eq. (19), the thermal energy entering the ice is captured by the subgrid-scale formulation, and no separate discretisation is needed to accurately capture it.

Comparing the two heat fluxes, Eqs. (15) and (19), provides an indication whether the thermal processes are dominated by freezing or melting:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{\displaystyle \frac{j_{\text{flow}}}{-j_{\text{ice}}}}=& \sqrt{{\displaystyle \frac{\mathit{\pi }}{k{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{p}}{\mathit{\rho }}_{\mathrm{w}}{f}_{\mathrm{0}}{k}_{\text{wall}}^{\mathrm{1}/\mathrm{3}}}}}\\ & \cdot {\displaystyle \frac{(t-t_{\mathrm{0}}{)}^{\mathrm{1}/\mathrm{2}}{h}^{\mathrm{5}/\mathrm{3}}{\left|\frac{\partial p}{\partial \mathit{\xi }}-{\mathit{\rho }}_{\mathrm{w}}\mathit{g}\cdot s\right|}^{\mathrm{3}/\mathrm{2}}}{T_{\mathrm{\infty }}}},\end{array}\end{array}$$

with the first term composed of physical constants and the second term of case-dependent variables. The melting process dominates when $j_{\text{flow}}/-j_{\text{ice}}>\mathrm{1}$, coinciding with the case of large opening heights or warm ice sheets. The relevance of the melting process will also increase over time as the surrounding ice warms up due to the presence of the water-filled crevasse. In contrast, for short timescales, the freezing process will be dominant for almost all glacial temperatures, which will impose a limit on the ability of the hydraulic fracture to develop.

2.3 Implementation

The ice sheet fracture problem is described by the momentum balance in the domain Ω, Eq. (1), and the mass balance within the fracture Γ_d, Eq. (11). These two equations are discretised using the finite-element method, using quadratic quadrilateral elements for the displacements u and using quadratic interface elements for the fluid pressure p. Near the expected fracture path, these elements have a characteristic size of 2.5 m, whereas away from the interface elements up to 20 m are used. The temporal discretisation is performed using a Newmark scheme for the solid acceleration and an implicit backward Euler scheme for the pressure, wall melting, and integration of quantities for post-processing (total water inflow and thermal fluxes). The only exception quantity not using an implicit time discretisation scheme is the viscous strains, which change slowly compared to all other variables within the system. While it is possible to include the strain increments in a time-implicit manner, we elect here to evaluate this term using an explicit Euler scheme and thus are only required to update this quantity once at the beginning of each time increment. The resulting systems of Eqs. (8), (10), and (14), including implementation details for the micro-scale problem at the scale of the fracture opening, are provided in the Supplement.

Figure 5 Depth dependence of the ice sheet temperature and resulting creep coefficient and tensile strength.

The macro-scale governing equations at the glacier length scale (Eqs. 1 and 11) are solved in a monolithic manner. This monolithic solver uses an energy-based convergence criterion [f_u f_p][du dp]^T<ϵ, placing equal importance on the convergence of the nonlinear CZM and the fluid pressure within the crack. Once the solution is converged, the stresses ahead of the crack along a pre-determined path (dashed line in Fig. 2) are calculated. The stresses are compared to the tensile strength to determine fracture propagation. If the crack propagates, a single new interface element is inserted and more iterations of the monolithic solver are performed (including checking for additional fracture propagation) to obtain a solution, where both unknown variables (i.e. displacement and pressure) are solved at the end of each time increment using the updated crevasse length. It is noteworthy that the path along which new interface elements are inserted is pre-determined, allowing the crevasse to only propagate straight down and then split into two basal cracks that can only propagate sideways in a straight line. While the pre-determination of crack path and insertion of cohesive interface elements only between continuum finite elements are limitations of our numerical approach, it is reasonably realistic given the 2D idealisation of the rectangular glacier domain. The requirement that the crack path be known a priori restricts the nucleation of crevasses elsewhere in the domain. Also, we do not model the surface hydrology associated with the formation of supraglacial lakes but rather assume a pre-existing lake with a known depth that intersects this initial crack.

To initialise the simulations, 1 d of time is simulated using time increments of Δt = 10 min, where the crevasse is not allowed to propagate. During this period, the viscous creep alters the stress state from that of a compressible linear-elastic material to that of a nearly incompressible material compatible with viscoelastic deformations. This 1 d initialisation period was long enough for the stresses within the ice to stabilise to the steady state, with further initialisation time not altering these stresses. After this initial time period, the crack is allowed to propagate, and the remainder of the simulation is performed using Δt = 2 s for a duration of 2 h.

Table 1 Material properties used within the simulations.

2.4 Material properties

The properties used to model the ice and rock layers are provided in Table 1, along with parameters for water flowing in the crack. The temperature of the ice sheet is approximated using a temperature profile from Jakobshavn Glacier (Ryser et al., 2014b, a), with the location of this profile ≈ 80 km away from the lake drainage observations used as comparison. This temperature profile, shown in Fig. 5, is directly used within the description of the melting and freezing process. Additionally, the creep coefficient of Glen's law, Eq. (3), is determined based on this temperature profile as (Weertman, 1983; Duddu and Waisman, 2012; Greve et al., 2014)

$$\begin{array}{}\text{(21)}& A=A_{\mathrm{0}}\text{exp}\left({\displaystyle \frac{-Q_{\mathrm{c}}}{R}}\left({\displaystyle \frac{\mathrm{1}}{T}}-{\displaystyle \frac{\mathrm{1}}{T_{\text{ref}}}}\right)\right),\end{array}$$

and the tensile strength of the ice as (Litwin et al., 2012)

$$\begin{array}{}\text{(22)}& f_t=f_{\mathrm{t}\mathrm{0}}-f_{\text{deg}}T.\end{array}$$

These equations use the ice temperature T and the reference temperature T_ref in Kelvin. They further use the activation energy related to the creep within ice, Q_c, the gas constant R, and the temperature degradation rate f_deg. The creep coefficient and tensile strength from these relations are shown in Fig. 5 for the temperature profile used.

3 Results

3.1 Importance of viscoelasticity

The deformations of the ice sheet using a viscoelastic and linear-elastic constitutive model are shown in Fig. 6. If the linear-elastic model is used, crack opening is determined by the balance between the elastic stresses and the pressure within the crevasse, Fig. 6a. This results in limited opening, and the vertical crevasse propagation does not consume a significant volume of lake water. The crevasse propagation rate is governed by this water flow, retaining a pressurised crevasse and reaching the bedrock layer after 13 min. Upon reaching the rock layer, once the pressure within the crevasse is significant enough to start the horizontal crack propagation, the ice sheet will begin to lift up solely due to the water pressure acting at the base of the vertical crevasse. As ice is assumed to be compressible elastic, the higher water pressure causes the crevasse to be wider at its base and more narrow towards the ice surface, owing to the lower water pressure. However, once the ice sheet starts to lift up significantly, the water pressure in the horizontal crack underneath the ice acts in the vertical direction, opposing vertical cryostatic pressure. At this moment, the water pressure required to further propagate the horizontal crack decreases substantially. This causes a significant burst in the propagation, as shown in Fig. 7, but the additional volume created with the horizontal crack leads to the reduction in water pressure at the inlet. As the ice sheet is lifted up, the deformation of the ice sheet resembles that of a floating ice beam with the depth-varying water pressure applying a triangularly distributed load on the vertical crevasse walls. The resultant bending moment causes the opening width at the upper portions of the crevasse to be restricted and eventually forces it to be fully closed, as illustrated by the solid blue line in Fig. 8. This prevents any further fluid from entering the crevasse, causing the only slight increases in pressure due to the freezing of fracture walls and thus eventually stopping the propagation of the horizontal crack at the ice–bedrock interface.

Figure 6 Overpressure of the water within the crevasse relative to the cryostatic ice pressure and deformations after 2 h. Deformations in (a) and (b) are magnified by 1000. Animations showing the deformations over the full simulation duration are provided in the Supplement.

Figure 7 Vertical crevasse depth (dashed lines) and horizontal crack length at the ice–bedrock interface (solid lines) following the definitions from Fig. 1. The vertical dotted line indicates the moment the vertical crevasse reaches the base of the ice sheet, after which point the horizontal crack propagation begins.

Figure 8 Crevasse mouth opening width due to viscous and elastic deformations in the surrounding ice and melting of the crevasse walls. Vertical dotted lines indicate the time when the vertical crevasse reaches the base of the ice sheet. Notably, the opening width continues to increase in the viscoelastic case, whereas it decreases to zero in the linear-elastic case. Crevasse opening due to melting is at least an order of magnitude less than that due to deformation.

In contrast, using the viscoelastic ice sheet rheology produces a crevasse that continuously widens during downwards propagation due to creep deformation induced by the overpressure within the crevasse. This means that part of the water entering the crevasse will need to be retained to compensate for the ever-increasing crevasse volume. As a result, the rate of downwards crevasse propagation is slower compared to the linear-elastic model, which does not account for this additional opening width. The viscoelastic creep strains also cause the opening of the crevasse to be ever-increasing while the crevasse is pressurised, Fig. 8, which allows more water to flow into the crevasse and downwards towards the bed. Once a crevasse penetrates the entire ice thickness and the crack tip reaches the bed, horizontal crack propagation commences, facilitated by the higher stresses due to the incompressible nonlinear viscoelastic response. Furthermore, as the ice sheet begins to lift up, even though the pressure within the crevasse suddenly drops due to the newly created volume, the viscoelastic creep deformations that occurred during the downward propagation allow the vertical crevasse to remain open (note that the crevasse closes at this stage in the linear-elastic model). As the crevasse remains open, water is able to flow through it, reach the bed, and lift up the ice sheet. The sustained transfer of water from the supraglacial lake to the bed is the main control on the rate of basal uplift and horizontal basal crack propagation. As the horizontal basal crack continues to propagate, the rate at which additional crack volume is created due to viscous deformations continues to increase, while the rate of water inflow through the vertical crevasse increases relatively slowly. This causes the pressure at the base of the vertical crevasse to decrease, as shown in Fig. 9. Eventually, this pressure becomes sufficiently low such that further crack propagation is paused, coinciding with a stress state where the viscous deformations allow for a crack volume enhancement equal to the water inflow. As the opening height of the vertical crevasse continues to increase due to viscous deformations even when the horizontal crack is halted, the rate of water inflow slowly starts to exceed the rate of volume increase, allowing the pressure within the horizontal crack to slowly recover. Once this pressure is sufficiently high, the horizontal crack resumes propagation. This alternation in pressure at the crack tip leads to episodic propagation. This process causes the undulations observed in the shape of the basal surface of the ice sheet (see Fig. 6) that are at most only ∼ 10 % of the total opening height and causes the stepped crack length plot in Fig. 7 due to episodic propagation and halting.

Figure 9 Pressure at the crack tip (black) and at the base of the crevasse (red) for the viscoelastic case. Shaded regions correspond to moments when the crack propagation is halted, as visible in Fig. 7.

3.2 Importance of melting

The thermal energies produced by friction, gained due to freezing, and lost due to conduction into the ice are given in Fig. 10. As the bedrock layer is at the melting point of ice (see Fig. 5), the only thermal energy loss due to conduction is through the vertical crevasse walls, whereas in the horizontal crack, no temperature difference exists to conduct heat into the ice. The conductive heat fluxes continue to increase proportionally with crevasse depth while the downward propagation is occurring, until the crevasse reaches the bedrock. After this time, thermal conduction is solely governed by the $t^{-\mathrm{3}/\mathrm{2}}$ dependency, causing the conduction rate to slowly decrease as the ice surrounding the crevasse warms. As this is unrelated to the rheological model used for the ice, the conduction of thermal energy follows the same trend for both models, with the conduction energy loss for the viscoelastic model lagging slightly behind due to the crevasse reaching the base at a slower rate. In contrast, frictional heating shows a strong dependence on the rheology model. In the linear-elastic model case, the smaller crevasse opening combined with the faster propagation causes significantly more frictional heating due to fluid flow at the initial stages. However, as the fluid flow stops at later stages due to the closing of the crevasse mouth, this frictional heating also halts. In the viscoelastic model case, the fluid has a wider crevasse to flow through, but only a limited amount of fluid flows into the horizontal crack at the glacier bed, which initially reduces the total thermal energy generated by friction (i.e. for t < 90 min in Fig. 10). However, because the fluid flow is continuously maintained due to wider crevasse opening, the length over which the fluid is transported increases as the horizontal crack propagates; therefore, the rate of heating increases throughout the simulation, eventually exceeding the heating produced by the linear-elastic model (i.e. for t > 90 min in Fig. 10).

Figure 10 Total thermal energy created due to water–ice friction, lost due to conduction into the ice, and consumed by freezing (positive) or melting (negative). Energies are normalised by the heat ρ_iℒ required to melt 1 m³ of ice, and the values are calculated per unit of out-of-plane length.

The difference between the conduction and frictional heat fluxes dictates the rate of freezing at the crevasse walls. In the linear-elastic model case, when the fluid flow stops the friction heat flux goes to zero, so freezing will occur at a steady rate at the vertical crevasse walls. Eventually, at a certain point in time, the vertical crevasse and the horizontal crack will fully refreeze, although this does not happen within the simulated time of 110 min shown in Fig. 10. In contrast, in the viscoelastic model case, the ever-increasing rate of frictional heating and the slowly reducing rate of conductive losses cause the freezing process to slow down and eventually cause some of the initially frozen crevasse walls to start melting. It should be noted, however, that the thermal energies produced due to friction and lost due to conduction are quite small, only sufficient to freeze ≈ 2 m³ m⁻¹ of ice over the full length of the crevasse (980 m downwards and up to 4 km sideways). This amounts to 1–2 mm of ice building up on the vertical crevasse walls over this time period compared to a crevasse opening width of ≈ 0.5 m, as shown in Fig. 8. This differs from the crevasses observed in reality, where opening widths of several metres are not uncommon. While melting is commonly credited with partly creating such large openings (Andrews et al., 2022), our results allude to processes such as rapid glacial sliding (due to reduced basal friction) occurring after the onset of the lake drainage event that could lead to large crevasse openings.

3.3 Application to rapid lake drainage

We conduct a comparative study with the observations from Das et al. (2008), as shown in Fig. 11. To obtain these changes in the water level from the 2D simulation results (producing fracture inflows in m³ m⁻¹), we assumed the estimates of the lake area as A_lake = 5.6 km² and the out-of-plane length of the crevasse as W_oop = 3.2 km (based on values estimated by Das et al., 2008 and used by Tsai and Rice, 2012), allowing the water level change to be obtained as $\mathrm{\Delta }{h}_{\text{lake}}=Q{W}_{\text{oop}}/A_{\text{lake}}$. As our model predicts horizontal cracks of 1.6 km in each direction compared to the out-of-plane width of W_oop = 3.2 km, assuming plane-strain conditions is reasonable for this case. We calculated fluid inflow rates using the lake height time series in Das et al., which has a 20 min sampling frequency. Due to this sample frequency, the water level changes consist of a linear interpolation between sample points, while the inflow rate is only piecewise continuous.

Figure 11 Vertical uplift of the ice sheet surface over time at set distances from the crevasse (x = 0 m is at the crevasse, and x = 500 m is 500 m to the left/right of it) using the viscoelastic rheology. For comparison, observational data from Das et al. (2008) are included. The vertical dotted line indicates the moment the crack reaches the base of the ice sheet. (Linear-elastic uplift is provided in the Supplement).

Figure 12 Model results using a linear-elastic (blue) and viscoelastic (black) ice rheology are compared to the reference observations (red circles) from Das et al. (2008): (a) lake water level versus time and (b) rate of water level change versus time. The vertical dotted line indicates the moment the vertical crevasse reaches the base of the ice sheet. For all lines, the vertical axis reports the values (in m³ m⁻¹) as directly resulting from simulations on the right and on the left (in m water level change), converted using A_lake = 5.6 km² and W_oop = 3.2 km.

Even though the 2D model considered here is highly idealised and does not represent the complex 3D crevasse and lake geometry in reality, the predicted surface uplift matches well with the measurements of a GPS station located 1–2 km from the crevasse for the first 90 min (Fig. 11). In contrast, the observed and simulated lake water level changes match the observations to a lesser extent, as shown in Fig. 12. One potential explanation for the larger mismatch in the fracture inflow compared to uplift is the observed development of secondary cracks during the hydraulic fracturing process, which could have significantly enhanced fluid transport to the ice sheet base. Additional assumptions for the numerical model are that the ice sheet is pristine (i.e. undamaged) and the ice–rock boundary is initially frozen, neither of which is strictly correct. Within the ice, pre-existing cracks, crevasses, and defects can link to the newly developed hydrofractures, which could significantly enhance the water inflow. Furthermore, fluid flow and movement at the ice–bed interface as it drains downglacier could influence the modelled fluid inflow and ice sheet uplift. These connections are also a potential explanation for the exponential increase in lake drainage, which differs from the behaviour observed within the simulations. One final point of potential mismatch is the conversion between water volumes resulting from our simulations to the lake drainage height reported by Das et al. (2008) (and inversely, from lake water level to volumes from Das et al.). We assume a simplified lake geometry with a constant area as water height decreases instead of taking into account lake bathymetry. This simplification is therefore likely responsible for the model's underestimation of lake water level change, particularly during the later stages of the lake drainage.

4 Discussion

4.1 Influence of rheology

With the linear-elastic rheology, the vertical crevasse is able to propagate through the 980 m thick ice sheet in 13 min, whereas with the elastic–viscoplastic rheology it takes about 25 min. Despite this time difference, only a small amount of lake water (< 100 m³ m⁻¹) is required to drive the fracture to the base. Upon reaching the base and lifting up the ice sheet, the fracture propagation process is able to consume a substantial amount of water (> 1000 ${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$) from the supraglacial lake. In contrast, with the viscoelastic rheology horizontal fracture propagation is episodic, characterised by periods of fast crack growth and stagnation. Specifically, horizontal fracturing stagnates when fluid inflow is fully accommodated by the space created due to viscoelastic creep deformation, limiting any increase in fluid pressure necessary to drive further propagation. However, as creep deformation widens the vertical crevasse, the increased fluid inflow leads to pressure increase and further fracture propagation.

A key difference is that with the viscoelastic rheology model, the method can capture the uplift of the ice sheet due to supraglacial lake drainage, consistent with observations, whereas with the linear-elastic rheology, uplift stagnates after approximately 1 km. This is because in the linear-elastic case, the depth-varying water pressure within the vertical crevasse and the resultant bending moment are sufficient to close the narrow crevasse opening at the top surface, thus preventing additional water from entering the vertical crevasse (see Figs. 12, 8, and 6). Our parametric studies (see Fig. SI1-SI2 in the Supplement) show that the bending effect in the linear-elastic case is independent of ice thickness over shorter timescales (hours) and will always close the crevasse near the ice surface, thereby stopping crevasse propagation.

While we only considered a long crevasse (using plane-strain assumptions), it is expected that these conclusions regarding the role of viscous strains hold if axisymmetric models were used or if the crevasse geometry is more complex; irrespective of the crack geometry, the overpressure within the crevasse will always result in viscous strains enhancing the opening, increasing the water transport towards the glacier bed. High stress concentrations in the area surrounding the crack tip will also still result in a short Maxwell timescale, such that viscous strains occur on timescales comparable to those relevant to hydrofracture. While the horizontal crack growth rate, surface elevation change (uplift), and water level change may be different, the main finding – viscous deformation exerts a much stronger control on hydrofracture propagation compared to thermal effects – would not change.

4.2 Influence of ice temperature

Thermal processes within the crevasse have a limited effect in our simulation studies (see Fig. 8). The change in crevasse opening due to freezing/melting (a few millimetres) is orders of magnitude smaller than the elastic/viscoplastic deformation (tens of centimetres). Consequently, frictional heating-induced melting is unable to prevent crevasse closure in the linear-elastic simulations. On the flip side, conduction loss-driven freezing at the crevasse walls is not enough to close the crevasse opening in the viscoelastic simulations under the thermal conditions considered here. This is not to say that thermal effects are insignificant in the complete process: the creation of supraglacial lakes is driven by melting, and results show that if the surface layers of ice are cold enough, crevasse propagation is halted before it reaches significant depths (see Supplement Sect. S3). We identify an ice temperature threshold of −8 °C at and below which crevasses freeze shut and propagation is prevented. Our results suggest that crevasse propagation models applied to rapid supraglacial lake drainage events in warmer ice regions do not need to account for thermal interactions (causing refreezing) on the short timescales of crevasse development. This finding can greatly simplify the parameterisation of coupled surficial and subglacial hydrology within models. For example, a simple parameterisation could involve identifying regions where ice surface temperatures are above −8 °C and in these regions introduce the influence of surface melt on subglacial hydrology and basal friction.

4.3 Relating to observational data

Advances in technology have been reflected in the production of observational data sets capturing transient lake drainage events (Andrews et al., 2018; Chudley et al., 2019; Das et al., 2008; Doyle et al., 2013; Hoffman et al., 2011; Mejía et al., 2021; Stevens et al., 2015, 2022; Lai et al., 2021). However, models seeking to represent hydraulic fracturing still use analytical linear-elastic fracture mechanics from the 1970s (Weertman, 1971, 1973) and 1990s (Desroches et al., 1994). The computational framework developed here advances the analysis of rapid lake drainage and fracture events. Our model results using a viscoelastic rheology produce realistic timescales for supraglacial lake drainage events on the Greenland Ice Sheet and are closely tied to lake volume. Based on the total inflow at 120 min (about 1300 m³ m⁻¹ in Fig. 12a), we estimate that for the North Lake in 2006 (Das et al., 2008), 0.0042 km³ of water drained in 2 h through a crevasse with an out-of-plane length of W_oop = 3.2 km. Notably, the inflow rate of lake drainage (black line in Fig. 12b) stabilised to 0.0032–0.0035 km³ h⁻¹ (1000–1100 ${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$) over the second hour of the simulation. Using this average inflow rate and assuming the same out-of-plane crevasse length reported in 2006, we now compare our results to observations of North Lake drainage from 2011 to 2013 (Stevens et al., 2015). Lake volumes and drainage duration over these 3 years were 0.0077 km³ over 3 h in 2011, 0.0077 km³ over 5 h in 2012, and 0.0057 km³ over 5 h in 2013. Assuming that the inflow rate at 2 h will remain constant for the next hour, Fig. 12a, we obtain a total inflow volume of 0.0074 km³ after 3 h, which matches the reported values for 2011, albeit overestimating the lake volumes for the 5 h drainage events from 2012 and 2013. The discrepancy in estimating lake volumes in 2012 and 2013 is likely related to the assumption of out-of-plane crevasse length, which leads to smaller inflow rates compared to 2011. Estimates of crevasse volume obtained from a network inversion filter for 2012 and 2013 are smaller than that for 2011 (Stevens et al., 2015).

We can also compare the rate of crevasse opening obtained within simulations to observed strain rates and slip velocities, as one of the assumptions in the computational model is that no basal slip occurs before basal uplifting occurs. Ice surface velocities that have been reported in the area surrounding North Lake are around 91 m yr⁻¹ (Ryser et al., 2014b), or around 1 cm h⁻¹ (in the absence of ongoing lake drainage events). Das et al. (2008) reported speed-ups of this velocity by 300 % during lake drainage due to reductions in basal friction. However, even at this increased rate, the enhancement to the crevasse width would be limited to only 6 cm compared to the 50 cm crevasse opening created due to the water pressure within the crevasse inducing local elastic and viscous deformations. For basal motion to have a significant influence on crevasse width, lake water would need to reach the bed and become distributed over a large enough area to influence basal water pressures, lift up the ice, and cause enhanced sliding. Because of these prerequisites, peak basal motion would have a delayed influence on crevasse width, only occurring after the crevasse has opened and transferred water to the bed. As this occurs well outside the considered time span within this work (2 h), the effect of basal motion on crevasse width can indeed be ignored in our study.

4.4 Future of the Greenland Ice Sheet

As supraglacial lakes expand inland and occupy higher-elevation areas in response to atmospheric warming (Howat et al., 2013), understanding the influence of these lake drainage events on subglacial hydrology is becoming more important. We foresee this model's application to cases where describing the fracture mechanics in high detail is desired, such as in constraining the formation and drainage of subglacial flood waves following rapid supraglacial lake drainage. While the model presented here is restricted in scope to only solve for a single hydrofracture and the resulting uplift, results can be coupled with data sets and other process-based models to fully investigate ice dynamic response to rapid supraglacial lake drainage. For example, crevasse formation models (e.g. Hoffman et al., 2018) or remote sensing data can constrain the timing of supraglacial lake drainage, with the later also defining pre-drainage lake conditions. Results from our subsequent model runs can be fed back into large-scale low-resolution glaciological models or smaller flood wave propagation models (e.g. Lai et al., 2021) to inform changes in basal conditions both locally and downstream along the resulting subglacial flood wave. Indeed, when interpreting GPS-derived ice motion of supraglacial lake drainage events, utilising a high-resolution model such as this would enable one to parse out the ice motion produced by fracture propagation and that produced by subglacial processes such as frictional sliding. Large volumes of water injected into the bed following rapid lake drainage create a subglacial flood wave that modulates ice velocity and can alter the subglacial drainage system as it moves downglacier. These flood waves can connect and drain previously hydrologically isolated regions of the subglacial drainage system that control minimum ice velocities, thereby having lasting effects on hydrodynamic coupling evident throughout the remainder of the melt season (Mejía et al., 2021). The area of the ice sheet's bed that can be modified is unconstrained but can extend tens of kilometres downglacier. While our model does not incorporate a basal drainage component or simulate the subglacial flood wave produced by the lake drainage event, these complex processes can be investigated in the future by coupling our model with a subglacial hydrology model to more fully assess the ice-dynamic and hydrological consequences of rapid lake drainage. The model presented here therefore provides not only a mathematical framework for interpreting in situ observations but also a mechanism to simulate detailed subglacial flooding, which can provide more accuracy when inferring subglacial transmissivity and establishing initial conditions of the subglacial flood wave produced as water drains downglacier after the lake drainage event.

4.5 Advancing ice sheet modelling

Crevasses play a role in two of the most poorly understood glaciological processes, subglacial hydrology and iceberg calving; consequently, they are also poorly represented in ice sheet models. The impact of subglacial hydrology on basal motion (Bueler and van Pelt, 2015) remains poorly or rarely represented in numerical ice sheet models, leading to potentially large model uncertainty (Aschwanden et al., 2021). Data assimilation techniques typically use observations of GrIS geometry and assume ice rheology is known (based on ice temperature, englacial properties, crystal orientation, and impurities) so that basal friction can be treated as the only unknown field parameter (Goelzer et al., 2017). While such techniques can provide modelled velocities close to observations in diagnostic simulations, nonphysical responses may be predicted in prognostic simulations (Seroussi et al., 2011). Existing ice sheet models (Lipscomb et al., 2019) do not incorporate changes in basal friction due to supraglacial lake drainage events, which can severely limit their applicability to future warmer scenarios. As we scale up our computational framework to glacier-scale simulations, we intend to use it develop simpler parameterisations linking surface hydrology to subglacial hydrology. Future studies should consider mapping the combinations of strain rates and surface temperatures that relate to lake formation and drainage, and introduce fluid volume and pressure in fractures as inputs into subglacial hydrology and basal friction models. Although this is not such a simple task, our two-scale modelling framework is a first step towards exploring the interactions between the thermal, hydraulic, and mechanical process controlling GrIS flow and fracture.

5 Conclusions

We have developed a two-scale computational method able to capture the hydraulic fracturing process responsible for rapid supraglacial lake drainage in high temporal and spatial resolution. By resolving ice sheet deformation in 2D, capturing the water within the crevasse in 1D, and using a sub-grid-scale formulation based on analytic expressions for thermal conduction, melting, and frictional heating, the relevant mechanisms surrounding crevasses in ice sheets are all captured even though they are relevant on drastically different length scales. By separating these mechanisms across scales, we achieve a highly (computationally) efficient scheme capable of capturing the dynamic processes of crevasse formation and subsequent ice uplifting. While no direct verification studies have been performed with this two-scale formulation for glacial fracturing, previous applications have demonstrated its accuracy for benchmark hydrofracture cases in rock media, albeit without the thermal model (Hageman and de Borst, 2022, 2019).

Our novel modelling framework allows us to explore the vertical and horizontal fracture propagation during rapid lake drainage events; discern the role of thermal effects; and more importantly, evaluate the repercussions of ignoring viscoelasticity even on the short timescales associated with these events. We find that viscoelastic creep deformation has a significant effect on the horizontal fracture propagation at the glacier bed and the vertical crevasse opening width, allowing for the greater fluid flow required for rapid lake drainage. Our findings reasonably illustrate the timescales and the order of magnitude of water volumes involved in the rapid drainage of supraglacial lakes, as observed by Das et al. (2008). This study demonstrates the utility of our modelling framework for developing a parameterisation of hydraulic fracture in a large-scale model of the Greenland Ice Sheet and ultimately for understanding the implications for sea level rise.