2.1. Viscoelastic rheology of undamaged ice

Assuming small elastic deformations, we additively decompose the total strain tensor ε into elastic and viscous components

(1) $${{\epsilon} _{{\rm kl}}} = {\epsilon} _{{\rm kl}}^{\rm e} + {\epsilon} _{{\rm kl}}^{\rm v}, $$

where ${\epsilon} _{{\rm kl}}^{\rm e} $ is the elastic strain (time independent and recoverable) component and ${\epsilon} _{{\rm kl}}^{\rm v} $ is the viscous (time dependent and irrecoverable) component. Making the usual assumption that glacier ice is isotropic, owing to its random polycrystalline microstructure, the elastic stress/strain relationship (multi-axial Hooke's law) is given by

(2) $${\epsilon} _{{\rm kl}}^{\rm e} = \displaystyle{1 \over E}[{\sigma _{{\rm kl}}} - \nu ({\sigma _{{\rm ii}}}{\delta _{{\rm kl}}} - {\sigma _{{\rm kl}}})],$$

where σ _kl denote components of the Cauchy stress tensor, E is Young's modulus, ν is Poisson's ratio, δ _kl is the Kronecker delta and repeated indices imply summation. The above equation can be rewritten in the form

(3) $${\sigma _{{\rm ij}}} = {C_{{\rm ijkl}}}{\epsilon} _{{\rm kl}}^{\rm e}, $$

where the fourth-order elasticity tensor C _ijkl is defined as

(4) $$\eqalign{{C_{{\rm ijkl}}} & = \displaystyle{E \over {2(1 + \nu )}}({\delta _{{\rm il}}}{\delta _{{\rm jk}}} + {\delta _{{\rm ik}}}{\delta _{{\rm jl}}}) \cr & \quad + \displaystyle{{E\nu} \over {(1 + \nu )(1 - 2\nu )}}{\delta _{{\rm ij}}}{\delta _{{\rm kl}}}.} $$

Denoting the deviatoric stress, $\sigma _{{\rm kl}}^{{\rm dev}} = {\sigma _{{\rm kl}}} - {\sigma _{{\rm ii}}}{\delta _{{\rm kl}}}/3$ , the viscous rheology can be expressed using the power-law creep relationship

(5) $${\dot \epsilon} _{{\rm kl}}^{\rm v} = A\tau _{\rm e}^{n - 1} \sigma _{{\rm kl}}^{{\rm dev}}, $$

where A is the temperature dependent viscosity coefficient, n is the flow law exponent, the dot decoration over a variable represents the material derivative and $\tau _{\rm e}^{\rm 2} = (3/2) \sigma _{{\rm ij}}^{{\rm dev}} \sigma _{{\rm ij}}^{{\rm dev}} $ is the (von Mises) stress invariant. Because in a Maxwell viscoelastic model the stress felt by the viscous and elastic elements is the same, provided that the elasticity tensor does not vanish, the stress can be computed from Eqn (3) and then used directly in Eqn (5) to compute the viscous strain rate.

2.2. Viscoelastic rheology of damaged ice

We assume isotropic damage of ice under tension to simplify the formulation. We introduce a scalar internal state variable D, such that its evolution from D = 0 to 1 represents the deterioration of ice from the fully intact undamaged state to the completely damaged state. For 0 < D < 1, we can define an effective stress ${\bar \sigma _{{\rm ij}}}$ in the material that is larger than the average stress σ _ij defined by

(6) $${\bar \sigma _{{\rm ij}}} = \displaystyle{{{\sigma _{{\rm ij}}}} \over {1 - D}}.$$

Following the hypothesis of equivalent strain (Lemaitre, Reference Lemaitre1992), the stress/strain relationship for damaged ice, in terms of the effective stress, can be expressed as

(7) $${\bar \sigma _{{\rm ij}}} = {C_{{\rm ijkl}}}{\epsilon} _{{\rm kl}}^{\rm e}. $$

The damage modified viscous strain rate tensor is now defined in terms of the effective stress as

(8) $${\dot \epsilon} _{{\rm kl}}^{\rm v} = A\bar \tau _{\rm e}^{(n - 1)} \bar \sigma _{{\rm kl}}^{{\rm dev}}, $$

with the effective deviatoric stress $\bar \sigma _{{\rm kl}}^{{\rm dev}} = {\bar \sigma _{{\rm kl}}} - {\bar \sigma _{{\rm ii}}}{\delta _{{\rm kl}}}/3$ and $\bar \tau _{\rm e}^{\rm 2} = (3/2)\bar \sigma _{{\rm ij}}^{{\rm dev}} \bar \sigma _{{\rm ij}}^{{\rm dev}} $ . Substituting these definitions into Eqn (8), the damage magnified viscous strain rate is

(9) $${\dot \epsilon} _{{\rm kl}}^{\rm v} = \displaystyle{A \over {{{(1 - D)}^n}}}\tau _{\rm e}^{(n - 1)} \sigma _{{\rm kl}}^{{\rm dev}}. $$

Thus, the presence of damage leads to a nonlinear strain rate enhancement factor of (1 − D)⁻ⁿ .

2.3. Effect of pore pressure on the rheology of damaged ice

The previous sections described a constitutive creep damage model for polycrystalline ice, analogous to those developed by Pralong and Funk (Reference Pralong and Funk2005), Duddu and Waisman (Reference Duddu and Waisman2012) and with the exception of the damage production law Krug and others (Reference Krug, Weiss, Gagliardini and Durand2014). In this section, we extend the continuum damage model to incorporate hydraulic fracture under wet (saturated) conditions where water can penetrate into microfractures (Fig. 1). Recalling Biot's theory of poroelasticity (Biot, Reference Biot1941, Reference Biot1955), the relationship between the homogenized Cauchy stress σ _ij and macroscopic solid effective stress ${\tilde \sigma _{{\rm ij}}}$ under saturated conditions can be defined as

(10) $${\sigma _{{\rm ij}}} = (1 - \phi ){\tilde \sigma _{{\rm ij}}} - \phi {P^{\rm h}}{\delta _{{\rm ij}}},$$

where P ^h represents the pressure of water filling the pores of a permeable medium, and ϕ is the average porosity of the medium. Assuming that damage and porosity are equivalent in isotropically damaged ice as a first approximation, we can extend the definition of the effective stress ${\bar \sigma _{{\rm ij}}}$ (defined in Eqns (6) and (7)) to express the homogenized Cauchy stress as

(11) $${\sigma _{{\rm ij}}} = (1 - D){\bar \sigma _{{\rm ij}}} - D{P^{\rm h}}{\delta _{{\rm ij}}}.$$

Fig. 1. Schematic showing the main features of the model. Surface and basal crevasses are present in the grounded ice.

Note that in the dry case P ^h = 0 and Eqn (11) reduces to the original definition of effective stress defined in Eqn (6). For simplicity, we further assume that water flow into pores is sufficiently rapid so that the water pressure P ^h in the microvoids and microcracks in damaged ice is hydrostatic. This implies

(12) $${P^{\rm h}} = {\rho _{\rm f}}g\langle z\rangle, $$

where ρ _f is the fluid density, g is the gravitational acceleration, 〈〉 denote Macaulay brackets defined such that 〈χ〉 = χ when χ > 0 and 〈χ〉 = 0 when χ < 0, and the hydraulic head z is the vertical distance between the water surface level z ₀ and the level of the material point z ₁ (i.e. z = z ₀ − z ₁).

Combining Eqns (12) and (11), we can now write the macroscopic stress/strain relationship as

(13) $${\sigma _{{\rm ij}}} = (1 - D){C_{{\rm ijkl}}}{\epsilon} _{{\rm kl}}^{\rm e} - {\rho _{\rm f}}g\langle z\rangle D{\delta _{{\rm ij}}}.$$

When D = 0 the above equation for Cauchy stress reduces to that of the undamaged material and when 〈z〉 = 0, that is, when the hydraulic head is below the material point, the equation reduces to that of the damaged material under dry conditions. Note that the effective solid stress ${\bar \sigma _{{\rm ij}}}$ increases under saturated conditions, as given by

(14) $${\bar \sigma _{{\rm ij}}} = {C_{{\rm ijkl}}}{\epsilon} _{{\rm kl}}^{\rm e} = \displaystyle{{{\sigma _{{\rm ij}}}} \over {1 - D}} + \displaystyle{D \over {1 - D}}{P^{\rm h}}{\delta _{{\rm ij}}}.$$

The damage modified viscous strain rate tensor under wet conditions is given by

(15) $${\dot \epsilon} _{{\rm kl}}^{\rm v} = A\bar \tau _{\rm e}^{(n - 1)} \bar \sigma _{{\rm kl}}^{{\rm dev}}, $$

with the effective deviatoric stress $\bar \sigma _{{\rm kl}}^{{\rm dev}} = {\bar \sigma _{{\rm kl}}} - {\bar \sigma _{{\rm ii}}}{\delta _{{\rm kl}}}/3$ and $\bar \tau _{\rm e}^{\rm 2} = (3/2)\bar \sigma _{{\rm ij}}^{{\rm dev}} \bar \sigma _{{\rm ij}}^{{\rm dev}} $ . Recalling that neither the von Mises stress, ${\bar \tau _{\rm e}}$ , nor the deviatoric stress, $\sigma _{{\rm kl}}^{{\rm dev}} $ , depend on pore pressure, Eqn (15) shows that unlike the elastic rheology, the viscous component of the rheology is invariant to the inclusion of pore pressure. This is illustrated more explicitly for the uniaxial example in Appendices A and B.

2.4. Damage evolution law

To complete the constitutive damage model description, we need to specify the damage evolution law. There is large uncertainty in the appropriate specification of a damage evolution law, but our formulation in theory is more general and does not depend on the particular choice of evolution law. Nonetheless, for definiteness we adopt a power-law isotropic damage rate function analogous to that was previously introduced by Pralong and Funk (Reference Pralong and Funk2005) and Duddu and Waisman (Reference Duddu and Waisman2012, Reference Duddu and Waisman2013):

(16) $$\dot D = \displaystyle{{B{{\langle \chi \rangle} ^r}} \over {{{(1 - D)}^{{k_\sigma}}}}}, $$

where the damage rate coefficient B is a (possibly temperature dependent) model parameter, the exponent r is a chosen constant, the exponent k _σ is a stress dependent parameter and χ is the Hayhurst equivalent stress given by Pralong and Funk (Reference Pralong and Funk2005):

(17) $$\chi = \left( {\matrix{ {\alpha {{\bar \sigma} ^{(1)}} + \beta {{\bar \tau} _{\rm e}} - (1 - \alpha - \beta ){{\bar \sigma} _{{\rm kk}}},} \hfill & {{\rm if}\quad {{\bar \sigma} ^{(1)}} \gt 0,} \hfill \cr {0,} \hfill & {{\rm if}\quad {{\bar \sigma} ^{(1)}} \le 0.} \hfill \cr}} \right.$$

In the above equation, α and β are material parameters that control the damage growth mechanism (a detailed discussion on their selection is provided by Pralong and Funk (Reference Pralong and Funk2005)), ${\bar \sigma ^{(1)}}$ is the largest effective principal stress and ${\bar \tau _{\rm e}}$ is the effective von Mises stress corresponding to the solid matrix of porous ice. In this paper, we consider that failure occurs due to tensile stresses only so that ice remains intact in compression. Hence, damage only accumulates when ${\bar \sigma ^{(1)}} \gt 0$ . Creep experiments on polycrystalline materials (including metals and ice at high temperatures) illustrate that the rate of creep damage growth increases drastically as we approach full collapse. To be consistent with previous studies (Pralong and Funk, Reference Pralong and Funk2005; Duddu and Waisman, Reference Duddu and Waisman2012) we introduce a stress dependent exponent of the form

(18) $${k_\sigma} = \left( {\matrix{ {{k_1} + {k_2}{{\bar \sigma} _{{\rm ii}}},} \hfill & {{\rm for}\quad 0 \le {{\bar \sigma} _{{\rm ii}}} \le 1\;{\rm MPa},} \hfill \cr {{k_1} + {k_2},} \hfill & {{\rm for}\quad {{\bar \sigma} _{{\rm ii}}} \gt 1\;{\rm MPa},} \hfill \cr {0,} \hfill & {{\rm for}\quad {{\bar \sigma} _{{\rm ii}}} \lt 0\;{\rm MPa}.} \hfill \cr}} \right.$$

2.5. Mechanical equilibrium

Assuming small deformations, the mechanical equilibrium can be described by the standard viscoelastic boundary value problem in the computational domain Ω as

(19) $${\sigma _{{\rm ij,j}}} + {b_{\rm i}} = 0,\quad {\rm in}\quad \Omega, $$

(20) $${\epsilon} _{{\rm mn}}^{\rm e} = \displaystyle{1 \over 2}({u_{{\rm m,n}}} + {u_{{\rm n,m}}}) - {\epsilon} _{{\rm mn}}^{\rm v} \quad {\rm in}\quad \Omega, $$

(21) $${\sigma _{{\rm ij}}} = C_{{\rm ijmn}}^{{\rm hd}} {\epsilon} _{{\rm mn}}^{\rm e}, \quad {\rm in}\quad \Omega, $$

(22) $${\sigma _{{\rm ij}}}{n_{\rm j}} = {\bar t_{\rm i}}\quad {\rm on}\quad {\Gamma _{\rm t}},$$

(23) $${u_{\rm i}} = {\bar u_{\rm i}}\quad {\rm on}\quad {\Gamma _{\rm u}},$$

where b _i is the body forces vector; ${\bar u_{\rm i}}$ denotes any prescribed displacement conditions corresponding to free slip or zero slip on the domain boundary Γ_u, ${\bar t_{\rm i}}$ denotes any prescribed traction conditions corresponding to seawater pressure on the domain boundary Γ_t; and n _j denotes the outward normal to the boundary Γ_t. Equation (19) is the static equilibrium equation in solid mechanics, which resembles the stationary Stokes approximation from the fluid mechanics point of view. The viscous strain ${\epsilon} _{{\rm mn}}^{\rm v} $ in Eqn (20) is calculated from the evolution law in Eqn (15), which defines ${\dot \epsilon} _{{\rm mn}}^{\rm v} $ . In Eqn (21), $C_{{\rm ijmn}}^{{\rm hd}} $ denotes the hydro-damage modified fourth-order elasticity tensor, defined as

(24) $$C_{{\rm ijmn}}^{{\rm hd}} = ((1 - D){\delta _{{\rm km}}}{\delta _{{\rm ln}}} - {D^{{\rm hyd}}}{\delta _{{\rm kl}}}{\delta _{{\rm mn}}}){C_{{\rm ijkl}}}.$$

Using the relations C _ijkl δ _kl = (E/(1 − 2ν))δ _ij = 3κδ _ij and ${\bar \sigma _{{\rm qq}}} = 3\kappa {\epsilon} _{{\rm pp}}^{\rm e} $ , the above equation can be derived from Eqn (13) as follows:

(25) $$\eqalign{{\sigma _{{\rm ij}}} &= (1 - D){C_{{\rm ijkl}}}{\epsilon} _{{\rm kl}}^{\rm e} - {\rho _{\rm f}}g\langle z\rangle D{\delta _{{\rm ij}}}, \cr & = (1 - D){C_{{\rm ijkl}}}{\epsilon} _{{\rm kl}}^{\rm e} - \displaystyle{{{\rho _f}g\langle z\rangle D} \over {3\kappa}} {C_{{\rm ijkl}}}{\delta _{{\rm kl}}} \cr & = (1 - D){C_{{\rm ijkl}}}{\epsilon} _{{\rm kl}}^{\rm e} - \displaystyle{{{\rho _f}g\langle z\rangle D} \over {3\kappa {\epsilon} _{{\rm pp}}^{\rm e}}} {\epsilon} _{{\rm pp}}^{\rm e} {C_{{\rm ijkl}}}{\delta _{{\rm kl}}}, \cr &= (1 - D){C_{{\rm ijkl}}}{\delta _{{\rm km}}}{\delta _{{\rm ln}}}{\epsilon} _{{\rm mn}}^{\rm e} - \displaystyle{{{\rho _{\rm f}}g\langle z\rangle D} \over {{{\bar \sigma} _{{\rm qq}}}}}{C_{{\rm ijkl}}}{\delta _{{\rm kl}}}{\delta _{{\rm mn}}}\epsilon _{{\rm mn}}^{\rm e} \cr &= \left( {(1 - D){\delta _{{\rm km}}}{\delta _{{\rm ln}}} - \displaystyle{{{\rho _{\rm f}}g\langle z\rangle D} \over {{{\bar \sigma} _{{\rm qq}}}}}{\delta _{{\rm kl}}}{\delta _{{\rm mn}}}} \right){C_{{\rm ijkl}}}{\epsilon} _{{\rm mn}}^{\rm e}, \cr & = ((1 - D){\delta _{{\rm km}}}{\delta _{{\rm ln}}} - {D^{{\rm hyd}}}{\delta _{{\rm kl}}}{\delta _{{\rm mn}}}){C_{{\rm ijkl}}}{\epsilon} _{{\rm mn}}^{\rm e},} $$

where κ denotes the bulk modulus of elasticity and the hydrostatic or hydraulic damage component D ^hyd is defined as

(26) $${D^{{\rm hyd}}} = \displaystyle{{{\rho _{\rm f}}g\langle z\rangle D} \over {{{\bar \sigma} _{{\rm qq}}}}}.$$

Evidently, the hydraulic damage D ^hyd is nonzero only when there is some existing damage and is proportional to the ratio of the effective fluid pore pressure P ^h D = ρ _f g〈z〉D and solid matrix pressure ${\bar \sigma _{{\rm qq}}} = 3\kappa {\epsilon} _{{\rm pp}}^{\rm e} $ . Subjected to plane strain assumptions, the solution to the nonlinear boundary value problem defined by Eqns (19–23) is obtained using the Galerkin FEM detailed by Duddu and Waisman (Reference Duddu and Waisman2012); Duddu and others (Reference Duddu, Bassis and Waisman2013) and Duddu and Waisman (Reference Duddu and Waisman2013). In the present finite element implementation, four-node bilinear quadrilateral elements were used to discretize the unknown displacement field and the four-point Gauss quadrature rule is used for integration. The internal state and history variables (e.g. damage, viscous strains) are stored at the quadrature points and an explicit forward Euler scheme is used to update these variables in time. We note that choosing higher order elements or finer resolutions is generally recommended in case higher accuracy is required.

3.1. Model geometry and parameters

We idealize the geometry of grounded marine-terminating glaciers as rectangular slabs of ice in contact with water, as shown in Figure 2. We apply a free slip boundary condition in the horizontal direction at the bottom edge of the slab (assuming minimal friction from the bed) and in the vertical direction on the left edge of the slab (where the ice slab is connected to the larger glacier). We apply a fixed (or zero displacement) boundary condition in the horizontal direction on the left edge of the slab. Assuming the influx of ice is independent of depth, this set of boundary conditions is translationally invariant in the horizontal direction and hence independent of the inflow velocity. Furthermore, the choice of zero displacement or velocity boundary condition is justified because we are interested in calculating the stress and deformation rates defined by the displacement or velocity gradients, thus, they are independent of the inflow velocity or displacement condition at the left edge.

Fig. 2. Schematic drawing of the idealized grounded ice slab with dimensions and boundary conditions.

We denote the depth of the surface and basal crevasses by d _s and d _b respectively, and the initial ice thickness by H. The initial notches play the role of pre-existing weaknesses or starter cracks in linear elastic fracture mechanics and provide the seeds for localized damage propagation. This assumption prevents the growth of non-physical damage areas on the top of the slab in areas where the FEM discretization may be coarse to reduce the computational burden. Alternative crack or damage initiation schemes can be employed (e.g. seeding the glacier with random defects), but our scheme (i.e. seeding crevasses using notches) allows us to easily perform model sensitivity studies (Duddu and Waisman, Reference Duddu and Waisman2013). Additionally, in all the following numerical examples, once damage localizes near the initial notches, we allow hydrostatic (or hydraulic) damage to occur only in the vicinity of the crack path. The slab length L = 2500 m is set to five times the ice thickness H = 500 m and the initial surface or basal crevasses (notches) are prescribed at midlength, to avoid edge effects at the inflow and outflow boundaries. The depths of the initial crevasses (notch) are set to 8% of the initial slab thickness H. The piezometric head or hydraulic head in surface crevasses h _s is defined as the height of the water column measured from the top edge of the slab, consequently, the water pressure at the bottom of the surface crevasse is proportional to (h _s + d _s). We specifically allow for h _s > 0 to examine the effect of a supra-glacial lake filling a crevasse, although it simulates an unphysical example because we do not model the lake nor the topography necessary to sustain the lake. The piezometric head at the bottom of the basal crevasse is assumed to be equal to the height of water level on the right edge of the slab denoted by h _w. All the relevant material and model parameters listed in Table 1 are assumed from Duddu and Waisman (Reference Duddu and Waisman2012) for ice at −10°C. The homogeneous ice density ρ _i = 910 kg m⁻³ and water density ρ _f = 1000 kg m⁻³. The value of the damage coefficient parameter B in Eqn (16), controlling the damage rate, is varied in the numerical studies so as to assess model sensitivity.

Table 1. Values of mechanical and damage parameters for ice at −10°C

The parameter A is the viscosity coefficient retrieved from Duddu and others (Reference Duddu, Bassis and Waisman2013) by taking A = (3/2)K _N; where K _N is the viscosity coefficient defined by Duddu and others (Reference Duddu, Bassis and Waisman2013).

3.2. Benchmark simulation

To demonstrate the capability of the proposed damage model to consistently calculate the hydraulic forces on the crevasse walls, we consider a rectangular ice slab (2500 m × 500 m) initialized with a surface and a basal crevasse in two different approaches. In the first approach, the crevasses are defined by notches and the hydraulic forces at the finite element nodes lying on the crevasse walls are calculated from the hydrostatic pressure distribution using the line integral, $\int_\Gamma {P^{\rm h}}d\Gamma $ , as indicated in Figure 3. Thus, in the first approach the geometrical features of the crevasses are explicitly meshed and the corresponding results provide a reference solution or benchmark. In the second approach, the crevasses are defined by fully damaged elements by specifying D = 1 inside the crevasse zone and the hydraulic forces at the finite element nodes lying on the crevasse walls are calculated from the stress distribution using the area integral, ${\int_\Omega} {\rho _{\rm f}}g\langle z\rangle D{\delta _{{\rm ij}}}d\Omega $ (as given by the second term in Eqn (13)). Thus, in the second approach the geometrical features of the crevasses are implicitly defined by the damage variable and the results are compared with those obtained from the first approach. The following parameters were used in this study: crevasse depths d _s = d _b = 25 m; the piezometric head h _s = 0 and h _w = 0.5H. The total hydraulic forces calculated from both approaches on either side of the crevasse are the same: 3126.9 kN for the surface crevasse and 59 411.8 kN for the basal crevasse. The horizontal stress distributions computed using the two approaches, shown in Figure 4, are identical to within numerical error. Thus, this benchmark investigation indicates that our damage model is capable of accurately describing the stress state of an ice slab with fully damaged crevasse (i.e. when D = 1). This example also demonstrates the main advantage of the continuum damage mechanics description that completely eliminates the need for remeshing as the fracture propagates.

Fig. 3. Hydraulic pressure distribution on the surface and basal crevasses, assuming complete ice material failure (D = 1) of the elements in the crevasse zone, which is represented by a notch in this figure. The values of the hydraulic pressure P ^h are given at the bottom and top of the surface and basal crevasses.

Fig. 4. Snapshot of horizontal stress, σ _xx , contours in the linear elastic configuration using two methods; (a) models crevasses as notches and hydraulic forces are applied as nodal forces; (b) represents the equivalent proposed damage mechanics technique and (c) the stress σ _xx profile along the slab centerline where the height of the material point is measured from the bottom of the slab.

3.3. Effect of hydrodamage on surface crevasse propagation

We first conducted several simulations to investigate the effect of hydrodamage on the depth d _s to which surface crevasses penetrate in relation to the seawater depth h _w. We consider three different values of h _w/H = {0, 0.5, 0.8} and take h _s = 0 to activate hydraulic damage under wet conditions. Figure 5 shows snapshots of damage contours (red zones indicate completely damage ice) for h _w/H = 0.8 and damage coefficient B = 10⁻⁵ MPa^−r s⁻¹. These simulation results emphasize the localized nature of crevasse propagation in glaciers driven by stress concentrations near pre-existing defects. Next, we conducted model sensitivity studies by varying the value of the damage coefficient B = {10⁻³, 10⁻⁴, 10⁻⁵} MPa^−r s⁻¹ and the corresponding time steps used in the analysis are dt = {0.1, 1, 10}s, respectively. The time step is chosen sufficiently small (on the order of seconds) to ensure accuracy and stability of the explicit time update scheme used for computing damage and viscous strain evolution. Crevasse depths were computed under dry (no hydrodamage – NHD) and wet (hydrodamage – HD) conditions. In Figure 6, the normalized surface crevasse depths (d _s/H) are plotted against simulation time (h) for different normalized seawater depths (h _w/H). Note that the depth of the surface crevasse d _s at a given time is measured as the vertical distance from the top of the slab to the farthest finite element node where the damage exceeds the critical damage value, that is, D >D _cr (Table 1). The following conclusions can be drawn from Figure 6:

1. (1) Water-filled surface crevasses experience more damage at any particular time and propagate to greater depths (for a given color compare the solid and dashed curves in Fig. 6), including full-depth fractures indicating calving events (i.e. d _s/H = 1); these results conform with Nye zero-stress results in Weertman (Reference Weertman1973) and linear elastic fracture mechanics (LEFM) results in Van der Veen (Reference Van der Veen1998b);

2. (2) The value of the damage coefficient B does not significantly change the final crevasse depth (compare the different colored dashed curves in Fig. 6), but it affects the rate at which crevasses propagates; consequently, it affects the time interval between calving events. However, model predicted damage (or crevasse) propagation rates are poorly calibrated by field or laboratory experiments. This result also demonstrates that the small strain assumption does not influence the conclusions drawn from our modeling studies; because similar crevasse penetration depths are retrieved for the case of B = 10⁻³ (where the accumulated viscous strains are small) and the case of B = 10⁻⁵ (where the accumulated viscous strains are larger).

Fig. 5. Snapshot of damage contours at different time steps for the isolated surface crevasse propagation model with h _w/H = 0.8. These results were simulated using B = 10⁻⁵ MPa^−r s⁻¹.

Fig. 6. The evolution of surface crevasse with time under different values of near terminus water depth h _w. The figures show simulation results for different values of B, the parameter in Eqn (16). The tag ‘HD’ in the legend means including Hydraulic Damage and ‘NHD’ means No Hydraulic Damage.

We next performed a sequence of simulations to investigate the stability of marine-terminating glaciers in relation to the piezometric head h _s in surface crevasses. We consider three different values of h _s = {0, 25, 50} m (where h _s >0 corresponds to the presence of a supraglacial lake) and recorded the temporal evolution of surface crevasses for different seawater depths h _w. In Figure 7, we plot the normalized maximum (or final) crevasse depth $d_{\rm s}^{{\rm max}} /H$ and the total time (h) elapsed till maximum crevasse depth is attained as a function of the normalized seawater depth h _w/H, under dry conditions and under wet conditions for three different values of piezometric head h _s. These simulations illustrate that:

1. (1) Under dry conditions, through-thickness surface crevasse propagation is not observed, regardless of the seawater level (blue curve in Fig. 7a); whereas, under wet conditions through thickness surface crevasse propagation always occurs except when the seawater level is sufficiently high (h _w > 0.8, as indicated by green, red and brown curves in Fig. 7). Thus, meltwater in surface crevasses destabilizes the glacier by driving through thickness crevasse propagation. These conclusions agree with the Nye zero-stress model in Weertman (Reference Weertman1973) that water-filled surface crevasses are highly likely to reach the bottom of the slab.

2. (2) An increase of seawater height generally decreases the rate of crevasse propagation (more pronounced when h _w > 0.5 in Figs 6, 7), consequently, it increases the total time elapsed till maximum crevasse depth is attained. Thus, the seawater level has a stabilizing effect on crevasse propagation as it applies a compressive crack-closing pressure. These results provide a qualitative measure of the conditions that lead to faster versus slower crevasse propagation, although the quantitative damage propagation rate remains poorly calibrated.

Fig. 7. Final crevasse depths ratios $d_{\rm s}^{{\rm max}} $ and corresponding simulation times for different values of h _s. The tag ‘NHD’ means No Hydraulic Damage. The points that are marked in the time plot with t → ∞ are those showing no crevasse propagation at all i.e. $d_{\rm s}^{{\rm max}} /H = 0$ on the left plot.

The numerical Nye-zero depth is calculated as the depth of the material point (from the top surface) at which the horizontal tensile stress vanishes (σ _xx  = 0) and the values were recorded prior to damage propagation for all the simulations. The final crevasse depths retrieved from our damage simulations were always within 10% to those predicted by the Nye zero-stress model.

3.4. Effect of hydrodamage on surface and basal crevasse propagation

We next investigated the effect of hydrodamage on the simultaneous propagation of surface and basal crevasses. Unless the ocean water level is sufficiently high (h _w/H > 0.8), surface crevasses always form in our simulation and it is not possible to have isolated basal crevasses without surface crevasses. This could be changed by setting a finite threshold for the largest principal stress as opposed to using a zero threshold as we did here. We find that, in accordance with the findings of Nick and others (Reference Nick, Van der Veen, Vieli and Benn2010) and Bassis and Walker (Reference Bassis and Walker2012), basal crevasses do not propagate unless they are water filled.

Through finite element simulation, we estimated surface and basal crevasse depths by varying the seawater height h _w/H = {0, 0.25, 0.6, 0.7, 0.8, 0.9} for two scenarios: (1) only basal crevasse is water filled and surface crevasse is dry, and (2) both surface and basal crevasses are water filled. The plots of normalized maximum (or final) crevasse depths ( $d_{\rm s}^{{\rm max}} $ and $d_{\rm b}^{{\rm max}} $ ) versus the normalized seawater height are shown in Figures 8, 9 for the two scenarios, respectively. In the case of a dry surface crevasse and water-filled basal crevasse (Fig. 8), our results indicate that the maximum total crevasse depth $(d_{\rm s}^{{\rm max}} + d_{\rm b}^{{\rm max}} )$ decreases as the seawater level increases until h _w/H = 0.8, but then increases as the seawater level further increases, h _w/H > 0.8. This is because the maximum surface crevasse depth decreases as seawater level increases, whereas, the maximum basal crevasse depth becomes significant only at high seawater levels, when the water pressure is sufficient to induce a tensile stress at the basal crack tips. Thus, our simulation results in Figure 8 are in good agreement with those published by Bassis and Walker (Reference Bassis and Walker2012). In Figures 8, 9, we plotted the time taken for the crevasses to reach the maximum (or final) depth as function of h _w/H. These results show that the crevasses propagation times are smallest for h _w/H = 0 and h _w/H = 0.9, indicating that the corresponding crevasse propagation rates are the largest for h _w/H = 0 and h _w/H = 0.9, which is attributed to the existence of higher tensile stresses at the surface and basal crack tips, respectively.

Fig. 8. Final crevasses' depths ( $d_{\rm s}^{{\rm max}} $ for surface and $d_{\rm b}^{{\rm max}} $ for basal) and corresponding simulation times under different values of near terminus water depth h _w; in these simulations, only the basal crevasses are water filled while the surface crevasses are dry. These results were simulated using B = 10⁻⁴ MPa^−r s⁻¹.

Fig. 9. Final crevasses' depths ( $d_{\rm s}^{{\rm max}} $ for surface and $d_{\rm b}^{{\rm max}} $ for basal) and corresponding simulation times under different values of side water pressure h _w; in these simulations, both surface and basal crevasses are water filled. These results were simulated using B = 10⁻⁴ MPa^−r s⁻¹. Note that final total crevasse depth is always larger when both basal and surfaces are water filled (compared blue dashed lines in Figures 8a and 9a), thus indicating a mutually positive effect.

An important point to note is that surface and basal crevasses propagate to a greater depth (or height) when both are water filled (compare blue dashed lines in Figures 8a and 9a), thus indicating a mutually positive effect. To further investigate the effect of water-filled surface crevasses on basal crevasse propagation and vice-versa, we plot the temporal evolution of surface and basal crevasses in Figure 10 for water level h _w/H = 0.25. The results in Figure 10 illustrate that the basal crevasse propagation is triggered when the surface crevasse approaches the bottom of the slab and alters the stress at the basal crevasse tip. The main conclusion of our study is consistent with previous studies: glaciers with dry surface crevasses are more stable than those with water-filled surface crevasses, and the situation worsens when the basal crevasses are water filled.

Fig. 10. Damage propagation for the case of h _w/H = 0.25 from Figure 9 (surface and basal crevasses are water filled). The top plot shows crevasses propagation with time and the following are snapshots of damage contours at different time steps. These results were simulated using B = 10⁻⁴ MPa^−r s⁻¹.

In this study, we did not exclusively investigate conditions that enable the propagation of only a basal crevasse without a surface crevasse. Because we assume that glacier ice has zero tensile strength, surface crevasses will always form while simulating an extending glacier and so it is not possible to have isolated basal crevasses without surface crevasses. However, the model can account for the tensile strength for ice by specifying a stress threshold for damage initiation (Pralong and Funk, Reference Pralong and Funk2005), which disables the formation of surface crevasses at very low deformation rates. Similarly, including a sliding law instead of a free-slip boundary condition would quantitatively alter our simulation results. Nonetheless, the hydraulic damage methodology we have developed is directly applicable to more complex geometries, boundary conditions and rheologies.