SeaRISE experiments revisited : potential sources of spread in multi-model projections of the Greenland ice sheet

The present paper revisits the future surfaceclimate experiments on the Greenland ice sheet proposed by the Sea-level Response to Ice Sheet Evolution (SeaRISE; Bindschadler et al., 2013) study. The projections of the different SeaRISE participants show dispersion, which has not been examined in detail to date. A series of sensitivity experiments are conducted and analyzed using the ice-sheet model for integrated Earth-system studies (IcIES) by replacing one or more formulations of the model parameters with those adopted in other model(s). The results show that large potential sources of the dispersion among the projections of the different SeaRISE participants are differences in the initialization methods and in the surface mass balance methods, and both aspects have almost equal impact on the results. The treatment of ice-sheet margins in the simulation has a secondary impact on the dispersion. We conclude that spinning up the model using fixed topography through the spin-up period while the temperature is allowed to evolve according to the surface temperature history is the preferred representation, at least for the experiment configuration examined in the present paper. A benchmark model experimental setup that most of the numerical models can perform is proposed for future intercomparison projects, in order to evaluate the uncertainties relating to pure ice-sheet model flow characteristics.


Introduction
Numerical modeling is an important technique for projecting the response of ice sheets to climate change (e.g., Huybrechts and de Wolde, 1999).Each of the processes simulated in ice-sheet experiments has a degree of uncertainty associated with it, and thus the final output may sometimes have significant dispersion among possible combinations of the methods used to represent them.Multi-model intercomparison over a standardized protocol of numerical experiments is a typical approach for evaluating the uncertainties in model projections.Several intercomparison experiments have been previously performed with focus on various topics, in particular on the behavior of the Greenland ice sheet under future climate changes.
As numerical models have become increasingly complex, it has become more difficult to examine the sensitivity to all uncertainties in all possible model formulations, both numerical and physical.Multi-model intercomparison is an effective, although not perfect, procedure for evaluation of model uncertainties.Greve and Herzfeld (2013) performed sensitivity studies of 500-year projections of the Greenland ice sheet under two scenarios, the AR4 climate scenario and doubled Published by Copernicus Publications on behalf of the European Geosciences Union.F. Saito et al.: SeaRISE revisited basal sliding, using two different numerical ice-sheet models.The models differ not only in the numerical and physical representation of ice-sheet dynamics but also in the method used to compute the surface mass balance from surface temperatures.Despite the differences, a common result is obtained, showing a larger sensitivity to climate warming than to a doubling of the basal sliding.Herzfeld et al. (2012) studied the sensitivity of Greenland ice-sheet projections to the regional updating of the bedrock topography for some glaciers, also using two different numerical ice-sheet models.Both models show significant impact in the response to the doubled sliding scenario by just changing a limited area of bedrock topography.Shannon et al. (2013) used four numerical icesheet models to evaluate the effect of enhanced basal sliding driven by surface runoff on 200 years of evolution of the Greenland ice sheet.Edwards et al. (2014) use six numerical ice-sheet models to evaluate three types of modeling uncertainties: climate model input, ice-sheet model choice, and the interaction of the two systems in terms of the surface mass balance-elevation feedback.While some common features from these papers can be extracted, some divergence in the results seems to be unavoidable.
SeaRISE (Sea-level Response to Ice Sheet Evolution) is a multi-model community effort to investigate the likely range of evolution of the Greenland and Antarctic ice sheets over the next few hundred years (Bindschadler et al., 2013).A total of eight models participated for the Greenland experiments (Nowicki et al., 2013).A series of century-scale sensitivity experiments to prescribed changes in surface climate, sub-ice-shelf melting, and basal sliding were performed.The results exhibit a large range in projected changes for the icesheet volume: projected Greenland ice-sheet contributions to the global sea level for the future-climate experiment under the A1B scenario range from 5.4 to 38.7 cm at 500 years from the present day.The projected ranges are larger for experiments where future-climate scenarios are amplified by a factor of 2, ranging from 8.5 to 142.6 cm.One of the objectives of the SeaRISE project is to show the possible range of uncertainties in the ice-sheet projection of current ice-sheet models, because no single model can be identified to be the best in every aspect (Bindschadler et al., 2013).The approach of the SeaRISE project is rather unrestricted: some aspects in the experiment protocol are standardized, while many others are left to the individual participants.The former includes part of boundary conditions of the ice-sheet model, such as the present-day surface temperature, surface accumulation and bedrock topography.Scenarios for future surface-climate changes such as a 100-year time series of surface temperature, precipitation, and surface melting are provided.The latter includes structural differences in ice-sheet models, such as model numerics or approximation level, and the treatment of some boundary conditions such as the surface mass balance scheme.Bindschadler et al. (2013) identified differences in the methods to compute the surface mass balance among the par-ticipants as the primary source of the dispersion in the results of future-climate experiments on the Greenland ice sheet.Nowicki et al. (2013) further concluded that variations in the initial ice volume, and thus the initialization of the ice-sheet topography, are another source of uncertainty.However, detailed quantitative evaluation of the reasons for the dispersion was beyond the scope of the two papers.The effects of some of the characteristics have already been argued in previous studies.Greve and Herzfeld (2013) compared 500-year future-climate experiments with three different grid spacings of 20, 10, and 5 km and concluded that the sensitivities in the simulated ice-sheet volume are insignificant.
The present paper performs a "one-model" approach to evaluate the relative impact of the various factors on Greenland ice-sheet projections under the SeaRISE protocol.The numerical model used in this paper is IcIES (ice-sheet model for integrated Earth-system studies), which also participated in the SeaRISE experiments.There are at least 10 characteristics that differ among the ice-sheet models participating in SeaRISE (see Table 2 in Bindschadler et al., 2013), and most have two or more variations.Some concern numerical aspects, such as grid resolution and time stepping, and others are physical aspects, such as ice-flow mechanics and surface mass balance.
This paper does not intend to cover the sensitivities of all of the aspects.The initialization methods and the surface mass balance methods, proposed in Bindschadler et al. (2013) as possible sources of variation, and three more characteristics, the bedrock topography boundary conditions, the basal-sliding methods, and the treatment of advance/retreat in the ice-sheet margin, are chosen to investigate sensitivities in the present paper.Of the four different sets of future scenarios under the SeaRISE protocol, the surface-climate experiment (C1 to C3), the basal-sliding experiment (S1 to S3), the ice-shelf melting experiment (M1 to M3), and a combination experiment, the present paper only revisits the surfaceclimate experiment.
In the next section, the five model setup characteristics of focus in this study are introduced to demonstrate the variety of choices among SeaRISE models used in Bindschadler et al. (2013).In Sect.3, a model description of IcIES is given to outline the setup adopted in the submission of Bindschadler et al. (2013).In Sect.4, we describe the setup of the five characteristics to replace the IcIES standard configuration in the present experimental design.Results and discussion follow to understand and compare the possible sources of spread among the results of the SeaRISE models.

Bedrock topography
SeaRISE provides several different versions of the present Greenland ice-sheet topography (available at http:// websrv.cs.umt.edu/isis/index.php/Present_Day_Greenland)."Greenland Developmental Data Set" (hereafter referred to as dev1.2) includes a Jakobshavn trough in the bedrock and bathymetry topography of Bamber et al. (2001, the second last version in the protocol).For the latest protocol, the bedrock topography including a new compilation of the subglacial troughs over Jakobshavn Isbrae, Helheim, Kangerlussuaq, and Petermann glaciers following Herzfeld et al. (2012) is proposed (hereafter referred to as JHKP).Although the differences between these data sets are localized, significant differences in the simulated global features are possible.Herzfeld et al. (2012) presented a significant difference in the present-day simulated topography and velocity field by using the JHKP data set and an older data set without inclusion of the above four glacier troughs (corresponding to a version before dev1.2 in SeaRISE).In addition, significant differences were found in the response of the Greenland ice sheet to doubled-sliding experiments over 500 years (i.e., equivalent to the S3 experiment in SeaRISE).

Basal-sliding formulation
The available methods to compute basal sliding have several degrees of freedom.One method applies a Heaviside function at the pressure-melting point of the basal temperature, i.e., the basal sliding is set to 0 when below the pressure-melting point.Others apply a smooth sliding transition around the pressure-melting point (Hindmarsh and Le Meur, 2001), i.e., the basal sliding gradually becomes close to 0 below the pressure-melting point, partly for numerical stability and partly for physical reasons to introduce subgrid-scale variation of the basal sliding.Some models in SeaRISE explicitly document such a smooth transition to implement melting at sub-melting point temperatures.

Initialization method
Obviously, the accuracy of the simulated present-day ice sheet is crucial for future projections.It is possible that small errors in the simulated present-day state may affect the short-term projections (Arthern and Gudmundsson, 2010;Yan et al., 2013).In addition, since the climate depends on the surface topography and ice extent, present-day climate forcing computed in the simulation may already have some bias.This bias occurs both for simulations with ice-sheet models coupled to sophisticated climate models, as well as in simulations using simple climate parameterizations.Some previous studies compute surface temperature by a combination of a reference field obtained from observation-based studies and their perturbation via the lapse rate and changes in surface topography relative to the present-day observed surface topography.This implies that the computed surface temperature field in the model is identical to the observation when the modeled surface topography is the same as the observation.
The choice of initialization method was left to participants in SeaRISE, and three different techniques were applied by the SeaRISE/Greenland participants.One method is called initialization by "tuning" in Bindschadler et al. (2013), which may be better termed "inversion" or "optimization".This method inverts given data fields, e.g., basal friction coefficients, to adjust present-day observation fields, e.g., surface velocity.Internal temperature fields are usually assumed to be in a steady state with computed velocity fields under the present-day conditions.The second and third methods are called initialization by "spinning up", wherein the model is run with the input of climate history of glacial/interglacial cycles, e.g., derived from the GRIP ice-core record.Although in principle these two initialization methods are not mutually exclusive (e.g., Edwards et al., 2014), the choice of the SeaRISE participants are either of the two.The first of these, hereafter referred to as "free spinning up", allows the icesheet topography to evolve freely under a prescribed climate history.A major disadvantage of free spinning up is that the present-day simulated topography often deviates from reality.The other initialization method is referred to as "fixed topography spinning up", in which the ice-sheet topography is fixed through the spin-up phase at a slightly smoothed measured present-day topography while ice-sheet temperatures freely evolve.The fixed topography spinning up is a hybrid of the two techniques in which the initial topography can be very close to the present-day observation while ice-sheet internal states include the influence of the long-term climate history.One major drawback is that the flow and temperature fields in the initial state are not in equilibrium (Goelzer et al., 2013), which leads to an artificial drift to restore the equilibrium after allowing evolution of the topography.
A number of studies have focused on the initialization methods and their impact on the simulation of the Greenland ice sheet.Rogozhina et al. (2011) compare the simulated present-day Greenland ice sheet obtained by several initialization methods including free transient spinning up.Pollard and DeConto (2012) presented a general and simple method to deduce spatial distribution of basal-sliding coefficients to reduce the errors in simulated surface topography that can be applied to any type of ice-sheet model.Morlighem et al. (2011) presented another approach in which uncertainties in the bedrock topography were also taken into account in the inversion method.Goelzer et al. (2013) presented a series of Greenland ice-sheet simulations with yet another hybrid technique to incorporate the influence of long-term climate history and obtain an initial ice-sheet topography close to the present-day conditions, by adjusting ice-temperature profiles and synthetic corrections over the surface mass balance.They concluded that the uncertainty arising from the surface mass balance methods and scenarios have a larger impact on the sensitivity of short-term projection of the Greenland ice sheet than those from the initialization methods, but the ex-F.Saito et al.: SeaRISE revisited perimental settings were not the same as the SeaRISE experiment.Aðalgeirsdóttir et al. (2014) presented a series of Greenland ice-sheet simulations using free transient spinning up as well as a flux-corrected initialization method, in which the surface mass balance during the initialization is modified such that the simulated present-day topography is close to observations (similar in principle to Goelzer et al., 2013 above).They concluded that the initialization methods are an important source of uncertainty.Yan et al. (2013) compared the evolution of the Greenland ice sheet to future-climate scenarios between two spin-up methods: free spinning up under transient and steady-state climate forcing.Both the simulated present-day ice-sheet topography and the simulated surface mass balance were different, thus the impact of the difference in the initialization method includes all of these components.Seroussi et al. (2013) found that the ice-sheet model is far more sensitive to changes in external forcing than its initial temperature for a 100-year-scale experiment, while futurescenario experiments from different initial conditions were not discussed.
So far, the influence of the fixed topography spinning up has not been discussed except for Goelzer et al. (2013), who showed an example in their configuration.This is a main target of the present paper.In addition, although Nowicki et al. (2013) concluded that variation of the initial ice volume may be a source of the uncertainties in the SeaRISE results, the influence of different choices for the initialization methods were not qualitatively evaluated.This paper extends their discussion and shows the relative significance to the short-term projection among other possible methods.

Treatment of advance of the ice-sheet margin
Precise simulation of the ice-sheet margins (ice-sheet extent) is a challenging issue.When ice-sheet topography and extent are allowed to evolve freely during future-warming experiments, it is possible to obtain sudden jumps in the position of the ice-sheet margin over many regions.Such changes reflect a strong flux imbalance near the margin in the simulated present-day state.Although detailed numerical implementation is not shown in Bindschadler et al. (2013), some participants in SeaRISE describe their methods as either fixing the ice-sheet margin (calving front) or limiting its advance (i.e., only retreat is allowed).While this is not necessarily true in reality because speeding up at the margin may result in advance before increased melting, some models just use this assumption.Previous studies have not demonstrated its influence on the sensitivity of the results, and so this issue is explored here.

Surface mass balance
The four aspects described above involve the structural (internal) rather than external (input) configuration.The method to compute the surface mass balance to drive ice-sheet mod-els instead affects the external configuration, and uncertainty relating to this aspect has a more direct impact on the simulated response of the Greenland ice sheet to climate warming.There has been a wide range of methods used to compute surface melting and/or surface mass balance in previous works including SeaRISE.
The method used to compute surface mass balance was left to individual choice in the SeaRISE project, which provided the future scenarios of precipitation, surface temperature, and surface melting, but whether or not to adopt unique parameterization of surface melting using the scenarios of precipitation and surface temperature was left to individual models.
Most participants adopted some form of the "positive degree-day" (PDD) scheme (Reeh, 1991) to compute surface melting.Even models using the PDD scheme, however, can vary in one or more parameters used in the scheme, e.g., the conversion coefficients from simulated degree-day to melting, the standard deviation (SD) of short-term statistical air temperature fluctuation (Gaussian noise added to parameterized monthly data), and so on.Previous studies showed how variation in PDD schemes and their coefficients can influence present-day and future simulation of the Greenland ice sheet (e.g., Stone et al., 2010).Bindschadler et al. (2013) argued that the variation of the surface mass balance method is the likely primary source of the dispersion in the results of future-climate experiments, although this assertion has not been quantitatively evaluated.This paper will demonstrate the relative significance of the surface mass balance method on the short-term projection compared to other model settings.

Aspects not tested in the present paper
The five aspects mentioned above are a subset of possible sources of the spread.As summarized in Table 2 in Bindschadler et al. (2013), there are at least 10 characteristics with different implementations among the participating ice-sheet models of SeaRISE.The remaining aspects are the numerical method (finite difference or finite element), the horizontal and vertical grid resolutions, the time step, the ice-flow mechanics (the shallow ice approximation, full Stokes), and the basal hydrology computation.The dependence on the stress in the basal-sliding formulation is also different among the models in addition to the sub-melt sliding formulation.It is possible for there to be other differences in model aspects not in the table, such as the ice enhancement factor, individual numerical schemes, and so on.Exploration of these remaining aspects was partly performed in previous studies (e.g., Greve and Herzfeld, 2013), and others are left for future studies.

Model description
The time-dependent, three-dimensional, thermodynamically coupled model used in this paper as well as in the SeaRISE project, called IcIES, is described in Saito and Abe-Ouchi (2005), Greve et al. (2011), andBindschadler et al. (2013).The model computes the evolution of the ice thickness, bedrock elevation, and ice temperature under a history of climate forcing, given in terms of surface mass balance and surface temperature, which may depend on the computed ice-sheet topography.The model parameters are the same as those described in Greve et al. (2011).In the present paper, the model domain spans 1500 km × 2800 km, with (151 × 281 grid points) corresponding to a horizontal resolution of 10 km.
The evolution of surface elevation is determined by the continuity equation for the local ice thickness with a history of the surface mass balance field.The temperature distribution is calculated by a thermodynamic equation with the surface temperature and geothermal heat flux given at the surface and base, respectively.Changes in the bedrock elevation are calculated by a linear model expressing local isostatic rebound with a prescribed time constant.
The shallow ice approximation is applied (Hutter, 1983) using Glen's flow law with an exponent of n = 3 (Paterson, 1994) for the velocity computation.The horizontal velocity vector v H is calculated for the given surface elevation h and bedrock topography b: where g is the acceleration of gravity, ρ I is the density of ice, and v B is the basal-sliding velocity.The rate factor A(T ), through which the velocity and temperature fields are coupled, follows Paterson (1994) and Huybrechts (1992).The formulation in Paterson (1994) is different from the one in Cuffey and Paterson (2010).We use the former in this study for a historical reason, to maintain consistency with the past numerical studies using IcIES including the submission to SeaRISE.Another reason is that the focus of this paper is on sensitivity to different external and technical configurations but not to "ice-flow" physics.The enhancement factor E in Eq. (1), which controls the softness of ice, implicitly reflects the effect of impurity and/or anisotropy of ice.It is used as a tuning parameter to improve the agreement between the measured and modeled surface topography.In the present paper, the constant value E = 3 is adopted in all experiments except where explicitly described.The basal-sliding velocity v B is computed using the Weertman sliding law, with an allowance for sub-melt sliding fol-lowing Hindmarsh and Le Meur (2001): where τ B , N B , and T B are the basal shear stress, basal normal stress, and basal temperature relative to the pressure-melting point, respectively.The function f (T B ) controls the occurrence of basal sliding (see Sect. 4).Following Huybrechts and de Wolde (1999), the exponents p, q and the coefficients C B are set to 3, 1, and 1.8 × 10 −10 N −3 yr −1 m 8 , respectively, for the standard configuration (v1, see Sect. 4).
The computation of the surface temperature follows Fausto et al. (2009): it linearly depends on the surface elevation, longitude, and latitude and an anomaly term that describes the paleoclimate temperature history or futureclimate-warming scenarios.The annual and summer mean surface temperatures are parameterized separately, and monthly mean temperatures are estimated from interpolation of the two fields using a sinusoidal function.The surface mass balance field is computed as the sum of the accumulation and ablation fields.The present-day mean annual precipitation (Ettema et al., 2009) is modified by a temperaturedependent function following Huybrechts et al. (2002).Conversion from the precipitation to the accumulation rate is computed statistically as in Huybrechts and de Wolde (1999), which is a function of the mean monthly temperature.Ablation (surface melting) is computed using the PDD method of Reeh (1991), which relates ablation to both air temperature and snow accumulation.The amount of melting is computed as the product of the number of positive degree days and PDD factors obtained by observations.It considers the possibility for melting even when the average daily temperature is below the freezing point, different melt rates for melt of snow and ice due to the albedo difference (Braithwaite and Olesen, 1989), and the production of superimposed ice and warming caused by the phase change.This method is adopted in most numerical studies with ice-sheet models (Ritz et al., 1997;Greve, 2000;Huybrechts et al., 2002).Four parameters control the surface melting in the PDD scheme in IcIES, the PDD factor for ice melt, β ice , the PDD factor for snow melt, β snow , the SD of short-term air temperature fluctuation, σ , and the saturation factor for the formation of superimposed ice, P max .The selection of the values of these parameters is described later.
All experiments in the present paper were performed with a newer revision of IcIES than that used for the SeaRISE project.To obtain stable simulations over all the experiments with a unique method, some modifications to the numerical representation were implemented.The physics and the mathematical formulation of the physics were not changed.The difference in the volumes of the simulated Greenland ice sheet for identical configurations varied at most by 0.3 %, which does not affect the conclusions of the present paper.Therefore, although the model itself is slightly modified, the  2013).The "A1B climate change" scenario, C1, over 500 years is now available, in which the first 100 years are obtained from climate model results, and the climate state of the final 400 years is kept constant at the 100-year climate.Two more "enhanced climate change" scenarios, C2 and C3, are defined in which the climate change of C1 with respect to the present day is amplified by factors of 1.5 and 2.0, respectively.In addition, a "constant present-day climate" scenario, C0, is defined for reference experiments.
One of the major uncertainties relating to ice-sheet dynamics stems from the basal-sliding processes because they are poorly understood due to the difficulties in direct observation (e.g., Nowicki et al., 2013).Often, the parameters relating to basal sliding are tuned to match present-day observed features such as ice-sheet topography and/or the surface velocity.Some models adopt spatially homogeneous parameters (e.g., Robinson et al., 2011), while others apply an inversion technique to compute spatially variable parameters (e.g., Seroussi et al., 2013).In the present paper, the impact of homogeneous changes in the basal-sliding coefficients are shown to interpret the results.Generally, the simulated icesheet thickness is too large, especially near the margin (Nowicki et al., 2013), and larger basal-sliding coefficients are required to reduce the error.In this paper, the cases of uniform doubled (v2) and quadrupled (v4) basal-sliding coefficients are examined.All of the experiments are repeated using these coefficients throughout the simulation.It is worth mentioning that the enhanced sliding experiments in the present paper differ from the "basal-sliding experiment" (e.g., S1) presented in SeaRISE.The former keeps the same value for the sliding coefficients over both the spin-up and the future, while the latter changes the coefficients for the future experiment only.
In addition to the experiment using uniform basal-sliding coefficients, some experiments are performed with a nonuniform basal-sliding coefficient field (case vm).Since the case vm partly relates to the initialization methods, it is described in detail in Sect.4.7.
Table 1 summarizes the sensitivity experiments in the present paper.The original IcIES submission, which is referred to as configuration O, adopts the following methods for the five characteristics: -the "Greenland Developmental Data Set" (dev1.2) for the bedrock topography; -basal sliding following the Weertman law without allowance of sub-melt sliding; -the "free" spinning-up method to initialize the presentday ice-sheet topography; -"free" advance/retreat of ice-sheet margin in response to the climate boundary condition; -the positive degree-day method for surface melting using a modification of Tarasov and Peltier (2002), where the standard deviation of the short-term statistical air temperature fluctuations to compute daily temperatures from monthly temperatures is set as 5.5 K in the IcIES original submission, which is slightly larger than the value of 5.2 K in Tarasov and Peltier (2002), Each of the five characteristics has two or more choices among the SeaRISE models.In the present paper, one method for each characteristic is chosen, to demonstrate structural uncertainties on the projection of Greenland ice sheet.A series of four experiments, A-B-D-E, is the sequence of one-by-one replacement in four methods: bedrock, sub-melt sliding, initialization, and margin advance, starting from the original configuration O. Experiment F is an additional sensitivity experiment which focuses on the impact of "fixed-topography" transient spin-up (will be described in Sect.4.6).D, E, and F have variation of fixed-topography steady-state spin-up, named as D s , E s , and F s , respectively.Configurations B , D , E , F , D s , E s , and F s were performed with an additional replacement in the surface mass balance computation.Finally, configuration E s is an additional experiment shown in the Appendix.The details of these replacements are described below.

Bedrock topography (A)
The bedrock topography dev1.2 used in the original configuration O is replaced by the JHKP data set in experiment A. All the procedures are then repeated with the new bedrock data.

Basal-sliding formulation (B)
The original IcIES submission adopts a Heaviside function at the pressure-melting point for the occurrence of basal slid- The Cryosphere, 10, 43-63, 2016 Table 1.Summary of numerical experiments in this paper.The bedrock column denotes the sources of bedrock topography as a boundary condition (see main text for interpretation of symbols).The column "sub-melt" denotes whether or not the sub-melt basal-sliding occurrence based on Eq. ( 3) is implemented.The initialization columns denote climate forcing used for initializing the ice-sheet topography, where "125 ky tr" stands for 125 kyr transient forcing based on ice-core record.Thickness columns denote how the ice thickness is computed during initialization phase, where "free" means that ice thickness is allowed to evolve freely, "fixed (obs)" means that ice thickness is kept fixed as the present-day observation through the initialization phase artificially, and "fixed (B 0 ka)" means that ice thickness kept fixed as the simulated topography at 0 ka obtained by experiments with configuration B. The margin column denotes whether the ice margin is allowed to advance freely (free) or limited to the initial condition (no advance) during future-climate experiments.The differences from the previous row are shown in bold.All the configurations are repeated with all the combinations of the basal-sliding coefficients (cases v1, v2, and v4) and the future-climate scenarios (C0, C1, C2, and C3).The experiment with suffix "s" (e.g., D s ) indicates steady-state initialization under the present-day conditions, which is denoted as "0 ka st" in the initialization column.The experiments denoted with prime (like B ) means switching the method to compute surface melting from PDD following Tarasov and Peltier (2002) (denoted as "T" in the surface melt column) to PDD of Huybrechts and de Wolde (1999) (denoted as "H").The table also includes an additional experiment E s shown in the Appendix, which uses another method of surface mass balance (indicated by symbol "S").Details are described in the Appendix.ing, which means that sub-melt sliding is prevented.It corresponds to the use of a binary operator with f = 1 if the bottom temperature is at the pressure-melting point and f = 0 otherwise, see Eq. ( 2).The Heaviside-function switch in A is replaced by an exponential function of the basal temperature to allow the occurrence of sub-melt sliding following Greve (2005) and Greve et al. (2011): where T B is the basal temperature relative to the pressuremelting point in • C, the parameter γ = 1 is used in the present paper.Formulation and/or the parameters of submelt sliding inclusion may vary among the SeaRISE models, but in the present paper the formulation above is chosen for demonstration of sub-melt sliding.The cases of uniform doubled (v2) and quadrupled (v4) basal-sliding coefficients, which are tests for model tuning and differs from the basalsliding experiment presented in SeaRISE, are examined with the allowance of sub-melt sliding occurrence.

Initialization method (D and D s )
For the original submission, IcIES used the "free spinningup" method.The background temperature history is based on the oxygen isotope record of the GRIP ice core (Dansgaard et al., 1993;Johnsen et al., 1997), which is provided by SeaRISE as a time series of temperature from 125 ka to the present.At the beginning, a steady-state simulation is performed under the climate field at 125 ka, and from this steady-state condition, the ice thickness and temperature and the bedrock topography are allowed to evolve freely until 0 ka.Two other methods are tested in the present paper: the "fixed topography transient spinning up" and the "fixed topography steady-state spinning up".The first is identical to the free spinning up except that the ice sheet and bedrock topographies are fixed to the present-day state and only the temperature can evolve.Thus the ice-sheet topography used as the initial condition for the future-climate experiment is In the fixed topography steady-state spinning-up method, a steady-state simulation is performed under present-day climate and topography fields with evolving temperature.This initialization method mimics the tuning method, where the ice-sheet topography is very close to present-day observations, while the influence of the long-term climate history is excluded.This initialization requires an inversion of, e.g., the coefficients of basal velocity, which is mimicked by different basal-sliding coefficients.In addition, experiments with spatially non-uniform basal-sliding coefficients are performed as another mimic of the tuning method, which will be described in the later section.

Treatment of advance of the ice-sheet margin (E and E s )
Both advance and retreat of the ice-sheet margin are freely allowed in the original configuration of IcIES.The thickness can be non-zero over the entire model domain during one step in the numerical time integration, but those grids that match a floating condition are immediately cut off.The configuration E (E s ) is equivalent to D (D s ) except that only retreat in the ice-sheet margin is allowed after the presentday simulation.There are some possibilities of how to implement the prohibition of ice-sheet advance numerically.In the present paper, the solution of the ice thickness beyond the present ice-sheet area is set to 0 during the time integration.1

Surface mass balance (e.g., B )
In the original IcIES submission, the PDD factor for ice melt is a cubic function of the local mean July surface temperature with a range between a minimum of 8.3 mm and a maximum of 17.22 mm ice equivalent per day per degree (Tarasov and Peltier, 2002).The factor for snow melt is a linear function of local mean July surface temperature with the range between a minimum of 2.65 mm and a maximum of 4.3 mm ice equivalent per day per degree.Some models in SeaRISE use a PDD scheme with different parameters, and others used other simplified schemes (Bindschadler et al., 2013).One variation of the PDD scheme is chosen in the present paper.Some models adopt constant (temperature-independent) coefficients, such as 3 and 8 mm ice equivalent per day per degree for snow and ice, respectively, following Huybrechts and de Wolde (1999).In the present paper, a variant of this PDD approach has been chosen with slightly larger SD of short-term statistical air temperature fluctuation as 5.5 K.

Impact of fixed-topography transient spin-up (e.g., F)
One aspect remaining to be discussed is the impact of non-equilibrium internal states originating from the fixedtopography transient spin-up.Since there is a feedback between climate and ice-sheet topography, the difference between free spin-up and fixed topography spin-up includes both the effect of internal temperature and of the initial topography.One way to minimize the initial discrepancy and to separate the impact of non-equilibrium internal states is to perform a free spinning-up simulation that ends with the same topography at the present-day.The impact of the internal non-equilibrium state is evaluated as follows: experiment F (or F ) is initialized by fixed-topography transient spin-up with the topography fixed through 125 kyr as the final state of the spin-up phase obtained by configuration B (or B ) instead of the present-day topography as D (or D ).Thus the difference between experiments B (B ) and F (F ) only stems from the internal thermal state due to the initialization methods, both having an identical initial topography.
To evaluate the impact of "no memory" of the transient past climate, further fixed-topography steady-state spin-up experiments are performed (F s and F s ).Instead of the topography being fixed at the present-day observation, as for configuration D s and D s , it is fixed at the final topography of the spin-up phases of experiments B and B , respectively.
The series of experiments outlined in this section is summarized as follows: -B = free topography + transient temperature,

Impact of non-uniform basal-sliding coefficient (vm, e1:vm)
Another aspect remaining to be discussed is the impact of initialization by tuning or inversion.There are three models in the SeaRISE Greenland experiment that use a form of inversion, and these all differ not only in the method and parameter tuned but also in other aspects such as basal-sliding formulation and surface mass balance.The results of the three models have already a dispersion as shown in Bindschadler et al. ( 2013), Fig. 3, due to partial or all combination of difference among the models.An inversion experiment could be performed using the same method as these three models or another method such as Pollard and DeConto (2012).
Generally, the inversion depends on the boundary conditions The Cryosphere, 10, 43-63, 2016 www.the-cryosphere.net/10/43/2016/such as surface mass balance as well as the ice-flow characteristics in individual models, which are different among the SeaRISE inversion models.Therefore, even if an experiment following one or all of the inversion methods in the SeaRISE is conducted, the many degrees of freedom mean that the results may not explain the dispersion in the SeaRISE results.However, "potential" explanations of the impact of an inversion are worthy of exploration.The essential difference between the inversion models and the others is the application of non-uniform parameter fields such as basalsliding coefficients (with a certain assumption for other fields such as ice temperature and enhancement factors).In order to demonstrate a potential impact of the inversion, we repeat some experiment configurations using a prescribed field of non-uniform basal-sliding coefficients kept constant throughout the simulation.Pollard and DeConto (2012) presented a general method to deduce spatial distribution of basal-sliding coefficients to reduce the errors in the simulated surface topography.In this method, the evolution of ice-sheet topography and temperature are computed for a prescribed surface mass balance, periodically adjusting the basal coefficient at each grid point according to the error of local surface elevation compared to the present-day observation.In the present paper, the method is applied with modification for the minimum and maximum limits of the basal-sliding coefficient, which are chosen as 10 −8 × C B,v2 and 10 4 × C B,v2 , respectively (see Eq. 2), after some trials.The same boundary condition as B is applied for this computation.As described in the model section, the standard enhancement factor in the present study is 3, but using this value never gives a reasonable coefficient field: even when no sliding is allowed as the lower limit of the coefficients, the interior part of the ice sheet is still lower by more than 400 m compared to the present-day observation.Following Pollard and DeConto (2012), we modify the enhancement factor to a smaller value.The enhancement factor set to 1 for this inversion procedure (configuration e1).The configuration vm is the run with the "inverted" non-uniform basalsliding coefficient fields that are computed under B :C0:e1.In addition, for comparison purposes, a subset of the uniform basal-sliding coefficient runs is repeated with the small enhancement factor E = 1.
The non-uniform basal-sliding coefficient cases vm are conducted as a variation of the uniform basal-sliding cases.All of the experiments are repeated using these coefficients fixed throughout the simulation, both over the spin-up and the future.

Results
Table 2 summarizes the simulated ice-sheet volumes at the end of the initialization phase (or at the beginning of futureclimate scenario experiments) compared to the present-day observations, and the root mean square of residuals of the thickness.The results of experiments except for O will be described later in this section.Under configuration O, the overestimation of the ice-sheet volume is within +6 % and with increased basal-sliding coefficient v4 within 0.5 % of the present-day observations.The good match of the simulated volumes can be explained by an overestimation around the ice-sheet margin and an underestimation over the interior regions (e.g., Bindschadler et al., 2013;Yan et al., 2013).The difference in surface elevation relative to the present-day observation is included in the Supplement (Fig. S1).Since generally the simulated thickness over the interior region is larger while that over the margin is smaller than the presentday observation, the root mean square of the residuals is as large as 500 m even when the total volume is close to the observation, which is a common feature among model studies, in particular using a spatially uniform basal velocity coefficient (e.g., Nowicki et al., 2013).The results of configuration O show volume losses of 34.1, 72.1, and 142.8 cm sea-level equivalent at 500 years under climate scenarios C1, C2, and C3, respectively.Standard basal-sliding cases v1 under all future-climate scenarios are within the range of original SeaRISE results.Simulated responses become larger with enhanced basal-sliding coefficient, and some cases are still within the original range of results, while some are above the range, for example, the simulated VAF response of C3:v4 is 17 cm more than the upper boundary of the original range.

C1:v1
In the following sections, the effects of replacement of the five model aspects are described in turn.The fractional changes of the effects of this series of experiments are summarized in Table S1 and Figs.S4-S6 in the Supplement.

Bedrock topography
Configuration A is equivalent to O (SeaRISE/IcIES configuration), except that the bedrock topography dev1.2 is replaced by the JHKP topography.Simulated VAF responses are affected by replacing the bed topography of a few regions but are less than +2.2 cm under all the combinations of climate and sliding coefficients (Fig. 2).3) at 500 years.The circles in the gray bars indicate the mean values of all the SeaRISE participants.

Basal-sliding formulation
Configuration B is equivalent to A (O with JHKP bedrock), except for the inclusion of sub-melt sliding during both the initialization and future scenario phases.Table 2 shows the simulated volumes at the end of the initialization with configuration B. The introduction of the sub-melt sliding results in a wider sliding area and therefore a smaller ice-sheet volume due to enhanced outward ice flow.The standard basalsliding coefficient case v1 shows ice-sheet volumes close to present-day observations (1.7 % overestimation).Similar to other configurations, such as O and A, the increase in the basal-sliding coefficient leads to smaller present-day icesheet volumes.In the case v4 with 4× coefficients, the resulting present-day ice volume underestimates observations by more than 10 %.For O and A, quadrupling the basalsliding coefficient varies the volume by around 5 % of observed, while it does the same for B by more than 12 %. Figure 3a-c show simulated present-day ice-sheet topography obtained by B:v1 to B:v4 cases, respectively.The Supplement includes a figure showing the difference in the surface elevation relative to the present-day observation (Fig. S2a-c).The interior part of the ice sheet becomes lower with an increasing basal-sliding coefficient.In addition, the ice-covered area around the northwest region is much reduced with a higher basal-sliding coefficient in particular in the B:v4 case, which partly contributes to the overall underestimation in the volume.
The sub-melt sliding treatment affects the VAF response more for larger sliding coefficients as shown in Fig. 2 (comparing B with A).For v4, the absolute increases in the VAF from A to B are similar between C2 and C3 scenarios (+26.4 and +24.9 cm, respectively), and the ratios of the increases in the VAF to the corresponding total VAF become smaller from the lower climate scenario C1 to the higher C3.Also, the case of v4 has proportionally less difference in the higher climate scenarios when comparing the change between B and A. The C1:v1 case results in a loss of 36.5 cm at 500 years (which is about 1 cm more than in case A).The largest difference between B and A is +26.4 cm for the C2:v4 case.

Initialization method
Configuration D is equivalent to B (free transient spin-up), except that the ice-sheet initial condition is obtained by a fixed topography spin-up given by the present-day observation.Because of the inconsistency in the internal temperature due to fixed topography spin-up, larger drifts are shown even under the constant-climate scenario run (C0) compared with those of the free spin-up configuration (B).Similar to Bindschadler et al. (2013), no configuration matches the observed rate of present-day volume change.These drifts are subtracted from the results under future-climate runs (C1 to C3), in order to isolate the response to the forcing alone.The simulated response of the VAF is 26.0 cm for D under the C1:v1 case, therefore it has −10.5 cm impact relative to B. This more than cancels the impacts of the treatment of bedrock topography and sub-melt sliding (Fig. 2).Under the C2 and C3 cases, VAF are 52.6 and 111.6 cm, which shows −24.5 and −39.3 cm impact, respectively.Thus, the impact of whether the topography is free or fixed within the spin-up to observed reaches around one-third of the range of the original SeaRISE experiments.Especially for under larger basal-sliding coefficients, cases v2 and v4, VAF are significantly reduced due to the different spin-up condition whether free or fixed, which are large enough to cancel the effect of including sub-melt sliding.Simulated responses in VAF are reduced to 50 % or less from B to D.
Figure 4 shows the changes in VAF relative to that under the constant-climate scenario C0 obtained by experiments B (free topography spin-up), F (fixed topography spin-up as for B), and D (fixed topography spin-up as for the observation), over all the combinations of climate scenarios and sliding coefficients.
Configuration F is equivalent to B (free transient spinup), except that the ice-sheet initial condition is obtained by a fixed topography spinning up as the final state of configuration B, which means that the initial topography for futureclimate runs are identical.Since internal thermal states are not in equilibrium under configuration F due to the artificial prohibition of topography evolution, the thermal conditions drift to restore the equilibrium during the future-climate run even under the constant-climate simulation.The effect of the non-equilibrium thermal state is not systemically larger from F to B. In the case of v1 basal sliding, F (fixed topography spin-up), shows a smaller response than B (free transient spin-up), under the C1 scenario, similar response under the C2, a larger response under C3, respectively, while in the case of v4 all F show a larger response under the three scenarios.Over all the combinations of climate and sliding coefficient examined in the present paper, the differences in the final states of VAF between B (free topography spin-up) and F (fixed topography spin-up as B) are smaller than the differences in VAF between B and D (fixed topography spin-up as observation).This means that the model sensitivity of the internal non-equilibrium thermal states is smaller than the sensitivity to the free or fixed topography options, when they are evaluated in terms of changes relative to the constant-climate experiment.The effect of the non-equilibrium thermal state is larger for larger VAF, because the elevation-ablation feedback amplifies the geometry changes.The maximum impact in the present paper is +14.5 cm sea-level equivalent for F under the C3:v4 case, which is 10.5 % of the variability of corresponding D cases.Configuration F s is equivalent to F (fixed topography transient spin-up as B) except that the initial condition of the ice sheet is obtained by a fixed topography steady-state spinning up as the final state of configuration B. All the experiments show almost identical sensitivity of VAF between steady-state and transient spin-up, in terms of relative changes in VAF to the corresponding constant-climate scenario cases.In other words, as long as the final topography is the same, it does not make much difference whether the spin-up used a transient climate or steady state.Therefore, if an initial state with free spinning-up methods ends at the observed topography, the time evolution of VAF is expected to come close to the one obtained by fixed spinning-up meth-ods, both under transient climate scenarios and under the constant present-day climate scenario imposed for the first 500 years.

Treatment of advance of the ice-sheet margin
The initialization phase of configuration E is identical to that of D (free margin, fixed-topography transient spin-up as observation), but advance of the ice-sheet margin is not allowed while retreat is freely allowed under future-climate runs.Prohibiting ice-margin advance has a smaller impact than the choice of spin up whether free or fixed (Fig. 2).The simulated response of VAF is 19.8 cm in experiment E and −6.2 cm relative to D under the C1:v1 case.Thus, under mild climate-warming scenarios like C1, the choice of spin-up whether free or fixed and the margin treatment has a larger effect on the response of Greenland ice sheet over 500 years compared with the effects of bedrock or sub-melt sliding.The impact of replacing the treatment of the margin is affected little by the choice of basal coefficients, but the larger basal coefficients tend to have slightly more impact from the replacement.This reflects the fact that higher velocity at the margin tends to result in more advance in the margin.Under higher climate scenarios such as C3, advance in the ice-sheet margin is not significant even in the freemargin experiments D, and thus less impacts are seen from the replacement of the margin treatment.

Surface mass balance
Figure 2 shows the simulated changes in VAF under all of the combinations of climate scenarios and basal-sliding coefficients by the series of experiment B (free topography spin-up), D (fixed topography spin-up as observation) and E (no advance in the margin).Surface mass balance is replaced from B (PDD of Tarasov and Peltier, 2002) to B (PDD of Huybrechts and de Wolde, 1999), and after that the same replacement sequences are followed as B to E (initialization and margin treatment).
Configuration B is equivalent to B, except that the surface melting parameterization of Tarasov and Peltier (2002), which was used in the IcIES original submission, is replaced by Huybrechts and de Wolde (1999), which was used by some of the SeaRISE participants.The future-climate runs C1 and C0 and the initializations are repeated using the new PDD methods.Table 2 shows the simulated initial volumes under the configuration of the B series and Fig. 2 shows the simulated changes in VAF under all of the combination of climate scenarios and basal-sliding coefficients by the series of experiment B , D , and E .
With the change of the surface mass balance method, the simulated present-day ice-sheet volumes become larger by about 4 %. Figure 3d-f show simulated present-day ice-sheet topographies obtained by experiments B :v1 to B :v4, respectively.The supplementary includes a figure showing the difference in the surface elevation relative to the present-day observation (Fig. S2d-f).The main difference between B and B is found in northwestern Greenland.The retreat of the icesheet margin over northwestern Greenland is not seen in the B cases (Fig. 3d-f).Changes over the interior region (around the summit) are small because the change in method primarily influences the ablation area near the ice-sheet margin.
Figure 2 shows a volume loss of 28.2 cm sea-level equivalent at 500 years for configuration B ; thus replacing the PDD methods in the C1:v1 case has an impact of ∼ −8.3 cm.This impact is slightly smaller than the impact of −10.5 cm by replacing the free vs. fixed topography methods from B (free topography spin-up) to D (fixed topography spin-up as observation).The smaller sensitivity partly stems from the overestimation in the present-day topography.Since the simulated initial volume is larger, less surface melting is expected because of the elevation-temperature feedback.Under stronger warming scenarios, the impact of the replacement of the surface melting method from B to B , which are −21.9 and −50.8 cm under C2 and C3, respectively, is similar or even larger than that of the free vs. fixed topography methods from B to D. Similar to the replacement of the free vs. fixed topography methods, the large impact due to different basal-sliding formulation is canceled by the replacement of the surface melting method and the results become closer among the three cases of basal-sliding coefficients under the same climate scenarios.Through all combinations of climate scenarios and basal-sliding coefficients, a significant influence in the simulated responses of VAF due to the different surface mass balance methods are shown.As shown in Fig. 2, the difference in the surface melting methods has similarly large influences on simulated responses as the free vs. fixed topography methods.
Similarly, configurations D (fixed topography spin-up as observation) and E (no advance in the margin) are equivalent to D and E, respectively, except for the surface melting parameterization.Under the lower future-climate scenario C1 (Fig. 2a), the influence of the replacement of surface mass parameterization is comparable to that of replacement of both the free vs. fixed topography methods or the treatment of icesheet margin (B vs. D; D vs. E).Under the higher futurescenario C3, the influence of the former becomes even larger than that of the latter.Simulated responses in VAF are reduced to around 60 % of those obtained using the original surface mass balance parameterization (B vs. B ; D vs. D ; E vs. E ) under C3 future-climate scenario.Similar results to the case of the other surface parameterization are obtained as shown in the comparison of F (fixed topography spin-up as B ), F s , and B (free topography spin-up).
A comparison between the results of F or F s (B plus different fixed topography) and B (B plus different surface mass balance) shows the relative influence of the internal inconsistency and the surface mass balance parameterization.Over all the combinations considered in the present paper, the impact of the internal non-equilibrium thermal state to the sim- ulated sensitivity of VAF is smaller than the impact of the difference in the surface melting methods for both a steadystate and transient spin up.

Non-uniform basal-sliding coefficient field
Figure 5 shows the difference in the simulated surface topography relative to the present-day observation and inverted basal-sliding coefficient field.The inversion procedure is performed using the surface mass balance method of Huybrechts and de Wolde (1999), with the ice enhancement factor E = 1, and prohibition of advance in the ice-sheet margin.The last constraint is somewhat arbitrary but is kept for simplicity.The inverted coefficients are smaller than the case v2 value by some orders of magnitude over the interior region, while larger around the margin.Although not perfect, the overestimation in surface elevation near the margin and the underestimation in the interior part are significantly reduced (see Fig. S2d-f in the Supplement for uniform basal-sliding coefficient cases).As mentioned, since the inverted field is a function of other aspects such as surface mass balance, a different distribution should be computed for each configuration.Since the experiment reported in this section is intended to demonstrate non-uniform basal-sliding coefficient fields, the same field is used through all the experiment in this section.Among the series of experiments, E s and E s (sub-melt sliding included; fixed topography steady-state spin-up as the present-day observation; no advance in the ice margin) are performed using the inverted field, with default enhancement factors E = 3 and E = 1. Figure 6 shows the simulated changes in VAF under all of the combinations of climate scenarios and basal-sliding coefficients by the series of experiments E s and E s .The case state in the present paper.Among the four coefficient cases v1-vm, under E s :e1 configuration, the non-uniform case vm shows the smallest changes in VAF under each climate scenario C1 to C3.Similar to the other configuration discussed above, simulated responses become larger with a uniformly enhanced basal-sliding coefficient (e1:v1 to e1:v4), while the case e1:vm is even smaller than the case e1:v1.As shown in Fig. 5b, most regions near margins have very large basal-sliding coefficients while most of the interior has a very small value (even smaller than the value of v1), which leads to larger and smaller changes in VAF over the margin and interior, respectively.Thus totally smaller change in VAF than v4 case is shown in the run with E s :e1 configuration.The same holds true for the configuration E s :e1 (E = 1, PDD following Tarasov and Peltier, 2002).With the default en-hancement factor E = 3, the relation among four coefficient cases, v1 to vm, varies with experiment: for example, the non-uniform case vm has a response between v2 and v4 under the E s case.For all the non-uniform basal-sliding coefficient cases, the changes in VAF at 500 years never exceed the changes obtained by the v4 (uniformly quadrupled) case and sometimes become slightly smaller than those obtained by the v1 case.

Discussion
The simulated response of the Greenland ice sheet is affected by the method options explored in the experiments presented in this paper and is explained partly by the difference in the initial state and partly by that in the initial drifts.
Replacement of the bedrock topography (O to A) as well as the sub-melt sliding treatment (A to B) leads to lower elevations in some regions and thus larger responses under futurewarming scenarios due to the elevation-ablation feedback.
Prohibiting ice-margin advance (D to E and so on) suppress ice-sheet growth under the constant future-climate scenario C0, and thus the relative responses under warming climate scenarios C1-C3 become smaller.The replacement of initialization method whether free or fixed to the observation (B to D) leads to smaller responses under future-warming scenarios, which cannot simply be explained by the difference in the initial volume for some particular cases.As shown in Table 2, the initial ice-sheet volumes of B:v1, B :v1, and B :v2 are larger than those with fixed topography spin-up experiments, while the simulated VAF response are larger.This is due to inconsistency in the internal temperature field due to fixed topography spin-up, which leads to larger drift under the constant-climate scenario to approach the steady more rapidly, in this case, to larger ice-sheet volume and leads to smaller response due to the elevation-ablation feedback.An inversion method following Pollard and DeConto (2012) is applied to compute an optimized basal coefficient field in order to emulate some of the SeaRISE models initialized with a tuning method.Since the inverted non-uniform field depends on all model properties such as bedrock topography, surface mass balance, the basal-sliding formulation, and ice-flow mechanics, it is not guaranteed that the result with the non-uniform field in the present paper explains the behavior of the SeaRISE models using a tuning method.However, since at least the SeaRISE models with free topography spin-up including IcIES have qualitatively similar results (overestimated thickness around the margin while underestimated in the interior), the computed basal-sliding coefficient field in the present paper may capture the general characteristics of the expected inverted basal-sliding coefficient field.The inverted field in the present paper generally shows larger values around the margin and smaller in the interior, which leads to a larger response around the margin and a smaller response in the interior, under the future-climate scenarios.The total response of Greenland ice sheet is determined how much both responses are compensated.At least throughout the configuration of the present paper, the simulated total responses in VAF do not significantly deviate from the uniform basal-sliding coefficient cases.
The four methods examined in the series of transient spinup, O-A-B-D-E with all cases of the basal-sliding coefficients v1 to v4 under all the future-climate scenarios C1 to C3 (bluish group in Fig. 2), are related to the ice flow but do not relate to the model inputs.Among these four aspects, the inclusion of sub-melt sliding enhances the ice-sheet response strongest (A to B), but using fixed-topography spinup cancels and even reduces this impact (B to D). Prohibition of ice-sheet advance is a secondary influence that can reduce the sensitivity (D to E).For the lower future-climate scenario case C1, the combination of all four aspects (Fig. 2a) affects the volume loss as much as 42 %, which leads to the response of 19.8 cm sea-level equivalent in experiment E. This value is very close to the average of SeaRISE participants (19.2 cm sea-level equivalent) presented in Bindschadler et al. (2013), regardless of the basal-sliding coefficient.For the higher future-climate scenario case C3 (Fig. 2c), the combination of all four aspects affects the volume loss by as much as 30 % of the total response, which is not enough to explain the large deviation of O from the average.The spread of the results due to different basal-sliding coefficients is similar between the C2 and C3 scenarios.Thus the source of spread in SeaRISE experiments can only partly be explained by variations in the experimental configuration of technical aspects of ice flow.The most influential of these is the specification of free or fixed geometry and slightly less, the treatment of the icesheet margin evolution.Using a non-uniform basal-sliding coefficient field and/or smaller enhancement factor has the potential to further reduce the volume loss (Fig. 6).Although significant changes in the volume loss are not shown using the inverted field in the present paper, it is still possible to have larger impacts on the changes using different basalsliding fields.
The uncertainty in the methods to compute surface melting can further influence the model sensitivity.Configuration E replaces all four technical aspects as well as the surface mass balance compared to the original configuration O. E results in a volume loss which is smaller than the average of the SeaRISE experiments for the C1 future-climate scenario.Even for the highest climate scenario, case C3, the volume response is slightly smaller than or close to the average of the SeaRISE experiments, regardless of the basal-sliding coefficient (Fig. 2c).Again, significant changes in the volume loss are not shown using the inverted field in the present paper (Fig. 6), but it still has the potential to explain the spread in the SeaRISE results.
In the series of the experiments in the present paper, the choices that have greatest effect on the simulated response are the method to compute the surface mass balance and the method to initialize the ice sheet, which have comparable ef-fects.This is consistent with the discussion of the possible reasons for spread in the SeaRISE results by Bindschadler et al. (2013) and Nowicki et al. (2013).The variation of the surface mass balance alone (B to B ) has some influence on the ice-sheet sensitivity, but not enough to completely cancel the large volume response obtained by the IcIES original configuration (i.e., configuration O with v1 basal sliding).The influence of the initialization methods (whether free or fixed topography) on the short-term ice-sheet sensitivity is comparable to the influence of uncertainties in the surface mass balance methods.Moreover, the influence of the artificial prohibition of the advance of ice-sheet margin, which is not discussed in the papers, is found to be secondary to the main two aspects but not negligible.
One drawback when using initialization methods, except for the "free" spin-up, is a drift due to inconsistency in simulated temperature fields.Comparison of the results between B (free topography spin-up) and F (fixed topography spin-up as B) or B and F , where the corresponding pairs have identical topography but different internal states, can show the influence of internal non-equilibrium thermal states.Over all the combinations in the present paper, the difference in the final states of VAF between B (B ) and F (F ) is smaller than the difference in that between B (B ) and D (D ).This implies that, at least in terms of changes relative to the constant-climate experiment, the influence of the internal non-equilibrium thermal states to the ice-sheet sensitivity is smaller than the influence of different initial states.The largest difference between B and F is found under the C3:v4 case, which shows a difference of +14.5 cm sea-level equivalent between the two different internal non-equilibrium thermal states.Since an expected counterpart of the D case, which has the identical topography to the present-day observation without artificial drifts, cannot be easily performed, an indirect evaluation is conducted as follows.This 14.5 cm effect is about 11 % of the simulated VAF response obtained by D C3:v4 case, and thus the effect of the internal non-equilibrium state is expected to remain minor relative to the total sensitivity.In other words, the initial topography has more effect on the future projection, in terms of relative to constant scenario runs, than the initial internal temperature field.Therefore, future-climate experiments initialized by fixed-topography spin-up are considered the preferable approaches for characteristic projections of the ice-sheet evolution by an icesheet model.In addition, in terms of changes relative to the constant-climate experiment, steady-state and transient spinup initializations show almost identical sensitivities during 500-year model runs.
Table 3 summarizes simulated changes in VAF of configurations B, F, D, and D relative to the corresponding constant future scenario experiments.Except for the lower sensitivity cases such as C1:v1 and C1:v2, the table shows that the effect of internal non-equilibrium states (B vs. F) is rather small compared to the effect of differences in surface mass balance methods (D vs. D ).Thus, the uncertainties due to surface mass balance must be another potential source of uncertainties in the simulated 500-year-scale future projections of the Greenland ice sheet rather than those due to ice-flow characteristics.
All the analysis in the present paper is performed using the anomaly relative to the result of the "constant" futureclimate experiment C0 ("experiment minus control"), following the discussion of the SeaRISE methods (Bindschadler et al., 2013;Nowicki et al., 2013).In other words, trends in the evolution of the ice-sheet volume at the present day, whether they are artificial or not or whether they are consistent with the present-day observation, are excluded from the discussion.Simulated trends vary among the configurations and range from −45 cm (E-v4) to +24 cm (D -v1) after 500 years among transient experiments.Steady-state experiments do not deviate much from the corresponding transient experiments.Simulated changes in VAF for some experiments are shown in the Supplement (Fig. S3).In reality, the trends arise as the result of long-term climate history.Since the trend is not necessarily zero, the actual future projection of the Greenland ice sheet should be evaluated as the sum of the trend and the anomalies.It is expected that such long-term memory has a smaller impact for the future changes in ice-sheet volume at least during the next 500 years, compared with the changes due to future surfaceclimate scenarios, because the results of transient spin-up (with long-term memory) and steady-state spin-up (without) show similar responses.In the present paper, only a part of the surface-climate experiments in SeaRISE has been revisited.The same procedures applied here can be followed for other series of experiments (e.g., basal-sliding experiments), which are left for the next study.
Therefore, although it cannot be confirmed, if a perfect spin-up (free evolution spin-up under transient climate ending with the present-day observed topography) could be obtained, then it can be expected that the VAF response of such an experiment would be close to that obtained using a fixedtopography spin-up with the present-day topography.Thus, a future-climate experiment initialized by fixed-topography spin-up (with the present-day topography) under either tran-sient climate history or steady-state climate can be considered a suitable approach for characteristic projection by an ice-sheet model in order to isolate the response to the prescribed climate scenario alone.While it cannot be fully confirmed, the analysis of the series of experiments in the present paper suggests that the large sensitivity of IcIES can be attributed to the use of a free topography during the spin-up, free evolving margin during the future experiment, and the difference in the surface melting parameterization.
Sensitivities due to different treatments of the margin advance need to be carefully interpreted, since marine boundaries are present for major Greenland outlet glaciers and thus marine-ice-sheet instabilities have been identified in numerical model studies (Nick et al., 2013).It is not mentioned explicitly, but most SeaRISE models determine the grounding line by a floating criterion (set H = 0 when the surface falls below flotation height) or fix the grounding line through time.Therefore marine-ice-sheet instabilities of the Greenland ice sheet are important in terms of future projection, but SeaRISE models do not have sufficient capability to represent grounding-line processes adequately.There are M1, M2, and M3 experiments in SeaRISE, which are called ice-shelf melting experiments.Since the SeaRISE Greenland models do not have explicit ice-shelf processes, the implementation of the "ice-shelf melting" varies greatly among the models, which is one of the reason why the spread of these results are very large (larger than C1, C2, and C3 spreads presented in this paper).Nowicki et al. (2013) state that "Thus, the current generation of Greenland whole ice-sheet models is not yet able to simulate the potential response to a warming ocean, and caution is needed when interpreting the SeaRISE response to this scenario, as the ensemble mean response likely underestimates the true potential response."For the same reason, the present paper focuses on atmospheric warming scenarios only, which means that the impact of margin retreat purely due to the surface mass balance is discussed.When marine-ice instability processes are included, the problem of margin advance/retreat may become more significant than those expected in the present paper.
Multiple combinations of changes in all of the aspects considered in the present paper (except for the bedrock topography) are tested in order to check for interactions between the uncertainties."One-at-a-time" effects are summarized in Table S1 and Figs.S4-S6.Although the detailed features vary among combinations, the general features in the results discussed in the present paper are also shown.Two aspects, free or fixed topography spin-up and the surface mass balance methods have larger influences than other aspects on the changes in VAF at 500 years over all the future-climate scenarios and the basal-sliding coefficients.Prohibition of ice-sheet advance has a large influence, in particular when the future-climate scenario is mild.The difference of the results by transient spin-up and those by steady-state spin-up are smaller among the other aspects throughout the combinations.Difference of the results by free transient spin-up and those by fixed transient spin-up (as free experiment) are always smaller than the difference of the former and those by fixed transient spin-up (as the observation).Except for C1:v4 and C2:v4 cases, inclusion of sub-melt sliding has less (or similar) influences than the two aspects of large impact (surface mass balance, fixed-topography spin-up).A large impact of sub-melt sliding inclusion is found when the surface mass balance follows Tarasov and Peltier (2002) and when the initial topography is the same as free topography spin-up (e.g., B, F).As described in Table 2 (see B and v4), simulated total volume at the present-day deviates most from the observation among the experiments, and the impact of switching off the sub-melt sliding inclusion (B to A) is as large as 10 %.Starting from such a small initial condition is considered to be a reason for the large impacts of changes in the sub-melt sliding formulation, through elevation-ablation feedback.
Since the ice-sheet models will become increasingly more complex, a one-model study such as the present paper cannot cover all possible variations among the existing models.It would be preferable that all participating models perform one common and highly controlled experiment that allows effective identification of the uncertainties due to specific variations in ice-sheet models.Such an experiment would not be an intercomparison for more realistic projections but rather an abstract test purely for model intercomparison purposes.The intercomparison experiments of the ice2sea projects (e.g., Edwards et al., 2014) mainly focus on model differences and therefore provide such controlled protocols except for the initialization methods.
The experiment in the present paper only covers a small part of the SeaRISE model choices, and thus there is insufficient comparison of the dependence of SeaRISE results on these choices.Nevertheless, it shows that structural and parametric uncertainties are just as important as initialization.In other words, it shows that if all the SeaRISE models repeated this study, the range of the results could widen beyond the current reported spread.Hence, it is important to systematically control and study uncertainties with such designed control experiments.
Here we propose a model intercomparison study to evaluate the uncertainties in modeled response that originate from modeled ice-flow characteristics such as ice-flow approximation level, basal-sliding formulation, and model resolution.The proposed experiment setup, which is referred to as the "benchmark" experiment, consists of a carefully controlled protocol to define the following characteristics.
-Initialization of the present-day condition using either -assimilation, -or fixed-topography spin-up; -preparation of "identical" model inputs in order to extract the influence by difference in ice-flow characteristics only, -(easier) not temperature but the spatial/temporal scenario of the surface mass balance with no topography or albedo feedback, -or provide an identical surface mass balance subroutine (not a scheme, in order to keep it really identical among the models) as well as scenarios, -with parameterization such as the PDD scheme, with a regional climate model, or with any methods used for ice-sheet future projections, as far as identical among the models; -performing two short-term future-climate experiments, a constant-climate experiment and a warming climate experiment, in order to subtract the influence of (artificial) drifts; -limitation of the advance of the ice-sheet margin to the present-day (initial) margin (Although the opposite approach is possible, this approach is much easier to implement in some models.Also in this case the treatment of boundary conditions over the ice-free grids does not need to be specified).
A demonstration of this type of experiment is presented in Appendix A. Since spinning-up methods are not specified, except for the ice-sheet topography, most types of ice-sheet models can easily perform this experiment, including computationally expensive full Stokes models, models using inversion techniques, and models using free evolution spinning up over a long climate history.This experiment configuration is a compromise to allow choice of initialization method by individual model but is, however, still proscribed enough to separate uncertainties and/or some feedbacks.The results of this benchmark would help to address the uncertainties obtained by other intercomparison experiments for more realistic projection with a large variety of model aspects like the SeaRISE experiments.
www.the-cryosphere.net/10/43/2016/The Cryosphere, 10, 43-63, 2016 The present paper revisits the future surface-climate experiments on the Greenland ice sheet proposed by the multimodel intercomparison SeaRISE (Bindschadler et al., 2013).A series of sensitivity experiments has been performed, using the ice-sheet model IcIES, to attempt to understand sources of the spread in the SeaRISE multi-model intercomparison.
Five aspects, surface balance parameterization, basal sliding, margin migration, initialization, and bed topography, are chosen to replace the standard formulation of IcIES by those adopted in other models, and all the experiments are conducted from spin-up to the simulation of future evolution.
The results show that the difference in the initialization methods as well as in the surface mass balance methods are large potential sources for the spread in the SeaRISE experimental results.In addition, the treatment of ice-sheet margin migration in the simulations also has a non-negligible impact on the spread among the multi-model projections.Performance of an initialization technique with fixed ice-sheet topography through time while temperature is allowed to evolve according to the surface temperature history or to the present-day condition is indirectly evaluated and found to provide an acceptable initial condition, at least for short-term projections.
The SeaRISE project, in which several ice-sheet models of different complexity participated to perform similar experiments, showed the degree to which current ice-sheet models and modeling choices diverge.Furthermore, Nowicki et al. (2013) show detail and careful analysis of all the results both globally and regionally to present how and where the models are similar or dissimilar.However, the SeaRISE protocol is not strictly controlled and most experimental configurations are left as the choice of the participants.Therefore, it is difficult to separate the effects of different choices by comparing only the submitted results.The present paper demonstrates that various implementations adopted in individual models can affect the simulated responses and how much they may contribute to the diversity in SeaRISE results.The analysis in the present paper is quite limited in terms of spin-up, and we propose a benchmark experiment to address this.If all models are used to perform a highly controlled experiment, it is easier to analyze the uncertainty due to model spin-up within the variation of current ice-sheet model structures.

Appendix A: Demonstration of the "benchmark" experiment
For a demonstration of the suggested benchmark experiment, configuration E s is performed by IcIES, which is the same as E s and E s except for the future surface mass balance scenarios.Steady-state initialization under fixed present-day topography is performed, and the future surface mass balance is imposed using the SeaRISE data sets without any correction.
Although most of the models did not use it, SeaRISE provided a transient future scenario of the surface mass balance computed by a variation of PDD method.The parameters for the PDD are described at http://websrv.cs.umt.edu/isis/index.php/Future_Climate_Data, where the standard deviation of the short-term statistical air temperature fluctuations is set as 4.5 K, and the PDD factors are set as 3 and 8 mm ice equivalent per day per degree for snow and ice, respectively.Actually, one participant, ISSM, in SeaRISE has a similar configuration to the benchmark: the surface mass balance is imposed with the SeaRISE data sets without any correction, initialization is based on inversion which enables initialization with a topography close to that of the present-day, and a fixed calving front is enforced (may correspond to prohibition of both advance/retreat).There is no explicit information about inclusion of the sub-melt sliding processes.The simulated response of VAF for this experiment is 5.4 cm sea-level equivalent at 500 years from the present under C1 scenario, which is the minimum response among the SeaRISE participants.
Figure A1 shows the simulated time series of VAF under C1 scenario with different uniform basal-sliding parameters v1 to v4, as well as runs using the inverted nonuniform basal-sliding field (Fig. 5b) with the default enhancement factor (vm) and a different enhancement factor E = 1 (e1:vm).The losses in VAF by IcIES are −10.8,−12.0, and −13.0 cm sea-level equivalent at 500 years with basal-sliding configuration of v1, v2, and v4, respectively; thus only 2.2 cm spread is attributable to the different basalsliding coefficient.Further, using the non-uniform basalsliding coefficient field leads to smaller losses in VAF: −9.0 and −6.7 cm sea-level equivalent for the vm and e1:vm cases, respectively.The smallest responses in the present paper are obtained under the E s configuration, which is even smaller than configuration E cases and is only 1.1 cm sealevel equivalent more than the smallest result of SeaRISE participants (ISSM, upper end of the gray bar in Fig. A1).Although the difference is very small, it is still possible that all the model aspects tested in the present paper are not sufficient to explain the SeaRISE spreads under future-climate scenarios.There are others differences in the properties such as higher-order physics, the numerical grid system, the basalsliding parameterization, and the distribution of basal-sliding coefficient field.Nevertheless, "net" uncertainties that stem from all the model properties except for those provided by external models (such as the surface mass balance) are expected be evaluated using this type of benchmark experiment.

-
F = fix to free topography + transient temperature, -F s = fix to free topography + steady-state temperature, -D = fix to observed topography + transient temperature, -D s = fix to observed topography + steady-state temperature.
Bindschadler et al. (2013)  presented their results in terms of the simulated time series of volume above flotation (VAF) under future-climate-warming scenarios, C1, C2, and C3, relative to that under the constant-climate scenario C0. Figure 1 shows the results of the present paper following the SeaRISE analysis under future-climate scenarios C1 with a standard basal-sliding coefficient v1. Figure 1 also shows the ranges of the results of the eight SeaRISE participants at 100, 200, and 500 years from the present, given in Table 3 of Bindschadler et al. (2013).The result of configuration O, which is a simulation corresponding to the original IcIES submission, is close to the largest response among the SeaRISE participants.

Figure 1 .
Figure 1.Simulated changes in VAF (volume above flotation, see the main text) obtained by future-climate runs under C1 (A1B climate forcing), with "standard" sliding coefficient (v1), in terms of the difference relative to the result of corresponding constantclimate experiments (C0).Each line is a different experimental configuration of O (IcIES SeaRISE compatible), A, B, B , D , and E .The vertical gray bars indicate the range of results summarized in the SeaRISE (Bindschadler et al., 2013, Table 3) at 100, 200, and 500 years.The circles in the gray bars indicate the mean values of all the SeaRISE participants.

Figure 2
Figure2shows simulated changes of VAF at 500 years obtained by a subset of the experiments in the present paper under the future-climate scenarios C1, C2, and C3 for the standard (v1), doubled (v2), and quadrupled basal-sliding coefficients (v4).The results of configuration O show volume losses of 34.1, 72.1, and 142.8 cm sea-level equivalent at 500 years under climate scenarios C1, C2, and C3, respectively.Standard basal-sliding cases v1 under all future-climate scenarios are within the range of original SeaRISE results.Simulated responses become larger with enhanced basal-sliding coefficient, and some cases are still within the original range of results, while some are above the range, for example, the simulated VAF response of C3:v4 is 17 cm more than the upper boundary of the original range.In the following sections, the effects of replacement of the five model aspects are described in turn.The fractional changes of the effects of this series of experiments are summarized in TableS1and Figs.S4-S6 in the Supplement.

Figure 2 .
Figure 2. Simulated changes in VAF at 500 years from the presentday obtained by future-climate runs in terms of the difference relative to the result of corresponding constant-climate experiments (C0).The top, middle, and lower panels are results of C1 (A1B climate forcing), C2 (1.5× A1B), and C3 (2× A1B), respectively.Each panel contains the results of experimental configuration of O (IcIES SeaRISE compatible), A, B, B , D, D , E, and E .The three bars from left to right in each configuration correspond to the results for v1 (using "standard" sliding coefficients), v2 (2×), and v4 (4×), respectively.The vertical gray bars at the right indicate the range of results summarized in SeaRISE (Bindschadler et al., 2013, Table 3) at 500 years.The circles in the gray bars indicate the mean values of all the SeaRISE participants.

Figure 3 .
Figure 3. Simulated present-day surface topography obtained by experiments with free spin-up initialization and sub-melt sliding: B (upper panels) and B (lower panels).Contour intervals are 200 and 1000 m for thin and thick lines, respectively.

FFigure 4 .
Figure 4.The same figures as Fig. 2 under experimental configuration of B, F, F s , D, D s , B , F , F s , D , D s .The left five experiments apply Tarasov and Peltier (2002) while the right five apply Huybrechts and de Wolde (1999) for the surface mass balance computation.

Figure 5 .
Figure 5.The result of the "inversion" procedure: (a) difference in the surface elevation (m) relative to the present-day observation and (b) "inverted" basal-sliding coefficient field in terms of fraction relative to the value of the case v2 in logarithmic scale.

Figure 6 .
Figure5shows the difference in the simulated surface topography relative to the present-day observation and inverted basal-sliding coefficient field.The inversion procedure is performed using the surface mass balance method of Huybrechts and de Wolde (1999), with the ice enhancement factor E = 1, and prohibition of advance in the ice-sheet margin.The last constraint is somewhat arbitrary but is kept for simplicity.The inverted coefficients are smaller than the case v2 value by some orders of magnitude over the interior region, while larger around the margin.Although not perfect, the overestimation in surface elevation near the margin and the underestimation in the interior part are significantly reduced (see Fig.S2d-f in the Supplement for uniform basal-sliding coefficient cases).As mentioned, since the inverted field is a function of other aspects such as surface mass balance, a different distribution should be computed for each configuration.Since the experiment reported in this section is intended to demonstrate non-uniform basal-sliding coefficient fields, the same field is used through all the experiment in this section.Among the series of experiments, E s and E s (sub-melt sliding included; fixed topography steady-state spin-up as the present-day observation; no advance in the ice margin) are performed using the inverted field, with default enhancement factors E = 3 and E = 1.Figure6shows the simulated changes in VAF under all of the combinations of climate scenarios and basal-sliding coefficients by the series of experiments E s and E s .The case E s :e1:vm is the simulation of most optimized present-day

Figure A1 .
Figure A1.Simulated changes in VAF obtained by future-climate C1 under experimental configuration of E s with uniform sliding coefficient cases v1, v2, and v4, with the inverted non-uniform sliding coefficient case vm, and that with ice enhancement factor E = 1 case e1:vm.
Four different future-climate experiments are presented inBindschadler et al. (2013): the surface-climate experiment, the basal-sliding experiment, the ice-shelf melting experiment, and a combination experiment.The present paper focuses on the surface-climate experiment, while the other three experiments are left for future studies.The surfaceclimate experiment leads to less abrupt changes after perturbation is applied than the other three, which is expected to emphasize the differences among various modeling approaches.In this future-climate experiment, changes in the climate conditions on the upper surface of the ice sheet are prescribed.Future scenarios of two fields, surface temperature and precipitation, are provided.The scenarios were calculated from the results of A1B scenario experiments by the mean of 18 climate models used in the Fourth Assessment Report, compiled by Bindschadler et al. ( www.the-cryosphere.net/10/43/2016/The Cryosphere, 10, 43-63, 2016 F. Saito et al.: SeaRISE revisited experiment design used for the submission is hereafter referred to as "IcIES" original configuration.

www.the-cryosphere.net/10/43/2016/ The Cryosphere, 10, 43-63, 2016 F. Saito et al.: SeaRISE revisited identical
to the present-day condition.Smoothing of the icesheet topography as used in some SeaRISE models is not applied for the present paper, in order to obtain the identical topography among runs with different model parameters.

Table 2 .
Simulated ice-sheet volume (×10 15 m 3 ), the percentage relative to the present-day observed volume 2.91 × 10 15 m 3 , and the root mean square of the difference in the thickness relative to the observation (m).Configurations correspond to the results for v1 (using "standard" sliding coefficients), v2 (2×), and v4 (4×) are shown.The volumes of other experiments such as D, E, etc. are identical to the observed value by definition.

Table 3 .
Simulated changes in VAF (cm) relative to corresponding constant future-climate experiments at 500 years from the present for the configurations B and F and their differences and the two configuration D and D and their differences.