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, Reference Meier2007; Moore and others, Reference Moore, Grinsted, Zwinger and Jevrejeva2013). 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, Reference Colgan2016). Specifically, calving occurs when the combination of basal and surface crevasses penetrates the entire thickness, thus isolating an iceberg (e.g. Benn and others, Reference Benn, Warren and Mottram2007b). The mechanisms associated with crevasse initiation and propagation, however, are complex and involve mechanical, thermal and hydraulic fracture processes (e.g. Alley and others, Reference Alley, Dupont, Parizek and Anandakrishnan2005; Bassis and Ma, Reference Bassis and Ma2015). 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, Reference Pollard, DeConto and Alley2015; DeConto and Pollard, Reference DeConto and Pollard2016). Scambos and others (Reference Scambos, Hulbe, Fahnestock and Bohlander2000, Reference Scambos2009) and Banwell and others (Reference Banwell2014) 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, Reference Nye1957) or linear elastic fracture mechanics (LEFM) models (Weertman, Reference Weertman1971, Reference Weertman1973; Smith, Reference Smith1976; van der Veen, Reference van der Veen1998a, Reference van der Veen1998b). The zero stress model, first proposed by Nye (Reference Nye1957), 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 (Reference Jezek1984), and most notably Benn and others (Reference Benn, Hulton and Mottram2007a) and Nick and others (Reference Nick, Van der Veen, Vieli and Benn2010), 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, Reference De Robin1974; Weertman, Reference Weertman1974).
To address the limitations of the zero stress model, Weertman (Reference Weertman1971, Reference Weertman1973) 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 (Reference Weertman1973) 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 (Reference Weertman1977) 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 (Reference Smith1976) 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 (Reference van der Veen1998a, Reference van der Veen1998b) 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 (Reference Tada, Paris and Irwin1973, Reference Tada, Paris and Irwin2000).
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, Reference Jiménez and Duddu2018b), so their relevance to iceberg calving in real glaciers can be somewhat limited. Moreover, van der Veen (Reference van der Veen1999) 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 (Reference Mottram and Benn2009) tested crevasse depth models using a field study at Breid. h. oamerkurjökull glacier in Iceland, and found that the van der Veen (Reference van der Veen1998b) LEFM model performs better than the Nye (Reference Nye1957) zero stress model when tuned for fracture toughness along with input data on crevasse spacing. However, Mottram and Benn (Reference Mottram and Benn2009) 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, Reference Pralong and Funk2005; Pralong and others, Reference Pralong, Hutter and Funk2006; Duddu and Waisman, Reference Duddu and Waisman2012, Reference Duddu and Waisman2013; Duddu and others, Reference Duddu, Bassis and Waisman2013; Mobasher and others, Reference Mobasher, Duddu, Bassis and Waisman2016; Jiménez and others, Reference Jiménez, Duddu and Bassis2017), 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 (Reference Pralong and Funk2005) 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 (Reference Borstad2012) 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, Reference Keller and Hutter2014a, Reference Keller and Hutter2014b). 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 (Reference Jiménez, Duddu and Bassis2017), based on poromechanics theory (Mobasher and others, Reference Mobasher, Duddu, Bassis and Waisman2016). 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, Reference Weertman1973).
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, Reference Kachanov1958; Rabotnov, Reference Rabotnov1963) and the hypothesis of strain equivalence (Lemaitre, Reference Lemaitre1971), the effective Cauchy stress tensor ${{\bar {\bi \sigma}}}$
in ice may be defined in terms of the ‘true’ Cauchy stress tensor σ as
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, Reference Murakami1983; Duddu and Waisman, Reference Duddu and Waisman2013). 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 σ 1 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
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 σ 1 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 (Reference Mobasher, Duddu, Bassis and Waisman2016, Reference Mobasher, Berger-Vergiat and Waisman2017)
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, Reference Biot1955; Coussy, Reference Coussy1995), which formalizes Terzaghi's principle (Terzaghi, Reference Terzaghi1951, Reference Terzaghi2007). The term ‘poro-damage’ was first used by Bary and others (Reference Bary, Bournazel and Bourdarot2000), 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, Reference Christmann, Plate, Müller and Humbert2016), 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, Reference Glen1955; Cuffey and Paterson, Reference Cuffey and Paterson2010). 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
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
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
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, Reference Glen1955; Nye, Reference Nye1957)
For polycrystalline ice, the parameter n is generally taken to be 3. We obtained the viscosity coefficient A = 7.156 × 10−7 MPa−3 s−1 at a temperature of − 10°C by taking the parameter A 0 from van der Veen (Reference van der Veen2013) and then applying the Arrhenius equation (Eqn (2.14) in van der Veen (Reference van der Veen2013)). 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
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
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, Reference Karr and Choi1989; Rist and others, Reference Rist, Sammonds, Oerter and Doake2002). 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
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 (Reference Jiménez, Duddu and Bassis2017) 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, Reference Murakami1983) 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
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
${\bar \chi }$
is the effective Hayhurst stress invariant (Hayhurst, Reference Hayhurst1972; Murakami and others, Reference Murakami, Kawai and Rong1988) defined as
α 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, Reference Weiss2004). 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 (Reference Pralong and Funk2005) 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 (Reference Duddu and Waisman2012) 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 (Reference Duddu and Waisman2012), except for α, β and l c, which were assumed from Pralong and Funk (Reference Pralong and Funk2005)
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, Reference Jiménez, Duddu and Bassis2017)
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, Reference Bazant1994), 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, Reference Duddu and Waisman2013), 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, Reference Hillerborg, Modéer and Petersson1976)
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, Reference Borstad and McClung2013). 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, Reference Paterson1994; van der Veen, Reference van der Veen1998b). The cohesive strength σ c is typically much lower in value than the tensile yield strength of ice (Schulson and Duval, Reference Schulson and Duval2009). Previously, van der Veen (Reference van der Veen1998b) 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, Reference Pralong and Funk2005; Krug and others, Reference Krug, Weiss, Gagliardini and Durand2014) 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 (Reference Weiss2003). 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
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 (Reference Jiménez, Duddu and Bassis2017). 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
The above equations are subject to appropriate Dirichlet and Neumann boundary conditions,
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
where $D^{\rm max} = 0.999$
is the maximum value of damage, and ξ = 10−16 ≈ 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 (Reference Jiménez, Duddu and Bassis2017)), 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 1, and k 2, 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 1 and k 2 do affect the rate of crevasse penetration, but they do not significantly influence the final crevasse depth (except for k 2 = 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 (Reference Mobasher, Berger-Vergiat and Waisman2017). 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 1, and k 2 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 1, k 2, we simply use the values for polycrystalline ice calibrated in Pralong and Funk (Reference Pralong and Funk2005) and Duddu and Waisman (Reference Duddu and Waisman2012) 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, Reference Anderson2005, Section A2.3.2) and superposing with the linearly varying overburden pressure as
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, Reference Hutchinson1968; Rice and Rosengren, Reference Rice and Rosengren1968) singularity, as given by (see Anderson, Reference Anderson2005, Section 3.2.3)
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, Reference Rice1968).
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, Reference Scambos2009). 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, Reference Benn2017). 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, Reference Bassis and Walker2012). These diverse calving styles were previously found to be related to glacier geometry and boundary conditions (Bassis and Jacobs, Reference Bassis and Jacobs2013; Duddu and others, Reference Duddu, Bassis and Waisman2013), 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, Reference Greve and Blatter2009). 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 (Reference van der Veen1998a, Reference van der Veen1998b) and Nick and others (Reference Nick, Van der Veen, Vieli and Benn2010). 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, Reference Karr and Choi1989). 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 (Reference Mobasher, Duddu, Bassis and Waisman2016) 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, Reference Benn, Hulton and Mottram2007a, Reference Benn, Warren and Mottram2007b; Nick and others, Reference Nick, Van der Veen, Vieli and Benn2010; Cook and others, Reference Cook, Zwinger, Rutt, O'Neel and Murray2012, Reference Cook2014; Krug and others, Reference Krug, Weiss, Gagliardini and Durand2014; Todd and Christoffersen, Reference Todd and Christoffersen2014; Pollard and others, Reference Pollard, DeConto and Alley2015; Ma and others, Reference Ma, Tripathy and Bassis2017; Sun and others, Reference Sun, Cornford, Moore, Gladstone and Zhao2017; Yu and others, Reference Yu, Rignot, Morlighem and Seroussi2017; Todd and others, Reference Todd2018). 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, Reference Jiménez and Duddu2018b, 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, Reference Weertman1973). 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, Reference Jiménez and Duddu2018b). 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, Reference Lo, Borden, Ravi-Chandar and Landis2019), 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, Reference Wu and Nguyen2018). 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, Reference Murakami1983). 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 (Reference Duddu and Waisman2013). 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, Reference Schulson and Duval2009). 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, Reference Jiménez, Duddu and Bassis2017). 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, Reference Nick, Van der Veen, Vieli and Benn2010) and the stress-intensity-factor-based LEFM model (Krug and others, Reference Krug, Weiss, Gagliardini and Durand2014; Jiménez and Duddu, Reference Jiménez and Duddu2018b). 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 (Reference Weertman1973). 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, Reference Mottram and Benn2009; Enderlin and Bartholomaus, Reference Enderlin and Bartholomaus2019).
Acknowledgments
We gratefully acknowledge the funding support provided by the National Science Foundation's Office of Polar Programs via grants #PLR-1341428, #PLR-1341568 and #PLR-1847173. We also thank scientific editor Alan Rempel and the two anonymous reviewers for their comments that motivated the discussion section in this article.
Appendices
APPENDIX A. Long wavelength approximation
Crevasse models predict penetration depth as a function of the horizontal Cauchy stress σ xx. In the long wavelength ‘shallow’ approximation for an incompressible fluid, σ xx varies linearly with depth (Weertman, Reference Weertman1957) and can be decomposed into a (constant) deviatoric stress τ xx and a linearly varying lithostatic stress (e.g. Bassis and Walker, Reference Bassis and Walker2012),
where x, z are the in-plane horizontal and vertical coordinates (see Fig. 2), ρ i = 917 kg/m3 is the density of ice, g is gravitational acceleration, H is the height of the glacier, and the Macaulay brackets $\langle x \rangle = {1\over 2}\lpar x + \vert x\vert \rpar$
. A force balance at the terminus shows that in the long wavelength approximation, the horizontal deviatoric stress τ xx is given by
where ρ w = 1020 kg/m3 is the density of seawater and h w is the seawater depth at the glacier terminus. The long wavelength approximation relies on four idealizations that: (i) far from the terminus (i.e. in the so-called ‘far-field’ region of the glacier), the stresses in ice only vary with the z coordinate; (ii) ice deformation under isothermal conditions is incompressible; (iii) there is no resistance to flow at the basal boundary (i.e. basal friction is negligible); and (iv) the out-of-plane strain rate is negligible (i.e. flow is confined ${\dot \epsilon }_{yy}=0$
). Under the above conditions, the longitudinal strain rate ${\dot \epsilon }_{xx}$
is constant with depth, and the far-field stress σ xx is independent of ice rheology, as noted in Bassis and Walker (Reference Bassis and Walker2012). We verified that the constant τ xx approximation is valid when the above four conditions are satisfied in our finite element simulations, which allows us to compare CPDM model with analytical crevasse models.
APPENDIX B. Zero stress model
Nye (Reference Nye1957) first proposed the zero stress model by assuming that crevasses will penetrate through a glacier to the depth where tensile stresses vanish. In practice, two assumptions are commonly made with the application of this model: (i) that ice has zero cohesive strength; and (ii) that the stress singularity at crevasse tips can be neglected, owing to interaction of stress fields in a uniform field of closely-spaced crevasses. Although the zero stress model is rheology independent, ice is generally assumed to be an incompressible material, and Eqn (A1) is solved such that σ xx(z s) = 0, where z s is the vertical position of the crevasse tip measured from the bottom of the glacier. The maximum penetration depth of a surface crevasse d s = H − z s in the long wavelength approximation away from the calving front is given by
Under idealized conditions, for a glacier that terminates on land (i.e. the seawater level h w = 0), the deviatoric stress $\tau _{xx}={1\over 4} \rho _{{\rm i}} g H$
, and the predicted surface crevasse depth $d_s = {1\over 2}H$
.
The original Nye zero stress model did not account for the presence of hydraulic pressure in water-filled crevasses, however, later works (Jezek, Reference Jezek1984; Nick and others, Reference Nick, Van der Veen, Vieli and Benn2010; Bassis and Walker, Reference Bassis and Walker2012) included an additional term into the horizontal stress balance that incorporates hydraulic pressure,
where σ w(z) is the depth-varying hydraulic pressure acting on the crevasse walls. Because the hydraulic pressure acts to open the crevasse, it is assumed in the literature that this pressure induces a positive (tensile) stress within the ice, and crevasses are expected to penetrate to the depth where the net stress vanishes to zero (Weertman, Reference Weertman1980; Jezek, Reference Jezek1984). The pressure term σ w can be expressed as
where h s is the hydraulic head within the surface crevasse. By solving Eqn (B2) for σ xx(z s) = 0, we obtain the depth of a surface crevasse in the long wavelength limit:
where h s = c sd s denotes the level of water in the crevasse and c s ∈ [0, 1] is the fraction of the crevasse depth filled by water. For $c_s \ge {{\rho _{\rm i}} / {\rho _{{\rm w}}}}$
, the horizontal stress σ xx(z s) > 0 for any z s, thus indicating full depth penetration of the surface crevasse. For $c_s \lt {\rho _{{\rm i}}} / {\rho _{{\rm w}}}$
the problem is nonlinear, so we use an iterative algorithm based on the bisection method to solve for d s. The zero stress model results are show in Fig. 9.
Previously, Weertman (Reference Weertman1977) derived an approximate analytical solution to show that the crevasse penetration depth is equal to the Nye depth when the cohesive (fracture) strength of ice is essentially zero and the crevasse spacing is much smaller than crevasse depths. Thus, the Nye zero stress model can be categorized as a fracture-mechanics-based model, but its applicability is limited to a perfectly or nearly uniform crevasse field (De Robin, Reference De Robin1974; Weertman, Reference Weertman1974). Due to its simplicity, researchers employed the zero stress model (Nick and others, Reference Nick, Van der Veen, Vieli and Benn2010; Bassis and Walker, Reference Bassis and Walker2012; DeConto and Pollard, Reference DeConto and Pollard2016; Sun and others, Reference Sun, Cornford, Moore, Gladstone and Zhao2017), although it would underestimate the penetration depth of isolated or widely spaced surface crevasses. The restriction to the long wavelength approximation is overly limiting, but the zero stress model can also be directly applied without resorting to this approximation when the full Cauchy stress is available from full Stokes 2D and 3D simulations (Todd and Christoffersen, Reference Todd and Christoffersen2014; Todd and others, Reference Todd2018). Here we focus on only the long wavelength approximation because it provides a simple test case for comparing different crevasse models and is used in shallow ice dynamics models.
APPENDIX C. LEFM model
To address the limitations of the zero stress model, LEFM models were proposed by Smith (Reference Smith1976); van der Veen (Reference van der Veen1998a, Reference van der Veenb) and Krug and others (Reference Krug, Weiss, Gagliardini and Durand2014). The models are ideally suited for a rectangular-plate-like glacier made of linear elastic ice with a single surface crack subjected to the far-field stress state given by Eqns (A1) and (B2). In the LEFM models, crack growth is governed by the stress intensity factor at the crack tip, which is proportional to the applied far-field stress and the square root of the initial crack length. The maximum penetration depth of surface crevasses can be determined by equating the net stress intensity factor $K_{\rm I}^{{\rm net}}$
at the crevasse tip to the experimentally measured critical stress intensity factor $K_{\rm Ic}$
(Rist and others, Reference Rist1999). The net stress intensity factor is evaluated as
where z′ = H − z is the vertical coordinate measured from the top of the ice slab, and M(z′, H, d) is an appropriate weight function (Tada and others, Reference Tada, Paris and Irwin1973, Reference Tada, Paris and Irwin2000).
The LEFM models of van der Veen (Reference van der Veen1998a) and Krug and others (Reference Krug, Weiss, Gagliardini and Durand2014) used weight functions corresponding to a finite-width plate with a single edge crack that requires the boundary opposite of the crack (i.e. the basal boundary) to be a traction-free surface. Therefore, they are not appropriate for a grounded glacier with free tangential slip at the basal boundary (Jiménez and Duddu, Reference Jiménez and Duddu2018b). Instead, we use the following weight function for a finite-width plate with symmetric double edge cracks:
where the functions f 1 and f 2 are given as
The penetration depth of surface crevasses is determined by equating the net stress intensity factor $K_{\rm I}^{\rm net}$
to the experimentally determined critical stress intensity factor $K_{\rm Ic}=0.1\, {\rm MPa}\, \sqrt {{\rm m}}$
. Because the net stress intensity factor is a nonlinear function of the surface crevasse depth d s, we use an iterative algorithm based on the bisection method to solve for d s. Thus, our LEFM results shown in Fig. 6 are appropriate for idealized grounded glaciers with free basal slip.
APPENDIX D. Finite element implementation
We follow the standard Galerkin method of weighted residuals to derive the weak forms for the incompressible nonlinear viscous rheology and damage evolution equations, which involves multiplying by test functions w, q, and w, integrating by parts, and then applying the divergence theorem. For brevity, we only show the weak form for the nonlinear viscous rheology, which can be stated as
Find ${\bi v} \in {\cal V}$
, ${\bar p} \in {\cal P}$
, and ${\dot D} \in {\cal D}$
such that $\forall$
${\bi w} \in {\cal V}$
, $q \in {\cal P}$
, $e \in {\cal D}$
:
where ${\cal V}$
is a vector function space and ${\cal P}\comma\; {\cal D}$
are scalar function spaces with appropriate smoothness for the bulk field and their variations. The above weak form is implemented in the open source finite element software FEniCS (Alnæs and others, Reference Alnæs2015). The advantage of FEniCS is that the weak form corresponding to any coupled boundary value problem can be directly specified using the Python application program interface, and FEniCS interprets, assembles and solves the discretized linear system.
We employ a decoupled solution procedure that is explicit in time and consists of two sequential computations at each time-step. First, we compute the flow velocity and effective pressure using the variational form of the full Stokes equations (D1) using a mixed finite element discretization, and apply the Picard iteration scheme while holding damage constant at its previous time-step value. We then use the converged velocity and pressure solutions to determine the effective Cauchy stress and the local damage (production) rate at the current time-step. Second, we compute the damage rate using the variational form of the nonlocal gradient damage evolution Eqn (D2) using the standard finite element discretization, and apply the forward Euler scheme to update damage at the current time-step. Because the solution to the nonlinear Stokes equations is time-independent (steady-state) and damage evolution equation is time-dependent and highly nonlinear, it is simpler to use an explicit-time forward Euler scheme to update the damage separately. In contrast, the implicit-time backward Euler scheme will require a monolithic solve of the Stokes and damage evolution equations, so it can become very cumbersome.
For numerical accuracy, efficiency and stability, we ensure that the damage increment is sufficiently small (i.e. ΔD < 0.05) at any integration point at each time step and determine the time-step size as
Thus, the full-Stokes CPDM model is computationally intensive because we use small time-step sizes (≤2 hours) and mesh sizes (≤2.5 m) to alleviate numerical instability and mesh sensitivity. Additionally, we use the P2-P1 (Taylor-Hood) element for the Stokes equations, wherein the velocity v is resolved using quadratic (P2) interpolation on a six-noded triangle and the effective pressure ${\bar p}$
is resolved using linear (P1) interpolation on a three-noded triangle. However, we use three Gauss integration points for both P2 and P1 element, so that we can interpolate the variables at the same integration points. To determine the effective Cauchy stress, we evaluate the effective deviatoric stress at the three integration points using the nodal velocity field of the P2 element and add it to the corresponding effective pressure (i.e. at the same three integration points) computed using the nodal pressure field of the P1 element. The damage rate ${\dot D}$
is also resolved using the P1 element with three integration points because its evolution is dictated by the Cauchy stress. The nonlinear viscosity parameter $\eta \lpar {\dot {\bi \epsilon }}\lpar {\bi v}\rpar \rpar$
is a function of the velocity, so we apply a Picard (i.e. fixed-point) iteration scheme to implicitly solve the finite element system of equations. The full details of the numerical implementation are described in Jiménez and others (Reference Jiménez, Duddu and Bassis2017) and sample codes are available for download from the project website (Jimenez and Duddu, Reference Jimenez and Duddu2018a).
Open access
