1. Introduction
Climatic variations are known to exert an external control on the dynamic evolution of glaciers and ice sheets. However, not all observed changes in ice movement and geometry evolution can be ascribed to changes in the climate system. Ice sheets evolve as a product of the interplay of geometry, thermodynamics and hydrology. Hence, accurate simulation of ice evolution requires proper treatment of both internal and external conditions. Irrespective of the climate forcing, periodic changes in ice geometry and flow velocities may occur.
Numerical modelling of periodic behaviour in ice masses has helped to provide insight into how internal instabilities may trigger oscillatory flow. Reference BuddBudd (1975) was the first to link basal sliding velocities to the thermal regime in a numerical model and simulated surging behaviour. Later work by Reference OerlemansOerlemans (1983) illustrated the relevance of the explicit treatment of basal water and including temperature advection in the heat equation. A ‘binge/purge’ mechanism, proposed by Reference MacAyealMacAyeal (1993), relates periods of excessive ice advance to changes in the thermal regime. Reference PaynePayne (1995) demonstrated that the interaction of basal dynamics and the basal heat budget may lead to limit cycles in the basal thermal regime. Similar multi-millennial oscillations were found in a three-dimensional (3-D) model, developed by Reference Marshall and ClarkeMarshall and Clarke (1997), in studies by Reference Calov, Ganopolski, Petoukhov, Claussen and GreveCalov and others (2002) and Reference Papa, Mysak and WangPapa and others (2006). Results of the Heinrich Event INtercOmparison (HEINO) as part of the Ice- Sheet Model Intercomparison Project (ISMIP), summarized by Reference CalovCalov and others (2010), demonstrate the ability of ice sheet models to simulate Heinrich-type oscillations. Heinrich events are periods of excessive iceberg calving from the Laurentide ice sheet during the last ice age, for which observational evidence originates from layers of ice-rafted debris in deep-sea sediment cores (Reference HeinrichHeinrich, 1988; Reference BondBond and others, 1992). They are thought to be related to periodic ice- stream surging in Hudson Strait (Reference Clarke, Marshall, Hillaire-Marcel, Bilodeau, Veiga- Pires, Clark, Webb and KeigwinClarke and others, 1999).
Surging of valley glaciers is an example of periodic ice movement, in which one cycle is characterized by a long period of quiescent flow, followed by a shorter period of high ice velocities. Theoretical work has emphasized the crucial role of basal conditions, facilitating sliding of the ice, in the occurrence of surging (Reference Clarke, Collins and ThompsonClarke and others, 1984; Reference ClarkeClarke, 1987). Several mechanisms have been proposed to explain surging, of which changes in the efficiency of the drainage system (Reference KambKamb and others, 1985) and a thermally controlled surging mechanism (Reference MurrayMurray and others, 2000; Reference Fowler, Straughan, Greve, Ehrentraut and WangFowler and others, 2001) are most widely accepted. Reference Murray, Strozzi, Luckman, Jiskoot and ChristakosMurray and others (2003) concluded, based on observed regional differences between surge dynamics, that at least two distinct surge mechanisms exist. Understanding the interaction between basal resistance to ice flow, the basal hydrological system and melting conditions near the base is required to address mechanisms causing surge-like behaviour.
Runaway feedback mechanisms arising from the interaction of basal thermodynamics, water production and basal shearing have been linked to flow variations of ice streams along the Siple Coast, West Antarctica (Reference Payne and DongelmansPayne and Dongel- mans, 1997; Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000b; Reference Anandakrishnan, Alley, Jacobel, Conway, Alley and BindschadlerAnandakrishnan and others, 2001), observed in satellite imagery of the Ross Ice Shelf and its tributaries (Reference Bindschadler and VornbergerBindschadler and Vornberger, 1998; Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others, 2000; Reference Jacobel, Scambos, Nereson and RaymondJacobel and others, 2000). Reference Payne and DongelmansPayne and Dongelmans (1997) and Reference PaynePayne (1999) address the role of feedback between internal shearing and ice viscosity (Reference Clarke, Nitsan and PatersonClarke and others, 1977) in ice-stream flow. Flow variability of connected ice streams is linked to changes in the drainage basin extent of the individual ice streams. Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others (2000b) use an undrained plastic-bed model (Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000a) to show that an instability arisingfrom the interaction of frictional heating, water storage and sliding velocities may lead to transient ice- stream behaviour. This instability has been proposed as a possible mechanism for the shutdown of Kamb Ice Stream (formerly Ice Stream C), which occurred ~1 60 years ago (Reference RoseRose, 1979; Reference Retzlaff and BentleyRetzlaff and Bentley, 1993; Reference Bougamont, Tulaczyk and JoughinBougamont and others, 2003). Reference Alley, Anandakrishnan, Bentley and LordAlley and others (1994) and Reference Anandakrishnan and AlleyAnandakrishnan and Alley (1997) argue that the stoppage of Kamb Ice Stream may have been caused by a topographic change leading to enhanced water flow to neighbouring Alley Ice Stream (formerly Ice Stream B2).
In this study, a 3-D thermomechanically coupled ice- sheet model, the Parallel Ice Sheet Model (PISM), is used to simulate ice flow on a synthetic valley-shaped basal topography. Within the model, combined stress balances, a detailed energy conservation scheme and an evolving basal hydrology interact within a single framework. Coupling of these model components in a 3-D framework allows for a detailed study of mutual interactions and leads to new insights into the processes involved in oscillations. This study aims to physically describe different types of oscillations that appear in the model output, in order to identify the role of internal destabilizing feedback mechanisms in maintaining oscillations. This involves a discussion of the interplay of stresses, heat fluxes and basal water. In simulations of Austfonna, Svalbard, using SICOPOLIS (Simulation COde for POLythermal Ice Sheets), the choice of specific sliding-law parameters was shown to facilitate characteristic dynamics of outlet glaciers, ranging from steady fast flow to cyclic surge behaviour (Reference Dunse, Greve, Schuler and HagenDunse and others, 2011). In line with Reference Dunse, Greve, Schuler and HagenDunse and others (2011), we study the dependence of the occurrence and characteristics of periodic behaviour on parameters controlling the basal resistance to sliding flow (Section 3.1). SICOPOLIS and PISM differ considerably with respect to the description of sliding. The use of a synthetic topography facilitates the analysis of oscillation features. The impact of modifying parameters defining the climate input on the cyclic behaviour is discussed in Section 3.2, and the role of the height dependence of the surface mass balance and temperature is studied separately (Section 3.3). Finally, the model robustness under grid refinement is verified using resolution experiments (Section 3.4).
2. Model and Set-Up
Numerical experiments are conducted with an open-source ice dynamical model, PISM. Within the model, ice flow is simulated as a product of mechanical, thermodynamical and hydrological processes. A more complete description of the model is given by Reference Bueler and BrownBueler and Brown (2009). Here, only model properties relevant to the simulated cyclic behaviour are discussed.
2.1. Model description
PISM combines two shallow (small depth-to-width ratio) stress balances, namely the shallow-ice approximation (SIA) and the shallow-shelf approximation (SSA), which are computationally efficient schemes to simulate ice flow by internal deformation and ice-stream flow, respectively. The widely used non-sliding SIA (Reference Morland and JohnsonMorland and Johnson, 1980; Reference HutterHutter, 1983) is a stress balance in which the driving stress is balanced solely by internal shearing. In cases where basal lubrication significantly adds to the surface velocity, stresses associated with basal resistance to the flow and longitudinal stretching can no longer be neglected. The SSA (Reference MacAyealMacAyeal, 1989; Reference Weis, Greve and HutterWeis and others, 1999) assumes that the driving stress is fully balanced by membrane stresses and friction at the ice/bedrock interface, and is often used to simulate ice flow in ice streams and ice shelves. In PISM, deformational velocities from the SIA and sliding velocities from the SSA are weighted and averaged to achieve a smooth transition from shearing flow to sliding flow. The hybrid implementation of stresses appears to be a powerful tool to deal with the fluxmatching problem at the boundary between grounded and sliding ice (Reference SchoofSchoof, 2006; Reference Bueler, Brown and LingleBueler and others, 2007; Reference Bueler and BrownBueler and Brown, 2009; Reference GoldbergGoldberg, 2011; Reference WinkelmannWinkelmann and others, 2011). Furthermore, in terms of computational costs, this shallow model is significantly more efficient than higher- order or full-Stokes alternatives.
An energy conservation scheme, as described by Reference Aschwanden and BlatterAsch-wanden and Blatter (2009), computes temperature and liquid water content fields using an ‘enthalpy-gradient method’ (Reference PhamPham, 1995). The scheme solves the enthalpy-balance equation simultaneously in temperate ice (ice at pressuremelting point) and cold ice (ice below pressure-melting point) and is therefore particularly suited to model the thermodynamic state of polythermal ice masses. Within the ice, the energy scheme accounts for enthalpy advection, diffusion and production of enthalpy by strain heating. Within the bedrock layer, enthalpy evolves diffusively, with a prescribed lower boundary flux equal to the geothermal heat flux. At the ice/bed interface, frictional heating occurs as ice slides over the bedrock. When the ice at the base is at pressure-melting point, excess energy is used for melting the ice. A simple diffusion relation describes the evolution of the effective thickness of the water layer, W, at the base (Reference Bueler and BrownBueler and Brown, 2009):
where S denotes the local water production by basal melting and K is a diffusion coefficient (200 m2 a-2). S may become negative, in which case refreezing of basal water occurs as the basal temperature is below pressure-melting point. The basal water diffusion relation provides a simplified description of water transport, in which water diffuses at a prescribed rate (controlled by K) from high to low effective thickness of the water layer. Clearly, the model does not include changes in the hydrological drainage efficiency, which are often related to Alaskan-type surging behaviour (Reference KambKamb and others, 1985). The basal water pressure, p w, is given by
With ρ the density of ice (910 kg m-3), g gravitational acceleration (9.81 ms-2) and H ice thickness. W 0 denotes the preset maximum value for W (of 2 m) at which till saturation is achieved. α is a factor defining the maximum ratio of pore-water pressure (p w) to overburden pressure (p l = ρgH), which is achieved in the case of till saturation. In this study, we used α = 0.97, which is in accordance with observations that α ~ 1 (Reference Lüthi, Funk, Iken, Gogineni and TrufferLüthi and others, 2002). The basal water effective thickness at the ice margin is set to zero.
Sliding velocities in PISM are derived by iteratively solving the SSA, in which the basal shear stress, τb , has a pseudoplastic form:
where τ c denotes the yield stress, u b is the basal velocity vector, u 0 is a constant threshold velocity (100 m a-1) and q is the pseudoplasticity exponent, which can take values ranging from 0 (perfectly plastic till) to 1 (linearly viscous till) (Reference ClarkeClarke, 2005). In this study, the till is assumed to behave plastically or nearly plastically, with values of q ranging from 0 to 0.3. Nearly plastic deformation implies that some, but limited, till deformation occurs when the applied stress is slightly smaller than the yield stress. The yield stress formulation is given by (Reference Bueler and BrownBueler and Brown, 2009)
Here c 0 is the till cohesion (set to zero) and ϕ is a material parameter describing the material strength of the till, which is assumed to be spatially invariant, in contrast to Reference Bueler and BrownBueler and Brown (2009). Values of ϕ range from 7.5° to 30° for the different experiments. Note that α = 1 would imply that the yield stress goes to zero in the case of till saturation.
In PISM, ice dynamics, thermodynamics and hydrology are coupled within a single framework. Sliding velocities are computed from the SSA stress equation, in which the driving stress is balanced by membrane stresses and basal friction. The basal shear stress is a function of the basal water pressure (Eqns (3) and (4)), while the basal water pressure itself is affected by basal melting (Eqns (1) and (2)). Basal melting follows from the enthalpy equation, in which frictional heat production depends on the magnitude of sliding velocities. This is an example showing how the different regimes of the model are connected. The interaction between dynamics, thermodynamics and hydrology is crucial in the description of cyclic flow.
2.2. Numerical set-up
Model experiments were conducted on a synthetic bottom topography (Fig. 1), using a 88 × 56 horizontal grid (grid spacing 1.25 km). The basal profile is inclined such that the central flowline is at a 1.0° inclination, thereby forcing the ice to flow to the left in Figure 1. An increasingly sloping bed topography towards the edges in the y-direction and upstream in the x-direction prevents ice flow out of the model domain. Consequently, all ice volume variations are due to changes in the mass-balance regime, related to fluctuations in the extent of the ice cap.
Contour maps of (a) the synthetic bedrock topography, (b) the prescribed surface temperature and (c) the prescribed surface mass balance (mw.e. a −1).
The vertical grid contains 150 unevenly distributed layers (more concentrated near the base). In the bedrock thermal model, a vertical grid containing ten evenly distributed layers extending 50 m into the bed is applied. An adaptive time-stepping scheme selects the largest time-step for which stability of the mass-continuity and energy-conservation scheme is assured, thereby significantly improving the computational efficiency (Reference Bueler and BrownBueler and Brown, 2009).
2.3. Experiments
To a large extent, parameters ϕ and q control the area and magnitude of sliding flow, as they directly affect the basal resistance experienced by the moving ice. In a first set of experiments, referred to as ‘basal parameter experiments’, a systematic exploration of the ϕ -q parameter space is carried out to study the relation between these parameters and the presence and characteristics of cyclic flow (Section 3.1).
A second set of model runs, referred to as ‘climate parameter experiments’, is performed to study the sensitivity of the oscillation characteristics to variations in parameters controlling the surface mass balance and temperature (Section 3.2).
Finally, in a third set of experiments, referred to as ‘climate feedback experiments’, the impact of including a surface-height-dependent mass balance and surface temperature formulation on the oscillation characteristics is studied (Section 3.3). (In the previous sets of experiments, a location-dependent surface mass balance and temperature are prescribed.) Feedback mechanisms associated with the inclusion of a surface-height dependence of the specific mass balance and ice surface temperature may interfere with otherfeedback mechanisms active during oscillations, which complicates the analysis of oscillatory flow. Hence, the time- dependent effect of surface-height variations on the surface mass balance and temperature is studied separately.
2.4. Climate set-up
The prescribed surface mass balance (Fig. 1) is assumed to increase linearly along the flowline with the height of the bed at a rate of 0.006 m w.e. a-1 m-1 and has an upper limit of 0.33 m w.e. a-1. A zero mass balance is prescribed at the domain boundaries to prevent ice formation by accumulation at the boundary of the model domain. Forced by the valley shape of the bed profile with steepest slopes near the domain boundaries, ice flow is always directed inward, i.e. away from the domain boundaries, and hence no ice will form at the domain edges. The imposed mass-balance distribution, together with the valley-shaped bedrock, forces the ice mass to form an ablating tongue downstream of a wide accumulation basin.
The prescribed surface temperature (Fig. 1) is set to increase linearly with the height of the bedrock at a rate of 0.007°Cm-1. In the basal parameter experiments (Section 3.1), the equilibrium-line altitude (ELA), E, and the 0°C surface-temperature height, z T0, correspond to a bedrock height of 1350m (along the centre flowline) and 400 m, respectively. In the climate parameter experiments (Section 3.2) values for E and z T0 range from 1250 to 1450 m and -200 to 1000 m, respectively.
In the climate feedback experiments (Section 3.3) the mass balance and surface temperature are a function of the time-dependent surface height instead of the fixed basal topography. Therefore, E and z T0 are corrected for the altitudinal difference, which is required to obtain an ice cap of an approximately similar size to that in the corresponding basal parameter experiments. More specifically, values of E and z T0 were corrected for the mean ice thickness derived from the corresponding run without the height feedback.
In all experiments, the chosen climate set-up results in a polythermal ice mass, with temperate ice at the base of the tongue surrounded by cold ice upstream. The experimental set-up is not designed to replicate one specific type of observed cyclic behaviour (like small-scale surging of valley glaciers or large-scale ice-stream flow behaviour) but rather serves the purpose of identifying the processes and interactions involved in multiple types of oscillatory flow. We purposely selected a climate set-up that produces an ice cap with a polythermal basal thermal structure, since previous work has indicated the dependence of a variety of flow instabilities on fluctuations in basal melt conditions. For example, changes in the basal thermal regime have been linked to unstable cessation of ice-stream flow (Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000b) and thermally controlled surging (Reference ClarkeClarke, 1976; Reference Fowler, Straughan, Greve, Ehrentraut and WangFowler and others, 2001) and have been proposed to explain Heinrich-type oscillations (Reference MacAyealMacAyeal, 1993). Additionally, the current set-up allows for varying sliding and climate parameters over a wide range of values, while the ice margins do not cross the domain boundaries.
3. Results
An underlying assumption in the derivation of the SIA and SSA and higher-order models, based on the Stokes model for a slow-flowing fluid (Reference FowlerFowler, 1997), is that inertia from flow accelerations can be neglected in the stress equations. This implies that the flow field is determined fully by the ice geometry (thickness and slope), boundary stresses and ice viscosity. Ice flow at a certain time is therefore not directly affected by the modelled ice velocities of the previous time- step. Consequently, oscillatory behaviour of masses based on inertia (e.g. a spring-mass oscillation) is intrinsically infeasible in a ‘slow’ flow model. In cases where the shape of an ice mass is forced out of its equilibrium state, ice velocities will change, such that the surface profile will relax back to its equilibrium shape in a nearly diffusive way following the mass-continuity equation (the SIA is fully diffusive, whereas longitudinal stresses in the SSA allow for ice flow in the across-slope direction). In the absence of feedback mechanisms enhancing flow or sliding variations, geometric perturbations will not initiate a sustained oscillation. Cyclic flow behaviour therefore shows up as internal or external feedback processes activated as the ice mass restores itself from an out-of-equilibrium situation. The aim of this study is to qualitatively describe possible feedback mechanisms in order to understand the oscillatory behaviour found in simulated ice velocities.
3.1. Basal parameter experiments
First, the impact of changes in the parameters controlling the basal shear stress on the presence and characteristics of oscillatory flow is discussed. Reference Jiskoot, Murray and BoyleJiskoot and others (2000) have illustrated that sediment properties exert a strong control on the occurrence of surging behaviour in Svalbard. In terms of computational costs, it is efficient to start the basal parameter experiments with a fully developed ice mass instead of prescribing ice-free conditions at t = 0. Therefore, a 5ka initialization run is performed with (ϕ, q) = (20°, 0.0), which serves as input for the experiments with varying values of ϕ and q. Starting from a fully developed ice mass, model runs with values of the till strength, ϕ, ranging from 7.5° to 22.5° and values of the sliding exponent, q, ranging from 0.0 to 0.3 are performed.
Figure 2 shows contour plots of temporal mean values of six variables, averaged over the run time, in this case for (ϕ, q) = (12.5°, 0.0). The maps illustrate the formation of an ice tongue with a sliding zone on a temperate base, partially saturated with water produced by basal melting. Mean basal ice velocities (up to 280 ma~1) dominate over internal deformation velocities in the sliding zone. In cold- based areas, sliding is absent and surface velocities are much lower. The water saturation level in Figure 2 is determined by a balance between generation of water by basal melting and diffusion of water from a high to low saturation level (Eqn (1)). Along the flowline, this leads to a saturated zone bounded downstream by an unsaturated zone as basal water goes to zero at the ice edge, and bounded upstream by an unsaturated zone under the influence of cold ice advection.
Contour plots of (a) mean surface height, h, (b) ice thickness, H, (c) surface velocity, vsurf, (d) basal velocity, vbase, (e) water saturation level, w, and (f) basal temperature, T base, in the basal parameter experiment with (ϕ, q) = (12.5°, 0.0).
Figure 3 gives an overview of modelled flow behaviour in the output of model runs with various combinations of ϕ and q. Three types of flow are distinguished: (1) stationary fast (sliding) flow, (2) high-frequency oscillatory fast flow (period 114-169 years) and (3) low-frequency oscillatory fast flow (period 1000+ years). The till strength is a measure of the material resistance of the bed to the ice flow near the base (Eqn (3)). Reducing the till strength directly enhances ice velocities in the sliding zone, as the gravitational driving stress experiences less resistance from basal friction. The sliding exponent, q, defines to what extent the basal shear stress is proportional to the ice velocity. Increasing the value of q enhances basal resistance and thus lowers the ice velocities in the sliding zone. Concerning the different types of flow, Figure 3 illustrates that increasing ice velocities (by lowering ϕ or q) at some point leads to a shift from stationary ice flow to high-frequency oscillating flow. Further enhancement of ice velocities at some point leads to a transition from high- to low-frequency cyclic flow. The different flow types, transitions and trends in the oscillation characteristics are analysed in the remainder of this section.
Overview of oscillation features in the basal parameter experiments. The coloured circles indicate the flow type. The value tags represent the mean maximum sliding velocity (ma−1) (black), the oscillation amplitude (m a−1) (brown) and the oscillation period in years of the maximum sliding velocity (green). Crosses indicate absence of oscillatory behaviour. The periodicity for the highfrequency oscillations ranges from 114 to 169 years. Low-frequency oscillations have a period of 1000+ years. The black circles mark the high- and low-frequency experiments used for further analysis.
3.1.1. Basal water feedback
As mentioned above, oscillations in the model can only arise when a feedback mechanism, causing flow enhancement, prevents relaxation to a steady-state geometry. To illustrate such a mechanism, we consider a perturbed initial surface slope in the sliding zone, which is somewhat steeper than the equilibrium slope. Larger driving stresses cause sliding velocities to increase instantaneously. Greater sliding velocities induce more heat production by frictional heating and consequently more basal meltwater production. A higher water saturation level reduces basal resistance to the flow and thus leads to further flow acceleration. This flowenhancing feedback mechanism, referred to as ‘basal water feedback’, was described by Reference RobinRobin (1955). Among others, Reference ClarkeClarke (1976) and Reference Fowler, Straughan, Greve, Ehrentraut and WangFowler and others (2001) also pointed out the importance of this feedback mechanism and related it to surging. The basal water feedback provides an additional ‘push’ to the ice flow when recovering from a perturbed geometry in regions with a nearly saturated base. In saturated areas, enhanced basal melting does not lead to a further reduction of the basal shear stress. Next, we discuss the two types of cyclic behaviour illustrated in Figure 3 in which the basal water feedback plays a decisive role.
3.1.2. High-frequency oscillations
The two types of oscillations found in the basal parameter experiments share the property that they exist by virtue of the basal water feedback. Nevertheless, some fundamental differences between the oscillations exist, which explain the observed discrepancies in oscillation characteristics (Fig. 3). First, we focus on one cycle of a high-frequency oscillation, shown for a run with (ϕ, q) = (12.5°, 0.0). Temporal mean values of the ice thickness, H, basal velocity, vbase, and water saturation, w, of this case are shown in Figure 2. Figure 4 shows deviations of the latter variables from the temporal mean atfour different times, t (years), covering one full cycle.
Contour maps of deviations from (a) the mean of the ice thickness, δH = H(t) − H mean, (b) the basal velocity, δvbase = vbase(t) − vbase,mean, and (c) the water saturation level, δw = w(t) − w mean, during one cycle of a high-frequency oscillation with (ϕ,q) = (12.5°, 0.0); t is given in years.
The four snapshots of the thickness deviation indicate oscillatory behaviour of the surface profile with a repetition period of 129 years. Steepening of the sliding zone leads to a maximum surface slope in the sliding zone at t = 35 years, whereas flattening of the surface leads to a minimum at t =100 years. At t = 0 and 70 years the surface slope in the sliding zone is close to the mean slope. Regarding the basal velocities, one may expect extrema to occur simultaneously with extrema in the surface slope as the driving stress is affected. However, Figure 4 shows that extrema in basal velocities coincide with surface slopes in the sliding zone close to the mean. While the surface is flattening in the sliding zone from t = 35 to 70 years, basal velocities continue to increase. This counter-intuitive behaviour is possible as additional basal water production by enhanced velocities leads to a strong enough reduction of the basal resistance near the front of the sliding zone, thereby causing basal velocities to increase and the sliding zone to expand. Conversely, surface steepening in the sliding zone after t = 100 years is accompanied by a reduction in ice velocities (by amplified reduction of basal water at the front of the sliding zone). Note that basal water feedback is only effectively changing basal resistance in regions with an unsaturated till, i.e. close to the boundaries of the sliding zone. The flow enhancement by basal water feedback results in a phase shift of about one-quarter of a period between extrema in the surface slope in the sliding zone and basal velocities and is a necessary ingredient for these oscillations to occur.
The feedback is further illustrated in Figure 5, showing deviations from the mean of the driving stress, τd, the basal shear stress, τb and the membrane stress, τb. The membrane stress comprises the combined effect of longitudinal stresses and transverse shear stresses. The driving stress is mainly controlled by slope variations, whereas the basal water saturation level controls, to a large extent, the basal shear stress. Between t = 35 and 70 years, τb is reduced by enhanced basal melting. Enhanced sliding velocities lead to flattening of the surface in the sliding area, causing an overall reduction of the driving stress. At the same time, at the up- and downstream boundary of the sliding zone, enhanced velocity gradients cause the surface to steepen locally and thus lead to a local increase of the driving stress. This local steepening in combination with reduced basal resistance results in enhanced longitudinal stresses at the down- and upstream boundaries of the sliding zone, which cause the sliding zone to expand, even when the driving stress is overall below its mean value. The expansion upstream of the sliding zone is small, since the transition to cold ice hinders growth of the area with a saturated till, due to refreezing of basal water.
Contour maps of deviations from (a) the mean of the driving stress, δτ d = τ d(t)−τd,mean (b) the basal shear stress, δτ b = τ b(t)−τb,mean, and (c) the membrane stress, δτ m = τ m(t)− τm,mean, during one cycle of a high-frequency oscillation with (φ, q) = (12.5°, 0.0); t is given in years. Values are only shown in regions with a water saturation level greater than 70%.
Figure 6 presents time series of the ice volume and maximum basal velocity for the basal parameter experiments with q = 0.1. Despite the impact of the high-frequency oscillations on the flow velocities, the ice volume is hardly affected. Since in these experiments the mass balance is only location-dependent (not a function of surface height), the total mass budget of the ice cap is only affected when the areal extent of the ice mass is altered. High-frequency oscillations tend to have a small impact on the extent of the ice cap and thus on the total mass budget. Termination of ice flow acceleration is primarily caused by reduction of the surface slope rather than advection of cold ice from higher areas. This is reflected in the absence of significant variations in the extent of the temperate base. It should be mentioned that Reference Dunse, Greve, Schuler and HagenDunse and others (2011) also found high-frequency oscillations in a model in which sliding is directly linked to basal temperatures. A more detailed comparison of high- frequency oscillations found by Reference Dunse, Greve, Schuler and HagenDunse and others (2011) and in this study is required to identify similarities and differences in the related cyclic mechanisms and to understand the exact importance of water treatment.
Time series of (a) ice volume and (b) maximum basal velocity for basal parameter experiments with q = 0.1 and φ ranging from 8.75° to 22.5°; LF, HF and S denote low-frequency, high-frequency and steady fast flow, respectively.
3.1.3. Low-frequency oscillations
Further reduction of the till strength, ϕ, in the basal parameter experiments at some point leads to a transition to a second type of oscillation, the low-frequency oscillation (Fig. 3). Figure 7 shows snapshots of the ice thickness, sliding velocities, water saturation level and basal temperatures at six different times, t (years), in a cycle for a run with (ϕ, q) = (7.5°, 0.3). At t = 0 years sliding is absent, while the surface recovers from a previous tongue advance. Excess energy from strain heating (by internal shearing) and geothermal heating in the temperate zone induces basal melting. Consequently, the water saturation level increases and basal sliding initiates as the yield stress weakens relative to the driving stress. As the tongue advances, the surface of the tongue flattens, thereby reducing the driving stress. Local steepening of the surface occurs down- and upstream of the sliding area. Enhanced driving stresses together with reduced basal resistance, amplified by the basal water feedback, cause flow enhancement and extension of the sliding zone in both the down- and upstream directions. Enhanced strain heating in steep regions promotes basal melting and thickening of the basal water layer. (Recall that water flows diffusively from high to low water thickness (Eqn (1)).) The saturated zone therefore expands in both the down- and upstream directions, despite refreezing at the transition from temperate to cold basal ice. Greatest sliding velocities are reached at t = 270 years. The tongue now extends far into the valley, causing the ice cap to have a negative mass budget. Further acceleration of the ice-stream flow is in part prevented by a limited inflow of ice coming from the accumulation zone and in part by the advection of cold ice into the sliding zone causing basal water to refreeze and the basal shear stress to increase. Consequently, the tongue of the ice cap will thin as surface melting is no longer compensated for by ice inflow. The reduced insulation and basal shear stress after thinning further reduce the available heat at the base and ultimately cause a shutdown of the ice stream. In a similar fashion, Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others (2000b) argued that the interplay of basal shearing, melting and sliding may lead to unstable stoppage of streaming flow. During the recovery phase, a positive mass budget enables ice-cap growth, while strain heating and geothermal heating gradually raise basal temperatures to melting point. This describes a full cycle with a periodicity of 1060 years.
Contour maps of (a) thickness, H, (b) basal velocity, vbase, (c) water saturation, w, and (d) basal temperature, T base, during one cycle of a low-frequency oscillation with (ϕ, q) = (7.5°, 0.3); t is given in years
Whereas in a high-frequency oscillation the oscillation period is determined by the response time of the ice mass to surface slope variations, in a low-frequency cycle the recovery time of the ice geometry and thermodynamics after a tongue advance is the decisive factor. Figure 6 illustrates the difference in volume fluctuations between high- and low-frequency oscillations. As in a high-frequency oscillation, the role of basal water feedback and local steepening in the low-frequency cycle is again to enhance acceleration and deceleration of ice flow, in that way keeping the ice cap out of a steady state. In contrast to high-frequency cyclicity, significant fluctuations in the extent of the temperate zone arise in low-frequency oscillations, thereby illustrating the important role of cold ice advection. The low-frequency cyclic mechanism described here is similar to the ‘binge/purge’ mechanism proposed by Reference MacAyealMacAyeal (1993) and the cyclic mechanisms described by Reference PaynePayne (1995) and Reference Dunse, Greve, Schuler and HagenDunse and others (2011).
3.1.4. Occurrence and magnitude of oscillations
The above description of the high- and low-frequency cycle mechanisms demonstrates the role of basal water feedback and local slope variations in maintaining oscillations. Figure 3 demonstrates that the presence and characteristics of the different types of oscillation depend strongly on the choice of parameters ϕ and q.
Lowering the till strength, ϕ, reduces the basal resistance experienced by the sliding ice at the ice/bed interface. As a result, sliding velocities increase and the sliding zone expands. In an environment with overall lower basal shear stresses, flow velocities are more sensitive to fluctuations of the driving stress. Under these conditions, perturbations of the surface slope have a larger impact on flow velocities. The sensitivity of flow velocity to local slope variations, further amplified by basal water feedback, is a measure of the amplitude of the oscillations. In this way, the transition from stationary fast flow to oscillatory flow for decreasing values of ϕ can be understood. Figure 3 also shows that high- frequency oscillation amplitudes are inversely proportional to the sliding exponent, q. This is explained by the fact that the basal shear stress, τb, is proportional to |u b| q , which implies that when q > 0 enhanced (reduced) sliding velocities experience enhanced (reduced) basal resistance. For larger values of q, this is a restrictive factor for the effectiveness of flow enhancement. Hence, a more plastic till promotes the occurrence and magnitude of high-frequency oscillations.
A transition to low-frequency cyclicity, in which the ice volume is strongly altered, occurs when local steepening and basal water feedback are strong enough to expand the temperate zone in the upstream direction. This may only arise in the case of intense flow enhancement, which occurs for very low values of ϕ. Expansion of the temperate zone facilitates a larger inflow of ice into the ablation area, which feeds the advancing tongue until cold-ice advection shuts down the enhanced inflow. The ‘activation’ of cold ice exposes a major difference between high- and low-frequency cyclic behaviour.
3.2. Climate parameter experiments
Next, the outcome of sensitivity experiments with a varying climate set-up is discussed. Starting from the output of the basal parameter experiment with (ϕ, q) = (12.5°, 0.0), 10 ka model runs with ELA, E, ranging from 1250 to 1450 m and 0°C surface-temperature height, z T0, ranging from —200 to 1000 m are performed. In all the experiments, the till strength and sliding exponent are fixed at 12.5° and 0.0, respectively. In reality, E and z T0 do not vary independently; however, for reasons of clarity the sensitivity to changes in E and z T0 is studied separately here.
Oscillation characteristics (period and amplitude of the basal velocity) are presented in Figure 8. The oscillation period does not show a clear trend with changes in E, but the oscillation amplitude is inversely proportional to E. Because the mass turnover depends on the height of the equilibrium line, the heat released by internal shearing becomes a more dominant component in the heat budget for low values of E. Mainly due to this effect, the available heat at the base increases for lower values of E, which leads to an increase in basal melting and subsequent expansion of the sliding zone. Lower values of E hence lead to an overall reduction of the basal resistance experienced by the ice at the ice/bedrock interface. As mentioned above, this favours oscillatory behaviour, since surface slope fluctuations have a stronger impact on flow velocities. In model runs with E ≥ 1400m, oscillations are absent. For values of E < 1250m, the ice tongue advances out of the model domain. Likely, a shift towards a low-frequency cycle will occur.
Oscillation characteristics (amplitude and period) as a function of (a) ELA, E, and (b) 0°C altitude, z T0, in the climate parameter experiments with (φ, q) = (12.5°, 0.0). The dotted line marks the standard set-up for E and z T0 used in the basal parameter experiments.
Changing the surface temperatures, through changes in z T0, indirectly affects multiple components in the heat budget. On the one hand, the surface temperature affects the ice viscosity and indirectly the ice volume, which has an impact on the amount of heat generation by internal shearing and trapping of geothermal heat. On the other hand, the surface temperature influences the amount of heat advection towards the base. The aforementioned effects seem to compensate for each other to a large extent for values of z T0 < 400 m, which implies that the amount of water production at the base, and thus the overall basal resistance, is rather constant. This is illustrated by a nearly invariant amplitude of the oscillations. For z T0 ≥ 400 m, more heat is available at the base, likely due to the reduced significance of cold-ice advection for higher surface temperatures. This reduces the overall basal resistance and favours stronger oscillations. From z T0 = 800 to 1000m, a transition from high- to low-frequency cyclic behaviour occurs.
3.3. Climate feedback experiments
In contrast to the previous experiments, in the following experiments the prescribed mass-balance and surface- temperature inputs are no longer time-independent. In the climate feedback experiments, the prescribed climate forcing is proportional to the evolving surface height, thereby introducing a time dependence. Starting without ice, 15 ka model runs were performed for multiple values of ϕ, while fixing q at 0.1. (Recall that in previous experiments the surface mass balance and temperature were a function of bedrock height rather than surface height.) Therefore, in all experiments a positive correction has been applied to the ELA, E, and 0°C altitude, z T0, equal to the spatial mean ice thickness in the corresponding basal parameter experiment, in which the surface height did not affect the climate forcing. In this way, the size of the ice cap and the mean surface temperature in the climate feedback experiments are not very different from those in the basal parameter experiments. Nevertheless, a direct (quantitative) comparison of oscillation amplitude and frequency with the outcome of the basal parameter experiments is not feasible.
3.3.1. Mass-balance/height feedback
The dependence of the mass balance on surface elevation provides an additional feedback mechanism which plays a role in the occurrence of oscillatory flow. Time series of ice volume and maximum basal velocities for five model runs with different values of ϕ (q = 0.1) are shown in Figure 9. Figure 10 shows maps of the evolution of the ice thickness and basal velocity at four points in time during one cycle in a run with ϕ = 15°. Only the tongue of the ice cap varies significantly in thickness during a cycle. As shown in Figure 9, t = t 2 and t 4 correspond to extrema in the ice volume. Extrema in the basal velocities are found during times of intermediate volume (t = t 1and t 3). At t = t 1, basal velocities are high since geothermal heating and internal shearing induce more basal melting after a period with large ice volume (Fig. 10). A large extent of the ice cap causes the total mass balance to be negative, inducing volume loss. Since the mass balance is surface-height dependent, thinning of the ice induces a more negative mass balance, causing more rapid thinning. In this way, the mass-balance/height feedback is ‘pushing’ the ice out of the mean geometric state. This feedback causes tongue thinning and terminus retreat to a minimum extent at t = t 2. As the tongue retreats, the total mass budget goes to zero, limiting further retreat of the terminus. Cooling of the thin ice and reduction of the sliding- zone extent result in lower sliding velocities and hence cause the mass budget to become positive after t = t 2. The resulting thickening is enhanced by the mass-balance/height feedback, and a delayed response of the ice area causes the mass budget to be positive at t = t 3. From t = t 3 to (4, terminus advance initiates and is further enhanced by higher sliding velocities as more heat is available for basal melting. The total mass budget reduces to zero at t = t 4, after which thinning sets in. This completes one cycle.
Time series of (a) ice volume and (b) maximum basal velocity in the mass-balance/height feedback experiments with q = 0.1 and φ ranging from 12.5 to 30.0°. At t 1 = 8000 years, t 2 = 8700 years, t 3 = 9400 years and t 4 = 10250 years, snapshots of the ice thickness in the run with φ = 15.0° (brown curves) are taken and shown in Figure 10.
Snapshots of (a) ice thickness and (b) maximum basal velocity in a run with (φ, q) = (15.0°, 0.1) in the massbalance/ height feedback experiments. The times t = t 1–t 4 refer to the times indicated in Figure 9. The dashed lines mark the position of the snout at t = t 1.
Figure 9 illustrates that high-frequency velocity oscillations may coexist with the low-frequency volume cycle described above. These high-frequency cycles seem to be of similar origin to the high-frequency oscillations in the basal parameter experiments. In line with the basal parameter experiments, the high-frequency oscillations only appear in experiments with low values of ϕ. Whether the low-frequency volume oscillations, forced by the mass-balance/height feedback, occur depends on whether a stable equilibrium state for the Fig. 10. Snapshots of (a) ice thickness and (b) maximum basal velocity in a run with ($, q) = (15.0°, 0.1) in the mass-balance/height feedback experiments. The times t = q-t4 refer to the times indicated in Figure 9. The dashed lines mark the position of the snout at t = t 1 ice cap exists. Under stable conditions, a perturbed surface profile will induce velocity variations that force the surface back to the equilibrium state. This relaxation effect is then stronger than the mass-balance/height feedback, which tends to enhance ice thickness perturbations. Unstable conditions may thus arise as the relaxation effect is dominated by the mass-balance effect. A strong dependence of the mass balance on altitude, i.e. a large balance gradient, directly enhances the strength of the mass-balance/height feedback and thus favours the occurrence of this type of oscillation. Furthermore, a low accumulation rate reduces the mass turnover, and thus the relaxation time, of the perturbed geometry, thereby increasing the likelihood of oscillations. Additional experiments with adjusted values of the balance gradient (not shown) demonstrate that the volume oscillations indeed vanish for small values of the balance gradient.
3.3.2 Surface-temperature/height feedback
Next, we introduce a surface height dependence (instead of a bedrock height dependence) of the surface temperature. Figure 11 shows time series of the ice volume and maximum sliding velocity for four model runs with a varying till strength. As mentioned above, a quantitative comparison with the basal parameter experiments is hindered by differences in the climate input (following the correction of E and z T0). Nevertheless, qualitatively comparing the results (cf. Figs 11 and 6) indicates very similar behaviour. When the till strength declines, a consecutive shift from stationary fast flow to high-frequency oscillatory flow and from high- to low-frequency cyclic flow occurs. For the chosen set-up, these results suggest that introducing a dependence of the surface temperature on the surface height does not affect the qualitative characteristics of oscillations.
Time series of (a) ice volume and (b) maximum basal velocity in the surface temperature/height feedback experiments with q = 0.1 and φ ranging from 10.0º to 22.5º.
3.4. Resolution experiments
To verify the robustness of the model under grid refinement, the mass-balance/height feedback run with (ϕ, q) = (15°, 0.1) was repeated with grid spacings of 0.83 km (standard/1.5) and 1.88 km (standard × 1.5). In the standard model run, we found two types of oscillatory flow, i.e. a low-frequency volume cycle and a high-frequency velocity oscillation (Fig. 9). Adaptive time-stepping assures stability for all experiments. Generally, the model time-step changes with grid resolution following Δt ~ Δx 2. In Table 1, mean values of the ice volume, ice thickness, ice area, sliding-zone extent and maximum basal velocity are given. Convergence towards a fine-grid result is clear for the mean ice volume and thickness. The mean volume shows a decreasing trend towards a higher resolution, with values in the 1.88 km run ~4% larger and in the 0.83 km run ~1% smaller than in the standard run. This mainly results from deviations of the mean ice thickness rather than the ice extent. The mean ice area does not converge; however, deviations from the standard run are small (up to 1.3% at 0.83 km resolution). The mean sliding-zone extent and maximum basal velocity converge under grid refinement, with absolute deviations from the standard run ranging from 2.4% (0.83 km) to 12% (1.88 km) and 1.6% (0.83 km) to 4.9% (1.88 km), respectively. Maps of the mean ice thickness and mean basal velocity for the three model runs are shown in Figure 12. Fluctuations of the ice volume show a decreasing but converging trend towards a finer grid, with an amplitude corresponding to 18% (1.88 km), 10% (1.25 km) and 8% (0.83 km) of the total volume. In terms of oscillation characteristics, the three model runs all show qualitatively similar behaviour, with a high-frequency surface oscillation superposed on a low-frequency volume cycle. Generally, we conclude that increasing the grid resolution leads to only minor changes in the model behaviour.
Mean volume, V mean, ice thickness, H mean, ice area, A mean, sliding-zone area, S mean, and maximum sliding velocity, v base,max, in the mass-balance/height feedback experiment with (φ, q) = (15.0º, 0.1) at resolutions of 0.83, 1.25 and 1.88 km. Values in parentheses indicate the positive (+) or negative (−) percentage deviations from the standard run
Contour maps of (a) mean ice thickness, H mean, and (b) maximum basal velocity, v base,max, in the mass-balance/height feedback experiment with (φ, q) = (15.0º, 0.1) at resolutions of 0.83 km (top), 1.25 km (middle) and 1.88 km (bottom).
4. Conclusion and Discussion
Model experiments were performed using PISM on a synthetic topography to study mechanisms causing oscillatory flow in a thermomechanically coupled polythermal ice- sheet model. Oscillatory behaviour in a slow-flow model occurs as internal feedback mechanisms, causing instability, force an ice mass out of an equilibrium state. One of the feedback mechanisms is basal water feedback, in which the interaction of basal sliding, melting and shearing may effectively enhance velocity fluctuations, thereby avoiding relaxation to a steady state. The basal parameter experiments illustrate that, for a fixed set-up, high-frequency oscillatory flow may initiate if the material till strength is lowered. Further reduction of the till strength amplifies oscillations and at some point leads to a transition to a low-frequency oscillatory flow type. High-frequency oscillations (period 114-169 years) and low-frequency cycles (period 1000+ years) differ significantly in terms of volume fluctuations. Variations in the extent of the basal area at pressure-melting point expose a control on the low- frequency cyclicity. The physical mechanism responsible for low-frequency oscillations in the basal parameter experiments is essentially similar to the oscillatory mechanism discussed by, for example, Reference Oerlemans and Van der VeenOerlemans and Van der Veen (1984), Reference MacAyealMacAyeal (1993), Reference PaynePayne (1995), Reference Calov, Ganopolski, Petoukhov, Claussen and GreveCalov and others (2002) and Reference Dunse, Greve, Schuler and HagenDunse and others (2011). In high-frequency oscillations, the basal water distribution fluctuates, while the temperate base is hardly affected. Variations in ice velocities induced by an oscillating surface slope are enhanced internally by basal water feedback. This is reflected in variations in the distribution of basal water within the temperate basal area. These oscillations therefore depend on the production and transport of basal water. Nevertheless, Reference Dunse, Greve, Schuler and HagenDunse and others (2011) also found high-frequency oscillations with SICOPOLIS in which sliding is directly linked to basal temperatures. A comparison of oscillations in SICOPOLIS and PISM on a similar set-up will give more insight into differences between simulated oscillation characteristics and the importance of linking sliding to evolving basal water. Low values of the sliding exponent, q, implying a nearly plastic till, favour the occurrence of oscillations driven by basal water feedback, as moderation of perturbed ice flow by related changes in the basal resistance is less pronounced.
Flow enhancement in the basal parameter experiments results from the interplay of basal shearing, water production and sliding. The high-frequency oscillations rely on the interaction of geometric variations with a redistribution of water on a temperate base, thereby affecting the till strength. This mechanism agrees, in a general sense, with the idea that changes in the ice thickness and movement of stored water in the subglacial till, affecting the rate of till failure, control surging behaviour of soft-bedded Alaskan-type glaciers (Reference Lingle and FatlandLingle and Fatland, 2003). We find that oscillatory flow is more likely to occur for lower material resistance of the bed, which is in line with findings that most, if not all, surge- type glaciers are underlain by a layer of easily deformable eroded material (Reference MurrayMurray and others, 2000; Reference Harrison and PostHarrison and Post, 2003). Due to the simplicity of the water transport model, which disregards the dependence of the drainage system’s transmissivity on water input, proposed mechanisms for surging based on transient behaviour of the drainage system’s efficiency are not modelled. Hence, simulating hard-bedded surging related to a transition between a slow distributed drainage system and a highly conductive conduit system (Reference KambKamb and others, 1985; Reference KambKamb, 1987), as well as soft- bedded surging associated with a destruction of the drainage system during acceleration (Reference Eisen, Harrison, Raymond, Echelmeyer, Bender and GordaEisen and others, 2005), would require a more detailed water transport model. However, thermally controlled soft-bedded surging (Reference Fowler, Straughan, Greve, Ehrentraut and WangFowler and others, 2001; Reference Murray, Strozzi, Luckman, Jiskoot and ChristakosMurray and others, 2003), often referred to as surges of Svalbard type, appears to be of similar origin to the low-frequency oscillations in this study. As in the modelled low-frequency oscillations, flow acceleration during a thermally controlled surge relies on a positive feedback between till weakening and sliding acceleration, whereas termination of flow acceleration is caused by basal refreezing. On a larger scale, Reference MacAyealMacAyeal (1993) used the same principles in a possible explanation for the Heinrich-type oscillations of the Laurentide ice sheet.
The frequency difference between the two oscillation types in the basal parameter experiments is mainly related to the impact of the oscillations on the ice volume and thermodynamic structure. In low-frequency oscillations, the ice volume and basal thermal regime vary significantly during an advance/retreat cycle. The length of the quiescent phase is inversely proportional to the accumulation rate (Reference Dowdeswell, Hodgkins, Nuttall, Hagen and HamiltonDowdeswell and others, 1995; Reference Eisen, Harrison and RaymondEisen and others, 2001), which controls the rate of recovery. A more dramatic advance requires a longer recovery time, thereby reducing the oscillation frequency. The periodicity of thermally controlled cyclicity therefore depends strongly on the chosen domain size. Since the experimental set-up in this study is not designed to study a specific type of observed cyclic behaviour, we are not able to compare simulated timescales of the low-frequency oscillations (1000+ years) with the observed periodicity of Svalbard-type surging (50-500 years) (Reference Dowdeswell, Hamilton and HagenDowdeswell and others, 1991) and the frequency of Heinrich-type oscillations (~7000 years). In low-frequency oscillations, the termination of sliding by refreezing of the bed is based on a similar feedback mechanism to that proposed for the stoppage of Kamb Ice Stream (Reference Tulaczyk, Kamb and EngelhardtTulaczyk and others, 2000b). Despite differences in the experimental setup, simulated timescales for cessation of streaming flow in this study correspond well with rates found by Reference Bougamont, Tulaczyk and JoughinBougamont and others (2003). In high-frequency oscillations, the period of oscillations is determined by the efficiency of flow enhancement in response to variations of the surface slope rather than recovery of the ice volume and thermodynamic structure. Frequency is therefore determined by the rate of flow enhancement, the area affected by flow enhancement and the mass turnover, affecting the relaxation time of a perturbed surface (Reference Jóhannesson, Raymond, Waddington and OerlemansJóhannesson and others, 1989).
Climate parameter tests demonstrate the high sensitivity of oscillatory behaviour to changes in the ELA and the surface temperature. Introducing a surface-height dependence of the mass balance may lead to another type of instability, in which surface perturbations are enhanced, causing sustained long-term fluctuations in the ice volume. These oscillations become more pronounced with a larger balance gradient and a low accumulation rate. Incorporating a height dependence of the surface temperature does not affect the qualitative aspects of cyclic flow in the chosen set-up of the experiments.
Numerical inconsistencies related to a discontinuous transition of ice flow at the ice-stream margin, as discussed by Reference Fowler, Straughan, Greve, Ehrentraut and WangFowler (2001) and Reference Bueler and BrownBueler and Brown (2009), are tackled in PISM by the hybrid implementation of stress equations, causing a smooth transition from non-sliding to sliding ice flow. Results of ISMIP-HEINO (Reference Greve, Takahama and CalovGreve and others, 2006; Reference CalovCalov and others, 2010) have demonstrated that oscillatory flow is not confined to models with a sharp transition from sliding to non-sliding flow. In this study, a physical description of the cyclic flow has been given and oscillations do not seem to arise from unreasonable physics or numerical artefacts. All experiments have been performed on a single bedrock topography with fixed dimensions. It would be worthwhile to investigate how a different basal set-up affects the oscillation characteristics. Furthermore, in all experiments in this study the base was polythermal, with a temperate zone downstream of a cold-based region. It would be interesting to also study flow patterns for different polythermal structures and for fully temperate ice. On a fully temperate base, low-frequency cycles are likely to be absent, since they rely on variations in the extent of the temperate base. Additionally, extending the analysis to tidewater glaciers would be a worthwhile experiment. Finally, itshould be stressed that the prescribed external climate forcing is constant in time. Hence, possible interactions of climatic fluctuations and internal flow instabilities are not considered.
The coupling of stress balances and the interplay between geometry evolution, thermodynamics and hydrology provide a good basis for analysis of cyclic flow. The diffusive basal water model used in this study does interact dynamically with changes in the heat budget near the base, the basal shear stress and sliding velocities, which makes it suitable for studying cyclic behaviour resulting from these interactions. Nevertheless, inclusion of a more realistic water transport model, simulating flow from high to low fluid potential with a transient hydraulic conductivity, would be an improvement (Reference Flowers and ClarkeFlowers and Clarke, 2002). Ultimately, incorporating a water transport scheme which accounts for fluctuations in the drainage efficiency associated with changes in the morphology of the hydrological system is required to study Alaskan-type cyclic behaviour.
Acknowledgements
This project is funded through the Royal Netherlands Academy of Arts and Sciences (KNAW) professorship of J. Oerlemans. Special thanks are due to C.H. Reijmer and E. Bueler for useful comments, which helped to improve the manuscript. We thank two anonymous reviewers for constructive and useful comments.
