Uncertainty in ice-sheet models

Uncertainty in ice-sheet models has two main sources. Firstly, the modelled ice sheet is a function of the initial decisions made in abstracting the real system to a conceptual model. By abstracting reality, we are accepting that uncertainty will be inherent in the results of our model – users of numerical modelswill seektominimalizethis. Forexamplea model which does not consider longitudinal stresses will fail to replicate particular aspects of the glacial system. Secondly, model uncertainty can arise from uncertainties in the input parameters, such as the choice of parameter values at model initialization. The first set of uncertainties can only be explored through comparisons of systems that use different levels of abstraction in constructing their model rules. Thus, for example Reference Hubbard, Hein, Kaplan, Hulton and GlasserHubbard and others (2005) showed how the inclusion of longitudinal stresses led to considerable variation in modelled ice-sheet extents in Patagonia compared to previous results (Reference Hulton, Sugden, Payne and ClappertonHulton and others, 1994, Reference Hulton, Purves, McCulloch, Sugden and Bentley2002) based around models using the shallow-ice approximation (Reference NyeNye, 1957).

Within the ice-sheet modelling community, intercomparison experiments have sought to explore the sensitivity of modelled results to differences in implementation (and in principle abstraction). For example, the European Ice Sheet Modelling Initiative (EISMINT) experiments (Reference Huybrechts and PayneHuybrechts and others, 1996; Reference PaynePayne and others, 2000) explored a variety of model types and their responses to a range of glaciological (e.g. steady-state Greenland ice sheet under present conditions) and hypothetical experiments (e.g. ice-stream formation on a hypothetical symmetric topography). In many modelling fields a similar approach, though motivated by a different philosophy, is applied in ensemble modelling. Here, a variety of models are usually run in parallel and the results are weighted to provide a probabilistic forecast (Reference AndersonAnderson, 1996). Such approaches are arguably better suited for use in policy formulation, where the focus of interest is centred on risk estimation (Reference Vaughan and SpougeVaughan and Spouge, 2002).

In general, such model intercomparisons allow exploration of differences between models for a given scenario, and are important in moving models towards a consensus view. However, since intercomparisons are also often used as a benchmark towards which modellers aim in development, they cannot be said to be a means of quantifying uncertainty. Additionally, sensitivity to ‘external’ (input) parameters is not usually the subject of those intercomparisons, leaving a gap in the understanding of uncertainties associated with ice-sheet modelling, which we aim to fill. This second set of uncertainties, the focus of this paper, can be explored in a variety of ways. The most common approach employed in glaciological modelling is to carry out a sensitivity study. In such studies, all but one parameter are typically held constant whilst the parameter under investigation is varied through some range, and the response of the modelled ice sheet is visualized or quantified (e.g. Reference Van de Wal and OerlemansVan de Wal and Oerlemans, 1994; Reference Fabre, Letréguilly, Ritz and MangeneyFabre and others, 1995; Reference Ritz, Fabre and LetréguillyRitz and others, 1997; Reference Huybrechts and de WoldeHuybrechts and de Wolde, 1999; Reference Purves and HultonPurves and Hulton, 2000; Reference PattynPattyn, 2003; Reference Essery and EtcheversEssery and Etchevers, 2004). Sensitivity studies are important because they allow modellers to quantify which parameters, in isolation, have the greatest influence on model results. However, since all but the most trivial ice-sheet models are inherently non-linear, it is often difficult to properly quantify uncertainty or allow assignment of probabilities to particular outcomes (Reference Van der VeenVan der Veen, 2002).

Parametric uncertainty analysis

Currently, despite the increased importance of numerical modelling in policy making, relatively little research has explored ways in which uncertainties might be consistently calculated in individual ice-sheet model results. Parametric uncertainty analysis (PUA), where all parameters are varied together, is one such approach. Marshall and others (2002) investigated the response of the North America ice sheet to a range of parameters, and completed a number of simulations where parameters were varied simultaneously. Reference Tarasov and PeltierTarasov and Peltier (2004) also explored the response of the North American ice sheet during the LGM to simultaneous variations in a range of parameters, whilst geophysically constraining the modelled ice extents. However, in both these experiments the variation of parameter values appears to have been based on a limited variation of parameter values over a predefined range. Such approaches allow the range of potential model responses to input parameters to be estimated, but these responses are not associated with probabilities.

Reference Vaughan and SpougeVaughan and Spouge (2002) carried out a full parametric uncertainty analysis of future ice-sheet mass balance for the Greenland ice sheet using parameterizations of accumulation and ablation. A key difference in this approach is that probability density functions (PDFs) are estimated for each parameter, and used in the model runs to allow calculation of a final PDF for a chosen model output parameter. PDFs are, in general, used to associate a probability with a certain parameter value or model outcome. For example using a normal distribution as a PDF, values close to the mean are more likely to occur than values that deviate more from the mean – in the case of a normal distribution, 66% of all values are within one standard deviation of the mean. The advantage of using PDFs in parametric uncertainty analysis is that by including a (realistic) likelihood for a certain parameter value to occur, probabilities can be assigned to model outcomes. This, in turn, allows modellers, policy makers and users to identify in a more meaningful way how each parameter might affect the result, and estimate, for instance, the risk of extreme events. Typically, PUA is carried out using Monte Carlo simulation (MCS) which is a brute-force technique, where large numbers of model runs are normally required to generate an ensemble of results from which PDFs can be calculated.

Topography in ice-sheet models

In ice-sheet modelling the bed topography on which models are run is generally treated as having no uncertainty associated with it. However, the data used in the generation of these bed topographies is often associated with large uncertainties, which vary as a function of locale (Hebeler and Purves, 2008). Furthermore, the interpolation methods used to resample these data to resolutions appropriate for typical ISMs introduce further uncertainties (Reference Hebeler and PurvesHebeler and Purves, 2004). Topography, as discussed by Reference KerrKerr (1993), plays a key role in influencing ice-sheet behaviour through its geometry, whereby the mass-balance profile is controlled by relief (e.g. through orographic effects) and the topographic geometry (e.g. hypsometry). Topography also influences ice-sheet configuration directly by facilitating or constraining ice flow, for example in a glacial trough (Reference Jamieson, Sugden, Cooper, Barrett, Stagg, Storey, Stump and WiseJamieson and others, 2008). The sensitivity of modelled ice-sheet extents to uncertainty in topography is, thus, relevant to furthering our understanding of ice sheets, and to refining numerical models that seek to reconstruct such systems.

Aims

The objectives of this paper are, therefore, threefold:

1. 1. to propose a generalizable set of methods for exploring the parametric uncertainty of ice-sheet model results;

2. 2. to explore the influence of uncertainty in bed topography on modelled ice extents and volumes;

3. 3. to compare the influence of topographic uncertainty to other sources of uncertainty (such as climatic or basal model inputs) on modelled ice extents and volumes.

The remainder of the paper is structured as follows. We briefly introduce the ice-sheet model used in our experiments, before we describe a method for exploring uncertainty in bed topography. We then review the parameters and associated PDFs derived for the parametric uncertainty experiments. The results of the experiments are presented as plots of ice extent and volume through time and by visualizations of ice extent and volume from ice sheets in equilibrium, using probability maps. Additionally histograms of the associated PDFs are given. Finally, we discuss the significance of the results and their implications for those constructing, using and interpreting ice-sheet models.

The GLIMMER ice-sheet model

Physical laws relating to the dynamics of glaciers have long been incorporated into numerical models of ice flow, in order to make predictions about past, present and future ice-sheet behaviour (e.g. Reference Hulton, Sugden, Payne and ClappertonHulton and others, 1994; Reference OerlemansOerlemans and others, 1998; Reference Huybrechts and Le MeurHuybrechts and Le Meur, 1999). We aim to test how uncertainty in various model input parameters can affect the results of such ice-sheet models. To achieve this, we employ the community ice-sheet model GLIMMER (General Land Ice Model for Multiply Enabled Regions; Reference PaynePayne, 1999) which builds on the foundations laid down by the modelling studies of Reference HuybrechtsHuybrechts (1986), Reference Boulton and PayneBoulton and Payne (1993) and Reference Payne and DongelmansPayne and Dongelmans (1997). Ice dynamics calculations implement the shallow-ice approximation (Reference HutterHutter, 1983), a widely used approach which assumes that bedrock and ice surface slopes are small enough that normal stress components can be neglected. The ice flow law parameter is handled by a three-dimensional thermomechanical model. A full outline of the numerics implemented in GLIMMER is provided by M.K.M. Hagdorn and others (http://glimmer.forge.nesc.ac.uk), and Reference PatersonPaterson (1994) provides additional context to the derivation of these mechanics.

The ability of a modelled ice sheet to grow and interact with its bed and climate is strongly dependent on basal ice velocities. In GLIMMER, the velocities of the basal ice, v _b, at any one point are determined as a function of the basal shear stress, τ_b, and the effective pressure (N = ice overburden pressure – basal water pressure) thus:

where k is a constant describing the thermomechanical properties of ice and p and q are positive integers (Reference PatersonPaterson, 1994). The basal slip coefficient, t _b (also known as the basal traction parameter), determines the dependence of the thermomechanical properties of ice on N and therefore assumes that k and N ^−q can be integrated as a single parameter. Equation (1) can therefore be simplified (M.K.M. Hagdorn and others, http://glimmer.forge.nesc.ac.uk) to:

Sliding is assumed to occur when ice temperature at the bed reaches pressure-melting point (taking into account both a user-defined geothermal heat flux, G _therm, and calculated frictional heat contributions), at which point the ice can detach from the bed and begin sliding due to the presence of water. Its ability to do so is a function of t _b, which is calculated as a linear function of basal melt rate, b _melt:

The slope of the function is given by t _bslope, the softness of the bed is given by a parameter, b _soft, and the basal traction is limited from becoming too large by the value of t _bmax. Variability in the basal traction parameter will therefore alter the behaviour and form of the resulting ice mass.

At each time-step, the model requires spatial parameterization of surface mass balance and air temperature as inputs. We follow the approach of Jamieson and Sugden (2008) to generate these inputs. The distinction between snowfall (which only occurs below a given temperature, t _snow) and rainfall across the model domain is determined using a positive degree-day (PDD) model (Reference ReehReeh, 1991). The PDD approach assumes that melt at the ice or land surface is proportional to the number of days (integrated through time) in which the air temperature, T, rises above freezing point. The number of PDDs is, therefore, proportional to the energy available for melting. Melt, w, is calculated by:

where DDF is the degree-day factor describing the density and albedo of snow or ice (DDF is different for each). The melt calculation further assumes that there is an annual sinusoidal cyclicity (i.e. seasonality) to air temperature, which is calculated by taking the mean annual temperature and an associated half-range, t _range. Further diurnal deviations from the sinusoidally shifted mean annual temperature are accounted for by assuming that this variability has a normal distribution with a standard deviation of 5°C. The mean annual air temperature, MAAT, at the land or ice surface is calculated given a MAAT at sea level which is then adjusted to the surface using an atmospheric lapse rate, lrate. The model also incorporates a firn model to account for the fraction, w _max, of melted snow that refreezes to become superimposed ice. The patterns of precipitation and MAAT used as inputs to this model are described under ‘Climate forcing’, below.

The Earth can be approximated as a thin elastic layer (the lithosphere) that floats above the asthenosphere. GLIMMER treats these two components separately and can do so in a number of ways. The flexural rigidity of the lithosphere is modelled as it responds to changes in load (Lambeck and Nakiboglu, 1980). The response of the asthenosphere to changes in overburden is modelled so mantle flow is approximated over time as an exponentially decaying hydrostatic response function (M.K.M. Hagdorn and others, http://glimmer.forge.nesc.ac.uk).

The model has been tested in the European Ice-Sheet Modelling Initiative (EISMINT) experiments (Reference Huybrechts and PayneHuybrechts and others, 1996; Reference PaynePayne and others, 2000) and it has been employed to gain better understanding of ice dynamics (e.g. Reference Le Brocq, Payne and SiegertLe Brocq and others, 2006). Furthermore, it has been used as a tool to reconstruct ice-sheet configurations over numerous regions and time periods (Reference PaynePayne, 1999; Reference Payne and BaldwinPayne and Baldwin, 1999; Jamieson and Sugden, 2008; Reference Lunt, Valdes, Haywood and RuttLunt and others, 2008).

Study area

For our experiments, we use basal topography of Fennoscandia extracted from 30 arcsec resolution GLOBE (Global Land One-Km Base Elevation) DEM data (http://www.ngdc.noaa.gov/mgg/topo/globe.html) between 50 and 72° N and between 14° W and 56° E. This is merged with ETOPO2 bathymetric data (http://www.ngdc.noaa.gov/mgg/global) at 2 arcmin resolution and transformed to Albers equal area projection. Thus the open water-bodies in GLOBE are filled by ETOPO2 data, and resampled to 10 km resolution. The model domain covers an area of 2580 km by 2370 km, with a minimum elevation of −3882 m and a maximum elevation of 1960 m.

Climate forcing

The primary objective of the experiments described in this paper is to investigate the influence of uncertainty in topography, climate and model parameters upon modelled ice-sheet extent and volumes. To this end, we are concerned with generating ice sheets in which parameters are varied and which can be compared to a benchmark simulation of a Fennoscandian ice sheet.

For precipitation rates, Climatic Research Unit (CRU; http://www.cru.uea.ac.uk) and Intergovernmental Panel on Climate Change (IPCC; http://www.ipcc-data.org/) observation data of recent climate were used to interpolate a precipitation scheme (Fig. 1), which can be scaled and varied over time. Notable features are high values of precipitation along the western coastal regions, with the highest precipitation of ∼2.6 m a⁻¹ occurring in the southwest fjords.

Fig. 1.

Precipitation scheme for Fennoscandia compiled from present-day CRU and IPCC observation data serving as baseline input for the GLIMMER ISM.

CRU and IPCC climate data from 30 year observation periods were also used to compile MAAT at sea level for the study area. From this data, latitudinal dependencies were derived to construct an input temperature field at sea level that matches the general trends from observation data and gets colder with increased latitude. For modelling steady-state ice sheets on Fennoscandia, recent mean temperature was lowered by a constant value over the whole modelling domain (cf. Reference Huybrechts and de WoldeHuybrechts and de Wolde, 1999; Reference Forsström and GreveForsström and Greve, 2004).

Preliminary sensitivity tests for different temperature and precipitation schemes were run to determine an optimal climate set-up for the DEM uncertainty analysis as well as the parametric uncertainty tests. Lowering of present-day temperature by 10–14 °C, approximately resembling the cooling during the Last Glacial Maximum (Reference Ritz, Fabre and LetréguillyRitz and others, 1997), proved to produce stable and sensible ice-sheet extents.

DEM uncertainty experiments

DEMs are subject to uncertainties generated from sources such as resampling and data accuracy. GLOBE DEM data, commonly used for modelling in Fennoscandia, are known to contain substantial uncertainties originating from numerous data sources, from compilation methods and from measurement errors. Despite claims as to the need for detailed error models as an integral part of digital elevation data (Reference Ehlschlaeger and GoodchildEhlschlaeger and Goodchild, 1994), common DEMs such as GLOBE are distributed with only global error or accuracy figures (Reference FisherFisher, 1998) as described by the root-meansquared error, rmse (D.A. Hastings and others, http://www.ngdc.noaa.gov/mgg/topo/globe.html). These global accuracy figures have proven unrealistic and are of limited use as they lack information on the spatial distribution of error which is often spatially correlated with topographic attributes such as altitude or slope (Reference Holmes, Chadwick and KyriakidisHolmes and others, 2000; Reference Oksanen and SarjakoskiOksanen and Sarjakoski, 2005).

Therefore, in order to realistically model DEM uncertainty including spatial dependencies, higher-order reference data are necessary where no explicit error model exists. Because of the limited availability of such reference data for large regions of Fennoscandia, a model to simulate GLOBE DEM uncertainty has been developed (Hebeler and Purves, in press a,b). This model assesses GLOBE DEM error properties over a number of mountainous regions where higher-order reference data (A. Jarvis and others, Shuttle Radar Topography Mission (SRTM), http://srtm.csi.cgiar.org) are available. Using these regions, characteristics of error magnitude and spatial configuration, including their correlation with DEM characteristics such as elevation, slope and roughness, have been assessed. Analysis showed GLOBE error to be best modelled using deterministic components of uncertainty, modelled using regression, combined with stochastic elements. While the deterministic components reproduced the amount and spatial configuration of uncertainty well, the stochastic elements allow the produced uncertainty surfaces to be used within Monte Carlo simulations. The resulting uncertainty surfaces are essentially surfaces of deviations from the original elevation data in metres, and are added to the GLOBE DEM.

Using this uncertainty model, derived from areas with existing reference data, a suite of 150 uncertainty surfaces at 1 km resolution for Fennoscandia were produced. These uncertainty surfaces were then added to the original GLOBE DEM, resampled to 10 km and used as input to the GLIMMER ISM for Monte Carlo simulations, thus allowing assessment of the impact of GLOBE DEM uncertainty on ISM results. The baseline model set-up was used in our simulations, with present MAAT lowered by −10, −12 and −14°C, producing three sets of ISM runs. In order to determine the number of runs necessary to obtain stable and reliable results, the standard deviation of modelled ice extent across all model runs, and the standard deviation of this standard deviation, were plotted against the increasing number of runs, until both measures stabilized towards a constant value (Reference Raaflaub and CollinsRaaflaub and Collins, 2006). For each set, the mean and standard deviation of modelled ice extent and volume after 30 kyr (model years) across all runs were evaluated, and probability maps (Hunter and Goodchild, 1995) were compiled, as a measure of the frequency with which each cell was glaciated across a suite of Monte Carlo simulation (MCS) runs.

Parameter sensitivity

Using the PDD approach for determining snowmelt in GLIMMER, five parameters influence the temperature forcing. These are the mean annual air temperature at sea level, the latitudinal dependency of temperature, the seasonal variation of temperature, the atmospheric lapse rate and the geothermal heat flux. Additionally, four parameters influence the ISM via their explicit or implicit dependence on temperature, namely the degree-day factors for ice and snow, the threshold temperature for precipitation to fall as snow and the meltwater refreezing fraction (Fig. 2).

Fig. 2.

Relational diagram of the ISM parameters used in the sensitivity and parametric uncertainty analyses.

The following parameters were thus selected for initial sensitivity testing in order to identify sensible parameter ranges for PUA:

• Mean annual air temperature and latitudinal dependency: MAAT

• Atmospheric lapse rate: lrate

• Annual temperature half-range (seasonal variation): t _range

• Precipitation: precip

• Degree-day factor for ice: DDF_ice

• Degree-day factor for snow: DDF_snow

• Geothermal heat flux: G _therm

• Threshold temperature for precipitation to fall as snow: t _snow

• Meltwater refreezing fraction: w _max.

The modelling of basal water and basal traction can influence ice-sheet configuration, because of their importance in controlling basal ice velocities. A flow enhancement factor is often used in ISMs as a tuning parameter to adjust ice-flow properties, for example, for layers of soft or warm ice. Thus, the following three additional parameters were tested, giving a total of 12:

• Flow enhancement factor: ffac

• Ice-thickness threshold for resolving physics: icelimit

• Basal traction rate: btrc.

The influence of each of these parameters on mass balance and ice dynamics and their interrelations within GLIMMER are schematically depicted in Figure 2.

For all parameters, initial values, as well as variation ranges, were derived from the literature and are described in detail in the following section. Additionally, for climate forcing, variation was derived from uncertainty in the input data (e.g. using differences between CRU and IPCC data) and by accounting for uncertainty given in the respective metadata. Using information from the literature along with details of input data uncertainty and sensitivity, a probability density function (PDF) was derived for each parameter, to be used in the PUA. For this analysis, a suite of 510 input configurations for GLIMMER was created, where each of the 12 parameters listed above was randomly assigned a value according to the corresponding PDF.

Parameter ranges and PDFs

For each parameter the value range and selection criteria used within each sensitivity test are given here. A PDF to be used within the PUA for each parameter is deduced using the input criteria as well as the model sensitivity. An overview of both the tested parameter range and the resulting model sensitivity is given at the end of this section.

Results and discussion of sensitivity tests

The sensitivity tests conducted prior to the PUA indicate that the modelled ISM extent and volumes are most sensitive to input parameters that influence mass balance (see Fig. 2; Table 1): the tested range of mean annual air temperature (Fig. 3a) and seasonal variation in air temperature (Fig. 3b) resulted in the largest variations in modelled extent (up to 200%) in the sensitivity tests. Up to a 100% change was recorded for variation in precipitation for the −12 and −14°C scenarios (Fig. 3e and f). The impact of lapse rate (Fig. 3c), precipitation at −10°C (Fig. 3d) and the DDFs of ice and snow (Fig. 3g and h) was in the range 40–75%.

Table 1.

Sensitivity test parameters: default values of the baseline approach, tested parameter range and impact on modelled ice extent and volume. Impact is measured as the deviation range from the baseline (default) value after 30 kyr in per cent. Ratio is the ratio of relative impact of parameter variation on extent when compared to volume. Precipitation (m a⁻¹) for each of the three temperature scenarios is varied in per cent relative to the baseline input scenario. Maximum basal traction is modelled using a slope of 0.02 for the basal melt function. Results of the DEM uncertainty test shown for comparison, with local maximum and minimum modelled uncertainty in metres, are discussed in the section ‘DEM uncertainty tests’

Model parameters that control ice dynamics have less influence over ice-sheet model results: variation of basal traction (Fig. 3l), the flow enhancement factor (Fig. 3m) and the refreezing factor (Fig. 3k) had an impact on modelled ice extent of ∼10–20%. Additionally, the variation of both the geothermal heat flux (Fig. 3i) and snow threshold temperature (Fig. 3j) parameters within their tested ranges was found to have limited influence on both modelled ice extent and volume (∼1–2%). Varying the ice-thickness threshold (Fig. 3n) for resolving of dynamics had a low to medium non-linear impact on model results, which varied by ∼5%. In general, the impact of changing the above parameters is slightly smaller on extent than volume.

While ‘external’ climate parameters have a greater influence on ice-sheet configurations, the ‘internal’ parameters that control ice dynamics can influence model stability, as shown by our sensitivity tests for different basal traction configurations. The effect these parameters have on model stability can be observed when examining ice-sheet behaviour. For example, particular selections of parameter values can generate oscillations in ice flow. Therefore parameter setups have been chosen that are unlikely to trigger model instabilities, and because of the use of these more conservative value ranges, the influence of internal ice-dynamics model parameters associated with internal ice dynamics might be low compared to the climate-forcing module.

While the strong influence of climate parameters concurs with previous findings (Reference Huybrechts and de WoldeHuybrechts and de Wolde, 1999), the sensitivity of modelled ice sheets to some other parameters is more controversial. For instance, Reference Huybrechts and de WoldeHuybrechts and de Wolde (1999) found basal melt rates had a large influence on the mass balance of the Antarctic and Greenland ice sheets, and melt rates of 10–30 m a⁻¹ have been found for ice shelves (Reference Huybrechts and de WoldeHuybrechts and de Wolde, 1999; Reference Rignot and ThomasRignot and Thomas, 2002). For grounded ice, basal melt rates are heavily determined by the geothermal heat flux, and high basal melt rates under the Antarctic and Greenland ice sheets are explained by multiples of the normally assumed fluxes of 56 mW m⁻². Reference Fahnestock, Abdalati, Joughin, Brozena and GogineniFahnestock and others (2001) report local basal melt rates of up to 1 m a⁻¹ for the Greenland ice sheet. In our experiments, however, maximum basal melt rates lie in the range 0.5–1 m a⁻¹, and variation of the geothermal heat flux showed little impact on modelled ice extent or volume, in agreement with the findings of Reference Ritz, Fabre and LetréguillyRitz and others (1997). This inconsistency could be explained by the conservative selection of parameter variation values, but since the applied values are averages over the whole modelling domain, higher rates appeared to be unreasonable and unsupported by the literature. Another possible explanation is that in our model set-ups, ice-thickness variation (and therefore topography) is a dominant factor influencing basal melt rates through overburden pressure. Because rapid variations in ice thickness are driving large differences in overburden pressure, they may overlay the impact of changes in geothermal heat flux. In this case, higher variations in G _therm than the ones applied in our studies (which we assumed to be unrealistic) would show a greater impact on modelled ice extents and volumes.

Our experiments showed that increasing flow enhancement factors had a considerable effect on reducing modelled ice volume and, to a lesser amount, ice extent. In contrast, Reference Ritz, Fabre and LetréguillyRitz and others (1997) found a slight increase in modelled ice extent for initial increases of the flow factor for their Greenland experiments. Similar behaviour can be observed for variation in basal sliding comparing the two sets of experiments. Possible explanations may lie both in the topography and the climate configurations. For example, the Greenland ice-sheet margin positions are not only limited through high ablation in lower areas, but also by mass transport from the accumulation areas. In this case, the lower ice viscosity simulated by a higher flow factor would allow ice to flow further without being ablated. The relatively large seasonal temperature variations in our model set-up make it likely for an area to experience ablation during some time of the year, especially for lower elevations. Therefore increases in ice extent during colder periods of the year may ablate during the summer, and not contribute to the long-term expansion of the ice sheet.

The baseline parameters and their range examined in our sensitivity tests have been carefully compiled. However, both parameter ranges and parameter PDFs are approximations or are still based on assumptions. The impacts on modelled ice extent and volume of varying parameters across the ranges given in Table 1 are therefore not applicable for a direct comparison with each other or results from other studies. Nevertheless, the impact and the impact ratio (impact on extent/volume, measured as deviation from the mean in per cent) for each parameter can be used to compare the relative influence of each of the factors on ISM results. The derived impact ratio (Table 1) can also be used to determine whether parameters have a predominant influence on either modelled ice extent or volume. This is an important consideration in understanding the mechanisms through which these parameters influence the ISM.

In general, the impact ratio of parameter variation on modelled ice extent and volume is ∼0.8 (Table 1). This slightly higher impact of uncertainty on modelled ice volume can be explained by the fact that ice extent is a two-dimensional measurement. By contrast, the calculation of ice volume is three-dimensional and the variation of ice thickness is explicitly included. Because of the influence of topographic uncertainty on the calculation of ice thickness, its direct impact is likely to be reflected in modelled ice volume. Parameters with a low impact ratio have a greater influence on modelled ice volume than extent, such as the flow enhancement factor and the geothermal heat flux, via its influence on basal sliding. Both parameters influence ice velocity, through changes in either ice viscosity or sliding behaviour. Increased ice velocities support the accelerated flow of ice to lower elevations, where it is more likely to be ablated. At the same time, because ice flows away from the accumulation areas faster, ice build-up is inhibited. Thus, if mass-balance gradients at the edges of the ice mass are very steeply negative (i.e. high ablation rates occur near the equilibrium-line altitude (ELA)) then ice extent may not change significantly with increasing ice velocities. At the same time ice thickness (and thus volume) decreases because the ice transports ice from the accumulation zone to the ablation zone more rapidly.

Conversely, the maximum basal traction rate and the threshold temperature for precipitation to fall as snow have impact ratios close to 1 (Table 1). This means their influence on modelled ice extent is larger than on volume, when compared to other parameters. While the overall impact of t _snow might be too small to support detailed analysis, it is interesting that higher rates of basal traction appear to favour changes in ice extent rather than volume. This is opposite to the effect of geothermal heat flux, a parameter that controls basal sliding via its dependence on basal melt. One explanation for this possible contradiction could lie in the fact that variation in geothermal heat flux affects the entire ice sheet, resulting in an overall change in the area of basal ice at its melting point. Changes in basal traction rate, however, only influence areas that already are at the pressure-melting point and are experiencing basal sliding. As this phenomenon occurs mostly in stream patterns and towards the edges of the ice sheet, where the ELA is positioned, changes in ice volume generated by the spatial pattern of ice velocity through basal traction may be relatively slow. At the same time, ice extent may be influenced more strongly through local advances of ice-stream outlets, because concentrated patches of thicker ice do not melt away as quickly as homogeneously distributed, but thinner, slow-moving ice advance. However, differences in, as well as overall values of, basal traction rate impact on modelled ice extent and volume are small, which means that small relative changes in model results can have strong influences on the associated impact ratio.

It is interesting to note the observable decrease of the impact ratio (Table 1) for the precipitation variation at −14°C. In comparison to the warmer climates tested, the impact upon volume versus extent increases. This may be due to the impact of air temperature on ice temperature and thus velocity. For the cold, −14°C, scenario, the ice-sheet boundary is strongly influenced by (low) ice velocities caused by the low temperatures as well as the high ice thickness dampening steeper slopes of the bedrock topography. This in turn results in areas of positive mass balance with increasing ice thickness, where ice cannot flow to ablation areas, so direct changes in mass balance mainly act on ice volume, rather than extent. Effectively, the stiffer, slower-flowing ice results in the build-up of a higher ice dome than the higher-temperature scenarios.