To account for gentle surface slopes, glacier surface elevations were taken from the SPIRIT Digital Elevation Model (DEM) V1 (whose correlation parameters are set for gentle slopes), based on SPOT5 Stereoscopic Survey of Polar Ice imagery acquired on 1 September 2008. The DEM has a 40 m resolution and a 30 m root-mean-square (RMS) absolute horizontal precision (http://polardali.spotimage.fr:8092/wstools/IPY/).

Bed topography was derived from ground-penetrating radar (GPR) data (Grabiec et al., 2012; Navarro et al., 2014) and depth soundings in the glacier forebay (Vieli et al., 2002). The resulting initial geometry of the modeled glacier is shown in Figure 2.

Surface velocities were measured daily at 16 stakes (Figure 1) between May 2005 and April 2011, using differential GPS (Puczko, 2012). In this study, we focus on a subset of eight stakes close to the modeled flowline (Figure 3). We additionally use velocities of the calving front measured in 2009 by terrestrial laser scanning (TLS; data generously provided by Jacek Jania [University of Silesia] from surveying and data processing by Jacek Krawiec [Laser 3D], Artur Adamek [Warsaw University of Technology], and Jacek Jania). Since a value for the velocity at the calving front is needed at each time step during the entire simulation period, but we only have TLS measurements for 2009, we assume a constant ratio between the frontal velocities and those at the closest stake.

The near-terminus surface velocities exhibit a seasonal pattern overlaid by strong interannual variability. Each year has a spring speed up, followed by a rapid slow down, followed sometimes by a gradual speed up through the winter (Figure 3A). In 2005–2006, the highest velocities occurred in summer, while in 2007 the highest velocities occurred in November–December. In 2010, velocities remained very high through winter and spring before dropping dramatically during the summer. In 2008–2011, velocities were much higher than in 2005–2008. The mean varied considerably between years (Figure 3B), discouraging the use of a single mean velocity in the model. Similarly, the very large interannual variability discourages the use of summer and winter means for summer and winter modeled periods (Figure 3C).

Interannual observations of Hansbreen's terminus position over the last decades reveal a generally smooth retreat with occasional abrupt changes (e.g., Vieli et al., 2002). Seasonal variations of the terminus position have also been observed (Blaszczyk et al., 2013). The “weekly” (assuming a week as a 1/48 of a year) terminus positions used in this paper were derived from time-lapse photographs of the Hansbreen terminus taken ca. every 3 h from December 2009 to September 2011 by three different cameras installed at surveyed positions on the eastern slope of Fugleberget (Figure 1). The two Canon EOS 1000D cameras (each equipped with a Canon EF-S 18–55 mm lens) were calibrated from images of a grid pattern using the Camera Calibration Toolbox for Matlab (www.vision.caltech.edu/bouguetj/calib_doc/), while the Canon Powershot A530 camera (which no longer existed at the time of the analysis) had to be calibrated and oriented simultaneously from multiple images of ground control points surveyed by L. Kolondra (unpublished report, 2011). Terminus positions were mapped by tracing the waterline along the terminus in each image, then projecting those pixel coordinates onto a horizontal plane at an altitude of 0 m (tidal heights were not available) to convert them to world coordinates. Standard errors of 0.47 m in water level due to tides (Zagórski et al., 2015; Michał Ciepły, personal communication, 2012), 0.68 pixels due to uncertainties in terminus tracing, and 3.23, 5.70, and 10.53 pixels for each of the three cameras due to uncertainties in camera calibration and orientation result in standard errors for width-averaged glacier length of 3.79, 3.37, and 14.54 m, respectively.

The ice mélange in the glacier forebay was qualitatively evaluated as either “complete,” “partial,” or “free” from the same time-lapse photographs used to measure terminus position. For the remaining period, we used the values of the nearest cell in a 25-km resolution time series of sea ice concentrations derived from Nimbus-7 SMMR and DMSP SSM/I-SSMIS passive microwave data (Cavalieri, 1996, updated yearly). These values were used taking into account the partial overlap of the grid cell with Hornsund fjord mouth, and their comparable sizes (Figure 4).

We apply a dynamical downscaling method (which uses a modified version of the Polar WRF 3.4.1 model) to produce—from a regional climate dataset consisting of meteorological, sea-surface temperature and sea-ice concentration data—input data for glacier thermomechanical modeling of Hansbreen (Finkelnburg et al., in preparation).

Surface mass balance (SMB) was obtained from European Arctic Reanalysis (EAR) data, with 2 km horizontal resolution and hourly temporal resolution, constrained by automatic weather stations (one in Hornsund and two in Hansbreen) and stake observations (Finkelnburg, 2013). First, ablation was calculated from the surface energy balance (SEB), which is resolved in the EAR by an optimized version of the unified NOAA Land Surface Model (Chen and Dudhia, 2001) of the Polar WRF 3.4.1 model. The algorithm solving for the SEB takes into account net radiation, sensible heat flux, latent heat flux, and ground heat flux, and encompasses all heat fluxes involved in melt and refreezing processes within the snowpack. Second, accumulation was obtained as the solid (frozen) precipitation of the Morrison bulk microphysics scheme for cloud physics used by the EAR. Finally, monthly mean SMB and surface meltwater (SMW) at each flowline point was calculated by applying bilinear interpolation to the available 2-km resolution hourly accumulation and ablation data (Figure 4).

Ice is treated as an incompressible viscous fluid. The Stokes system of equations describing the dynamical model is composed of equations describing the steady conservation of linear momentum and the conservation of mass of an incompressible continuous medium:

where σ is the Cauchy stress tensor, u is the velocity vector, g is the gravity force vector and ρ is the density. The body force f is added to account, in our 2-dimensional (2-D) model, for the friction on the lateral side of the glacier. To this end, the concept of shape factor (Nye, 1965) is here extended to the full-Stokes formulation by defining the body force f as (Jay-Allemand et al., 2011):

where the shape factor f = f(x) is a scalar function of the transversal shape of the glacier and t is the unit vector tangent to the upper surface. Jay-Allemand et al. (2011) evaluated f(x) by assuming that the transverse shape of the bedrock is a parabola, and they found an empirical estimate of the shape factor:

where h(x) = z_s(x)−b(x) is the ice thickness (expressed as the difference between the surface and bed elevations) and w(x) is the half-width at the glacier surface.

As the constitutive relation, we adopt Nye's generalization of Glen's flow law (Glen, 1955; Nye, 1957):

This equation links the deviatoric stress τ to the strain rate ε˙. The effective viscosity η is written as

where Iε˙22 represents the square of the second invariant of the strain rate tensor, A is the softness parameter in Glen's flow law and E is an enhancement factor.

The constitutive relation (4) is expressed in terms of deviatoric stresses, while the conservation of linear momentum (1) is given in terms of full (Cauchy) stresses. Both stresses are linked through the equation

where p is the pressure (compressive mean stress). Typical values for the flow-law exponent (n = 3) and softness (A = 0.1 bar⁻³a⁻¹) were used in the model (Albrecht et al., 2000; Vieli et al., 2002). This value of softness is adequate (see e.g., Cuffey and Paterson, 2010) for a polythermal glacier like Hansbreen, composed mostly of temperate ice except for a thin upper layer of cold ice in the ablation zone (Jania et al., 1996).

We introduce a scalar damage variable D that quantifies the loss of load-bearing surface area due to fractures, known as fracture-induced softening (Borstad et al., 2012). This softening is taken into account through the introduction of the damage within Glen's law. Following Borstad et al. (2012) and Krug et al. (2014), the enhancement factor E can be linked to the damage D as

For undamaged ice (D = 0), E = 1 and the flow regime is unchanged. As damage increases (D > 0), E > 1, ice viscosity decreases, and flow velocity increases.

In this study, we use a very simple function for D. Damage is nonzero only in the lowermost 2 km of the glacier (near-terminus heavily crevassed area) and increases linearly toward the terminus, where it reaches a maximum value of 0.4 (Krug et al., 2015). This value provided a good fit to observed velocities in preliminary experiments (not shown here) of the sensitivity of the modeled velocities to changes in damage.

The time evolution of the glacier surface is calculated by solving the free-surface evolution equation

where z_S is the surface elevation, t is time, u_S and v_S are the horizontal and vertical components of the flow velocity at the surface, respectively, and b is the surface mass balance.

The upper surface of the glacier is a traction-free zone with unconstrained velocities. At the ice divide at the head of the glacier, horizontal velocity and shear stresses are set to zero.

For boundary conditions at the bed, we use a friction law that relates the sliding velocity to the basal shear stress in such a way that the latter is not set as an external condition but part of the solution:

where C, the friction coefficient, is determined using the inverse Robin method proposed by Arthern and Gudmundsson (2010) and modified by Jay-Allemand et al. (2011). The latter study includes a regularization parameter, λ, for which we have adopted a value of 0.4 determined from preliminary tests. The inverse method infers the basal friction parameter by reducing the mismatch between observed and modeled surface velocities using a cost function.

Since the inversion procedure requires a continuous function for the surface velocity, we calculate it as a sixth-degree polynomial regression for each 4-week period.

At the glacier terminus, we set backstress to zero above sea level and equal to the water-depth-dependent hydrostatic pressure below sea level. In model runs with ice mélange, the additional backstress is applied to the calving face, in the opposite direction of ice flow (negative X). In the absence of further data, we assumed the range of 30–60 kPa estimated by Walter et al. (2012) for Store Gletscher as indicative of the order of magnitude that we could expect for Hansbreen, despite their different settings. While Store Gletscher is buttressed by a rigid proglacial mélange, with a thickness reaching 75 m (Todd and Christoffersen, 2014), Hansbreen presents a thinner layer of ice mélange made up of a mixture of growlers, bergy bits and small icebergs bonded by sea ice. As no measurements of ice mélange thickness in Hansbreen forebay are available, but the sea-ice maximum thickness in Hornsund fjord is known to be around 1 m (Kruszewski, 2012), we assume a mean ice mélange freeboard height of 0.5 m and a mean thickness of 4.5 m.

The CDw calving criterion (Benn et al., 2007) assumes that calving is triggered by the downward propagation of transverse surface crevasses near the calving front as a result of the extensional stress regime. Following Nye (1957), crevasse depth is calculated as the depth where the longitudinal tensile strain rate tending to open the crevasse equals the creep closure resulting from the ice overburden pressure. This procedure incorporates the full stress solution into the crevasse depth criterion. In Benn's model, calving is assumed to occur when surface crevasses reach the waterline.

Following Todd and Christoffersen (2014), crevasse depths are calculated from the balance of forces:

where σ_n, the “net stress,” is positive for extension and negative for compression. The first term on the right-hand side of Equation (10) represents the opening force of longitudinal stretching, adapted by Todd and Christoffersen (2014) from Otero et al. (2010); τ_e represents the effective stress, τe2=τxx2+τzx2. τ_e is multiplied by the sign function of the longitudinal deviatoric stress, τ_xx, to ensure that crevasse opening is only produced under longitudinal extension (τ_xx > 0). The second term on the right-hand side is the ice overburden pressure, which leads to creep closure, where ρ_i is the density of glacier ice, g is the acceleration of gravity and d is the crevasse depth. The last term represents the water pressure which contributes to open the crevasse, which is a function of the depth of water filling the crevasse.

At each time step, the glacier is divided into a rectangular mesh with 10 vertical layers and a horizontal grid size of ca. 50 m in the upper glacier and ca. 25 m near the terminus. The Stokes system of equations (1) is solved by a finite element method using Elmer/Ice and the 2-D stress and velocity fields are computed along the central flowline (Figure 1). The new surface elevations are computed from the surface mass-balance input and the solved surface velocities using the free-surface evolution equation and the grid nodes are shifted vertically to fit the new geometry.

At the terminus, the grid nodes are shifted down-glacier according to the velocity vector and the length of the time step and the terminus position is updated according to the calving criterion.

Prognostic model runs were carried out with a 1-week (1/48 of a year) time step, starting from the 2008 glacier geometry. Every 4 weeks (four time steps), we ran an initialization process which consisted of solving the Robin problem (Jay-Allemand et al., 2011) to force a best-fit friction coefficient to be used for those four model runs. This forcing was done to minimize the misfit between the observed and modeled velocities. The choice of the initialization time step was made as a compromise between the time resolution needed for capturing the sudden changes in velocity and an acceptable computational cost.

Given the difficulties of measuring the depth of water in crevasses, we ran the model for a range of crevasse water depths (from 6 to 12 m) to evaluate the sensitivity of the model to this parameter. We found that calving rate is highly dependent on the depth of water in crevasses, with an increase of just a few meters causing the glacier to switch from advance to retreat (Figure 5). If a constant water depth is prescribed, the model is unable to reproduce the observed terminus position fluctuations (Figure 5), although the results allow us to select the water depth which, on average, best fits the observations (10 m, as illustrated by Figure 5).

To test the effect of ice mélange on calving rate and terminus position, we varied ice mélange backstress from 0 to 70 kPa (based on Walter et al., 2012, as discussed in the Introduction) in steps of 10 kPa and multiplied by the fraction of ice mélange coverage (weekly average) in the glacier forebay, varying from 0 when no ice mélange is present to 1 when ice mélange completely fills the glacier forebay. When there is only partial coverage of ice mélange, ice mélange remains concentrated near the margins of the glacier and the fjord and therefore continues to exert some backstress on the glacier front. This supports the prescription of backstress even at low ice mélange concentrations. In this experiment, we prescribed a constant crevasse water depth of 10 m, the best-fit value from the previous experiment.

We found that the effect of ice mélange backstress on glacier front position was significant, even under low stresses (Figure 6), and that 50 kPa (for full mélange coverage), applied to the front of the glacier, yielded the best fit to observations during the winter when ice mélange was present. The effect of a backstress of 50 kPa on the modeled longitudinal deviatoric stress profile is shown in Figure 7.

As discussed by Todd and Christoffersen (2014), the relationship between surface melt rate and crevasse water depth depends on many factors, including the distribution, shape, and depth of crevasses, the melting and refreezing on crevasse walls, and the potential drainage of water from crevasses into englacial, subglacial, and proglacial water bodies. Since observations of these processes are very scarce and the water in crevasses starts as surface meltwater, we chose to parameterize crevasse water depth in terms of the available surface meltwater predicted by the SEB model (Figure 8).

Use of such a parameterized time-varying crevasse water depth, in combination with the best-fit ice mélange backstress from the previous experiment, yielded terminus positions in very good agreement with the observations (Figure 8), reproducing the winter advance and the subsequent summer retreat.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.