3.1. A brief overview of calving laws

Calving rate is defined as the difference between the ice velocity at the glacier terminus and the change of glacier length over time, i.e.

(1)

where u _c is the calving rate, <$> is the vertically averaged ice velocity at the terminus, L is glacier length and t is time. This equation can be interpreted in two ways (Reference Van der VeenVan der Veen, 1996): (1) a forward approach, in which the calving rate is estimated from independent environmental variables and then the changes in terminus position are determined from calving rate and ice velocity (e.g. Reference Siegert and DowdeswellSiegert and Dowdeswell, 2004) and (2) an inverse approach, in which the calving rate is determined from ice velocity and changes in front position (e.g. Reference Van der VeenVan der Veen, 1996, Reference Van der Veen2002). Various ‘calving laws’ have been proposed following either approach. In what follows, we outline the major points of these calving laws. A more detailed account on the subject is given in the thorough review by Reference Benn, Warren and MottramBenn and others (2007b).

Most of the early calving laws were based on the forward approach, focusing on the direct estimate of the calving rate, relating it to some independent variable and then using it, together with the glacier velocity, to infer changes in ice-front position. Two main approaches were followed in the choice for the independent variable. The first (Reference Brown, Meier and PostBrown and others, 1982; Reference Pelto and WarrenPelto and Warren, 1991) considers the calving rate to depend linearly on the water depth at the calving front, with the particular linear relationship being based on the fit to field observations. The main problem with this approach is that the empirical relationships between water depth and calving rate vary greatly from one glacier to another (Reference HaresignHaresign, 2004), are quite different for tidewater and freshwater glaciers (Reference Funk and RöthlisbergerFunk and Röthlisberger, 1989) and also vary with time (Reference Van der VeenVan der Veen, 1996), which makes them of limited use for prognostic models. The second main approach for the choice of independent variables is that of Reference SikoniaSikonia (1982), later used, for example, by Reference VenterisVenteris (1999), which relates the calving rate to a height-above-buoyancy at the glacier terminus.

The inverse approach of Reference Van der VeenVan der Veen (1996) also uses the height-above-buoyancy criterion, but inverts the problem, focusing on the factors controlling the terminus position rather than those controlling the calving rate. His calving criterion was based on the observation that, at Columbia Glacier, Alaska, the terminus is located where the height of the terminal cliff is ∼50 m above buoyancy. Reference Vieli, Funk and BlatterVieli and others (2000, Reference Vieli, Funk and Blatter2001, Reference Vieli, Jania and Kolondra2002) adopted a modified version of the Reference Van der VeenVan der Veen (1996) criterion, in which the fixed heightabove-buoyancy is replaced by a fraction of the flotation thickness, the fraction being an adjustable model parameter. The main drawback of any of the height-above-buoyancy calving laws is that they do not allow for the development of floating ice tongues, which restricts their application to glaciers with grounded calving fronts.

3.2. Benn and others’ (2007a) calving criterion

A major advance in calving models came with the calving criterion proposed by Reference Benn, Hulton and MottramBenn and others (2007a). This criterion was preceded by a full analysis of the calving problem (Reference Benn, Warren and MottramBenn and others, 2007b), which considered the most prominent calving processes: the stretching associated with surface velocity gradients; the force imbalance at terminal ice cliffs; the undercutting by underwater melting; and the torque arising from buoyant forces. Reference Benn, Warren and MottramBenn and others (2007b) then established a hierarchy of calving mechanisms, concluding that longitudinal stretching in the large-scale velocity field of the glacier near the terminus can be considered the first-order control on calving. The other mechanisms are second-order processes, superimposed on the first-order mechanism.

Consequently, Reference Benn, Hulton and MottramBenn and others’ (2007a) criterion assumes that the calving is triggered by the downward propagation of transverse surface crevasses, near the calving front, as a result of the extensional stress regime. The crevasse depth is calculated following Reference NyeNye (1955, Reference Nye1957), assuming that the base of a field of closely spaced crevasses lies at a depth where the longitudinal tensile strain rate tending to open the crevasse equals the creep closure resulting from the ice overburden pressure. Crevasses partially or totally filled with water will penetrate deeper, because the water pressure contributes to the opening of the crevasse. These arguments lead to the following equation for crevasse depth, d:

(2)

We note that the factor 2 is located at a different position in the corresponding equations of Reference Benn, Hulton and MottramBenn and others (2007a,Reference Benn, Warren and Mottramb), which are now recognized to be in error (personal communication from D. Benn, 2009). g is the acceleration due to gravity, ρ _i and ρ _w are the densities of ice and water, A and n are Glen’s flow-law parameters, d _w is the water depth in the crevasse (Fig. 4) and is the longitudinal strain rate (∂u/∂x) minus the threshold strain rate required for crevasse initiation, . Reference Benn, Hulton and MottramBenn and others (2007a) adopted the simplifying assumption that , and hence . This could lead to a slight overestimate of the crevasse depth, as confirmed by some field measurements (e.g. Reference HoldsworthHoldsworth, 1969). In more recent work, Reference Mottram and BennMottram and Benn (2009) considered and used it as a tuning parameter to fit computed and observed crevasse depths.

Fig. 4. (a) Schematics for the model of calving by crevasse depth, adapted from figure 1 of Reference Benn, Hulton and MottramBenn and others (2007a). (b) Close-up showing some of the variables in greater detail.

Reference NyeNye’s (1957) formulation does not account for stress concentrations at the tip of the fracture. This is admissible, however, because crevasses near the terminus appear as fields of closely spaced crevasses where stress-concentration effects are reduced by the presence of nearby fractures. Formulations such as those of Reference WeertmanWeertman (1973) and Reference SmithSmith (1976), which take into account stress-concentration effects for approximating the penetration depth of isolated crevasses, do not seem appropriate for crevasses near the terminus, because they rely on the assumption of isolated crevasses. This has been confirmed by field measurements (Reference MottramMottram, 2007) which show that Reference WeertmanWeertman’s (1973) formulation consistently overestimates the depths of crevasses. There exist more rigorous frameworks for the calculation of crevasse depths, such as the linear elastic fracture mechanics (LEFM) used by Reference Van der VeenVan der Veen (1998) and Reference RistRist and others (1999). However, there are a number of difficulties associated with the use of LEFM, most notably that the assumption of linear rheology is not suitable for glacier ice. LEFM is also very sensitive to crevasse spacing, which is often unknown. Finally, field observations (Reference MottramMottram, 2007; Reference Mottram and BennMottram and Benn, 2009) have shown that Nye’s approach performs as well as the more complicated LEFM approach of Reference Van der VeenVan der Veen (1998) in predicting observed crevasse depths.

Reference Benn, Hulton and MottramBenn and others’ (2007a) criterion provides that calving occurs when the base of the crevasses reaches sea level. The terminus position is then located where the crevasse depth equals the glacier freeboard above sea level, h, i.e.

(3)

where h = H − D _w, H being the ice thickness and D _w the sea (or lake) water depth.

The justification for using sea level instead of the glacier bed as the relevant crevasse penetration for triggering calving is as follows. Observations of many calving glaciers show that surface crevasses near the terminus penetrate close to the waterline (Reference Benn and EvansBenn and Evans, 2010), and it has been noted that calving often develops as collapse of the subaerial part of the calving front followed, after some delay, by buoyant calving of the subaqueous part (e.g. Reference O’Neel, Marshall, McNamara and PfefferO‘Neel and others, 2007). An alternative explanation is that, near the terminus, longitudinal crevasses intersecting the ice front are commonly observed in addition to the transverse ones. This is especially true for fjord glaciers. If such crevasses extend below the waterline, then there is a free hydraulic connection between crevasses and sea or lake water, implying a continuous supply of water to the crevasses and associated deepening, by hydrofracturing, until the crevasse depth reaches the bed (Reference Benn and EvansBenn and Evans, 2010). Both these arguments justify, as a first approximation, taking sea level as the critical depth of crevasse penetration for triggering calving.

Equations (2) and (3) show that the first-order calving processes can be parameterized using longitudinal strain rates, ice thickness and depth of water filling the crevasse. The longitudinal strain-rate term provides the link between the calving processes and ice dynamics, and allows an easy implementation of this calving criterion in prognostic models of tidewater glacier dynamics, no matter whether the tongues are grounded or floating, or are ice shelves.

Reference AlleyAlley and others (2008) proposed a simple law for ice-shelf calving which relies on the same basic assumption as Reference Benn, Hulton and MottramBenn and others’ (2007a) model, that calving is dominated by cracks transverse to the flow. They differ in that Reference AlleyAlley and others’ (2008) calving law is empirical, and the calving rate is parameterized in terms of the product of stretching rate, ice thickness and half-width of the ice shelf.

3.3. Three-dimensional extension of, and our improvements to, Benn and others’ (2007a) model

Limitations of the published implementations of Reference Benn, Hulton and MottramBenn and others’ (2007a) calving model are that they are two-dimensional and that the dynamic models employed are simplistic (driving stress supported either by basal drag, or by drag at the lateral margins, or by a combination of both; i.e. longitudinal-stress gradients are not taken into account). To overcome such limitations, we have developed a three-dimensional extension of Reference Benn, Hulton and MottramBenn and others’ (2007a) calving criterion which uses a full-Stokes model of glacier dynamics. We refer to the application of this model to Johnsons Glacier as experiment 1. We detected a further limitation of Reference Benn, Hulton and MottramBenn and others’ (2007a) model, related to the use of Reference NyeNye’s (1955, Reference Nye1957) formula for calculating crevasse depths. Consequently, we developed a further improvement to Reference Benn, Hulton and MottramBenn and others’ (2007a) model, which we refer to as experiment 2. Experiment 3, which further departs from Reference Benn, Hulton and MottramBenn and others’ (2007a) model, consists of an additional improvement in the calculation of crevasse depth. Finally, the introduction of a ‘yield strain rate’ is considered in experiment 4. All of these successively refined calving models share the basic features described below, and their implementation steps are very similar. We first present their common aspects and then give the details for each particular experiment. In the following, we use d ₀ to denote the crevasse depth when the crevasse does not have any filling water. In the case of crevasse depths calculated using Reference NyeNye’s (1955, Reference Nye1957) formula, we therefore have

(4)

where B is the stiffness parameter related to the softness parameter used in Equation (2), A, by B = A ^−1/n . The common features of all of our experiments, corresponding to their three-dimensional model geometry, are:

1. 1. We consider the glacier length, L, as a function of x and y, i.e. L = L(x, y), which represents the glacier length following each flowline.

2. 2. The longitudinal strain rate, , calculated as ∂u/∂x in Reference Benn, Hulton and MottramBenn and others’ (2007a) model, now becomes the strain rate following the ice-flow direction at each point, calculated as

(5)

1. where u₁ = (u ₁, v ₁) and u₂ = (u ₂, v ₂) are the velocity vectors at two positions, x₁ = (x ₁, y ₁) and x₂ = (x ₂, y ₂), with x₂ = x₁ + u₁ Δt. We have neglected the vertical component of velocity, which is very small near the calving front. Equation (5) gives an approximation to the rate of change of the velocity field, u, along the direction of flow in which changes in both magnitude and direction of u along the flowline are taken into account. A more accurate estimate of the strain rate in the direction of flow could be obtained by computing the norm of the directional derivative of the vector field, u, in the direction of flow. The latter, however, involves a rather cumbersome computation which gives no appreciable differences in the calculated crevasse depths. The procedure described is approximately equivalent to rotating the traction vector to align it in the direction of ice flow.

2. 3. is computed from the velocity field resulting from the solution of a full-Stokes model of ice dynamics described in section 4.

The steps for implementing all the models considered can be summarized as follows:

1. 1. Let L ₀(x, y) be the terminus position at a starting time, t ₀.

2. 2. We solve the dynamical full-Stokes model to obtain the velocity field, u(x, y). The velocity field under consideration (whether at the surface or at different depths) depends on the particular model employed.

3. 3. We compute from u(x, y) using Equation (5).

4. 4. Using , and perhaps a given (or modelled) height of water partially filling the crevasse, we calculate the crevasse depth, d(x, y). The method used for calculating the crevasse depth depends on the particular model under consideration.

5. 5. The new terminus position, L ₁(x, y), will be located where d − h = 0.

6. 6. We use the computed velocity field at the surface, u_s(x, y), and a given accumulation/ablation field, a(x, y), to calculate the new glacier geometry after a time, Δt, calving the glacier front, as determined by L ₁(x, y). Δt can be assigned any value, provided a(x, y) is available at such a temporal resolution. The choice of Δt determines whether or not the model reproduces seasonal variations. This is important if one wishes to consider the contribution of surface meltwater to crevasse deepening.

7. 7. We repeat the above processes until we reach a steady-state configuration (surface geometry and front position) consistent with the prescribed a(x, y).

Experiment 2

The model considered in this experiment involves an important conceptual difference to Reference Benn, Hulton and MottramBenn and others’ (2007a) calving model. The derivation of Reference NyeNye’s (1957) formula for calculating crevasse depth assumes plane strain, i.e. v = 0 (v being the velocity component transverse to the direction of flow) and τ_xy = τ_yz = 0 (deviatoric shear-stress components acting on the plane of flow), from which the normal stress component, σ_yy , also equals zero. Reference NyeNye’s (1957) analysis also assumes ∂w/∂x = 0. As Nye points out, ‘The main feature of a real glacier that is omitted in our model is the variation in the y [z in Nye’s terminology] direction – that is to say, the influence of the valley walls’. Consequently, the calculation of crevasse depth using Reference NyeNye’s (1957) equation is an approximation which becomes better for wider glaciers, and probably very good for an ice sheet. Reference Benn, Hulton and MottramBenn and others (2007a) rightly stress that their analysis focuses on flow along the centre line of an outlet glacier, because at exactly the centre line (but not if we move away from it) the same assumptions of plane strain hold. Therefore, the use of Reference NyeNye’s (1957) formula for calculating crevasse depths is valid only for crevasses near the centre line, and is expected to worsen as we move away from it.

To reduce this problem, instead of using Nye’s equation as given by Equation (4), we have derived the crevasse depth equation by balancing two terms: (1) the tensile deviatoric stress tending to open the crevasse, calculated directly from Reference NyeNye’s (1957) generalization of Reference GlennGlen’s (1955) constitutive equation and the strain field produced as output of our full-Stokes model, and (2) the ice overburden pressure tending to close the crevasse, approximated as ρ _i gd ₀. Reference NyeNye’s (1957) generalization of Reference GlennGlen’s (1955) flow law is given by

(6)

where τ_ij is the deviatoric stress tensor, and

(7)

are the strain-rate tensor and the effective strain rate, respectively. If, in order to allow direct comparison with Equation (4), we consider that the ice flows in the direction of x, we find

(8)

with . The crevasse depth is thus given by

(9)

which can be compared with Equation (4) with . We note the absence in Equation (9) of the factor 2 and the presence of the effective strain-rate factor. In the general case of ice flow in any direction, we use

(10)

with the strain rate following the ice-flow direction at each point, calculated using Equation (5). In this experiment, as in experiment 1, we use the velocity field at the surface, u_s(x, y), to compute the along-flow strain rate.

4.1. Basic equations

The equations describing the dynamical model are the usual ones for steady conservation of linear momentum and conservation of mass for an incompressible continuous medium:

(11)

(12)

where σ_ij and u_i represent the stress tensor and velocity vector components, respectively, and ρ is ice density and g_i are the components of the vector accelaration of gravity. We follow Einstein’s convention of summation on repeated indexes. As a constitutive relation, we adopt Glen’s flow law, given by Equation (6). The constitutive relation is expressed in terms of deviatoric stresses, while the conservation of linear momentum is given in terms of full stresses. Both stresses are linked through the equation

(13)

where p is the pressure (compressive mean stress). Conservation of angular momentum implies the symmetry of the stress tensor, i.e. σ_ij = σ_ji . We have set n = 3 and taken B as a free parameter to be tuned to fit the observed and computed velocities at the glacier surface.

4.2. Boundary conditions

The domain boundary can be divided into portions that have boundary conditions of different types. The upper surface is a traction-free area with velocities unconstrained. At the ice divides the horizontal velocities and the shear stresses are expected to be small, so we set null horizontal velocities and shear stresses and leave the vertical velocity unconstrained. Note that, in a three-dimensional model, the horizontal velocities at the divides are not necessarily zero. What should be zero are their components normal to the divide plane, but the components tangent to the latter could be non-zero, though they will usually be small. At the basal nodes the horizontal velocities are specified according to a Weertman-type sliding law (e.g. Reference PatersonPaterson, 1994, chapter 7), while the vertical component is derived from the horizontal ones and the mass-continuity condition, assuming a rigid bed:

(14a)

(14b)

(14c)

where u, v and w are the components of velocity and the subscript ‘b’ denotes evaluation at the glacier bed. The variables s and b represent the vertical coordinates of the glacier surface and bed, respectively, and H _w is the height of the basal water column. As a simplifying approximation, we computed the basal shear stress using the shallow-ice approximation (SIA) and calculated the ice pressure at the bed hydrostatically. This is a rough approximation for a model that solves the full-Stokes system of equations in the interior of the domain. It could be improved by introducing an iteration loop for the computation of the basal shear stress. The SIA estimate would be the initial iteration, successively refined by the basal shear stresses computed using the full-Stokes model. However, this would be extremely costly computationally. Other workers using finite-element three-dimensional models of glacier dynamics have used the same approach (e.g. Reference HansonHanson, 1995). In the above equations, we used p = 2 and q = 1, though different values were tested during the tuning procedure.

As there are no basal water-pressure measurements available for Johnsons Glacier, we derived a functional form for the height of the basal water column. At the glacier front, we assume that it equals the sea-water depth (Reference Vieli, Funk and BlatterVieli and others, 2000). For the interior of the glacier, we use different laws for the ablation and accumulation zones. We parameterize them in terms of the distances, measured along flowlines, between the point and the ice front (d _front,point)_, the point and the equilibrium line (d _ela,point)_, the front and the equilibrium line (d _front,ela) and the equilibrium line and the summit (d _ela,summit) (Fig. 2):

(15)

(16)

where the superscripts ‘ab’ and ‘ac’ denote ablation and accumulation zones. The constant, C, defining the slope of the straight line was roughly fitted through a pre-tuning procedure for which we assigned to the constants B in the flow law and K in the sliding law values from a previous paper on Johnsons Glacier dynamics that did not include calving (Reference Martín, Navarro, Otero, Cuadrado and CorcueraMartín and others, 2004). Different trials for C led to the choice C = 40. The height of the basal water column given by Equations (15) and (16) is plotted in Figure 5. The assumption that, at the glacier front, it equals the sea-water depth is a reasonable one. The assumption that it smoothly increases as we move up-glacier through the ablation area is also quite reasonable, and similar linear variations (as a simple choice for the functional form) have been used by others (e.g. Reference Vieli, Funk and BlatterVieli and others, 2000). Finally, the assumption that it decreases as we move from the equilibrium-line altitude (ELA) to the glacier head has also been used by other workers (e.g. Reference HansonHanson, 1995) and is consistent with the equipotential surfaces of the hydraulic potential dipping up-glacier (with a slope of ∼11 times the slope of the glacier surface), resulting in englacial water conduits dipping down-glacier, approximately perpendicular to the equipotential surfaces (Reference ShreveShreve, 1972).

Fig. 5. Height of the basal water-saturated ice (dark blue), as given by Equations (15) and (16), along the flowline of Johnsons Glacier indicated by a red line in Figure 2. Light blue indicates unsaturated glacier ice above the basal water-saturated ice.

At the glacier front, we have set stress boundary conditions, taking them as null above sea level and equal to the hydrostatic pressure below sea level.

4.3. Numerical solution

Equations (11) and (12) are a Stokes system of equations, which is solved without discarding any stress components, i.e. using a full-Stokes solution. It is non-linear because of the constitutive relation (Equation (6)) employed. The unknowns are the velocity components and the pressure. This system is reformulated in a weak form, whose solution is approximated by finite-element methods with Galerkin formulation, using a mixed velocity/pressure scheme (e.g. Reference Quarteroni and ValliQuarteroni and Valli, 1994). The basis functions of the approximating spaces for velocity and pressure are taken as quadratic and linear, respectively. This choice of polynomial spaces of different degree is dictated by considerations of convergence and stability of the numerical solution (e.g. Reference Carey and OdenCarey and Oden, 1986). The ice density and the rheological parameters, B and n, as well as the density, have been taken as constant across the entire glacier. The non-linear system of equations resulting from the finite-element discretization was solved iteratively, using a direct procedure based on fixed-point iteration (Reference Martín, Navarro, Otero, Cuadrado and CorcueraMartín and others, 2004). The linear system associated with each step of this iterative procedure was solved by lower–upper decomposition. The iteration procedure stops when the norm of the vector difference between the solutions for successive iterations falls below a prescribed tolerance. The system actually solved is a non-dimensional one. The details of the dimensionless variables are given by Reference Corcuera, Navarro, Martín, Calvet and XimenisCorcuera and others (2001). Reference Martín, Navarro, Otero, Cuadrado and CorcueraMartín and others (2004) give further details of the numerical-solution procedure, including the full set of finite-element equations.

Experiment 1

From the velocity field, we compute the surface strain rate using Equation (5) and then the depth of the crevasses using Equation (4). As we have no a priori knowledge of the depth of water filling the crevasses, d _w, we set it to zero and therefore use Equation (4) instead of Equation (2). The values of d ₀(x, y) − h(x, y) are plotted in Figure 10a. The contour line d ₀ − h = 0 gives the terminus location in our three-dimensional extension of Reference Benn, Hulton and MottramBenn and others’ (2007a) calving model. According to Figure 10a, Johnsons Glacier front is located at a position behind (up-glacier from) that predicted by the calving criterion. This could be due to our assumption of crevasses empty of water, i.e. d _w = 0. To establish whether this is a reasonable hypothesis, we recomputed d(x, y) − h(x, y) for different amounts of water filling the crevasses, using Equation (2). The results for the case in which d _w equals, at each point, one-half of the ice thickness, H, are shown in Figure 10b. We observe that d − h = 0 for most of the ice front. In other words, we would need a height of water d _w = H/2 filling the crevasses near the terminus to properly calculate, using the three-dimensional extension of Reference Benn, Hulton and MottramBenn and others’ (2007a) calving criterion, the location of the presently observed calving front.

Fig. 10. (a) Plot of contour lines for d ₀(x, y) − h(x, y) in experiment 1; contour level interval is 20 m. The zero line should correspond to the location of the calving front, but it does not. (b) Contour lines for d(x, y) − h(x, y) in experiment 1 with a depth of water filling the crevasses equal to half the ice thickness at each point. In this case, the zero contour line gives a good fit to the position of the calving front. (c) Contour lines for d ₀(x, y) − h(x, y) in experiment 2. Contour level interval is 20 m. (d) Plot of contour lines for d ₀(x, y) − h(x, y) in experiment 3. (e) Contour lines for d(x, y) − h(x, y) in experiment 3, with a depth of water filling the crevasses equal to one-tenth of the ice thickness at each point. (f) Contour lines for d ₀(x, y) − h(x, y) in experiment 4. (g) Plot of contour lines for d(x, y) − h(x, y) in experiment 4, with a depth of water filling the crevasses equal to one-sixth of the ice thickness at each point. Contour line interval is 10 m unless otherwise stated. (h) Location of the area considered.

Experiment 2

In this case we compute the crevasse depth using Equation (10) instead of Equation (4). The corresponding contour lines of d ₀(x, y) − h(x, y) are shown in Figure 10c. The differences between the results of experiments 1 and 2 are quite remarkable. In the latter case, there is no need for any water filling the crevasses. In fact, the model solution shows that the d ₀ − h = 0 contour line is located within the glacier domain, meaning that the glacier should be calved along this line. In other words, the crevasse depths are slightly overestimated.

Experiment 3

In this experiment we consider the model-computed tensile deviatoric stress opening the crevasse as a function of depth, and determine the depth at which it is balanced by the ice overburden pressure tending to close the crevasse. The contour lines for d ₀(x, y) − h(x, y) are shown in Figure 10d. The −10 m contour line lies very close to the calving front, indicating that an additional 10 m of crevasse depth near the terminus would be required to have the glacier calving at the presently observed front position. If we fill the frontal crevasses with a height of water, d _w, equal to one-tenth of the ice thickness, the contribution of the water pressure, ρ_wgd_w , to the deepening of the crevasse makes the contour line d − h = 0 lie very close to the observed front position, as shown in Figure 10e. In other words, this experiment, with such an amount of water, accurately models the presently observed calving front position.

Experiment 4

This experiment examines the influence of yield strain rate, , in the computation of the crevasse depth and, consequently, in the calculation of the position of the calving front. We used , which is the yield strain rate corresponding to the yield stress of 60 kPa that Reference Mottram and BennMottram and Benn (2009) found produced the best fit in their experiment. We introduced such a yield strain rate in the most physically plausible of the above experiments, i.e. experiment 3, but with crevasses empty of water. The results are shown in Figure 10f. For experiment 3, we found the contour line d ₀ − h = −10 m lay very close to the observed calving front position. Adding a yield stress has the effect of reducing the crevasse depth and, correspondingly, moving a given d ₀ − h contour line farther down-glacier, so we could not match the observed front position by varying the value of (any value of , no matter how small, worsens the fit). However, the idea of a critical-stress intensity factor necessary to overcome the fracture toughness of the ice is physically plausible, so we kept the and searched for the amount of water filling the crevasses needed to move the d − h = 0 contour line as close as possible to the calving front. The required height of water, d _w, was one-sixth of the ice thickness, as shown in Figure 10g.

7.1. Limitations of Benn and others’ (2007a) calving criterion and its implementation

1. 1. The most important limitation of Reference Benn, Hulton and MottramBenn and others’ (2007a) calving criterion is the plane strain assumption inherent in the use of Reference NyeNye’s (1957) formula for calculating crevasse depths. (We discussed this when we introduced experiment 2 in section 3.3.) The key issue is that the accuracy of the approximation is better for wider glaciers, and probably will be best for ice sheets. For a given glacier, the estimate will be best near the centre line. Glaciers showing narrowing or widening of their tongues (especially if these are sharp) will not be good candidates for such a crevasse-depth model, because of the expected influence of lateral (transverse to flow) compression or extension, unless we strictly focus on the centre line. This is probably the reason why the crevasse-depth model behaved so poorly for Johnsons Glacier. This glacier fills a basin that drains through a narrow calving front, implying fast flow acceleration and transverse-to-flow compression in the vicinity of the calving front.

2. 2. Our second criticism does not refer to the calving criterion itself but to the simplistic models of glacier dynamics employed by Reference Benn, Hulton and MottramBenn and others (2007a). They used three different models of glacier dynamics: (1) driving stress entirely balanced by basal drag; (2) driving stress supported by drag at the lateral margins; and (3) resistance to flow provided by a combination of basal and lateral drag. Therefore, longitudinal-stress gradients were not taken into account in the force balance. This oversimplification strongly contrasts with the use of a calving criterion based on the penetration depth of crevasses developed near the calving front as a result of the extensional-stress regime, because one might expect longitudinal-stress gradients. However, such models, though simplistic, are not inconsistent with the calving criterion because the crevasses open in response to longitudinal stresses, not longitudinal-stress gradients. In the original model formulation, Reference Benn, Hulton and MottramBenn and others (2007a) assume that the driving stress is everywhere equal to the sum of the basal and lateral drag. In such a situation, one can have longitudinal stresses (if the driving stress changes down-glacier) without longitudinal-stress gradients. Nevertheless, the longitudinal-stress gradients are important near the terminus of calving glaciers, and their inclusion is an important feature of the present modelling study.

7.2. The role of other calving mechanisms

Reference Benn, Warren and MottramBenn and others’ (2007b) analysis of the main mechanisms responsible for calving concludes that the stretching associated with surface velocity gradients near the terminus is the first-order mechanism for triggering calving. Other mechanisms, such as the force imbalance at terminal ice cliffs, the undercutting by underwater melting and the torque arising from buoyant forces (which may be oscillatory, as happens in the case of the variable torque due to ocean tides), are second-order processes, superimposed on the first-order mechanism. These contributions need to be accounted for to get a better representation of calving, though the necessary field data are rarely available. The relative importance of each mechanism depends on the particular case study. For Johnsons Glacier, considering the minute portion of the glacier with the bed below sea level, the torque arising from buoyant forces is certainly small. Because of the closed embayment, waves hitting the glacier front are rather weak, so wave erosion at sea level is small. There is, however, some melting at the underwater part, attributed to warm waters in the embayment due to the very shallow water depth (a few metres). Finally, we note that the force imbalance at the terminal ice cliff is accounted for in our model, because the stress boundary condition set at the terminal cliff (null stresses in the emerged part, and stresses equal to hydrostatic pressure in the submerged part) are, by themselves, the expression of such force imbalance.

7.3. Difficulties and weaknesses of our modelling

During the tuning of free parameters of the model, we had difficulties in closely fitting the relatively large velocities observed at the glacier surface near the terminus. The very small surface slopes in this area imply, by virtue of the Weertman-type sliding law, small sliding velocities. We failed to substantially increase the sliding velocity by decreasing the value of the effective pressure term in the sliding law. Moreover, the small velocity gradients near the terminus also implied small crevasse depths, with the consequence that they never reached sea level, preventing calving from occurring. We finally obtained good agreement with observed velocities, as well as deeper crevasse depths, by setting different values of the sliding parameter, K, for the ablation and accumulation zones.

An evident weakness of the full-Stokes model employed in this paper is that it uses, for reasons of computational cost, a basal boundary condition in terms of velocities that make use of basal shear stresses computed by balancing the driving stress. Also our choice of the parameterization of height of the basal water pressure, contributing to sliding, involved a certain amount of speculation, because of the lack of basal water-pressure measurements on Johnsons Glacier. All of these arguments diminish the accuracy of the sliding model employed.

7.4. Comments on the results

The results of experiment 1 show that our straightforward three-dimensional extension of Reference Benn, Hulton and MottramBenn and others’ (2007a) calving criterion requires a large amount of water to fill the crevasses near the calving front (water height about half the ice thickness) to properly reproduce the presently observed location of Johnsons Glacier calving front. For this particular glacier, there are no field observations to verify whether the assumed amount of water filling the crevasses is realistic.

Experiment 2 had markedly better results than experiment 1. No water filling the crevasses was required to fit the observed front position. We attribute this improvement to the fact that the effective strain rate, with all of its components, comes into play for computing the crevasse depth. Nevertheless, this model slightly overestimated the crevasse depth necessary for triggering calving.

The model involved in experiment 3 has a yet more solid physical foundation than that of experiment 2. In experiment 3, a very small amount of water (height of water filling the crevasse about one-tenth of the ice thickness at the point under consideration) was required to provide a good fit to the observed calving front position. We believe that a model using strain rate and stresses dependent on depth provides a better approximation to real glaciers.

Introducing the effect of the yield stress rate, , (experiment 4) did not improve the fit to the observed front position of Johnsons Glacier. In spite of this, including in the modelling should not be discarded, because of its interpretation as a heuristic device to fulfil the role of a critical stress intensity factor required to overcome the fracture toughness of the ice.

In experiments 1, 3 and 4 we used a certain amount of water filling the near-front crevasses in order to accurately reproduce the observed front position of Johnsons Glacier calving front. We thus used the water filling the crevasses as a free parameter for tuning our model to fit the observations. This is true of the particular application in this paper, but it is not inherent to the use of Reference Benn, Hulton and MottramBenn and others’ (2007a) calving criterion, i.e. water filling the crevasses is not necessarily used as a tuning parameter in Reference Benn, Hulton and MottramBenn and others’ (2007a) calving criterion. It could, instead, be an input parameter supplied by a model of melting at the glacier surface.

Our three-dimensional extension of Reference Benn, Hulton and MottramBenn and others’ (2007a) calving criterion, with Reference NyeNye’s (1957) formula (experiment 1), failed to adequately model the crevasse depths near the calving front of Johnsons Glacier, and thus to reproduce its current position. This, however, brings no discredit to either criterion or formula. Simply, the basic assumption of plane strain did not apply in our case study. Although the models considered in experiments 2, 3 and 4 represented clear improvements over those of experiment 1, all of them rely on the same fundamental assumptions of Reference Benn, Hulton and MottramBenn and others’ (2007a) calving criterion: (1) that calving is triggered by the downward propagation of transverse crevasses near the calving front resulting from the extensional stress regime and (2) that the crevasse depth is calculated by balancing the tensile deviatoric stress opening the crevasse with the ice overburden pressure tending to close it. The differences between the models presented here are more in the details of the particular methods than in the underlying physical ideas.

An interesting and appealing feature of the modelling results shown here is that they produce a ‘calving bay’ at the glacier front. Intuitively, this appears to be the result of pure extensional strain near the glacier centre line and simple shear at the glacier margins, with associated rotation of the principal strain axes. This provides an additional source of model validation.

7.5. The role of water as a mechanism favouring calving

Reference Benn, Hulton and MottramBenn and others’ (2007a) calving criterion allows the effect on crevasse deepening of water filling the crevasses to be included. To our knowledge, there are few reported observations of water-filled crevasses near the front of tidewater glaciers (e.g. Reference Mottram and BennMottram and Benn, 2009). In contrast, there are well-known examples of ice shelves where water filling the crevasses was indeed a contributing factor to the crevasse deepening and subsequent calving. The collapse of the northern part of Larsen B ice shelf, Antarctica, in 2002 is a clear example; in this case, however, other factors also contributed to its disintegration, such as ice-shelf thinning by bottom melting due to warming ocean waters (Reference Shepherd, Wingham, Payne and SkvarcaShepherd and others, 2003).

Even if there is no water filling the crevasses, we hypothesize that the effect of surface water on triggering calving can still be important. During the melt season, it is clear that surface meltwater drains into the crevasses. If there is a good hydraulic connection between the crevasse and the glacier bed, such water will not remain in the crevasse but will penetrate to the glacier bed. However, on its way to the bed this water can contribute to crevasse widening and deepening by melting the crevasse walls and by enlarging the conduits to the bed. No matter whether this flow is mainly accomplished through conduits or fractures, it will certainly contribute to the weakening of the material as a whole and, consequently, favour the occurrence of calving. Of course, the state of the englacial drainage system will play a role in the effectiveness of such mechanisms. For instance, we can think of situations where a very active drainage system allows water to drain to the bed, and eventually the ocean, from very shallow depths within the ice, without supplying any important amount of water to the near-front crevasses. There are also regions, such as many Antarctic ice shelves, where there is no evidence of a free hydraulic connection between fractures and the ocean. In fact, hot-water drilling often indicates just the opposite, as a connection with the ocean is not made until the drill penetrates near the ice-shelf bed or into the marine ice layer.

In light of the above arguments, investigation of the contribution of water to the triggering of calving is a research subject that deserves further attention.

7.6. Polythermal structure of Johnsons Glacier

The value of the stiffness parameter that best fitted the computed and observed velocities was B = 0.23 Mpa a^1/3. This is close to the value obtained for the neighbouring (but ending on land) Hurd Glacier (0.22 Mpa a^1/3; Reference OteroOtero, 2008), for which ground-penetrating radar measurements and geomorphological evidence suggest a polythermal structure (Reference NavarroNavarro and others, 2009). It is also quite similar to, though slightly higher than, the value of 0.20 Mpa a^1/3 for Storglaciären, Sweden (Reference HansonHanson, 1995), which is also known to be polythermal. The higher stiffness of Johnsons Glacier could be interpreted as indicative of a larger relative amount of cold ice. To give a rough idea of the degree of variation of B with temperature, two examples extracted from Reference PatersonPaterson (1994) are: B = 0.167 Mpa a^1/3 at 0°C; and B = 0.236 Mpa a^1/3 at −2°C. This subject is of interest because the South Shetland Islands glaciers have commonly been assumed to be temperate (e.g. Reference Qin, Zielinski, Germani, Ren, Wang and WangQin and others, 1994; Reference Furdada, Pourchet and VilaplanaFurdada and others, 1999; Reference Benjumea, Macheret, Navarro and TeixidóBenjumea and others, 2003), while more recent work (Reference Molina, Navarro, Calver, García-Sellés and LapazaranMolina and others, 2007; Reference NavarroNavarro and others, 2009) has suggested that at least some Hurd Peninsula glaciers ending on land are polythermal. Our results suggest that the tidewater Johnsons Glacier could also be polythermal. In addition to its dynamic implications, the polythermal structure has an impact on the glacier sensitivity to climate changes. Temperate glaciers are especially sensitive to climate changes, as almost all heat supplied is used for melting. If, however, the glaciers are polythermal, heat is first used to raise the temperature to the melting point. Comparing the specific heat capacity with the latent heat of fusion (e.g. Reference PatersonPaterson, 1994, p.205, table 10.1) shows that this effect is only marginal, because the energy needed to melt a given amount of ice at the pressure-melting point is ∼150 times the energy needed to raise the temperature of the same amount of ice by 1°C. However, other effects, such as the insulating effect of the cold ice layer, should be taken into account. As discussed by Reference Blatter and HutterBlatter and Hutter (1991), the temperature gradient through the cold ice layer is a dominant control on the capability to transport away from the cold/temperate transition surface the latent heat released when the top of the temperate layer freezes, which has important implications on the hydrothermal regime of polythermal glaciers.

In any case, the cold layer of these glaciers would have temperatures slightly below the melting point, so the assumption of isothermal ice mass made in our modelling is reasonable.