Permanent fast flow versus cyclic surge behaviour: numerical simulations of the Austfonna ice cap, Svalbard

Abstract

A large part of the ice flux within ice caps occurs through spatially limited fast-flowing units. Some of them permanently maintain fast flow, whereas others operate in an oscillatory mode, characterized by short-lived active phases followed by long quiescent phases. This surge-type behaviour results from intrinsic rather than external factors, thus complicating estimates of glacier response to climate change. Here we present numerical model results from Austfonna, an ice cap on Svalbard that comprises several surge-type basins. Previous studies have suggested a thermally controlled soft-bed surge mechanism for Svalbard. We systematically change the parameters that govern the nature of basal motion and thereby control the transition between permanent and oscillatory fast flow. Surge-type behaviour is realized by a relatively abrupt onset of basal sliding when basal temperatures approach the pressure-melting point and enhanced sliding of marine grounded ice. Irrespective of the dynamic regime, the absence of considerable volumes of temperate ice, both in the observed and simulated ice cap, indicates that fast flow is accomplished by basal motion over a temperate bed. Given an idealized present-day climate, the equilibrium ice-cap size varies significantly, depending on the chosen parameters.

1. Introduction

The role of glaciers as indicators of climate change and their major contribution to sea-level rise is widely acknowledged. However, extraction of climate signals from glaciers is not straightforward, because the history, current state and future evolution of glaciers result from the interplay of external factors (e.g. changes in surface air temperature or precipitation) and also intrinsic glacier dynamics (Hagen and others, 2005; Yde and Paasche, 2010). Surge-type glaciers illustrate this affinity in a drastic way. Periodically, they experience significant changes in glacier dynamics and geometry that are largely decoupled from climate variability, challenging both observational glaciologists and modellers to deliver estimates about their response to climate change. This is especially true for predictions on decadal timescales, in great demand by the public in general and decision makers in particular.

Surface velocities of the order of several hundred to 1000 m a⁻¹ are typical for glaciers during the surge phase. Such high flow velocities are achieved through basal motion (Clarke, 1987) and require a temperate glacier bed, i.e. basal temperatures at the pressure-melting point (pmp). For a frozen bed, ice deformation is the exclusive mode of glacier motion, typically yielding moderate surface velocities of the order of 1–10 m a⁻¹, depending on shear stress and ice temperature. If conditions for basal motion are maintained, the dynamics of a glacier may be characterized by permanent fast flow, such as observed for ice streams, otherwise a temporal pattern of alternating fast and slow glacier flow may evolve. These different dynamic regimes have an impact on the glacier’s net mass balance. Surge-type glaciers are characterized by an oscillatory mode of equilibrium and do not maintain a steady mass flux that equals the theoretical balance flux to maintain a steady-state surface profile (Clarke, 1987). Instead, their dynamics are characterized by a long quiescent phase of inefficient ice flow, undershooting the balance flux, and a short-lived surge phase of super-efficient ice flow, greatly overshooting the balance flux. During the quiescent phase, the glacier can be divided into an active thickening zone, the ‘reservoir zone’, and an almost stagnant depleted zone, the ‘receiving zone’, separated by the dynamic balance line (DBL) (Meier and Post, 1969; Dolgoushin and Osipova, 1975). The reservoir and receiving zones do not usually coincide with the glacier’s accumulation and ablation zones, which are separated by the equilibrium-line altitude (ELA); instead the DBL marks a boundary zone where glacier outflow is restricted (Clarke and others, 1984). The distribution of mass from the reservoir zone into a receiving zone during the surge may be accompanied by a significant advance of the terminus (Meier and Post, 1969).

Observations on Variegated Glacier, Alaska, suggest that 95% of the motion during the surge in 1982–83 was due to basal sliding and only ∼5% due to internal deformation (Kamb and others, 1985). Basal processes are therefore considered fundamental in initiating and maintaining glacier surges. Factors determining the absolute contribution of basal motion to the overall ice flow include the thermal regime at the glacier base, basal shear stress, basal water pressure and the lithology of the underlying bedrock, as well as the presence of deformable sediments. Kamb (1987) and Raymond (1987) discussed implications of different basal hydraulic-drainage systems and pointed out that a highly pressurized system of linked cavities may facilitate fast flow and prevail during a surge, while a switch to a hydrologically efficient channel system reduces the basal water pressure and may cause surge termination. Clarke and others (1984) suggested the presence of highly deformable, water-saturated sediments as an alternative explanation of highly enhanced basal motion. Basal motion of polythermal ice bodies is spatially and temporally restricted to basal areas at the pmp. The cold basal areas are frozen to the ground and basal motion is negligible on timescales smaller than 100 ka (Shreve, 1984). Therefore, the thermal evolution of the bed plays a decisive role (Clarke, 1976).

The occurrence of surge-type glaciers worldwide is spatially limited to several clusters, suggesting that certain regional conditions favour surge behaviour (Meier and Post, 1969; Dolgoushin and Osipova, 1975). One such cluster is found on Svalbard, a highly glacierized archipelago in the Eurasian Arctic (Fig. 1). Statistical model results indicate a thermally controlled soft-bed surge mechanism with largest surge probability for long polythermal glaciers resting on unlithified or easily erodible beds (Hamilton and Dowdeswell, 1996; Jiskoot and others, 2000). The sparse observations on glacier surges that are available conform with these theoretical results. Murray and others (2000) found that the advance of the surge front of Bakaninbreen, a 17 km long surge-type glacier in southern Svalbard, has been associated with the down-glacier expansion of a temperate ice base. Furthermore, they found evidence for subglacial sediments within the area affected by the surge. Observations of soft beds underlying polythermal surge-type glaciers are also available from other areas of Svalbard (e.g. Glasser and Hambrey, 2001).

Fig. 1. Map of the Austfonna ice cap and its location within the Svalbard archipelago (inset). Contours show surface elevations (m a.s.l.), blue solid lines indicate drainage basins (Hagen and others, 1993) and the red solid line the main ice divide. Basins of interest are marked with bold letters: the three known surge-type basins Etonbreen (ET), terminating into Wahlenbergfjorden (WBF), and Bråsvellbreen (BR) and Basin-3 (B3), terminating into the Barents Sea. Other basins of interest are Duvebreen (DU), Leighbreen (LE) and Basin-5 (B5).

Large-scale surge behaviour bears the potential for abrupt ice-sheet disintegration. Ice-rafted debris found in marine sediment cores from the North Atlantic provides a record of exceptional iceberg-calving events. These events, known as ‘Heinrich events’, occurred at time intervals of ∼7 ka. It has been proposed they are associated with periodic surges (‘binge–purge oscillations’) of the prehistoric Laurentide ice sheet (MacAyeal, 1993) that covered large parts of North America between ∼95 and 20 ka BP. Oscillatory ice-sheet behaviour similar to Heinrich events has been generated by a number of models for simplified set-ups (MacAyeal, 1993; Payne, 1995; Hindmarsh and Le Meur, 2001; Calov and others, 2002, 2010; Greve and others, 2006). Calov and others (2010) compare nine thermomechanical ice-sheet models, one model utilizing a combination of the shallow-shelf (SSA) and shallow-ice approximation (SIA), the others being purely SIA type, and showed that seven models, including the SSA/SIA model, are capable of producing oscillatory behaviour. The oscillations were characterized either by single- or multi-peaked spectra; however, differences in periodicity are large. Although a discrete switch from no-slip conditions to full basal sliding is not required for the occurrence of oscillations (Calov and others, 2010), the smoothing of the transition over a broad subfreezing temperature range may prohibit oscillatory behaviour (Greve and others, 2006). The amplitude of the oscillations increased with increased enhancement of basal sliding (Calov and others, 2010). Furthermore, the influence of changed climate conditions was tested. No clear relationship between oscillatory behaviour and air temperature was found, but an increase in accumulation, leading to an increased ice thickness, forced some of the models from an oscillatory mode to a mode of steady fast flow. For the remaining models that exhibited oscillatory behaviour, the period of the oscillations decreased, associated with a decrease in time required for the build-up of a critical ice thickness in the reservoir zone at which the surge is released (Calov and others, 2010).

Here we present numerical simulations of Austfonna, a large polythermal ice cap in Svalbard (Fig. 1). Several of the drainage basins of this ice cap are known to have surged in the past. The aim is to generate an ice cap that is similar in size and dynamic behaviour to present-day Austfonna and that provides suitable initial conditions for prognostic model runs in future studies. We investigate a number of model parameters that exert controls on the nature of the dynamics of the simulated ice cap. In light of the suggested surge mechanism for Svalbard, we focus, in particular, on the description of basal sliding and the role of the basal thermal regime.

2. The Austfonna Ice Cap

Austfonna (the Norwegian name means ‘eastern ice cap’), is located on the island Nordaustlandet (80° N) in the northeast of Svalbard (Fig. 1). Covering an area of >8000 km², it is by far the largest individual ice body on the highly glacierized archipelago.

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 (Dowdeswell 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 (Schytt, 1969; Solheim, 1991), considerably longer than for smaller glaciers in western Svalbard or other areas in general (Dowdeswell and others, 1991). Over 28% of the entire area covered by Austfonna lies below sea level (Dowdeswell, 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 (Dowdeswell 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⁻¹ (Moholdt 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 (Dowdeswell, 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 (Dowdeswell, 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⁻¹ (Dowdeswell and others, 1999). Bevan 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 (Bamber 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 (Pinglot and others, 2001), although large interannual variations in the total amount of snow accumulation have been observed in recent years (Taurisano and others, 2007; Dunse and others, 2009). These observations also indicate a clear asymmetry in snow accumulation that is reflected in the net SMB (Schuler and others, 2007) and, hence, the ELA: ∼450 m a.s.l. in the northwest and 300 m a.s.l. in the southeast (Moholdt and others, 2010). Geodetically derived mass balance over the period 2002–08 suggests a net SMB close to zero (Moholdt 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 (Dowdeswell and others, 2008; Moholdt and others, 2010).

2.2. Thermal regime

Austfonna’s thermal structure is described as polythermal (Dowdeswell 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. Zagorodnov 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 (Zagorodnov 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 (Zagorodnov and others, 1989b). The temperature gradient in the lower 400 m suggests a moderate geothermal heat flux of ∼40 mW m⁻² (Ignatieva and Macheret, 1991). Information on the thermal structure of the ice cap is also contained in the airborne RES data (Macheret 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 (Björnsson and others, 1996). The Institute of Geography, USSR Academy of Sciences, collected such data in 1984 (Macheret and Vasilenko, 1988). Similar to the RES data of Dowdeswell and others (1986), the absence of a continuous IRH indicates that significant volumes of temperate ice do not exist. Kotlyakov and Macheret (1987) and Macheret 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 (Macheret and Vasilenko, 1988; Dowdeswell and others, 1999). Furthermore, Solheim and Pfirman (1985) and Solheim (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 (Hjort and others, 1995; Blake, 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 (Lauritzen 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. Model

We employ the ice-sheet model SICOPOLIS (SImulation COde for POLythermal Ice Sheets; http://sicopolis.greveweb.net/), which solves the thermomechanically coupled field equations: the equations of mass, momentum and energy, along with Glen’s flow law as constitutive equation and boundary conditions for the free surface and the bed. In addition, a jump condition is solved for the cold/temperate transition surface that separates cold and temperate ice regions (Greve, 1997). The model takes advantage of the simplification of the field equations known as the shallow-ice approximation (Hutter, 1983; Morland, 1984). Under the assumption that typical length scales exceed vertical scales, the stress in the ice reduces to the shear stress, τ, in the horizontal plane, and the hydrostatic pressure, P = ρg(h − z), where ρ is the ice density, g the gravitational acceleration, h the surface elevation and z the height within the ice column counted positive upward. The resulting effective stress, σ = P|∇h|, reveals that ice flow only depends on local ice thickness and surface slope and is aligned with the surface gradient. The total ice flow is the sum of the ice flow by internal deformation and basal sliding; the latter is discussed in more detail in section 3.1.

The thermomechanical coupling of the model equations stems from the dependence of the rate factor on ice temperature and water content. The temperature evolution in the cold ice is balanced by horizontal and vertical advection, vertical heat conduction and dissipative strain heating. In the temperate ice the temperature follows directly from the pressure field and the energy balance is provided by an evolution equation for the water content. Basal melting is computed by balancing the heat fluxes from the ice and the lithosphere and the heat generated by basal sliding. The geothermal heat flux is applied as a boundary condition 5 km below the ice base to account for the thermal inertia of the lithosphere (Ritz, 1987). Isostatic adjustment of the lithosphere due to changes in ice load is computed by the elastic-lithosphere/relaxing-asthenosphere approach (see Greve and Blatter, 2009, and references therein).

Simulations are run on a fixed rectangular grid at 2 km horizontal resolution encompassing 69 × 67 gridpoints in the horizontal plane. The vertical resolution comprises 81 gridpoints in the cold ice column, with increased spatial resolution towards the base, and 11 equally distributed gridpoints each in the temperate part of the ice column, if applicable, and the lithosphere. A fixed time-step of 0.025 years is used for all simulations in this study, except when a different value is stated. The values of relevant physical parameters are listed in Table 1. The model is forced by monthly mean air temperature, monthly mean precipitation, sea-level stand and geothermal heat flux. See Greve (1997) for more details on how SICOPOLIS computes the surface mass balance from the given input fields and accounts for firn warming and superimposed ice formation.

Table 1. Physical parameters used in all model experiments

3.1. Basal sliding

The dynamic boundary condition at the glacier bed is given by a Weertman-type sliding law (Greve 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 Hindmarsh and Le Meur (2001) and Greve (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 (Pfeffer, 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 Marshall 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 (Tarasov and Peltier, 1997) or a function of sea level (Zweck 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 (Albrecht and others, 2010). This approach provides a promising basis for implementation of physically based calving laws, such as the waterline crevasse-depth model of Benn 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 (Jakobsson 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 Dowdeswell 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 (Vasilenko and others, 2010). The profiles follow the routes investigated by the University of Oslo and NPI for mass-balance purposes since 1998 (Schuler and others, 2007; Taurisano and others, 2007; Dunse 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⁻¹ (Schuler 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 (Smith and Barstad, 2004; Schuler 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 (Moholdt and others, 2010).

The vertical temperature profile through the ice at the summit presented by Zagorodnov and others (1989b) indicates a geothermal heat flux of ∼40 mW m⁻² (Ignatieva 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; Tulaczyk 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. Results

In section 4.1 we provide an overview of all experiments (experiments 1–5) over the full range of γ. Based on the temporal variability of ice-cap variables, in the following sections we then present particular runs that generate an ice cap characterized by permanent fast flow (section 4.2) and oscillatory fast flow, including cyclic surge-type behaviour (section 4.3), and show contour maps of selected field variables as well as the temporal evolution of the ice thickness at specific locations to provide indicators of individual basin activity.

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

We have shown that by choosing distinct combinations of the parameters in the basal-sliding law, representing a relatively abrupt onset of basal motion and enhanced sliding of marine grounded ice, SICOPOLIS generates an ice cap with similar dynamics to Austfonna. With regard to the observational evidence, the experiment describing soft sediments beneath marine grounded ice is considered the most realistic scenario. However, a perfect match of the simulated and modern observed ice cap in terms of glacier outline, geometry and surface velocities is not achieved. Each experiment represents homogeneous idealizations of the basal conditions, whereas actual basal conditions, such as bedrock lithology, thickness and water saturation of subglacial sediments, vary both spatially and temporally. Therefore, sliding parameters should ideally be functions of space and time, rather than constants. Furthermore, variations in basal water pressure associated with changes in the basal hydraulic system will occur beneath the real ice cap, but are not taken into account in the model. Nevertheless, the evolution of the subglacial hydrology is expected to go hand in hand with the basal thermal regime: once the basal temperature reaches the pmp, basal melt can occur. If the hydraulic conductivity of the bedrock and an emerging basal hydraulic drainage system are low, local water storage increases the basal water pressure. This weakens subglacial sediments, if present, facilitating enhanced basal motion as represented by experiment 5 and providing a potential surge mechanism as suggested by Fowler and others (2001). In the following, we describe the evolution of a surge cycle and the underlying mechanisms.

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. Bjö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 Payne (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 Calov and others (2002) and agree with observations of polythermal surge-type glaciers in other parts of Svalbard (Sund 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 (Schytt, 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 (Kamb 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 Clarke (1976) and Raymond (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. Bueler and others (2007) and Bueler 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 (Payne 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 (Bueler and Brown, 2009). Submelt sliding removes the singularity in the velocity field. Greve and others (2006) and Calov 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.

6. Concluding Remarks

The simulated dynamic behaviour of Austfonna is sensitive to the parameter combination of the applied sliding law. Submelt sliding in a defined temperature range and sliding enhancement of marine grounded ice determine whether spatially limited flow units operate in a mode of permanent or oscillatory fast flow. These results are in line with previous model experiments for the idealized Heinrich Event INtercOmparison (HEINO) set-up by Greve and others (2006) and Calov and others (2010). In the case of cyclic surge behaviour, enhanced sliding during the active phase is required to draw down the ice surface sufficiently and to reestablish a cold base, thereby initiating the quiescent phase. Considerable volumes of temperate ice do not occur in any of the model runs, indicating that fast flow is accomplished exclusively by basal motion over a temperate bed. The dynamic regime is shown to have a significant impact on the equilibrium size of Austfonna. Ice volume and areal extent decrease with both enhanced sliding and an increased temperature range at which submelt sliding is allowed to occur. The oscillatory behaviour generated by the model conforms with observational evidence of previous surge activity. However, none of the experiments result in an exact match of the simulated and modern observed ice cap, since the latter follows from an external-climatic and intrinsic-dynamic history with basal conditions being functions of space and time. Nevertheless, we have achieved a set of idealized initial conditions for prognostic model runs that allow us to assess the uncertainties in the dynamic response of the ice cap to climate change.