# Probabilistic parameterisation of the surface mass balance-elevation feedback in regional climate model simulations of the Greenland ice sheet

T. Edwards, X. Fettweis, O. Gagliardini, F. Gillet-Chaulet, H. Goelzer, J.

Gregory, M. Hoffman, P. Huybrechts, A. Payne, M. Perego, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

T. Edwards, X. Fettweis, O. Gagliardini, F. Gillet-Chaulet, H. Goelzer, et al. Probabilistic parame-
terisation of the surface mass balance-elevation feedback in regional climate model simulations of the
Greenland ice sheet. The Cryosphere, 2014, 8 (1), pp.181-194. 10.5194/tc-8-181-2014. hal-03778266

HAL Id: hal-03778266

https://hal.science/hal-03778266

Submitted on 15 Sep 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

# Probabilistic parameterisation of the surface mass balance–elevation feedback in regional climate model simulations of the Greenland ice sheet

T. L. Edwards¹, X. Fettweis², O. Gagliardini³,⁴, F. Gillet-Chaulet³, H. Goelzer⁵, J. M. Gregory⁶,⁷, M. Hoffman⁸, P. Huybrechts⁵, A. J. Payne¹, M. Perego⁹, S. Price⁸, A. Quiquet³, and C. Ritz³

¹Department of Geographical Sciences, University of Bristol, Bristol BS8 1SS, UK

²Department of Geography, University of Liege, Laboratory of Climatology (Bat. B11), Allée du 6 Août, 2, 4000 Liège, Belgium

³Laboratoire de Glaciologie et Géophysique de l'Environnement, UJF – Grenoble 1/CNRS, 54, rue Molière BP 96, 38402 Saint-Martin-d'Hères Cedex, France

⁴Institut Universitaire de France, Paris, France

⁵Earth System Sciences & Departement Geografie, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussels, Belgium

⁶NCAS-Climate, Department of Meteorology, University of Reading, Reading, UK

⁷Met Office Hadley Centre, Exeter, UK

⁸Fluid Dynamics and Solid Mechanics Group, Los Alamos National Laboratory, T3 MS B216, Los Alamos, NM 87545, USA

⁹Department of Scientific Computing, Florida State University, 400 Dirac Science Library, Tallahassee, FL 32306, USA

*Correspondence to:* T. L. Edwards (tamsin.edwards@bristol.ac.uk)

Received: 30 December 2012 – Published in The Cryosphere Discuss.: 27 February 2013

Revised: 3 November 2013 – Accepted: 2 December 2013 – Published: 30 January 2014

**Abstract.** We present a new parameterisation that relates surface mass balance (SMB: the sum of surface accumulation and surface ablation) to changes in surface elevation of the Greenland ice sheet (GrIS) for the MAR (Modèle Atmosphérique Régional: Fettweis, 2007) regional climate model. The motivation is to dynamically adjust SMB as the GrIS evolves, allowing us to force ice sheet models with SMB simulated by MAR while incorporating the SMB–elevation feedback, without the substantial technical challenges of coupling ice sheet and climate models. This also allows us to assess the effect of elevation feedback uncertainty on the GrIS contribution to sea level, using multiple global climate and ice sheet models, without the need for additional, expensive MAR simulations.

We estimate this relationship separately below and above the equilibrium line altitude (ELA, separating negative and positive SMB) and for regions north and south of 77° N, from a set of MAR simulations in which we alter the ice sheet surface elevation. These give four “SMB lapse rates”, gradients that relate SMB changes to elevation changes. We assess un-

certainties within a Bayesian framework, estimating probability distributions for each gradient from which we present best estimates and credibility intervals (CI) that bound 95 % of the probability. Below the ELA our gradient estimates are mostly positive, because SMB usually increases with elevation: 0.56 (95 % CI: –0.22 to 1.33) kg m⁻³ a⁻¹ for the north, and 1.91 (1.03 to 2.61) kg m⁻³ a⁻¹ for the south. Above the ELA, the gradients are much smaller in magnitude: 0.09 (–0.03 to 0.23) kg m⁻³ a⁻¹ in the north, and 0.07 (–0.07 to 0.59) kg m⁻³ a⁻¹ in the south, because SMB can either increase or decrease in response to increased elevation.

Our statistically founded approach allows us to make probabilistic assessments for the effect of elevation feedback uncertainty on sea level projections (Edwards et al., 2014).
# 1 Introduction

Over the past two decades the Greenland ice sheet (GrIS) has been losing mass at an increasing rate, on average $142 \pm 49 \text{ Gt a}^{-1}$ with a total contribution to global sea level of about 8 mm (Shepherd et al., 2012). The GrIS has the potential to raise global sea level by several centimetres this century, and more in the next, with larger regional changes. The sensitivity of the GrIS to climate change is not well-known (IPCC, 2007), so it is important to improve estimates of its response and make projections of the resulting contribution to sea level over the next one to two centuries to inform policy and planning. Underestimating sea level rise would leave coastal cities around the globe at risk, while overestimating it could result in unwarranted expenditure on coastal defence. Projections should therefore include probabilistic assessments of uncertainty if they are to provide the most robust and complete information for making decisions.

Predictions of the GrIS response to projections of future climate change are made with physically based ice sheet models (ISMs) forced with climate model simulations. ISMs simulate both parts of ice sheet response: the flow of ice subject to its boundary conditions (dynamic); and surface mass balance (SMB), which is the sum of surface accumulation and surface ablation (broadly speaking, the balance of snowfall versus meltwater runoff). However, SMB models included in ISMs are usually rather simple. Most often they use an empirically derived positive degree-day (PDD) scheme, in which melting is parameterised as a function of the sum of daily air temperatures above melting point, and runoff is usually modelled from temperature and precipitation with a simple snow pack model (e.g. Janssens and Huybrechts, 2000). Daily climate means are often approximated from seasonal means to reduce the input data set size.

At the other end of the spectrum of model complexity are regional climate models (RCMs). These simulate the atmosphere and surface over a limited spatial domain, with higher spatial and temporal resolution than global climate models (GCMs), and are forced at their boundaries with GCM simulations or reanalysis data such as ERA-40. Some RCMs, such as MAR (Modèle Atmosphérique Régional: Fettweis, 2007) and RACMO2/GR (e.g. Ettema et al., 2009), include complex snow/ice schemes that represent many of the physical processes that govern SMB. Such RCMs have been shown to be quite successful in reproducing the current SMB of the GrIS (e.g. Ettema et al., 2009; Fettweis et al., 2011; Vernon et al., 2013). RCMs are computationally expensive so only short and/or a small number of simulations can be performed.

Some have suggested that PDD descriptions of ice sheet response are too sensitive to climate change (van de Wal, 1996; van de Berg, 2011). In contrast, comparisons made between RACMO2/GR and the Janssens and Huybrechts (2000) PDD model by Vernon et al. (2013) and Hanna et al. (2011) find the RCM is more sensitive. In an attempt to make the most robust comparison (e.g. using the same ice

sheet extent and forcing from the same RCM), Goelzer et al. (2013) find that a PDD model underestimates sea level rise by 14–31% compared to MAR. These large variations in response relative to RCMs may reflect the simplicity of the PDD scheme.

Ideally, then, we would prefer future projections of GrIS SMB to be made with the more complete representations in RCMs rather than simple parameterisations such as the PDD model (for example Rae et al., 2012; Fettweis et al., 2013). But the ice flow component of an ISM is still needed to simulate the dynamical response of the GrIS. ISMs are run at higher resolution than RCMs (kilometres rather than tens of kilometres), to better represent glacier flow at the ice sheet margin.

As the ice sheet evolves in response to climate change, it also affects the local climate through feedback processes. Some, like the ice albedo feedback, may be simulated within the RCM. Others relating to the dynamical response, including the evolving geometry of the ice sheet, can only be simulated by coupling the RCM and ISM, or else parameterising the feedback to adjust the input climate forcing throughout the simulation.

One important ice–climate feedback is the set of interactions between the atmosphere and the ice sheet surface elevation; here we focus on the feedback between the atmosphere and ice surface/snow pack. The two main parts of this SMB–elevation feedback are (i) temperature, where an initial increase in air temperature that leads to ice melting lowers the surface elevation and exposes the ice to warmer temperatures through the atmospheric lapse rate; and (ii) precipitation, where surface elevation changes affect air temperature and atmospheric circulation and therefore the location and amount of precipitation. Surface topography in RCMs is usually held constant, so they do not incorporate the elevation feedback at all. PDD schemes include a parameterisation of the temperature aspect of the feedback, using an atmospheric lapse rate to adjust the input temperature forcing as the ice sheet surface evolves. They do not represent the precipitation aspect of the feedback except, in some cases, through a scaling factor for temperature. Most PDD schemes assume constant feedbacks (temperature–elevation, i.e. atmospheric lapse rate correction; precipitation–elevation, i.e. scaling correction; and ice albedo) that do not vary across the ice sheet or with climate change (discussed by Robinson et al., 2010; Helsen et al., 2011; Stone et al., 2010), though there are exceptions (Tarasov and Peltier, 2002).

If we are to simulate SMB with an RCM and how that SMB is affected by ice topography changes (unlike Rae et al., 2012; Fettweis et al., 2013, who use RCMs with constant ice sheet topography), we must either couple an ISM to an RCM, or else force an ISM with RCM output using a parameterisation of the relationship in terms of an “SMB lapse rate”. Coupling an ISM to an RCM or GCM is rarely done because it is technically challenging (one example is Ridley et al., 2005), and because the climate models, particularly RCMs,
are too computationally expensive to simulate the timescales
of long-term ice sheet response. The computational expense
also drastically limits opportunities to perform multiple sim-
ulations to sample uncertainties in modelling choices.

The pragmatic solution is therefore to parameterise the
SMB–elevation feedback. This allows us to explicitly sim-
ulate the SMB and dynamical responses without the tech-
nical challenges and substantial computational expense of
coupling ISMs to RCMs. Provided the parameterisation ad-
equately represents the feedback in MAR, this allows us
to perform many simulations that we otherwise could not,
because we can force ISMs with MAR that have not yet
been coupled to it, and sample uncertainties in the feedback
and ice sheet modelling with additional simulations that we
would not otherwise have computational resources to per-
form.

Helsen et al. (2011) provide the first such parameter-
isation, for the relationship between SMB and height in
RACMO2/GR, and use this to adjust the SMB forcing ap-
plied to an ISM. Franco et al. (2012) derive relationships
between the individual components of SMB (snowfall, rain-
fall, meltwater runoff, and loss by sublimation and evapo-
ration) and elevation changes in MAR, to correct low res-
olution SMB simulations onto a higher resolution ice sheet
topography. Hakuba et al. (2012) study the SMB response to
surface elevation changes in a version of the ECHAM5 GCM
(Roeckner et al., 2003) by lowering the ice sheet topography
to 75%, 50% and 25% of the present day, though do not
parameterise the relationship. We develop on these studies
in method (presented here) and application (Edwards et al.,
2014).

We derive a new parameterisation for the elevation feed-
back in MAR using a suite of simulations in which the MAR
GrIS surface height is altered. The parameterisation is a set
of four gradients that relate SMB changes to height changes.
These can be used to adjust the input SMB forcing as the ice
sheet geometry evolves (Edwards et al., 2014). The four gra-
dients are used according to whether the adjusted mean SMB
of the previous decade is positive or negative, and whether
the grid cell is in the north or south of the ice sheet. Eleva-
tion feedback uncertainty can be sampled with different SMB
lapse rates; with careful experimental design this can give a
probabilistic assessment of the effect of elevation feedback
uncertainty on sea level. ISM and GCM uncertainty can also
be sampled by using different models. We present these re-
sults in a companion paper (Edwards et al., 2014).

## 2 Method

We derive the parameterisation from a set of MAR simu-
lations in which the surface elevation is altered (Sect. 2.1).
We try various choices for the parameterisation structure,
judging them by their success in reproducing the SMB re-
sponse in MAR and their flexibility and ease of implementa-

tion (Sect. 2.2). After deciding on the structure, we estimate
probability distributions for the four SMB–elevation gradi-
ents (Sect. 2.3).

## 2.1 Climate simulations

The regional climate model MAR (Fettweis, 2007) has been
adapted for simulating the climate over Greenland, with
full coupling to a complex snow/ice model and relatively
high horizontal resolution (25 km). MAR is one of the few
RCMs (another is RACMO2/GR) that includes the positive
feedback between ice surface albedo and melting (Fettweis,
2007), though this is only partially included because the ice
sheet extent and elevation are constant (there is no change in
the ice–tundra boundary). MAR has been shown to simulate
GrIS SMB quite successfully (e.g. Fettweis et al., 2011).

We use a set of eight simulations, each 20 yr long, in which
MAR is forced at the boundaries by the ECHAM5 GCM
(Roeckner et al., 2003) under the SRES A1B emissions sce-
nario (Nakićenović et al., 2000). Two are control simulations,
using the default ice surface topography based on Bamber
et al. (2001): they are the first two decades (2000–2019, t₁)
and last two decades (2080–2099, t₂) of the MAR ECHAM5-
A1B simulation described by Rae et al. (2012) and Fettweis
et al. (2013). The other six are perturbation experiments,
three for each time period, in which we alter the GrIS sur-
face height. We use three types of height change: uniform
lowering by 50 m (“−50 m”), uniform lowering by 100 m
(“−100 m”), and NonUniform changes (“NonUn”) derived
from a GrIS simulation by Ridley et al. (2005). Ridley et al.
(2005) couple GISM (Greenland Ice Sheet Model: Huy-
brechts and de Wolde, 1999) to the HadCM3 GCM (Gordon
et al., 2000) so that the elevation feedback is included, and
quadruple the atmospheric CO₂ concentrations from prein-
dustrial values. We use the resulting GrIS surface height
change after 140 yr, at which point the ice sheet has lost 10 %
of its original volume. We interpolate these height changes
from the GISM 20 km polar stereographic grid to the MAR
grid, and add them to the default topography over ice sheet
grid cells. The ice sheet area is not changed: no cells are
changed from ice to tundra or vice versa. Any negative height
values that result after applying the changes are set to zero,
to avoid the ice surface being specified below the bedrock.
Our analysis uses the mean of each two-decade simulation,
over which the SMB time series is approximately stationary
(Rae et al., 2012).

Figure 1 shows the default (control) topography and the
height difference between the NonUn and control experi-
ments. Figure 2 shows the SMB changes for the NonUn ex-
periments and Fig. 3 the uniform height change experiments.
These figures show that large decreases in elevation gener-
ally decrease SMB, due to increased melting and decreased
snowfall (Franco et al., 2012). There are two main excep-
tions to this that arise from the complex effects of topogra-
phy on local air circulation and precipitation. In the NonUn
<figure><img src="image_5.png" /></figure>

**Fig. 1.** Left: ice sheet surface elevation in the control experiments. Right: elevation change in the NonUn experiments (NonUn – control). Red dashed line is 77° N.

<figure><img src="image_6.png" /></figure>

**Fig. 2.** Mean SMB change (perturbed minus control) in the NonUn experiments, 2000–2019 (left) and 2080–2099 (right). Red dashed line is 77° N.

experiments, there is an increase in SMB along the western ice sheet margin while there is a thinning of the ice sheet (Fig. 2). Here the lowering of the ice sheet surface dampens the “barrier wind” that brings warm air from the tundra along the ice sheet margin and enhances melting (van den Broeke and Gallée, 1996). In the uniform height change experiments, surface lowering can lead to either a decrease or increase in SMB (Fig. 3): a decrease in elevation exposes ice to warmer air temperatures, which can increase the moisture content of the air and enhance precipitation, but conversely an increase in elevation may cause air to rise and cool, also encouraging precipitation (Fettweis et al., 2005; Franco et al., 2012). These aspects show the importance of using a surface energy balance based RCM, rather than simpler models, to account for such phenomena. The consequences of this complexity for the parameterisation are discussed in Sect. 3.2.

Figure 4 shows SMB responses versus height changes for the two NonUn experiments, with arrows pointing from control to NonUn result, separated into regions north and south

<figure><img src="image_7.png" /></figure>

**Fig. 3.** Mean SMB change (perturbed minus control) in the –100 m (top row) and –50 m (bottom) experiments for the 2000–2019 (left column) and 2080–2099 (right) simulations. Red dashed line is 77° N.

of latitude 77° N (this choice of latitude is explained later). The structure of the data is somewhat similar to that found by Helsen et al. (2011) for RACMO2/GR, with a broadly linear positive relationship below the equilibrium line altitude (ELA: the line at which SMB equals zero) and a negative, weaker relationship above the ELA. The behaviour is linear below the ELA within each time period because we use the simulation mean: in a constant climate, the average melting is approximately proportional to the average temperature, which is approximately proportional to elevation. The behaviour above the ELA, particularly south of 77° N (the majority of the ice sheet), is reminiscent of the complex relationship found between precipitation and height in MAR by Franco et al. (2012).

There is a clear offset between the beginning and end of the century. At a given height, particularly below the ELA, the SMB is lower in the warmer climate at the end of the century. This is partly due to the linear dependence of melting on local temperature in a constant climate (described above), but also to two mechanisms that accelerate the melting and runoff as the climate warms. The first is the positive ice albedo feedback. Bare ice appears each summer after the accumulated winter snowpack melts, and it has a lower
<figure><img src="image_8.png" /><figcaption>**Fig. 4.** Changes in SMB when perturbing the height of the MAR ice sheet from the control topography to the NonUn-based topography for grid cells north (left) and south (right) of 77° N. Arrows point from control to NonUn experiment. Data with height change $|\Delta h| < 25$ m are excluded.</figcaption></figure>

albedo than snow. So in a warming climate the maximum area of bare ice (the ablation zone) increases, and a positive albedo feedback amplifies the warming. MAR has a more realistic, lower albedo for bare ice (around 0.45) than most RCMs and therefore a greater sensitivity to warming; Fetweis et al. (2013) estimate that surface melting increases exponentially with rising temperatures. The second mechanism is the type of precipitation falling on the ice sheet. In the latter part of the century most summer precipitation falls as rain rather than snow, and most of this runs off directly to the ocean rather than accumulating as ice. Both mechanisms accelerate the decrease in SMB as the A1B scenario progresses.

Figure 5 shows the SMB responses versus height changes for all of the perturbation experiments except NonUn $t_2$ (reserved as a test: Sect. 2.3), divided into four partitions of SMB (negative and positive) and region (north and south of 77° N). Each data point shows the SMB response ($\Delta S_i = S_i^{pert} - S_i^{cont}$) versus the height perturbation ($\Delta h_i = h_i^{pert} - h_i^{cont}$) for a given grid cell $i$, so each grid cell can appear up to five times. We exclude the 906 grid cells of the NonUn $t_1$ experiment that have $|\Delta h_i| < 25$ m (see Sect. 2.3). We also exclude cells in which the SMB crosses the ELA between the control and perturbed experiments (i.e. in which the perturbed and control SMB have opposite signs) to make distinct data sets for positive and negative SMB.

Most of the variation in Fig. 5 is from the NonUn simulation, because this has the widest range of height perturbations. The south has a steeper slope, a stronger relationship between $\Delta S$ and $\Delta h$, than the north.

The uniform perturbation experiments are the short vertical bands at $\Delta h = (-50 \text{ m}, -100 \text{ m})$. Most of the uniform experiment data in the south (bottom two subfigures of Fig. 5) show the behaviour we expect: when elevation decreases, SMB decreases (most points are in the bottom left quadrant, where both $\Delta S$ and $\Delta h$ are negative; $\Delta S/\Delta h$ is therefore

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 5.** Scatter plots of SMB change $\Delta S$ versus height change $\Delta h$ for grid cells north (top row) and south (bottom) of 77° N, divided into grid cells with SMB in both the control and perturbed experiments less than zero (left column) and greater than or equal to zero (right). Data with height change $|\Delta h| < 25$ m are excluded.</figcaption></figure>

positively signed). But in the north (top two subfigures), the response is the opposite: when elevation decreases, SMB increases (most points are in the top left quadrant, with positive $\Delta S$ and negative $\Delta h$; $\Delta S/\Delta h$ is negative). This change in response from north to south can be seen in the maps in Fig. 3, particularly for the –100 m experiment (top two subfigures) where the SMB change along the margin is positive (red) north of 77° N and mostly negative (blue) in the south. However, most of the uniform experiment data do lie within the range of the NonUn results.

## 2.2 Parameterisation structure

The parameterisation comprises four “SMB lapse rates”, gradients that characterise a linear relationship between SMB change and surface elevation change. When testing the parameterisation, we use the gradients to adjust the control SMB using the NonUn height changes and compare with the actual NonUn SMB results. In a companion paper we use the parameterisation with several ice sheet models to dynamically adjust projections of future SMB as the GrIS shape evolves (Edwards et al., 2014). The four gradients correspond to the four possible combinations of the grid cell adjusted mean SMB over the past decade being positive or
negative and the grid cell latitude being north or south of 77° N. We estimate these gradients from the ratios of SMB changes to height changes ($\Delta S/\Delta h$) in the surface elevation perturbation experiments. This parameterisation structure is determined by a combination of a priori choices and informal tests.

We choose the structure of our parameterisation with the following aims: to preserve as much of the SMB–elevation relationship in the MAR simulations as possible; to use as few assumptions as possible; to be applicable to any SMB forcing from MAR; and to be simple for the ice sheet modeler to implement. We test the ability of the parameterisation to reproduce the SMB field in the NonUn $t_2$ simulation when applied to the control $t_2$ simulation using the NonUn height changes.

We parameterise the relationship between elevation and mean SMB, using total SMB rather than its individual components (as in Franco et al., 2012) so it is easier to implement in ISMs and requires only one simulated variable as the input forcing. We also choose to parameterise changes in SMB as a function of changes in elevation (in common with Franco et al., 2012), rather than absolute values (as in Helsen et al., 2011). If we were to parameterise the relationship between absolute SMB $S$ and absolute height $h$, using a linear model $S = a + bh$ (e.g. in Fig. 4), we would force the adjusted SMB to lie along a single line lying somewhere between the data from the two time periods $t_1$ and $t_2$, with large uncertainty in the intercept due to the climate dependence of SMB at a given height. Instead, we can parameterise the relationship between SMB changes and height changes, $\Delta S = b\Delta h$, estimating only the gradient $b$. This way SMB can be adjusted up or down the slope apparent in the data, rather than onto a single line with constant intercept. Eliminating the intercept in this way preserves the climate dependence of the SMB–elevation relationship in the MAR simulations, and removes half the unknown parameters. Working with anomalies rather than absolute values is also a standard approach in climate modelling, because the former are thought to be simulated more reliably than the latter.

The adjusted SMB is $S^{\text{adj}} = S^{\text{RCM}} + b\Delta h$, where $S^{\text{RCM}}$ is the original SMB (kg m⁻² a⁻¹), $\Delta h$ the height change (m), and $b$ the SMB–elevation gradient $b = \Delta S/\Delta h$ (kg m⁻³ a⁻¹). More specifically, for a given MAR grid cell in a year $t$ a gradient $b_t$ is used to adjust the control SMB $S_t^{\text{RCM}}$ using the height difference between the NonUn and control experiments, $S_t^{\text{adj}} = S_t^{\text{RCM}} + b_t (h^{\text{NonUn}} - h^{\text{cont}})$. The gradient $b_t$ is selected according to the “reference” SMB and latitude of the grid cell, where the reference is the mean of the adjusted SMB over the previous 10 yr (see Edwards et al., 2014, for more details). Using the adjusted SMB for the reference means the gradient selection evolves as part of the feedback, which helps to make the method more robust with changing climate.

Our height perturbation simulations allow us to derive the gradients directly from SMB responses to height changes for each grid cell, rather than the difference in SMB between grid cells in different locations on the ice sheet (as in Helsen et al., 2011). This is important because the SMB response may be determined by different physical processes due to local topography and atmospheric circulation patterns. Each grid cell $i$ provides an estimate of the gradient $b = \Delta S/\Delta h$ from the SMB change (perturbed SMB minus control SMB, $\Delta S_i = S_i^{\text{pert}} - S_i^{\text{cont}}$) versus the elevation change (perturbed height minus control height, $\Delta h_i = h_i^{\text{pert}} - h_i^{\text{cont}}$). In Fig. 4 these correspond to the arrow slopes; in Fig. 5 they are the y axis values divided by the x axis values.

We choose not to make the gradients a function of grid cell location (Helsen et al., 2011; Franco et al., 2012), other than the north–south divide, to avoid dependence on the MAR grid resolution (Franco et al., 2012) and make the parameterisation as generic as possible. A spatially varying parameterisation would be tailored to the current shape of the ice sheet and the gradients would need to be interpolated for the ISM grid, which could lead to distorting edge effects at discontinuities such as the margin, ELA, and grid cell boundaries (e.g. Franco et al., 2012).

We do not make the gradients a function of climate or time (beginning versus end of the century), because this would restrict our ability to apply the parameterisation to other MAR simulations: for the missing years of the A1B scenario (2020–2079) we could interpolate or otherwise scale the results, but this would be less reliable or applicable for other emissions scenarios and simulations forced by other GCMs.

To guide our choices for other aspects of the parameterisation structure, we consider various methods of estimating and applying the gradients and quantify their relative success in reproducing SMB changes in one of the perturbation experiments. We estimate the gradients from the SMB responses in the two NonUn simulations, then use them to adjust the SMB in the control $t_2$ experiment according to the NonUn height changes. We quantify success by comparing the parameterised cumulative SMB change with the actual results in the NonUn $t_2$ simulation, in terms of both the root mean square error in the spatial pattern and the error in the GrIS total (not shown). We base our decisions on a combination of practical considerations (such as ease of implementation) and these informal sensitivity tests, rather than a systematic optimisation across all possible choices.

Our final gradients are a function of SMB sign (positive/negative) and region (north/south), because these divisions make substantial improvements to the parameterisation while not introducing much complexity when implementing in ISMs. The clear difference in SMB response above and below the ELA has already been discussed (Sect. 2.1). We also choose to divide by region because of the distinct regimes in Fig. 4 in which the north has a shallower gradient and larger intercept than the south. The uniform height
change simulations also indicate that the northern mar-
gin behaves differently (Fig. 3). We test the performance of
north–south divisions in half degree intervals in the range
74–79° N, and also compare with using no division, and find
that 77° N gives the best result.

We test two other functional dependencies for the gradi-
ents: eight divisions in SMB rather than two, and height de-
pendence as well as SMB dependence. The improvements
are not marked enough to justify the extra complexity.

We try three methods for estimating gradients: (a) a lin-
ear model of S versus h, (b) a linear model of ΔS versus
Δh with zero intercept, and (c) a non-parametric method.
We apply each to the four data sets (positive/negative SMB,
north/south), and grid cells with | Δh | < 25 m are excluded.
In method (a), a linear fit of S vs. h estimates the gradi-
ent b in Fig. 4; this is a similar approach to Helsen et al.
(2011), except that we then make the SMB adjustment with
our anomaly method rather than an intercept. In method (b), a
linear fit of ΔS versus Δh estimates the gradient b in Fig. 5; a
zero intercept reflects our expectation that mean SMB change
is zero if there is no height change. In method (c), we use
a non-parametric approach instead of a linear model. This
takes the median of ΔS/Δh ratios (y/x in Fig. 5) as an es-
timate of b. Method (a) is the least successful, and (b) is the
most successful. But we judge that (b) is not an appropri-
ate method, because the fit residuals for grid cells above the
ELA vary systematically as a function of height change. Part
of this may be due to our constraint of a zero intercept, but
the data also clearly have non-linear structure (Fig. 5). The
non-parametric method avoids model assumptions such as
normally distributed fit residuals, allowing us to capture all
the aspects of the MAR response. Our final method is there-
fore based on (c), though we use the full distribution rather
than the median (Sect. 2.3).

Our final parameterisation of the SMB–elevation feedback
is therefore a set of four gradients, **b** = (b<sub>p</sub><sup>N</sup>, b<sub>n</sub><sup>N</sup>, b<sub>p</sub><sup>S</sup>, b<sub>n</sub><sup>S</sup>), that
are used to adjust SMB with a linear model of SMB change
versus elevation change. A gradient is selected from the set
of four according to whether the mean of the adjusted SMB
in the previous decade is positive (p) or negative (n) and
whether the grid cell is north or south of 77° N (N, S). The
gradients are estimated from the ratios ΔS/Δh of grid cells
in the MAR perturbation experiments.

The full SMB–elevation relationship is complex, but our
aim is to create a parameterisation that is straightforward to
implement: we have therefore partially linearised it, by par-
titioning the data four ways and by using a linear model of
SMB adjustment. In the following section our aim is to ac-
count for this approximation with non-parametric assessment
of uncertainties in those linear parameters. With this struc-
ture it is easy to implement the parameterisation, to use it in
forcing an ice sheet model, and to assess the impacts of the
parametric uncertainty (arising in large part from this lineari-
sation) using additional ice sheet model simulations.

2.3 Parameter estimation

We now turn to formal statistical inference to obtain the fi-
nal gradient values. We wish to assess the full uncertainty
in the SMB–elevation relationship rather than using only
tuned (“best estimate”) values or performing ad hoc sensi-
tivity tests. This is particularly important given there are op-
posite sign SMB responses to elevation changes in the simu-
lations. So we estimate full probability distributions for each
of the four gradients (b = b_p^N, b_n^N, b_p^S, b_n^S) using a Bayesian
approach. This also allows us to propagate the probabilis-
tic SMB–elevation feedback uncertainty to predictions of the
GrIS contribution to sea level (Edwards et al., 2014).

We derive initial (“prior”) distributions for the four gradi-
ents using SMB responses from five of the six perturbation
simulations: –100 m, t₁ and t₂; –50 m, t₁ and t₂; and NonUn
t₁. These five experiments thus include parameter estimates
under different climates (t₁ and t₂), different height changes
(–50 m, –100 m, and NonUn), and different locations for a
uniform height change (–50 m, –100 m).

We reserve the final simulation (NonUn t2) as a test of
the parameterisation, reweighting the prior distributions us-
ing the degree of success in reproducing the cumulative sea
level change to obtain updated (“posterior”) distributions. We
choose NonUn t2 because the NonUn height changes span a
wider range and are closer in spatial pattern to those expected
in a warmer climate than the uniform height changes, and be-
cause the SMB signal is larger for t2 than for t1; we are more
concerned that the parameterisation is valid under a warmer
climate than the present day.

This division of simulations allows us to try a wide range
of candidates for parameter values but assign larger weights
to those that match the target we wish to reproduce: the ag-
gregate behaviour of the whole ice sheet.

We use histograms of the ratio of SMB changes to height
changes, ΔS/Δh (Fig. 6), as a basis for our prior distribu-
tions for the four linear gradient values. These are the same
data as in Fig. 5, which shows ΔS versus Δh. Our mini-
mum threshold for the denominator, | Δh | ≥ 25 m, removes
extreme values from the tails of these distributions which sta-
bilises estimation of the ratios. All four distributions show
that SMB is sometimes positively correlated with height,
sometimes negatively correlated. Above the ELA (Figs. 6b
and d) the histograms for b_p^N and b_p^S are very narrow: the
SMB responses for a given height change are small in magni-
tude with little variation. Below the ELA (Figs. 6a and c) the
histograms for b_n^N, b_n^S are much broader, showing the wide
variation in response for different regions of the ice sheet.
These histograms are dominated by the four uniform pertur-
bation simulations.

Each of the four histograms has a different number of grid
cells, so we take equally sized subsets of each to obtain a joint
sample of the gradient set **b**: for each histogram we order the
values of ΔS/Δh and take the 0.5 % to 99.5 % quantile val-
ues in 0.5 % steps, giving 199 samples of the four gradients
<figure><img src="image_11.png" /><figcaption>Fig. 6. Histograms of the ratio ΔS/Δh for grid cells north (top row) and south (bottom) of 77° N, divided into grid cells with SMB in both the control and perturbed experiments less than zero (left column) and greater than or equal to zero (right). Data with height change | Δh | < 25 m are excluded.</figcaption></figure>

(b<sub>p</sub><sup>N</sup>, b<sub>n</sub><sup>N</sup>, b<sub>p</sub><sup>S</sup>, b<sub>n</sub><sup>S</sup>). These prior distributions are shown in light grey in Fig. 7.

We use each of these 199 prior estimates of the gradient set to adjust the control SMB in 2080–2099 according to the NonUn height change, and assess their success in reproducing the target NonUn t<sub>2</sub> experiment. Each gradient set is used to calculate a spatial pattern of cumulative SMB change and the corresponding total GrIS cumulative sea level contribution.

We simplify the statistical modelling by choosing comparisons so that the differences (discrepancies) between the adjusted and target SMB at each location are approximately i.i.d. (independent and identically distributed) in space. We make the comparisons approximately independent by “thinning” (Rougier and Beven, 2013), using only every 5th grid cell (125 km spacing). This spacing removes spatial correlation: an empirical variogram of the thinned discrepancies is flat for all lengths up to around 800 km (the width of the ice sheet). We assume the discrepancies are identically distributed in space, that is, that the model is equally likely to match the target at every location. We also assume that the discrepancies are normally distributed. In the absence of further information and as a first attempt to describe parameterisation uncertainty, these choices and assumptions allow us

<figure><img src="image_12.png" /><figcaption>Fig. 7. Prior (light grey) and posterior (dark grey) distributions of the four value gradient set, b = (b<sub>p</sub><sup>N</sup>, b<sub>n</sub><sup>N</sup>, b<sub>p</sub><sup>S</sup>, b<sub>n</sub><sup>S</sup>) for regions north (b<sup>N</sup>, top row) and south (b<sup>S</sup>, bottom) of 77° N, and SMB less than zero (b<sub>n</sub>, left column) and greater than or equal to zero (b<sub>p</sub>, right).</figcaption></figure>

to avoid the difficult task of modelling the spatial correlation and variation of the discrepancies.

These assumptions translate to a simple metric for assessing the gradient estimates. The scoring, “likelihood”, function is a multivariate (for multiple locations) independent Gaussian with constant variance; the exponent is the sum of squared differences between the adjusted SMB and the target SMB over the subsampled grid cells (independent: a product of Gaussians) divided by the “discrepancy variance” σ<sup>2</sup> (identically distributed: constant variance). The multiplicative constant is discarded due to normalisation later. So the score s<sub>j</sub> for the jth of 199 samples of b is

$$s_j = \exp \left[ -\frac{1}{2\sigma^2} \sum_i (f_i^j - z_i)^2 \right], \quad (1)$$

where f is the adjusted SMB, z the target SMB, and i the grid cell index. The discrepancy variance is a parameter that represents how closely we expect the parameterised SMB to match the target; our choice is discussed below.

The weight given to each gradient set is the normalised score, w<sub>j</sub> = s<sub>j</sub> / Σs<sub>j</sub>. Note that a single weight is calculated for each gradient set b, rather than individual weights for each of the four components. The most successful
<table><caption>**Table 1.** The 2.5% quantile, best estimate, and 97.5% quantile estimates of the SMB–elevation gradients in kg m⁻³ a⁻¹, below (SMB < 0) and above (SMB ≥ 0) the ELA, for regions north and south of 77° N.</caption><thead><tr><th></th><th>Region</th><th>2.5%</th><th>Best estimate</th><th>97.5%</th></tr></thead><tbody><tr><td rowspan="2">SMB &lt; 0</td><td>North</td><td>-0.22</td><td>0.56</td><td>1.33</td></tr><tr><td>South</td><td>1.03</td><td>1.91</td><td>2.61</td></tr><tr><td rowspan="2">SMB ≥ 0</td><td>North</td><td>-0.03</td><td>0.09</td><td>0.23</td></tr><tr><td>South</td><td>-0.07</td><td>0.07</td><td>0.59</td></tr></tbody></table>

("maximum likelihood") gradient set $\tilde{b}$ has the smallest sum of squared differences and therefore the largest weight. We calculate posterior distributions for the four components of $b$ by reweighting the prior distributions with the normalised weights. We estimate probability densities from the histograms with kernel density estimates and use these to estimate the modes of the posterior distributions, which are our best estimates of the gradients. As we are in a Bayesian framework, our uncertainties are expressed as "credibility intervals" rather than confidence intervals. We estimate 95% credibility intervals with bootstrapping: we resample 100 000 times from the 199 gradient values (with replacement, using the normalised weights), smooth these with the same bandwidth, and estimate the 2.5% and 97.5% quantiles.

Our statistical framework requires minimal choices: the form of the likelihood function; the spacing for the subsampling; and a value for the discrepancy variance. We also choose to set the bandwidth (standard deviation of the smoothing) for the kernel density estimation because the automatically chosen value (Silverman, 1986) does not seem to adequately resolve the distribution shapes. We test various options and make our final choices with the following considerations: thinning so that the discrepancies appear approximately uncorrelated in space; the variance $\sigma^2$ chosen such that the weights are not concentrated on a small number of gradient estimates and most or all of the discrepancies for the maximum likelihood parameterisation ($\tilde{b}$) are in the range $\pm 3\sigma$ (Pukelsheim, 1994); and the posterior distribution of total GrIS sea level contribution is close to the target. We choose the smoothing bandwidth so that the density profile captures the main features of the histogram. Our final choices are a Gaussian likelihood function; subsampling distance 5 grid cells (125 km); discrepancy variance $\sigma^2 = (20 \times 10^3 \text{ Gt})^2$; and bandwidths 0.15 kg m⁻³ a⁻¹ for gradients below the ELA and 0.05 kg m⁻³ a⁻¹ above the ELA. Sensitivity tests for these choices are described in the next section.

## 2.4 Results

Figure 8 shows the adjusted cumulative SMB from the maximum likelihood parameterisation ($\tilde{b}$) and the target. The maximum likelihood gradient set reproduces the target well

<figure><img src="image_13.png" /><figcaption>**Fig. 8.** Cumulative SMB change at the end of the NonUn 2080–2099 simulation: (left) target MAR simulation (perturbed minus control) and (right) result from maximum likelihood gradient set $\tilde{b}$ applied to the NonUn height change (adjusted minus control). Red dashed line is 77° N.</figcaption></figure>

in most areas, but cannot reproduce the SMB increases with decreasing elevation along the western and southeastern ice sheet margins. Figure 9 shows the discrepancies between the two for all grid cells and the subset used for the likelihood calculation. Most of the discrepancies are small over the ice sheet interior and larger at the margin.

Figure 7 shows the posterior distributions (dark grey) for the four gradients; Table 1 gives the best estimates and 95% credibility intervals. The posterior distributions are mostly positive, with much larger gradients below the ELA, particularly in the south, than above. Most of the distributions are fairly symmetric, except the south above the ELA which has a low best estimate and a long tail of larger values. The weighting has a particularly strong effect for grid cells below the ELA, drastically narrowing the distributions: effectively the likelihood scoring gives high weights to the gradient estimates derived from the NonUn 2000–2019 experiment (large, positive values), rather than the uniform height change experiments (small, positive and negative values), because these are most successful in reproducing the patterns of change in the NonUn 2080–2099 experiment.

We can apply the same weights to the total GrIS cumulative sea level contributions for each sample of the gradient set (Fig. 10). The prior distribution is centred close to zero: the prior estimate of the elevation feedback is that it has no net effect. The update narrows and shifts the posterior distribution so that it is centred on the target, a positive contribution from the feedback.

We test the sensitivity of the results to the elevation threshold and statistical modelling choices. Varying the threshold (from the default 25 m) between 10 m and 50 m in 5 m intervals changes the results by no more than 0.02 kg m⁻³ a⁻¹, and in most cases 0.01 kg m⁻³ a⁻¹ or zero, for the majority of the gradient best estimates and CI bounds. The exceptions are the best estimates in the south ($b_S^n, b_S^p$) and the upper CI
<figure><img src="image_14.png" /></figure>

**Fig. 9.** Cumulative SMB change at the end of the NonUn 2080–2099 simulation: (left) error in the maximum likelihood gradient set $\tilde{b}$ applied to the NonUn height change (adjusted minus perturbed), and (right) the subset of those grid cells used in the calculation of the weights (discrepancies $f_i^j - z_i$ in Eq. 1). Red dashed line is 77° N.

bound for the latter ($b_S^P$), where the changes are in the range 0.06–0.12 kg m$^{-3}$ a$^{-1}$. These correspond to a 5 % change for $b_S^P$; for the smaller gradient $b_S^P$ the fractional changes are larger, but the changes plateau (i.e. the gradient estimates stabilise) above 25–30 m, as intended with the use of the threshold.

We try substituting the Gaussian likelihood with a Cauchy (Student's t distribution with one degree of freedom; very heavy tailed), scaled to match a Gaussian at the 25 % and 75 % quantiles. Our motivation is that the histogram of discrepancies for the maximum likelihood gradient set is fairly sharply peaked. The effect of a Cauchy likelihood is to distribute the weights over a much smaller number of gradient sets, which drastically narrows the posterior distributions. If we also reduce the bandwidths to match these narrower distributions (from 0.15 to 0.05 kg m$^{-3}$ a$^{-1}$ below the ELA and from 0.05 to 0.03 kg m$^{-3}$ a$^{-1}$ above), the CI widths decrease by 54–84 %, and the best estimates increase by 13–38 % for three of the gradients and 186 % (from 0.07 to 0.20 kg m$^{-3}$ a$^{-1}$) for $b_p^S$. Because the weights are so concentrated, and we wish to be conservative with uncertainty estimates, we choose the Gaussian likelihood. An alternative approach would be to set a larger discrepancy variance for the ice sheet margin grid cells than the interior, though one might be less confident in assigning the value of two uncertain parameters rather than one.

The discrepancies for the maximum likelihood parameterisation are all within ±1.5σ, which indicates that our discrepancy variance is too large; on the other hand, reducing σ concentrates the weights on a smaller number of gradient estimates, leading to narrower posterior distributions and 95 % CIs. Changing σ from 20 to 15 or 25 Gt does not affect the best estimates much (0–11 %) except for the small-valued $b_p^S$ (43 %). Increasing or decreasing σ by 5 Gt has the effect of

<figure><img src="image_15.png" /></figure>

**Fig. 10.** Prior (blue) and posterior (red) distributions of cumulative change in sea level at 2099 using parameterised elevation feedback for height changes in the NonUn simulation. The target is the result from the NonUn 2080–2099 experiment (vertical black line).

increasing or decreasing the CI widths by 12–24 %. Decreasing σ to 15 or 10 Gt broadens the discrepancies to about 2σ, but concentrates the weights rather more. Again, we err on the side of conservatism in our choice.

In the grid cell sampling (required for independence), using different spacing does not have a monotonic effect on the results. Decreasing the spacing from 5 grid cells to 4 (100 km) or increasing it to 6 (150 km) both have the effect of decreasing most best estimates and CI widths. This indicates it is not a problem of using too short a correlation length (violating the independence assumption) but of sensitivity to the grid cell sampling, most likely at the ice sheet margin. Of these three choices, the 5 cell spacing produces the best match to the cumulative sea level change in Fig. 10; in other words, both the 4 and 6 cell spacings concentrate the weights on smaller gradients (smaller SMB adjustments), which match the target spatial pattern well for the particular sampled cells but perform poorly for the ice sheet total using all grid cells. We alter the offset of the sampling, which also has a non-monotonic effect. Shifting both the longitudinal and latitudinal offsets by –3 cells gives a small decrease in the best estimate (0 to –4 %), while offsets of –2, –1 and +1 all give higher best estimates (15–34 %, except $b_p^S$ 71–100 %). The effect on CI width is also mixed; the largest effect is on $b_p^S$, up to 26 %.

Using a larger discrepancy variance for the margin than the interior would reduce the sensitivity of the results to sampling, because the margin grid cells would have less effect on the likelihood value.
Using automatically set bandwidths gives wider CIs, but oversmooths the distributions, for SMB < 0 and narrower CIs, but undersmooths the distributions, for SMB ≥ 0. Changing our bandwidths from 0.15 to 0.1 or 0.2 kg m⁻³ a⁻¹ below the ELA and from 0.05 to 0.03 or 0.07 kg m⁻³ a⁻¹ above affects the CI widths by small amounts ±0–6%, except for $b_p^N$ (–23 to 27%).

## 3 Discussion

### 3.1 Advantages and strengths

The main advantage of parameterising the GrIS SMB–elevation feedback is that it allows us to force ISMs with SMB simulated by the MAR RCM, which has a more physically realistic representation of the processes than simple schemes such as the PDD, while incorporating the feedback.

The parameterisation can also be used without ISMs to make a first order adjustment to SMB, improving projections such as Rae et al. (2012) and Fettweis et al. (2013) by incorporating the elevation feedback (in effect omitting only the dynamical ice response): for this, the SMB in a given year for each grid cell can be converted to an ice-equivalent height change. A third use is adjusting low resolution SMB fields to the observed surface elevation, for better comparisons with observations or inputs to ISMs (as do Franco et al., 2012).

We have confidence in our parameterisation due, for example, to the similarity in patterns between the maximum likelihood result and the target (Fig. 8), and the centring of the posterior sea level distribution on the target (Fig. 10), and have quantified the effects of the complex non-linear responses in MAR on the feedback uncertainty.

There are several advantages to our approach relative to the parameterisations by Helsen et al. (2011) and Franco et al. (2012). The first relate to our RCM simulations, in which the relationship between SMB and height appears to be more complex (e.g. Fig. 4). For Helsen et al. (2011), this may be partly due to the different schemes in MAR and RACMO2/GR, but in general it is due to our use of simulations in which both the surface elevation and climate boundary conditions are altered. Franco et al. (2012) alter the grid resolution, which produces local changes to elevation; Helsen et al. (2011) do not use altered topography to estimate their SMB gradients, though they do to assess their performance. Neither force the RCM with a global climate different to the present day. Altering the elevation means there is no need for a “space-for-time” substitution. This improves the relevance and robustness of the parameterisation because it is based on results from height changes at a given location, rather than height changes across different spatial locations. We found it is also important to use a wide range of height perturbations, that is, that 50–100 m changes are not sufficient to explore the relationship. It may also be important to apply height changes with the spatial pattern expected

under climate change (NonUn) rather than a uniform or fractional (Hakuba et al., 2012) lowering, because the effects on local atmospheric circulation are potentially quite different. We see it is important, particularly for RCMs that include the albedo feedback, to assess the elevation feedback under different global climate conditions, rather than studying one climate era (Franco et al., 2012) or correcting the ice sheet elevation for other climates using a temperature lapse rate (Helsen et al., 2011).

The second set of advantages relate to our parameterisation structure. Using only a gradient (in common with Franco et al., 2012), rather than a gradient and intercept (Helsen et al., 2011), is more robust because it minimises the problem of the climate-dependent offset. In other words, parameterising the relationship between SMB changes and height changes, rather than absolute values, retains more information about the response. A further aspect of flexibility is our choice to estimate the gradients with a non-parametric method (no assumed functional form) rather than a linear model as both Helsen et al. (2011) and Franco et al. (2012) do. Furthermore our parameterisation is very flexible because, unlike the previous studies, it does not depend on spatial location (other than the north–south divide) so it does not depend on the RCM resolution or require interpolation to the ISM grid, and is easy to implement.

The third advantage relates to parameter assessment. We estimate the gradients within a formal probabilistic framework. This allows us to provide not only a best estimate parameterisation but the full probability distributions, so that ISMs can be used to explore the effect of this uncertainty on the GrIS contribution to sea level and express these as credibility intervals.

### 3.2 Limitations and further work

This is a parameterisation of the SMB-feedback response in an RCM, not of the real world. We have not attempted to estimate the parametric or structural uncertainty of MAR, and have not used an observational constraint. One approach to explore structural uncertainty would be to compare with parameterisations derived for other RCMs. Assessing MAR parametric uncertainty would require a perturbed parameter ensemble such as the 11-member HadRM3 ensemble of Murphy et al. (2009), which is very computationally expensive. We could incorporate observations into the elevation feedback by using them, rather than the target simulation, to calculate the likelihood. We would have to take care that a parameterisation based on observed SMB changes would give a coherent result when applied to RCM simulations. However, the effect of a first-order SMB adjustment on sea level is negligible for a present day ERA-INTERIM forced simulation (not shown), so observational constraints might in any case be of limited use.

We use MAR because it is the most successful of the three RCMs presented by Rae et al. (2012) at reproducing
the current SMB of the GrIS. If we were to study additional RCMs, we would derive the parameterisation separately for each. Our preliminary assessment of another RCM used by Rae et al. (2012), HIRHAM, indicates that the SMB response to height is much more linear and less variable than in MAR, most likely because fewer processes are incorporated such as the albedo feedback. Uncertainty in SMB projections are generally thought to be dominated by the choice of GCM, rather than RCM (e.g. Rae et al., 2012); two GCMs are used for the projections presented by Edwards et al. (2014).

We note that the parameterisation may not be robust for
changes in elevation and SMB much greater than those in
the NonUn A1B 2080–2099 simulation.

We could parameterise each component of SMB sepa-
rately (Franco et al., 2012) or make the parameterisation
structure more complicated in other ways described, but this
would have hindered our aim to test the results from the pa-
rameterisation in several ISMs (Edwards et al., 2014).

Our estimation of the gradients is non-parametric, but our
adjustment of SMB with these gradients uses a linear model
with zero intercept. This linear model is an approximation
we use to reduce the complexity of the parameterisation. Fig-
ure 4 indicates that the relationship is not quite linear above
or below the ELA: for example, the gradient in the north is
slightly shallower at the lowest elevations. Figure 5 indicates
that the intercept may be non-zero below the ELA. A more
complex parameterisation could account for these departures
from our model, though it would be harder to implement.

Figures 2 and 3 show that a surface lowering can often lead
to an increase in SMB. This is particularly the case in the uni-
form elevation change simulations (Fig. 3) for the north and
east, for small elevation changes, and for the beginning of the
century. It is also apparent in the NonUn simulations (Fig. 2)
along the western margin. This behaviour is likely to derive
from the precipitation component of SMB, which has a com-
plex, non-linear relationship with surface elevation (Franco
et al., 2012). The maximum likelihood gradient set does not
reproduce the SMB increases with decreasing height in the
west and southeast (Fig. 8). The probability distributions do
incorporate this variation by including the full range of re-
sponses: in other words, other samples from **b** give differ-
ent correction patterns depending on whether the individual
components are positive or negative (Edwards et al., 2014).
But within an individual ISM simulation the four compo-
nents of **b** are held constant. One way to represent this com-
plex behaviour more fully would be with a stochastic pa-
rameterisation, in which the gradients are randomly sampled
from the distributions through the simulation rather than held
constant. This would incorporate the effect of both positive
and negative values of each gradient within a single simu-
lation rather than separate simulations (as in Edwards et al.,
2014). However, this would require much more complex sta-
tistical modelling to describe the spatial and temporal corre-
lation structure of the gradients, and more complex imple-
mentation.

We exclude grid cells with opposite sign SMB in the con-
trol and perturbation simulations when estimating the gradi-
ents, because of our division at the ELA. (This exclusion only
applies to estimating the gradients, not to applying them:
when adjusting the MAR SMB, a grid cell may be above
the ELA before and below after). The arrows that cross the
ELA in Fig. 4 indicate that this filtering may tend to remove
smaller than average gradients from the below ELA sample
and larger gradients from the above ELA sample. There is
no significant change in albedo in grid cells that remain on
one side of the ELA or the other, so these show a linear re-
lationship between SMB and elevation. Non-linearity occurs
mainly for grid cells that are above the ELA (where no bare
ice appear in summer) in the control simulation and move be-
low the ELA with a new elevation, or vice versa, but these are
not included in the analysis. This exclusion might therefore
lead to an underestimate of the uncertainty.

Boundaries can lead to unexpected edge effects. If the ice
sheet were to retreat past 77° N then the parameterisation
would shift to using only the south values (i.e. a larger re-
sponse). However, as we have discussed, we would not rec-
ommend using this parameterisation for elevation and SMB
changes far beyond our simulations.

The gradients change stepwise across the north–south
boundary at 77° N. In principle, this could be smoothed out
with a soft transition in a slightly more complex implemen-
tation. We include only one boundary rather than several
to minimise edge effects. It might be useful to use a fur-
ther regional division, west–east at around 40° W, because
MAR projects different precipitation responses to a warming
climate: along the eastern coast snowfall tends to increase,
while along the western coast summer precipitation begins
to fall as rain.

Our choice of a 25 m threshold was made to stabilise
the gradient estimates. Most results are not sensitive to this
choice, but we note that a different threshold would in partic-
ular alter the best estimate and CI upper bound of $b_S^0$.

Using the mean of each 20 yr simulation might lead to
an underestimate of uncertainty by averaging over temporal
variability. However, this variability is incorporated when the
parameterisation is applied to an annual SMB time series, so
care would have to be taken to avoid double counting.

Finally, the choices of the structure and parameter estima-
tion depend on the approximations and the prioritisation of
aspects described in Sects. 2.2 and 2.3. Different choices for
the statistical modelling might be justifiable, and future work
could explore this more thoroughly. As mentioned, an alter-
native choice might be to set a larger discrepancy variance
for the ice sheet margin grid cells than the interior. The sub-
set of data chosen to ensure independent discrepancies could
also be selected by hand rather than regularly spaced, with
the aim of choosing the most informative cells: for example,
picking cells along the margin and a smaller number from the
interior. Aggregation of data is another possible method for
removing correlation (Rougier and Beven, 2013). Different
statistical modelling choices would certainly be appropriate if parameterising the elevation feedback for a different RCM.

## 4 Conclusions

Surface mass balance of the GrIS can be modelled with the sophisticated, physically based energy balance schemes available in some RCMs, but this is usually at the expense of including the elevation feedback. To include the feedback requires coupling the RCM with an ISM, but this is computationally expensive and technically challenging, which effectively precludes the exploration of uncertainties in the structure and parameter values of the ISM, and in the elevation feedback in the RCM. The only way to incorporate the physical modelling of SMB processes and elevation feedback while also exploring these model uncertainties is with a parameterisation such as the one presented here.

We estimate the SMB–elevation feedback separately below and above the ELA and for regions north and south of 77°N from a set of MAR simulations in which we alter the ice sheet surface elevation. We make advances on previous parameterisations, in particular by including probabilistic assessment of the parameters. This demands statistical modelling choices that must balance adequate representation of the processes in MAR with ease of implementation for widespread use by ice sheet modellers. In general we approach this with linearisation of the processes (for simplicity of use) combined with Bayesian inference of parameter probability distributions (to capture the range of uncertainties arising from linearisation). The maximum likelihood parameterisation is successful in representing the cumulative sum, and most of the spatial pattern, of MAR SMB changes at the end of the century under the A1B scenario and non-uniform elevation changes. Best estimates and 95% credibility intervals for the parameters provide ice sheet modellers with the opportunity to fully explore the non-linear behaviour. In Edwards et al. (2014) we propagate the probabilistic uncertainties presented here to future sea level projections for the GrIS.

*Acknowledgements.* This work was supported by funding from the ice2sea programme from the European Union 7th Framework Programme, grant number 226375. Ice2sea contribution number 120. S. Price, M. Hoffman, and M. Perego were supported by the US Department of Energy (DOE) Office of Science, Advanced Scientific Computing Research and Biological and Environmental Research programs. M. Hoffman was partially supported by the Center for Remote Sensing of Ice Sheets at the University of Kansas through US National Science Foundation grant ANT-0424589. Thanks to Jonty Rougier for helpful discussions. We thank Fuyuki Saito and two anonymous referees for their helpful reviews.

Edited by: E. Larour

## References

Bamber, J. L., Layberry, R. L., and Gogenini, S. P.: A new ice thickness and bed data set for the Greenland ice sheet 1: measurement, data reduction, and errors, J. Geophys. Res., 106, 33773–33780, 2001.

Edwards, T. L., Fettweis, X., Gagliardini, O., Gillet-Chaulet, F., Goelzer, H., Gregory, J. M., Hoffman, M., Huybrechts, P., Payne, A. J., Perego, M., Price, S., Quiquet, A., and Ritz, C.: Effect of uncertainty in surface mass balance-elevation feedback on projections of the future sea level contribution of the Greenland ice sheet, The Cryosphere, 8, 195–208, doi:10.5194/tc-8-195-2014, 2014.

Ettema, J., van den Broeke, M. R., van Meijgaard, E., van de Berg, W. J., Bamber, J. L., Box, J. E., and Bales, R. C.: Higher surface mass balance of the Greenland ice sheet revealed by high-resolution climate modeling, Geophys. Res. Lett., 36, L12501, doi:10.1029/2009GL038110, 2009.

Fettweis, X.: Reconstruction of the 1979–2006 Greenland ice sheet surface mass balance using the regional climate model MAR, The Cryosphere, 1, 21–40, doi:10.5194/tc-1-21-2007, 2007.

Fettweis, X., Gallée, H., Lefebre, L., and van Ypersele, J.-P.: Greenland surface mass balance simulated by a regional climate model and comparison with satellite derived data in 1990–1991, Clim. Dynam., 24, 623–640, 2005.

Fettweis, X., Tedesco, M., van den Broeke, M., and Ettema, J.: Melting trends over the Greenland ice sheet (1958–2009) from spaceborne microwave data and regional climate models, The Cryosphere, 5, 359–375, doi:10.5194/tc-5-359-2011, 2011.

Fettweis, X., Franco, B., Tedesco, M., van Angelen, J. H., Lenaerts, J. T. M., van den Broeke, M. R., and Gallée, H.: Estimating the Greenland ice sheet surface mass balance contribution to future sea level rise using the regional atmospheric climate model MAR, The Cryosphere, 7, 469–489, doi:10.5194/tc-7-469-2013, 2013.

Franco, B., Fettweis, X., Lang, C., and Erpicum, M.: Impact of spatial resolution on the modelling of the Greenland ice sheet surface mass balance between 1990–2010, using the regional climate model MAR, The Cryosphere, 6, 695–711, doi:10.5194/tc-6-695-2012, 2012.

Goelzer, H., Huybrechts, P., Fürst, J. J., Nick, F. M., Andersen, M. L., Edwards, T. L., Fettweis, X., Payne, A. J., and Shannon, S. R.: Sensitivity of Greenland ice sheet projections to model formulations, J. Glaciol., 59, 733–749, doi:10.3189/2013JoG12J182, 2013.

Gordon, C., Cooper, C., Senior, C. A., Banks, H., Gregory, J. M., Johns, T. C., Mitchell, J., and Wood, R. A.: The simulation of SST, sea ice extents and ocean heat transports in a version of the Hadley Centre coupled model without flux adjustments, Clim. Dynam., 16, 147–168, 2000.

Hakuba, M. Z., Folini, D., Wild, M., and Schär, C.: Impact of Greenland’s topographic height on precipitation and snow accumulation in idealized simulations, J. Geophys. Res., 117, D09107, doi:10.1029/2011JD017052, 2012.

Hanna, E., Huybrechts, P., Cappelen, J., Steffen, K., Bales, R. C., Burgess, E., McConnell, J. R., Peder Steffensen, J., van den Broeke, M., Wake, L., Bigg, G., Griffiths, M., and Deniz Savas, D.: Greenland Ice Sheet surface mass balance 1870 to 2010 based on Twentieth Century Reanalysis, and links
with global climate forcing, J. Geophys. Res., 116, D24121, doi:10.1029/2011JD016387, 2011.

Helsen, M. M., van de Wal, R. S. W., van den Broeke, M. R., van de Berg, W. J., and Oerlemans, J.: Coupling of climate models and ice sheet models by surface mass balance gradients: application to the Greenland Ice Sheet, The Cryosphere, 6, 255–272, doi:10.5194/tc-6-255-2012, 2012.

Huybrechts, P. and de Wolde, J.: The dynamic response of the Greenland and Antarctic ice sheets to multiple-century climatic warming, J. Climate, 12, 2169–2188, 1999.

IPCC: Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change, edited by: Solomon, S., Qin, D., Manning, M., Chen, Z., Marquis, M., Avery, K. B., Tingor, M., and Miller, H. L., Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 2007.

Janssens, I. and Huybrechts, P.: The treatment of meltwater retention in mass-balance parameterizations of the Greenland ice sheet, Ann. Glaciol., 31, 133–140, 2000.

Murphy, J. M., Sexton, D. M. H., Jenkins, G. J., Booth, B. B. B., Brown, C. C., Clark, R. T., Collins, M., Harris, G. R., Kendon, E. J., Betts, R. A., Brown, S. J., Humphrey, K. A., McCarthy, M. P., McDonald, R. E., Stephens, A., Wallace, C., Warren, R., Wilby, R., and Wood, R. A.: Climate change projections, in: UK Climate Projections, no. 2, Met Office Hadley Centre, Exeter, UK, 1–194, 2009.

Nakićenović, N., Alcamo, J., Davis, G., and de Vries, B.: Special report on emissions scenarios: a special report of Working Group III of the Intergovernmental Panel on Climate Change, Intergovernmental Panel on Climate Change, 2000.

Pukelsheim: The 3-sigma-rule, The American Statistician, 48, 88–91, 1994.

Rae, J. G. L., Aðalgeirsdóttir, G., Edwards, T. L., Fettweis, X., Gregory, J. M., Hewitt, H. T., Lowe, J. A., Lucas-Picher, P., Mottram, R. H., Payne, A. J., Ridley, J. K., Shannon, S. R., van de Berg, W. J., van de Wal, R. S. W., and van den Broeke, M. R.: Greenland ice sheet surface mass balance: evaluating simulations and making projections with regional climate models, The Cryosphere, 6, 1275–1294, doi:10.5194/tc-6-1275-2012, 2012.

Ridley, J., Huybrechts, P., Gregory, J., and Lowe, J.: Elimination of the Greenland Ice Sheet in a High CO₂ Climate, J. Clim., 18, 3409–3427, 2005.

Robinson, A., Calov, R., and Ganopolski, A.: An efficient regional energy-moisture balance model for simulation of the Greenland Ice Sheet response to climate change, The Cryosphere, 4, 129–144, doi:10.5194/tc-4-129-2010, 2010.

Roeckner, E., Bäuml, G., Bonaventura, L., Brokopf, R., Esch, M., Giorgetta, M., Hagemann, S., Kirchner, I., Kornblueh, L., Manzini, E., Rhodin, A., Schlese, U., Schulzweida, U., and Tompkins, A.: The atmospheric general circulation model ECHAM 5, Part I: Model description, Max-Planck-Institut für Meteorologie, Hamburg, Germany, 126 pp., 2003.

Rougier, J. C. and Beven, K. J.: Model limitations: the sources and implications of epistemic uncertainty, in: Risk and Uncertainty Assessment for Natural Hazards, edited by: Rougier, J. C., Sparks, R. S. J., and Hill, L. J., Cambridge University Press, Cambridge, 40–63, ISBN 978-1-107-00619-5, 2013.

Shepherd, A., Ivins, E. R., Geruo A., Barletta, V. R., Bentley, M. J., Bettadpur, S., Briggs, K. H., Bromwich, D. H., Forsberg, R., Galin, N., Horwath, M., Jacobs, S., Joughin, I., King, M. A., Lenaerts, J. T. M., Li, J., Ligtenberg, S. R. M., Luckman, A., Luthcke, S. B., McMillan, M., Meister, R., Milne, G., Mouginot, J., Muir, A., Nicolas, J. P., Paden, J., Payne, A. J., Pritchard, H., Rignot, E., Rott, H., Sandberg Sørensen, L., Scambos, T. A., Scheuchl, B., Schrama, E. J. O., Smith, B., Sundal, A. V., van Angelen, J. H., van de Berg, W. J., van den Broeke, M. R., Vaughan, D. G., Velicogna, I., Wahr, J., Whitehouse, P. L., Wingham, D. J., Yi, D., Young, D., and Zwally, H. J.: A Reconciled Estimate of Ice-Sheet Mass Balance, Science, 338, 1183–1189, 2012.

Silverman, B. W.: Density Estimation, Chapman and Hall, London, 1986.

Stone, E. J., Lunt, D. J., Rutt, I. C., and Hanna, E.: Investigating the sensitivity of numerical model simulations of the modern state of the Greenland ice-sheet and its future response to climate change, The Cryosphere, 4, 397–417, doi:10.5194/tc-4-397-2010, 2010.

Tarasov, L. and Peltier, W. R.: Greenland Glacial History and Local Geodynamic Consequences, Geophys. J. Int., 150, 198–229, 2002.

van de Berg, W. J., van den Broeke, M., Ettema, J., van Meijgaard, E., and Kaspar, F.: Significant contribution of insolation to Eemian melting of the Greenland ice sheet, Nat. Geosci., 4, 679–683, 2011.

van de Wal, R. S. W.: Mass-balance modelling of the Greenland ice sheet: a comparison of an energy-balance and a degree-day model, Ann. Glaciol., 23, 36–45, 1996.

van den Broeke, M. R. and Gallée, H.: Observation and simulation of barrier winds at the western margin of the Greenland ice sheet, Q. J. Roy. Meteor. Soc., 122, 1365–1383, 1996.

Vernon, C. L., Bamber, J. L., Box, J. E., van den Broeke, M. R., Fettweis, X., Hanna, E., and Huybrechts, P.: Surface mass balance model intercomparison for the Greenland ice sheet, The Cryosphere, 7, 599–614, doi:10.5194/tc-7-599-2013, 2013.