Effect of uncertainty in surface mass balance--elevation feedback on projections of the future sea level contribution of the Greenland ice sheet

Effect of uncertainty in surface mass balance--elevation feedback on projections of the future sea level contribution of the Greenland ice sheet


Introduction
The Greenland ice sheet (GrIS) response to climate change has two parts: surface mass balance (SMB), which is the sum of surface accumulation and surface ablation (broadly speaking, the balance of snowfall versus meltwater runoff); and dynamic, the changes in ice flow and discharge from the ice sheet.Various approaches to simulating these have been taken for making projections of the GrIS contribution to sea level.As with all simulation problems, there is a trade-off between representing more processes, with the aim of increasing the physical realism of the simulations, and technical and computational resources, which limit implementation and the number of simulations.
The dynamic response is simulated with ice sheet models (ISMs), which solve the Stokes equations in complete or approximate form.SMB can be simulated with sophisticated, physically based energy balance schemes in regional climate models (RCMs) such as MAR (Modèle Atmosphérique Régional: Fettweis, 2007) and RACMO2/GR (e.g.Ettema et al., 2009).The two aspects of the GrIS response can thus be combined by forcing an ISM with SMB simulated by an RCM.But RCMs usually use a fixed surface topography, neglecting the important effects of ice sheet surface elevation changes on the atmosphere (Edwards et al., 2014) by omitting the SMB-elevation feedback.An alternative to using SMB from an RCM is to simulate it within the ISM, so that the evolving ice sheet topography can dynamically alter the SMB.The most common approach to this is with a positive degree day (PDD) scheme, which parameterises SMB as a function of temperature and precipitation (supplied from observations for the present day, or a regional or global climate model for future projections) and possibly applies a simple snowpack model (e.g.Janssens and Huybrechts, 2000).PDD models incorporate the temperature aspect of the SMB-elevation feedback through a lapse rate correction, but not the precipitation aspect (except, in some cases, through a scaling factor for temperature), and also represent SMB much more simply than RCMs such as MAR and RACMO2/GR (Edwards et al., 2014).
To make the best use of physically based simulations of the GrIS response to climate change -simulating ice flow with an ISM and SMB with an RCM, while also including the SMB-elevation feedback -one can couple an ISM to an RCM.However, this is technically challenging and computationally expensive, effectively precluding exploration of climate and ice sheet modelling uncertainties.
The only way, therefore, to incorporate physical modelling of ice flow, SMB processes, and the SMB-elevation feedback while also exploring model uncertainties is with a parameterisation such as the one we present in a companion paper (Edwards et al., 2014), where we characterise the SMB response to elevation in MAR using a suite of simulations in which the MAR GrIS surface height is altered.The parameterisation is a set of four gradients, or "SMB lapse rates", that relate SMB changes to height changes below and above the equilibrium line altitude (ELA) and for regions north and south of 77 • N.Here we apply the parameterisation in five ISMs to adjust MAR projections of SMB under the SRES A1B scenario (Nakićenović et al., 2000) as the ice sheet geometry evolves.
The climate community have been attempting to quantify uncertainty in global climate model (GCM) predictions for some time, focusing on uncertainty in their parameter values with perturbed parameter ensembles (PPEs) but also attempting to estimate structural uncertainty, which is uncertainty about the remaining discrepancy between a model and reality at the model's most successful parameter values (e.g.Sexton et al., 2011;Sexton and Murphy, 2011).But probabilistic quantification of uncertainties in ice sheet model predictions has hardly yet been attempted.There is also an urgent need to propagate uncertainties along the causal chain from greenhouse gas forcing scenarios to the impacts of climate change.ISM projections are only beginning to tackle these challenges.Multi-model comparisons such as MIS-MIP (Pattyn et al., 2012) and SeaRISE (Bindschadler et al., 2013) compare ISMs with each other using standardised experiments, but the ice2sea project (http://www.ice2sea.eu) is among the first to explore systematically emissions scenarios and climate and ice sheet model uncertainties within a coherent framework (e.g.Rae et al., 2012;Fettweis et al., 2013).
In Edwards et al. (2014), we make uncertainty assessment an integral part of the parameterisation by estimating probability distributions for the four parameters (elevation feedback gradients).Here we propagate these uncertainties to future projections by sampling values from the four distributions.We also use five different ISMs to assess the effect of ISM structural uncertainty that arises from different representations of ice flow and initialisation procedures, and explore GCM structural uncertainty by forcing MAR with two different GCMs.Our coherent experimental design and statistical framework, unusual in sea level projections, allow us to propagate the three types of model uncertainty along the causal chain from SRES scenario to sea level contribution; assess the relative importance of these through time; and present probabilistic assessments of the effect of elevation feedback parameterisation uncertainty on the projected GrIS sea level contribution under the A1B scenario.

Climate projections
The regional climate model MAR (Fettweis, 2007) has been adapted for simulating the climate over Greenland, with full coupling to a complex snow/ice energy balance model and relatively high horizontal resolution (25 km).Unlike most RCMs, MAR includes the positive feedback between ice surface albedo and melting (Franco et al., 2013) et al. (2014).We use five climate simulations performed for the ice2sea project.The first is a twenty year simulation of 1989-2008 in which MAR is forced at the boundaries by the ERA-INTERIM reanalysis (Dee et al., 2011), which we use as a baseline for initialisation and projections.The next two are twenty year simulations of 1980-1999 in which MAR is forced by two GCMs, ECHAM5 (Roeckner et al., 2003) and HadCM3 (Gordon et al., 2000), under 20th century climate forcings, which we use to calculate projection anomalies.The final two are one hundred year projections  forced by the two GCMs under the SRES A1B emissions scenario (Nakićenović et al., 2000).The A1B scenario is only defined until 2100, but we wish to make GrIS projections for two hundred years .We did not have resources to simulate the second century (the climate responding to the A1B scenario fixed at 2100 values), so we extend the MAR simulations by repeating the final decade (2090-2099) ten times.This method is likely to underestimate the longer-term climate response, as discussed later.
MAR simulates SMB values only over grid cells categorised as permanent ice, so the low spatial resolution (relative to ISMs) leads to missing values for parts of the GrIS margin.As part of the data processing, we extrapolate the SMB field outwards from the ice sheet using a simpler representation of the SMB-elevation relationship than in Edwards et al. (2014).For each land grid cell, we use the ice grid cells within 100 km to fit a linear model of SMB versus surface height, and use this relationship with the land grid cell elevation to assign an SMB value.This was performed before the parameterisation of Edwards et al. (2014) was developed.

Ice sheet models
We implement the parameterisation in five ISMs with varying structure and complexity: Elmer/Ice, which solves the full Stokes equations; GISM, MPAS and CISM, which use a higher order approximation that reduces computational expense (HO: e.g.Pattyn, 2003); and GRISLI, which uses a hybrid of the first order shallow ice and shallow shelf approximations (SIA and SSA: Bueler and Brown, 2009) to further reduce computational expense.We also use a SIA version of GISM, and call the two versions GISM-HO and GISM-SIA.Elmer/Ice and MPAS use finite element numerical methods on unstructured grids, while the others use finite difference methods on regular grids at 5 km resolution.
Elmer/Ice builds on Elmer, the open-source finite element, partial differential equation solving parallel code mainly developed by the CSC-IT Center for Science Ltd in Finland.The unstructured mesh allows a variable grid resolution to focus computational resources at the ice sheet margin; here we use a minimum horizontal grid size of less than 1 km.GISM-SIA is a thermomechanical ISM (Huybrechts and de Wolde, 1999) that has been modified and extended for pro-jections on centennial timescales using a new higher-order approximation of the force balance (GISM-HO: Fürst et al., 2011Fürst et al., , 2013)).The MPAS-Land Ice model is based on the MPAS (Model for Prediction Across Scales) climate modelling framework of Ringler et al. (2008); here we use a regular 5 km resolution hexagonal mesh.GRISLI (GRenoble Ice Shelves and Land Ice model) is a thermomechanically coupled ISM; it is a hybrid model that for grounded ice uses the SIA for vertical shearing and the SSA as a sliding law (Ritz et al., 2001;Bueler and Brown, 2009).The Community Ice Sheet Model (CISM) version 2.0 includes improvements to all components of the Glimmer-CISM SIA model (Rutt et al., 2009).The models are summarised in Table 1; more detailed information is given elsewhere (Elmer/Ice: Gillet-Chaulet et al., 2012; GISM-HO: Goelzer et al., 2013;MPAS: Perego et al., 2012;CISM: Price et al., 2011;Lemieux et al., 2011;Evans et al., 2012;GRISLI: Quiquet et al., 2012).

Initialisation
Determining initial conditions for ISMs involves finding a balance between observations of the present-day ice sheet, reconstructions of past climate changes (to which the ice sheet is still responding), and the physical laws and parameterised processes incorporated in the ISM, while accounting for uncertainties and limitations in all of these.ISM initialisation mostly uses ad hoc tuning methods rather than formal data assimilation as in numerical weather forecasting.Our initialisation procedures use observations of present-day ice sheet geometry (surface elevation, bedrock elevation and ice thickness; e.g.Bamber et al., 2013), ice velocities (Joughin et al., 2010;Bamber et al., 2000), and geothermal heat fluxes (Shapiro and Ritzwoller, 2004).
Methods vary between ISMs but in general there are three stages.First, we estimate the 3-D ice-temperature field by solving the heat equation, in some cases accounting for the response to atmospheric temperature changes over one or more glacial-interglacial cycles.Second, we infer the spatial pattern of the basal drag coefficient that leads to the best agreement with observed ice velocities, given observational uncertainties and model limitations.During these first two stages, the geometry of the ice sheet is held constant to the observed values (Bamber et al., 2013).For the third stage, we finally allow the ice sheet geometry to evolve or "relax" so that it is internally consistent with the ice temperature and flow fields.
We derive ice temperatures for Elmer/Ice using the computationally cheaper SIA model SICOPOLIS (Greve, 1997).We derive them for CISM, MPAS and GRISLI from a quasisteady-state CISM simulation forced with a constant surface temperature field simulated by the RACMO/GR regional climate model (Ettema et al., 2009, as  of several glacial-interglacial cycles, rescaling the ice temperature field to the observed ice thickness (Goelzer et al., 2013).Geothermal heat flux fields in all spin-up simulations are from Shapiro and Ritzwoller (2004).Ice temperature spin-up simulations are summarised in Table 1.
We infer the basal drag coefficient from observed velocities using the control method (Elmer/Ice: Morlighem et al., 2010), an iterative inverse method (GRISLI), or tuning (CISM, MPAS); the methods and velocity targets are also summarised in Table 1.
We obtain the initial ice sheet geometry by forcing the ISM with a present-day mean SMB field.For most models this is the 1989-2008 mean (the reference period for the ice2sea project) of the ERA-INTERIM forced MAR simulation; for GISM we use the 1960-1990 mean SMB calculated with the GISM PDD scheme from ERA-40 data (Uppala et al., 2005) and observed surface elevation (similar to Hanna et al., 2011).We obtain the initial geometry for Elmer/Ice by allowing the upper surface to evolve for 55 yr under the present-day SMB forcing, and use the same initial geometry for GRISLI, allowing it to evolve for a further 145 yr to obtain a coherent initial state for this model.We initialise the GISM models by allowing the geometry to evolve for 1000 yr while the ice sheet margin is fixed to observed values under present-day SMB forcing and ice thickness changes are limited to 0.2 m a −1 .These relaxation periods are determined arbitrarily or by available computational resource, though most of the flux divergence anomalies that arise from observational and modelling uncertainties die away after a few years (e.g.Gillet-Chaulet et al., 2012).We do not relax the geometry for CISM and MPAS.However, these simulations are already closer to the present-day state than for the other models because they use a quasi-steadystate spin-up forced with modern-day surface temperatures, rather than glacial-interglacial cycles, and calibrate the basal drag with balance velocities (which have complete coverage and are in equilibrium with the geometry), rather than surface velocities.
We use two methods to correct any remaining drift that arises from model imbalances.For Elmer/Ice and GRISLI, we perform a control simulation and subtract this from the results (Sect.2.5).For the others we diagnose and apply a "synthetic SMB" correction, which is the additional SMB required to keep the ice sheet close to the present-day observed geometry under present-day SMB forcing.The synthetic SMB correction is applied unaltered in perturbed simulations (i.e. the parameterisation adjusts only the MAR SMB forcing and not the synthetic SMB correction).This is similar, in principle, to the flux corrections that were formerly in common use in atmosphere-ocean GCMs.In CISM and MPAS there is no relaxation in the initialisation procedure, so the synthetic SMB accounts for the entire initial mass imbalance generated by the initial geometry.The mean synthetic SMB correction for both GISM models is 9 cm a −1 , with less than 0.3 % of the total ice sheet area having a correction of greater than 10 m a −1 , mainly at marine-terminating ice margins (Goelzer et al., 2013).For MPAS the mean synthetic SMB correction is −7.3 cm a −1 , with 1.6 % of the total ice sheet area having a correction of greater than 10 m a −1 , and The Cryosphere, 8, 195-208, 2014 www.the-cryosphere.net/8/195/2014/1989-2008 1989-2008 1960-1990 1960-1990 1989-2008 1989-2008 for CISM it is 2.4 cm a −1 , with less than 2 % of the total ice sheet area having a correction of greater than 10 m a −1 .We exclude Greenland's glaciers and ice caps from the projection results using a mask produced for ice2sea (Rastner et al., 2012).
The basal drag calibration and initialisation methods are summarised in Table 2.More details can be found in the references given in the previous section and the table.

Boundary conditions
The boundary conditions are geothermal heat flux, ice temperature, bedrock elevation, and projections of SMB.Geothermal heat fluxes are from Shapiro and Ritzwoller (2004) and held fixed in the simulations.We also keep the ice temperature fields (determined during initialisation; Sect.2.3 and Table 1) held fixed because there is negligible ice temperature response to changing atmospheric forcing over two centuries.We use bedrock elevations from the ice2sea data set (Bamber et al., 2013) and hold these fixed because there is negligible isostatic adjustment on this timescale (e.g.Goelzer et al., 2013).
The SMB forcings are therefore the only time-varying boundary conditions.We force the ISMs with "anomalycorrected" SMB projections, S RCM , to remove the mean discrepancy in MAR between the GCM-forced and ERA-INTERIM reanalysis-forced simulations (Fettweis et al., 2013).For this we calculate SMB anomalies with respect to the present day by subtracting the 1989-2008 mean SMB (1989-1999: 20th   S ERA40  1960−1990 ).If used, the synthetic SMB field is also applied.Our initialisation and anomaly-correction methods use the approximation that SMB trends during the period 1989-2008 are small relative to the A1B projections (Rae et al., 2012).

Parameterisation
We use the SMB-elevation feedback parameterisation to adjust the anomaly-corrected SMB forcing as the GrIS surface height in each ISM evolves.The adjustment is made each year, for each grid cell, using one of four parameters (gradients) selected according to the current "reference" SMB and the region of the grid cell (Edwards et al., 2014).
For a given ISM grid cell in a given year (t), a gradient (b t ) (kg m −3 a −1 ) is used to adjust the SMB forcing (kg m −2 a −1 ) using the height difference (m) between the previous year and the start of the simulation: where In Edwards et al. (2014), we estimate probability distributions for each of the four parameters (gradients).Our best estimate values of the parameters are the modes of these four distributions, and our credibility intervals (CIs) encompass the central 95 % of each.We propagate the parameter uncertainties to sea level projections by performing simulations using different values sampled from the distributions: (a) the four modes, i.e. best estimates; (b) the four 2.5 % quantiles, i.e. lower bounds of the CIs; and (c) the four 97.5 % quantiles, i.e. upper bounds.The parameter estimates are given in Table 3.They differ slightly to those in Edwards et al. (2014) due to changes made to the parameter estimation code for an update to that study: a small error was corrected, unnecessary rounding was removed, and the code was restructured (thus changing the random seed for bootstrapping).One difference is 0.03 (SMB ≥ 0, 97.5 % bound), and the other differences are 0.02 or less.
With one model, GISM-SIA, we also sample the four distributions more thoroughly by using 99 percentiles: i.e. the four 1 % quantiles in one simulation; the four 2 % quantiles in the next, and so on up to the four 99 % quantiles.
We therefore generate five simulations for each ISM: -Control: forced with S init (Elmer/Ice, GRISLI) or S init + S syn (other models), to check for or subtract model drift from the projections; -No feedback: forced with S RCM with no adjustment, to estimate the response of the GrIS without elevation feedback; -Best estimate: forced with S adj , using the "best estimate" values of the four parameters (Table 3); -2.5 % and 97.5 % quantiles: as for best estimate, but using the parameter values that correspond to the bounds of the 95 % CI (Table 3).
We perform these for the two MAR simulations forced by ECHAM5 (all ISMs) and HadCM3 (Elmer/Ice, GISM-HO, GISM-SIA and GRISLI) under the A1B scenario over 2000-2200.As described above we also use GISM-SIA to perform a simulation with each of the 99 percentile estimates of the four parameters, forming a 99 member PPE; we do this only for the ECHAM5 projection.
The drift in the control simulations is very small (0.03 %, or less, of the cumulative projected sea level contribution at 2200) for all models except Elmer/Ice, which has a drift of 2-2.5 % (−4 mm).For this model, drift is not constrained during the initialisation procedure and the free-surface elevation has been allowed to diverge from observations for a relaxation period of only 55 yr.The applied SMB is not corrected by a synthetic SMB, and the remaining drift shows the drift of the model when directly applying the 1989-2008 mean SMB given by MAR forced under ERA-INTERIM.This drift is corrected by subtracting the control simulation from the projections.

Results
The additional cumulative sea level contribution due to the SMB-elevation feedback, under the A1B scenario and averaged over five ISM projections for ECHAM5 and three for HadCM3, is 4.3 % (best estimate; 95 % credibility interval 1.8-6.9%) at 2100, and 9.6 % (best estimate; 95 % credibility interval 3.6-16.0%) at 2200 (Figs. 1 and 2; Tables 4 and  5).We exclude GISM-SIA from summary statistics because it is so similar in structure, and identical in set-up, to GISM-HO that it would effectively give double weighting to GISM.In all ISM predictions the elevation feedback is positive: the lower bounds of the 95 % CIs (the 2.5 % quantile simulations) are always larger than the "no feedback" case.This amplifying effect is expected, because mean GrIS SMB becomes negative at around 2050 (Rae et al., 2012) and the gradient estimates for negative SMB are generally positive (Edwards et al., 2014) thus adjusting the SMB in the same direction: more negative, giving greater sea level contribution.However, this was not guaranteed, because there are both positive and negative values in the mean SMB forcing, in the SMB changes from 2000-2199, and in the gradient sets (Edwards et al., 2014).
The larger feedback adjustment at 2200 arises from the negative SMB forcing at the end of the first century which is sustained (by repetition of the decade 2090-2099 from the MAR simulation: Sect.2.1) during the second century.
The Elmer/Ice and GRISLI total sea level contributions are virtually identical, but this is a coincidence: the changes in ice discharge and SMB are different, and the differences just compensate each other so that the sea level contributions are the same.The two models do have similar initialisation procedures, but the differing dynamic and SMB responses lead us to believe this is not the reason for the similar results.CISM seems to be an outlier in the ISM ensemble, particularly at 2200 (Fig. 2; Table 5), in both the no feedback result (low) and the magnitude of the feedback (small), though the latter may be a result of the former.The GISM-SIA results are very similar to those of GISM-HO.This can be explained by the identical model set-up and their similarity in response to perturbations; the only difference between the two models is in the dynamic response, which is a much smaller contribution than the SMB response (Fürst et al., 2013).We can compare ECHAM5 and HadCM3 using the three ISMs that used both projections (Elmer/Ice, GISM-HO and GRISLI, from here on called the "starred" models).The HadCM3-forced projection gives consistently higher sea level contributions than the ECHAM5-forced, because the HadCM3-forced MAR simulation projects greater meltwater runoff (Rae et al., 2012).Without elevation feedback, the mean GrIS sea level contributions at 2100 are 57 mm for ECHAM5 and 64 mm for HadCM3; with feedback, these increase by 3 mm (5 %).At 2200, the mean contributions without feedback are 164 mm for ECHAM5 and 176 mm for HadCM3; the feedback contributes an additional 16 mm (10 %) and 19 mm (11 %), respectively.
Not only does the contribution due to the feedback increase with time, but so does the uncertainty.Our experimental design allows us to inspect the changes and relative importance of the different modelling uncertainties.At 2100 (Fig. 1; Table 4), the GCM differences are, respectively, 2.6 and 1.7 times those of the elevation feedback parameterisation and ISM uncertainties: the mean difference between the two GCMs (averaged over the no feedback and three feedback estimates) for the three starred models is 7.8 mm; the mean spread of all five ISMs for ECHAM5 is 4.5 mm; and the mean 95 % CI (over all ISMs and GCMs) is 3.0 mm.But by 2200 (Fig. 2; Table 5), the relative importance has changed.The ISM spread and elevation feedback parameterisation uncertainty overtake the difference between the two GCMs: the mean difference between the two GCMs is 13.7 mm; the mean spread of ISMs for ECHAM5 is 25.8 mm; and the mean 95 % CI (over all ISMs and GCMs) is 20.8 mm.In other words, the parameterisation and ISM uncertainties increase by a factor of seven and six, respectively, between 2100 and 2200, while the difference between the two GCMs barely doubles.The relative contributions to our uncertainty thus depend on the timescale of interest.This can be seen clearly in Fig. 3, which shows the fractional uncertainties as a function of time with a 15 yr running mean applied.Although we have a small number of ice sheet and climate models, it is useful to examine the relative importance of different uncertainties in our study to inform again (the range increases faster than the signal).The time dependence of the GCM structural uncertainty explored (the difference between ECHAM5 and HadCM3 results divided by their mean, averaged over the best estimates for the three starred ISMs) appears at first glance to be similar, but the underlying reason is different: it rapidly falls to a minimum when the two GCMs projections happen to coincide for a short period (see also Rae et al., 2012).After a short increase due to the GCMs diverging again, the fractional uncertainty slowly falls because the signal increases faster than the range.The fractional uncertainty due to the parameterisation (95 % CI divided by the best estimate, averaged over all ISMs and GCMs) is initially very small, but increases fairly linearly with time: the width of the CI scales with the signal.In the latter part of the second century the parameterisation uncertainty overtakes that of the GCMs.These results are discussed further in the next section.
The 99 member GISM-SIA PPE gives the entire sea level probability distribution (Fig. 4) because it jointly samples the full probability distributions of each of the four gradient parameters rather than just the bounds of the 95 % CIs.We can inspect any other interval, such as the 90 % CI .We can also test whether the projected distribution of sea level contributions is symmetric about the best estimate of 186.1 mm.The 95 % CI is symmetric (±10.4 mm), and the 90 % CI nearly so (−8.8 mm, +8.6 mm), but the 50 % interval (−3.8 mm, +3.5 mm) indicates the probability density leans slightly towards the higher sea level contributions.Gradient (kg m −3 a −1 ) Sea level contribution at 2200 (mm) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q North, SMB < 0 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q South, SMB < 0 North, SMB >= 0 South, SMB >= 0 The unimodal and (broadly) symmetric shape of the PPE sea level distribution is probably a result of the unimodal and (in most cases) symmetric shapes of the four posterior probability distributions of the gradient values (Edwards et al., 2014), combined with the linear model of adjustment.Figure 5 shows the dependence of the sea level contribution on the gradient sets.However, this cannot show which of the four gradients is most important because they are varied jointly.But it illustrates some characteristics of the parameterisation structure and uncertainties: first, the relationship between sea level contribution and SMB lapse rates is fairly linear, save for a clear change in slope for the south SMB ≥ 0 lapse rate.This is due to the linear adjustment model, but is perhaps surprisingly robust considering the SMB ≥ 0 lapse rates range from negative to positive values, and considering there is a positive feedback in adjusting SMB.In fact the latter effect can be seen in the SMB < 0 slopes, which become slightly steeper at larger SMB lapse rates.Second, the magnitude of the slopes illustrate the characteristics of the SMB lapse rate probability distributions (Edwards et al., 2014): shallower slopes for the two SMB < 0 lapse rates are a result of their broader probability distributions, i.e. larger absolute uncertainties, relative to the SMB ≥ 0 lapse rates.Similarly, the south SMB ≥ 0 slope change is due to the asymmetric probability distribution.

Elevation feedback
It is important to quantify the SMB-elevation feedback, relative to projections without it (such as Rae et al., 2012), to improve our assessments of GrIS sensitivity to climate change.Using RCMs to simulate SMB represents the processes more fully than with ISMs; using ISMs to evaluate the feedback itself gives a more complete picture than with RCMs alone (by modifying their topography), because they simulate the entire GrIS response rather than only the SMB aspect.Parameterising the SMB-elevation feedback gives the best of both worlds, allowing ISMs to be forced with SMB from RCMs and include the feedback, while being computationally cheap enough to explore model uncertainties.The response of SMB to changes in ice surface elevation is complex and nonlinear, and even with a probabilistic uncertainty assessment, a linear parameterisation cannot fully represent the SMB changes that would be modelled by a fully coupled RCM-ISM.However, the advantages of our particular parameterisation over other studies are described in the companion paper (Edwards et al., 2014).
We present the first application of an SMB-elevation parameterisation to RCM projections of future sea level rise.The feedback always increases sea level rise, even for the lower bounds of our 95 % CIs (despite the negative gradient values: Table 3).But the additional contribution is relatively small, with best estimates 4.3 % at 2100 and 9.6 % at 2200.
We can compare these results with estimates from the RCM alone.When MAR is forced by the GCM MIROC5 under the Representative Concentration Pathway RCP8.5 (Moss et al., 2010), changing the ice sheet height by an amount equivalent to cumulative projected SMB changes from 2000-2079 gives a much greater feedback contribution than our result: 5-15 % over just one century (Fettweis et al., 2013).At first glance this might appear to be due to the use of a high-end scenario (relative to the medium scenario A1B), giving larger SMB changes and therefore a larger feedback.
However, the greater feedback contribution in the RCM most likely arises from the absence of ice dynamics.Perhaps counter-intuitively, sea level rise from the total GrIS response (SMB and dynamics) is smaller than from the SMB response alone.Recent ice2sea studies (e.g.Gillet-Chaulet et al., 2012;Goelzer et al., 2013) have highlighted the complexity of the interactions between SMB and dynamic responses.Dynamic mass loss and SMB-related mass loss are not additive: mass loss by one process decreases loss by the other, so we cannot neatly divide sea level contributions into, for example, 60 % SMB and 40 % dynamic.SMB removes ice (e.g. through melting) before it can reach the ice sheet margin to be discharged.Conversely, dynamics remove ice before it can melt and also partly compensate surface height decreases (from increased melting) by redistributing mass gained at the top of the ice sheet, thickening the ablation zone and therefore www.the-cryosphere.net/8/195/2014/The Cryosphere, 8, 195-208, 2014 reducing melting.Further discussion can be found in Gillet-Chaulet et al. (2012) and Goelzer et al. (2013), though this is not a recent finding: Huybrechts and de Wolde (1999), for example, find a 24 % lower sea level contribution by 2100 from Greenland when dynamics are included.In our study the time integrated SMB changes for the ECHAM5forced MAR simulation contribute 62.7 mm to sea level by 2100 and 211.4 mm by 2200 (Goelzer et al., 2013), while the total responses from the ice sheet models are substantially smaller: 53.7-57.4mm at 2100 and 150.8-170.5 mm at 2200 (Tables 4 and 5).Therefore, it is not surprising that SMB-only projections, without coupling to an ice sheet model, show larger SMBelevation feedback contributions.The results from Fettweis et al. (2013) are likely overestimates.Only a coupled RCM-ISM would allow full evaluation of this response, but our results suggest that the feedback effect is small enough that using such models for one-or two-century simulations may not be justified given the additional computational expense, technical difficulty, and severe restrictions on exploring modelling uncertainties.
We can also compare these results with estimates from an ISM alone.Some have previously suggested that the simpler PDD descriptions of ice sheet response are too sensitive to climate change (van de Wal, 1996;van de Berg, 2011), but recent studies find they are less sensitive than MAR and RACMO/GR (Goelzer et al., 2013;Helsen et al., 2011;Vernon et al., 2013;Hanna et al., 2011).Goelzer et al. (2013) compare projections from our parameterisation versus the GISM PDD scheme, using the same ice sheet extent and RCM, and Helsen et al. (2011) compare their parameterisation of the RACMO2/GR SMB-elevation feedback with PDD results.In both cases the PDD schemes predict a smaller sea level contribution in a warmer climate than the RCMs: for example, Goelzer et al. (2013) find the PDD model gives a 31 % lower sea level contribution than MAR.There are many differences between physically based energy balance RCM schemes and empirical PDD schemes, including the presence of albedo feedback in the RCMs, higher spatial resolution in the PDD schemes, and different treatment of refreezing (Goelzer et al., 2013).Our parameterisation of the SMB-elevation feedback is more complex than a PDD scheme, in that it estimates probability distributions rather than a fixed lapse rate, includes nonlinear precipitation aspects, and incorporates variation of the feedback with climate, topography and region (through the use of four probability distributions).
This parameterisation could be used to correct low resolution MAR simulations onto a high resolution digital elevation map; the difference in elevation between the two topographies is analogous to the geometry evolution in ISM simulations.This correction might improve the realism of MAR SMB simulations, which would be useful not only for forcing ISMs but also comparing with observations.

Modelling uncertainties
It is a research priority to improve quantification of ice sheet model uncertainties, so as to make meaningful and complete statements about the difference between our projections and the real world.Research in this area is sparse, but the ice sheet modelling community can make rapid improvements by applying established uncertainty quantification techniques from climate and other environmental modelling research areas, as we have done here.
Our method is novel because we propagate three types of model uncertainty -GCM structural uncertainty, ISM structural uncertainty, and elevation feedback parameterisation uncertainty -along the causal chain from SRES scenario to sea level within a coherent experimental design and statistical framework.In particular, parameter uncertainty is estimated with a probabilistic method, which gives well-defined credibility intervals (CIs) rather than simple sensitivity analysis; all ISMs use the same parameterisation and parameter sampling, and are forced with the same GCMs and RCM, enabling comparisons across ISMs; and MAR is forced with two GCM forcings under the same SRES scenario, enabling comparisons between ECHAM5 and HadCM3.A similar approach is taken by Shannon et al. (2013) for a new parameterisation of GrIS basal lubrication.This experimental design gives us the opportunity to assess the relative importance of various modelling uncertainties and how these vary with time.We discuss parametric uncertainties, structural uncertainties, time dependence, and future directions.
The feedback parameter uncertainties translate to a unimodal and fairly symmetric probability distribution for sea level contribution (Fig. 4), because of the linear SMB adjustment and the sampling from unimodal and (in most cases) symmetric probability distributions for the gradients (Edwards et al., 2014).If the gradient distributions derived from MAR were strongly asymmetric or multimodal, or the SMB adjustment nonlinear, the sea level distribution might be a different shape.Here is the real value of Bayesian parameter estimation (Edwards et al., 2014) and sampling these uncertainties with perturbed parameter ensembles, because such things are generally not known in advance.The slight asymmetry in the sea level probability distribution described in the previous section might become more pronounced for timescales longer than two centuries.
If substantially greater computational and time resources were available, a natural next step might be a more thorough exploration of parametric uncertainties, by simultaneously perturbing additional parameters of the ISMs and, if possible, parameters of the climate models.Climate model PPEs are used extensively (e.g.Murphy et al., 2009); RCM PPEs are relatively rare due to their computational expense (e.g.Sexton and Murphy, 2011).Two ISM PPEs are presented by Stone et al. (2010) and Applegate et al. (2011), who vary multiple parameters simultaneously in order to explore their interactions.For example, surface lowering driven by Structural uncertainties are much more challenging to quantify.Using multi-model ensembles is a necessary first step.But there is no well-defined space from which to sample (as there is for parameters).Models have common components and errors that are likely to reduce the ensemble range relative to our uncertainty about the climate and GrIS.It is extremely difficult to express this probabilistically -to quantify the degree to which the spread of a model ensemble encompasses our true uncertainty about the system -and there is no consensus about how to do so.Given these difficulties, and our small ensemble sizes, we do not claim to quantify ISM and GCM structural uncertainties, only to sample from them.It is better to sample from this space, to explore it partially, than to ignore it.Coherent experimental design such as this -systematic sampling of models, parameterisation, and parameter values -is essential to better understanding of uncertainties.The structural uncertainty we explore most is in the ISMs.This arises from structural differences in ice sheet modelling, including representation of the Stokes equations, numerical methods, grid resolution, initialisation procedures, and drift correction.Given these differences, it might seem surprising that the results from different ice sheet models are so similar.Goelzer et al. (2013) compare GrIS sea level projections with different plausible modelling choices for a single ice sheet model (GISM-HO) and the same RCM (MAR) and GCM (ECHAM5) as here, and obtain a wider range of projections.But our results confirm the importance of SMB versus dynamic mass loss.The large spread in results in Goelzer et al. (2013) arises from differing SMB projections: an RCM, a PDD model forced with an RCM, and a PDD model forced with a GCM.In our study all models (by definition) use the same projected SMB changes, and also use the same mask to define ice sheet extent (Sect.2.3).The dynamical response of the GrIS is smaller than the SMB response, so differences in model physics, numerical methods, grid and even initialisation cannot substantially alter the sea level contribution: SMB is the dominant driver of sea level rise.This confirms results by others such as Gillet-Chaulet et al. (2012) and Goelzer et al. (2013).
We use only two GCMs and do not explore RCM structural uncertainty at all.We do incorporate uncertainty in the structure of our parameterisation of MAR, by including a discrepancy variance term when estimating the gradients (Edwards et al., 2014), but this refers to our parameterisation of the MAR elevation feedback, not the MAR simulation of the feedback in the real world; the latter is the RCM structural uncertainty.We could sample RCM structural uncertainty using multiple RCMs, but this is extremely challenging due to their computational expense.An ensemble of RCMs would also require a new parameterisation of the SMB-elevation feedback for each.For further exploration of MAR and GCM uncertainties, we refer to Fettweis et al. (2013) who evaluate the effects of changing the MAR tundra/ice mask and ice sheet topography, and using a range of GCMs for forcing, on simulations of current climate and projections of future climate change.
We therefore do not claim to make a full assessment of modelling uncertainty in sea level projections.Such an assessment would require perturbing other parameters of the ice sheet models, perturbing parameters of the climate models, and quantifying ice sheet model and climate model structural uncertainties through, for example, large multimodel ensembles of ISMs, RCMs and GCMs.This is beyond present capabilities, though we believe it should be our aim.These additional uncertainties have potentially large effects on predictions of future mass loss, which would result in wider credibility intervals than those presented here.We reiterate that our uncertainty assessment pertains mainly to uncertainties from the elevation feedback parameterisation, with additional exploration of the effects of ice sheet and global climate model uncertainties.
Our quantification of the relative magnitude of modelling uncertainties through time (Fig. 3) is specific to our particular choices and opportunities: the RCM, GCMs and ISMs; the emissions scenario; our method of extending the scenario from 2100 to 2199; and (trivially) our chosen CI width.If we used a scenario in which the SMB projections were more negative than for A1B, such as RCP8.5, or less negative, such as E1 (Rae et al., 2012), the elevation feedback would be correspondingly larger or smaller and so would the uncertainty.Similarly, if we were to use 200 yr GCM and RCM simulations, that is, drive MAR with GCM simulations that were forced in the second century with the A1B scenario at 2100, rather than repeating the decade 2090-2099 from the MAR simulations ten times, the SMB forcing would likely continue to decrease rather than remain fixed.This would lead to a more rapid increase in the sea level projections with elevation feedback.The ISM fractional uncertainty might then not increase in the second century, or not by as much (Fig. 3); this would depend on the relative responses of the ISMs.The GCM range might increase more quickly than the mean, giving a time-dependent fractional uncertainty that is more similar in shape to the ISM fractional uncertainty (Fig. 3).
Such choices and limitations are unavoidable given the substantial computational expense of a three stage model chain in which two stages are climate models.Nevertheless, this work begins the process of determining where to focus computational resources and model development for different decision-making timescales.We present the first application of a parameterisation of the GrIS SMB-elevation feedback to projections of future climate change.This is currently the only practicable method for making projections for the GrIS contribution to sea level that can incorporate physical modelling of ice flow, SMB processes, and the SMB-elevation feedback while also exploring modelling uncertainties.
We use a coherent experimental design and Bayesian framework to make probabilistic statements (credibility intervals) about the effect of feedback parameterisation uncertainty on sea level projections and explore GCM and ISM structural uncertainties.Probabilistic assessments are more meaningful than single projections or sensitivity analyses: they are easier to interpret, more straightforward to compare across studies, and provide more robust and complete information for decision-making.
We make projections for the SRES A1B emissions scenario with the MAR RCM, ECHAM5 and HadCM3 GCMs, and five ISMs.The experimental design allows us to estimate the relative importance of these model uncertainties, and how they vary over time, for the GrIS contribution to sea level over the next two centuries.The additional sea level contribution due to the SMB-elevation feedback, averaged over five ISM projections for ECHAM5 and three for HadCM3 is 4.4 % (best estimate; 95 % credibility interval 1.8-6.9%) at 2100, and 9.6 % (best estimate; 95 % credibility interval 3.6-16.0%) at 2200.In all results, the lower bounds of our 95 % credibility intervals for sea level contributions with the SMB-elevation feedback are larger than for the "no feedback" case.
The relative contributions to uncertainty in our study vary with time.In 2100, the GCM differences are larger than the ISM range and the parameterisation uncertainty, but by 2200, the ISM and parameterisation uncertainties are greater.This work indicates that the areas where future research should be directed -such as understanding climate or ice sheet model differences, or performing perturbed parameter ensemblesdepend on the timescale of interest.For sufficiently large ensemble sizes, this method can be used to give an indication where future sea level research should be directed.
The relatively small sea level rise contributed by the elevation feedback suggests that century-long simulations with coupled RCM-ISM models may not be justified for the additional computational expense and technical difficulty of coupling, at least for medium to low emissions scenarios.

Fig. 5 .
Fig. 5. GISM-SIA projections for sea level contribution in 2200 as a function of elevation feedback gradient values.

www.the-cryosphere.net/8/195/2014/ The Cryosphere, 8, 195-208, 2014Table 1 .
Shapiro and Ritzwoller (2004)op) treatment of Stokes equations and (bottom) derivation of ice temperature (IT) boundary conditions.HO, SIA and SSA are higher order, shallow ice and shallow shelf approximations; SS is steady state; SAT is surface air temperature data used for the spin-up simulation; g-ig is glacial-interglacial.All models use bedrock elevation data fromBamber et al. (2013)and geothermal heat flux fromShapiro and Ritzwoller (2004).

Table 2 .
Summary of ice sheet model basal drag calibration and initialisation.PDD is positive degree day.
century climate forcings; 2000-2008: A1B scenario) from the A1B projections.We add these anomalies to the present-day SMB used for initialisation, to give S RCM 2000,...,2199 = S A1B 2000,...,2199 − S 20C,A1B 1989−2008 + S init , where S init is generally S ERAI 1989−2008 , the 1989-2008 mean of the ERA-INTERIM forced simulation (except GISM, is selected from one of four values according to the reference SMB (S ref < 0 or S ref ≥ 0) and region (north or south of 77 • N) of the grid cell, where S ref t is the mean S adj of the previous decade or, for the first decade, all available years.

Table 3 .
The 2.5 % quantile, best estimate, and 97.5 % quantile estimates of the SMB-elevation gradients in kg m −3 a −1 , below (SMB < 0) and above (SMB ≥ 0) the ELA, for regions north and south of 77 • N.

Table 5 .
As for Table 4, but at 2200.

The Cryosphere, 8, 195-208, 2014 www
Shannon et al., 2013) processes (such as basal lubrication:Shannon et al., 2013)could interact differently with the elevation feedback depending on the values of the relevant parameters.Performing multiple parameter perturbations might affect the shape of the distribution of projected sea level contributions. .the-cryosphere.net/8/195/2014/