A three-dimensional calving model: numerical experiments on Johnsons Glacier, Livingston Island, Antarctica

Abstract

Calving from tidewater glaciers and ice shelves accounts for around half the mass loss from both polar ice sheets, yet the process is not well represented in prognostic models of ice dynamics. Benn and others proposed a calving criterion appropriate for both grounded and floating glacier tongues or ice shelves, based on the penetration depth of transverse crevasses near the calving front, computed using Nye’s formula. The criterion is readily incorporated into glacier and ice-sheet models, but has not been fully validated with observations. We apply a three-dimensional extension of Benn and others’ criterion, incorporated into a full-Stokes model of glacier dynamics, to estimate the current position of the calving front of Johnsons Glacier, Antarctica. We find that two improvements to the original model are necessary to accurately reproduce the observed calving front: (1) computation of the tensile deviatoric stress opening the crevasse using the full-stress solution and (2) consideration of such a tensile stress as a function of depth. Our modelling results also suggest that Johnsons Glacier has a polythermal structure, rather than the temperate structure suggested by earlier studies.

1. Introduction

Iceberg calving is an important mass-loss mechanism from ice shelves and tidewater glaciers for many mid- and high-latitude glaciers and ice caps (e.g. Dowdeswell and others, 2008) and for the polar ice sheets. It accounts for about half the losses from the Greenland ice sheet (e.g. Thomas, 2004) and it has long been believed to be the dominant mechanism of ice loss from the Antarctic ice sheet (e.g. Bentley, 2004). Recent work, however, has suggested that the ice loss from Antarctica is ∼50% due to calving and ∼50% due to basal melting, although there is a wide variability (10–90%) among different drainage basins (personal communication from E. Rignot, 2009). Despite its important contribution to the mass budget of glaciers and ice sheets, an adequate representation of calving is still missing from prognostic models of ice dynamics. As Benn and others’ (2007b) review of the calving problem pointed out, understanding the processes controlling calving rates has long been a major unsolved problem in glaciology. The calving process is further complicated by its close link to two other longstanding problems of glaciology: the realistic computation of stresses near the grounding line (e.g. Hindmarsh, 2006; Schoof, 2007a,b) and the appropriate representation of basal sliding (e.g. Fowler, 1981, 1986; Schoof, 2005, in press; Gagliardini and others, 2007). Recent catastrophic ice-shelf break-up (MacAyeal and others, 2003; Shepherd and others, 2003) followed by acceleration of outlet glaciers that feed them (Rott and others, 1996; Rignot and others, 2004; Scambos and others, 2004), and acceleration of Greenland outlet glaciers (Howat and others, 2007, 2008; Holland and others, 2008; Joughin and others, 2008) have increased the interest in obtaining a proper understanding of calving processes. Recent work has further stressed the need for such an understanding; for example, Nick and others (2009) showed that the recent ice acceleration, thinning and retreat of Greenland’s large Helheim Glacier began at the calving terminus, then propagated extremely rapidly upstream through dynamic coupling along the glacier. This suggests that such changes are unlikely to be caused by basal lubrication through surface melt propagating to the glacier bed, a flow acceleration mechanism proposed for the Greenland ice sheet by Zwally and others (2002) and supported, amongst others, by Parizek and Alley (2004).

Many so-called ‘calving laws’ have been proposed to incorporate the calving processes into prognostic models of glacier and ice-sheet dynamics (see section 3 for a discussion of the main ones). Most of these laws are empirical (i.e. observation-based) rather than theoretical (i.e. physically based). Moreover, some are valid only for particular types of glaciers (e.g. tidewater glaciers with grounded calving front). These characteristics make many of them inappropriate for use in prognostic models. The calving criterion recently proposed by Benn and others (2007a) is based on the penetration depth of transverse crevasses near the calving front, computed from the strain field using Nye’s (1955, 1957) formula; i.e. it is theoretical. A great advantage of this criterion is flexibility: it can be applied to both floating and grounded tidewater glaciers, and to ice shelves. However, the models of ice dynamics used by Benn and others (2007a) are limited in that they are two-dimensional, do not incorporate longitudinal-stress gradients and are strictly valid only for crevasses near the centre line of glaciers. In this paper, we present a three-dimensional extension of Benn and others’ (2007a) calving criterion, which uses a full-Stokes model of glacier dynamics, thus providing a more realistic representation of the velocity gradient producing the transverse crevasses. Additionally, we have computed the crevasse depths using not only Nye’s (1955, 1957) simplified formula but also more advanced methods based on the full-Stokes stress field derived from the model solution. We have applied the improved model to Johnsons Glacier, a grounded tidewater glacier on Livingston Island, Antarctica. The record of the front positions of Johnsons Glacier spans only a few years during the last decade, and during this observation period the front has remained at a nearly constant position, so a full modelling exercise of time evolution to follow the front-position changes of the glacier has not been possible. Instead, our modelling experiment is a diagnostic one, aimed at establishing whether the model adequately reproduces the current front position of Johnsons Glacier. Our results validate fundamental assumptions made in Benn and others’ (2007a) calving criterion, although the above-mentioned improvements to the model were necessary to accurately reproduce the current position of the calving front.

2. Geographical Setting, Glacier Geometry and Glaciological Data

Johnsons Glacier is a small (∼5.6 km²) tidewater glacier on Livingston Island, South Shetland Islands, Antarctica (Fig. 1), that terminates in a 50 m high ice cliff extending 500 m along the coast. A local ice divide, with altitudes between 200 and 330 m a.s.l., defines it as a separate glacier basin within the Hurd Peninsula ice cap (Fig. 2). The northern part of the glacier has steeper slopes (typical values ∼10°) than those in the southern part (typical values ∼6°). The confluence of the northern and southern flows of ice results in a folded and highly fractured terminal zone (Ximenis and others, 2000).

Fig. 1. Location and map of Livingston Island, South Shetland Islands. The right panel shows details of Hurd Peninsula, where Johnsons and Hurd Glaciers are located.

Fig. 2. Surface topography of Johnsons Glacier derived from geodetic measurements (total station and differential GPS) in 1999–2000; contour level interval is 20 m. The black dots indicate stakes for ice-velocity and mass-balance measurements. The red curve indicates an example flowline on which the parameters involved in the parameterization of the height of the basal water column are shown.

Ice surface velocities of Johnsons Glacier reach values up to 44 m a⁻¹ near the calving front (Ximenis, 2001). Accumulation and ablation rates show a large spatial and temporal (yearly) variability, with maximum accumulation rates ∼1 m w.e. a⁻¹ (reached in the northern sector because of the topography and the prevailing northeasterly wind direction) and maximum ablation rates up to −4 m w.e. a⁻¹ measured over the past 10 years (Ximenis, 2001; Otero, 2008). The equilibrium line (approximate location shown in Fig. 2) lies close to the 180 m altitude contour line.

Johnsons Glacier has traditionally been considered a temperate glacier, on the basis of limited temperature–depth profiles measured at some shallow/intermediate boreholes (Furdada and others, 1999), which have been said to be consistent with measurements at other locations in the South Shetland Islands (e.g. Qin and others, 1994), and also on the basis of measured radio-wave velocities in ice (Benjumea and others, 2003). However, more recent investigations (Molina and others, 2007; Navarro and others, 2009) have suggested that at least the land-terminating glaciers on Hurd Peninsula are polythermal. Ground-penetrating radar sections of Johnsons Glacier also show, in the ablation zone, some layers and patches of cold ice.

The surface topography of Johnsons Glacier, determined from geodetic measurements (total station and differential GPS; Molina and others, 2007), is shown in Figure 2, while Figure 3 shows the bedrock topography retrieved from low-frequency (20 MHz) radio-echo sounding measurements (Benjumea and others, 2003; Navarro and others, 2005, 2009). Figure 3 also shows the location of the radar profiles and, in the proglacial embayment, the bathymetric profiles. The radar data show a maximum ice thickness of ∼160 ± 3 m and an average thickness of 97 ± 3 m. Total ice volume in 2000 was estimated at 0.545 ± 0.014 km³ (Molina and others, 2007). Note, in Figure 3, that most of the glacier bed is above sea level, reaching values slightly below sea level only next to the calving front and in the central part of the basin, further up-glacier from the terminus.

Fig. 3. Bedrock topography of Johnsons Glacier determined by subtracting the ice thickness (retrieved from low-frequency (20 MHz) radio-echo sounding measurements) from the surface topography (Fig. 2); contour line interval is 20 m. The red lines on the glacier surfaces indicate the radar profiles, and the curves in the proglacial embayment indicate the bathymetric profiles. The ice thickness for the highly crevassed terminal area (down-glacier from the dashed curve in the figure) was determined by interpolation between the glacier-bed topography up-glacier from the dashed curve and the sea-bed topography in the neighbourhood of the terminal cliff.

Although Johnsons Glacier has been losing mass for at least the past 50 years, the geodetic mass balance during the period 1956–2000 has been moderately negative at −0.23 m w.e. a⁻¹ for the ensemble Johnsons–Hurd Glaciers (Molina and others, 2007), and the Johnsons Glacier calving front has remained at a nearly constant position during the past decade (D. García and J. Calvet, unpublished data). The latter is consistent with the small water depth in the proglacial embayment (just a few metres) and the absence of a reverse-slope bed near the present calving front (Meier and Post, 1987; Van der Veen, 1996; Vieli and others, 2001). However, it is difficult to quantify, a priori, how much the present surface geometry differs from a steady-state one, because this depends on the availability of an adequate representation of the calving processes.

We have previously published two modelling experiments for Johnsons Glacier. One uses a two-dimensional model (Corcuera and others, 2001), the other uses a three-dimensional model (Martín and others, 2004); however, neither of these models includes a calving front. Instead, we excluded the glacier terminus region from the model domain and defined an artificial boundary in the lower part of the glacier. There we established velocity boundary conditions based on measured ice velocities at the Johnsons Glacier network of stakes, which is dense near the terminus. In this paper, we extend our domain to include the calving front.

3. Calving Laws

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 (Van 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. Siegert 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. Van der Veen, 1996, 2002). 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 Benn 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 (Brown and others, 1982; Pelto 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 (Haresign, 2004), are quite different for tidewater and freshwater glaciers (Funk and Röthlisberger, 1989) and also vary with time (Van 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 Sikonia (1982), later used, for example, by Venteris (1999), which relates the calving rate to a height-above-buoyancy at the glacier terminus.

The inverse approach of Van 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. Vieli and others (2000, 2001, 2002) adopted a modified version of the Van 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 Benn and others (2007a). This criterion was preceded by a full analysis of the calving problem (Benn 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. Benn 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, Benn 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 Nye (1955, 1957), 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 Benn and others (2007a,b), 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, . Benn 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. Holdsworth, 1969). In more recent work, Mottram 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 Benn and others (2007a). (b) Close-up showing some of the variables in greater detail.

Nye’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 Weertman (1973) and Smith (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 (Mottram, 2007) which show that Weertman’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 Van der Veen (1998) and Rist 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 (Mottram, 2007; Mottram and Benn, 2009) have shown that Nye’s approach performs as well as the more complicated LEFM approach of Van der Veen (1998) in predicting observed crevasse depths.

Benn 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 (Benn 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. O‘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 (Benn 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.

Alley and others (2008) proposed a simple law for ice-shelf calving which relies on the same basic assumption as Benn and others’ (2007a) model, that calving is dominated by cracks transverse to the flow. They differ in that Alley 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 Benn 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 Benn 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 Benn and others’ (2007a) model, related to the use of Nye’s (1955, 1957) formula for calculating crevasse depths. Consequently, we developed a further improvement to Benn and others’ (2007a) model, which we refer to as experiment 2. Experiment 3, which further departs from Benn 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 Nye’s (1955, 1957) 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 Benn 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 1

This experiment is a three-dimensional extension of Benn and others’ (2007a) calving model that parallels their implementation steps. In this experiment, we use the model-computed velocity at the surface, u_s(x, y), to calculate , and then Equation (4) to estimate the crevasse depth. Before proceeding we note that, in Benn and others’ (2007a) model, the velocity used to calculate crevasse depth using Equation (4) was the vertically averaged velocity. Benn and others’ (2007a) model assumes that the longitudinal strain rate is dominated by the along-flow variations of basal sliding, i.e. the creep component is small compared with basal sliding, as usually occurs in fast-flowing calving glaciers. Therefore, it does not matter whether the velocity used for computing the crevasse depth is that at the surface, at the bed, or the vertically averaged one, because their along-flow derivatives are nearly equal. However, this may not be the case, as discussed below in experiment 3.

Experiment 2

The model considered in this experiment involves an important conceptual difference to Benn and others’ (2007a) calving model. The derivation of Nye’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. Nye’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 Nye’s (1957) equation is an approximation which becomes better for wider glaciers, and probably very good for an ice sheet. Benn 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 Nye’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 Nye’s (1957) generalization of Glen’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 ₀. Nye’s (1957) generalization of Glen’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.

Experiment 3

This experiment involves another important conceptual difference, compared to Benn and others’ (2007a) calving model and the two previous experiments, all of which are underpinned by the assumption that the creep component is small compared with basal sliding, and that the along-flow deviatoric stress relevant to the opening of the crevasse will be very close to its value on the surface. However, strain rates, and corresponding stresses, are functions of depth within the ice. Consequently, the new improvement to the model considered in experiment 3 is based on finding the depth at which the model-computed tensile deviatoric stress, considered as a function of depth, equals the ice overburden closure pressure.

Experiment 4

A final improvement to our modelling of Johnsons Glacier calving is to use as the longitudinal strain rate in the flow direction minus the threshold strain rate required for crevasse initiation, . This was introduced in the theory of Benn and others (2007a), but their implementation into models of glacier dynamics assumes . Mottram and Benn (2009) incorporated the effect of such a threshold strain rate (referred to as ‘yield strain rate’) in the computation of crevasse depths from the strain rates measured at the surface of Breiðamerkurjökull, Iceland, using Equation (4). As they discuss, ‘the idea of a yield strain rate is introduced as a heuristic device to fulfil the role of a critical-stress intensity factor required to overcome the fracture toughness of the ice. The inclusion of allows the Nye model to be tuned, to allow for the observation that crevasses form only when the applied stresses exceed some (variable) value (Vaughan, 1993; Van der Veen, 1999)’. They used as a tuning parameter to adjust the predicted depths to those measured in the field. They tested a range of values of yield stresses (converted to yield strain rate in the model) of 10, 30, 50, 60 and 100 kPa, concluding that 60 kPa produced the best fit between observations and model calculations. We use this value in our computations, to test whether it improves the modelling of the Johnsons Glacier calving front position.

4. Full-Stokes Model of Glacier Dynamics

We consider the ice mass to be an incompressible and isotropic non-linear viscous fluid. Because of the extremely low Reynolds number (e.g. Re ∼ 10⁻¹³ for ice sheets; Schiavi, 1997), we consider a stationary quasi-static flow regime. Our glacier-dynamics model can be separated into three submodels (dynamical, thermal and free-surface evolution) whose unknowns are functions of space and time, and are solved separately, through an iterative procedure that uncouples the equations. A limited set of front positions is available for Johnsons Glacier, spanning a period of only 10 years. Moreover, the bedrock in the terminus area and the proglacial embayment is very flat and has a shallow water depth (just a few metres), resulting in a nearly constant glacier front position. Additionally, Johnsons Glacier is mostly temperate, though our radar studies reveal some patches of cold ice. Therefore, we restrict ourselves to the dynamical submodel, described below.

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. Paterson, 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. Hanson, 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 (Vieli 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 (Martí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. Vieli 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. Hanson, 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 (Shreve, 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. Quarteroni 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. Carey 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 (Martí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 Corcuera and others (2001). Martín and others (2004) give further details of the numerical-solution procedure, including the full set of finite-element equations.

5. Domain Gridding and Tuning of Model Parameters

Our finite-element grid is an unstructured grid made up of tetrahedral elements constructed by the Voronoi method, following Delaunay triangulation rules. Details of the procedure can be found in chapter 3 of Otero (2008). The grid used in the computations is shown in Figure 6, and consists of 5845 nodes and 3731 elements. This grid is a compromise solution between grid size and computation time, based on the results of a test of the sensitivity of the model to grid size.

Fig. 6. Finite-element grid of tetrahedra used for the model computations. Note that this figure has a different orientation than the maps of Johnsons Glacier in this paper, in order to properly show the location of the calving front in the grid. The latter corresponds to the area surrounded by the red curve in the foreground. The boundaries contoured with a blue curve correspond to ice divides. Those without any coloured curve represent margins where the glacier has contact with lateral walls.

The free parameters in our model are the stiffness parameter, B, in the constitutive relation and the coefficient, K, in the sliding law. Tuning of these parameters aimed to minimize the differences between computed and observed surface velocities. This misfit was calculated as

(17)

where the superscripts ‘c’ and ‘m’ denote computed and measured, respectively, and the summation extends over the N points at which measurements of ice velocity at the surface are available. Note that Equation (17) takes into account the misfit in both magnitude and direction of the velocity vectors, averaged over the total number of measurements. In the tuning procedure, B values were increased from 0.18 to 0.33 MPa a^1/3, in steps of 0.01 Mpa a^1/3, and K from 0 to 2 m a⁻¹ Pa⁻¹, in steps of 0.1 m a⁻¹ Pa⁻¹. The results of the tuning are shown in Figure 7a. Although an absolute minimum of the misfit (E ∼ 5.2 m a⁻¹) is found at K = 0.9 m a⁻¹ Pa⁻¹, B = 0.22 Mpa a^1/3, it lies within an elongated ‘valley’ defined by the misfit contour line E = 5.4 m a⁻¹, i.e. there is a wide range of values of K and B that produce nearly equal misfits. Within this valley, the misfits computation/observation corresponding to the highest values of both B and K (upper right side of the figure) are physically meaningless. They imply such a large amount of sliding that velocities corresponding to internal deformation are forced to become negative in order to fit the computed velocities to the observed velocities at the glacier surface. The model velocities computed for the choice of free parameters corresponding to the absolute minimum of the misfit mentioned above, unfortunately, do not provide an appropriate fit to the observed velocities, especially for locations close to the calving terminus, as Figure 7b clearly shows.

Fig. 7. (a) Tuning of the free parameters of the model, B and K, for the case in which a single value of K is used for the whole glacier. The plot shows the magnitude of the differences between the computed and observed velocities at the glacier surface, calculated using Equation (17). The exact location of the absolute minimum (E ∼ 5.2 m a⁻¹) is indicated by the red dot. (b) Comparison of computed (red) and measured (green) velocities at the glacier surface, for the choice of model parameters B = 0.22 MPa a^1/3 and K = 0.9 m a⁻¹ Pa⁻¹, with K constant across the whole glacier.

To minimize the differences between computation and observation near the terminus, we made a number of tests, selecting different powers for the basal shear stress and the effective pressure terms in the Weertman-type sliding law (p and q in Equations (14a) and (14b)), with unsatisfactory results. Finally, we succeeded in getting a good fit using p = 2 and q = 1 and setting different K values for the accumulation and ablation zones. (K is used for the former and 3K for the latter.) The results of the tuning are shown in Figure 8a, from which the best choice (misfit of E ∼ 4.1 m a⁻¹) for the model parameters is B = 0.23 Mpa a^1/3 and K = 0.8 m a⁻¹ Pa⁻¹ (this K value corresponds to the accumulation area; that for the ablation area would be three times as large). The corresponding velocity field is shown in Figure 8b, which manifests a clear improvement in the fit computation/observation as compared with that obtained with the earlier choice of model parameters. Differences still exist, especially near the calving front, which are more pronounced in direction than in magnitude of the velocity vector field.

Fig. 8. (a) Tuning of the free parameters of the model, B and K, for the case in which different values of K are used for the accumulation and ablation areas. In the latter, a value of 3K was used instead of K in the sliding law. The exact location of the absolute minimum of the misfit between computed and observed surface velocities (E ∼ 4.1 m a⁻¹) is indicated by the red dot. (b) Comparison of computed (red) and measured (green) velocities at the glacier surface, for the choice of model parameters B = 0.23 MPa a^1/3 and K = 0.8 m a⁻¹ Pa⁻¹, the latter corresponding to the accumulation zone, and the corresponding constant taken as 3K in the ablation zone.

For a highly crevassed terminus with substantial void spaces, one expects a flow law describing a continuous material to be inadequate. ‘Softening’ the material by assigning a lower B value in the crevassed region may partially allow for this. However, in our tuning we did not attempt to assign different values for B at the crevassed and non-crevassed areas, because the ‘boundary’ between them is not well defined.

Note that the tuning of the stiffness parameter, B, has implications for the computed crevasse depths. Crevasses will penetrate deeper in stiffer ice, due to slower creep closure rates. Equations (4) and (10) show that an increase in B implies a corresponding increase in d ₀. However, an increase in B also implies slower flow, and possibly a decrease in strain rate, thus counteracting (at least partly) the effects of the increase in B on crevasse deepening.

6. Results for the Predicted Front Position

The field of computed velocities corresponding to the optimal choice of model parameters is shown in Figure 9. We use this field for all the computations discussed in this section.

Fig. 9. Model-computed velocities corresponding to the optimal choice of model parameters. We use this velocity field to calculate the strains.

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 Benn 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 Benn 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 Mottram 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. Discussion

We first discuss the main hypotheses of our modelling and then its results and implications.

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

1. 1. The most important limitation of Benn and others’ (2007a) calving criterion is the plane strain assumption inherent in the use of Nye’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 Benn 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, Benn 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

Benn 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 Benn 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 Benn and others’ (2007a) calving criterion, i.e. water filling the crevasses is not necessarily used as a tuning parameter in Benn 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 Benn and others’ (2007a) calving criterion, with Nye’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 Benn 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

Benn 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. Mottram 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 (Shepherd 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; Otero, 2008), for which ground-penetrating radar measurements and geomorphological evidence suggest a polythermal structure (Navarro 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 (Hanson, 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 Paterson (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. Qin and others, 1994; Furdada and others, 1999; Benjumea and others, 2003), while more recent work (Molina and others, 2007; Navarro 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. Paterson, 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 Blatter 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.

8. Conclusions and Outlook

We draw the following main conclusions from this study:

1. 1. Benn and others’ (2007a) calving criterion is well founded in physical principles and is easy to implement. Being based on the crevasse penetration depth, computed directly from the strain-rate field derived from the velocity field, it has a direct link with glacier dynamics which allows its use in prognostic models.

2. 2. Our three-dimensional extension of Benn and others’ (2007a) calving criterion, combined with the use of a full-Stokes model of ice dynamics, was able to reproduce the presently observed position of Johnsons Glacier calving front only if a large amount of water filling the near-front crevasses is hypothesized. We attribute the failure to reproduce the observed front position to the fact that the hypothesis of plane strain inherent in Nye’s (1957) formula is not fulfilled for Johnsons Glacier.

3. 3. Using a modified criterion for crevasse depth, which computes the deviatoric longitudinal stress tending to open the crevasse using the full-stress solution, we substantially improved the results, to the point of avoiding the need to introduce water into the crevasses for the model to fit the observed front position. Nevertheless, this model slightly overestimated the crevasse depth necessary for triggering calving.

4. 4. The model that considers the tensile deviatoric stress opening the crevasse as a function of depth, and computes the crevasse depth by finding the depth at which this stress exactly balances the ice overburden pressure tending to close the crevasse, is the most physically plausible model and provides a good fit to observations if a small amount of water filling the near-front crevasses is hypothesized.

5. 5. Introducing a yield strain rate does not improve, for Johnsons Glacier, the fit between model computations and observations, unless we assume a certain (but reasonable) amount of water filling the crevasses. The existence of such a yield strain rate is physically plausible as a mechanism to overcome the fracture toughness of the ice.

6. 6. Our main difficulties in model implementation were related to the Weertman-type sliding law. The practical difficulties concerned estimating the height of the basal water column for a glacier lacking field measurements of basal water pressure.

In future work, we intend to apply our model to a glacier with a good record of front position changes, which is not available for Johnsons Glacier. The main weakness of the present implementation of the full-Stokes model lies in its inadequate representation of the basal boundary condition; further applications should improve this aspect.