2.1. Ice-cap geometry and geometric changes

Surface elevation and ice thickness, and thereby bed topography, were extensively mapped by airborne radio-echo sounding (RES) in 1983 (Reference Dowdeswell, Drewry, Cooper, Gorman, Liestøl and OrheimDowdeswell and others, 1986). Austfonna consists of a main central dome merging with a somewhat smaller dome to its south that together feed a number of drainage basins (Figs 1 and 2a). A main ice divide oriented southwest/northeast forms a natural border between the northwestern basins, predominately terminating on land or in fjords, and the southeastern basins that are, to a large extent, grounded below sea level and form a continuous calving front towards the Barents Sea. Several of the drainage basins are known to have surged in the past, specifically Etonbreen (1938 or earlier) in the west, and Bråsvellbreen (1937–38) and Basin-3 (between 1850 and 1873; Fig. 1) in the southeast. Observations do not span one entire surge cycle of any of these basins. Based on the approximate mass transfer from the reservoir into the receiving zone, the duration of the quiescent phase has been estimated as ∼150–500 years (Reference SchyttSchytt, 1969; Reference SolheimSolheim, 1991), considerably longer than for smaller glaciers in western Svalbard or other areas in general (Reference Dowdeswell, Hamilton and HagenDowdeswell and others, 1991). Over 28% of the entire area covered by Austfonna lies below sea level (Reference DowdeswellDowdeswell, 1986), while for individual basins, headed by Bråsvellbreen, the proportion is as large as 57%. The lowest bedrock elevations, >150 m below sea level, are found in depressions below the lower reaches of Basin-3 and Leighbreen (Fig. 2). The calving margins of the marine grounded ice add up to a total length of ∼230 km (Reference Dowdeswell, Benham, Strozzi and HagenDowdeswell and others, 2008).

Fig. 2. Maps of (a) the bedrock topography and (b) modern observed ice thickness. The white solid line in (a) indicates the modern observed glacier outline. Contours are spaced at 50 m.

The ice cap is undergoing surface elevation changes that are potentially attributed to a build-up during a quiescent phase with interior thickening of up to 0.5 m a⁻¹ and marginal thinning of 1–3 m a⁻¹ (Reference Moholdt, Hagen, Eiken and SchulerMoholdt and others, 2010). Marginal thinning encompasses the stagnant receiving zones of the three known surge-type basins, which are characterized by depleted ice-surface profiles and low driving stresses, both indicative of glaciers in an early quiescent phase (Reference DowdeswellDowdeswell, 1986). Basins with steep ice-surface profiles and relatively high marginal shear stresses are interpreted to be frozen to their beds and may also be surge-type, particularly if they contain ice regions grounded below sea level (Reference DowdeswellDowdeswell, 1986). Spaceborne radar interferometry of Austfonna from the mid-1990s revealed a velocity structure of a slow-moving polar ice cap (<10 m a⁻¹), interrupted by distinct fast-flow units having velocities in excess of 100 m a⁻¹ (Reference Dowdeswell, Unwin, Nuttall and WinghamDowdeswell and others, 1999). Reference Bevan, Luckman, Murray, Sykes and KohlerBevan and others (2007) suggested that slow ice motion is the key factor for the interior thickening, since the actual ice flux across the equilibrium line appeared to be only half the balance flux. Changes in the accumulation/ablation pattern have been suggested as an alternative mechanism for the observed geometric changes (Reference Bamber, Krabill, Raper and DowdeswellBamber and others, 2004). This is not supported by observations of surface mass balance (SMB). In situ data from shallow ice cores do not indicate a significant trend in the net SMB of the accumulation area for the periods 1963–86 and 1986–98/99 (Reference Pinglot, Hagen, Melvold, Eiken and VincentPinglot and others, 2001), although large interannual variations in the total amount of snow accumulation have been observed in recent years (Reference TaurisanoTaurisano and others, 2007; Reference Dunse, Schuler, Hagen, Eiken, Brandt and HøgdaDunse and others, 2009). These observations also indicate a clear asymmetry in snow accumulation that is reflected in the net SMB (Reference Schuler, Loe, Taurisano, Eiken, Hagen and KohlerSchuler and others, 2007) and, hence, the ELA: ∼450 m a.s.l. in the northwest and 300 m a.s.l. in the southeast (Reference Moholdt, Hagen, Eiken and SchulerMoholdt and others, 2010). Geodetically derived mass balance over the period 2002–08 suggests a net SMB close to zero (Reference Moholdt, Hagen, Eiken and SchulerMoholdt and others, 2010). However, the net mass balance of Austfonna is negative, −1.3 ± 0.5 km³ a⁻¹, due to significant retreat of the marine margins (Reference Dowdeswell, Benham, Strozzi and HagenDowdeswell and others, 2008; Reference Moholdt, Hagen, Eiken and SchulerMoholdt and others, 2010).

2.2. Thermal regime

Austfonna’s thermal structure is described as polythermal (Reference Dowdeswell, Drewry, Cooper, Gorman, Liestøl and OrheimDowdeswell and others, 1986), although no direct evidence of a considerable temperate ice volume is available. In the firn area, latent heat release from internal refreezing of surface meltwater contributes significantly to warming the near-surface firn and ice. Current near-surface (10 m) temperatures in the firn area exceed those in the ablation area, despite lower air temperatures (e.g. Reference Zagorodnov, Sin’kevich and ArkhipovZagorodnov and others, 1989a; and as evident in our thermistor measurements in the ablation area of Etonbreen and at the summit from 2006 to 10; T.V. Schuler, unpublished data). Direct measurements of the temperature distribution throughout an ice column are available from a borehole in the summit area, drilled in 1987 (Reference Zagorodnov, Sin’kevich and ArkhipovZagorodnov and others, 1989a). The vertical temperature profile shows that the bulk of the 567 m thick ice column is cold, underlain by an isothermal basal layer ∼30 m thick that is close to the pmp. A temperature minimum of −8°C at ∼160 m depth at the drilling site is associated with lower temperatures during the Little Ice Age, when firn warming was presumably insignificant (Reference Zagorodnov, Sin’kevich and ArkhipovZagorodnov and others, 1989b). The temperature gradient in the lower 400 m suggests a moderate geothermal heat flux of ∼40 mW m⁻² (Reference Ignatieva and MacheretIgnatieva and Macheret, 1991). Information on the thermal structure of the ice cap is also contained in the airborne RES data (Reference Macheret, Bobrova and SankinaMacheret and others, 1991). An internal reflection horizon (IRH), separating a transparent zone above and a zone of high backscatter below, may be associated with the cold/temperate transition surface (CTS). Such IRHs are a characteristic of many glaciers in Spitsbergen/western Svalbard (Reference BjörnssonBjörnsson and others, 1996). The Institute of Geography, USSR Academy of Sciences, collected such data in 1984 (Reference Macheret and VasilenkoMacheret and Vasilenko, 1988). Similar to the RES data of Reference Dowdeswell, Drewry, Cooper, Gorman, Liestøl and OrheimDowdeswell and others (1986), the absence of a continuous IRH indicates that significant volumes of temperate ice do not exist. Reference Kotlyakov and MacheretKotlyakov and Macheret (1987) and Reference Macheret and VasilenkoMacheret and Vasilenko (1988) noted that the absence of continuous IRHs does not exclude the presence of a thin bottom layer of temperate ice, nor the occurrence of water at the ice/bed interface, because a thin layer of temperate ice may not be resolved by the radar system. The vanishing of bedrock returns in some regions (e.g. Leighbreen) may indicate the presence of a localized basal temperate layer. A simple theoretical consideration leads to a similar result: given a constant temperature gradient determined by the geothermal heat flux, an ice thickness more than 500 m, which is about the maximum ice thickness of the ice cap, is required to approach the pmp when internal strain heating is neglected. Without a significant contribution from strain heating or latent heat release by infiltrating surface meltwater, the bulk of the ice cap is expected to be cold, but potentially underlain by a temperate base.

2.3. Basal properties

Turbid meltwater plumes arising from several points beneath the marine-ice margin of Austfonna provide evidence that the marine grounded areas are, at least to some extent, underlain by soft sediments and subjected to basal water flow (Reference Macheret and VasilenkoMacheret and Vasilenko, 1988; Reference Dowdeswell, Unwin, Nuttall and WinghamDowdeswell and others, 1999). Furthermore, Reference Solheim and PfirmanSolheim and Pfirman (1985) and Reference SolheimSolheim (1986) mapped continuous submarine sediment ridges in front of Bråsvellbreen and Basin-3 and conclude that these ridges were partly pushed up by the advancing surge front and partly deposited by vigorous meltwater drainage along the entire surge front during the period of maximum extension. A rhombohedral pattern of linear, discontinuous sediment ridges and mounds is interpreted as being formed by a squeeze-up process, filling bottom crevasses that opened during surge advance. After surge termination, the surge lobe was stagnant and subglacial features have been preserved. These subglacial sediments may originate from glacial erosion of the underlying bedrock and/or marine sedimentation during periods when the glacier extent over the seabed was considerably smaller then today (e.g. the Holocene Optimum ∼9 to 5 ka BP (Reference Hjort, Mangerud, Adrielsson, Bondevik, Landvik and SalvigsenHjort and others, 1995; Reference BlakeBlake, 2006)). A geological map of Nordaustlandet published by the Norwegian Polar Institute (NPI) indicates that an east–west line along Wahlenbergfjorden separates hard bedrock types in the north, i.e. Late Precambrian Hecla Hoek rocks and Caledonian granite/ migmatite complexes, from soft bedrock types in the south, i.e. flat lying young sediments of Middle Carboniferous to Lower Cretaceous age (Reference Lauritzen and OhtaLauritzen and Ohta, 1984). Information on bedrock type for presently glacierized areas is not available, but it is noteworthy that the known surge-type basins, Etonbreen and Bråsvellbreen, are located adjacent to the soft, easily erodible bedrock types.

3.1. Basal sliding

The dynamic boundary condition at the glacier bed is given by a Weertman-type sliding law (Reference Greve and BlatterGreve and Blatter, 2009). Typically, sliding is enabled for a temperate ice base and no-slip conditions prevail for a cold base frozen to the ground. Here we use a modified expression of the basal sliding velocity, v _b, following Reference Hindmarsh and Le MeurHindmarsh and Le Meur (2001) and Reference GreveGreve (2005), that allows for partial sliding as the basal temperature approaches the pmp, termed ‘submelt sliding’:

(1)

where τ _b is the basal shear stress, N _b the basal pressure, C _b, p and q are the sliding-law coefficient and exponents and T’_b the basal temperature relative to pmp. Full sliding is enabled for a temperate ice base and exponentially decreases when basal temperatures drop below the pmp. The transition between no-slip and slip conditions at the glacier base is prescribed by a submelt-sliding parameter, γ (°C). It defines the temperature deviation from pmp for which full sliding is reduced by a factor e⁻¹. A large positive value of γ allows for initiation of sliding well below the pmp in a smooth way. A small positive value leads to abrupt onset of sliding close to the pmp. The basal pressure is given either by the full hydrostatic ice-overburden pressure, P _f,

(2)

(where H _i is the thickness of the ice), or by the reduced basal pressure, P _r, the difference between the ice-overburden pressure and the basal water pressure. To account for the buoyancy force exerted upon the ice volume submerged below sea level, the basal water pressure is assumed to equal the sea-water pressure,

(3)

where ρ _sw is the density of the sea water and D _w the water depth. Note that drainage of basal water from beneath the marine terminus requires an additional pressure head in order to overcome the sea-water pressure (Reference PfefferPfeffer, 2007). To assure numerical stability and to keep sliding velocities given by Equation (1) within a reasonable range, P _r is constrained to be ≥0.2 × P _f.

3.2. Calving

Grounded marine-ice margins typically constitute a vertical calving front that fulfils a certain stability criterion, i.e. the actual ice thickness must exceed the flotation thickness at all times (H _i ≥ ρ _sw ρ ⁻¹ D _w). Otherwise the grounded margin will either disintegrate or form a floating ice tongue or ice shelf. Models using a fixed grid in the map plane usually lack the ability to track the position of the calving front on a subgrid scale. Nevertheless, the models can simulate a retreat of a marine ice margin by setting H _i at particular gridpoints to zero. Advances of the calving front, however, require the marine-ice margin to protrude at least one gridcell in one single time-step while overcoming the stability criterion. This cannot be realized at reasonable flow velocities. To bypass this shortcoming, models usually allow the ice to extend underwater without becoming afloat. Here this ‘underwater ice’ is treated as regular ice, but is further subjected to calving that acts on its surface as an additional ablation term. We employ the modified version of the calving parameterization of Reference Marshall, Tarasov, Clarke and PeltierMarshall and others (2000) that relates the calving flux, Q _c, to the water depth, D _w, to some power k and the ice thickness, H _i, to some power l:

(4)

where k _c is a calving parameter that may include effects of margin geometry, longitudinal stress gradients, ice temperature or hydraulic fracturing. Here we keep it simple and set k _c to 10⁻⁴ m⁻¹ a⁻¹, and k = l = 1. These parameter choices enable retreat and advance cycles of marine grounded ice margins within a bathymetric range comparable with modern observed values. Increasing the calving parameter to k _c = 10⁻³ m⁻¹ a⁻¹ leads to retreat of the marine grounded basins to shallow waters, as are characteristic for the fjord heads in the northwest, while the deeper southeastern areas become completely and permanently deglacierized.

This treatment of marine ice represents an improvement over earlier approaches in large-scale ice-sheet models that limit the marine-ice extent by defining a bathymetric contour at which the ice is simply cut off, either described by a constant value (Reference Tarasov and PeltierTarasov and Peltier, 1997) or a function of sea level (Reference Zweck and HuybrechtsZweck and Huybrechts, 2005). A new concept of subgrid parameterization of ice-front motion that maintains a vertical calving front has recently been applied to another large-scale finite-difference ice-sheet model (Reference Albrecht, Martin, Haseloff, Winkelmann and LevermannAlbrecht and others, 2010). This approach provides a promising basis for implementation of physically based calving laws, such as the waterline crevasse-depth model of Reference Benn, Hulton and MottramBenn and others (2007).

3.3. Input data and initialization

The required model input consists of surface and bedrock topography, providing the initial ice-cap geometry, the climate forcing in terms of precipitation and air-temperature fields, sea-level stand and geothermal heat flux.

The surface elevations above sea level are based on a 1 :250 000 map of Svalbard (Sheet 3) published by NPI in 1998. The bathymetry is represented by the publicly available International Bathymetric Chart of the Arctic Ocean, Version 2.0 (Reference JakobssonJakobsson and others, 2008) having a horizontal resolution of 2 km (Figs 1 and 2a). The maximum ice thickness of ∼600 m coincides with the central dome (Fig. 2b). The bedrock beneath the ice cap is derived by subtracting scattered ice-thickness data from the surface topography. We supplement the airborne RES data published by Reference Dowdeswell, Drewry, Cooper, Gorman, Liestøl and OrheimDowdeswell and others (1986) with two new datasets from spring 2008. The Danish National Space Centre acquired 60 MHz airborne RES data, while a ground team collected 20 MHz ground-penetrating radar data along profiles ∼800 km long (Reference Vasilenko, Navarro, Dunse, Eiken and HagenVasilenko and others, 2010). The profiles follow the routes investigated by the University of Oslo and NPI for mass-balance purposes since 1998 (Reference Schuler, Loe, Taurisano, Eiken, Hagen and KohlerSchuler and others, 2007; Reference TaurisanoTaurisano and others, 2007; Reference Dunse, Schuler, Hagen, Eiken, Brandt and HøgdaDunse and others, 2009) and cover special areas of interest, such as Duvebreen and Basin-3.

We have derived an idealized present-day climate based on field data acquired since 1998 and European Centre for Medium-Range Weather Forecasts (ECMWF) ERA-40 and ERA-Interim Re-analysis. The temperature input is based on a 5 year record (2004–09) collected by an automatic weather station on Etonbreen at 510 m a.s.l., 50–100 m above the recent mean ELA. Monthly mean temperatures are distributed over the entire model domain using a mean lapse rate of −0.0045 K m⁻¹ (Reference Schuler, Loe, Taurisano, Eiken, Hagen and KohlerSchuler and others, 2007). The precipitation field is based on ERA-40/ERA-Interim Reanalysis at a 6 hour temporal resolution for the period 1960–2006. Monthly mean precipitation fields have been downscaled using a simple precipitation model (Reference Smith and BarstadSmith and Barstad, 2004; Reference Schuler, Crochet, Hock, Jackson, Barstad and JohannessonSchuler and others, 2008). The precipitation fields are kept constant over the entire model period while we allow the surface temperature to adjust for changes in ice-surface elevation given the constant lapse rate. Application of this idealized present-day climate on the present-day topography leads to a net SMB close to zero (+0.07 m w.e. a⁻¹),in agreement with recent observations for the period 2002–08 (Reference Moholdt, Hagen, Eiken and SchulerMoholdt and others, 2010).

The vertical temperature profile through the ice at the summit presented by Reference Zagorodnov, Sin’kevich and ArkhipovZagorodnov and others (1989b) indicates a geothermal heat flux of ∼40 mW m⁻² (Reference Ignatieva and MacheretIgnatieva and Macheret, 1991). The ice temperature for the initial state assumes a steady-state temperature gradient according to the geothermal heat flux, while the annual mean air temperature constrains the temperature at the ice surface. The temperature within the ice volume is limited to the basal pmp determined by the ice-overburden pressure. At the summit, good agreement of the initial temperature profile with the observed one is achieved for the lower 400 m. In the upper layer, the initial ice temperatures are too low. However, after ∼100 model years, the firn-warming routine of SICOPOLIS has established a temperature profile in the summit area similar to the observed one, throughout the ice column.

3.4. Model experiments

We define a series of model experiments (experiments 1–5), each characterized by a unique combination of the sliding-law parameters C _b, p, q and N _b. The individual experiments represent a broad range of potential basal environments, facilitating basal motion to a varying degree, from hard rock or completely immobile sediment/permafrost to soft water-saturated sediments. Table 2 lists relevant parameters used in all model experiments. Experiment 1 considers sliding over hard rock everywhere with basal pressure given by the hydrostatic ice-overburden pressure. Experiment 2 also considers sliding over hard rock everywhere, but reduced basal pressure applied for the marine grounded parts of the ice cap. Experiments 3 and 4 differ from experiments 1 and 2 by a threefold enhancement of basal sliding of the marine grounded ice. Experiment 5 applies a linear sliding law to the marine grounded ice and represents fast basal motion by plastic deformation of water-saturated sediments (Hind-marsh, 1997; Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000). Each experiment was run four times using different values of γ: 0.05, 0.2, 0.5 and 1.0°C, describing a progressively smoother transition from no-slip to full basal sliding.

Table 2. Sliding-law parameter combinations of the five model experiments. The sliding-law coefficient, C _b(C _r for hard rock and C _s for soft sediments as specified in Table 1), and the exponents p and q are first given for ice areas resting on land, then for the ice areas grounded below sea level. N _b is either the full (P _f) or reduced basal pressure (P _r). Note that in the case of a linear sliding law (marine grounded ice in experiment 5) basal sliding is independent of the basal pressure (exponent q = 0)

All model experiments take the present-day geometry of Austfonna as the initial condition. The present-day reference climate is applied as external forcing. We run the model for 30 ka and evaluate the last 5 ka, when all runs have reached a state with negligible long-term trends, i.e. the temporal evolution of a particular variable is described by a stationary signal.

4.1. Mean ice-cap size and variability

Total ice volume, areal extent and basal area at the pmp are presented in Figure 3. The left panels (a, b, c) give mean values while the associated variability is shown to the right (d, e, f). The temporal variability is expressed in terms of two standard deviations normalized by their respective mean value to capture the amplitudes of the stationary signal (∼95% of the respective output variable is then contained, assuming a normal distribution). Both progressively increasing enhancement of basal sliding (experiments 1–5) and an increasing value of γ lead to a smaller ice cap in terms of volume and areal extent, a consequence of increasing mass transport from the interior of the ice cap towards the margins. The dashed lines (Fig. 3a and b) indicate a ±10% region around the current ice volume and areal extent, marked by solid lines. Assuming that Austfonna is close to equilibrium, model runs using γ = 0.2°C generally provide the best fit of simulated to current ice-cap size. For experiment 5, best fit is achieved using γ = 0:05°C. The simulated temperate basal area is largest for experiment 1, sliding over hard rock, and smallest for experiment 5, deformation of soft sediments, coinciding with the largest and smallest ice volume, thickest and thinnest ice, respectively.

Fig. 3. Comparison of the five model experiments with varied submelt sliding parameter, γ. Mean values of (a) total volume, (b) total glacierized area and (c) percentage of temperate basal area relative to the total area for the time period 25–30 ka and their associated twofold standard deviation normalized with respect to the mean (d, e, f). The grey dashed lines envelope model results within 10% of the modern observed values (a, b) marked by a solid grey line. The modern basal area at the pmp (dashed line in c) follows from theoretical considerations using ice-cap geometry and geothermal heat flux.

Concerning the temporal variability, an abrupt onset of basal sliding, represented by a small value of γ, leads to large fluctuations in the presented quantities (Fig.3d–f). Comparing the different experiments for a given γ, experiment 5, the experiment with most enhanced basal motion, yields largest variability. The fluctuations rapidly decrease for progressively smoother onset of basal motion, and for γ = 1°C the variability has vanished in all experiments. For experiment 5, significant variability also exists for γ = 0.2 and 0.5°C. In terms of normalized standard deviations, the fluctuations in areal extent are up to two times larger than those in volume, indicating that mass redistribution dominates over mass changes. This redistribution is effected by an ice flux far in excess of the balance flux. Fluctuations in basal areal extent at the pmp are even larger than fluctuations in total areal extent, as high as 32.4% in the case of experiment 5 with γ = 0:05°C.

Figure 4 demonstrates the effect of abruptness in sliding onset (described by different values of γ) on the temporal evolution of the ice cap by means of time series of the areal extent, for two values of γ. Figure 4a presents results for model runs using a smooth transition (γ = 1.0°C), and Figure 4b presents those for an abrupt transition from no sliding to full sliding (γ = 0:2°C). When γ = 1:0°C (Fig. 4a) the areal extent is constant over time and lacks temporal variability, whereas for γ = 0.2°C (Fig. 4b), temporal variability in areal extent exists for all experiments. For experiments 1–4 fluctuations are clustered between 0.7% (experiment 3) and 1.2% (experiment 4), while for experiment 5 they are significantly higher, at ∼3.5% (Fig. 3e).

Fig. 4. Time series of total glacierized area for a quasi-steady-state time period of model years 25–30 ka for all model experiments. Model runs with (a) γ set to 1°C and (b) γ set to 0.2°C.

4.2. Permanent fast flow

Based on the absence of temporal variability in scalar variables, experiment 1 with γ = 1:0°C is classified as permanent fast flow. Figure 5 shows the basal velocity field and the basal temperature relative to pmp in the steady state. This experiment did not yield temporal fluctuations, hence only one arbitrary time slice is shown. Basal sliding dominates the overall flow everywhere and the surface velocity field (not shown) looks very similar. The bulk of the ice cap is slow-moving (<10 m a⁻¹), interrupted by several fast-flow units with moderate velocities of 50–150 m a⁻¹ . The flow units coincide with the warmest basal areas, but sliding is already initiated at basal temperatures ∼2°C below the pmp. The flow units also coincide with valleys and troughs of the bed topography where ice thickness is greatest, facilitating enhanced basal shear and basal temperatures.

Fig. 5. (a) Steady-state basal velocity and (b) basal temperature relative to pmp for experiment 1 with γ set to 1°C. The modelled glacier outline is indicated by the grey solid curve.

An intriguing result is the absence of significant volumes of temperate ice, indicating that regional fast flow is accomplished exclusively by basal motion in regions where the glacier base is near the pmp. This result agrees with the observations of the thermal structure of Austfonna.

4.3. Oscillatory fast flow

Experiment 5 reveals the largest variability in ice-cap geometry, suggesting that spatial and temporal variability exist in ice-cap dynamics. Using γ = 0:2°C leads to a variability similar in magnitude to that found using the extreme case of γ = 0.05°C. We discuss the case of γ = 0.2°C as this model run illustrates oscillatory fast flow including cyclic surge behaviour very well (Fig. 6).

Fig. 6. Snapshots of the basal velocity (left panels) and basal temperature relative to pmp (right panels) at indicated model years for experiment 5 with γ set to 0.2°C.

Some drainage basins (e.g. Leighbreen and Basin-3 (Fig. 1)) are active at all given instances, but with variable flow intensity. Basin-3 seems capable of generating large-scale surges (Fig. 6d; 29.67 ka). Other basins are only active for short periods of time when they show very high basal velocities (>500 m a⁻¹ up to several km a⁻¹), otherwise they have very low activity (e.g. Etonbreen (Fig. 6a; 29.2 ka) and Basin-5 (Fig. 6b; 29.43 ka)). As in the runs that produce permanent fast flow (section 4.2), no significant volumes of temperate ice exist. Fast flow is now strictly confined to the basal areas at, or very close to, the pmp, and large variations in basal temperature are observed for the surging basins. The active phase is associated with a bed at, or close to, the pmp, and entails massive mass fluxes along the flowlines, whereas basal temperatures are below the pmp during the quiescent phase. Figure 6b and f represent a time when the ice cap is characterized by increased fast flow and surges (29.43 ka) before the activity reduces again and the surge lobe of Basin-5 has become stagnant (Fig. 6c and g; 29.5 ka).

To investigate the behaviour of individual drainage basins and the periodicity of their oscillations, we present temporal ice-thickness variations at specific locations, lined up in 2 km intervals along the modern ELA, i.e. the elevation contour of 450 m a.s.l. in the northwest and 300 m a.s.l. in the southeast. These points are fixed in space and change their relative position in the course of the simulation, with respect to both the simulated ELA and DBL. Groups of 7–12 specific points are selected to represent a particular basin. Individual points within the group may show different behaviour (e.g. some points may become ice-free, periodically, and only glacierized in the case of a considerable advance of the terminus during a large-scale surge event). Temporal variations in the ice thickness of individual basins differ in both period and amplitude (Fig. 7). Etonbreen and Basin-5 show very marked oscillations of regular frequency and large amplitudes of up to ∼100 m. Etonbreen shows gradual thickening followed by abrupt thinning, whereas Basin-5 is characterized by abrupt thickening followed by gradual thinning. The contrasting behaviour of periodic thickening and thinning can be explained by the migration of the simulated ice divide towards the northwest, moving the specific points from their location along the modern ELA well into the reservoir zone (Etonbreen) or receiving zone (Basin-5). Leighbreen, Basin-3 and Bråsvellbreen show fluctuations of smaller amplitude and at higher frequencies, while Duvebreen shows very little change in time. Spectral analysis of the individual time series reveals principal periods in the range 200–500 years for all basins, in agreement with the estimated periods of large-scale surge events (150– 500 years). While Etonbreen and Basin-5 show typical surge-type behaviour with short-lived surge events following long quiescent phases, Leighbreen, Basin-3 and Bråsvellbreen may also be surge-type, but they maintain spatially limited fast flow with minor oscillations on a decadal timescale between large-scale surge events.

Fig. 7. Mean ice-thickness evolution at specific points fixed in space on individual basins during experiment 5 with γ set to 0.2°C. The dashed green lines indicate the times represented in Figure 6.

5.1. Through a surge cycle

In agreement with field observations and theory, fast flow is achieved primarily by basal motion. Abrupt flow oscillations are thus governed by temporal variations in basal motion and closely related to the basal thermal regime. Reference BjörnssonBjörnsson and others (1996) have shown that the basal thermal regime of Svalbard glaciers is closely linked to ice thickness.

During the quiescent phase of a surge-type basin, the basal temperature is below pmp over a large fraction of the bed. The reservoir zone experiences substantial thickening, while the receiving zone is subjected to severe thinning and retreat. The increased ice thickness leads to increased insulation from the cold atmosphere and increased overburden pressure, the latter resulting in increasing basal shear and, hence, dissipative strain heating. The bed eventually approaches the pmp in the interior of the ice cap. Sliding sets in according to the chosen γ and initiates the active surge phase. Friction generated by basal motion generates additional heat. Compressive flow in a boundary zone where fast-flowing ice (basal motion enabled) meets slowly flowing ice (basal motion disabled or highly reduced) leads to local thickening and surface steepening, provoking local enhancement of basal shear stress and strain heating. This mechanism was described as a trigger of flow instabilities by Reference PaynePayne (1995). The warm-based area expands towards the margins and provides a mechanism for forward propagation of the surge front. The advancing surge lobe is subjected to dynamic thinning and surface flattening, thereby reducing both thermal insulation and generation of heat, and the area affected by fast basal motion shrinks back upstream. Once the reservoir zone is tapped and the surface profile along the flowline flattened sufficiently, conditions for a temperate bed can no longer be maintained, basal motion becomes insignificant and brings the surge to a halt. Another quiescent phase is initiated, characterized by insufficient ice flux from the accumulation towards the ablation area and hence a gradual build-up of a reservoir zone. These distinct phases in the simulated surge cycle agree with those reported by Reference Calov, Ganopolski, Petoukhov, Claussen and GreveCalov and others (2002) and agree with observations of polythermal surge-type glaciers in other parts of Svalbard (Reference Sund, Eiken, Hagen and KääbSund and others, 2009). The simulated surge advances of ∼10–15 km over widths of ∼15–25 km (Fig. 6) are of similar size to the 1930s surge of Bråsvellbreen, during which the 30 km wide front advanced by ∼20 km (Reference SchyttSchytt, 1969).

To describe the evolution of the surge cycles in detail, we extract time-series data at one specific location on Etonbreen. Given the parameter combinations defined by experiment 5, the temporal evolution of ice thickness, basal temperature relative to pmp and relative contribution of basal motion to the overall ice flow at this location is clearly influenced by periodic surge behaviour (Fig. 8). The contribution of basal sliding vanishes if the basal temperature relative to pmp falls below −1°C (Fig. 8b and c). Basal motion becomes the dominant mechanism of glacier flow, i.e. >70% of the total (surface) motion, if basal temperatures are within 0.7°C of the pmp (experiments 1–4). The point at which basal motion becomes the principal flow mechanism does not coincide with a constant basal temperature, but rather falls within a critical temperature range that allows for large variations in basal motion (experiment 5). This indicates that additional factors (e.g. the basal shear stress) determine the exact onset of fast basal sliding that initiates the surge phase. The sawtooth-shaped signal shows cycles of gradual thickening over 250–300 years followed by abrupt thinning over <50 years. Every ice-thickness peak precedes a corresponding peak in basal temperature (Fig. 8a, experiment 5). No clear relation exists between the basal temperature and the phase of gradual thickening (quiescent phase); however, the mean ice thickness allows maintenance of basal temperatures within the critical temperature range. The basal shear stress increases with increasing ice thickness and increasing surface slope (not shown in Fig. 8) and eventually reaches a level at which basal sliding becomes dominant. Frictional heating (v _b · τ _b) provides a positive feedback mechanism to further increase basal temperatures and, hence, basal motion. Each peak in basal temperature coincides with a corresponding peak in the contribution of basal motion, the latter reaching ∼95%, the same fraction as observed during the 1982–83 surge of Variegated Glacier (Reference KambKamb and others, 1985). Basal motion facilitates fast flow, and results in rapid thinning and surface flattening (not shown). The resulting reduction in basal shear stress limits internal heat sources and lowers the contribution of basal motion to the overall ice flow. It should be noted that the selected point has a bedrock elevation of 240 m a.s.l., i.e. the sliding law parameters at this location remain constant during all experiments. The differences between the experiments are thus of a purely geometric nature and only depend on the efficiency of the ice flux further down-stream, where the ice is grounded below sea level and sliding parameters subjected to changes. Enhanced basal motion of the marine grounded ice also entails geometric changes in the interior ice cap, grounded above sea level, and is conducive to the generation of surge-type behaviour. The simulated surge-type behaviour therefore appears to be ultimately controlled by the geometric evolution of the glacier, as mentioned by Reference ClarkeClarke (1976) and Reference RaymondRaymond (1987).

Fig. 8. (a) Ice thickness evolution, (b) basal temperature relative to pmp, and (c) relative contribution of basal motion to the overall ice flow, i.e. the surface motion, at one specific location on Etonbreen. Time-series data are given for the period 25–30 ka for all five experiments using γ = 0.2°C.

5.2. Numerical robustness

The mechanism initiating oscillatory behaviour of ice sheets evokes some controversy among theoretical glaciologists. Reference Bueler, Brown and LingleBueler and others (2007) and Reference Bueler and BrownBueler and Brown (2009) pointed out a numerical error arising from the abrupt transition from no-slip conditions at the ice base to full basal sliding that is employed by models such as those applied to EISMINT II (European Ice-Sheet Modeling Initiative II) experiment H (Reference PaynePayne and others, 2000). The numerical artefact arises from a singularity in the horizontal velocity field that leads to infinite vertical velocities and hence, unrealistic strain heating. The effect is smeared out in the case of a large horizontal grid spacing of the order of 10 km, typical for large-scale ice-sheet models, but becomes significant at finer grids (Reference Bueler and BrownBueler and Brown, 2009). Submelt sliding removes the singularity in the velocity field. Reference Greve, Takahama and CalovGreve and others (2006) and Reference CalovCalov and others (2010) demonstrated that large-scale oscillations can also occur for a smooth transition of the warm- and cold-based sliding regimes and also for a higher-order model that incorporates longitudinal stresses through a combination of SIA and SSA. The implementation of submelt sliding is, in addition, motivated by physical reasoning. Sliding is likely to occur on a subgrid scale, before the entire area represented by a single gridpoint has reached the pmp.

To test the robustness of the model experiments, we ran experiment 5 using γ = 0.05°C at different spatial (1, 2 and 4 km) and temporal (1, 0.1, 0.025 and 0.00625 years) resolutions. This model run was expected to be the most sensitive to numerical inaccuracies as it represents maximum sliding enhancement and the most abrupt transition from no-slip conditions to full basal sliding. With a time-step of 1 year, only the 4 km resolution run performed stably. However, the generated ice-cap size is ∼10% smaller than for the higher temporal resolutions for which the results converge. The runs at 2 and 1 km spatial resolution converge at a time-step of 0.025 years, although the result at 0.1 year differs only by ∼1%. We conclude that a time-step of 0.025 years is sufficiently small for all three spatial resolutions. Using this, the total volume ranges from 2471 km³ at 1 km resolution to 2490 km³ at 2 km resolution (<0.8%) and the total area from 8257 to 8437 km² (<2.2%). The resulting mean ice thickness varies from 293 m at 1 km resolution to 300 m at 4 km resolution (<2.2%). Slightly larger differences (12.5%) exist for the area at the pmp, ranging from 1257 km² at 1 km up to 1424 km² at 4 km resolution, probably associated with a larger mean ice thickness in the case of the 4 km resolution. The oscillatory behaviour is generally unaffected by the spatial resolution. Fluctuations in volume range from 1.2% (4 km) to 1.7% (1 km), for areal extent from 1.9% (4 km) to 2.2% (1 km), for mean ice thickness from 1.9% (2 km) to 2.2% (1 km) and from 16.0% (1 km) to 19.7% (4 km) for the basal area at the pmp. The runs at 1 km spatial resolution give the largest temporal variation in most cases, but there is no unambiguous trend of the scalar variable values with resolution, and the differences do not indicate numerical instabilities. The slight differences in the mean and standard deviations of the presented variables can be explained by the fact that individual flow units are represented differently at the various spatial resolutions: some small and narrow flow units may be present at high and absent at low spatial resolution.