Introduction

Iceberg calving is an important mass loss process, not only for the many smaller tidewater glaciers around the world, but also for larger ice sheets. Calving accounts for 50% of ice flux from the Greenland ice sheet and the majority of output from Antarctica (Reference Jacobs, Hellmer, Doake, Jenkins and FrolichJacobs and others, 1992; Reference Parizek and AlleyParizek and Alley, 2004). Recently, marine-terminating glaciers have received attention as dramatic retreats of calving fronts have been observed, particularly in Greenland (Reference Luckman, Murray, de Lange and HannaLuckman and others, 2006; Reference Howat, Joughin, Fahnestock, Smith and ScambosHowat and others, 2008). Wide speculation surrounding the cause of tidewater glacier retreats includes links to higher ocean temperatures (Reference Holland, Thomas, de Young, Ribergaard and LyberthHolland and others, 2008; Reference MurrayMurray and others, 2010), with decreased back-stress from reduced sea-ice coverage (Reference JoughinJoughin and others, 2008; Reference Amundson, Fahnestock, Truffer, Brown, Lüthi and MotykaAmundson and others, 2010) being proposed as a possible mechanism. Ice-flow models incorporating the representation of calving processes at the ice margin have potential to improve our understanding of the mechanisms surrounding glacier retreat; however, calving is a complex process which has proven difficult to model adequately.

The calving rate, U _C, of a glacier is defined as the difference between the velocity at the calving front, U _F, and the rate at which its length, L, changes:

Until recently, the two calving models adopted most widely for use in ice-flow models were simple empirical relationships, one relating the calving rate linearly to water depth at the front of the glacier (Reference Brown, Meier and PostBrown and others, 1982) and the other calculating changes in terminus position using a fixed ratio between flotation height and ice thickness (Reference Van der VeenVan der Veen, 1996). One disadvantage of these methods is that they have to be tuned to each glacier individually, using field measurements which are not always available. The flotation model also breaks down as a glacier moves through areas of basal overdeepening, where treatment of floating ice may be required (Reference Nick, Van der Veen, Vieli and BennNick and others, 2010). Recently, a new calving criterion has been developed that uses the location at which surface crevasses penetrate to the waterline (sea level) to define a new terminus position (Reference Benn, Hulton and MottramBenn and others, 2007). Although this crevasse-depth criterion does not represent the full complexity of processes thought to be involved in calving, it may provide a first-order control on predicted terminus location sufficient for realistic modelling.

Validation of the crevasse-depth criterion has been explored previously. A modified crevasse-depth criterion has been tested in a diagnostic (non-time-evolving) model by Reference Otero, Navarro, Martin, Cuadrado and CorcueraOtero and others (2010), who found that it successfully predicted the terminus position of Johnsons Glacier, a small glacier on Livingston Island, Antarctica. However, Johnsons Glacier is slow-moving and terminates in shallow water, so it does not provide an adequate test of the model’s performance in a fast-flowing outlet glacier context. It is also important to test the calving criterion in a prognostic (time-evolving) model. This was performed in a one-dimensional (1-D) idealized outlet glacier model and it was found that the crevasse-depth criterion is able to reproduce seasonal cycles of retreat and advance of the type observed for Greenland marine outlet glaciers (Reference Nick, Van der Veen, Vieli and BennNick and others, 2010). The work has been extended (Reference Vieli and NickVieli and Nick, 2011) to a model of Jakobshavn Isbræ, Greenland, uncovering high sensitivity and rapid adjustment of marine outlet glaciers to perturbations at their marine boundary.

Significantly, the sensitivity of the model to the entire range of its input variables has not yet been tested thoroughly. In this work we aim to investigate the sensitivity of the crevasse-depth criterion coupled to an ice-flow model to one key input variable, namely the depth of water in crevasses. In order to apply the crevasse-depth criterion with a realistic stress distribution, we use flowline data from Columbia Glacier, Alaska. However, rather than attempting to reproduce the behaviour of this glacier, we emphasize that our work is aimed at understanding the sensitivity of the crevasse-depth criterion to a key controlling parameter, the depth of water in crevasses.

Columbia Glacier

Columbia Glacier (Fig. 1) is a temperate tidewater glacier with a history of rapid retreat, measured at about 18km between the early 1980s and the present, which has been well observed and recorded. This comprehensive data record with its extensive bed elevation measurements makes Columbia Glacier particularly suitable as a test case for a tidewater glacier flow model. Our work uses a dataset composed of repeat aerial photography from which measurements of glacier length, surface elevation and speed exist (Reference KrimmelKrimmel, 2001; Reference O’Neel, Pfeffer, Krimmel and MeierO’Neel and others, 2005). We focus on the near-terminus region, extending ∼5 km upstream from the calving front and on short sub-annual timescales. Bathymetric measurements using side-scan sonar provide the topography downstream of approximately the 2003 terminus on a 5m grid (surveyed by the US National Oceanic and Atmospheric Administration (NOAA) in September 2005). Upstream of this region a continuity-based model (Reference RasmussenRasmussen, 1989; Reference O’Neel, Pfeffer, Krimmel and MeierO’Neel and others, 2005; Reference EngelEngel, 2008), constrained by airborne radar measurements provides basal topography. In the downstream regions covered by bathymetry, basal elevations are accurate to ±20 m, while in the region of continuity estimates errors are on the order of 25% (100 m). The bed topography and surface geometries show that the glacier remained grounded for the majority of its retreat, although it was observed to have a floating tongue temporarily between 2007 and 2009 (Reference Walter, O’Neel, McNamara, Pfeffer, Bassis and FrickerWalter and others, 2010), meaning that for the majority of the retreat we do not need to consider the effects of floating ice in the model.

Fig. 1. (a) Map of Columbia Glacier (source: US Geological Survey). Thick arrows and dotted line show the location of the central flowline used in this model. Location of 1993 calving front marked. (b) Oblique photograph of Columbia Glacier, with front marked to show retreat (source: R.M. Krimmel).

The Model

We set up a two-dimensional (2-D) model of Columbia Glacier using the observed bed topography and surface elevation from 1993 on an idealized flowline along the approximate centre of the glacier. The path of the flowline is shown in Figure 1 and the 1993 bed and surface elevation in Figure 2b. Flowline locations used in this paper reflect the horizontal distance (km) from the ice divide.

Fig. 2. (a) Comparison of observed and modelled horizontal flow velocity. (b) 1993 surface and bed elevation data along the central flowline of Columbia Glacier, showing different data sources.

The 1993 start time was chosen as it characterizes the dynamics of Columbia Glacier during its retreat, which is a period of particular interest. The associated length change and surface lowering in a time of retreat prohibits the use of an initialization (spin-up) period before model runs begin. Full Stokes computations result in long model run times, which also prevent the execution of an initialization period in the glacier’s steady-state position 15 years prior to 1993. This lack of initialization makes the model sensitive to surface elevation uncertainty in the digital elevation model (DEM) (Reference Zwinger and MooreZwinger and Moore, 2009), which could lead to inaccuracy in the simulated surface evolution. Despite this source of inaccuracy, the mid-retreat set-up is sufficient to provide a realistic glacier geometry in order to perform sensitivity testing.

There is some precedent for using 2-D ice-flow models for calving studies (Reference Vieli, Funk and BlatterVieli and others, 2001; Reference Oerlemans and NickOerlemans and Nick, 2005). Adaptations have been made using finite-difference models (e.g. Reference Nick, Vieli, Howat and JoughinNick and others, 2009) to include representations of lateral stresses and the cross-flow flux divergence as the valley width varies. However, three-dimensional (3-D) flow processes (lateral drag and horizontal spreading) are complicated to incorporate using finite-element methods and we have not attempted to do so here. Figure 1 shows several regions where we would expect this assumption to fail: convergent flow at km 35 and 50 where tributary glaciers join the main trunk and ∼km53 where the valley narrows. Lateral drag is also considered to have a stabilizing effect on glacier termini (Reference O’Neel, Pfeffer, Krimmel and MeierO’Neel and others, 2005; Reference Benn, Hulton and MottramBenn and others, 2007), but should only have a significant effect on the stress balance around the terminus if there are large along-flow gradients. Neglecting these 3-D processes is possible because of the broad and constant channel width (5–6 km) upstream of the 1993 terminus in our region of focus (km 55–60; Fig. 1).

The numerical model is implemented using the open-source finite-element software Elmer (http://www.csc.fi/elmer).

The velocity and stress profile within the glacier is calculated using Stokes flow for an incompressible fluid:

where _u is the 2-D velocity vector, _p is pressure, ρ _i = 910 kgm^–3 is the ice density (assumed constant throughout) and g is gravitational acceleration. The deviatoric stress tensor τ is related to the strain-rate tensor by Glen’s flow law:

where D is the strain-rate tensor and δ is its second invariant:

The Arrhenius factor A=56×10^–16 s^–1 (kPa)^–3 is an average of two experimentally observed values for temperate ice both measured in situ on glaciers (Reference PatersonPaterson, 1994).

Boundary Conditions

Measured flow speeds at Columbia Glacier demonstrate that sliding is the dominant component of surface motion (Reference Meier and PostMeier and Post, 1987; Reference O’Neel, Pfeffer, Krimmel and MeierO’Neel and others 2005), and borehole water pressures confirm low values of effective pressure (Reference MeierMeier and others, 1994). We thus choose to follow Reference Nick, Van der Veen and OerlemansNick and others (2007) in using a Weertman-type sliding law, modified to include the effect of basal water pressure on sliding velocity. Hence sliding velocity U _b is defined with respect to basal shear stress by the relation:

where N _eff is the difference between ice overburden pressure and basal water pressure, τ _b is basal shear stress determined from the glacier geometry and A _s is a sliding parameter. Basal water pressure is calculated assuming a direct connection between the glacier bed and the sea, such that

where P _w, the basal water pressure, depends on b, the depth of the bed below sea level, and ρ _w = 1000 kgm^–3.

The basal sliding parameter A _s has been chosen from a range around the value of 2.9×10¹⁴ (m kg)^1/2 a^–2 used by Nick and others (2007) to find a best-fit solution of modelled-to-observed surface velocity from 1993. The best result was found using two regions of different sliding coefficient. As can be seen in Figure 2b, there is a region of reversed bed slope towards the front of the glacier (km 52– 62). In this area, modelled velocities tend to be lower than observed, with the discrepancy probably caused by 3-D effects neglected in the flowline model. To account for the difference, in this region (down-glacier from km 55) we enhance sliding by decreasing the sliding coefficient to 9.1×10¹⁴ (m kg)^1/2 a^–2. For the rest of the glacier, the sliding parameter is 7.2×10¹⁶ (m kg)^1/2 a^–2. The values produce an average error of 249 ma^–1 over the km 50–60 region, which is relatively small compared with absolute velocity values of 7000 ma^–1 around the front. The fit of modelled-to-observed velocity can be seen in Figure 2a.

Accumulation and ablation data are taken from Reference TangbornTangborn (1997), who parameterized surface mass balance as a function of altitude using low-altitude precipitation and temperature observations and the area–altitude distribution of the glacier. The annually averaged surface mass-balance profile over the period 1949–96 is approximated by the empirically fit equation:

where B is the surface mass balance and S(x) is elevation in ma.s.l. The summer of 1993 had a mass balance very close to this average. Using this mass balance and the surface velocity, a new surface elevation profile is calculated at each time-step by the kinematic boundary condition:

The available geometry data reach back nearly to the ice divide, so the upstream boundary condition is zero horizontal ice flux. The glacier’s upper surface is a moving boundary, where node positions are updated accordingly to the solution of Equation (9). At the calving front, between calving events, the boundary moves freely, updated according to the calculated velocity profile, with a stress boundary condition accounting for water pressure acting against the flow below the waterline. Within the glacier, the mesh nodes are moved at each time-step in order to maintain an even distribution. For a more thorough discussion of the numerical methods used in Elmer to solve these equations see Reference Gagliardini and ZwingerGagliardini and Zwinger (2008).

Calving

In our model, calving is considered to occur where crevasses penetrate to sea level. Crevasse depth is calculated using the method suggested by Reference NyeNye (1955, Reference Nye1957), which assumes that crevasses will penetrate to the point at which the longitudinal tensile strain rate is balanced by the creep closure rate due to ice overburden pressure. Our model neglects any interaction between crevasses and the stress field in the surrounding ice. Linear elastic fracture mechanics (LEFM) can accommodate these conditions, but Reference Mottram and BennMottram and Benn (2009) showed that in practice the two methods produce similar results; therefore, we adopt the Nye method for simplicity. We also omit the effect of critical yield strain rate (or fracture toughness) on the crevasse depth. As a result, modelled crevasse depths may be larger than if fracture toughness were included as a parameter.

Following Reference Otero, Navarro, Martin, Cuadrado and CorcueraOtero and others (2010), strain rates are calculated along the flowline (x-direction) rather than in the direction of the vector of maximum longitudinal strain rate. The latter method is more accurate, but was found to give no appreciable difference in calculated crevasse depth and suffers from a higher computational cost. Crevasses can be said to penetrate to the point at which the horizontal component of the Cauchy stress tensor (σ_xx ) is zero. This depth calculation can be adapted to include the effects of water (depth D _w, measured from the base of the crevasse) within crevasses; hence, the base of the crevasse field occurs where

At each time-step, the crevasse-depth profile is examined to determine whether it crosses the waterline at any point (or multiple points) along the glacier’s length. If it does, the glacier is cut vertically at the point furthest upstream to form a new terminus position. A new mesh is then created using the updated surface profile and terminus position.

Model Experiments

The implementation of the crevasse-depth calving criterion in this model allows the prediction of discrete calving events. Whether the glacier calves depends on the stress profile around the front, and a sufficiently short time-step will ensure multiple time-steps between calving events. This is important, as in such a scheme calving at every time-step would produce a spurious dependence of calving rate on the time-step chosen. Therefore, we select the time-step of model runs in order to ensure that calving does not occur at every time-step.

Prognostic model runs were carried out with durations of 0.5 and 1.0 year, starting in each case from the 1993 geometry. To test the sensitivity of the model to the input variable, a variety of different crevasse water depths (D _w) were used, ranging from 0 to 10m. Results using crevasse water depths >10m were rejected as they led to calving occurring at every time-step. For most of the model runs, a time-step of 0.005 year (1.8 days) was found to produce reliable results, with multiple time-steps occurring between calving events. To test if the calculated calving rate was dependent on this time-step, the experiment with 9 m water depth was run again with time-steps of 0.002 year (0.7 day) and 0.01 year (3.7 days). For the time-step of 0.002 year, the simulated calving rate was found to differ by only 2% from the original result, although for 0.01 year calving starts to occur at every time-step, leading to an unreliable result. The experiment with 10m water depth used a 0.001 year (0.4 day) time-step, which was found to be the maximum time-step still to produce periods with no calving. Model runs on such short time-steps take a significant time to complete, so this experiment was only run over 0.5 year compared with 1 year for the others. Sensitivity to mesh grid size was also tested, with a range of grid sizes between 10 and 50m being used around the terminus, widening to 100 m further upstream. It was found that the varying mesh size produced a difference of ∼10m in the crevasse-depth profile, which is small compared with the size of calving events produced. Therefore a grid size of 50m around the front was used for most experiments.

Results

Overall, the behaviour of the calving model depends strongly on the crevasse water depth (Fig. 3). With no surface water in crevasses, the model glacier does not calve at all but steadily advances. Over the length of the model runs, for small crevasse water depths (5–7 m) the model glacier advances significantly, while for larger depths (8– 10 m), after an initial advance, the model glacier retreats (Table 1). The observed retreat rate in 1993 of 0.64 kma^–1 (Reference KrimmelKrimmel, 2001) falls between the values modelled with crevasse water depths of 9 and 10 m. At the end of the experiment, the different water-depth scenarios reach significantly different end points (Fig. 4).

Fig. 3. Computed evolution of terminus position over time for different depths of water in crevasses. Pale grey lines show the size of calving events and dashed lines show the position after each calving event to highlight the overall trend of retreat/advance (for clarity these lines exclude any events <20m in size). Lines for 9 m water depth are not shown as they do not differ significantly from those for 8 m.

Fig. 4. Glacier surface profiles for the initial model geometry in 1993 and at the end of each model run for different values of crevasse water depth (for 10m the model run is only 0.5 year compared with 1 year for others).

Throughout the experiment the models exhibit significantly different calving behaviour. Examining the calving rates (Table 1; Fig. 5) we see that calving rate increases with crevasse water depth, as would be expected given the calving criterion used. Model runs with 5 and 7m water depths exhibit infrequent relatively small calving events. For greater crevasse water depths the calving events are not only more frequent, but also tend to be larger, generally 200– 300 m, whereas the 5 and 7 m experiments show an average calving event size of 54 m.

Fig. 5. Modelled calving rates and final terminus position. Calving rates calculated using the average size of, and time between, each calving event in the model run. Glacier length change calculated over length of each run and in case of retreat after initial ten events.

Discussion

Our results show that modelled calving rate is highly sensitive to changes in crevasse water depth. Results were only reliable (with multiple time-steps between calving events) for crevasse water depths ranging between 0 and 10 m. Within this range we see a wide disparity in terminus behaviour, ranging from a 3.5 kma^–1 advance to 1.9 kma^–1 retreat. This effect may be partly due to the low surface gradient around the glacier terminus, as seen in Figure 4. These areas of low surface slope also exhibit a shallow crevasse-depth profile, so that small changes in crevasse depth can lead to large horizontal shifts in calving position. This shallow terminus is more pronounced in the modelled surface elevation than the observed and may be partly due to the lack of a lateral drag term, as lateral drag will tend to impede horizontal spreading around the glacier terminus. However, a flat tongue is a feature of many glaciers, although less pronounced at Columbia Glacier in 1993, so the result may still be considered interesting from a general modelling perspective.

One particular feature of the sensitivity to meltwater is the change in behaviour between 7 and 8m crevasse water depth, suggesting that the calving rate can be extremely sensitive to relatively small changes in this variable. This sensitivity to crevasse water depth may make the crevasse-depth calving criterion difficult to implement in a predictive ice-flow model. A 1 m change in crevasse water depth will produce a significant effect on the stress, in turn affecting crevasse depth and calving rate, but it is likely to be small compared with any potential measurement errors. Crevasse water depth is an inherently difficult property to measure, as the crevassed regions of any glacier are a dangerous and difficult working environment and the steep-sided narrow shape of crevasses makes them difficult to observe remotely. It is also a difficult property to estimate using surface mass-balance models, which would require estimates of crevasse width and length and also an idea of the rate at which water drains from the crevasse. Drainage rates are likely to be high in the extensively fractured regions around a calving face.

Sensitivity to crevasse water depth had been indicated in previous modelling work (Reference Vieli and NickVieli and Nick, 2011) and is an effect that may be expected, as many calving glaciers exhibit a strong seasonal cycle of winter advance and summer retreat, which may be explained by the presence of more surface meltwater in the summer months. As many studies have previously focused on ocean warming as a potential trigger for glacier retreat (e.g. Reference MurrayMurray and others, 2010) the result may also be interesting in highlighting the potential for atmospheric warming to cause glacier retreat.

In contrast to previous calving models, the model we have presented represents glacier retreat as a sequence of individual calving events. We find that generally a time-step can be chosen which ensures that there are multiple time-steps of glacier advance between each calving event. It seems likely that this is due to the fact that the maximum horizontal deviatoric stress is located some distance behind the calving face (Reference Hanson, Hooke and LeHanson and Hooke, 2003), which is not the area of lowest surface elevation. Hence a period of surface lowering is required before calving occurs again. The modelled increase in calving rate with crevasse water depth is characterized by a general increase in both frequency and size of calving events. This is the opposite behaviour to that expected from observation (Reference Walter, O’Neel, McNamara, Pfeffer, Bassis and FrickerWalter and others, 2010) and theory (Reference Amundson and TrufferAmundson and Truffer, 2010), where small icebergs calve frequently while large calving events happen only infrequently. This apparent contradiction can be explained by considering the simplifications made in formulating the calving model. Rather than being a representation of a specific physical event, each modelled calving event can be considered to represent a rapid change of terminus position, which could occur by the release of an individual large iceberg, or a disintegration into many smaller blocks. In many cases the region downstream of the calving point is also crevassed below sea level, so may be expected to disintegrate rather than break off as a single berg. It should also be noted that we would not expect a uniform crevasse field in a real glacier and stochastic variations in crevasse depth will control the exact timing at which an iceberg is released. Consequently, the model simulates average calving behaviour rather than being able to identify individual calving events. Nevertheless, the frequency of modelled calving events is physically meaningful. The model establishes event-driven models as a potential method to investigate short-timescale physics at the tidewater margin.

This model has the potential to provide a new method to investigate the effects of external controls on different styles of calving from tidewater glaciers. Methodological improvements are needed to switch from the kind of sensitivity experiment presented in this paper to an accurate representation of glacier behaviour. Future improvements should accommodate lateral drag in the modelled stress distribution, as well as seasonality in the input variables (such as crevasse water depth, which may be expected to vary with surface ablation). A realistic calving model should also represent other processes affecting calving, such as submarine melt and resistive stress from sea ice that may alter the stress profile around the terminus.

Conclusion

Recent debates have focused on causes for observed acceleration and retreat of calving glaciers, particularly those in Greenland. Previous modelling work (Reference Nick, Vieli, Howat and JoughinNick and others, 2009) suggests that penetration of surface meltwater to the bed of a calving glacier (causing increased basal water pressure and therefore increased sliding velocity) is unlikely to force observed levels of acceleration. Our results suggest that acceleration and retreat could be triggered by surface meltwater via enhanced fracturing and deepening of crevasses. This feedback increases calving rates, allowing higher air temperatures to play a key role in initiating retreat and dynamic instability.

The strong dependence of the modelled calving rate on crevasse water depth is likely to cause difficulty in applying the method to predictive glacier and ice-sheet models. Results may be significantly influenced by poorly constrained input data. A successful physically based crevasse-depth calving model must be constrained by better observations and include more thorough development of detailed physics. However, the ability to simulate glacier calving as a sequence of individual physically meaningful events means that it has potential to investigate calving style through different phases of advance and retreat and/or given variable geometry. Potential exists to greatly enhance our understanding of the complex interactions between calving and ice dynamics.