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 (Reference Hagen, Eiken, Kohler and MelvoldHagen and others, 2005; Reference Yde and PaascheYde 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−1 are typical for glaciers during the surge phase. Such high flow velocities are achieved through basal motion (Reference ClarkeClarke, 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−1, 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 (Reference ClarkeClarke, 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) (Reference Meier and PostMeier and Post, 1969; Reference Dolgoushin and OsipovaDolgoushin 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 (Reference Clarke, Collins and ThompsonClarke 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 (Reference Meier and PostMeier 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 (Reference KambKamb 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. Reference KambKamb (1987) and Reference RaymondRaymond (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. Reference Clarke, Collins and ThompsonClarke 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 (Reference ShreveShreve, 1984). Therefore, the thermal evolution of the bed plays a decisive role (Reference ClarkeClarke, 1976).
The occurrence of surge-type glaciers worldwide is spatially limited to several clusters, suggesting that certain regional conditions favour surge behaviour (Reference Meier and PostMeier and Post, 1969; Reference Dolgoushin and OsipovaDolgoushin 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 (Reference Hamilton and DowdeswellHamilton and Dowdeswell, 1996; Reference Jiskoot, Murray and BoyleJiskoot and others, 2000). The sparse observations on glacier surges that are available conform with these theoretical results. Reference MurrayMurray 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. Reference Glasser and HambreyGlasser and Hambrey, 2001).
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 (Reference MacAyealMacAyeal, 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 (Reference MacAyealMacAyeal, 1993; Reference PaynePayne, 1995; Reference Hindmarsh and Le MeurHindmarsh and Le Meur, 2001; Reference Calov, Ganopolski, Petoukhov, Claussen and GreveCalov and others, 2002, Reference Calov2010; Reference Greve, Takahama and CalovGreve and others, 2006). Reference CalovCalov 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 (Reference CalovCalov and others, 2010), the smoothing of the transition over a broad subfreezing temperature range may prohibit oscillatory behaviour (Reference Greve, Takahama and CalovGreve and others, 2006). The amplitude of the oscillations increased with increased enhancement of basal sliding (Reference CalovCalov 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 (Reference CalovCalov 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 km2, 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 (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).
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−1 and marginal thinning of 1–3 m a−1 (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−1), interrupted by distinct fast-flow units having velocities in excess of 100 m a−1 (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 km3 a−1, 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−2 (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.
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 (Reference GreveGreve, 1997). The model takes advantage of the simplification of the field equations known as the shallow-ice approximation (Reference HutterHutter, 1983; Reference MorlandMorland, 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 (Reference RitzRitz, 1987). Isostatic adjustment of the lithosphere due to changes in ice load is computed by the elastic-lithosphere/relaxing-asthenosphere approach (see Reference Greve and BlatterGreve 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 Reference GreveGreve (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.
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’:
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−1. 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,
(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,
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.
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 ρ −1 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:
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−4 m−1 a−1, 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−3 m−1 a−1 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−1 (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−1),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−2 (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.
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.
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.
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).
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−1), interrupted by several fast-flow units with moderate velocities of 50–150 m a−1 . 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.
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).
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−1 up to several km a−1), 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.
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 Reference Fowler, Murray and NgFowler 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. 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).
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 km3 at 1 km resolution to 2490 km3 at 2 km resolution (<0.8%) and the total area from 8257 to 8437 km2 (<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 km2 at 1 km up to 1424 km2 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 Reference Greve, Takahama and CalovGreve and others (2006) and Reference CalovCalov 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.
We thank H. Jiskoot and an anonymous reviewer for thorough and constructive comments that helped improve this paper. This study is a contribution to the International Polar Year project GLACIODYN funded by the Norwegian Research Council (grant 176076/S30). The final stage was supported by funding to the ice2sea project from the European Union 7th Framework Programme, grant No. 226375, ice2sea contribution No. 024. T.D. was supported by a short-term predoctoral fellowship through the Japan Society for the Promotion of Science (research stay at the Institute of Low Temperature Science, Sapporo) and an Arctic field grant through the Svalbard Science Forum (fieldwork). T.D. thanks E. Vasilenko for guidance through relevant Russian publications, and O. Salvigsen for providing insight into the geology of Nordaustlandet.