Vulnerability of firn to hydrofracture: poromechanics modeling

Abstract

On the Greenland Ice Sheet, hydrofracture connects the supraglacial and subglacial hydrologic systems, coupling surface runoff dynamics and ice velocity. In recent decades, the growth of low-permeability ice slabs in the wet snow zone has expanded Greenland's runoff zone, but observations suggest that surface-to-bed connections are rare, because meltwater drains through crevasses into the porous firn beneath ice slabs. However, there is little quantitative evidence confirming the absence of surface-to-bed fracture propagation. Here, we use poromechanics to investigate whether water-filled crevasses in ice slabs can propagate vertically through an underlying porous firn layer. Based on numerical simulations, we develop an analytical estimate of the water injection-induced effective stress in the firn given the water level in the crevasse, ice slab thickness, and firn properties. We find that the firn layer substantially reduces the system's vulnerability to hydrofracture because much of the hydrostatic stress is accommodated by a change in pore pressure, rather than being transmitted to the solid skeleton. This result suggests that surface-to-bed hydrofracture will not occur in ice slab regions until all pore space proximal to the initial flaw has been filled with solid ice.

1. Introduction

Over the last two decades, around 55% of mass loss from the Greenland Ice Sheet has come from the runoff of surface meltwater, with the remainder driven by ice dynamics (van den Broeke and others, 2009; Mouginot and others, 2019). Passive microwave observations and regional climate models also show a long-term increase in the area of the ice sheet experiencing surface melt, the maximum elevation of where melting occurs, and the total length of the annual melt season (Fettweis and others, 2011; Colosio and others, 2021). Therefore, understanding how much and how quickly surface meltwater can be transported through the supraglacial and englacial hydrologic systems and how those systems are evolving with time is critical for assessing current and future sea level contributions from the Greenland Ice Sheet.

Water transport processes vary significantly across the ice sheet. In the bare ice ablation zone, surface meltwater flows efficiently over the impermeable ice surface in streams or river and forms lakes in closed basins (Smith and others, 2015; Yang and Smith, 2016). Particularly in Southwest Greenland, this supraglacial system is connected to the ice sheet bed through fractures, moulins, and rapid lake drainage events, and most melt eventually enters the subglacial system (Das and others, 2008; Dow and others, 2015; Koziol and others, 2017; Andrews and others, 2018; Hoffman and others, 2018; Lai and others, 2021; Poinar and Andrews, 2021). These englacial transport pathways are primarily formed by hydrofracture (Stevens and others, 2015; Poinar and others, 2017) and lead to a coupling between surface melting and ice dynamics, where meltwater delivery to the bed can cause transient, seasonal increases in ice velocity that may temporarily increase ice discharge (Zwally and others, 2002; Schoof, 2010; Moon and others, 2014).

In contrast, in the accumulation zone, meltwater percolates in the porous near-surface firn layer, where it may refreeze locally (Harper and others, 2012; Machguth and others, 2016) or be stored in buried liquid water aquifers (Forster and others, 2014). These processes buffer runoff and prevent water from reaching the subglacial system as long as pore space remains for storage. The processes by which the percolation zone may transition to a bare ice ablation zone under persistent atmospheric warming are not yet fully understood, particularly the timescales over which the hydrologic system evolves from local retention to efficient runoff and drainage.

The development of multi-meter thick ice slabs in the near-surface firn of the wet snow zone appears to be one mechanism by which the ice sheet may transition rapidly from retention to runoff. These continuous, low-permeability layers of refrozen ice block vertical percolation and allow water to flow freely over the surface, despite the presence of a relict porous firn layer at depth (MacFerrin and others, 2019; Tedstone and Machguth, 2022). This produces a surface hydrologic network that qualitatively resembles that of the bare ice ablation zone (Yang and Smith, 2016; Tedstone and Machguth, 2022). However, similar to the data shown in Figure 1, Culberg and others (2022) presented observational evidence from Northwest Greenland that when surface meltwater drains into surface crevasses in ice slabs, it is stored in the underlying relict firn layer where it refreezes into ‘ice blobs’ and does not reach the ice sheet bed. Therefore, they suggest that meltwater in ice slab regions is unlikely to influence the subglacial system, and that surface-to-bed connections may not develop until all local firn pore space is depleted (Culberg and others, 2022).

Figure 1. Ice-penetrating radar observations of water storage beneath an ice slab in Northwest Greenland. (a) Radar observations in 2016 show an ice blob that has refrozen in the porous firn beneath the ice slab. (b) An idealized schematic of an ice blob, after Culberg and others (2022), where the initial crevasse that allowed for drainage is thought to be relatively shallow and to have healed as the ice blob refroze.

While Culberg and others (2022) present strong evidence for large-scale water storage in ice blobs, their evidence that fractures do not propagate through the firn layer and that no water drains deeper into the englacial system is at best circumstantial. The ice-penetrating radar data they present do not directly resolve any surface crevasses (see Fig. 1 for example) and the ice blobs were imaged anywhere from 6 months to 5 years after the initial drainage (Culberg and others, 2022), by which time any deep englacial pathways would have closed by creep. The authors speculated that leak-off into the permeable firn might reduce the crack tip stress enough to prevent unstable propagation, based on literature from the hydraulic fracking community (Culberg and others, 2022). However, to date there have been no investigations of hydrofracture processes in firn, leaving no direct quantitative support for this assumption. In principle, it is possible that ice blobs represent leak-off from the walls of much deeper crevasses that were at one time fully filled with water. Even if ice blobs only form when fractures cannot propagate through the firn, there is currently no way to investigate whether these drainage dynamics are particular to Northwest Greenland, or if this response should generalize to the entire Greenland ice sheet or potentially even to Antarctic ice shelves.

Current models of ice sheet hydrofracture are not suited to investigating this question. The most common approach is to use linear elastic fracture mechanics (LEFM). For example, Lai and others (2020) treated the effects of a near-surface firn layer on dry fracture propagation on Antarctic ice shelves. They assumed that due to its lower density, the presence of firn leads to a lower overburden stress and reduced viscosity. However, their analysis for firn focuses only on water-free crevasse. For analysis of hydrofracture in firn, leakage of water into the firn needs to be considered, whereas LEFM analyses typically assume that the medium is impermeable and incompressible. In fact, work in other fields on the hydraulic fracturing of permeable reservoir suggests that it is the compressibility and permeability of the medium that can lead to resilience to hydrofracture (Bunger and others, 2005; Detournay, 2016; Meng and others, 2020, 2022, 2023; Chen and others, 2022). The observations of Culberg and others (2022) also highlight that even if a crevasse is filled with water, the presence of firn should instead limit fracture propagation. Therefore, the compressibility and permeability of firn need to be considered.

To address this challenge, here we develop a poromechanical model to predict the effective stress in the firn layer beneath a water-filled fracture in an ice slab. This approach can describe both fluid flow out of the fracture and the solid-fluid coupling that impacts stresses. Based on these simulations, we propose an analytical model for the maximum effective stress in the firn layer for both constant water pressure and constant injection rate conditions. We then apply this model to assess the vulnerability of Greenland's ice slab regions to hydrofracture and analyze the behavior of the system as a function of water availability, ice slab thickness, and the mechanical and hydraulic properties of the firn.

2. Methods

In regions with high velocity gradients, dry surface fractures may form in ice slabs. If meltwater flows into these crevasses, they may continue to propagate until they reach the underneath permeable firn layer as shown in Figure 2. The meltwater then penetrates into the firn layer either by infiltration or fracturing. Previous research has used two-phase continuum models to study meltwater flow through snow without considering flow-induced deformation in the porous snow layer (Meyer and Hewitt, 2017; Moure and others, 2023). We develop a two-dimensional, poroelastic continuum model to quantify the stress and pressure changes in the firn layer during meltwater penetration (Biot, 1941; Wang, 2000; Coussy, 2004). Here, we consider two scenarios of water infiltration into the porous firn layer:

1. 1. Constant pressure boundary condition : a constant water height (H _w) in the surface crevasse.

2. 2. Constant flow rate boundary condition: a constant water injection velocity (V _inj) at the crevasse tip.

Figure 2. A schematic showing stresses acting on the porous firn layer during water infiltration from the crevasse.

When stress is applied to porous media, part of the stress is transmitted through the pore fluid and part of the stress is transmitted through the solid skeleton. Effective stress (σ′)–the fraction of the total stress (σ) that is transmitted through the solid skeleton–controls the mechanical behavior of porous media (Terzaghi, 1943). We assume the porous firn layer to be an isotropic, linear elastic continuum. Figure 2 shows the stresses acting on the firn layer from initial stresses and hydrostatic water pressure. We assume plane-strain condition for this 2D model ($\epsilon _{yy} = 0,\; \, {\partial \over \partial y} = 0$). To rationalize the crossover from infiltration to fracturing regimes quantitatively, we adopt a fracturing criterion for cohesive porous media: the horizontal tensile effective stress (σ′_xx) should exceed the material tensile strength (σ′_t) to generate vertical fractures (Wang, 2000; Coussy, 2004). The criterion has been verified experimentally with hydrostone (Haimson and Fairhurst, 1969) and granular packs made of polyurethane (Meng and others, 2023). The fracture criterion at the crevasse tip (point A in Fig. 2) is written as follows:

(1)$$\sigma'_{xx} = \sigma'_{xx, 0} + \delta\sigma'_{xx}\geq\sigma'_t,\; $$

where σ′_xx,0 is the initial horizontal effective stress before water infiltration, δσ′_xx is the infiltration-induced effective stress, and σ′_t is the tensile strength of firn that ranges from 45 kPa to 1 MPa with the density increasing from 500 kg m⁻³ to 900 kg m⁻³ (Petrovic, 2003; Shimaki and Arakawa, 2021).

2.1 Initial effective stresses before water infiltration

Before water penetrates into the surface crevasses, the firn layer has atmospheric pore pressure (p ₀ = 0), and the effective stress is equal to the total stress. We first derive the expression for the initial effective stress in the firn layer without the presence of water, σ′_xx,0. Terzaghi's effective stress tensor σ′ is the portion of the stress supported through deformation of the solid skeleton, and here we adopt the convention of tension being positive. The stress–strain constitutive equation for the linear elastic firn layer is:

(2)$${\boldsymbol \sigma'} = {3K\nu\over 1 + \nu}\epsilon_{kk}{\boldsymbol I} + {3K( 1-2\nu) \over 1 + \nu}{\boldsymbol \epsilon},\; $$

where K, ν, ${\boldsymbol \epsilon }$ are the drained bulk modulus, the drained Poisson ratio of the firn layer (Biot, 1941; Terzaghi, 1943), and the strain tensor, respectively. The constitutive equations for plane-strain are obtained by inserting the constraint that $\epsilon _{yy} = 0$ into Eqn. (2) and noting that $\epsilon _{kk} = \epsilon _{xx} + \epsilon _{zz}$:

(3)$$\sigma_{xx, 0}' = {3K\nu\over 1 + \nu}( \epsilon_{xx, 0} + \epsilon_{zz, 0}) + {3K( 1-2\nu) \over 1 + \nu}\epsilon_{xx, 0},\; $$

(4)$$\sigma_{zz, 0}' = {3K\nu\over 1 + \nu}( \epsilon_{xx, 0} + \epsilon_{zz, 0}) + {3K( 1-2\nu) \over 1 + \nu}\epsilon_{zz, 0}.$$

Solving Eqn. (4) for $\epsilon _{zz, 0}$ and substituting into Eqn. (3) yields

(5)$$\sigma'_{xx, 0} = {\nu\over 1-\nu}\sigma'_{zz, 0} + {3K( 1-2\nu) \over ( 1 + \nu) ( 1-\nu) }\epsilon_{xx, 0}.$$

Eqn. (5) gives the initial effective lateral stresses at the crevasse tip (point A in Fig. 2):

(6)$$\sigma'_{xx, 0} = -{\nu\over 1-\nu}\rho_i gH_i + {3K( 1-2\nu) \over ( 1 + \nu) ( 1-\nu) }\epsilon_{xx, 0},\; $$

where ρ _i, H _i are the density and height of the ice slab above the firn layer, respectively. By assuming linear elasticity for the porous firn layer under plane-strain conditions, the initial horizontal effective stress is partitioned into the lithostatic stress and the extensional/compressional stress induced by the horizontal strain due to ice flow.

2.2 Water infiltration into the porous firn layer

Following Coussy (2004), the poroelasticity equation states that

(7)$$\delta{\boldsymbol \sigma} = \delta{\boldsymbol \sigma'}-b\delta p{\boldsymbol I},\; $$

where δ σ, δ σ′, δp are the infiltration-induced changes in the Cauchy total stress tensor, the effective stress tensor, and the pore pressure, respectively, and b ∈ [0,  1] is the Biot coefficient of the porous medium. We then derive the expression for the infiltration-induced horizontal effective stress, δσ′_xx. Figure 2 shows the stresses acting on the firn layer with water injection through a crevasse with an opening width of 2L _crev and either a constant water height in the crevasse (H _w) or constant water injection velocity (V _inj). To quantify the stresses and pressure changes during the water infiltration into the dry cohesive firn layer, we develop a two-dimensional, two-phase poroelastic continuum model. In the following, we present the extension of Biot's theory to two-phase flow (Jha and Juanes, 2014; Bjørnarå and others, 2016), where we consider small deformations.

2.2.1 Geomechanical equations

Under quasi-static conditions, the balance of linear momentum leads to the equation of equilibrium as follows:

(8)$$\nabla\cdot( {\boldsymbol \sigma}_0 + \delta{\boldsymbol \sigma}) + ( \rho_{b, 0} + \delta\rho_b) {\boldsymbol g} = {\bf 0},\; $$

where g is gravitational acceleration, σ₀ is the initial total stress before water infiltration, δ σ is the infiltration-induced change in the total stress, and ρ _b,0 is the initial bulk density before water infiltration, ρ _b,0 = (1 − ϕ)ρ _s + ϕρ _a. The infiltration-induced change in the bulk density for the solid-fluid system is δρ _b = (ρ _w − ρ _a)ϕS _w (Bjørnarå and others, 2016), where ρ _s is the solid ice density, ϕ is the porosity, and $\rho _\alpha$ and $S_\alpha \in [ 0,\; \, 1]$ are the density and saturation of the fluid phase α (water w or air a), respectively.

We obtain the differential form by subtracting static initial conditions from Eqn. (8):

(9)$$\nabla\cdot( \delta{\boldsymbol \sigma}) + ( \delta\rho_b) {\boldsymbol g} = {\bf 0}.$$

The 2D stress balance equation becomes:

(10)$$\matrix{{\partial( \delta\sigma_{xx}) \over \partial x} + {\partial( \delta\sigma_{zx}) \over \partial z} = 0,\; \quad{\rm in\, x\, direction,\; }\cr {\partial( \delta\sigma_{xz}) \over \partial x} + {\partial( \delta\sigma_{zz}) \over \partial z} + ( \delta\rho_b) g = 0,\; \quad{\rm in\, z\, direction.} }$$

The infiltration-induced strain tensor is defined as $\delta {\boldsymbol \epsilon } = {1\over 2}[ \nabla {\boldsymbol u} + ( \nabla {\boldsymbol u}) ^T]$, where u = [u,  y,  w] is the infiltration-induced displacement vector, and u,  y, and w are the displacements in x, y, z directions, respectively. For 2D deformation in plane-strain condition, the strains are written as:

(11)$$\eqalign{& \delta\epsilon_{xx} = {\partial u\over \partial x},\; \quad\delta\epsilon_{zz} = {\partial w\over \partial z},\; \quad\delta\epsilon_{xz} = {1\over 2}\left({\partial u\over \partial z} + {\partial w\over \partial x}\right),\; \\ & \quad\delta\epsilon_{kk} = \delta\epsilon_{xx} + \delta\epsilon_{zz}.}$$

Using equations (7), (2), and (11), the stress balance equation (10) can be expressed as a function of infiltration-induced displacements u(x,  z,  t), w(x,  z,  t), the infiltration-induced change in the pore pressure δp(x,  z,  t), and the saturation of the water phase S _w(x,  z,  t).

2.2.2 Fluid flow equations

For the two-phase immiscible flow system, the conservation of fluid mass can be written as follows:

(12)$${\partial( \phi\rho_\alpha S_\alpha) \over \partial t} + \nabla\cdot( \rho_\alpha\phi S_\alpha {\boldsymbol v_\alpha}) = 0.$$

The phase velocity ${\boldsymbol v_\alpha }$ is related to the Darcy flux ${\boldsymbol q_\alpha }$ in a deforming medium by the following relation:

(13)$${\boldsymbol q_\alpha} = \phi S_\alpha( {\boldsymbol v_\alpha}-{\boldsymbol v_s}) = -{k_0\over \eta_\alpha}k_{r\alpha}( \nabla ( \delta p_\alpha) -\rho_\alpha{\boldsymbol g}) ,\; $$

where v_s is the velocity of the solid phase, k ₀ is the intrinsic permeability of the porous firn layer, g is the gravity vector, and $\eta _\alpha$, $k_{r\alpha } = k_{r\alpha }( S_\alpha )$ and $\delta p_\alpha$ are the dynamic viscosity, relative permeability, and infiltration-induced fluid pressure for phase α (water w or air a), respectively. As the initial fluid pressure in the dry firn layer is p ₀ = 0 before water infiltration, the total fluid pressure equals to infiltration-induced fluid pressure, $p_\alpha = \delta p_\alpha$. Since capillary pressure is negligible here, (p _c = p _w − p _a = 0) the two phases have the same fluid pressure p. The relative permeability functions are given as Corey-type power law functions (Meyer and Hewitt, 2017; Moure and others, 2023):

(14)$$k_{rw} = S_w^{a_w} \quad {\rm and} \quad k_{ra} = ( 1-S_w) ^{a_a},\; $$

where the fitting parameters are the exponents a _w = 3 and a _a = 2 (Bjørnarå and others, 2016).

The mass conservation equation for the solid phase is given below:

(15)$${\partial[ \rho_s( 1-\phi) ] \over \partial t} + \nabla\cdot[ \rho_s( 1-\phi) {\boldsymbol v_s}] = 0.$$

Assuming isothermal conditions and using the equation of state for the solid, the following expression for the change in porosity is obtained (Lewis and Schrefler, 1999):

(16)$${\partial\phi\over \partial t} = ( b-\phi) \left(c_s{\partial ( \delta p) \over \partial t} + \nabla\cdot{\boldsymbol v_s}\right),\; $$

where c _s is the compressibility of the solid phase, b ∈ [0,  1] is the Biot coefficient of the porous medium. Combining equations (13) and (16), we expand equation (12) and arrive at the governing equation for two phase fluids (see details of algebra in the Supplementary Material):

(17)$${\phi{\partial S_\alpha\over \partial t} + S_\alpha\left(b{\partial( \delta\epsilon_{kk}) \over \partial t} + \left(\phi c_\alpha + ( b-\phi) c_s\right){\partial ( \delta p) \over \partial t}\right) + \nabla\cdot{\boldsymbol q}_\alpha = 0,\;} $$

where $c_\alpha$ (Pa⁻¹) is the compressibility of fluid phase and $\delta \epsilon _{kk}$ is the infiltration-induced volumetric strain of the solid phase.

Adding equation (17) for water and air phases, and imposing that S _a + S _w ≡ 1 for the porous firn layer, we obtain the pressure diffusion equation:

(18)$$b{\partial( \delta\epsilon_{kk}) \over \partial t} + {1\over M}{\partial ( \delta p) \over \partial t} + \nabla\cdot{\boldsymbol q}_t = 0,\; $$

where q_t = q_w + q_a is the total Darcy flux for water and air phases. The Biot modulus of the porous firn, M, is related to fluid and firn properties as ${1\over M} = \phi S_w c_w + \phi ( 1-S_w) c_a + ( b-\phi ) c_s$, where c _w, c _a are the water and air compressibility, respectively (Coussy, 2004).

2.2.3 Summary of governing equations

We use a 2D, two-phase poroelastic continuum model to solve the infiltration-induced stress and pressure changes. The model has four governing equations, two derived from conservation of fluid mass [Eqn. (18) for the water–air fluid mixture and Eqn. (17) for the water phase] and two derived from conservation of linear momentum [Eqn. (10)]. After applying the constitutive law [Eqn. (2)] and the definition of the effective stress [Eqn. (7)], stresses in Eqn. (10) are expressed in terms of displacements and pore pressure. The model solves the time evolution of four unknowns: (1) pore pressure field δp(x,  z,  t); (2) water saturation field S _w(x,  z,  t); (3) horizontal displacement field u(x,  z,  t), and (4) vertical displacement field w(x,  z,  t) of the porous firn layer. The governing equations are summarized and written in x,  z coordinates as follows:

(19)$$\matrix{b{\partial( \delta\epsilon_{kk}) \over \partial t} + {1\over M}{\partial ( \delta p) \over \partial t}-k_0{\partial\over \partial x}\left(\left({k_{rw}\over \eta_w} + {k_{ra}\over \eta_a}\right){\partial ( \delta p) \over \partial x}\right)\cr \quad- k_0{\partial\over \partial z}\left(\left({k_{rw}\over \eta_w} + {k_{ra}\over \eta_a}\right){\partial ( \delta p) \over \partial z}-\left({\rho_w k_{rw}\over \eta_w} + {\rho_a k_{ra}\over \eta_a}\right)g\right) = 0,\; }$$

(20)$$\eqalignno{{\matrix{\phi{\partial S_w\over \partial t} + S_w\bigg(b{\partial ( \delta\epsilon_{kk}) \over \partial t} + \bigg(\phi c_w + ( b-\phi) c_s\bigg){\partial ( \delta p) \over \partial t}\bigg)-{k_0\over \eta_w}{\partial\over \partial x}\bigg(k_{rw}{\partial ( \delta p) \over \partial x}\bigg)\cr \quad- {k_0\over \eta_w}{\partial\over \partial z}\bigg(k_{rw}\big({\partial ( \delta p) \over \partial z}-\rho_wg\bigg)\bigg) = 0,\; }}}$$

(21)$$\eqalign{\textstyle{\partial\over \partial x}\bigg({3K\nu\over 1 + \nu}\delta\epsilon_{kk} + {3K( 1-2\nu) \over 1 + \nu}\delta\epsilon_{xx}-b\delta p\bigg) \cr+ \textstyle{\partial\over \partial z}\bigg({3K( 1-2\nu) \over 1 + \nu}\delta\epsilon_{xz}\bigg) = 0,\;}$$

(22)$$\eqalignno{& \textstyle{\partial\over \partial x}\bigg({3K( 1-2\nu) \over 1 + \nu}\delta\epsilon_{xz}\bigg) + {\partial\over \partial z}\bigg({3K\nu\over 1 + \nu}\delta\epsilon_{kk} + {3K( 1-2\nu) \over 1 + \nu}\delta\epsilon_{zz}-b\delta p\bigg) \cr& + ( \delta\rho_b) g = 0,\;}$$

where ${\delta \epsilon _{xx} = {\partial u\over \partial x},\; \, \delta \epsilon _{zz} = {\partial w\over \partial z},\; \, \delta \epsilon _{xz} = {1\over 2}( {\partial u\over \partial z} + {\partial w\over \partial x}) ,\; \, \delta \epsilon _{kk} = \delta \epsilon _{xx} + \delta \epsilon _{zz}}$.

Solving the four coupled governing equations [Eqns. (19), (20), (21), (22)], we obtain the spatiotemporal evolution of the saturation, displacements and pressure field. To quantify the vulnerability of firn to hydrofractures based on Eqn. (1), the model outputs the infiltration-induced change of horizontal effective stress as follows:

(23)$$\delta\sigma'_{xx} = {3K\nu\over 1 + \nu}\delta\epsilon_{kk} + {3K( 1-2\nu) \over 1 + \nu}\delta\epsilon_{xx}.$$

2.2.4 Initial and boundary conditions

The model solves the infiltration-induced pressure and stress changes in the porous firn layer, where 0 ≤ x ≤ L, and 0 ≤ z ≤ H as shown in Figure 2. Water flows into the crevasse, the tip of which has an opening L _crev. Water infiltrates into the porous firn layer via the crevasse tip at either a constant height, H _w, or a constant velocity, V _inj. We initialize the model by specifying u(x,  z,  0) = w(x,  z,  0) = δp(x,  z,  0) = 0. The water saturation is zero everywhere except at the crevasse tip, where it is kept at S _w = 1 as follows:

(24)$${S_w( x\leq L_{\rm{crev}},\; 0,\; t) = 1,\; S_w( x> L_{\rm{crev}},\; 0,\; t) = S_w( x,\; z> 0,\; 0) = 0.}$$

We consider the stress (or displacement) and pressure (or flow rate) boundary conditions on the domain area. On the left boundary (x = 0), axis of symmetry requires that ${\partial \over \partial x} = 0$, and horizontal displacement equals zero. On the right boundary, we assume it is far from the crevasse tip and unaffected by the infiltration. On the bottom boundary where the firn layer touches the impermeable, rigid ice column, the displacements and vertical water flow rate are zero. On the top boundary, when the water surface height exceeds the ice slab height (e.g. when a lake overlies the ice slab), the lake depth adds to the initial lithostatic stresses. Otherwise the overlying ice slab provides constant lithostatic stresses. The vertical water flowrate is zero everywhere except at the crevasses, where either δp = ρ _w gH _w or q _w,z = V _inj. The boundary conditions are summarized as follows:

(25)$$\matrix{& u\vert _{x = 0} = {\partial w\over \partial x}\vert _{x = 0} = {\partial p\over \partial x}\vert _{x = 0} = {\partial S_w\over \partial x}\vert _{x = 0} = 0,\; \cr & u\vert _{x = L} = w\vert _{x = L} = \delta p\vert _{x = L} = 0,\; \cr & u\vert _{z = H} = w\vert _{z = H} = q_{w, z}\vert _{z = H} = 0,\; \cr & \delta\sigma'_{xx}\vert _{z = 0} = 0,\; \cr & \delta\sigma'_{zz}\vert _{z = 0} = 0,\; \cr & q_{w, z}\vert _{z = 0}( x> L_{\rm{crev}}) = 0,\; \cr & \delta p\vert _{z = 0}( x\leq L_{\rm{crev}}) = \rho_w gH_w\quad {\rm or} \quad q_{w, z}\vert _{z = 0}( x\leq L_{\rm{crev}}) = V_{\rm{inj}}. }$$

2.2.5 Model parameters

The four poroelastic constants in the model are the drained bulk modulus K, the drained Poisson ratio ν, the Biot coefficient b, and the Biot modulus M of the firn layer. We calculate the Biot coefficient from the relationship $b = 1-{K\over K_s}$ (Coussy, 2004), where K _s is the bulk modulus of the ice grain. The Biot modulus of the porous firn, M, is related to fluid and firn properties as ${1\over M} = \phi S_w c_w + \phi ( 1-S_w) c_a + ( b-\phi ) c_s$ (Coussy, 2004), which is a spatiotemporal variable as water penetrates into the firn layer. A summary of the modeling parameters is given in Table 1.

Table 1. Modeling parameters for the 2D, two-phase poroelastic continuum model

2.2.6 Numerical implementation

We use a finite volume numerical scheme to solve the four coupled governing equations [Eqns. (19), (20), (21), (22)]. We sequentially solve time-discretized equations within a timestep. We update the flow and mechanics simultaneously [Eqns. (19), (21), (22)] using implicit time discretization, and separately from the saturation by treating the saturation as a fixed variable. Then we solve the water transport equation [Eqn. (20)] with prescribed pressure and displacement fields. By the end of each timestep, the four unknown quantities are all updated. See the Supplementary Material for the convergence and mesh independence analysis. The modeling results reach convergence with the mesh size dx=dz=0.5 m, which is adopted for all the simulation presented here.

3. Modeling results

3.1 Spatiotemporal evolution of pressures and stresses

Figures 3 and 4 compare our modeling results for water infiltration with two different boundary conditions at the crevasse tip: constant pressure (H _w = 10 m), and constant flow rate (V _inj = 0.05 m s⁻¹). In both cases, water infiltrates into the porous firn layer due to the pressure gradient and gravity, resulting in an elliptical-like shape for the water saturation profile as shown in Figure 3a. Figure 3b presents the temporal evolution of the pore pressure field, highlighting the pressure diffusion within the water phase. The viscous dissipation is constrained within a certain depth, below which the water flow becomes purely gravity-driven. Figure 3c shows the spatiotemporal evolution of the infiltration-induced horizontal effective stress change (δσ′_xx(x,  z,  t)). At the crevasse tip, the firn is under the maximum tensile effective stress, which makes it the most vulnerable place for hydrofracturing (see the fracture criterion in Eqn. (1)).

Figure 3. Modeling results for the water infiltration with H _w = 10 m (left panel), V _inj = 0.05 m s⁻¹ (right panel) in the domain 0 < x < L, 0 < z < H, where L = H = 30 m. A sequence of snapshots shows the spatiotemporal evolution of (a) water saturation field, S _w(x,  z,  t), (b) pore pressure field, δp(x,  z,  t), and (c) infiltration-induced horizontal effective stress change, δσ′_xx(x,  z,  t). Infiltration time t = 20 s, 40 s, 180 s from snapshot (i), (ii) to (iii).

Figure 4. Modeling results for the water infiltration with H _w = 10 m (black line) and V _inj = 0.05 m s⁻¹ (blue line). Time evolution of (a) injection pressure P _inj(t) at the crevasse tip, and (b) infiltration-induced horizontal effective stress change at the crevasse tip δσ′_xx(t). We use their maximum values (δp, δσ′_xx,max) to evaluate the vulnerability of the porous firn layer to hydrofracturing. The markers indicate times for the snapshots shown in Figure 3: t = 20 s, 40 s, 180 s in sequence.

Figure 4 presents the temporal evolution of the infiltration-induced pore pressure (P _inj(t)) and the effective stress (δσ′_xx(t)) at the crevasse tip. For the constant velocity injection condition, Figure 4a shows that it takes some time for the injection pressure to build up, and thus the water invading front is slightly delayed compared with the constant pressure condition. With either constant pressure or constant flow rate as the boundary condition, Figure 4b shows that the tensile effective stress at the crevasse tip experiences a logarithmic-like growth in time with a fast increase in the first 30 seconds. To avoid boundary effects, we terminate the simulation when water infiltrates into half of the domain depth, and take the corresponding pressure (δp) and tensile effective stress at the crevasse tip (δσ′_xx,max) for the following scaling analysis. To confirm that the effective stress does not increase noticeably as infiltration keeps going, we have extended the simulation time and find that the effective stress value we choose effectively captures its maximum value developed during infiltration (see Supplementary Material).

3.2 Analytical model

3.2.1 Scaling between δσ′_xx,max and δp

To check whether fracture initiates in the porous firn layer during meltwater infiltration, we focus on the stress and pressure changes at the crevasse tip, which has been shown to be the most vulnerable place. To implement the fracture criterion in Eqn. (1) more efficiently, we develop a scaling relationship between δσ′_xx,max and δp, which dictates how much of the infiltration-induced pore pressure change is transmitted to the solid skeleton. We recall the poroelasticity equation for the effective stress [Eqn. (7)], and propose a scaling law for the infiltration-induced horizontal effective stress change at the crevasse tip as follows:

(26)$$\delta\sigma'_{xx, \rm{max}} = \delta\sigma_{xx, \rm{max}} + b\delta p = -\gamma ( b\delta p) + ( b\delta p) = \beta ( b\delta p) ,\; $$

where we assume the horizontal total stress change is linearly proportional to the pore pressure change with a numerical pre-factor 0 < γ < 1, and is negative as it is compressive. We then conduct a series of simulations to determine the numerical pre-factor 0 < β < 1.

We conduct a series of simulation under constant injection pressure or constant injection velocity boundary conditions by varying the modeling parameters, including b,  H _w,  V _inj,  L _crev and k ₀. We set H _w < H _i in all simulations with a constant pressure boundary condition. Figure 5 shows that when we combine all the numerical data of δσ′_xx,max and δp onto a single plot, our proposed scaling is robust across a range of parameters and both boundary conditions: δσ′_xx,max = β(bδp), where β = 0.22.

Figure 5. The scaling between δσ′_xx,max and bδp retrieved from simulations. Markers represent all simulation data in Figure 7 from water infiltration with a constant injection pressure (cross markers) or constant injection velocity (circular markers). We use the same colors and markers as Figure 7 for each simulation. The black, blue, red, green marker color represents a range of b, H _w or V _inj, L _crev, and k ₀, respectively. The dashed red line represents the analytical prediction: δσ′_xx,max = βδp (Eqn. (26)), with the prefactor β fitted to be 0.22.

3.2.2 The analytical expression of δp

To find the modeling parameters that impact the stresses, we develop analytical predictions of δp under the two boundary conditions. For the constant pressure injection at the crevasse tip, the pressure change is equal to the hydrostatic water pressure, δp = ρ _wgH _w. For the constant velocity injection at the crevasse tip, we derive δp from fluid continuity and Darcy's law. We assume that above a certain water infiltration depth H ₀, whose expression is derived later, water infiltrates much faster than the hydraulic conductivity of the firn. Therefore gravity is negligible and water invades in an approximately radially symmetric manner as shown in Figure 3(i)(ii). We approximate the flow near the tip as spreading out from the initial crack width L _crev at an angle θ. From the fluid continuity, we obtain:

(27)$$V_{\rm{inj}}L_{\rm{crev}} = V( r) ( L_{\rm{crev}} + \theta r) \rightarrow V( r) = {V_{\rm{inj}}L_{\rm{crev}}\over L_{\rm{crev}} + \theta r},\; $$

where V(r) is the radial velocity for water at distance r from the crevasse tip. At r = H ₀, the water velocity decays to the gravity-driven flow rate, which is also the hydraulic conductivity of water flow in porous firn. We derive the expression for H ₀ as follows:

(28)$$V( H_0) = V_{\rm{grav}} = {\rho_wgk_0\over \eta_w} \rightarrow H_0 = {L_{\rm{crev}}\over \theta}\left({\eta_wV_{\rm{inj}}\over \rho_wgk_0}-1\right).$$

Pore pressure diffuses from the crevasse tip (r = 0) to H ₀, below which the flow becomes gravity-driven, resulting in the observed elliptical-like shape of the water invading front. We calculate the pressure diffusion from r = 0 to H ₀ by Darcy's law:

(29)$$\matrix{V( r) & = -{k_0\over \eta_w}{\partial p\over \partial r}\rightarrow \displaystyle\int_{0}^{H_0}-{\eta_w\over k_0}V( r) \partial r = \displaystyle\int_{\delta p}^{0}\partial p,\; \cr & \rightarrow\delta p = {1\over \theta}{\eta_wV_{\rm{inj}}L_{\rm{crev}}\over k_0}ln\left({\eta_wV_{\rm{inj}}\over \rho_wgk_0}\right). }$$

For all constant injection velocity simulations (Figs. 7e–h), we plot the simulated δp as a function of ${\eta _wV_{\text {inj}}L_{\text {crev}}\over k_0}ln( {\eta _wV_{\text {inj}}\over \rho _wgk_0})$ as markers in Figure 6, the slope of which implies the prefactor θ in Eqn. (29). We conclude that simulation results agree well with the theory [Eqn. (29)] when θ = 3π/5. The inferred spreading angle θ is larger than π/2 because the infiltration front slightly deviates from a quarter circle under the effect of gravity. Note that Eqn. (29) applies to the condition when the water velocity decays to the gravity-driven flow rate before the invading front reaches approximately half of the domain depth. For an unrealistically large crevasse opening or water injection velocity at the crevasse tip, the invading front keeps expanding in a quarter circular-like shape, and H ₀ in Eqn. (29) is replaced by the depth when we terminate the simulation (${H\over 2}$ in this case). We include the analysis of large injection velocity or crevasse opening in the Supplementary Material. However, the large water pressure induced there (in the range of MPa) is not physical as it should have been capped by hydrostatic water pressure, ρ _wgH _w.

Figure 6. The scaling between the infiltration-induced pore pressure change at the crevasse tip (δp) and viscous pressure (${\eta _wV_{\text {inj}}L_{\text {crev}}\over k_0}ln( {\eta _wV_{\text {inj}}\over \rho _wgk_0})$) from the water infiltration with a constant injection velocity. Circular markers represent simulation data in Figure 7 from water infiltration with a constant injection velocity. We use the same colors as Figure 7 for each simulation. The blue, red, green marker color represents a range of V _inj, L _crev, and k ₀, respectively. The red dashed line represents the analytical prediction from Eqn. (29) with θ = 3π/5.

Figure 7. The dependence of δσ′_xx,max on modeling parameters (b,  H _w,  V _inj,  L _crev, k ₀) for the water infiltration with a constant injection pressure (top panels, (a)~(d)) and a constant injection velocity (bottom panels, (e)~(h)). On the top panel, H _w = 10 m except in (b), and on the bottom panel, V _inj = 0.05 m s⁻¹ except in (f). The markers represent simulation results, and the dashed red line represents analytical prediction from Eqn. (30) with β = 0.22.

3.2.3 The analytical expression of δσ′_xx,max

Now that we have the scaling between δσ′_xx,max and δp (Eqn. (26)), and the analytical expression of δp (Eqn. (29) with θ = 3π/5), we arrive at the final analytical expression of infiltration induced maximum effective stress as follows:

(30)$$\matrix{& \delta\sigma'_{xx, {\max}} = \beta( b\delta p) ,\; \cr & \delta p = \left\{\matrix{ \rho_wgH_w, \hbox{ with a constant }H_w,\; \cr {5\over 3\pi}{\eta_wV_{\rm{inj}}L_{{\rm crev}}\over k_0}ln\left({\eta_wV_{\rm{inj}}\over \rho_wgk_0}\right),\;\hbox{with a constant }V_{\rm{inj}}. }\right. }$$

To check that Eqn. (30) holds for the conducted simulations, Figure 7 presents the dependence of δσ′_xx,max on the modeling parameters, where the red dashed line represents the analytical prediction from Eqn. (30) with the numerical pre-factor β fitted to be 0.22.

Note that to cause fracture we need the maximum effective stress σ′_xx,max ≡ σ′_xx,0 + δσ′_xx,max to exceed the tensile strength σ′_t. Combining Eqn. (1), (6), and (30), we arrive at the final expression for the fracture criterion:

(31)$$\eqalign{\sigma'_{xx, \rm{max}} & = \sigma'_{xx, 0} + \delta\sigma'_{xx, \rm{max}} = -{\nu\over 1-\nu}\rho_i gH_i \cr& + {3K( 1-2\nu) \over ( 1 + \nu) ( 1-\nu) }\epsilon_{xx, 0}-{\nu\over 1-\nu}\rho_wg \langle H_w-H_i\rangle + \beta( b\delta p) \geq\sigma'_t,\; \cr \delta p & = \left\{\matrix{\rho_wgH_w ,\hbox{ with a constant} H_w, \cr {5\over 3\pi}{\eta_wV_{\rm{inj}}L_{\rm{crev}}\over k_0}ln\left({\eta_wV_{\rm{inj}}\over \rho_wgk_0}\right),\hbox{ with a constant $V_{\rm{inj}}$, } }\right.}$$

where 〈x〉 = max(x,  0) follows the rule of Macaulay bracket. The third term in the expression of σ′_xx,max accounts for additional background hydrostatic stress when there is ponding of the water on top of the ice slab. Figure 5 shows that the numerical pre-factor β is fitted to be 0.22.

3.3 Physical limits

Eqn. (31) provides different forms for the maximum effective stress in the firn depending on whether a constant pressure or constant injection velocity is assumed. Before applying this model to study the vulnerability of ice slab regions to hydrofracture, it is important to consider which boundary condition is most consistent with physical conditions on the ice sheet.

The constant pressure boundary condition straightforwardly represents a static water load in a partially or fully water-filled crevasse. It does not directly account for transient processes, such as water level fluctuations within a crevasses as water flows in from surface runoff or out through the permeable firn. However, by exploring the induced stresses for a range of water levels, we can constrain the plausible range of the maximum effective stress in the firn layer.

It is tempting to think of the constant injection velocity boundary condition as representing the transient case where water is flowing into the crevasse at a constant rate. However, this is not a good analogy. As Figure 4 shows, a constant injection rate leads to a roughly logarithmic-like increase in pressure with time, as more water is forced into the firn per time step than can be evacuated from the injection point due to the relatively low intrinsic permeability of the firn. However, a crevasse is not a closed system and the top remains open to the atmosphere. Therefore, when the rate of water injection into the crevasse exceeds the rate at which water can flow out through the firn, continuity of mass and pressure require that the water will start to fill the crevasse, rather than increase pressure in the firn layer. As a result, the constant injection rate boundary condition leads to artificially high estimates of firn pore pressure and, by extension, maximum effective stress, because our simulations assume a closed system and do not allow for backflow into the open crevasse. Therefore, we caution that the constant injection rate equations should not be used to calculate firn stresses. However, the constant injection rate model does provide an important relationship between pressure and injection velocity that we will later exploit to qualitatively estimate whether crevasses may fill with water, given typical stream flow rates into fractures.

4. Applications to the Greenland Ice Sheet

We now apply the analytical model developed in Section 3 to assess the vulnerability of Greenland's ice slab regions to hydrofracture. To do this, we seek to answer the following questions:

1. 1. Given typical firn conditions in Greenland, will fractures in ice slabs fill with water?

2. 2. If so, does hydrostatic loading of an ice slab crevasse induce sufficient stress in the firn layer to initiate fracture?

The first question is important because, as water flows into the top of an ice slab crevasse, either from distributed hillslope flow or where a stream intersects the fracture, it will also flow out of the fracture tip into the permeable firn. If water can be evacuated into the firn about as quickly as it enters the crevasse, the crevasse will not fill with water and there will be no additional hydrostatic stress that might drive hydrofracture. However, if the rate of infiltration into the firn is smaller than the rate of injection into the crevasse, the water level will rise within the crevasse. In this scenario, the second question becomes relevant, and we can apply Eqn 31 to determine whether, or under what conditions, the maximum effective stress induced in the firn layer would be sufficient to initiate fracture.

4.1 Mechanical and hydraulic parameters

To answer these questions by applying the analytical model developed in Section 3, we must first define reasonable values for the physical, mechanical, and hydraulic parameters of the ice slab-firn system in our area of interest. Unfortunately, given the large spatial extent of ice slab regions, the sparsity of subsurface observations within them, and the uncertainty in the few available measurements, it is difficult to choose a single representative value for any of these parameters. Therefore, we take a Monte Carlo simulation approach to this problem. For each variable, we define an empirical distribution of reasonable values using a compilation of in situ, laboratory, and remote sensing measurements reported in the literature. For the hydraulic and mechanical properties, which have never been measured directly in these regions, we use various empirical relations to define these properties as a function of firn density. Table 2 lists these unknown variables, the distributions we assign to them or the relation from which we calculate them, and the sources on which we base these distributions or relations.

Table 2. Monte carlo simulation parameters

* Note that ${\scr N}( \mu ,\; \, \sigma )$ denotes a normal distribution with mean μ and standard deviation σ and ${\scr U}[ a,\; \, b]$ denotes a uniform distribution over the values from a to b (inclusive).

The relation between open porosity and firn permeability is taken from Adolph and Albert (2014), which defined a power law relation between firn density and air permeability based on measurements from firn core samples collected at the North GRIP ice core drill site. We define our own relations between firn density and the mechanical parameters–Poisson ratio (ν) and Biot coefficient (b) – due to the lack of reports in the literature. The Poisson ratio has been measured with ultrasonic wave velocities at the laboratory scale and active seismic investigations at the field scale. We collate data sets from Schlegel and others (2019), King and Jarvis (2007), and Smith (1965) and use Monte Carlo simulation to build an expanded set of data points that cover the reported uncertainty for each measured data point. We then calculate the best linear fit between firn density (ρ _f) and ν using this expanded data set [Supplementary Figure 4] and define the uncertainty to be the one half of the 68% prediction interval on the measurements (reported as $\sigma _\nu$ in Table 2). The Biot coefficient is defined as a function of the ratio between the bulk modulus of the firn (K) and single grain elastic stiffness of ice (K _s). We define a relation between ρ _f and K in the same way as we did for ν, but this time using a quadratic fit to bulk modulus data compiled from the same sources [Supplementary Figure 5].

We run two sets of Monte Carlo simulations, the first where H _w is always less than or equal to H _i and the second where H _w is always greater than H _i. At each iteration in each simulation, we first draw ρ _i, ρ _f, L _crev, K _s, and H _i from the distributions defined in Table 2. If we are simulating a scenario in which we desire that H _w ≤ H _i, we then draw H _w from the distribution ${\scr U}[ 0,\; \, H_i]$. Alternately, if we desire that H _w > H _i, we draw H _w from the distribution ${\scr U}[ H_i + 0.1,\; \, H_i + 40]$. We use the randomly selected value of ρ _f to calculate K, ν, or k ₀ as appropriate. For example, K is drawn from a distribution defined as ${\scr N}( K_\mu ( \rho _f) ,\; \, \sigma _K)$ and then used to calculate b directly. For all the analyses that follow, we run 1,000,000 iterations of the Monte Carlo simulation, equivalent to solving the equations summarized in Figure 9 for a million different plausible configurations of an ice slab-firn system. The output of this simulation process is a distribution of the physically plausible range of solutions to the equation we are solving, given what we know about conditions in Greenland's ice slab regions.

4.2 Rate of crevasse filling

To determine whether fractures in ice slabs can fill with water, we start with the constant injection velocity solution in Eqn. (31) which provides a relationship between pressure and firn infiltration velocity. This equation shows that the change in pore pressure is related to the injection velocity, crevasse opening width, and firn permeability. Since the crevasse is open to the atmosphere at its top, we know that the maximum possible pressure change in the firn would be the hydrostatic pressure induced by a water-filled crevasse. Therefore, we set Eqn 31 equal to this hydrostatic pressure to estimate the maximum rate at which water can infiltrate into the firn.

(32)$$\rho_{\rm{w}}gH_{\rm{w}} = {5\over 3\pi}{\eta_wV_{\rm{inj}}L_{\rm{crev}}\over k_0}ln\left({\eta_wV_{\rm{inj}}\over \rho_wgk_0}\right).$$

Using the Monte Carlo approach described in Section 4.1, we numerically solve Eqn. (32) to compute a distribution of plausible infiltration rates (V _inj) into the firn beneath ice slabs.

We compare this distribution of infiltration rates to field measurements of supraglacial river and stream discharge to assess the balance between water flowing into and out of the crevasse. For this purpose, we reduce the system to two-dimensions and calculate discharge into the firn as follows:

(33)$$Q_{\rm{\,firn}} = V_{\rm{inj}}L_{\rm{crev}}.$$

We compare Q _firn (m² s⁻¹) to field observations of supraglacial stream and river discharge collected in the ablation zone of Southwest Greenland. Gleason and others (2016) measured width, depth, and discharge at 38 locations on a series of small streams – generally, less than 3 m wide and 0.3 m deep. We calculate an equivalent 2D discharge by dividing the measured volumetric discharge by the stream width. Smith and others (2015) also developed a power law relation between stream width (w) and discharge (Q) shown in Eqn. (34) based on field measurements in the Southwest Greenland ablation zone.

(34)$$Q = 0.099w^{1.85}.$$

To further estimate the range of plausible discharge rates for streams in ice slab regions, we measure the width of a series of representative streams in the Northwest and Southwest Greenland ice slab areas using WorldView imagery, making 3–5 measurements per stream (see Supplementary Fig. 6). Measured stream widths vary from 0.7-17 m, values which fall comfortably within the range of stream widths originally used to calibrate Eqn. (34) (Smith and others, 2015). We use Eqn. (34) to calculate the estimate discharge of these streams and rescale by the measured width to obtain equivalent 2D discharge. To ensure a conservative estimate of stream discharge, we limit this calculation to stream widths between 0.7-8 m.

Figure 8 shows the results of this analysis. The blue histogram shows the rate at which water would drain out of a crevasse in an ice slab into the underlying firn Q _firn, and the red histogram shows the rate at which water can drain into the crevasse from the surface based on the Gleason and others (2016) measurements. The gray shaded region shows the range of discharge rates calculated from our measured stream widths using Eqn. (34). We do not plot these values as a histogram, since we do not sample a sufficient number of streams to accurately capture the true distribution of stream widths in the region. However, we can bound the reasonable range of stream discharge from these selected observations. This analysis demonstrates two possible regimes. Q _surf and Q _firn are similar for the streams with the lowest discharge rates and crevasses with the largest openings, highest pressures, and most permeable underlying firn. This suggests that where crevasses are fed by small streams or hillslope flow, no hydrofracture will occur because water drains into the firn as quickly as it enters the crevasse, and therefore no hydrostatic stress is induced by water within the crevasse itself. However, we also see that the discharge from larger streams (>3 m wide) may be significantly greater than the rate at which water can drain out of a crevasse into the firn. In this second regime, the crevasse will fill with water, leading to an additional hydrostatic load that might be sufficient to vertically propagate the crevasse into the firn layer. We note that it is possible that we overestimate the typical discharge for large streams, since the Smith and others (2015) empirical relation was calibrated during a particularly high melt year and in the ablation zone where melt rates are higher. However, given that our firn infiltration rates are also overestimated, since they are based on the pressure of a fully water-filled crevasse, we assess that it is still plausible that some crevasses in ice slab regions may fill with water. Therefore, we address the second question – does hydrostatic loading of an ice slab crevasse induce sufficient stress in the firn layer to initiate fracture?

Figure 8. Comparison of the rate of water infiltration into the firn from the crevasse tip versus the rate at which surface streams may feed water into the crevasse. Blue bars show the distribution of firn water infiltration rates in response to a range of pressures equivalent to that induced by a water-filled crevasse. Red bars show small stream discharge values measured in the ablation zone of Southwest Greenland. The grey shaded region shows the range of plausible discharge values estimated from the empirical relation between stream width and discharge from Smith et al. (2015), given stream widths measured from WorldView imagery in the Northwest and Southwest Greenland ice slab regions. For larger stream, water infiltrates out of the crevasse tip more slowly than water enters the crevasse opening by stream flow, suggesting that some crevasses may partially fill with water.

4.3 Maximum effective stress in the firn layer

4.3.1 Effects of water depth

We first explore how the maximum effective stress in the firn layer changes as a function of water height in the crevasse. We consider two scenarios: a partially or fully water-filled crevasse (0 ≤ H _w ≤ H _i) and a supraglacial lake overtop a crevasse (H _w > H _i), where the surface area of the fracture is assumed to be negligible when compared to the total extent of the lake. Unfortunately, it is difficult to directly quantify if or when the magnitude of the maximum effective stress meets the fracture criterion defined in Eqn. (31) for two reasons. First, the tensile strength of firn (σ′_t) is not well constrained, particularly in Greenland where the mechanical properties may be significantly altered by water infiltration and refreezing. Second, we are not aware of a reliable way to estimate the background glaciological stress due to ice flow in a solid mechanics framework. This would require measurements of the background strain ($\epsilon _{xx, 0}$, rather than strain rate) at an appropriate time scale, which is not well-defined. Therefore, instead of directly comparing the total maximum effective stress in the firn to the tensile strength of firn, we instead compare the maximum effective stress in the firn to the total stress at the crevasse tip in an equivalent system composed of solid ice. In this way, we can evaluate whether the presence of firn will increase or decrease the likelihood of full-depth hydrofracture without assuming a particular value of tensile strength. To further simplify the analysis, we neglect the poorly-defined stress component due to the background glaciological strain ($\epsilon _{xx, 0}$), since this term is additive and should be of a similar magnitude for both the ice slab-firn case and solid ice case. Therefore, it should have minimal impact on the relative difference between the maximum effective stress in the firn and total stress in solid ice. We then define a non-dimensional maximum effective stress and non-dimensional water height as follows and rescale Eqn. (31).

(35)$$\tilde{\sigma}'_{xx, \rm{max}} \equiv {\sigma'_{xx, \rm{max}}\over \rho_{i}gH_{i}},\; $$

(36)$$\tilde{H} \equiv {H_w\over H_i}.$$

Figure 9 summarizes the final dimensional and non-dimensional equations for each scenario that we use in the remaining analysis, given the assumptions and simplification described above.

Figure 9. The theoretical prediction of the maximum effective stress at the crevasse tip σ′_xx,max. We present σ′_xx,max in both dimensional and dimensionless forms.

Given the equations in Figure 9, we use our Monte Carlo simulation approach to estimate the physically plausible distribution of maximum effective stresses that might be induced in the firn by a water-filled crevasse in an ice slab, given a range of water heights and the inherent uncertainty in the mechanical properties of the firn. Figure 10 shows the results of this non-dimensional analysis, with the average behavior of a ice slab-firn system shown in the white line, overlaid on a two-dimensional histogram that shows the spread of the possible solutions that would be consistent with the state of the firn in Greenland.

Figure 10. Non-dimensional analysis of maximum effective stress ($\tilde {\sigma }'_\equiv {\sigma '_{xx, {\max}}\over \rho _{i}gH_{i}}$) as a function of water height within the crevasse ($\tilde {H}\equiv {H_w\over H_i}$). The average behavior of an ice slab-firn system in Greenland is shown in the white line. This is overlain on a 2D histogram showing the full distribution of plausible solutions derived from the Monte Carlo analysis, given our uncertainty in the mechanical properties of the firn layer. The behavior for a solid ice system is shown in the blue line. When a firn layer is present, the maximum effective stress never becomes tensile, as most of the hydrostatic load is accommodated by a change in pore pressure, rather than being transmitted to the solid skeleton. Note that here we evaluate $\tilde {\sigma }'$ without the contribution of the background glaciological strain $\epsilon _{xx, 0}$. If $\epsilon _{xx, 0}$ is known it can be included via Eqn. (2).

We find that when a firn layer is present, the maximum effective stress is always compressive, even accounting for uncertainty in the mechanical properties, and in the case of an overlying supraglacial lake, actually becomes more compressive as the lake deepens. Typically, maximum effective stress in the ice slab-firn system is less than in a solid ice system, the exception being when $\tilde {H} = {H_w\over H_i}\lessapprox 0.6$. The non-dimensional form of the equations provides a clear explanation for the physics underlying this behavior. We can think of the first term in each equation – some constant multiplied by $\tilde {H}$ – as the hydrostatic term that describes how the maximum effective stress changes as the water pressure in the crevasse changes. The second term is a lithostatic term that describes the background state of stress in the system. Before water is added to the crevasse, the maximum effective stress in the firn is greater than in solid ice, because the firn's lower Poisson ratio reduces the proportion of the overburden stress is transmitted horizontally and can act to close the fracture. However, as water begins to fill the fracture, the maximum effective stress increases more slowly in the ice slab-firn system compared to the solid ice system, because only a portion of the hydrostatic stress is transferred to the solid skeleton, with the remainder being accommodated by an increase in pore pressure. As a result, once the water level in the crevasse exceeds $H_w\gtrapprox 0.6H_i$, the effective stress at the fracture tip in solid ice exceeds the maximum effective stress experienced by the firn. The exact point of this crossover can be calculated as a function of ν and b as shown in Eqn. (37).

(37)$$\tilde{H} = \left({1-2\nu\over 1-\nu}\right)\left({1\over 1-\beta b}\right)\left({\rho_i\over \rho_w}\right)\approx 0.6.$$

Once $\tilde {H}\geq {\rho _i\over \rho _w}$, the effective stress in the solid ice system will always be tensile. This is the critical conclusion of classical hydrofracture analyses in glaciology – that, due to the density difference between water and ice, a water-filled crevasse will always be in net tension and can propagate unstably (van der Veen, 2007; Lai and others, 2020). However, in the presence of an underlying firn layer, this transition is never reached and the maximum effective stress always remains compressive, due to the mitigating effect of pore pressure.

In the case of a supraglacial lake, where H _w > H _i, the classical solution for solid ice hydrofracture is identical to the solution for a completely water-filled crevasse and is independent of the depth of the lake. This is because the additional overburden from the lake contributes equally to the lithostatic and hydrostatic terms, so that the effective stress simply remains a function of the density difference between the ice and water in the crevasse. However, in the presence of an underlying firn layer, the maximum effective stress is significantly reduced by the presence of an overlying supraglacial lake and in fact decreases linearly as lake depth increases. The system behaves in this way because the overlying lake water contributes more to increasing the lithostatic stress acting to close the fracture than it does to increasing the hydrostatic stress, due to the competing influence of the Poisson ratio and Biot coefficient on the transfer of stress to the firn skeleton.

4.3.2 Effects of firn porosity

Since ν and b are both a function of firn porosity, we also explore the change in non-dimensional maximum effective stress as a function of firn porosity and non-dimensional water height. In Figure 11, we plot the same Monte Carlo data points shown in the histogram in Figure 10 in firn porosity versus $\tilde {H}\equiv {H_w\over H_i}$ space, taking the median simulation values in each 2D bin. For a partially or fully water-filled crevasse, the effective stress increases as the water height increases, due to the greater hydrostatic pressure. Effective stress also increases a firn porosity increases. Softer, more porous firn has a higher Biot coefficient and therefore a stronger fluid-solid coupling, so more of the hydrostatic stress is felt by the solid skeleton. More porous firn is also more compressible and has a lower Poisson ratio, so less of the lithostatic stress is transmitted horizontally and can act to close the crevasse. Instead, the firn compacts vertically under the overburden. In the case of a supraglacial lake over a crevasse, we instead find that an increase in lake depth reduces the effective stress due to the increasing influence of the lithostatic stress component. As expected, effective stress still increases as firn porosity increases, but this influence is more significant at greater lake depths, since it reflects how the coupling between hydrostatic stress and maximum effective stress is modulated by the Biot coefficient. Overall, we find that a low porosity, stiff firn matrix is the most stable.

Figure 11. Non-dimensional maximum effective stress as a function of firn porosity and non-dimensional water height in the crevasse. (a) Water-filled crevasses. The non-dimensional effective stress remains compressive but increases as firn porosity increases and water height increases due to the increasing water pressure, stronger fluid-solid coupling, and reduced lithostatic stress due to a lower Poisson ratio. (b) Supraglacial lake over a crevasse. Non-dimensional effective stress becomes more compressive as the water level increases, due to the additional lithostatic stress from the overlying lake. Firn porosity plays a great role in determining the non-dimensional effective stress as the water level increases, since it modulates both the hydrostatic stress transmitted to the solid skeleton, and the portion of the lithostatic stress transmitted horizontally.

5. Discussion

Our results demonstrate that the presence of a porous firn layer underneath Greenland's ice slabs leads to a significant resilience to hydrofracture in these regions. Where water drains into crevasses through hillslope flow or smalls streams, the rate of water injection into the crevasse is closely balanced by the rate at which water infiltrates into the permeable firn layer. As a result, water will not rise within the fracture and the hydrostatic pressure on the firn remains low and insufficient to cause hydrofracture in the firn layer. Where large streams, rivers, or lakes intersect fractures, the rate of surface water inflow may be significantly greater than the firn infiltration rate, leading to crevasse filling. However, even in this case, the firn largely stabilizes the system and prevents hydrofracture.

Specifically, we find that at water levels below the critical transition point defined in Eqn. (37), the underlying firn layer experiences a slightly greater maximum effective stress than if the system consisted of solid ice as shown in Figure 10. However, in both cases, the effective stress remains compressive. Once the water level exceeds this transition point, the firn becomes a stabilizing mechanism that significantly reduces the maximum effective stress as 78% of the total stress is accommodated by a change in pore pressure, rather than being transmitted to the solid skeleton. In the case of a supraglacial lake, a firn layer is particular stabilizing as the lithostatic stress increases more quickly with lake depth than the hydrostatic stress. Therefore, even in the case of a fully water-filled crevasse, hydrofracture of a porous firn layer appears to be highly unlikely. It is important to note that we do not include the effect of glaciological stresses due to ice flow in our analysis. Therefore, while it is clear that an ice slab-firn system is less likely to hydrofracture than a solid ice system under the same conditions, given sufficient backgrounds stress, hydrofracture in firn would certainly also be possible. A complete analysis of the conditions under which this might occur would require improved estimates of the tensile strength of icy firn, as well as a validated method for measuring ice strain over the relevant time scale within this solid mechanics approach.

5.1 Implications for Greenland

These results are consistent with the observations of ice blob formation in the Northwest Greenland ice slab region discussed in the introduction. Where a porous firn layer exists, the system is infiltration-dominated and leak-off of water from the fracture into the porous firn significantly reduces the likelihood of further crevasse propagation. Our analysis shows that this holds for all observed ice slab configurations on the Greenland Ice Sheet, not particular to the northwest. More generally, our work quantitatively validates the speculation in Culberg and others (2022) that a solid ice column is needed for hydrofracture and that therefore, as Greenland warms, there will be a time lag between the development of ice slabs and the formation of surface-to-bed connections that can couple the supraglacial and subglacial hydrology.

5.2 Implications for Antarctica

These results also have important implications for the future stability of Antarctica's ice shelves. Hydrofracture has been implicated in the breakup of the Larsen B and other ice shelves (Scambos and others, 2000, 2003; Banwell and others, 2013), leading to a loss of buttressing and significant accelerations in inland ice flow that increase mass loss from the continent (Scambos and others, 2004; Rignot and others, 2004). However, most ice shelves still retain some firn layer, and previous work hypothesized that all pore space in the firn would need to be filled with refrozen meltwater before hydrofracture could occur (Kuipers Munneke and others, 2014; Alley and others, 2018). This hypothesis was based on the assumption that surface ponds could not form until the firn layer was completely removed (Kuipers Munneke and others, 2014). The discovery of ice slabs in Greenland has since demonstrated that supraglacial hydrology may develop without complete filling of firn pore space (MacFerrin and others, 2019; Tedstone and Machguth, 2022), suggesting a potential mode for more rapid destabilization of ice shelves under ongoing atmospheric warming. Our results now quantitatively demonstrate that even if ice shelves were to develop ice slabs and rapidly transition from firn storage to supraglacial runoff in a similar way to Greenland, this alone would not be sufficient to prime them for immediate hydrofracture-induced disintegration. Instead, complete filling off all local firn pore space with solid ice may be necessary because of the firn's resilience to hydrofracture. This would require a longer period of sustained warming to achieve than the formation of ice slabs alone.

5.3 Assumptions and future work

While we have derived an idealized description of maximum effective tensile stress in the firn under a static water load, our results rest on a number of modeling assumptions that should be tested in future work. For example, since capillary pressure is insufficient to drive firn deformation or fracturing, we have focused on water infiltration rates that are much larger than the firn hydraulic conductivity, and thus can safely neglect capillarity in the model. Therefore our model cannot capture the gravity fingering instability under unsaturated flow conditions (Cueto-Felgueroso and Juanes, 2009). The effect of capillarity is beyond the scope of the current study, but might be important for studying the formation of ice pipes or ice lenses. With the large water infiltration rates in the model, it takes only several minutes to penetrate the depth of the firn layer. Therefore we neglect meltwater refreezing that takes hours (Moure and others, 2023) and snow compaction that takes years (Meyer and Hewitt, 2017).

The current poroelastic fracture criteria does not converge to the solid fracture criteria as the porosity or permeability is reduced. The difference comes from the boundary condition at the fracture tip. For solid ice, the medium is impermeable, and therefore the change of horizontal effective stress at the fracture tip is equal to the hydrostatic water pressure, δσ′_xx|_x=z=0 = ρ _wgH _w. For porous firn, we assume the firn layer underneath the crevasse tip is fully permeable (even as the porosity decreases) and leaves the skeleton stress-free at the crevasse tip, δσ′_xx|_x=z=0 = 0. Such assumption aligns with the fact that the effective stress always vanishes at the free surface of the solid skeleton in a porous medium (MacMinn and others, 2015; Auton and MacMinn, 2019; Meng and others, 2023). To span the transition from porous and permeable firn to non-porous and impermeable solid ice, we could introduce a “permeability load parameter” to modulate the boundary condition, as suggested in Auton and MacMinn (2019).

In terms of fracture dynamics, our model only predicts the conditions needed for fracture initiation in the firn layer and does not consider the dynamics of fracture propagation. Future work might consider the behavior of deeper crevasses that partially penetrate the firn layer, or the full transient propagation path of a shallower water-filled crevasse that initially is entirely within the ice slab. Similarly, here we have considered a static water load, but a fully transient model could be used to study the effect of diurnal fluctuations in water levels (or other transient filling processes) on the evolution of effective stress within the firn layer. As previously discussed, future work that incorporates background stresses due to the flow of ice, as well as the estimates of the tensile strength of firn, would also be of significant value.

6. Conclusions

Understanding the vulnerability of Greenland's ice slab regions to hydrofracture is critical for assessing where and how quickly the supraglacial and subglacial hydrologic systems can become coupled as the equilibrium line retreats inland. Previous observational work suggested that when meltwater drains into fractures in ice slabs, the water is largely be trapped in the underlying porous firn layer and does not reach the ice sheet bed (Culberg and others, 2022). However, this work presented little theoretical argument to support the idea that fractures could not propagate through a firn layer. Here, we developed a poromechanical model to analyze the maximum effective stress in the firn layer beneath a water-filled fracture in an ice slab. Our results show that the firn layer stabilizes the system in two ways. First, for low rates of water flow into a crevasse, this water can quickly leak-off into the firn and prevent the crevasse from filling in the first place. Second, even if water can fill the crevasse, a significant portion of the hydrostatic stress is accommodated by changes in pore pressure, reducing the effective stress felt by the solid skeleton and preventing fracture propagation in firn. Indeed, we find that the maximum effective stress in the firn remains compressive for all reasonable values of the firn material and hydraulic parameters and ice slab and crevasse geometries. This demonstrates that firn is broadly resilient to hydrofracture across a wide range of ice sheet conditions and confirms that surface-to-bed connections are unlikely to develop until all local pore space has been filled with refrozen ice, with the possible exception of regions with very high background stresses due to ice flow. Our model now provides a clear quantitative backing and physical explanation for the apparent lack of surface-to-bed connections in ice slab regions.