2.1 Ice flow model and rheology

We used the finite-element library libMesh (Kirk et al., 2006) to implement the Stokes equations for continuum momentum and mass conservation:

where σ is the Cauchy stress tensor, ρ_i the ice density, g the gravitational force vector and u the velocity vector. As we assume ice to be incompressible and isotropic, the Cauchy stress tensor can be decomposed into an isotropic and a deviatoric part σ^′:

where ${\mathit{\sigma }}_{\mathrm{m}}=\frac{\mathrm{1}}{\mathrm{3}}\text{tr}\left(\mathit{\sigma }\right)=\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\sigma }}_{ii}$ is the isotropic mean stress and I the identity matrix. The ice rheology is described as viscous power-law fluid (Glen's flow law), linking the deviatoric stress tensor σ^′ to the strain rate tensor $\dot{\mathit{\epsilon }}$:

The effective shear viscosity η is defined as

where ${\dot{\mathit{\epsilon }}}_{\mathrm{e}}=(\frac{\mathrm{1}}{\mathrm{2}}{\dot{\mathit{\epsilon }}}_{ij}{\dot{\mathit{\epsilon }}}_{ij}{)}^{\frac{\mathrm{1}}{\mathrm{2}}}$ is the effective strain rate, A the fluidity parameter, n=3 the power-law exponent and κ_ε is a finite strain rate parameter included to avoid infinite viscosity at low stresses (Greve and Blatter, 2009).

The model domain was discretized with second-order nine-node quadrangle elements with Galerkin weighting. Model variables are approximated with a second-order approximation for the velocities u and w and a first-order approximation for the mean stress σ_m (forming a LBB-stable set). The accuracy of the solution was improved with adaptive mesh refinement near the calving front. The Stokes equations with the nonlinear rheology were solved with the PETSc nonlinear solver SNES to a relative accuracy of 10⁻⁴ (Balay et al., 2008).

2.2 Model geometry and scaling

We used a two-dimensional version of the model to conduct the geometrical tests, as illustrated in Fig. 2. The geometry is defined in a Cartesian coordinate system with horizontal axis x and vertical axis z with origin at sea level at the calving front (where x=0). The ice moves from right to left. The idealized glacier geometry used in all model experiments consists of a block of ice resting on a flat bed with a characteristic length L=2000 m and a characteristic ice thickness H=200 m. The domain was discretized with 20 elements in the vertical and 200 elements in the horizontal which, after mesh refinement, led to a spatial resolution of 2.5 m in the terminus area.

All numerical results are scalable with reference values for ice thickness H_ref and overburden stress σ_ref and are therefore independent of the geometrical extent. This validity of the scaling was tested by running the model for different ice thicknesses, which recovered identical flow and stress results. The velocity scale u_ref was chosen as the vertical surface velocity caused by uniaxial confined compression in pure shear of an ice block under its own weight (Cuffey and Paterson, 2010):

The coordinates and the water depth at the calving front H_w are scaled by the ice thickness H_ref:

All stress and velocity components are scaled according to

Figure 2Geometry of the idealized grounded glacier. α is the slope angle of the calving front above the vertical cliff.

Download

2.3 Boundary conditions

The upper surface of the glacier was described as a traction-free surface boundary. Basal motion was parametrized with a slipperiness coefficient C, which relates the basal velocity u_b with basal shear stress τ_b (Gudmundsson and Raymond, 2008; Ryser et al., 2014):

This boundary condition was implemented as a two-element layer with constant viscosity ${\mathit{\eta }}_{\mathrm{s}}=h_{\mathrm{s}}/C$, which was added at the bottom of the model domain representing the glacier. At the lower boundary of this “sediment layer”, a Dirichlet boundary condition with zero velocity ($\mathit{u}=v=\mathrm{0}$) was imposed. A layer thickness of h_s=10 m was chosen, although tests with varying h_s showed no significant differences. This simple approach allowed us to capture the physical processes that are relevant to this study. In the case of vanishing basal motion the two-element layer was not used for the computation, and Dirichlet boundary conditions ($\mathit{u}=v=\mathrm{0}$) were imposed directly at the bottom of the model domain representing the glacier.

At the calving front a normal stress boundary condition was imposed below the water level, while the surface above water was kept stress-free. The stress boundary condition thus reads

where σ_nn and ${\mathit{\sigma }}_{n{t}_i}$ are the normal and tangential tractions applied on the calving front (σ_nn is negative, i.e., compressive since z<0 below water) and ρ_w is water density (Table 1).

At the upstream boundary of the glacier domain velocities were fixed to zero. Additional modeling experiments showed that different values for this upstream boundary condition do not affect the results of the analysis.

2.4 Sensitivity analysis strategy

The stress state and flow field near the calving front is analyzed in three suites of numerical experiments that investigate the effect of variations in relative water level ω, the slope of the calving front and basal motion.

The water level sensitivity experiments were performed for relative water levels $\mathit{\omega }=\mathrm{0},\mathrm{0.25},\mathrm{0.5},\mathrm{0.75},\mathrm{0.85}$ and ${\mathit{\omega }}_{\mathrm{f}}=\frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}$, where the last value is the relative water level at flotation. The calving front for this experiment was vertical and the bottom boundary without sliding (i.e., zero velocity Dirichlet boundary condition). All these experiments were undertaken with both the density of ocean water (ρ_w=1028 kg m⁻³) and freshwater (ρ_w=1000 kg m⁻³).

The calving front slope sensitivity experiments were performed on a geometry with the upper part of the calving front reclining at various angles. The lower 25 % of the calving front height was set vertical, and the upper part inclined at angles from 90, 75, 60 and 45^∘, until it reached the maximum surface height (see Fig. 2 for illustration). This particular geometrical setup was chosen to represent a simplified geometry of Eqip Sermia, which has a 50 m high vertical cliff at the bottom with a 45^∘ inclined slope up to the top at 200 m. For this experiment, the relative water level was set to ω=0 and the sliding velocity was set to zero.

The bed slipperiness sensitivity experiments were performed on a block geometry with a vertical calving front and a relative water level ω=0.5. The basal slipperiness coefficient C was varied from 0 to 1000 $\mathrm{m}{\mathrm{MPa}}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$ with 333 $\mathrm{m}{\mathrm{MPa}}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$ increments. A slipperiness of 1000 $\mathrm{m}{\mathrm{MPa}}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$ corresponds to a sliding speed of 300 m a⁻¹ for a typical tidewater outlet glacier in Greenland with a driving stress of 0.3 MPa.

2.5 Stress invariant combinations

Any criterion for fracture propagation or damage evolution should be independent of the choice of coordinate system and can therefore be expressed as a function of the invariants and eigenvalues of the stress tensor. Hayhurst (1972) proposed a linear combination of three stress invariants to describe the creep rupture of ductile and brittle materials under multi-axial states of stress. The invariants chosen were maximum principal stress σ₁, first stress invariant $I_{\mathrm{1}}={\mathit{\sigma }}_{\mathrm{m}}=\frac{\mathrm{1}}{\mathrm{3}}{\mathit{\sigma }}_{ii}$ and the von Mises stress $J_{\mathrm{2}}={\mathit{\sigma }}_{\mathrm{e}}=(\frac{\mathrm{3}}{\mathrm{2}}{\mathit{\sigma }}_{ij}^{\prime }{\mathit{\sigma }}_{ij}^{\prime }{)}^{\frac{\mathrm{1}}{\mathrm{2}}}$ to form the stress combination

where the weights α, β and γ fulfill the conditions

The Hayhurst stress χ_H has been used as a criterion for the initiation and evolution of damage in several glaciological studies (Duddu and Waisman, 2012, 2013; Duddu et al., 2013; Mobasher et al., 2016; Pralong and Funk, 2005; Pralong et al., 2003).

Figure 3Sensitivity experiment results for varying water depth. (a) Scaled Hayhurst stress distribution. (b) Scaled horizontal velocity distributions. (c) Scaled Hayhurst stress along the surface. (d) Scaled horizontal velocity magnitude along the surface. In panels (a, b), the subplots show increasing water depth from the bottom to the top (water level at $\widehat{z}=\mathrm{0}$). Solid and dashed lines in panels (c, d) correspond to experiments with sea- and freshwater densities, respectively.

Download

Figure 4Scaled velocities along the vertical face of the calving front (solid lines) for different relative water levels. Horizontal line markers show the relative water level for each curve.

Download

Figure 5Sensitivity experiment results for varying calving front slope. Panels (a–d) are the same as Fig. 3. In panels (a, b), the subplots show decreasing calving front slopes from the bottom to the top. In panel (c), the local minimum of stress close to the calving front is located where the front reaches its maximal height. In panel (d), vertical lines on the curves for inclined fronts mark the distance at which the maximal surface height is reached.

Download

To investigate the full spectrum of possible stress states that lead to the initiation of damage, we investigated linear combinations of five stress invariants: ${\mathit{\sigma }}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{e}},{\mathit{\sigma }}_{\mathrm{m}}$ and additionally the third invariant of the stress tensor I₃=det(σ) and third invariant of the deviatoric stress tensor $J_{\mathrm{3}}=\text{det}\left({\mathit{\sigma }}^{\mathbf{\prime }}\right)$. This extended linear combination reads

with weights $\mathit{\alpha },\mathit{\beta },\mathit{\gamma },\mathit{\varphi }$ and μ that fulfill the conditions

We performed a sensitivity analysis based on the five stress invariants of Eq. (14) by systematically varying the weights with 0.1 increments (Eq. 15).

Figure 6Sensitivity experiment results for varying bed slipperiness C. Panels (a–d) are the same as Fig. 3. In panels (a, b), the subplots show increasing bed slipperiness from the bottom to the top. Units for bed slipperiness C are m MPa a⁻¹.

Download

3.1 Sensitivity analyses

All sensitivity experiment results shown in Figs. 3, 5 and 6 exhibit similar velocity and stress patterns. The effects of varying water level, basal slipperiness and calving front slope on the stress field are displayed as Hayhurst stress with parameters chosen according to Pralong and Funk (2005) (Table 1). In general, the modeled velocities and stresses increase towards the calving front, with a local stress maximum at the surface that is located less than one ice thickness upstream of the calving front. This zone of high stress extends diagonally down towards the calving front where it has a second local maximum closely above the water level. For experiments with a relatively low water level, the absolute maxima in stress are found at the bottom of the calving face.

3.1.1 Water level height

The depth of the water at the calving front significantly impacts the stress regime and consequently the ice flow pattern and magnitude near the terminus. The effect of different water depths on the stress field is displayed as Hayhurst stress in Fig. 3a and c.

For a reduction in the relative water level from ω=ω_f to ω=0 the maximum Hayhurst stress at the surface increases from 0.08 to 0.42 σ_ref and the location of the stress peak at the surface moves from 0.1 to 0.5 H upstream of the front, whereas the Hayhurst stress at the vertical calving front increases from 0.15 to 0.81 σ_ref. Interestingly, the local maxima at the front are always located near the water level. Further, the position of the global stress maximum for low water levels (below 0.25) is found at the bottom of the calving front instead of the surface (Table 2).

Figures 3b, d and 4 illustrate how velocities close to the calving front increase by more than 1 order of magnitude when the water level is decreased from near flotation (ω=0.85) to shallow water (ω=0.25). Note that for all water depths the velocities are only affected up to approximately 2.5 ice thicknesses upstream from the front.

Extrusion flow, a velocity pattern for which maximum horizontal velocity occurs below the surface (Waddington, 2010), is clearly visible in Fig. 4 in the vicinity of the calving front for the low water level cases. This pattern of extrusion flow near the terminus was also observed by Hanson and Hooke (2000) and Leysinger-Vieli and Gudmundsson (2004).

In summary, increasing relative water depth leads to decreased flow velocities and lower stresses and moves the peak of the Hayhurst stress at the surface closer to the front.

Table 2Maximum scaled Hayhurst stress and velocity for water depth experiments. The s and f letters indicate whether the scaled Hayhurst stress maxima were found at the surface or at the bottom of the calving front, respectively.

Download Print Version | Download XLSX

Table 3Maximum Hayhurst stress and velocity for calving front slope experiments.

Download Print Version | Download XLSX

Table 4Maximum Hayhurst stress and velocity for bed slipperiness coefficient experiments.

Download Print Version | Download XLSX

3.1.2 Calving front slope

Results from the sensitivity experiment on calving front slope displayed in Fig. 5 show large variations in stresses and flow speeds. Maximum Hayhurst stresses are found at the bottom of the calving front for all cases ranging from 0.57 to 0.81 σ_ref for slope angles between 45 and 90^∘ (Table 3). A second, local maximum occurs at the surface behind the end of the slope, but the magnitude strongly decreases with decreasing slope. The maximum velocity for a 45^∘ slope is ∼4 times smaller than for a vertical calving front (Fig. 5d). Thus, as the calving front gets steeper, the stresses and the velocities increase. Again, the peak in Hayhurst stress at the surface moves further upstream as the calving front is becoming more gentle and a further local stress maximum occurs along the sloped surface. Moreover, the velocities along the surface peak not towards the front corner as in the vertical front case but rather towards the bottom of the sloped surface, which is another sign of extrusion flow.

Figure 7Sensitivity experiments result for an inclined surface (c, f), reverse bed (a, d) and the simple rectangular geometry (b, e). The left and right panels show the scaled Hayhurst stress and velocity distributions, respectively.

Download

3.1.4 Bed and surface slope

In the modeling presented so far we used a glacier geometry with horizontal surface and bed. Consequently the driving stress and hence velocities and stresses far upstream from the calving front are close to zero. In reality glaciers have a sloping surface. Therefore, we repeated some of the above experiments on a simple glacier geometry with a sloped bed and surface, a fixed cliff height and no sliding. Bed and surface slopes were chosen as −5 and 5^∘, respectively. Figure 7 illustrates the results: a reclining slope at the surface (i.e., surface height increasing towards the inland) with a flat bed leads to higher stresses and velocities upstream of the calving front as compared to the flat surface. However, the stress maximum and its location in the vicinity of the calving front remains almost identical (Fig. 7b, c). Similar results are obtained for a reverse bed slope with a flat surface (Fig. 7a, b).

To summarize, the stress and velocity fields in the vicinity of the calving front are only slightly altered for sloping bed and surface. It is, however, noteworthy that the reclining surface slope induces higher stresses near the surface, which could potentially induce crevassing and thus advect pre-damaged ice to the calving front.

Figure 8Combinations of five stress tensor invariants at the surface of an idealized glacier with a vertical calving front without water pressure and zero basal motion. Each black line represents a linear combination of five stress invariants. The blue envelope contains the maxima of all stress invariant combinations. The green triangle, red square and purple circle represent the maximum of the scaled Hayhurst stress ${\widehat{\mathit{\chi }}}_{\mathrm{H}}$, von Mises stress ${\widehat{\mathit{\sigma }}}_{\mathrm{e}}$ and maximum principal stress ${\widehat{\mathit{\sigma }}}_{\mathrm{1}}$, respectively.

Download

3.2 Stress invariant combinations

The Hayhurst stress, typically used as the driving force for damage evolution (Duddu and Waisman, 2012, 2013; Duddu et al., 2013; Mobasher et al., 2016; Pralong and Funk, 2005; Pralong et al., 2003), is not the only possible combination of objective stress measures. Here we attempt a systematic analysis of the possible stress invariant combinations (Eq. 14) and the corresponding locations of the stress maxima along the glacier surface. We illustrate this analysis at the example of the block geometry without any water pressure (ω=0) in Fig. 8, where all possible linear combinations of five stress invariants along the surface are displayed. While the stress combinations show a wide variety of curves, the maximum achievable stress states are dominated by the von Mises stress σ_e and the maximum principal stress σ₁, both of which contribute to the Hayhurst stress. Hence, these two stress invariants are likely the driving factors for material failure in the vicinity of the calving front. An important aspect illustrated in Fig. 8 is the horizontal position of the stress maximum, which is limited to x_max≃0.7 H_ref. This analysis thus suggests that a zone with maximum crevasse opening cannot be located in greater distance from the calving front than x_max for an idealized glacier without pre-damaged ice.

Figure 9Envelopes of stress invariant combinations at the surface of the idealized glacier with zero basal motion for varying relative water level ω. The green triangle, red square and purple circle represent the maximum of the scaled Hayhurst stress, von Mises stress and maximum principal stress, respectively, for each water depth.

Download

The magnitudes and positions of the maximum stress invariant combinations for different relative water levels ω are shown in Fig. 9 (blue area corresponds to Fig. 8). The maximum stress for dry conditions (ω=0) is located ∼0.7 H_ref from the calving front with a maximum von Mises stress of ∼ 0.45 σ_ref, whereas a water level close to flotation (ω=0.85) leads to a stress maximum of ∼0.15 σ_ref at ∼0.25H_ref from the calving front. Figure 9 clearly illustrates that water pressure at the calving front exerts a stabilizing effect on the calving front by both lowering the stresses and decreasing the distance from the calving front at which the stress maximum is located, as argued earlier by Bassis and Walker (2012).

The Hayhurst, maximum principal and von Mises stress distributions are shown in Figs. 3a, B1a and B1b, respectively.

4.1 Stress parametrization

The distribution along the glacier surface of scaled maximum principal stress ${\widehat{\mathit{\sigma }}}_{\mathrm{1}\mathrm{s}}$ is shown for different relative water levels in Fig. 10 (this is approximately the tensile stress along the surface, whereas Fig. 3c shows Hayhurst stress). The similarity in shape of these stress curves allows for an approximate representation by a function that depends on the relative water level ω alone:

where $\widehat{\widehat{x}}$ is a stretched and shifted version of the scaled (by ice thickness) horizontal coordinate $\widehat{x}$. This stretching function is somewhat cumbersome and is given in Appendix A. The extensional stress reaches the maximum at $\widehat{\widehat{x}}=\mathrm{1}$ (setting the derivative of Eq. 16 to 0) with magnitude ${\widehat{\mathit{\sigma }}}_{\mathrm{1},\mathrm{m}}=a\left(\mathit{\omega }\right)\exp (-\mathrm{1})$ and can be approximated by

It is interesting to note that the maximum extensional stress at the surface has a similar form as the mean deviatoric stress in Eq. (1) but is ∼60 % higher. The scaled horizontal position of the stress maximum can be approximated by

4.2 Analytical calving relation

Using the parametrizations of magnitude and position of the maximum extensional stress at the surface (Eqs. 17a and 18) the calving rate can be estimated under simple assumptions on crevasse formation.

One major assumption is that a large crevasse forms at the location of the maximum tensile surface stress where the ice is weakened until failure. Such crevassing seems realistic as both observations and model results show the formation of huge crevasses. When failure of the surface ice is complete, we assume that all ice in front of the crevasse is removed and a new calving front forms at the location of the crevasse. Here, we do not consider explicitly which processes are responsible for downward propagation of the crevasse. Several processes could be considered, such as bottom crevassing, hydro-fracturing by ponding water in surface crevasses, rapid elastic crevasse propagation (Krug et al., 2014), ice break-off in multiple steps (e.g., a surface slump, followed by subaqueous buoyant calving) or continued material fatigue due to tidal forcing. The proposed calving relation relies on the major assumption that processes responsible for ice break-off act on faster timescales than the formation of the surface crevasse and, therefore, that the calving process is uniquely determined by the time to failure at the surface stress maximum. Thus, the average calving rate ${\overline{u}}_{\mathrm{c}}$ can be calculated as the distance of the stress maximum divided by the time to failure T_f. In dimensional coordinates this is

Assuming further that crevasse formation can be described by isotropic damage formation with damage variable D, the stress in the damaged material is $\stackrel{\mathrm{̃}}{\mathit{\sigma }}=(\mathrm{1}-D{)}^{-\mathrm{1}}\mathit{\sigma }$ (Pralong, 2006; Pralong et al., 2003). The isotropic damage evolution relationship employed here is

where B is the rate factor for damage evolution, r and k are constants, σ₀ is the stress in the work zone and σ_th a stress threshold for damage creation. Integrating this relation over time, the time to failure, i.e., the time required for damage to evolve from an initial value D₀ to a critical value D_c, reads

We further assume that the stress in the work zone is the maximum tensile stress ${\mathit{\sigma }}_{\mathrm{0}}={\mathit{\sigma }}_{\mathrm{1},\mathrm{m}}$. Inserting the parametrizations for maximum tensile stress and stress maximum position (Eqs. 17c and 18) in the above relation yields

The term in square brackets is constant, and after renaming it the effective damage rate $\stackrel{\mathrm{̃}}{B}$ the expression reads

with parameter values $\stackrel{\mathrm{̃}}{B}=\mathrm{65}{\mathrm{MPa}}^{-r}{\mathrm{a}}^{-\mathrm{1}}$ and σ_th=0.17 MPa, which were determined from a calibration with data discussed below in Sect. 5.3. The parameter value r=0.43 is chosen according to Pralong and Funk (2005).

5.1 Sensitivity analyses

The stress intensity and therefore ice deformation rates are decreasing as the relative water level increases due to the pressure exerted by the water at the calving front. This feature is already captured by the depth-integrated extensional stress at the front (Eq. 1) and, in more detail, in the parametrized maximum extensional stress (Eq. 17a), illustrated in Fig. 10. In both cases the square dependence of the horizontal stress on relative water level controls fracture or damaging processes, the magnitude and rate of which depend linearly on the stress intensity.

In addition, the detailed modeling shows that the stress peak at the glacier surface moves upstream for lowering relative water level (Figs. 3 and 10), implying that crevasses are likely to open in greater distance from the calving front and leading to detachment of larger masses during calving.

A higher relative water level results in a more stable calving front (Bassis and Walker, 2012), which seems to be in contrast with the often-used relations which predict that calving rates increase with water depth (Brown et al., 1982; Hanson and Hooke, 2000; Meier and Post, 1987). In nature, however, glaciers terminating in deeper waters are also thicker and calve at higher rates as they experience higher absolute (unscaled) stresses. Furthermore, submarine frontal melting is likely to lead to higher calving rates by over-steepening of the front (O'Leary and Christoffersen, 2013), although the melt undercutting effect on calving rates seems to be limited (Cook et al., 2014; Krug et al., 2015).

Using freshwater instead of seawater at the calving front yields slightly higher stresses and velocities (Fig. 3c, d). This difference can be explained by the reduced back pressure applied by freshwater on the calving front, which results from a lower water density.

The model results demonstrate that reclining calving fronts lead to lower velocities and stresses and thereby implicitly confirm that inclined calving fronts should reach larger stable heights than vertical cliffs, as observed for example at Eqip Sermia (200 m high at 45^∘). This sensitivity analysis on front slope may, together with observational data on non-vertical calving fronts, provide constraints on parameters of ice resistance to failure. Further, the presence of extrusion flow along the reclining calving face of an idealized glacier was demonstrated. Such a velocity pattern has been observed and measured on an inclined slope at the northern front of Eqip Sermia (Lüthi et al., 2016) but is rarely discussed in modeling studies (Hanson and Hooke, 2000; Leysinger-Vieli and Gudmundsson, 2004).

Basal sliding leads to increased stresses at the surface throughout the computational domain. Thus, basal sliding may cause an onset of ice damaging and crevasse opening in a greater distance from the calving front (Fig. 6c). The velocity patterns in Fig. 6b show that the influence of bed slipperiness is only apparent in the proximity of the calving front, even for high sliding coefficients. Moreover, stress distributions are almost identical for all bed slipperiness experiments, which implies that basal sliding has a negligible effect on the stability of the calving front. Basal sliding adds a constant velocity at the bottom of the domain rather than affecting the velocity gradients. This result does not include any spatial variation in bed slipperiness, which would likely be caused by including a water-pressure-dependent sliding relation. Effective pressure (the difference between ice normal stress at the bottom and water pressure) typically decreases towards the calving front for real glaciers with sloping surface and may cause additional sliding towards the front, an effect that is not considered in this modeling effort.

For a sloping glacier surface the location and magnitude of the stress maximum in the vicinity of the calving front remain almost identical, as shown in Fig. 7. Similar results are obtained for a reverse bed slope with a flat surface, with a smaller influence on stresses and velocities than for the reclining surface. However, the effect on stresses and velocities upstream of the calving front is not visible for the reverse bed slope with a flat surface. This indicates that, for a glacier with a reclining surface slope, ice can potentially start damaging and forming crevasses at the surface far upstream from the calving front.

5.2 Calving relation

The proposed calving rate parametrization (Eq. 22) is simple and only requires two geometrical quantities: frontal ice thickness H and water depth H_w. The assumptions about the failure process are lumped into three parameters – $\stackrel{\mathrm{̃}}{B}$, σ_th and r – which can be determined by data calibration (Sect. 5.3). The parametrization exhibits many similarities with established calving relations but is formulated in terms of two quantities that are calculated by any ice flow model. It therefore is a drop-in replacement for other calving relations used in glacier models of different complexity.

The calving rate parametrization (Eq. 22) has some interesting properties, which are illustrated in Figs. 11 and 12. Holding constant the relative water depth, the absolute water depth or the ice thickness results in different calving laws:

• for constant relative water level ω the calving rate grows roughly like ${\overline{u}}_{\mathrm{c}}\propto H^{\mathrm{1}+r}$ (black and gray lines in Fig. 11);

• for constant absolute water depth H_w=ωH a fit shows that roughly ${\overline{u}}_{\mathrm{c}}\propto H^{\mathrm{1.25}}$ (red and orange lines in Fig. 11);

• for constant ice thickness the calving rate decreases with increasing relative water level (Fig. 12) roughly like

  $$\begin{array}{ll}{\displaystyle }{\overline{u}}_c& {\displaystyle }\propto \left(\mathrm{1}-{\mathit{\omega }}^{\mathrm{2.8}}\right){\left(\mathrm{1}-\mathrm{1.3}{\mathit{\omega }}^{\mathrm{2}}\right)}^r\\ {\displaystyle }& {\displaystyle }\simeq \left(\mathrm{1}-{\mathit{\omega }}^{\mathrm{2.8}}\right){\left(\mathrm{1}-{\left({\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{{\mathit{\rho }}_{\mathrm{i}}}}\mathit{\omega }\right)}^{\mathrm{2}}\right)}^r.\end{array}$$

The predicted calving rate for a given water depth depends on the thickness of the glacier, which is the result of the mass fluxes in the terminus area. Thus, calving rates depend on the surface evolution and hence the upstream dynamics of the glacier. The semi-empirical calving rate parametrization is therefore, in the sense of inclusion of upstream dynamics, similar to the position based calving models (Benn et al., 2007a, 2017; Nick et al., 2010; Todd and Christoffersen, 2014). The formulation as a calving rate also makes this parametrization relatively easy to use in larger-scale fixed grid models.

Figure 11Calving rates predicted by the parametrization in relation to ice thickness. Calving rates increase with increasing total ice thickness for a given water depth H_w=ωH (red and orange lines), relative water level ω (black and gray lines) or freeboard H−H_w (blue lines). Note that the gray line refers to a front at flotation.

Download

Figure 12Calving rates predicted by the parametrization as a function of relative water level. Calving rate decreases under increase of the relative water level ω for constant total ice thickness H.

Download

Figure 13Calving rates (m d⁻¹) predicted by the parametrization are shown as contours in dependence of H and ω. The hatched region indicates the states excluded by the maximum calving front criterion (Bassis and Walker, 2012). The gray area indicates states where the stress threshold σ_th precludes calving. Blue dots with numbers indicate calving rates determined from measurements, shown in Table 5.

Download

Figure 14Comparison of measured calving rates with predictions from the calving parametrization. The glacier names are abbreviated according to Table 5.

Download

5.3 Calibration of the parametrization

The calving rate parametrization (Eq. 22) contains three empirical parameters: $\stackrel{\mathrm{̃}}{B}=\mathrm{65}{\mathrm{MPa}}^{-r}{\mathrm{a}}^{-\mathrm{1}}$ and σ_th=0.17 MPa, which were obtained by calibration with data, and r=0.43, which is taken from Pralong and Funk (2005). Calving rates thus obtained are not very sensitive to the exact choice of parameter values, which are within the range of previous studies (Duddu and Waisman, 2012; Lliboutry, 2002; Vaughan, 1993).

Table 5Values of calving front height, water depth and calving rate for different glaciers.

Download Print Version | Download XLSX

To calibrate these parameter choices, data on calving rate, ice thickness and water depth for a wide variety of tidewater glaciers in the Arctic were collected. The data set covers the full range of water levels (relative and absolute), velocities and ice thicknesses that are found in Arctic tidewater glaciers. Unfortunately, many studies report only width-averaged data on calving front geometry and calving rate, which are not suitable for our proposed relation which relies on local stresses on a flowline. Only a limited set of point data on calving front geometries are available from the published literature from which total ice cliff thickness, water depth and calving rate can be obtained. For the calibration, we used the values shown in Table 5 from diverse data sources.

Contours of calving rates calculated with Eq. (22) are shown in Fig. 13 together with the maximum theoretical calving front height predicted by Bassis and Walker (2012). Figure 14 plots the same calving rate data against results from the parametrization. While a sizable spread of the data is visible, especially for low calving rates, the general agreement shows that the parametrization is well suited to estimate calving rates for this set of tidewater glaciers in the Arctic.

Note that the derivation of the parametrization is independent of the specific geometry or location of a tidewater glacier and thus the calibration is expected to be “global” and valid for any tidewater glacier.