the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Ice viscosity governs hydraulic fracture that causes rapid drainage of supraglacial lakes
Tim Hageman
Jessica Mejía
Ravindra Duddu
Emilio MartínezPañeda
Fullthickness crevasses can transport water from the glacier surface to the bedrock where high water pressures can open kilometrelong cracks along the basal interface, which can accelerate glacier flow. We present a first computational modelling study that describes timedependent fracture propagation in an idealised glacier causing rapid supraglacial lake drainage. A novel twoscale numerical method is developed to capture the elastic and viscoelastic deformations of ice along with crevasse propagation. The fluidconserving 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 rapidlakedrainage 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.
 Article
(3668 KB)  Fulltext XML

Supplement
(7777 KB)  BibTeX
 EndNote
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 physicsbased 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 waterfilled 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 linearelasticfracturemechanicsbased (LEFMbased) models were used to estimate the depth of waterfilled crevasses, namely, the Nye zero stress (Jezek, 1984; Benn et al., 2007; Nick et al., 2010), dislocationbased LEFM (Weertman, 1971), and stressintensityfactorbased LEFM models (Smith, 1976; Van Der Veen, 1998). Recently, stressintensityfactorbased LEFM models for estimating waterfilled crevasse depth have been proposed using finiteelementbased (Jimenez and Duddu, 2018) and boundaryelementbased methods (Zarrinderakht et al., 2022). Alternatively, continuum damage mechanics and phase field models for waterfilled 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 linearelastic 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 longerterm processes such as the slow propagation of crevasses over weeks to months (Poinar et al., 2017) or moulin formation posthydrofracture (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 largescale mass and momentum balances (Rice et al., 2015; Andrews et al., 2022). Within numerical simulations, reducedorder models have been used to include the interactions between basal water and basal friction (Pimentel and Flowers, 2011). Some studies utilised the finiteelement 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 elasticviscoplastic 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 linearelastic 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.
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 outofplane 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 outofplane 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 outofplane 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 planestrain 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 planestrain 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.
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 zeroverticaldisplacement 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 noslip 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 frozenbed 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 leakoff from the base of the ice sheet. Imposing no leakoff from the waterfilled crevasse deviates from realistic conditions, where water may slowly seep into the bedrock layer at the base. This leakoff 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 finiteelement formulation is used to accurately capture the fracturing and subsequent uplift process. This formulation includes the linearelastic and nonlinear viscoelastic deformations of the ice.^{1} This allows for the inclusion of both shorterterm (seconds to minutes) elastic deformations and longerterm (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 smallscale 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 planestrain 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 historydependent 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 fluidfilled fracture is described through its pressure p, and its history is included through the timesincefracture t_{0} 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 nonNewtonian viscous liquid for longterm 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:
where the density ρ_{π} corresponds to either the bedrock density ρ_{r} or the ice density ρ_{i}; $\mathit{g}=[\mathrm{0}\phantom{\rule{0.25em}{0ex}}\mathrm{9.81}\phantom{\rule{0.125em}{0ex}}\mathrm{m}\phantom{\rule{0.125em}{0ex}}{\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}\phantom{\rule{0.25em}{0ex}}{\mathit{\epsilon}}_{yy}\phantom{\rule{0.25em}{0ex}}{\mathit{\epsilon}}_{zz}\phantom{\rule{0.25em}{0ex}}{\mathit{\epsilon}}_{xy}{]}^{\mathrm{T}}=\mathbf{L}\mathit{u}$. The rock is assumed to be a linearelastic solid, whereas ice is assumed to be a Norton–Hoff type elasticviscoplastic material. Thus, the stress in ice is defined by the linearelastic strain, which is obtained by offsetting the total strain with the viscous strain history. The constitutive law for ice and rock is defined by
where the isotropic linearelastic 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
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 waterfilled crevasse to be localised close to the crevasse, so the temperature field and the creep coefficient A are timeinvariant within the bulk of the ice sheet.
Throughout this work, ice is considered a viscoelastic solid, which encompasses the special case of a linearelastic 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 linearelastic 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 linearelastic 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)
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 temperaturedependent 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 linearelastic 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.
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
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 (tensiondriven fracture). When modelling sharp cracks within a linearelastic material, the stresses at the crack tip become singular. Propagation criteria in linearelastic 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 elastoplastic 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 finiteelement 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
This cohesive traction depends on the fracture release energy G_{c} and the tensile strength of the ice f_{t}, unlike linearelastic 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(Hy)$. 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 noslip 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 glacierslidinginduced crevasse widening would be significant if the basal friction were weak.
We use an extrinsictype CZM implementation, in which the interface elements are inserted dynamically into the finiteelement mesh as the crack propagates. As a result, the cohesive zone model is only applied postcracking, while the notyetfractured 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 tractionseparation 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 tractionfree 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 twoscale scheme: resolving the fluid mass conservation within the fracture at the macroscale (implemented through a standard finiteelement scheme) while resolving the thermal processes as microscale effects (numerically resolved on a perintegrationpoint basis).
2.2.1 Pressuredriven 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–Stricklertype flow law (Gauckler, 1867; Strickler, 1981), resulting in a volume flux q produced by turbulent flow as (Tsai and Rice, 2010, 2012)
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}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{h}}^{\mathrm{1}}$ (see results section and Fig. 12), providing a Reynolds number for this flow as $\mathit{Re}={\mathit{\rho}}_{\mathrm{w}}\leftq\right/{\mathit{\mu}}_{\mathrm{w}}\approx $ 3 × 10^{5}, well above the typical transition point from laminar to turbulent flow (Re>𝒪(10^{3})). 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 finetuning these friction parameters.
The total aperture of the fracture is obtained based on the displacement jump across the discontinuity and melting layer thickness as
and the friction factor f related to the (constant) wall roughness k_{wall} and the reference friction factor f_{0} as (Strickler, 1981)
Together, Eqs. (7) and (9) allow for the fluid transport within the crack to be described explicitly as
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):
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 nearnegligible in practice (especially with the properties used in this paper), it provides a dampinglike 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 penaltylike approach:
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 postprocessing purposes, this flux is additionally integrated over time to obtain the total volume of water flowing into the crevasse as
As a 2D domain is considered, this produces the total volume of fluid that has entered the crevasse per unit of outofplane width (given in m^{3} m^{−1}). 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 outofplane width of the crevasse.
Equation (11) is solved to capture the macroscale 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 microscale behaviour on a pointbypoint basis. This novel and efficient twoscale approach is needed to account for the dynamic feedback between fluid flux and wall melting, which precludes a closedform 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 nearfreezing 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
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 selfreinforcing: 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 flowproduced heat flux is given by (Andrews et al., 2022)
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:
using the heat capacity c_{p} and the thermal conductivity k, and using the surfacenormal coordinate η. Because the ice–water boundary is consistently at T_{w} = 0 °C, an analytic solution for this equation is given by (White, 2006)
using erfc to indicate the complementary error function, T_{∞} the temperature of the ice and t_{0} 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
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):
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 timesincefracture 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 microscale from the mechanical problem, which is considered the macroscale. 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 twoscale 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 centimetresized 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 subgridscale 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:
with the first term composed of physical constants and the second term of casedependent 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 waterfilled 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 finiteelement 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 postprocessing (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 timeimplicit 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 microscale problem at the scale of the fracture opening, are provided in the Supplement.
The macroscale governing equations at the glacier length scale (Eqs. 1 and 11) are solved in a monolithic manner. This monolithic solver uses an energybased 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 predetermined 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 predetermined, 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 predetermination 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 preexisting 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 linearelastic 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.
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)
and the tensile strength of the ice as (Litwin et al., 2012)
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.1 Importance of viscoelasticity
The deformations of the ice sheet using a viscoelastic and linearelastic constitutive model are shown in Fig. 6. If the linearelastic 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 depthvarying 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.
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 everincreasing crevasse volume. As a result, the rate of downwards crevasse propagation is slower compared to the linearelastic model, which does not account for this additional opening width. The viscoelastic creep strains also cause the opening of the crevasse to be everincreasing 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 linearelastic 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.
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 linearelastic 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 linearelastic model (i.e. for t > 90 min in Fig. 10).
The difference between the conduction and frictional heat fluxes dictates the rate of freezing at the crevasse walls. In the linearelastic 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 everincreasing 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^{3} m^{−1} 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^{3} m^{−1}), we assumed the estimates of the lake area as A_{lake} = 5.6 km^{2} and the outofplane 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 outofplane width of W_{oop} = 3.2 km, assuming planestrain 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.
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, preexisting 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.1 Influence of rheology
With the linearelastic 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^{3} m^{−1}) 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}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\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 linearelastic rheology, uplift stagnates after approximately 1 km. This is because in the linearelastic case, the depthvarying 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. SI1SI2 in the Supplement) show that the bending effect in the linearelastic 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 planestrain 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 heatinginduced melting is unable to prevent crevasse closure in the linearelastic simulations. On the flip side, conduction lossdriven 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 linearelastic 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^{3} m^{−1} in Fig. 12a), we estimate that for the North Lake in 2006 (Das et al., 2008), 0.0042 km^{3} of water drained in 2 h through a crevasse with an outofplane 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^{3} h^{−1} (1000–1100 ${\mathrm{m}}^{\mathrm{3}}\phantom{\rule{0.125em}{0ex}}{\mathrm{m}}^{\mathrm{1}}\phantom{\rule{0.125em}{0ex}}{\mathrm{h}}^{\mathrm{1}}$) over the second hour of the simulation. Using this average inflow rate and assuming the same outofplane 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^{3} over 3 h in 2011, 0.0077 km^{3} over 5 h in 2012, and 0.0057 km^{3} 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^{3} 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 outofplane 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^{−1} (Ryser et al., 2014b), or around 1 cm h^{−1} (in the absence of ongoing lake drainage events). Das et al. (2008) reported speedups 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 higherelevation 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 processbased 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 predrainage lake conditions. Results from our subsequent model runs can be fed back into largescale lowresolution 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 GPSderived ice motion of supraglacial lake drainage events, utilising a highresolution 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 icedynamic 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 glacierscale 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 twoscale modelling framework is a first step towards exploring the interactions between the thermal, hydraulic, and mechanical process controlling GrIS flow and fracture.
We have developed a twoscale 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 subgridscale 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 twoscale 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 largescale model of the Greenland Ice Sheet and ultimately for understanding the implications for sea level rise.
The MATLAB code described in the methods and used to generate the results for this research is available from https://doi.org/10.5281/zenodo.13629447 (Hageman, 2024), where documentation for this code and postprocessing scripts are also provided.
Animations at https://doi.org/10.5446/68350 (Hageman et al., 2024) show the deformations and pressure over time for Fig. 6.
(1) Implementation details of the finiteelement scheme are given. (2) A parametric study regarding the effect of ice sheet thickness on the closing of crevasses using linearelastic and viscoelastic rheologies is given. (3) A parametric study on the role of ice sheet temperatures on stagnating crevasse propagation due to freezing is shown. (4) A figure shows the uplift over time using a linearelastic rheology (the linearelastic counterpart of Fig. 11). (5) Animations show the deformations and pressure over time for Fig. 6. The supplement related to this article is available online at: https://doi.org/10.5194/tc1839912024supplement.
TH: software, formal analysis, data curation, investigation, writing – original draft, visualisation, and methodology. JM: formal analysis, investigation, and writing – original draft. RD: investigation, writing – original draft, writing – review and editing, conceptualisation, and methodology. EMP: writing – review and editing, and conceptualisation.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
Tim Hageman acknowledges financial support through the research fellowship scheme of the Royal Commission for the Exhibition of 1851. Emilio MartínezPañeda acknowledges financial support from UKRI's Future Leaders Fellowship programme (grant no. MR/V024124/1). Ravindra Duddu acknowledges funding support from the NSF Office of Polar Programs via CAREER grant no. PLR1847173, and the Royal Society via the International Exchanges programme, grant no. IES/R1/211032. Jessica Mejía acknowledges support by the Heising–Simons Foundation, grant no. 20201910. The authors also acknowledge computational resources and support provided by the Imperial College London research computing service (2024).
This research has been supported by the Royal Commission for the Exhibition of 1851, the UK Research and Innovation programme (grant no. MR/V024124/1), the Office of Polar Programs (grant no. CAREER PLR1847173), the Royal Society (grant no. IES/R1/211032), and the Heising–Simons Foundation (grant no. 20201910).
This paper was edited by Nanna Bjørnholt Karlsson and reviewed by Douglas Benn and one anonymous referee.
Andrews, L. C., Hoffman, M. J., Neumann, T. A., and Catania, G. A.: Seasonal Evolution of the Subglacial Hydrologic System Modified by Supraglacial Lake Drainage in Western Greenland, J. Geophys. Res.Earth, 123, 1479–1496, https://doi.org/10.1029/2017JF004585, 2018. a, b
Andrews, L. C., Poinar, K., and Trunz, C.: Controls on Greenland moulin geometry and evolution from the Moulin Shape model, The Cryosphere, 16, 2421–2448, https://doi.org/10.5194/tc1624212022, 2022. a, b, c, d
Aschwanden, A., Bartholomaus, T. C., Brinkerhoff, D. J., and Truffer, M.: Brief communication: A roadmap towards credible projections of ice sheet contribution to sea level, The Cryosphere, 15, 5705–5715, https://doi.org/10.5194/tc1557052021, 2021. a
Bamber, J. L., Oppenheimer, M., Kopp, R. E., Aspinall, W. P., and Cooke, R. M.: Ice sheet contributions to future sealevel rise from structured expert judgment, P. Natl. Acad. Sci. USA, 166, 11195–11200, https://doi.org/10.1073/pnas.1817205116, 2019. a
Benn, D. I., Hulton, N. R., and Mottram, R. H.: 'Calving laws', 'sliding laws' and the stability of tidewater glaciers, Ann. Glaciol., 46, 123–130, https://doi.org/10.3189/172756407782871161, 2007. a
Bevis, M., Harig, C., Khan, S. A., Brown, A., Simons, F. J., Willis, M. J., Fettweis, X., van den Broeke, M. R., Madsen, F. B., Kendrick, E., Caccamise, D. J., van Dam, T., Knudsen, P., and Nylen, T.: Accelerating changes in ice mass within Greenland, and the ice sheet's sensitivity to atmospheric forcing, P. Natl. Acad. Sci. USA, 116, 1934–1939, https://doi.org/10.1073/pnas.1806562116, 2019. a
Boone, T. J. and Ingraffea, A. R.: A numerical procedure for simulation of hydraulicallydriven fracture propagation in poroelastic media, Int. J. Numer. Anal. Met., 14, 27–47, https://doi.org/10.1002/nag.1610140103, 1990. a
Bueler, E. and van Pelt, W.: Massconserving subglacial hydrology in the Parallel Ice Sheet Model version 0.6, Geosci. Model Dev., 8, 1613–1635, https://doi.org/10.5194/gmd816132015, 2015. a
Carrier, B. and Granet, S.: Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model, Eng. Fract. Mech., 79, 312–328, https://doi.org/10.1016/j.engfracmech.2011.11.012, 2012. a
Christoffersen, P., Bougamont, M., Hubbard, A., Doyle, S. H., Grigsby, S., and Pettersson, R.: Cascading lake drainage on the Greenland Ice Sheet triggered by tensile shock and fracture, Nat. Commun., 9, 1–12, https://doi.org/10.1038/s41467018034208, 2018. a
Chudley, T. R., Christoffersen, P., Doyle, S. H., Bougamont, M., Schoonman, C. M., Hubbard, B., and James, M. R.: Supraglacial lake drainage at a fastflowing Greenlandic outlet glacier, P. Natl. Acad. Sci. USA, 116, 25468–25477, https://doi.org/10.1073/PNAS.1913685116, 2019. a, b, c
Clayton, T., Duddu, R., Siegert, M., and MartínezPañeda, E.: A stressbased porodamage phase field model for hydrofracturing of creeping glaciers and ice shelves, Eng. Fract. Mech., 272, 108693, https://doi.org/10.1016/J.ENGFRACMECH.2022.108693, 2022. a, b
Crawford, A. J., Benn, D. I., Todd, J., Åström, J. A., Bassis, J. N., and Zwinger, T.: Marine IceCliff Instability Modeling Shows MixedMode IceCliff Failure and Yields Calving Rate Parameterization, Nat. Commun., 12, 1–9, https://doi.org/10.1038/s41467021230707, 2021. a
Das, S. B., Joughin, I., Behn, M. D., Howat, I. M., King, M. A., Lizarralde, D., and Bhatia, M. P.: Fracture propagation to the base of the Greenland ice sheet during supraglacial lake drainage, Science, 320, 778–781, https://doi.org/10.1126/SCIENCE.1153360, 2008. a, b, c, d, e, f, g, h, i, j, k, l
de Borst, R.: Fluid flow in fractured and fracturing porous media: A unified view, Mech. Res. Commun., 80, 47–57, https://doi.org/10.1016/j.mechrescom.2016.05.004, 2017. a
de Fleurian, B., Gagliardini, O., Zwinger, T., Durand, G., Le Meur, E., Mair, D., and Råback, P.: A double continuum hydrological model for glacier applications, The Cryosphere, 8, 137–153, https://doi.org/10.5194/tc81372014, 2014. a
Desroches, J., Detournay, E., Lenoach, B., Papanastasiou, P., Pearson, J. R. A., Thiercelin, M., and Cheng, A.: The crack tip region in hydraulic fracturing, P. Roy. Soc. Lond. AMat., 447, 39–48, 1994. a
Doyle, S. H., Hubbard, A. L., Dow, C. F., Jones, G. A., Fitzpatrick, A., Gusmeroli, A., Kulessa, B., Lindback, K., Pettersson, R., and Box, J. E.: Ice tectonic deformation during the rapid in situ drainage of a supraglacial lake on the Greenland Ice Sheet, The Cryosphere, 7, 129–140, https://doi.org/10.5194/tc71292013, 2013. a, b, c
Duddu, R. and Waisman, H.: A temperature dependent creep damage model for polycrystalline ice, Mech. Mater., 46, 23–41, https://doi.org/10.1016/J.MECHMAT.2011.11.007, 2012. a, b
Duddu, R., Jiménez, S., and Bassis, J.: A nonlocal continuum porodamage mechanics model for hydrofracturing of surface crevasses in grounded glaciers, J. Glaciol., 66, 415–429, 2020. a
Gauckler, P.: Etudes Théoriques et Pratiques sur l'Ecoulement et le Mouvement des Eaux, C. R. Acad. Sci. Paris, 64, 818–822, 1867. a
Ghosh, G., Duddu, R., and Annavarapu, C.: A stabilized finite element method for enforcing stiff anisotropic cohesive laws using interface elements, Computer Method. Appl. M., 348, 1013–1038, 2019. a
Glen, J. W.: The creep of polycrystalline ice, P. Roy. Soc. Lond. AMat., 228, 519–538, https://doi.org/10.1098/rspa.1955.0066, 1955. a, b
Goelzer, H., Robinson, A., Seroussi, H., and van de Wal, R. S.: Recent Progress in Greenland Ice Sheet Modelling, Current Climate Change Reports, 3, 291–302, https://doi.org/10.1007/s406410170073y, 2017. a
Goelzer, H., Nowicki, S., Payne, A., Larour, E., Seroussi, H., Lipscomb, W. H., Gregory, J., AbeOuchi, A., Shepherd, A., Simon, E., Agosta, C., Alexander, P., Aschwanden, A., Barthel, A., Calov, R., Chambers, C., Choi, Y., Cuzzone, J., Dumas, C., Edwards, T., Felikson, D., Fettweis, X., Golledge, N. R., Greve, R., Humbert, A., Huybrechts, P., Le clec'h, S., Lee, V., Leguy, G., Little, C., Lowry, D. P., Morlighem, M., Nias, I., Quiquet, A., Rückamp, M., Schlegel, N.J., Slater, D. A., Smith, R. S., Straneo, F., Tarasov, L., van de Wal, R., and van den Broeke, M.: The future sealevel contribution of the Greenland ice sheet: a multimodel ensemble study of ISMIP6, The Cryosphere, 14, 3071–3096, https://doi.org/10.5194/tc1430712020, 2020. a
Greve, R., Zwinger, T., and Gong, Y.: On the pressure dependence of the rate factor in Glen's flow law, J. Glaciol., 60, 397–399, https://doi.org/10.3189/2014JOG14J019, 2014. a
Hageman, T.: THageman/MATLAB_IceHydroFrac: Final version (final), Zenodo [code], https://doi.org/10.5281/zenodo.13629447, 2024. a
Hageman, T. and de Borst, R.: Flow of nonNewtonian fluids in fractured porous media: Isogeometric vs standard finite element discretisation, Int. J. Numer. Anal. Met., 43, 2020–2037, https://doi.org/10.1002/nag.2948, 2019. a
Hageman, T. and de Borst, R.: A refined twoscale model for Newtonian and nonNewtonian fluids in fractured poroelastic media, J. Comput. Phys., 441, 110424, https://doi.org/10.1016/j.jcp.2021.110424, 2021. a, b
Hageman, T. and de Borst, R.: Direct simulation vs subgrid scale modelling of fluid flow in fractured or fracturing porous media, Computat. Geosci., 26, 503–515, https://doi.org/10.1007/s10596022101386, 2022. a
Hageman, T., Pervaiz Fathima, K. M., and de Borst, R.: Isogeometric analysis of fracture propagation in saturated porous media due to a pressurised nonNewtonian fluid, Comput. Geotech., 112, 272–283, https://doi.org/10.1016/j.compgeo.2019.04.030, 2019. a
Hageman, T., Mejía, J., Duddu, R., and MartínezPañeda, E.: Ice viscosity governs hydraulic fracture causing rapid drainage of supraglacial lakes, TIB AVPortal [video], https://doi.org/10.5446/68350, 2024. a
Hewitt, I. J., Schoof, C., and Werder, M. A.: Flotation and free surface flow in a model for subglacial drainage. Part 2. Channel flow, J. Fluid Mech., 702, 157–187, https://doi.org/10.1017/JFM.2012.166, 2012. a
Hofer, S., Lang, C., Amory, C., Kittel, C., Delhasse, A., Tedstone, A., and Fettweis, X.: Greater Greenland Ice Sheet contribution to global sea level rise in CMIP6, Nat. Commun., 11, 6289, https://doi.org/10.1038/s41467020200118, 2020. a
Hoffman, M. J., Catania, G. A., Neumann, T. A., Andrews, L. C., and Rumrill, J. A.: Links between acceleration, melting, and supraglacial lake drainage of the western Greenland Ice Sheet, J. Geophys. Res.Earth, 116, 1–16, https://doi.org/10.1029/2010JF001934, 2011. a, b
Hoffman, M. J., Perego, M., Andrews, L. C., Catania, G. A., Price, S. F., Lüthi, M. P., Neumann, T. A., and Johnson, J. V.: Widespread Moulin Formation During Supraglacial Lake Drainages in Greenland, Geophys. Res. Lett., 45, 778–788, https://doi.org/10.1002/2017GL075659, 2018. a
Howat, I. M., de la Peña, S., van Angelen, J. H., Lenaerts, J. T. M., and van den Broeke, M. R.: Brief Communication “Expansion of meltwater lakes on the Greenland Ice Sheet”, The Cryosphere, 7, 201–204, https://doi.org/10.5194/tc72012013, 2013. a
Imperial College London research computing service: https://doi.org/10.14469/hpc/2232, 2024. a
Jellinek, H. H. G. and Brill, R.: Viscoelastic Properties of Ice, J. Appl. Phys., 27, 1198–1209, https://doi.org/10.1063/1.1722231, 1956. a
Jezek, K. C.: A modified theory of bottom crevasses used as a means for measuring the buttressing effect of ice shelves on inland ice sheets, J. Geophys. Res.Sol. Ea., 89, 1925–1931, https://doi.org/10.1029/JB089IB03P01925, 1984. a
Jimenez, S. and Duddu, R.: On the evaluation of the stress intensity factor in calving models using linear elastic fracture mechanics, J. Glaciol., 64, 759–770, 2018. a
Krawczynski, M. J., Behn, M. D., Das, S. B., and Joughin, I.: Constraints on the lake volume required for hydrofracture through ice sheets, Geophys. Res. Lett., 36, 10501, https://doi.org/10.1029/2008GL036765, 2009. a
Lai, C. Y., Stevens, L. A., Chase, D. L., Creyts, T. T., Behn, M. D., Das, S. B., and Stone, H. A.: Hydraulic Transmissivity Inferred from IceSheet Relaxation Following Greenland Supraglacial Lake Drainages, Nat. Commun., 12, 1–10, https://doi.org/10.1038/s41467021241866, 2021. a, b, c
Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res.Earth, 117, 1022, https://doi.org/10.1029/2011JF002140, 2012. a
Liang, Y. L., Colgan, W., Lv, Q., Steffen, K., Abdalati, W., Stroeve, J., Gallaher, D., and Bayou, N.: A decadal investigation of supraglacial lakes in West Greenland using a fully automatic detection and tracking algorithm, Remote Sens. Environ., 123, 127–138, https://doi.org/10.1016/J.RSE.2012.03.020, 2012. a
Lipscomb, W. H., Price, S. F., Hoffman, M. J., Leguy, G. R., Bennett, A. R., Bradley, S. L., Evans, K. J., Fyke, J. G., Kennedy, J. H., Perego, M., Ranken, D. M., Sacks, W. J., Salinger, A. G., Vargo, L. J., and Worley, P. H.: Description and evaluation of the Community Ice Sheet Model (CISM) v2.1, Geosci. Model Dev., 12, 387–424, https://doi.org/10.5194/gmd123872019, 2019. a, b
Litwin, K. L., Zygielbaum, B. R., Polito, P. J., Sklar, L. S., and Collins, G. C.: Influence of temperature, composition, and grain size on the tensile failure of water ice: Implications for erosion on Titan, J. Geophys. Res.Planets, 117, 8013, https://doi.org/10.1029/2012JE004101, 2012. a
Mejía, J. Z., Gulley, J. D., Trunz, C., Covington, M. D., Bartholomaus, T. C., Xie, S., and Dixon, T.: Isolated cavities dominate Greenland Ice Sheet dynamic response to lake drainage, Geophys. Res. Lett., 48, 1–11, https://doi.org/10.1029/2021gl094762, 2021. a, b, c
Nick, F. M., Van Der Veen, C. J., Vieli, A., and Benn, D. I.: A physically based calving model applied to marine outlet glaciers and implications for the glacier dynamics, J. Glaciol., 56, 781–794, https://doi.org/10.3189/002214310794457344, 2010. a
Pattyn, F.: The Paradigm Shift in Antarctic Ice Sheet Modelling, Nat. Commun., 9, 1–3, https://doi.org/10.1038/s4146701805003z, 2018. a
Pimentel, S. and Flowers, G. E.: A numerical study of hydrologically driven glacier dynamics and subglacial flooding, P. Roy. Soc. AMath. Phy., 467, 537–558, https://doi.org/10.1098/RSPA.2010.0211, 2011. a
Poinar, K., Joughin, I., Lilien, D., Brucker, L., Kehrl, L., and Nowicki, S.: Drainage of southeast Greenland firn aquifer water through crevasses to the bed, Front. Earth Sci., 5, https://doi.org/10.3389/feart.2017.00005, 2017. a
Réthoré, J., de Borst, R., and Abellan, M.A.: A twoscale approach for fluid flow in fractured porous media, Int. J. Numer. Meth. Eng., 71, 780–800, https://doi.org/10.1002/nme.1962, 2006. a
Rice, J. R., Tsai, V. C., Fernandes, M. C., and Platt, J. D.: Time scale for rapid draining of a surficial lake into the Greenland ice sheet, J. Appl. Mech.T. ASME, 82, 071001, https://doi.org/10.1115/1.4030325, 2015. a
Ryser, C., Lüthi, M. P., Andrews, L. C., Catania, G. A., Funk, M., Hawley, R., Hoffman, M., and Neumann, T. A.: Caterpillarlike ice motion in the ablation zone of the Greenland ice sheet, J. Geophys. Res.Earth, 119, 2258–2271, https://doi.org/10.1002/2013JF003067, 2014a. a
Ryser, C., Luthi, M. P., Andrews, L. C., Hoffman, M. J., Catania, G. A., Hawley, R. L., Neumann, T. A., and Kristensen, S. S.: Sustained high basal motion of the Greenland ice sheet revealed by borehole deformation, J. Glaciol., 60, 647–660, https://doi.org/10.3189/2014JOG13J196, 2014b. a, b
Selmes, N., Murray, T., and James, T. D.: Fast draining lakes on the Greenland Ice Sheet, Geophys. Res. Lett., 38, 15501, https://doi.org/10.1029/2011GL047872, 2011. a, b
Seroussi, H., Morlighem, M., Rignot, E., Larour, E., Aubry, D., Ben Dhia, H., and Kristensen, S. S.: Ice Flux Divergence Anomalies on 79north Glacier, Greenland, Geophys. Res. Lett., 38, https://doi.org/10.1029/2011GL047338, 2011. a
Shannon, S. R., Payne, A. J., Bartholomew, I. D., Broeke, M. R. V. D., Edwards, T. L., Fettweis, X., Gagliardini, O., GilletChaulet, F., Goelzer, H., Hoffman, M. J., Huybrechts, P., Mair, D. W., Nienow, P. W., Perego, M., Price, S. F., Smeets, C. J. P., Sole, A. J., Wal, R. S. V. D., and Zwinger, T.: Enhanced basal lubrication and the contribution of the Greenland ice sheet to future sealevel rise, P. Natl. Acad. Sci. USA, 110, 14156–14161, https://doi.org/10.1073/pnas.1212647110, 2013. a
Smith, L. C., Chu, V. W., Yang, K., Gleason, C. J., Pitcher, L. H., Rennermalm, A. K., Legleiter, C. J., Behar, A. E., Overstreet, B. T., Moustafa, S. E., Tedesco, M., Forster, R. R., LeWinter, A. L., Finnegan, D. C., Sheng, Y., and Balog, J.: Efficient meltwater drainage through supraglacial streams and rivers on the southwest Greenland ice sheet, P. Natl. Acad. Sci. USA, 112, 1001–1006, https://doi.org/10.1073/pnas.1413024112, 2015. a
Smith, R. A.: The Application of Fracture Mechanics to the Problem of Crevasse Penetration, J. Glaciol., 17, 223–228, https://doi.org/10.3189/S0022143000013563, 1976. a
Stevens, L. A., Behn, M. D., McGuire, J. J., Das, S. B., Joughin, I., Herring, T., Shean, D. E., and King, M. A.: Greenland supraglacial lake drainages triggered by hydrologically induced basal slip, Nature, 522, 73–76, https://doi.org/10.1038/nature14480, 2015. a, b, c, d
Stevens, L. A., Nettles, M., Davis, J. L., Creyts, T. T., Kingslake, J., Hewitt, I. J., and Stubblefield, A.: Tidewaterglacier response to supraglacial lake drainage, Nat. Commun., 13, 6065, https://doi.org/10.1038/s41467022337632, 2022. a
Strickler, A.: Contributions to the Question of a Velocity Formula and Roughness Data for Streams, Channels and Closed Pipelines, Tech. rep., W. M. Keck Laboratory of Hydraulics and Water Resources, 1981. a, b
Sun, X., Duddu, R., and Hirshikesh: A porodamage phase field model for hydrofracturing of glacier crevasses, Extreme Mechanics Letters, 45, 101277, https://doi.org/10.1016/j.eml.2021.101277, 2021. a
Trunz, C., Covington, M. D., Poinar, K., Andrews, L. C., Mejia, J., and Gulley, J.: Modeling the Influence of Moulin Shape on Subglacial Hydrology, J. Geophys. Res.Earth, 127, e2022JF006674, https://doi.org/10.1029/2022JF006674, 2022. a
Tsai, V. C. and Rice, J. R.: A model for turbulent hydraulic fracture and application to crack propagation at glacier beds, J. Geophys. Res.Earth, 115, 3007, https://doi.org/10.1029/2009JF001474, 2010. a, b, c
Tsai, V. C. and Rice, J. R.: Modeling Turbulent Hydraulic Fracture Near a Free Surface, J. Appl. Mech., 79, 031003, https://doi.org/10.1115/1.4005879, 2012. a, b, c
Van Der Veen, C. J.: Fracture mechanics approach to penetration of surface crevasses on glaciers, Cold Regions Science and Technology, 27, 31–47, https://doi.org/10.1016/S0165232X(97)000220, 1998. a
van der Veen, C. J.: Fracture propagation as means of rapidly transferring surface meltwater to the base of glaciers, Geophys. Res. Lett., 34, L01501, https://doi.org/10.1029/2006GL028385, 2007. a, b
Weertman, J.: Theory of waterfilled crevasses in glaciers applied to vertical magma transport beneath oceanic ridges, J. Geophys. Res., 76, 1171–1183, https://doi.org/10.1029/JB076I005P01171, 1971. a, b
Weertman, J.: Can a water filled crevasse reach the bottom surface of a glacier?, IASH Publ., 95, 139–145, https://doi.org/10.1017/CBO9781107415324.004, 1973. a, b
Weertman, J.: Creep deformation of ice, Annu. Rev. Earth Pl. Sc., 11, 215–240, https://doi.org/10.1146/annurev.ea.11.050183.001243, 1983. a
White, F. M.: Viscous Fluid Flow, 3 edn., McGrawHill, NewYork, ISBN 007124493X, 2006. a, b
Witherspoon, P. A., Wang, J. S. Y., Iwai, K., and Gale, J. E.: Validity of Cubic Law for fluid flow in a deformable rock fracture, Water Resour. Res., 16, 1016–1024, https://doi.org/10.1029/WR016i006p01016, 1980. a
Zarrinderakht, M., Schoof, C., and Peirce, A.: The effect of hydrology and crevasse wall contact on calving, The Cryosphere, 16, 4491–4512, https://doi.org/10.5194/tc1644912022, 2022. a, b
The term “viscoelastic” is used to refer to the combination of reversible elastic deformations and irreversible viscous/plastic deformations resulting from the use of Hooke's and Glen's laws. While this terminology is more common in the nonNewtonian fluids/rheology and glaciology community, in the solid mechanics community these models are referred to as viscoplastic models, namely the Norton–Hoff and Bingham–Maxwell models.