A non-local continuum poro-damage mechanics model for hydrofracturing of surface crevasses in grounded glaciers

Abstract

Hydrofracturing can enhance the depth to which crevasses propagate and, in some cases, allow full depth crevasse penetration and iceberg detachment. However, many existing crevasse models either do not fully account for the stress field driving the hydrofracture process and/or treat glacier ice as elastic, neglecting the non-linear viscous rheology. Here, we present a non-local continuum poro-damage mechanics (CPDM) model for hydrofracturing and implement it within a full Stokes finite element formulation. We use the CPDM model to simulate the propagation of water-filled crevasses in idealized grounded glaciers, and compare crevasse depths predicted by this model with those from linear elastic fracture mechanics (LEFM) and zero stress models. We find that the CPDM model is in good agreement with the LEFM model for isolated crevasses and with the zero stress model for closely-spaced crevasses, until the glacier approaches buoyancy. When the glacier approaches buoyancy, we find that the CPDM model does not allow the propagation of water-filled crevasses due to the much smaller size of the tensile stress region concentrated near the crevasse tip. Our study suggests that the combination of non-linear viscous and damage processes in ice near the tip of a water-filled crevasse can alter calving outcomes.

1. Introduction

The mass loss from glaciers and ice sheets represents the largest contribution to eustatic sea level rise in the 21st century (Meier and others, 2007; Moore and others, 2013). Crevasses can influence both the mass balance and the dynamics of glaciers and ice sheets by enhancing surface ablation, meltwater retention, basal sliding, viscous deformation and iceberg calving (e.g. Colgan and others, 2016). Specifically, calving occurs when the combination of basal and surface crevasses penetrates the entire thickness, thus isolating an iceberg (e.g. Benn and others, 2007b). The mechanisms associated with crevasse initiation and propagation, however, are complex and involve mechanical, thermal and hydraulic fracture processes (e.g. Alley and others, 2005; Bassis and Ma, 2015). For instance, meltwater in surface crevasses or seawater in basal crevasses apply pressure on the crevasse walls and promote crevasse growth deeper into the glacier; this hydraulic-pressure-driven fracture is commonly referred to as hydrofracture. Thus, calving is intricately linked to climate dynamics through processes such as hydrofracturing, and can potentially contribute to rapid sea level rise (Pollard and others, 2015; DeConto and Pollard, 2016). Scambos and others (2000, 2009) and Banwell and others (2014) proposed that meltwater-driven hydrofracture of surface crevasses can lead to rapid fracturing of ice-shelf fronts and accelerate the collapse of certain portions of the Antarctic ice sheet. To gain a better understanding of the hydrofracturing process in glaciers and the conditions enabling full thickness crevasse propagation or calving, we formulate a physically-consistent, continuum poro-damage mechanics model and compare it with existing fracture mechanics models.

Crevasses have historically been thought to be opening (mode I) fractures formed under the action of tensile normal stress. Therefore, researchers sought to estimate crevasse penetration depths in glaciers and ice shelves analytically using either the zero stress model (Nye, 1957) or linear elastic fracture mechanics (LEFM) models (Weertman, 1971, 1973; Smith, 1976; van der Veen, 1998a, 1998b). The zero stress model, first proposed by Nye (1957), estimates the penetration depth of air-filled crevasses assuming that ice has zero fracture resistance or cohesive strength, and that crevasses will propagate to the depth where the longitudinal tensile stress in ice vanishes. Later, Jezek (1984), and most notably Benn and others (2007a) and Nick and others (2010), extended the zero stress model for estimating the depth of water-filled crevasses by assuming that they will propagate to the depth where the longitudinal tensile stress becomes equal to the (compressive) water pressure. The zero stress model neglects the stress singularity at the crevasse tip as predicted by the mathematical theory of elasticity; however, it may be valid for a quasi-uniform field of closely spaced crevasses in an idealized rectangular glacier (De Robin, 1974; Weertman, 1974).

To address the limitations of the zero stress model, Weertman (1971, 1973) proposed a dislocation-based LEFM model that accounts for the stress singularity at the crevasse tip and estimated the penetration depths of isolated water-free and water-filled crevasses such that the crack opening displacement is positive along the fracture surface. However, Weertman (1973) assumed that the ice thickness is infinite, and so this model is only applicable when the crevasse depth is much smaller than the ice thickness. Later, Weertman (1977) showed that the penetration depth of a uniform field of closely-spaced water-free crevasses, calculated using dislocation-based fracture mechanics, is equal to that found by the Nye zero stress model provided the fracture strength of ice is taken to be zero. Smith (1976) first proposed a stress-intensity-factor-based LEFM model for estimating the penetration depth of water-filled crevasses based on the assumption that the ice thickness is infinite. To account for the finite ice thickness of real glaciers, van der Veen (1998a, 1998b) proposed stress-intensity-factor-based LEFM models that used weight functions corresponding to a single edge crack in a finite width plate, as given by Tada and others (1973, 2000).

The aforementioned analytical crevasse depth models for mode I fracture, including zero stress and LEFM models, consider idealized or synthetic geometries and boundary conditions. Existing crevasse models are strictly applicable for mode I fracture of thin plates with rectangular geometry and with specific displacement or traction boundary conditions applied on their edges (Jiménez and Duddu, 2018b), so their relevance to iceberg calving in real glaciers can be somewhat limited. Moreover, van der Veen (1999) has argued that crevasses are the manifestation of mixed-mode fracture due to combined opening (mode I), sliding (mode II) and tearing (mode III) fracture. Mottram and Benn (2009) tested crevasse depth models using a field study at Breid. h. oamerkurjökull glacier in Iceland, and found that the van der Veen (1998b) LEFM model performs better than the Nye (1957) zero stress model when tuned for fracture toughness along with input data on crevasse spacing. However, Mottram and Benn (2009) also remarked that ‘measuring crevasse depths accurately is difficult, dangerous and time-consuming,’ and it is likely that the reported crevasse depths from field studies systematically underestimate the true depths, which makes both model calibration and validation challenging.

An alternative approach for calculating crevasse depth involves the use of continuum damage mechanics (CDM) within numerical ice flow models (Pralong and Funk, 2005; Pralong and others, 2006; Duddu and Waisman, 2012, 2013; Duddu and others, 2013; Mobasher and others, 2016; Jiménez and others, 2017), which offers certain advantages. First, a CDM model can be easily incorporated into numerical ice sheet models without requiring complicated algorithms for tracking the propagation of crevasses. Second, in a CDM model, damage is evolved based on the Cauchy (true) stress evaluated by the ice flow model and is valid for any arbitrary glacier geometry or boundary conditions. Third, a CDM model can account for crevasse initiation without requiring the presence of pre-existing crack or damage (i.e. subcritical damage), unlike the stress-intensity-factor-based LEFM model. Pralong and Funk (2005) were the first to implement a local creep damage model within an incompressible non-linear Stokes flow formulation to simulate crevasse propagation and glacier calving. Borstad and others (2012) used the CDM concept to estimate the damage in ice shelves based on satellite-observed surface velocities, but this study focused on assessment (diagnostic) and not on prediction (prognostic) of damage evolution.

Nonetheless, the CDM approach has its limitations. First, the CDM model within a full Stokes numerical formulation is computationally expensive for investigating crevasse propagation in real glaciers or ice shelves. Second, the CDM model involves several empirical parameters that may not be uniquely calibrated from existing experiments or observations, and this parametric uncertainty can affect its predictive capability. Third, the CDM model will need a tensorial damage variable to treat mixed-mode fracture under multi-axial tension-compression stress state, which can exacerbate the issues of computational cost and parametric uncertainty (Keller and Hutter, 2014a, 2014b). These limitations can be addressed by implementing shallow shelf/ice or higher-order Stokes approximations (instead of full Stokes formulations) using parallel computing, and calibrating model parameters with new experiments and satellite observations to reduce uncertainty. Thus, the CDM approach is a promising technique that can enable high-fidelity simulation and provide a better understanding of the mechanisms of thermal, mechanical and hydraulic fracture processes driving crevasse propagation. Furthermore, CDM can be advantageous when examining the propagation of fractures in different glacier geometries with spatially varying profiles of ice temperature, viscosity, density and boundary conditions.

Here we present a non-local continuum poro-damage mechanics (CPDM) model to simulate hydrofracturing by extending the CDM model of Jiménez and others (2017), based on poromechanics theory (Mobasher and others, 2016). The CPDM alleviates spurious mesh-size sensitivity and artificial diffusion of damage in crevasse propagation simulations and accounts for the feedback between viscous (or elastic) and damage processes at the crevasse tip. We compare the CPDM model with more traditional LEFM and Nye zero stress models to better understand how the predictions of these models vary under different idealized scenarios. We first illustrate the parametric and mesh-size insensitivity of the CPDM model for estimating the crevasse penetration depth. We next investigate the role of ice rheology, crack tip acuity and fracture/damage process zone size on the propagation of water-filled crevasses. Thus, through idealized simulation studies, we address two important modeling questions related to glacier calving: (1) Are the crevasse depths predicted by the CDPM model consistent with those predicted by zero stress and LEFM models? (2) What is the role of ice rheology and fracture process zone size on crevasse propagation and calving outcomes? In addition, we also address a fundamental glacier mechanics question: ‘Can a water-filled crevasse reach the bottom surface of a glacier?’ (Weertman, 1973).

The rest of the article is organized as follows: in Section 2, we present the strong form of the governing equations of the non-local CPDM model for hydromechanical fracture, including the notion of poro-damage mechanics, constitutive models for ice, and creep damage evolution law; in Section 3, we present parametric and sensitivity studies before comparing the depths of water-filled crevasses predicted by CPDM, LEFM and zero-stress models; in Section 4, we discuss the implications of the rheology and damage model assumptions and idealized simulation results for ice fracture and calving for real glaciers; and in Section 5, we offer some concluding remarks. In the Appendices, we briefly review the formulation of the analytical zero-stress and LEFM models, and the finite element implementation of the CPDM model.

2. Model formulation and implementation

2.1. Notion of poro-damage mechanics

We represent damage using an isotropic scalar variable D ∈ [0, 1], where D = 0 and D = 1 represent the undamaged and fully damaged state, respectively, at a material point in the continuum. The continuum region where 0 < D ≤ 1 describes a finite thickness damage zone, which may be interpreted as partially degraded material with distributed microcracks or microvoids before failure (D < 1) and an open crack after failure (D = 1). Thus, from a computational modeling perspective, CDM can be viewed as a diffuse interface alternative to LEFM that assumes a sharp (zero-thickness) crack interface. Typically, in the CDM the fully damaged zones (D = 1) are physically interpreted as air-filled or open cracks. Following the principle of effective stress (Kachanov, 1958; Rabotnov, 1963) and the hypothesis of strain equivalence (Lemaitre, 1971), the effective Cauchy stress tensor ${{\bar {\bi \sigma}}}$ in ice may be defined in terms of the ‘true’ Cauchy stress tensor σ as

(1)$${\bar{\bi \sigma}} = \displaystyle{{\bi \sigma } \over {\left( {1-D} \right)}} = \displaystyle{1 \over {\left( {1-D} \right)}}\left[ {\matrix{ {\sigma _{xx}} & {\sigma _{xy}} & {\sigma _{xz}} \cr {\sigma _{xy}} & {\sigma _{yy}} & {\sigma _{yz}} \cr {\sigma _{xz}} & {\sigma _{yz}} & {\sigma _{zz}} \cr } } \right].$$

We can extend the CDM formulation to incorporate hydraulic fracture in the quasi-static regime, where water can permeate the damage zone and exert hydrostatic pressure. In a finite-thickness damage zone (0 < D < 1) saturated with water within an otherwise undamaged ice slab, as shown in Figure 1, the isotropic damage variable D may be physically interpreted as the ratio of the area of microvoids and microcracks to the total area on a planar surface (e.g. a principal plane) through a representative volume element (RVE) (Murakami, 1983; Duddu and Waisman, 2013). Thus, the damage variable D is intricately related to porosity ϕ, which is defined as the ratio of the volume of microvoids to the total volume within the RVE. The physical RVE of the damage zone contains microcracks and microvoids within which there is hydraulic pressure, while the surrounding intact ice sustains an effective Cauchy stress. The equivalent RVE homogenizes the mechanical response such that the macroscopic maximum principal stress σ ₁ maintains force balance with the microscopic solid stress ${\bar {\sigma}}_{1}$ in ice and water pressure p _w within the voids or cracks. Thus, we can define a simple relation for the macroscopic maximum principal stress as

(2)$$\sigma_{1} = \lpar 1-D\rpar {\bar{\sigma}}_{1} - D p_{\rm w}\comma\; $$

where p _w is the hydraulic pressure acting within the damage zone.

Fig. 1. Schematic illustration of continuum poro-damage mechanics: (a) in a damage zone saturated with water, we assume that the microvoids and microcracks in the representative volume element (RVE) are completely filled with water; (b) on a principal plane of the physical RVE, hydraulic pressure p _w acts in the regions of the microcracks and microvoids, whereas the microscopic solid stress ${\bar {\sigma}}_{1}$ acts on the surrounding intact ice. Thus, the mechanical stress response is homogenized, so that the averaged macroscopic maximum principal stress σ ₁ on a principal plane in the equivalent RVE maintains force balance according to Eqn (2).

Extending the definition to three dimensions, the macroscopic Cauchy stress in saturated damaged ice is given by Mobasher and others (2016, 2017)

(3)$${\bi \sigma} = \lpar 1-D\rpar {\bar {\bi \sigma}} - D p_{\rm w} {\bi I}.$$

The above definition of the macroscopic Cauchy stress resembles that of the effective stress in a saturated porous medium according to Biot's theory of poroelasticity (Biot, 1955; Coussy, 1995), which formalizes Terzaghi's principle (Terzaghi, 1951, 2007). The term ‘poro-damage’ was first used by Bary and others (2000), although it is less commonly used in the literature.

2.2. Constitutive model for damaged ice

The mechanical response of ice over shorter time scales (seconds to hours) is often described by linear elastic or visco-elastic constitutive models (Christmann and others, 2016), whereas glacier and ice-sheet flow over longer times scales (days to centuries) is best described by a non-linearly viscous constitutive model known as Glen's law (Glen, 1955; Cuffey and Paterson, 2010). To better understand the role of ice rheology and fracture process zone on crevasse propagation, we perform two sets of numerical simulations, one with a linear elastic rheology and one with a non-linear viscous rheology, and compare the non-local CPDM model against LEFM and zero-stress models. Assuming the deformation or flow of glacier ice is incompressible, the effective stress in ice can be decomposed into deviatoric and volumetric parts as

(4)$${\bar {\bi \sigma}} = {\bar {\bi \tau}} - {\bar p} {\bi I}\comma\; $$

where ${{\bar {\bi \tau}}}$ is the effective deviatoric stress tensor, ${\bar p}=-{1\over 3}{\rm trace}\lsqb {{\bar {\bi \sigma} }}\rsqb$ is the effective pressure, and I is the identity tensor. Because pressure is constitutively indeterminate for incompressible deformation, we use a two-field velocity–pressure formulation with non-linear viscous rheology to define the deviatoric stress as a function of strain-rate. Similarly, we use a two-field displacement–pressure formulation with the linear elastic rheology to define the deviatoric stress as a function of strain.

2.2.1. Non-linear viscous rheology

Assuming polycrystalline ice to be an isotropic and incompressible non-Newtonian fluid, the effective deviatoric stress can be defined as

(5)$$\bar{{\bi \tau}} = 2 \eta \lpar {\dot {\bi \epsilon}}\rpar {\dot {\bi \epsilon}}\comma\; $$

where ${\dot {\bi \epsilon }} = {1\over 2} \lpar \nabla {\bi v} + \nabla ^\top {\bi v} \rpar$ is the viscous strain rate defined by the symmetric gradient of the velocity field v, and the non-linear viscosity coefficient is given by

(6)$$\eta ({\dot {\rm \epsilon }}{\rm )} = \displaystyle{1 \over 2}A^{-1/n}\left( {{{\dot {\rm \epsilon }}}^{\rm eq}} \right)^{1/n-1}.$$

In the above equation, A is a temperature-dependent viscosity coefficient, n is a viscosity exponent controlling the non-linearity of the flow model, ${\dot \epsilon }^{\rm eq} = \sqrt {{1\over 2} {\dot {\bi \epsilon }} : {\dot {\bi \epsilon }} }$ is the equivalent strain rate (i.e. the second strain-rate invariant) and the colon (:) denotes inner product between two tensors. The values of the model parameters used in our study are listed in Table 1. We note that this rheology model is consistent with the Glen–Nye law defined as (Glen, 1955; Nye, 1957)

(7)$${\dot {\bi \epsilon}} = A \left(\bar{\tau}^{{\rm eq}} \right)^{n - 1} \bar{{\bi \tau}}.$$

For polycrystalline ice, the parameter n is generally taken to be 3. We obtained the viscosity coefficient A = 7.156 × 10⁻⁷ MPa⁻³ s⁻¹ at a temperature of − 10°C by taking the parameter A ₀ from van der Veen (2013) and then applying the Arrhenius equation (Eqn (2.14) in van der Veen (2013)). Combining Eqns (4) and (5) and substituting them into Eqn (3), we obtain the constitutive equation for the non-linearly viscous ice rheology that incorporates hydraulic pressure within the damage zone as

(8)$${\bi \sigma} = \lpar 1-D\rpar \left[2\eta\lpar {\dot {\bi \epsilon}}\rpar {\dot {\bi \epsilon}} - \bar{\,p}{\bi I} \right]- D p_{\rm w} {\bi I}.$$

Table 1. Material properties of ice at − 10°C assumed in this study

2.2.2. Linear elastic rheology

Assuming polycrystalline ice to be an isotropic and incompressible elastic solid, the effective deviatoric stress can be defined as

(9)$$\bar{{\bi \tau}} = {E\over \lpar 1+\nu\rpar } {\bi \epsilon}\comma\; $$

where E is the Young's modulus and ν = 0.5 is the Poisson's ratio, and ${\bi \epsilon } = {1\over 2} \lpar \nabla {\bi u} + \nabla ^\top {\bi u} \rpar$ is the small strain tensor defined by the symmetric gradient of the displacement field u. For polycrystalline ice, we assume E = 9500 MPa at a temperature of − 10°C (Karr and Choi, 1989; Rist and others, 2002). Combining Eqns (4) and (9) and substituting them into Eqn (3), we obtain the constitutive equation for the linear elastic ice rheology that incorporates hydraulic pressure within the damage zone as

(10)$${\bi \sigma} = \lpar 1-D\rpar \left[{E\over \lpar 1+\nu\rpar } {\bi \epsilon} - \bar{\,p}{\bi I} \right]- D p_{\rm w} {\bi I}.$$

Equations (8) and (10) describe the rheology of intact ice for D = 0 whilst for D = 1 they describe the stress state in water under hydrostatic conditions. For 0 < D < 1 the stress is defined by a combination of the effective solid (ice) stress and fluid (water) stress, so that it satisfies microscale force balance. Thus, Eqns (8) and (10) can be interpreted as a stress decomposition that is valid for both linear poro-elastic and non-linear poro-viscoplastic constitutive models, respectively. In contrast, the superposition principle used for obtaining net stress by adding the glaciological stress ice with the water pressure in crevasses in the zero-stress or LEFM models (see Appendix B) is strictly valid only for linear elastic/viscous media.

2.3. Creep damage evolution law

We employ the gradient non-local creep damage mechanics formulation originally presented in Jiménez and others (2017) to simulate time-dependent crevasse propagation. Although the failure of ice is described by the progressive accumulation of microcracks and microvoids, the creep damage evolution law (Murakami, 1983) is phenomenologically formulated and does not explicitly identify the micromechanical mechanisms, such as void/crack growth or coalescence. Because damage is represented using a scalar variable, this damage evolution law cannot account for damage-induced anisotropy depending on microcrack orientation. We assume that isotropic damage will only increase under a tensile stress state (i.e. wherever the pressure is negative) to describe predominantly mode I brittle fracture. Therefore, we define the material time-derivative of local damage in the Lagrangian description as

(11)$${\dot D}^{{\rm loc}} = \left\{ {\matrix{ {B\displaystyle{{{\langle \bar{\chi }\rangle }^r} \over {{\left( {1-D} \right)}^{k_\sigma }}}} & {{\rm if}\;\bar{\,p} \le 0,} \cr 0\hfill & {{\rm if}\;\bar{\,p} \le 0,}}} \right. $$

where B is a damage rate coefficient, r is a damage rate exponent, $k_\sigma$ is an experimentally calibrated parameter accounting for the local damage rate enhancement due to prior damage and is defined as

(12)$$k_\sigma = k_1 + k_2 \, {\rm trace}\lsqb {\bi \sigma}\rsqb \comma\; $$

${\bar \chi }$ is the effective Hayhurst stress invariant (Hayhurst, 1972; Murakami and others, 1988) defined as

(13)$$\bar{\chi} = \alpha \bar{\sigma}^{\lpar {\rm I}\rpar } + \beta \bar{\sigma}^{\rm v} + \lpar 1 - \alpha - \beta\rpar \, {\rm trace}\lsqb \bar{{\bi \sigma}}\rsqb \comma\; $$

α and β are parameters describing the brittle and ductile nature of the fracture process, respectively, ${\bar {\sigma}}^{\lpar {\rm I}\rpar }$ is the effective maximum principal stress, and ${\bar {\sigma}}^{\rm v} = ( {3\over 2} {{\bar {\bi \tau}}} : {{\bar {\bi \tau}}} ) ^{1/2}$ is the effective von Mises stress.

In Eqn (11), we enforce that the local damage rate is zero if the (compressive) pressure p > 0 or the Hayhurst stress χ < 0. This condition still allows for tensile fracture in regions where macroscopic shear stress is the primary cause of damage, but arrests damage growth in regions that are confined or subject to multiaxial compression. Alternatively, we can enforce that the local damage rate is zero if the effective maximum principal stress ${\bar {\sigma}}^{\lpar {\rm I}\rpar }\lt 0$, but we note that this did not significantly affect the final crevasse depth in our simulation studies. At initial stages, when D ≪ 1 the term $\lpar 1-D \rpar ^{k_\sigma }$ is not dominant and the damage rate is determined by the Hayhurst stress χ, which describes subcritical crevasse nucleation and propagation (Weiss, 2004). We also define critical damage $D^{\rm cr}=0.6$ to capture the transition from slow damage growth in subcritical stages to fast damage growth and rupture at the later stages of creep fracture. In Eqn (13), α weights the effective Hayhurst stress toward the maximum principal stress, describing brittle failure behavior; whereas β weights the Hayhurst stress toward the von Mises stress, describing ductile failure behavior. The sum of these parameter values is constrained by α + β ≤ 1. Pralong and Funk (2005) previously estimated α = 0.21 and β = 0.63 based on limited laboratory experimental data, and so these values are not well calibrated. A more detailed discussion of the physical meaning and significance of all damage model parameters can be found in Duddu and Waisman (2012) and the parameter values used in this study are listed in Table 2.

Table 2. Damage law parameters are all assumed from Duddu and Waisman (2012), except for α, β and l _c, which were assumed from Pralong and Funk (2005)

To maintain thermodynamic consistency and alleviate mesh-size sensitivity, we implement a non-local implicit gradient formulation for the creep damage rate (Jiménez and others, 2017)

(14)$${\dot D}-\displaystyle{1 \over 2}l_{\rm c}^2 \nabla ^2{\dot D} = {\dot D}^{{\rm loc}}, $$

where l _c is an assumed non-local length scale. In a local damage model, the width of the damage zone (and the amount of strain energy dissipated) is dependent on the finite element mesh size (Bazant, 1994), which leads to pathological mesh-size dependence. By ‘smearing’ damage in a regularized manner within the damage zone, the non-local damage model alleviates both mesh-size dependence and directional mesh bias so long as the finite element mesh size is sufficiently smaller than length scale l _c (Duddu and Waisman, 2013), thus ensuring thermodynamic consistency. However, the length scale parameter is indicative of the length of the fracture process zone ahead of the crack tip and an approximate estimate is given by (e.g. Hillerborg and others, 1976)

(15)$$l_{\rm c}\approx \displaystyle{{EG_{{\rm Ic}}} \over {\sigma _{\rm c}^2 }} = \displaystyle{{K_{{\rm Ic}}^2 {\rm (}1-\nu ^2{\rm )}} \over {\sigma _{\rm c}^2 }}, $$

where $G_{\rm Ic}$ is the Griffith fracture energy. Typically, the quasi-brittle tensile fracture of earth materials such as snow, ice and rocks is characterized by two parameters, namely, the fracture toughness and the fracture process zone size (Borstad and McClung, 2013). However, only the fracture toughness of glacier ice is well characterized with the critical stress intensity factor in the range $K_{\rm Ic} = 0.1 - 0.4 \, {\rm MPa}\, \sqrt {{\rm m}}$ (Paterson, 1994; van der Veen, 1998b). The cohesive strength σ _c is typically much lower in value than the tensile yield strength of ice (Schulson and Duval, 2009). Previously, van der Veen (1998b) estimated that a stress of 30–80 kPa is needed for a single crevasse to form, which suggests the cohesive strength could be in that range. Other studies (Pralong and Funk, 2005; Krug and others, 2014) allude to the cohesive strength parameter based on a stress threshold for damage initiation and estimate it to be in the range of 10–200 kPa. Assuming the cohesive strength σ _c = 0.1 MPa for ice, we get the length scale in the range l _c ≈ 0.75 −12 m. Here we take l _c = 10 m based on the suggestion in Weiss (2003). Note that as l _c → 0, the non-local CPDM mimics the LEFM assumption of a sharp (or zero-thickness) crack, but its implementation becomes more computationally expensive because the finite element mesh size must be reduced correspondingly.

2.4. Strong form of governing equations

We specify the strong form of the governing equilibrium (momentum balance) and continuity (mass balance) equations for the incompressible flow or deformation of ice. For the incompressible non-linear viscous rheology, these equations are commonly known as the Stokes equations. Incorporating the constitutive Eqn (8) and the definition of the strain-rate, the strong form in terms of the unknown vector velocity field v and scalar effective pressure field ${\bar p}$ in the updated reference domain Ω is given by

(16)$$\left. \matrix{\hfill \nabla \cdot \left\{ {\left[ {1-D} \right]\eta \,\left( {{\dot {\bf \epsilon }}{\rm (}{\bi v}{\rm )}} \right)\,\left[ {\nabla {\bi v} + \nabla ^{\rm \top }{\bi v}} \right]} \right\} \cr \hfill -\left( {1-D} \right)\nabla \bar{\,p}-\nabla (Dp_{\rm w}{\rm )} + {\bi b} = 0 \cr \hfill \psi {\rm (}D{\rm )}\nabla \cdot {\bi v} = 0} \right\}\quad \quad {\rm on}\;\Omega $$

where b = {0, ρ _ig}^T is the external body force vector due to gravity and g is the acceleration due to gravity. If we take the water pressure p _w = 0, then the above poro-damage formulation becomes equivalent to the damage formulation of Jiménez and others (2017). Thus, we can parametrize hydrofracture by including one additional external force term $\nabla \lpar D p_{\rm w}\rpar$ in the CDM formulation, which is elegant and computationally attractive. For the incompressible linear elastic rheology, using the constitutive Eqn (10) and the definition of the strain, the strong form in terms of the unknown vector displacement field u and scalar effective pressure field ${\bar p}$ in the reference domain Ω is given by

(17)$$\left. \matrix{\hfill \nabla \cdot \left\{ {\left[ {1-D} \right]\displaystyle{E \over {(1 + \nu {\rm )}}}\left[ {\nabla {\bi u} + \nabla ^{\rm \top }{\bi u}} \right]} \right\} \cr \hfill -\left( {1-D} \right)\nabla \bar{\,p}-\nabla (Dp_{\rm w}{\rm )} + {\bi b} = 0 \cr \hfill \psi {\rm (}D)\nabla \cdot {\bi u} = 0} \right\}\quad \quad {\rm on}\;\Omega. $$

The above equations are subject to appropriate Dirichlet and Neumann boundary conditions,

(18)$$\eqalign{& v_i = {\bar{v}}_i\quad {\rm or}\quad u_i = {\bar{u}}_i\quad {\rm on}\;\Gamma ^{\rm D}, \cr & {\sigma {\hat {\bi n}}} = {\bar{\bi t}}\quad {\rm on}\;\Gamma^{\rm N},} $$

where ${\bar v}_i$ or ${\bar u}_i$ are prescribed velocities or displacements, respectively, on the Dirichlet boundary $\Gamma ^{\rm D}$, ${\bar {\bi t}}$ is the applied traction vector, and $\hat {{\bi n}}$ is the outward unit normal vector on the Neumann boundary $\Gamma ^{\rm N}$.

The function ψ(D) is incorporated in the formulation to relax the incompressibility constraint in fully damaged material points and is written as

(19)$$\psi (D) = \left\{ {\matrix{ {1,} \hfill & {{\rm if}\quad D < D^{{\rm max}},} \hfill \cr {\xi}, \hfill & {{\rm if}\quad D = D^{{\rm max}},} \hfill \cr } } \right.$$

where $D^{\rm max} = 0.999$ is the maximum value of damage, and ξ = 10⁻¹⁶ ≈ 0 is chosen as a small number close to machine precision. We restrict damage to a maximum value $D^{\rm max}=0.999$ to prevent ill-conditioning and rank deficiency in the tangent stiffness matrix, which allows the simulation of crevasse propagation to proceed until the final depth or full-depth penetration is achieved. At partially damaged material points (i.e. $0 \lt D \lt D^{\rm max}$), we assume the local density is equal to that of undamaged ice, but the viscosity is degraded. Once a material point fails (i.e. $D^{\rm max}=0.999$), we set the density to zero in the case of a dry crevasse, thus indicating open space, or to the density of water in the case that the material point is below the water level within a water-filled crevasse. Thus, we can describe a two-phase ice–water medium based on the damage variable and the crevasse is described by a finite-thickness zone of fully failed material. The weak form and the finite element implementation of the CPDM model is briefly discussed in Appendix D.

3. Numerical results

In this section, we simulate the propagation of water-filled surface crevasses in idealized grounded glaciers. We first illustrate the sensitivity of crevasse depth predictions to mesh size and model parameters. We then compare the penetration depths of isolated and closely spaced water-filled surface crevasses with those predicted by the LEFM and zero stress models, respectively. In all the simulations, we consider an idealized rectangular glacier of length L = 1000 m and height H = 125 m. For simplicity, we neglect lateral shear and restrict the domain to a flow line near the terminus of a grounded glacier with x and z representing the along-flow and vertical coordinates. The glacier is grounded on a rigid, frictionless (free-slip) bed and terminates at the ocean with a seawater depth h _w, as depicted in Figure 2. To discount the free translation motion of the glacier, we apply a boundary condition to enforce zero horizontal flow velocity at the left edge of the domain. The hydrostatic pressure from seawater is applied as a traction (Neumann) boundary condition normal to the right edge of the ice slab with magnitude −ρ _wg〈h _w − z〉. The hydraulic pressure p _w within the crevasse (damage zone) is applied using the poro-mechanics formulation with magnitude p _w = −ρ_wg〈h _s − (z − z _s)〉. To keep the ratio of h _s/d _s constant in the simulation, we detect the crevasse depth d _s at each time step and then adjust the water level h _s accordingly. To be consistent with the LEFM model, damage is only permitted to nucleate beneath the initial damage zone within a structured mesh with a constant element size $l_{\rm elem}$.

Fig. 2. Schematic of the glacier with height H, length L, seawater level h _w, surface crevasse height d _s and water level h _s within an isolated surface crevasse. The origin is set at the lower-left corner of the glacier with x and z as the horizontal and vertical coordinates, respectively, and y is the out-of-plane coordinate forming a right-handed system.

3.1. Isolated water-filled surface crevasse: sensitivity studies

In our previous work (see Section 4.1 in Jiménez and others (2017)), we verified the accuracy of the finite element formulation of the non-linear Stokes equations describing ice-sheet flow by showing that numerical results converge to a known analytical solution with progressive mesh refinement. In this section, we examine the sensitivity of crevasse propagation rates and the steady-state or final depth with meltwater in the crevasse and with seawater at the nearby terminus ($h_{\rm s}/d_{\rm s}=50\percnt$, $h_{\rm w}/H=50\percnt$). For comparison, we also examine an air-filled crevasse and without seawater pressure at the ice terminus (h _s = 0, h _w = 0).

3.1.1. Mesh size sensitivity

We first conducted a series of crevasse growth simulations using progressively refined finite element meshes. For each simulation, we use a structured mesh in the anticipated damage zone with element sizes $l_{\rm elem}=5$, 2.5 and 1.25 m. An initial 10 m × 10 m damage zone is centered along the top of the ice slab, and we only permit damage to nucleate beneath the initial damage zone to simulate the propagation of an isolated crevasse. This is done by setting the local damage rate ${\dot D}^{\rm loc}=0$ at integration points beyond the distance l _c from the vertical line centered at the initial notch (i.e. at x = L/2). The results of this study are shown in Figures 3a, c, which show the evolution of the crevasse depth ratio d _s/H with time t for an air-filled and water-filled surface crevasse, respectively. The crevasse depth d _s at a given time is obtained by determining the smallest vertical coordinate of fully-damaged material points and then subtracting it from the glacier height. This simple post-processing scheme for identifying crevasse depth within the damage mechanics approach obviates the need for implementing crack-tip tracking algorithms. For the air-filled crevasse, the CDM model predicts almost the same crevasse depth d _s versus time, provided that the finite element size $l_{\rm elem}$ is smaller than the characteristic length scale l _c. For the water-filled crevasse, the CPDM model predicts the crevasse depth d _s versus time with less than $8\percnt$ difference between the coarse mesh size ($l_{\rm elem}=5\, {\rm m}$) and the fine mesh size ($l_{\rm elem}=1.25\, {\rm m}$), and with less than $2\percnt$ difference between the medium mesh size ($l_{\rm elem}=2.5\, {\rm m}$) and the fine mesh size. We observed less than $1\percnt$ difference between the results from the fine ($l_{\rm elem}=1.25\, {\rm m}$) and finest mesh sizes ($l_{\rm elem}= 0.0625\, {\rm m}$).

Fig. 3. Surface crevasse depth d _s normalized with the domain height H = 125 m versus time for varying mesh size l _elem and varying non-local length scale size l _c. In subfigures (a) and (b), we consider a dry crevasse ($h_{\rm s}/d_{\rm s}=0\percnt$) within a land terminating glacier ($h_{\rm w}/H=0\percnt$); in subfigures (c) and (d), we consider a water-filled crevasse ($h_{\rm s}/d_{\rm s}=50\percnt$) within a marine terminating glacier ($h_{\rm w}/H=50\percnt$). The crevasse depth prediction from the non-local CDM approach using the poro-damage assumption is reasonably insensitive to finite element mesh size l _elem and the length scale parameter l _c, so long as the mesh size and length scale are sufficiently small.

3.1.2. Length scale sensitivity

We next conducted a series of crevasse growth simulations by progressively reducing the value of the non-local damage length scale l _c to determine if the uncertainty in l _c affects the predicted crevasse depth. For each simulation, we use a structured mesh in the anticipated damage zone with length scales l _c = 10, 5, 2.5 and 0.0625 m. An initial l _c × l _c m damage zone is centered along the top of the ice slab, and we only permit damage to nucleate beneath the initial damage zone, as in the previous study. We take the finite element size as $l_{\rm elem}=l_{\rm c}/4$ to ensure the numerical error is sufficiently small (${\lt } 2 \percnt$). Figures 3b, d show the length scale sensitivity of crevasse propagation with time (d _s/H versus t) for an air-filled and water-filled surface crevasse, respectively. Although the size of l _c affects the rate of crevasse propagation, it does not affect the final crevasse penetration depth so long as l _c is sufficiently small (i.e. l _c ≤ 10 m). Note that as the length scale l _c is reduced, the creep damage rate is smeared in a much thinner fracture process zone and the crevasse propagates more rapidly, thus describing brittle (instantaneous) crack propagation similar to LEFM in the limit when l _c → 0.

3.1.3. Parametric sensitivity

We next conducted a series of crevasse growth simulations by varying the damage model parameters $D^{\rm max}$, r, k ₁, and k ₂, to determine the parametric sensitivity of the predicted crevasse depth. The parameter $D^{\rm max}$ is used to specify a small residual viscosity and alleviate numerical convergence issues. The value of this parameter does not affect the final crevasse depth if it is chosen to be as close to unity as possible. Figure 4a illustrates that it is adequate to take $D^{\rm max}=0.999$ to obtain the final crevasse depth that is insensitive to this parameter. The sensitivity studies shown in Figures 4b–d indicate that values of r, k ₁ and k ₂ do affect the rate of crevasse penetration, but they do not significantly influence the final crevasse depth (except for k ₂ = 3). The damage parameter B also only affects the rate of crevasse propagation and not the final crevasse depth, as previously demonstrated in Mobasher and others (2017). We also studied the sensitivity to parameters α and β in a limited scope. We find that the predicted final crevasse depths showed negligible difference for α = 1, β = 0 (purely brittle failure) and α = 0.21, β = 0.64 (brittle-ductile failure). This indicates that for the ranges of values of α ∈ [0.21, 1] and β ∈ [0, 0.64], these parameters appear to have no influence on the steady-state crevasse depth. However, we emphasize that the material parameters α and β cannot be chosen arbitrarily and need to be calibrated using experimental data from fracture tests under multiaxial loading, which is currently unavailable.

Fig. 4. Parametric sensitivity study of damage law parameters $D^{\rm max}$, r, k ₁, and k ₂ shown in subfigures (a), (b), (c), and (d), respectively. We consider an isolated water-filled crevasse ($h_{\rm s}/d_{\rm s}=50\percnt$) in a marine terminating glacier (H = 125 m and $h_{\rm w}/H=50\percnt$). For each study, we plot the normalized surface crevasse depth d _s/H versus time. These simulation studies suggest that crevasse propagation rate is sensitive to damage model parameters, but the final crevasse depth is reasonably insensitive, thus reducing the uncertainty in crevasse depth predictions.

3.1.4. Summary

These studies indicate that the final crevasse penetration depth predicted by the non-local CPDM model is insensitive to the non-local length scale for l _c < 20 m, finite element mesh size for l _elem ≤ l _c/4 and damage control parameter for $D^{\rm max}\leq 0.999$. To be consistent with the LEFM and zero stress theories, henceforth we consider only purely brittle tensile fracture. Therefore, in all the following simulation studies we chose l _elem ≤ 2.5 m, l _c ≤ 10 m, α = 1 and β = 0 and $D^{\rm max}=0.999$. Because the final crevasse depth is reasonably insensitive to the choice of creep damage parameters B, r, k ₁, k ₂, we simply use the values for polycrystalline ice calibrated in Pralong and Funk (2005) and Duddu and Waisman (2012) using available laboratory test data.

3.2. Isolated water-filled surface crevasse: CPDM versus LEFM models

In this section, we compare the depths of isolated, water-filled surface crevasses predicted by the CPDM model with those from the LEFM model (see Appendix C). We consider the same geometry used previously, and we prescribe a 10 m × 10 m initial defect as a pre-damaged zone centered along the top surface of the slab. We first conducted crevasse growth simulations using the incompressible non-linear viscous ice and linear elastic rheologies. Figure 5 shows the evolution of the damage illustrating the propagation of an isolated water-filled surface crevasse in an rectangular ice slab of height H = 125 m for meltwater level ratio $h_{\rm s}/d_{\rm s}=50\percnt$ and seawater level ratio $h_{\rm w}/H=50\percnt$.

Fig. 5. The evolution of the damage profile showing the propagation of an isolated water-filled surface crevasse in an rectangular ice slab of height H = 125 m for meltwater level ratio $h_{\rm s}/d_{\rm s}=50\percnt$ and seawater level ratio $h_{\rm w}/H=50\percnt$. These finite element method (FEM) results are obtained using the non-linearly viscous (Stokes flow) rheological model and the continuum poro-damage mechanics (CPDM) approach with a purely brittle fracture criterion by setting α = 1 and β = 0 in Eqn (13).

3.2.1. Simulations with non-linear viscous ice rheology

Figure 6a shows the normalized crevasse depths d _s/H for isolated crevasses filled with water to varying levels, h _s/d _s = 0, 12.5, 25, 37.5, 50, 62.5, 75, 87.5 and 100%, within ice slabs terminating at the ocean with varying sea levels, h _w/H = 0, 50 and 90% (near-floating case). We assumed a purely brittle fracture criterion (i.e. α = 1 and β = 0) in the CPDM model and compared the results against the ‘double edge cracks’ LEFM model (see Appendix C). The CPDM model shows better agreement with the LEFM model when h _w/H = 0 and $50\percnt$ because the damage mechanics approach is able to account for stress concentration at the isolated crevasse tip. However, when $h_{\rm w}/H = 90\percnt$ (i.e. a near-floating glacier) and $h_{\rm s}/d_{\rm s} \gt 50\percnt$, the CPDM model deviates from the LEFM model. Our initial suspicion was that this discrepancy in the near-floating case arises due to difference in the assumption of ice rheology (i.e. linear elastic versus non-linear viscous) and this led to our next set of experiments.

Fig. 6. Surface crevasse depth d _s normalized with the domain height H = 125 m for varying water levels h _s filling the surface crevasse. The solid, dashed and dotted lines depict the ‘double edge cracks’ LEFM model result for different seawater depths h _w at the terminus. The markers (i.e. blue squares, orange triangles and black dots) represent finite element method (FEM) results using (a) non-linear viscous and (b) linear elastic rheological models. We employ the continuum poro-damage mechanics (CPDM) approach with the maximum principal stress (MPS)-based damage criterion by setting α = 1 and β = 0 in Eqn (13).

3.2.2. Simulations with linear elastic ice rheology

To test if the non-linear rheology of ice was responsible for the discrepancy between the CPDM and LEFM models, we used the incompressible linear elastic rheology to estimate the penetration depths of isolated water-filled crevasses, by assuming purely brittle failure (i.e. α = 1 and β = 0). Figure 6b shows the corresponding normalized crevasse depths d _s/H for isolated crevasses filled with water to varying levels h _s/d _s. We now observe crevasse growth in the near-floatation glacier case (i.e. $h_{\rm w}/H = 90 \percnt$), which shows a better agreement between the CPDM and the ‘double edge cracks’ LEFM model results. This indicates that the assumptions of the elastic ice rheology and purely brittle tensile fracture in the CPDM model are, not surprisingly, more consistent with the LEFM model. However, the CPDM-predicted crevasse depths still lag behind the LEFM depths. This is a consequence of the non-local damage length scale l _c = 10 m, which blunts the crevasse tip and decreases the crack driving stress, which is investigated next.

3.2.3. Comparison of stress at the crevasse tip in the near-floating glacier

To better understand the role of stress singularity and crack tip acuity, we conducted finite element simulations by representing the crevasse with a non-local damage zone and a zero-thickness (or sharp) crack by introducing double nodes along the crevasse walls. We refine the mesh size to a few centimeters near the crack tip to accurately capture the singularity in the stress field and use adaptive mesh refinement to coarsen the mesh away from the crack tip. We compared the maximum principal stress $\sigma ^{\lpar {\rm I}\rpar }$ in the vicinity of the crevasse tip (i.e. at the nearest integration point) by considering different ice rheology (linearly elastic and non-linearly viscous) for each case. In the sharp or zero-thickness crack case, hydraulic pressure on the crevasse walls is applied as a traction boundary condition; whereas, in the finite-width damage cases, hydraulic pressure is incorporated using the CPDM formulation. We evaluated the approximate stress field near the sharp crack tip in a linear elastic material by using Westergaard's stress function in polar coordinates with the crack tip as the origin (see Anderson, 2005, Section A2.3.2) and superposing with the linearly varying overburden pressure as

(20)$$\sigma ^{({\rm I})}(r) = \displaystyle{{K_I} \over {\sqrt {2\pi r} }}-\rho _{\rm i}g(d_s + r),$$

where r is the distance away from the crack tip along the crack length for a crevasse with depth d _s = 25 m. As shown in Figure 7, in the near-floating case ($h_{\rm w}/H = 90 \percnt \rpar$ the maximum principal stress at the tip of a fully water-filled crevasse ($d_{\rm s}/H = 100 \percnt$) is positive only within a small region near the crevasse tip. For linear elastic ice in Figure 7a, the stress becomes zero at about r ≈ 1.2 m for a sharp crack or damage zone with l _c = 1 m, and r ≈ 0.5 m for a damage zone with l _c = 10 m. For non-linear viscous ice in Figure 7b, the stress becomes zero at about r ≈ 0.7 m for a sharp crack, r ≈ 0.2 m for damage zone with l _c = 1 m and r ≈ 0.1 m for a damage zone with l _c = 10 m. We find that the stress decay with the non-linear viscous rheology (red line decays as ≈ r ^−1/4) is steeper than with the linear elastic rheology (black line decays ≃ r ^−1/2) in Figure 7b. This result matches well with the stress field near a sharp crack tip in a power-law strain hardening elasto-viscoplastic material defined by the so-called HRR (Hutchinson, 1968; Rice and Rosengren, 1968) singularity, as given by (see Anderson, 2005, Section 3.2.3)

(21)$$\sigma ^{({\rm I})}(r) = {k}^{\prime}\left( {\displaystyle{J \over r}} \right)^{{1// ( n + 1)}},$$

where k′ is a proportionality constant, r is the distance away from the crack tip along the crack length, the viscous exponent n = 3 and J is the path-independent integral (Rice, 1968).

Fig. 7. The largest value of the maximum principal stress $\sigma ^{\lpar {\rm I}\rpar }$ in the vicinity of the crack tip is plotted for a surface crevasse of depth d _s = 25 m. We take ice thickness H = 125 m, seawater level $h_{\rm w}/H=90\percnt$ (i.e. near-floatation glacier) and water level within the surface crevasse $h_{\rm s}/d_{\rm s}=100\percnt$ (i.e. the crevasse is fully filled). The crevasse is represented in three different ways: as a sharp (zero-thickness) crack and as damage zones with l _c = 1 m and l _c = 10 m. The stress $\sigma ^{\lpar {\rm I}\rpar }$ rapidly decays away from the crack tip in the non-linear viscous rheology case and becomes negative within a short distance, which explains the lack of crevasse propagation.

3.2.4. Summary

These studies indicate that the CPDM model is generally consistent with the ‘double edge crack’ LEFM model. The discrepancy between the CPDM and LEFM models in the near-floating case can be attributed to the differences in the assumption of ice rheology and crack tip acuity. Because the stress is positive (tensile) in a very small region, the CPDM model with non-linear viscous rheology in Figure 6b with l _elem = 2.5 m and l _c = 10 m does not predict crevasse propagation. In contrast, the CPDM model with linear elastic rheology in Figure 6a predicts crevasse propagation because the tensile stress region near the crevasse tip is twice as larger. Moreover, the CPDM does not exactly describe the sharp-crack LEFM model due to the effects of crack tip blunting introduced by finite size of the non-local damage zone. This study suggests that the combination of viscous and damage processes in ice near the tip of a water-filled crevasse in a near-floating glacier can alter crevasse penetration depths.

3.3. Closely-spaced water-filled surface crevasses: CPDM versus zero stress model

In this section, we compare the depths of closely-spaced, water-filled surface crevasses predicted by the CPDM model with those from the zero stress model (see Appendix B). We consider the same glacier geometry used previously, but prescribe several 10 m × 20 m initial defects spaced 50 m apart as pre-damaged zones along the top surface of the slab. To decrease the computational cost, we exploit the periodicity of crevasses and reduce the length of the computational domain to L _r = 200 m. In all simulations, gravity loading is applied as a body force, and a Dirichlet condition on the flow velocity ${\bar v}_x= {\dot \epsilon }_{xx}L_r$ is prescribed on the right-edge, whilst ${\bar v}_x=0$ is prescribed on the left-edge of the reduced domain. The horizontal strain rate ${\dot \epsilon }_{xx}$ is calculated using Eqn (7) by substituting the appropriate τ _xx calculated from Eqn (A2) with the glacier height H = 125 m. With this setup, we conducted several crevasse growth simulations and determined the final penetration depths for both incompressible non-linear viscous and linear elastic ice rheologies. Figure 8 shows the evolution of the damage illustrating the propagation of closely-spaced water-filled surface crevasse in an rectangular ice slab of height H = 125 m for meltwater level ratio $h_{\rm s}/d_{\rm s}=50\percnt$ and seawater level ratio $h_{\rm w}/H=50\percnt$.

Fig. 8. The evolution of the damage profile showing the propagation of closely-spaced, water-filled surface crevasses in an rectangular ice slab of height H = 125 m for meltwater level ratio $h_{\rm s}/d_{\rm s}=50\percnt$ and seawater level ratio $h_{\rm w}/H=50\percnt$. These finite element method (FEM) results are obtained using the non-linearly viscous (Stokes flow) rheological model and the continuum poro-damage mechanics (CPDM) approach with a purely brittle fracture criterion by setting α = 1 and β = 0 in Eqn (13). To decrease the computational cost, we exploit the periodicity of crevasses and reduce the length of the computational domain to L _r = 200 m (dark shaded region). We then reconstruct the full glacier domain by stitching together several periodic domains (light shaded regions).

3.3.1. Simulations with elastic and viscous rheology

Similar to the studies in the previous section, we estimate the penetration depths of surface crevasses in relation to sea water height and water level in the crevasse. In Figure 9, we show the normalized crevasse depths d _s/H of a uniform field of crevasses filled with water to varying levels, h _s/d _s = 0, 12.5, 25, 37.5, 50, 62.5, 75, 87.5 and 100%, within ice slabs terminating at the ocean with varying sea levels, h _w/H = 0, 50 and 90% (near-floating case). The CPDM model shows better agreement with the zero stress model for h _w/H = 0 and $50\percnt$, but deviates for $h_{\rm w}/H = 90\percnt$ (i.e. a near-floating glacier) and $h_{\rm s}/d_{\rm s} \gt 50\percnt$. Unlike in the isolated crevasse studies with the linear elastic ice rheology, the CPDM model does not predict any crevasse propagation in the near-floating glacier case for closely-spaced crevasses. Moreover, the stress concentration at the crevasse tip is diminished by the interaction with neighboring crevasses, so the final crevasse depths are smaller for closely-spaced crevasses than those obtained for the isolated crevasses. Comparing Figures 9a, b, we notice that the crevasse depths predicted by the CPDM model using linear elastic and non-linear viscous ice rheologies are generally the same, except for a few cases with $d_{\rm s}/H \gt 50 \percnt$ and $h_{\rm s}/d_{\rm s} \gt 75 \percnt$. The anomalies (or non-monotonic behavior) in predicted crevasse depths in Figure 9a arise from the prescribed velocity boundary condition ${\bar v}_x= {\dot \epsilon }_{xx}L_r$ on the right edge of the reduced domain of length L _r = 200 m. As crevasses propagate deeper, the ice slab stretches or flows faster, so calculating ${\bar v}_x$ with strain rate ${\dot \epsilon }_{xx}$ from Eqn 7 is not consistent anymore. Specifically, in non-linear viscous rheology simulations, this boundary condition induces additional compressive pressure for crevasse depth ratios $d_{\rm s}/H\gt 50\percnt$, which leads to crack arrest.

Fig. 9. Surface crevasse depth d _s normalized with the domain height H = 125 m for varying water levels h _s filling the closely-spaced field of crevasses. The solid, dashed and dotted lines depict the zero stress model results for different seawater depths h _w at the terminus. The markers (i.e. blue squares, orange triangles and black dots) represent finite element method (FEM) results using (a) non-linear viscous and (b) linear elastic rheological models. We employ the continuum poro-damage mechanics (CPDM) approach with the maximum principal stress (MPS)-based damage criterion by setting α = 1 and β = 0 in Eqn (13).

3.3.2. Summary

The zero stress model is independent of ice rheology parameters, if we assume incompressibility of ice flow or deformation. The CPDM model is generally consistent with zero stress model in that the depths of closely-spaced crevasses are insensitive to ice rheology and is much less than that predicted by the LEFM model for similar boundary conditions. This is because the singularity at the crevasse tips vanishes due to stress redistribution and the stress at the crevasse tip becomes equal to that given by the long wavelength approximation. The prescribed velocity boundary condition on the right edge of the reduced domain is necessary to facilitate the growth of a uniform field of closely-spaced crevasses; however, as crevasses propagate deeper we find that this boundary condition leads to an inconsistency with the non-linear viscous rheology, so the predicted crevasse depths do not agree well with the Nye zero stress model. We expect that the implementation of periodic boundary conditions will mitigate this issue, but this would only be an appropriate assumption away from the glacier calving front. Near the calving front, a traction boundary is more appropriate and closely-spaced crevasses will not propagate uniformly because the stress distribution is non-uniform in the longitudinal (x) direction.

4. Discussion

4.1. Vulnerability of near-floating grounded glaciers to hydrofracture

The CPDM model with non-linear viscous rheology predicts that a fully water-filled crevasse may not penetrate the entire thickness of a near-floating grounded glacier, which contradicts the prediction of the LEFM and zero stress models. We illustrated that this difference mainly arises due to stress redistribution in non-linearly viscous ice drastically that reduces the size of the tensile region at the crevasse tip. This leads us to the conclusion that non-linear viscous and non-local damage processes can reduce the vulnerability of near-floating grounded glaciers to hydrofracture; this has implications for Greenland glaciers with abundant meltwater. It is possible that the CPDM model implemented within a large-scale ice-sheet simulation would result in a more stable calving front and predict smaller calving rates from grounded tidewater glaciers, when compared to LEFM and zero stress models. However, our simulations neglect the additional effect of bending once calving fronts become buoyant and may not be representative of floating Antarctic ice shelves (see, e.g. Scambos and others, 2009). Moreover, a combination of surface and basal crevasse propagation can cause calving in floating ice shelves, which is influenced by variations in ice-front buoyancy (Benn and others, 2017). Observations of iceberg calving indicate diverse range of styles, for example, floating ice shelves calve episodic tabular icebergs and grounded tidewater glaciers frequently calve smaller icebergs (Bassis and Walker, 2012). These diverse calving styles were previously found to be related to glacier geometry and boundary conditions (Bassis and Jacobs, 2013; Duddu and others, 2013), but this study implicates that ice rheology and crack-tip acuity can also play an important role.

4.2. Compressible elastic and incompressible viscoplastic rheology of ice

Generally, glacier or ice-sheet flow models assume incompressible non-linear viscoplastic rheology, which allows for efficient numerical simulation over longer geophysical timescales using Stokes-based approximations (Greve and Blatter, 2009). As given by Eqns (A1) and (A2), the long wavelength stress approximation assumes incompressibility, which is used in the analytical LEFM and zero stress models of van der Veen (1998a, 1998b) and Nick and others (2010). For the sake of comparison, in this paper we also considered the incompressible linear elastic rheology so that the far-field stress state in the glacier (sufficiently away from the terminus) is the same across all models. However, on shorter timescales ice exhibits compressible viscoelastic and incompressible viscoplastic behavior with the Young's modulus E = 9500 MPa and the Poisson's ratio ν = 0.35, as calibrated from laboratory experiments (Karr and Choi, 1989). For compressible elasticity, the bulk modulus is approximately three times the shear modulus (i.e. same order of magnitude), but for (nearly) incompressible elasticity the bulk modulus is several orders of magnitude larger than the shear modulus. It is important to note that the longitudinal stress variation with depth predicted by the compressible linear elastic ice rheology is different from that given by Eqns (A1) and (A2). We previously considered this compressible elastic and incompressible viscoplastic rheology along with local creep damage evolution law in Mobasher and others (2016) and find that surface crevasses can propagate to a greater depth with this rheology compared to purely incompressible viscous/elastic rheology used in this study. Therefore, it may be important to include and understand the influence of compressible elastic nature of ice on crevasse and rift propagation in glaciers and ice shelves over short and long time scales.

4.3. Hydrofracture parameterization in ice-sheet models

In recent years, several studies implemented quasi-analytic calving criteria based on zero stress or LEFM models in Stokes-flow-based formulations (Benn and others, 2007a, 2007b; Nick and others, 2010; Cook and others, 2012, 2014; Krug and others, 2014; Todd and Christoffersen, 2014; Pollard and others, 2015; Ma and others, 2017; Sun and others, 2017; Yu and others, 2017; Todd and others, 2018). These calving criteria incorporated a simple hydrofracture parametrization by defining a net longitudinal stress as the sum of the longitudinal stress generated by the extensional flow of the glacier (glaciological stress) and the hydrostatic pressure of water filling the crevasse. It is important to note that the net stress is not exactly the same as the Cauchy (true) stress at the crevasse tip, but only an equivalent applied stress for estimating the stress intensity factor at the crevasse tip in a linear elastic ice slab using the weight function method (see Jiménez and Duddu, 2018b, for more details). Although the net stress definition is not strictly applicable with viscous ice flow models due to non-linearity, this simple parametrization still accounts for the fact that a water-filled crevasse propagates deeper than a water-free crevasse (Weertman, 1973). For a specific ice-shelf or glacier this parameterization can be well-calibrated, but it is not clear whether it can accurately represent hydrofracture in a continental-scale simulation. As further investigation is warranted, we urge caution when using quasi-analytic calving laws relying on LEFM and zero stress theories in viscous Stokes-flow-based ice-sheet/glacier models. In contrast, the CPDM model uses the decomposition of stress in a saturated damaged zone into effective stress in ice and water pressure within crevasses. The effective stress is then computed from mass and momentum balance equations along with the ice constitutive law. Thus, CPDM provides an alternative approach for representing hydrofracturing that is physically-consistent and simpler to implement in ice-sheet models.

4.4. Realistic simulation of glacier calving process

The CPDM model accounts for the gradual degradation of ice viscosity due to damage accumulation and represents open crevasses as fully damaged zones, so it captures the influence of subcritical damage and crevasse propagation on glacier and ice-sheet dynamics. The CPDM model can handle arbitrary glacier geometries and frictional basal boundary conditions, whereas the weight functions in LEFM are specific to rectangular geometries and idealized boundary conditions, namely, free slip or no slip at the glacier base (Jiménez and Duddu, 2018b). In this paper, however, we considered only an idealized case – a rectangular grounded glacier with free basal slip and initialized it with a pre-crack or notch – to make a direct comparison between the CPDM model and the quasi-analytic LEFM and zero stress models. Note that the ‘double edged crack’ weight function used in the LEFM model is specific to an idealized rectangular plate with free tangential slip and zero normal displacement at the basal surface (see Appendix C). The CPDM approach proposed here can be used with any ice rheology or damage evolution law. For example, if we use the linear elastic ice rheology and a strain-energy-based phase field damage evolution law (e.g. Lo and others, 2019), then CPDM can emulate LEFM without having to determine stress intensity factors or appropriate weight functions. However, the drawback is that it is computationally more expensive due to mesh size and length scale restrictions, which may be overcome through length scale insensitive formulations (e.g. Wu and Nguyen, 2018). Thus, the CPDM approach is versatile and can enable realistic simulation of glacier calving process, even when the damage zone ahead of crevasses tips is not vanishingly small and the viscous flow of the ice cannot be neglected. However, to fully assess the role of viscous and non-local damage processes on fracture propagation, more detailed experiments and observations are required to better constrain the damage evolution law and other physical aspects, such as ice compressibility or fabric evolution.

4.5. Mixed-mode fracture propagation under multi-axial stress

Usually, brittle and quasi-brittle materials like glacier ice, concrete and rocks are weaker in tension and are susceptible to mode I fracture. The creep damage model accounts for brittle-ductile nature of failure in polycrystalline materials in the form of microvoids and microcracks in a phenomenological way, but it is mainly suited for mode I (opening) fracture under tension–torsion loadings (Murakami, 1983). Under macroscopic shear, torsion or biaxial tensile loading, brittle fracture in ice can still be determined by maximum tensile (principal) stress and the creep damage model based on scalar stress invariants can capture the propagation of inclined or curved fractures, as demonstrated in Duddu and Waisman (2013). Therefore, we assume that isotropic damage will evolve only under a tensile stress state, which is also consistent with the assumptions of LEFM and zero stress models for mode I hydrofracture. However, mode II (sliding) fracture is plausible in ice under compressive stress states, leading to the propagation of wing or comb cracks near existing inclined cracks (Schulson and Duval, 2009). We note that the creep damage model does not describe mode II fracture due to sliding and interfacial shear stress at the crack tip. Although LEFM or cohesive zone models can resolve mode I and mode II fracture, a key challenge is that mode II fracture toughness of brittle glacier ice is difficult to determine directly from laboratory experiments, and data are scarce in the existing literature. Further work is necessary for developing a mixed-mode fracture criteria through integrated laboratory experiments and numerical simulations.

5. Conclusion

In this paper, we presented a non-local CPDM model for hydrofracturing and used it to simulate the propagation of water-filled surface crevasses. When crevasses are air-filled, this model reduces to the standard non-local CDM model (Jiménez and others, 2017). We first illustrated that the CPDM model-predicted final crevasse depths are independent of the mesh size and the damage length scale, and insensitive to creep damage parameters. We then compared the CPDM-predicted penetration depths of surface crevasses with two existing models in the literature: the zero-stress model (Nick and others, 2010) and the stress-intensity-factor-based LEFM model (Krug and others, 2014; Jiménez and Duddu, 2018b). The CPDM model simulates the time-dependent propagation of water-filled crevasses; whereas, the LEFM or Nye zero stress models predict if a crevasse with a specified depth will propagate or not. Therefore, we compared the final (or steady-state) crevasse depth obtained from the CPDM model (after 2 years of real-time, when crevasses cease to propagate in all simulated cases) with the depths obtained from LEFM or Nye zero stress models. The CPDM model illustrates that the presence of water pressure within crevasses results in deeper penetration into the glacier, which is in agreement with Weertman (1973). The CPDM model results show good agreement with those of the LEFM model for isolated water-filled crevasses and the zero stress model for closely-spaced water-filled crevasses, so long as the glacier does not approach buoyancy. However, as the glacier approaches buoyancy (i.e. near-floating case $h_{\rm w}/H=90\percnt$), the CPDM model does not allow the propagation of water-filled crevasses, contradicting the zero stress and LEFM models that allow full depth propagation.

In conclusion, the non-local CPDM approach combining poro-mechanics and damage-mechanics theories is suited for simulating hydrofracturing. The broader applicability of the CPDM model can be realized through its ability to describe crevasse formation and propagation in relation to complex glacier/ice-shelf geometries and boundary conditions, thus improving our understanding of the calving process. This is because the CPDM model consistently accounts for the effects of water pressure within crevasses (i.e. damage zone) on the local stress field and ice rheology. The CPDM model within a full Stokes flow formulation is only viable for simulating fracture propagation in small, idealized cases and not for real glacier simulations, owing to the restrictions on time-step size and mesh resolution (see discussion in Appendix D). However, using shallow shelf or depth-integrated viscosity approximations and higher order Runge–Kutta time-update schemes for damage can reduce the computational cost and enable realistic ice-shelf fracture simulation. The CPDM model can also be useful in developing simpler calving criteria that can be incorporated into large-scale ice-sheet models, as an alternative to LEFM or Nye zero stress models. Finally, we acknowledge the need to validate crevasse models, including the proposed CPDM model, using crevasse depth and spacing measurements, although there is some uncertainty associated with current observational data on crevasse depths (Mottram and Benn, 2009; Enderlin and Bartholomaus, 2019).