Modelling environmental influences on calving at Helheim Glacier in eastern Greenland

Calving is an important mass-loss process for many glaciers worldwide, and has been assumed to respond to a variety of environmental influences. We present a grounded, flowline tidewater glacier model using a physically-based calving mechanism, applied to Helheim Glacier, eastern Greenland. By qualitatively examining both modelled size and frequency of calving events, and the subsequent dynamic response, the model is found to realistically reproduce key aspects of observed calving behaviour. Experiments explore four environmental variables which have been suggested to affect calving rates: water depth in crevasses, basal water pressure, undercutting of the calving face by submarine melt and backstress from ice mélange. Of the four variables, only crevasse water depth and basal water pressure were found to have a significant effect on terminus behaviour when applied at a realistic magnitude. These results are in contrast to previous modelling studies, which have suggested that ocean temperatures could strongly influence the calving front. The results raise the possibility that Greenland outlet glaciers could respond to the recent trend of increased surface melt observed in Greenland more strongly than previously thought, as surface ablation can strongly affect water depth in crevasses and water pressure at the glacier bed.


Introduction
Calving is an important mass loss process for both ice sheets and smaller tidewater glaciers (e.g., Church et al., 2001;Rignot and Thomas, 2002;Blaszczyk et al., 2009) and plays a strong role in tidewater glacier dynamics (Meier and Post, 1987).In the past there has been debate about whether changes in glacier dynamics and terminus retreat, observed at a number of locations worldwide, were triggered by changes in calving rate or changes in glacier velocity, a debate made more complicated by the potential feedback between the two mechanisms (van der Veen, 2002).Modelling work has since shown that changes at the terminus are able to trigger retreat (Nick et al., 2009) and it is now generally thought that calving rates are affected by a number of external environmental variables.The four main hypothesized mechanisms for tidewater glacier retreat are as follows: 1. increased water in crevasses, which promotes calving by increasing the depth of fractures around the calving front (Benn et al., 2007a); 2. increased basal water pressure, leading to faster sliding, which can also enhance fracturing by increasing longitudinal strain rates (van der Veen, 2002); 3. undercutting of the terminus by submarine melt, which alters the stress field around the terminus and is hypothesized to enhance calving (Motyka et al., 2003;O'Leary and Christofferson, 2013); linked, as increased air and ocean temperatures are typically strongly correlated (Deser and Blackmon, 1993).In addition, the physical mechanisms relating ice dynamics to oceanic and atmospheric warming are also interrelated.For example, modelling of the proglacial fjord by O' Leary (2011) and Xu et al. (2012) shows that undercutting rates at the terminus are affected not only by the surrounding water temperature, but also by the rate of subglacial discharge, which is directly linked to glacier surface ablation and hence air temperatures.The interaction between different forcing variables and mechanisms makes it difficult to distinguish the precise causes of tidewater glacier retreat by observation alone.By including calving in an ice flow model it is possible to experiment with applying these environmental forcing factors separately and deduce the likely sensitivity of tidewater glaciers to them.Such a study has previously been performed by Vieli and Nick (2011) using a physically-based calving criterion based on the modelled depth of crevasses in the ice.However, the vertically averaged model used in previous studies (e.g., Nick et al., 2010;Vieli and Nick, 2011) is likely to lead to inaccuracies in modelled crevasse penetration, as earlier models have shown strong vertical variation in longitudinal stresses around the calving front (Hanson and Hooke, 2000;Hanson and Hooke., 2003).A previous publication presented a two-dimensional, vertically-resolved flowline model of Columbia Glacier, Alaska, using a similar calving criterion (Cook et al., 2012).That study was limited in scope due to the time-consuming, manual implementation of calving.In this paper we present an extension of the previous model which now automates the calving process, making model runs significantly faster and allowing a greater scope of experiments.The model has also been extended to include a parameterization of lateral processes and a non-linear sliding law, although it is still currently limited to representing grounded ice.
The chosen study location is Helheim Glacier, in eastern Greenland (Fig. 1a), which was selected due to its recent interesting behaviour.The mean annual terminus location of Helheim Glacier was stable throughout the late 1980s and early 1990s, although sub-annual variations of between 2 and 4 km are typical (Bevan et al., 2012).After a period of thinning in the 1990s (Abdalati et al., 2001), the front began to retreat in 2001.Over the period 2001-2005 the front retreated by 6 km, the main trunk sped up by around 3.0 km a −1 and the surface lowered by 40 m (Howat et al., 2005;Luckman et al., 2006).By 2006, thinning had stopped near the front and the glacier slowed, advancing approximately 4 km from its most retreated position.Since 2007 the glacier has once again shown a stable mean annual terminus position, although somewhat further retreated than the pre-2004 position (Howat et al., 2007;Bevan et al., 2012).Although large annual variation in terminus position is typical during some periods at Helheim Glacier, the rapid and consistent retreat of 2001-2005 was exceptional and a similar retreat was also observed in other glaciers in the region, leading to the con- clusion that it was likely the result of a single environmental event (Howat et al., 2008).The cause of this synchronous retreat has been widely debated, with high air and sea temperatures hypothesized as potential triggers (e.g., Joughin et al., 2008b;Murray et al., 2010;Christofferson et al., 2012).The bed elevation data available for Helheim Glacier are limited to ice-covered areas, where airborne radar measurements can be made; therefore, this retreat has also provided a modelling opportunity with data available to set up a model which has sufficient bed data available to experiment with a full range of terminus behaviour: retreating, steady and advancing.In this paper our time-evolving, flowline model of Helheim Glacier is tested for sensitivity to the four main environmental inputs outlined above.

Model setup
This study uses a vertical two-dimensional (2-D) ice flow model to analyse the stresses along a flowline within the glacier, which can then be used to determine whether and where calving occurs.The model is allowed to evolve through time, with calving applied according to the crevassedepth calving criterion proposed by Benn et al. (2007b) and used in previous modelling studies (Nick et al., 2010;Otero et al., 2010;Cook et al., 2012).The criterion states that calving will occur where the surface crevasse field penetrates below sea level.The depth of surface crevasses is calculated by the Nye formulation (Nye, 1955(Nye, , 1957)).Allowing for the effect of water in crevasses this means that the base of the crevasse field occurs where where σ xx is the longitudinal component of the Cauchy stress, D w the depth of water in the crevasse, ρ w = 1000 kg m −3 is the density of the water and g = 9.81 m s −2 is gravitational acceleration.The ice flow solution is calculated using the open source Elmer/Ice finite element modelling software (Gagliardini et al., 2013).Calving is applied using additional code which examines the output from the Elmer/Ice software and updates the geometry to allow for any calving events.If the model is run with a sufficiently short time-step, multiple time-steps will occur between each calving event, thus allowing individual calving events to be distinguished.In this case the model also provides data on the modelled size and frequency of calving events which can be compared to observed data.The calving criterion is not intended to precisely represent individual calving events, but to provide a way to link calving rates to ice dynamics in a physically-based manner, so that only the bulk calving behaviour of the model is physically meaningful.It should also be noted that the calving events do not necessarily represent individual icebergs, but instead a rapid change of terminus position which could also encompass the release of numerous smaller ice volumes.

Data sources
The study used data from Helheim Glacier at its point of furthest retreat in July 2005, since the terminus position at this date allowed 2 km of bed elevation data for the model to advance into; this enabled an investigation of a full range of glacier behaviour.Elevation and velocity were measured along a flowline close to the centre of the glacier, identified using surface debris features and velocity data (Fig. 1a).A flowline model is not well suited to representing areas of convergent flow, such as where the outlet glacier meets the Greenland Ice Sheet.Therefore, a point was identified at 152 km along the flowline where the behaviour of the glacier changes from ice-sheet-type, funnelling ice from a wide range of angles, to an outlet glacier where the flow is channelized and the velocity vectors are roughly parallel; the model was limited to locations downstream of this point.
Bed elevation was taken from a gridded data product available from the Center for Remote Sensing of Ice Sheets (CRe-SIS) (Gogineni, 2012); this product is derived from radar  (Leuschen et al., 2011).These data provide a horizontal resolution of 500 m and vertical error of around 20 m in areas of good bed signal, while areas of poor radar return use an interpolated bed with higher associated error.Surface elevation data were taken from a pair of ASTER images (July 2005), with photogrammetry providing a digital elevation model (DEM) with a vertical accuracy of ∼ 10 m (Murray et al., 2010).Surface velocity was measured by feature tracking between consecutive Landsat-7 ETM+ images (July, August and September 2005), the method and errors are described in detail in Bevan et al. (2012).The model output was compared to observed terminus position; this was measured at the intersection between the flowline and a manually digitized ice front, mainly using images from Landsat-7 ETM+, but also supplemented by ASTER, SPOT 5 and airborne LiDAR images where available (Fig. 2b).Errors on these terminus positions have been estimated to be less than 100 m (Bevan et al., 2012).The resulting initial geometry of the model is shown in Fig. 1c.

Ice flow model
The geometry is used to find the velocity and stress fields within the glacier using Stokes flow for an incompressible where u is the 2-D velocity vector in the vertical plane, τ is the deviatoric stress tensor, p is pressure, ρ i = 910 kg m −3 is the ice density and g = (0,0,-9.81)m s −2 is gravitational acceleration.The deviatoric stress is related to the strain rate (˙ ij ) by the constitutive relation: where the effective viscosity µ is given by where ˙ 2 e = tr(˙ 2 )/2 is the square of the second invariant of the strain rate.We use n = 3, as it is the most commonly used exponent in glaciology.The Arrhenius factor A depends on the temperature of the ice relative to the pressure melting point (T ) according to the Arrhenius equation: where R is the universal gas constant, and Q, the creep activation energy, and the constant A 0 depend on the ice temperature.For this study standard values of these constants were taken from Paterson (1994).In the absence of observational data, the ice temperature profile with depth was approximated from a thermomechanical modelling study of Jakobshavn Isbrae, a similar outlet glacier on the west coast of Greenland (Funk et al., 1994).The ice temperature against relative depth was estimated by fitting a quadratic function to the modelled profile in Funk et al. (1994) where the ice temperature T is measured in • C, z is the height above the bed and H is the ice thickness.This temperature profile means that the basal ice never reaches the pressure melting point.

Boundary conditions
The glacier surface is assumed to be traction-free, and elevation changes are calculated by the equation for the kinematic boundary condition: where s is the surface elevation and M the surface mass balance.Mass balance values were estimated from Andersen et al. (2010), who directly observed melt from stakes placed along the glacier in two consecutive melt seasons in 2007 and 2008.These measurements were assumed to approximate the entire summer melt.Using a polynomial fit of ablation rates against elevation this leads to a summer ablation rate of: which places the entire model domain in the ablation area.In winter, surface mass balance was assumed to be zero since the dynamic effect of any winter snow is negligible if it all melts in the succeeding summer.
Since the model uses only a two-dimensional geometry, technically there are no lateral boundary conditions to apply.However, this neglects important 3-D processes such as the confluence of tributaries, the effect on volume conservation of changes in channel width, and lateral drag from the glacier margins, which has been suggested to have a stabilising effect on glacier termini (O'Neel et al., 2005).Little can be done to address the effect of tributaries without switching to a 3-D geometry, but by measuring the channel width at the surface and making assumptions about the channel geometry we are able to estimate the lateral drag and flux contribution from changes in channel width.We assume a trapezoidal channel shape with a width at the base two thirds of that at the surface, a constant lateral velocity gradient and sliding at the margins of half the centreline velocity.This leads to a volume contribution from lateral spreading or constriction of ice (φ): where H is the ice thickness and U x and W are the velocity and channel width at the surface.This was added as an additional surface mass balance term.The geometry can also be used to derive a vertically varying lateral drag term: where w(z) is the channel width at an elevation z and u x the velocity along flowline.The vertically varying force (F ) was imposed on each element in the glacier body, opposing flow.
To initialize the model and remove the destabilizing effects of errors in the DEMs used, a spin-up run was performed allowing the surface elevation to evolve until its rate of change came into line with annual mass balance values.This process provides a stable starting point for experiments, although it also means that there are some areas of difference between the modelled and observed geometry and surface velocity profile (Fig. 1b and c).These differences occur where the model fails to reproduce the full flow behaviour of the ice, such as regions where 3-D flow features are significant, or where there are errors in the bed DEM, which have been shown to have a significant effect on modelled velocities (Zwinger and Moore, 2009;Seroussi et al., 2011).
At the bed, a no-penetration boundary condition was applied.Basal drag was implemented in the model using the non-linear scheme proposed by Gagliardini et al. (2007).The drag depends on the effective pressure (N), which is defined as the difference between ice overburden pressure and water pressure at the bed (P w ) which was taken to be equivalent to depth below sea level.
In this model the basal drag is given by where A s is the sliding parameter without caviation (as used in Weertman-type sliding laws), C is the maximum value reached by τ b /N and u b is the basal sliding velocity.The constants C, A s and q were tuned to find the best fit to the observed surface velocity when the model reached its relaxed state.The values used were C = 1.0, q = 2.0 and a good fit required two different sliding parameters, A s = 1.66×10 −58 (m kg) 0.5 a −2 where x <166.1 km and A s = 2.76×10 −57 (m kg) 0.5 a −2 where x >166.1 km.
At the inflow boundary a constant velocity of 3500 m a −1 was applied (an average of the observed surface velocity of 4006 m a −1 and the modelled sliding velocity of approximately 3000 m a −1 ).See Fig. 1b for the surface velocity profile after spin-up.
At the terminus, water pressure was applied to the calving face, which was allowed to evolve according to the modelled velocity profile.In experiments where an additional backstress was applied to the calving face to replicate the effect of ice mélange, this acted parallel to the x axis, in the opposite direction to ice flow, over an appropriate vertical range.In experiments designed to replicate submarine melt, nodes on the calving face were negatively displaced, parallel to the x axis and according to the assumed melt rate.Because the applied water pressure acts perpendicular to the boundary, in these experiments the model was able to fully represent the buoyancy forces arising from the undercut cavity.However, the model is unable to represent floating ice caused by detachment of the basal ice boundary from the underlying bedrock, and model runs in which the glacier reached the flotation point were terminated.

Model experiments
The first experiment performed was a control run with no external forcing, in order to produce a base set of terminus position and calving event size data.This was used for sensitivity analysis and model validation.The model was then tested using the four external variables hypothesized to affect terminus behaviour: depth of water in crevasses, basal water pressure, undercutting of the terminus by submarine melt and backstress from ice mélange.The results from the crevasse water depth experiments were also used to identify cases with steady and retreating terminus behaviour; these settings were then used to test the other variables in scenarios with differing terminus behaviour.Each variable was applied in a seasonally varying manner, with maximum values occurring in a summer season, assumed to last from May to September in most cases.The only exception was experiments using backstress from ice mélange, which was applied over a reduced winter season of January-May, effectively creating an additional autumn season (September-January) in which no surface ablation or undercutting occurred and no backstress was applied.The precise values used for each variable are discussed separately in the relevant results section for the sake of clarity, and a list of all experiments performed is provided in Table 1.
In each experiment, the time-step was chosen to ensure that multiple time-steps occurred between calving events.This choice is important: if calving were allowed to occur at every time-step there would be a spurious dependence of calving rate on the time-step.The time-step necessary depends on the environmental forcing used; in cases where the model produced more frequent calving events, a time-step of 0.001 yr (8.8 h) was used, while for cases with less frequent events (such as those with low crevasse water depth) a maximum of 0.003 yr was chosen.The model used a variable mesh size of 100 m near the upstream boundary and 40 m at the calving front.

Control run results
In an initial control experiment the model was allowed to run freely for a 5 yr period, with no water in crevasses or other external forcing variables.Over the 5 yr period the model produced 130 calving events (or one every 14 days, on average), with the size distribution roughly normally distributed around a mean length of 479 m (Fig. 2a).Comparing the modelled and observed terminus behaviour over this period, we see that both the model and the real glacier advanced substantially from their position in July 2005, although the real glacier advanced somewhat sooner and further than the model (Fig. 2b).After 1 yr the model's terminus position stabilized, whereas the true glacier continued to fluctuate in length, probably in response to external climate factors.
www.the-cryosphere.net/8/827/2014/The Cryosphere, 8, 827-841, 2014 Table 1.Full list of experiments performed, arranged by results section.S and W refer to summer and winter values.D w is the crevasse water depth (cwd), dP is the change in hydrostatic head applied to the basal water pressure, M f is the maximum melt rate applied at the calving face and F f is the force applied to the calving front over a vertical range of −120 to 15 m a.s.l.In ice mélange experiments where both the cwd and the backstress vary, the length of seasons in cwd and backstress do not match, meaning that there is also an autumn season where both the cwd and the backstress are set at zero.

Cwd
Basal water pressure Undercutting Ice mélange

Sensitivity analysis
To test if the chosen mesh density and time-step affected the results, terminus position and calving event size were examined using a range of time-steps (0.001-0.01 a) and grids (10-50 m).A grid size of 40 m was selected on the basis that it was able to resolve all calving events identified with the 10 m mesh, and produced only a 0.2 % change in mean calving event size, and 1.3 % change in terminus position at the end of a 5 yr run, compared to a grid size of 30 m (the maximum density which could be run in serial due to memory constraints).Changes in timestep size of only 0.001 a had a significant effect on the mean size of calving events produced (up to 35 % in some cases), as smaller time-steps allowed smaller calving events to be resolved.However, the 0.003 a timestep chosen had a less than 1.0 % change in final terminus position compared to a 0.001 a timestep, giving us confidence in the terminus position as a robust model result.
The model was also tested for sensitivity to errors in the input data used, by examining the terminus position over a 5 yr run in advancing, steady and retreating experiments (Cook, 2012).No observable effect on the evolution of terminus position arose from changes in the length of the surface relaxation period, the precise formulation of the lateral drag term, the applied surface mass balance or minor changes in the bed elevation profile.Changes to the inflow velocity boundary condition (±500 m a −1 ) and the mass balance term used to correct for 3-D flux properties (50-110 %) were found to result in minor differences in terminus evolution -the final terminus position remained the same, but the timing of the advance or retreat differed by an average of 0.75 a with a maximum of 2.3 a. Changes in the basal water pressure input to sliding of up to ±100 m showed similar results, though in one case (with an altered hydrostatic head of 100 m) a change in final terminus position of 2.5 km occurred.These three input parameters have the potential to affect results, but the impact of small errors is likely to be minor.The only input variable to have a substantial effect on model results was the ice temperature profile used.A change of ±5 • C in the glacier ice was able to alter final terminus position by up to 4.0 km.To reduce the impact of uncertainty in the ice temperature profile on future modelling work, an extension of the model to include full thermomechanical coupling should be used, preferably constrained by observations of ice temperature through borehole measurement.

Model validation
The nature of these experiments makes it difficult to use normal validation methods to compare the model's performance against observed measurements.The terminus evolution of the model may be tuned to a wide range of behaviour by altering the crevasse water depth, which is a poorly known input variable; therefore, comparing terminus behaviour on its own was not an appropriate validation method.Instead we used the control run, which provided a reasonable fit of modelled to observed terminus position, and chose two methods of comparing the model output to known glacier behaviour: comparing calving event size distribution and the glacier's dynamic response to calving.
The model dynamics may be tested by examining the effect of calving events on glacier velocity.It has been reported by Nettles et al. (2008) that large calving events at Helheim Glacier occur at the same time as rapid increases in flow speed of up to 15 % near the calving front, decreasing further upglacier.When the modelled velocity was examined immediately after large calving events (400-450 m), a similar increase in velocity was observed of up to 10 % near the front, decreasing to around 1 % 20 km upstream.This behaviour provides evidence that the model provides a reasonable representation of the interaction between calving and glacier dynamics.
To test the performance of the calving law, we compared modelled calving events with similar observations from satellite imagery.A calving event is a retreat in the terminus position over a single timestep, which is directly comparable to the high frequency observations of terminus position which can be obtained from satellite imagery.During its retreat, Helheim Glacier was observed to produce distinct calving episodes days to weeks apart, during which the terminus retreated by 0.5 to 1 km (Joughin et al., 2008b).This qualitatively agrees with the style of calving produced by the results presented here, with a typical calving length of 479 m and frequency of every 14 days.A similar calving event size was also produced by a recent particle-based model of Helheim Glacier (Bassis and Jacobs, 2013).The approximately normal calving event size distribution produced by the model disagrees with other published studies which observe a power-law type distribution of iceberg sizes, where small ice losses are the most common (Dowdeswell, 1989;Chapuis, 2011;Åström et al., 2013).This difference is to be expected given the mesh resolution of 40 m and time-step of 0.003 yr used in most experiments -this is insufficient to model very small calving events, but able to capture large changes in terminus position which are likely to dominate long-term behaviour.The discrepancy is also affected by the fact that each calving event predicted by the model does not necessarily occur by the release of a single iceberg but could also indicate a rapid disintegration into smaller ice masses.In fact, this second possibility may be the most likely; in most of the modelled calving events, areas of ice downstream of the predicted calving point also become crevassed below the water line in the same timestep.

Crevasse water depth
Few data are available on the depths of water in crevasses; therefore, a range of crevasse water depths (henceforth abbreviated as cwd) from 0 to 50 m was selected to represent a wide range of terminus behaviour (see Table 1).As would be expected, larger water depths had the effect of deepening the base of the surface crevasse field (Fig. 3a), causing more frequent calving events to occur.The mean calving rate increased from 10.3 km a −1 for the experiment with zero cwd to 38.2 km a −1 with 50 m cwd.The higher calving rate had a significant effect on the terminus evolution (Fig. 3b), with greater crevasse water depths causing the terminus to retreat up to 17.5 km, while for lower water depths the terminus advanced by up to 3.9 km.
In some experiments the terminus reached the flotation point, after which point the experiment was terminated as the model cannot represent floating ice. Figure 3b also shows that the terminus position at times became fixed at a particular point.These locations typically coincided with a local bedrock maximum, on which the terminus was pinned as advance of the model into deeper water triggered calving.This behaviour echoes results using empirical calving laws which also inhibit advancement into deeper water (Nick and Oerlemans, 2006).However, unlike the experiments using empirical calving laws, as the terminus geometry evolved our model www.the-cryosphere.net/8/827/2014/The Cryosphere, 8, 827-841, 2014 was able to advance past bedrock overdeepenings, although this often coincided with the model reaching flotation point.This behaviour mirrors the observation by Joughin et al. (2008b) that Helheim Glacier advanced past an overdeepening by forming a floating ice tongue to "bridge" the gap.

Basal water pressure
Observations at Helheim Glacier have linked basal water pressure and velocity; a study by Andersen et al. (2010) used a surface mass balance model to find a correlation between observed surface velocity and estimates of melt water availability.The relationship between surface meltwater and velocity is expected to be complicated; velocities tend to increase with increasing basal water pressure, however, the basal water pressure depends not only on the flux of water to the bed but also on the type of drainage system at the bed.If sufficient water is added to the system it may trigger a change from a distributed to a more efficient channelized drainage system, thus decreasing basal water pressure (Bartholomew et al., 2010;Howat et al., 2010).Generally, the percentage of seasonal changes in velocity appear to be higher on the ice sheet than on outlet glaciers, but even outlet glaciers can show seasonal variations in speed of up to 15 %, thought to be caused by changes in basal water pressure (Joughin et al., 2008a).Satellite observations of Helheim Glacier have shown sub-annual velocity changes of approximately 20-30 % (Bevan et al., 2012); however, it is unclear whether these arise mainly from changes in basal sliding conditions or coincident changes in terminus position.
The model was tested with a summer increase in basal water pressure of 50 m (measured as change in hydrostatic head), equivalent to a change in velocity of approximately 15 % as observed by Joughin et al. (2008a), which was designed to represent a spring speed-up event.Since the observed velocities were derived during the period of highest sliding, this meant that a new set of sliding parameters had to be used for this experiment.Compared to the control experiment the model with variable basal water pressure has a lower average sliding velocity, due to the lower sliding rate in the winter period.
Although the seasonal change in basal sliding had relatively little effect on the instantaneous longitudinal deviatoric stress, causing a mean percentage change of only 1.2 % (Fig. 4a), it decreased both the mass flux reaching the calving front and longitudinal extension, thus having a long-term effect on the predicted crevasse field.Two sets of experiments were performed, first using cwd fixed at 0 m and second using a seasonally varying cwd, with a summer value of 30 m. Seasonal variation in basal water pressure had a significant effect on the modelled terminus behaviour in both sets of experiments (Fig. 4), with lower winter sliding velocity causing a decrease in calving rates and hence a relative advance of the terminus.For the experiment with 0 m cwd the change in mean calving rate was not statistically significant, although the final terminus position was altered by 1.9 km.For the experiment with 30 m cwd, the final terminus position changed by 7.3 km and the mean winter calving rate was 5.1 km a −1 lower than during the summer high sliding period (Fig. 4e).

Undercutting by submarine melt
Undercutting of the calving face by submarine melt is hypothesized to affect the surrounding stress balance by removing supporting ice, and hence affect calving rates (O' Leary and Christofferson, 2013).Recent observations around Greenland suggest that submarine melt rates in this region may be sufficient to have an effect on calving; studies have used either measurements of fjord temperature and water velocity to estimate the heat flux transported to the glacier terminus, or used the difference between ice flux to the grounding line and surface ablation to infer average The Cryosphere, 8, 827-841, 2014 www.the-cryosphere.net/8/827/2014/melt rates across the calving face.Estimated mean melt rates range from 11 to 1424 m a −1 for glaciers across Greenland (Rignot et al., 2010;Enderlin and Howat, 2013) and between 205 and 1238 m a −1 for Helheim Glacier (Sutherland and Straneo, 2012;Enderlin and Howat, 2013).Based on these estimates, a range of maximum undercutting rates from 1000 to 5000 m a −1 were applied parallel to the flowline, with a vertical profile of zero melt at the waterline, increasing linearly towards the bed.During the winter season a fixed maximum melt rate (m.m.r.) of 150 m a −1 was used, based on plume modelling results by O'Leary (2011).A full list of values used is given in Table 1.
The effect of undercutting on the longitudinal deviatoric stress distribution around the front is shown in Fig. 5a.There are areas where the stress is clearly affected, particularly at the bed around a kilometre behind the front, where the stress is affected by increased buoyancy forces acting on the ice.However, the mean change in the stress profile is only 0.6 % from the original profile, and the change in crevasse field depth is relatively small, with a maximum of 12 m.As might be expected from this, the undercutting applied had relatively little effect on the modelled terminus behaviour, with all experiments reaching the same position after the 5 yr model run (Fig. 5b).The experiments did show some difference in the timing of the advance, which was delayed by up to a maximum of 0.44 a for the 5000 m a −1 experiment.In order to test this result further, a wider range of undercutting rates was tested to see at what point the terminus would begin to retreat.As shown in Fig. 6, the rates required to cause terminus retreat were between 10 and 20 km a −1 , around an orderof-magnitude higher than published estimates of melt rates.

Ice mélange
Ice mélange is a mixture of calved ice and sea ice, often observed in front of Greenland glaciers.In winter months this mélange can freeze solid, and provide a stress opposing the flow of the glacier, which has been hypothesized to inhibit calving.Relatively little is known about the mélange due to the difficulties in performing studies in the proglacial fjord environment.In order to include the effects of the ice mélange in the calving model both its thickness and the stress it applies on the calving face must be known.Using Li-DAR data from July 2007, the ice mélange in front of Helheim Glacier was measured to have a typical freeboard of 10-20 m, discounting isolated larger icebergs.We use a mean freeboard value of 15 m, equivalent to a full thickness of approximately 135 m.The only study to measure the stress exerted on the calving face used a coincident change in velocity of 15 % upon the break up of the ice mélange at Store Gletscher, in western Greenland, to estimate a backstress of 30-60 kPa over the full calving face, equivalent to a force of 1.8-3.6 × 10 7 Nm −1 (Walter et al., 2012).Initial experiments showed that values lower than this had no observable effect on the model, therefore experiments were performed using a range of backstress from 0.0 to 50.0 × 10 7 Nm −1 over the vertical range of 120 m below to 15 m above sea level, applied during the winter season (see Table 1).
The effect of two of the backstress values on the modelled longitudinal deviatoric stress profile can be seen in Fig. 7.A backstress of 5.0 × 10 7 Nm −1 caused very little difference to the stress field around the front (mean absolute change of 8.5 %), while for a backstress of 25.0 × 10 7 Nm −1 the change in the stress field was 60.0 %.The effect on the modelled crevasses was also more significant, showing a clear reduction in longitudinal crevassed extent of 490 m in the 25.0 × 10 7 Nm −1 experiment, while the 5.0 × 10 7 Nm −1 experiment showed a change in crevassed extent of 25 m (Fig. 7).When the modelled terminus evolution was examined, only the very highest backstress was seen to be able to affect terminus behaviour (Fig. 8).For both sets of experiments, only a backstress of 50.0 × 10 7 Nm −1 caused the final terminus position to change, although in the 0 m cwd experiment the lower values changed the timing of the glacier's advance by up to 0.31 a.

Discussion
The modelled glacier behaviour is shown to be significantly more sensitive to environmental factors related to air temperature and surface run-off than to oceanic forcing.This is of particular interest given that surface mass balance models generally agree that surface ablation has been steadily increasing in Greenland over the past two decades (Vernon et al., 2013).The strong effect of water in crevasses may be expected as this relationship arises axiomatically from the calving model used.Selecting the values of cwd to use was challenging as no measurements of observed water depth have been published, and likely values are difficult to predict given the changing volume of crevasses as ice advects downstream, and the probable high rate of drainage in the fractured area around the terminus.Observations have shown that water-filled crevasses can occur even near the calving margin (Danielson and Sharp, 2013), so it is at least known that high water depths are possible in this region, and may become more frequent if the trend of increasing surface ablation continues.
If a hydraulic connection exists between the glacier surface and the bed, surface meltwater can also affect the terminus behaviour by altering basal water pressure, and hence sliding velocity.Observations of seasonal velocity variations of 15 % on Greenland outlet glaciers had previously been concluded to be insufficient to have a large effect on terminus behaviour (Joughin et al., 2008a).The results presented here indicate that changes in basal water pressure leading to velocity changes of this magnitude can in fact have a strong influence on modelled calving rates.The finding supports previous tidewater glacier modelling work (Vieli et al., 2000) and observations at other tidewater glaciers (Kamb et al., 1994;Sugiyama et al., 2011;Danielson and Sharp, 2013) which had already observed a significant dynamic response to changes in basal water pressure, and highlights that changes in basal conditions should not be excluded as an important influence on tidewater glaciers.The link between changes in air temperature and run-off, and changes in velocity may, however, be complicated.Adaptations in the basal hydrology system can mean that an increased input of water to the bed results in a more efficient drainage system and hence reduces basal water pressure (Bartholomew et al., 2010).In the context of previous debate over whether dynamic changes in tidewater glaciers are triggered at the terminus or by changes in velocity (e.g., van der Veen, 2002), the results from these two sets of experiments (crevasse water depth and basal water pressure) suggest that either mechanism may have a strong influence on glacier dynamics.
In contrast to previous results, the insensitivity to undercutting by submarine melt shown in these experiments was surprising.Previous research has indicated that a link between glacier retreat in southeast Greenland and ocean warming was likely (Murray et al., 2010;Vieli and Nick, 2011), with retreat coinciding with warm ocean temperatures.Previous modelling work has also indicated that undercutting has the potential to cause tidewater glacier retreat (Vieli and Nick, 2011), and in modelling work using a similar two-dimensional, high-resolution solution of stresses around the calving front, cavity formation has been shown to have a significant effect on the stresses around the terminus (O'Leary and Christofferson, 2013).The study by O'Leary and Christofferson (2013) was diagnostic, not time-evolving, and we hypothesize that this was the reason for the difference between their results and ours.
To test this hypothesis, we set up two experiments: one applied an undercut cavity of varying length to a fixed terminus while the second applied a fixed undercutting rate, but also allowed the ice around the terminus to evolve as the undercut cavity size increased.The results in Fig. 9 show that when the terminus geometry was fixed, our model reproduced the linear relationship between cavity length and stress displacement previously reported (O'Leary and Christofferson, 2013).However, when the ice was allowed to deform The insensitivity to undercutting would be significantly different if the undercutting were to be concentrated at the waterline (i.e.notching), as has previously been observed in lake-terminating glaciers (Kirkbride and Warren, 1997;Benn et al., 2001) and tidewater glaciers in Svalbard (Vieli et al., 2002).In that case the remaining ice block would be unsupported and highly likely to shear off from the glacier, leaving behind a submarine ice foot.This style of calving was not considered in our model.The terminus of the model remained on a downward sloping bed throughout the experiments, which may also have contributed to the observed insensitivity to undercutting.Previous work on ice sheets has shown that glaciers may be stable even under grounding-line retreat on downward facing slopes, but much more sensitive on reverse bed slopes (Schoof, 2007).To fully determine the sensitivity of tidewater glaciers to undercutting, future work should include experiments on a wider range of bed topography.
The experiments presented also showed little relationship between backstress and terminus behaviour, except in experiments in which the backstress was likely to be unrealistically high.Although little is known about the likely magnitude of backstress provided by a proglacial mélange, it is unlikely to exceed the shear strength of sea ice, which has been measured at 550 ± 120 kPa (Frederking and Timco, 1984), while our experiments used stresses of up to 3.7 MPa.Although these results show that backstress from ice mélange is unlikely to have a direct effect on the crevasse field around the terminus, it has been hypothesized by Amundson et al. (2010) to affect calving rates by exerting enough force to prevent a calved block of ice from tipping, thus keeping it pinned in contact with the main body of the glacier.Such large blocks of ice could feasibly have a much greater effect on the stress at the front of the glacier, so the winter ice mélange should not be discounted as having a possible effect on calving rates, and future modelling work should include this type of process.
Model validation using terminus evolution, calving event size and dynamic response to calving demonstrates that the model reproduces some key aspects of tidewater glacier behaviour.However, there are areas in which it could be improved; in its current form the model has a number of shortcomings, particularly the approximated ice temperature profile, and the oversimplified parameterization of lateral drag.The simplified lateral drag means that the model cannot account for lateral variations in bed topography and sliding rate, which could affect the results.This problem could be addressed by scaling the model up to three dimensions.The current temperature field could be improved by using  thermodynamically coupled equations to solve for the ice temperature.The current model presents a thermodynamic contradiction in that basal temperature is below the pressure melting point, while our sliding law assumes a non-zero basal water pressure.This was necessary to perform our basal sliding experiments, however an improved solution for ice temperature would allow us to apply basal sliding in a more rigorous manner.Aside from these shortcomings, the fact that the terminus behaviour depends strongly on poorly constrained input variables such as the crevasse water depth means that at present it is not suitable for use in predicting future tidewater glacier behaviour.If the model could be more thoroughly validated it would be a powerful tool to investigate the behaviour of tidewater glaciers, capable of including a wide range of environmental inputs and producing a detailed set of model outputs for comparison to observation.First, it would also have to be tested using a better constrained data set with reliable basal elevation, ice temperature and crevasse depth data.Of these, the data on crevasse depths would be hardest to procure, and are unlikely ever to be available on a wide scale.Although this is a significant limitation of the model, it is possible that with improved data sets for validation and refined models of crevasse depth it could be overcome.
Despite this limitation in the model presented, by using consistent crevasse water depths we have produced a set of properly controlled experiments which allow the impact of other environmental variables on calving to be investigated.The results show that both undercutting of the calving face by subaqueous melt and backstress from ice mélange have little impact on calving in the model.This result does not depend on the crevasse water depth, and is observed in experiments using both 0 m cwd and a seasonally applied depth of 30 m.In the case of basal sliding the model was observed to respond more strongly to changes in basal water pressure when combined with a larger crevasse water depth.Although we can conclude from both sets of experiments (0 m and 30 m cwd) that a change in basal sliding caused a change in terminus evolution, the strength of the response seems to be affected by the crevasse water depth used and further experiments using a better validated model of crevasse depth would be required to assess the magnitude of the response more rigorously.

Conclusions
This paper presents a flowline, tidewater glacier model, capable of representing the calving of grounded ice using a physically-based calving criterion.The results demonstrate The Cryosphere, 8, 827-841, 2014 www.the-cryosphere.net/8/827/2014/that it reproduces key aspects of observed calving behaviour, examined qualitatively by comparing the size and frequency of calving events, and the subsequent dynamic response.However, it would benefit from better data availability for more rigorous validation.In response to previous debate over whether glacier retreat is triggered first by glacier acceleration or increases in calving rate, the model results indicate that either process has the potential to cause terminus retreat and it may be more enlightening to consider the specific mechanisms involved.
The results presented show a surprising insensitivity to oceanic influences and a noticeable sensitivity to the effects of surface ablation.The lack of sensitivity to simulated ice mélange may be because of the relatively simple representation of its effect on the calving front used in the model.However, the weak effect of submarine melt on calving rates appears to arise from the ability of ice deformation in the model to counteract the stress displacement caused by a cavity in the calving front.The strong effect of both crevasse water depth and basal water pressure on the model's behaviour raises the possibility that outlet glaciers may be strongly affected by recent trends in surface ablation.The results demonstrate the potential for this type of model to improve understanding of tidewater glacier dynamics, their response to external forcing, and to diagnose the causes of glacier retreat.

Fig. 1 .
Fig. 1.Model set-up.(a) Location of Helheim Glacier, with path of chosen flowline and two indicative terminus positions.Numbers indicate distance along flowline in km.(b) Comparison of measured and modelled surface velocity at end of spin-up.(c) Comparison of observed surface profile with relaxed surface profile used to initialize model.

Fig. 2 .
Fig. 2. Model output for 5 yr run with no seasonal inputs.(a) Histogram of calving event sizes (b) Comparison of modelled and observed terminus evolution.

Fig. 3 .
Fig. 3. Results of crevasse water depth (cwd) experiments.(a) Longitudinal deviatoric stress near terminus (first timestep), showing the base of the crevasse field (white lines) with differing depths of water in crevasses from 0 to 50 m.Lower profiles correspond to greater cwd.(b) Terminus position over the 5 yr model run, grey shading indicates winter periods where there is no water in crevasses, white shading indicates summer where a range in cwd was applied.For clarity the terminus position shown is measured after each calving event.Experiments were terminated where the model reached the flotation point.

Fig. 4 .
Fig. 4. Results of basal water pressure experiments.(a) Change in longitudinal deviatoric stress when a change in hydrostatic head at the bed of 50 m is applied (first timestep).The base of the crevasse field is almost identical in both the initial and perturbed experiments.(b) Seasonally applied change in basal water pressure.(c) Resultant change in surface velocity at calving front.(d) Evolution in terminus position with fixed and seasonally varying basal water pressure.(e) Repeat of results shown in (d) but with an additional seasonal change in cwd of 30 m. Grey lines show the continuous terminus position, while darker lines show the terminus position after each calving event to highlight overall trends.

Fig. 5 .
Fig. 5. Results of submarine melt experiments.(a) Change in longitudinal deviatoric stress caused by an applied cavity of 200 m, the resulting displacement in the base of the crevasse field (solid and dashed lines) was small.(b) Terminus evolution over 5 yr model run (measured after each calving event).Grey shaded areas indicate winter season with fixed m.m.r. of 150 m a −1 , white areas indicate summer with varying m.m.r.applied.The 1000 and 2000 m a −1 lines coincide.(c) Terminus evolution repeated with an additional seasonal change in cwd of 30 m applied (all lines overlie each other).

Fig. 6 .
Fig. 6.Terminus location after a one year model run, forced with varying rates of submarine melt.(top) No water in crevasses.(bottom) Seasonally varying cwd with 0 m in winter and 30 m in summer.

Fig. 7 .
Fig. 7. Change in longitudinal deviatoric stress caused by applying a backstress of both 5.0 × 10 7 Nm −1 and 25.0 × 10 7 Nm −1 to the terminus, over a vertical range of +15 to −120 m (first timestep).Solid lines show the depth of the crevasse field in the initial state and dashed lines show crevasse depth with backstress applied.The lower backstress value causes little change in the modelled crevasse depth.

Fig. 8 .
Fig. 8. Terminus evolution of modelled glacier with varying backstress applied.(a) No water in crevasses, (b) seasonally varying cwd applied with a winter value of 0 m and summer value of 30 m.The two subfigures have different vertical scales.In (a) and (b), blue shading indicates winter when backstress is applied and there is no water in crevasses.White shading indicates summer when backstress is zero and there may be water in crevasses.In (b), yellow shading indicates autumn when there is no backstress and the crevasses are dry.In (b) all lines overlie each other, with the exception of 50.0 × 10 7 Nm −1 , which shows different terminus behaviour.

Fig. 9 .
Fig. 9. Displacement in stress field caused by undercut cavities of varying length, applied with both a fixed and evolving terminus.When the terminus is allowed to evolve, the ice displacement counterbalances the change in the stress field caused by the applied cavity.