Simulating the Greenland ice sheet under present-day and palaeo constraints including a new discharge parameterization

In this paper, we propose a new sub-grid scale parameterization for the ice discharge into the ocean through outlet glaciers and inspect the role of different observational and palaeo constraints for the choice of an optimal set of model parameters. This parameterization was introduced into the polythermal ice-sheet model SICOPOLIS, which is coupled to the regional climate model of intermediate complexity REMBO. Using the coupled model, we performed large ensemble simulations over the last two glacial cycles by varying two major parameters: a melt parameter in the surface melt scheme of REMBO and a discharge scaling parameter in our parameterization of ice discharge. Our empirical constraints are the present-day Greenland ice sheet surface elevation, the surface mass balance partition (ratio between total ice discharge and total precipitation) and the Eemian interglacial elevation drop relative to present day in the vicinity of the NEEM ice core. We show that the ice discharge parameterization enables us to simulate both the correct ice-sheet shape and mass balance partition at the same time without explicitly resolving the Greenland outlet glaciers. For model verification, we compare the simulated total and sectoral ice discharge with other estimates. For the model versions that are consistent with the range of observational and palaeo constraints, our simulated Greenland ice sheet contribution to Eemian sea-level rise relative to present-day amounts to 1.4 m on average (in the range of 0.6 and 2.5 m).


Introduction
Modelling the response of the Greenland ice sheet (GrIS) to anthropogenic warming has already been undertaken for more than 2 decades (Huybrechts et al., 1991;van de Wal and Oerlemans, 1997;Huybrechts and de Wolde, 1999;Greve, 2000) and attracted considerable attention in recent years (Vizcaíno et al., 2010;Goelzer et al., 2011;Graversen et al., 2011;Applegate et al., 2012;Lipscomb et al., 2013;Stone et al., 2013), including higher-order and full Stokes modelling approaches (Price et al., 2011;Seddik et al., 2012;Goelzer et al., 2013).The recent SeaRISE ice sheet modelling project (Nowicki et al., 2013) highlighted the importance of treatment of processes at the ice-ocean interface for the response of models of the Greenland ice sheet to future climate change.
Observational data indicate that during the past decade mass loss by the GrIS, both through surface melt and enhanced ice discharge, has contributed appreciably to global sea level rise (Shepherd et al., 2012).The latest projections suggest that the GrIS will contribute notably to sea level rise during the next century (Hanna et al., 2013).In the longerterm perspective, the GrIS can become even more important, because even if the global temperature is stabilized at the level of 2 • C above preindustrial, the GrIS would still continue to melt and, in the long-term perspective, can lose a significant fraction of its mass even for moderate warming (Ridley et al., 2005;Robinson et al., 2012).
Models of the GrIS contain a number of parameters that can be used for tuning the model using observational constraints.Present-day extent and surface elevation of the GrIS Published by Copernicus Publications on behalf of the European Geosciences Union.
are accurately known and it is natural to use them as such constraints (e.g.Stone et al., 2010;Greve et al., 2011;Lipscomb et al., 2013).At the same time, it is known that coarseresolution ice-sheet models have problems simulating the correct margins of the GrIS and they systematically overestimate its volume.One reason is that under present-day conditions, most ice discharge into the ocean occurs through relatively narrow outlet glaciers.As a result, although ice discharge into the ocean currently accounts for more than half of surface accumulation, the majority of Greenland the ice sheet margin is located several tens of kilometres away from the ocean.Since current ice-sheet models do not resolve outlet glaciers and their interaction with the ocean in the fjords, the modelled GrIS needs too much contact with the ocean to produce realistic discharge.This leads to systematic overestimation of the ice area and volume and makes observational "geometrical" constraints difficult to apply.
However, the observed shape of the GrIS is not the only characteristic that can serve as a constraint for the GrIS models.Recently, Robinson et al. (2011) introduced the mass balance partition (defined as the ratio between total ice discharge into the ocean and total precipitation over the GrIS, MBP) as another constraint on the GrIS.The MBP is an important characteristic for short-term as well as for longterm (future) behaviour of the GrIS, because it determines the GrIS mass balance sensitivity to climate change.In particular, the present-day MBP is related to the long-term stability properties of the GrIS (Robinson et al., 2012), i.e. for low MBP values (large surface melt) the modelled GrIS is more susceptible to warming than for high ones.For the longterm stability of the GrIS, the MBP is a more important characteristic than the present-day shape of the GrIS.However for short-term (centennial time scale) future global warming simulations, such an ice sheet would be an unfavourable initial condition.This is because during a considerable portion of time of such a future simulation, the modelled ice would melt in areas where in reality no ice exists (Goelzer et al., 2013).
Since standard coarse-resolution GrIS ice sheet models cannot simulate a realistic present-day surface orography of the GrIS and at the same time have the correct mass balance partition, we developed a novel approach, which allows us to circumvent this problem without resolving individual ice streams and outlet glaciers -and without an increase in computational cost.This approach is in the spirit of our previous modelling work (Robinson et al., 2010(Robinson et al., , 2011(Robinson et al., , 2012) ) and is based on a rather simple semi-empirical parameterization of ice discharge through the outlet glaciers.We propose the usage of this approach until a more complete representation of fast processes is available.Considering the above concerns, this approach is feasible for short-term as well as for longterm simulations.
In addition to the present-day constraints on the ice-sheet shape and the MBP, we use the Eemian as a palaeo constraint.Eemian conditions have already been recognized earlier as an important palaeo constraint for GrIS model parameters (Tarasov and Peltier, 2003) and have been applied more recently in several studies (Robinson et al., 2011;Born and Nisancioglu, 2012;Stone et al., 2013).While the Greenland Summit position might be located too far in the interior of the GrIS to serve as a strong palaeo constraint, the position of the new NEEM ice core appears more promising, because it is located rather near the ice margin, a location which is very sensitive to the climate changes during the Eemian (Born and Nisancioglu, 2012).In the present paper, we make use of the recently published estimate of the Eemian elevation drop at a position of about 200 km upstream of NEEM (NEEM community members, 2013), where the borehole ice sampled today was deposited during the Eemian.
The paper is summarized as follows.First, we give a short description of the ice-sheet model SICOPOLIS and the regional energy-moisture balance model REMBO (Sect.2).Our new discharge parameterization is comprehensively explained in Sect.3. Different metrics of the ice-sheet shape are introduced in Sect. 4 and the the model setup is given in Section 5. We discuss the behaviour of the shape metrics in Sect.6, while we present constraints on our model parameters in Section 7. Further simulations inspect Eemian climate and GrIS stability (Sect.8) and compare our findings with those of others (Sect.9).We close with a discussion and finally our conclusions.

Model description
For this study, we used the three-dimensional polythermal ice-sheet model SICOPOLIS (version 2.9) coupled to the regional energy-moisture balance model (REMBO).SICOPO-LIS treats the evolution of ice thickness, ice temperature and water content (Greve, 1997) based on the shallow ice approximation (Hutter, 1983).The dependence of the ice velocities on the ice temperature and water is introduced via the rate factor.SICOPOLIS enables a free and easy choice of several parameters including resolution.In our paper, Greenland is mapped onto a stereographic plane with 76 × 141 grid points (20 km grid spacing) using the topographic data set by Bamber et al. (2001).The vertical is resolved by 90 layers with decreasing layers thickness towards the bed of the ice sheet.A 10-layer thermal rock bed is coupled to the overlying ice sheet via heat fluxes.The geothermal heat flux is prescribed at the lower border of the thermal bedrock.We use the Weertman-type sliding law with powers as described in Robinson et al. (2011).The bedrock adjusts to the load caused by the ice sheet's weight using a local lithosphere relaxing asthenosphere model with a time delay of 3000 years.
The regional energy-moisture balance model REMBO is a climate model of intermediate complexity and it is described in detail by Robinson et al. (2010).REMBO uses diffusiontype equations for surface air temperature and atmospheric water content.For temperature, the well-known Budyko-Sellers energy balance approach is implemented.Planetary albedo is related to surface albedo via a linear parameterization based on empirical data.The lateral boundary conditions for temperature and relative humidity are taken from climatology (Uppala et al., 2005) for the 1958-2001 period, which in this paper is referred to as "present day".REMBO includes a 1-layer snowpack model with a simple parameterization of refreezing.Surface albedo depends on snow thickness and the melt rate.Surface melt is computed using a simple parameterization of van den Berg et al. (2008) and depends on both temperature and insolation.The formula for surface melt contains a free parameter c m (melt parameter), which is one of the major parameters determining the sensitivity of the ice sheet to climate change (Robinson et al., 2010).It reads which relates the surface melt m with the free melt parameter c m .The remaining variables ρ w , L m , τ s , α s , S, λ and T , are the density of water, the latent heat of ice melting, the total transmissivity, the surface albedo, the insolation at the top of the atmosphere, the long wave radiation coefficient and the surface temperature, respectively.It is important to note that REMBO resolves the seasonal cycle of temperature and insolation, which are both important for properly modelling surface melt (van de Berg et al., 2011).Please, refer to Robinson et al. (2010) for more details.The coupling between the models is bi-directional, i.e.SICOPOLIS provides the climate model with information about surface elevation and spatial extent of the ice sheet.In turn, REMBO provides SICOPOLIS with surface mass balance and mean annual surface temperature.

Ice discharge parameterization
Most ice discharge of the GrIS is brought into the ocean via fast-flowing narrow outlet glaciers (Rignot and Mouginot, 2012), most of which cannot be resolved by the coarseresolution GrIS models used for palaeoclimate simulations and long-term sea level change projections.This motivates us to include this portion of ice discharge in the mass balance of the ice sheet, which reads where h is the ice thickness, ∇ • q is the divergence of largescale lateral ice flow explicitly resolved by the model, d is the divergence of the ice flow associated with non-resolved ice streams/outlet glaciers which is parameterized (see below), b is the surface mass balance and b * is the basal mass balance, which is rather small.Note that, although b and d have the same units (m s −1 ) and are treated as surface fluxes, they represent completely different physical processes.The shows the ice-covered cells, while the dark gray shaded area R indicates the region over the ice sheet where the discharge parameterization applies.The (half) circle illustrates how that band R with the width of about R is determined in our scheme.Namely, the centre of such a circle is applied to every land point (open circles over the brown area).The smallest distance to the ocean l ij is depicted here for one example ice grid point.It is determined for every grid point inside the band R. For the discharge parameterization, the ice thickness h ij and the smallest distance to the ocean l ij are evaluated at grid points (i, j ), see Eq. (3).
integral of d over the entire area of the GrIS represents the total solid ice discharge into the ocean through outlet glaciers.
Total solid ice discharge into the ocean is equal to the sum of parameterized ice discharge and the lateral ice discharge into the ocean explicitly computed by SICOPOLIS.Under present day climate conditions the latter is much smaller than the parameterized solid discharge.This is explained by the fact that simulated modern GrIS has rather limited direct contact with the ocean.Note that Eq. (2) per definition conserves mass.
The divergence of the explicitly non-resolved fast lateral ice flow d in each grid cell (i, j ) is parameterized in a simple heuristic statistical approach via the local thickness of ice h and the distance from the actual grid cell to the nearest ocean grid cell l as where c, p and q are constant model parameters.R is defined as the area which consists of grid points located within a distance not more than R from the nearest ice-free land surface point or the nearest ocean point in case if the ice sheet has direct contact with the ocean (see Fig. 1).The rationale for this parameterization is the following.The ice thickness near the ice margins represents the amount of ice available to the outlet glaciers, while the distance to the coast can be regarded as a statistical measure of the outlet glacier density: if the ice margin is far away from the coast, it is very unlikely that any outlet glacier has contact with the ocean and there is only minor calving flux into the ocean, while one would expect a large calving flux for small distances to the coast.It is assumed that significant influence of outlet glaciers does not penetrate inside of the ice sheet (distance more than R).
Additionally, we assume that ice discharge into the ocean via fast moving ice streams only occurs from the area where ice surface is descending toward the coast.This is enforced by setting a maximum value of α 0 = 60 • to the angle between the gradients of surface elevation ∇ z s and distance field ∇ l (see Fig. 2 for the characteristics of the field l).Using the definition of the scalar product for the angle between above two gradients, this reads ∇z s • ∇l |∇z s ||∇l| ≥ cos (α 0 ) , evaluated at (i, j ). (4) This prevents our parameterization from simulating ice discharge into the ocean from an ice margin, which is oriented towards the interior of the Greenland island, which could happen when the GrIS is retreating under warm climates, and when only a few small ice sheets remain; see the small land bridge between the two ice caps in Fig. 3b for an example.The discharge parameterization is applied only to the icecovered grid cells that are located not more than R (here R = 120 km) from the ice margin (see Fig. 1).The resulting belt encloses the regions of ice with rather high velocities as found by satellite measurements (Joughin et al., 2010;Rignot and Mouginot, 2012).This reflects that our parameterization primarily accounts for the ice discharge via the near margin fast-flow and outlet glaciers.
The value of parameter c depends strongly on the choice of the powers p and q to maintain total ice discharge known form observational data.For convenience, we normalize c as i.e. for any fixed value of p and q we selected c 0 such that c d has a value of about one.In practice, after selecting p and q we chose c 0 so that the parameterization applied to observed Greenland elevation matched observed total discharge, for which we used 350 Gt yr −1 .The latter value is just about the average of the totals of ice discharge as found by Reeh (1994) and Rignot and Kanagaratnam (2006).Although the discharge for the modelled present-day GrIS will not be precisely 350 Gt yr −1 , such an approach guarantees that all valid values of the parameter c d keep the order of magnitude of about one for any power p and q.We thus have three free parameters c d , p and q in our ice discharge parameterization (Eqs.3 and 5).
As seen in Fig. 2, the minimal distance to the coast l can reach up to about 400 km in the centre of Greenland.In our parameterization, the inverse dependence of the ice discharge on the distance to the coast results in high ice discharge over regions very near to the coast (dark brown colour) and low discharge further inland (lighter brown colours).For a distance of around 400 km, the parameterized divergence of fast flow would nearly vanish.In our discharge parameterization, this is only relevant for an ice sheet under warmer climates, which can be much smaller than the present-day GrIS.In fact, because our parameterization only applies for R < 120 km around the ice margins, most parameterized divergence of fast flow occurs at the near-margin ice area for present-day Greenland (Fig. 3a).These are the regions where a high present-day ice discharge due to many outlet glaciers is observed (e.g.Moon et al., 2012), e.g. over the north-western region of the GrIS.For regions with fewer marine outlet glaciers, e.g. in the south-west, the observed ice margin resides rather in regions further from the coast.
In the model, the ice-marginal ring R and the distance field are computed every time step.The former is necessary, because the ice sheet changes its shape in time and because the coastline can potentially change its shape due to sea level change.

Measures of geometrical characteristics of an ice sheet
There are numerous possibilities to define a measure of the performance of a model based on the comparison of simulated geometrical characteristics of an ice sheet with observational data.The simplest is arguably to use the error in simulated total ice area and ice volume, which we define as where A obs , A mod , V obs and V mod are the observed ice area, the modelled ice area, the observed ice volume and the modelled ice volume, respectively.These errors can approach zero in principle, but this does not guarantee accurate simulation of the ice-sheet geometry, since regional errors can compensate each other.Therefore, we choose a stronger constraint based on the error in ice thickness expressed relatively to the total ice thickness.This reads where H obs ij and H mod ij are the observed and modelled ice thickness at the horizontal grid position (i, j ), respectively.The indices i, j run over the entire domain of the computational area and assume that the ice thickness is zero outside the ice-covered area.This error only approaches zero when the ice-sheet thickness is correctly simulated in each grid cell.Simulations assessed with the different error measures are presented in Sect.6.

Simulation setup
Following Robinson et al. (2011), we run the coupled REMBO-SICOPOLIS model through two glacial cycles starting at 250 kyr BP.These simulations serve a dual purpose: to perform a model spin-up necessary to simulate the present-day state of the GrIS and to apply palaeoclimate constraints (see below) to additionally reduce the range of model parameters.To drive the model through two glacial cycles, we apply variations in insolation due to changes in orbital parameters, equivalent CO 2 concentration and regional temperature anomaly obtained from the CLIMBER-2 model (Petoukhov et al., 2000;Calov et al., 2005;Ganopolski et al., 2010).We took these anomalies from the standard simulation as in Ganopolski and Calov (2011).The applied forcing is illustrated in Fig. 1 by Robinson et al. (2011).To generate an ensemble of model realisations, we vary two parameters: the discharge parameter c d (Sect.3) and the melt parameter c m (Sect.2).The discharge parameter c d is varied in steps of 0.2 and the melt parameter c m in steps of 5 W m −2 .The geothermal heat flux is set to 50 m W m −2 and the sliding coefficient to 15 m (yr Pa) −2 .All other parameters are the same as in Robinson et al. (2010).
We fix the powers p and q by minimizing the relative error in present-day ice thickness (see Sect. 4).Based on ensemble simulations over the parameter space (c d , c m ) for different powers p and q (not explicitly displayed here), we found that decreasing p and increasing q reduce the error in ice thickness err(H ).For simplicity, we chose the integers p = 1 and q = 3.To normalize the parameter c d (see Sect. 3 and Eq. ( 5)), we set the dependent parameter c 0 = 2.61 × 10 4 m 3 s −1 .All ensemble simulations presented in this paper are performed with these values for the powers and for the parameter c 0 in the ice discharge parameterization.

Simulations of the GrIS in the entire parameter space
Figure 4 compares the three error measures introduced in Sect. 4 for the ice-sheet shape in the phase space of the melt, c m , and discharge, c d , parameters.At first sight, the error fields in the three panels look similar.The smallest errors appear approximately along the descending diagonal, which is characterized by decreasing values of c m and increasing values of c d .One can also see that the parameter combinations with small errors are not limited to our sampled space: these regions expand in the direction of the descending diagonal.This underlines the need for more constraints (Sect.7). Figure 4 illustrates that our discharge parameterization allows us to reduce the errors in total area and volume practically to zero.However, we found no parameter combination for which the error in ice thickness was much lower then 20 %.Still, these are considerable improvements compared www.the-cryosphere.net/9/179/2015/The Cryosphere, 9, 179-196, 2015 to the standard version of the model without ice discharge parameterization constrained only by the mass balance partition (Robinson et al., 2011), which overestimates total ice volume and area by ca.20 % and has a relative thickness error of ca. 30 %.When the mass balance partition constraint is ignored, one can improve model performance in simulation of ice sheet shape by increasing surface melt.By choosing c m = −40 W m −2 , one can make all three errors comparable with those of the best model version with ice discharge parameterization.However, such a high melt factor practically eliminates ice discharge into the ocean and, as shown below, drastically affects the GrIS stability.In fact, it causes the GrIS to be unstable even under near-present-day climate conditions.

Constraining the model parameters
In addition to the relative error in present-day ice thickness, we use the following as further empirical constraints on the ensemble of the model realizations: the present-day surface mass balance partition and the Eemian drop in surface elevation relative to present day at the upstream position of the NEEM ice borehole.Figure 5a-c illustrates all constraints used in this paper.We accept a value of 20 % for the relative error in ice thickness.This choice is not totally arbitrary, because a closer inspection of the error field shows a minimum error in ice thickness of 18.2 %, i.e. there is indeed a plateau defined by ice thickness error values ≤ 20 %, as illustrated in Fig. 5a by the medium green shading.Within the parameter space, the error in ice thickness varies much more strongly for values higher than 20 %.This plateau structure of the error field motivated us to choose 20 % as the error limit.
As mentioned in the introduction, the mass balance partition is the amount of total ice discharge compared to total precipitation.In our work, we always refer to MBP as a characteristic of the ice sheet defined in its present-day state.Its practical definition is the total ice discharge divided by the total precipitation for the simulated present-day ice sheet.In Robinson et al. (2011), the MBP was diagnosed by REMBO from simulations with prescribed observed present-day icesheet topography.This was done because of systematic (and regionally significant) deviations of the simulated presentday GrIS from observational data.With the ice discharge parameterization, we can now safely operate with the simulated present day ice-sheet topography for determining MBP, because of the better match between simulated and observed topography.This means that our new MBP values simulated with ice discharge parameterization slightly differ from our former approach (c d = 0), and the valid MBP range (Fig. 5b) corresponds to somewhat lower values of c m compared to those by Robinson et al. (2011).However, the MBP has large inherent uncertainty, which we derived from results of regional climate models, following Robinson et al. (2011) and yielding a range of 45 to 65 %.
From measurement of air content in the NEEM borehole samples, the NEEM community members (2013) found that the surface elevation of the source area was 130 ± 300 m lower during the Eemian than at present day, which we exploit as our third constraint.Following these findings, we assume the maximal surface elevation drop at this location during the Eemian (130 to 115 kyr BP) did not exceed 430 m (compared to present day).Accounting for the trajectory tracing results by NEEM community members (2013), we defined the deposition position of Eemian ice in the NEEM ice core at the location about 200 km upstream from the NEEM drilling site at (45 • W, 76 • N) (see Fig. 7d), denoted NEEMup hereafter.We use the Eemian elevation drop at the NEEM source location as additional empirical constraint.
Figure 5 illustrates that none of the constraints is redundant, because the regions of valid simulations for all three constraints intersect each other and there is plenty of space without a common crossover for every constraint.discharge parameter, while the mass balance partition constrains the range of the melt parameter, the upper bound of which is then further constrained by the NEEMup elevation data.While the valid region of all three constraints cover an approximately equally large part of the parameter space (Fig. 5a-c), only a relatively small subset of model parameters (Fig. 5d) is consistent with all of these constraints simultaneously.
Figure 6 shows the range of ice margins which are consistent with the different constraints shown in Fig. 5.It can be seen that the MBP constraint alone (Fig. 6a) gives a rather broad band of valid ice margins, while the NEEMup constraint (Fig. 6b) alone results in an ice margin range, which is quite comparable with that of the err(H ) constraint (Fig. 6c).Finally, all three constraints together give a pronounced reduction of spread of ice margins (Fig. 6d) compared to the single-constraint cases.Simulated ice margins in the south, mid-west and north-east of Greenland compare well with observations.However, there are regions with rather strong mis-match in the south-west and in the north.Parts of this mismatch can be attributed to our model biases in precipitation.For example, REMBO simulates too much snow accumulation in the north-east and south-west of Greenland compared to the compilation by Bales et al. (2009).
Figure 5d depicts the simulated difference in GrIS volume between Eemian and present day expressed in units of global sea level equivalent.Compared to the other figure panels, we show here results from simulations with refined melt parameter spacing of 1 W m −2 .We enhanced the resolution of the melt parameter sampling, because the region of valid simulations appears rather elongated in the parameter space.The estimated Eemian sea-level contribution increases with increasing (less negative) c m .This is understandable, because surface melting increases with a higher (less negative) melt parameter.Nevertheless, there is also an increase of the GrIS contribution to the Eemian sea-level highstand for increasing discharge parameter values.Obviously, there is an interplay between ice discharge and surface melt, because the ice www.the-cryosphere.net/9/179/2015/The Cryosphere, 9, 179-196, 2015 discharge removes ice from the ice sheet and brings the ice surface into lower regions of the atmosphere, where stronger surface melt can occur.Averaged over the parameter space of valid simulations, we have a contribution of the GrIS to Eemian sea-level rise (above present-day value) of 1.4 m.The minimum contribution of the GrIS sea level rise among all valid simulations is 0.6 m and the corresponding maximum is 2.5 m.

Eemian versus present day and GrIS stability
Figure 7 shows the simulated present-day and Eemian ice sheet distributions from model versions with high, medium and low sea-level rise contributions of the Eemian compared to present day.While all fields look rather similar for the present day, there is a considerable difference between the corresponding Eemian fields.However, the present-day surface elevations for the different valid parameters sets still show slight differences.Naturally, these differences appear mainly near the ice margin, while the interior of the ice sheet remains almost unchanged for any valid parameter set.As is often the case in such optimization problems, there is a tradeoff concerning agreement with observations in certain regions (see Fig. 9a for the observed surface elevation).While the simulation with c d = 0.8, c m = −53 W m −2 (Fig. 7a) better resembles the ice-free south-western region, the northern region around Petermann Glacier matches the observation less well.This situation is opposite for the simulation with c d = 1.2, c m = −66 W m −2 (Fig. 7c).For all valid parameter sets, our simulated reduction in Eemian ice volume is accompanied by a strong retreat of ice in Greenland in particular in its northern part, see Fig. 7df, which spans the simulated lowest and highest Eemian to present-day GrIS contribution to sea level rise.For model versions with high sensitivity to climate forcing, the GrIS splits into two parts: a small ice cap in southern Greenland and a larger ice sheet in central Greenland (Fig. 7d).For the intermediate and low sensitivity model versions, the GrIS remains in one piece (Fig. 7e and f).In all valid model versions, there is a strong retreat of ice, mainly in western and northern Greenland.Interestingly, the NEEM location almost becomes ice free at 121 kyr BP in our most sensitive model version (see Fig. 7d).Nonetheless, an ice-free NEEM position during the Eemian would not contradict the existence of Eemian ice and most probably pre-Eemian ice in the NEEM ice core at present day, as reported by NEEM community members (2013), since the Eemian ice was accumulated farther upstream of NEEM.Similar argumentation would hold for Camp Century as well.
Figure 8 shows time series of ice volume and the NEEMup surface elevation for simulations over the last two glacial cycles from previous work with the same model (Robinson et al., 2011) and our present approach, which includes the sub-grid scale discharge parameterization.At all times, the valid model versions with the discharge parameterization simulate less ice volume than that without the discharge parameterization (Fig. 8a).This has two reasons: (i) previously, the model was not tuned for agreement with present-day surface elevation (ice volume).The present-day surface mass balance partition was (and is here) regarded as the more adequate characteristic to capture the sensitivity of the GrIS to long-term climate change.(ii) In our present approach, the inclusion of the discharge parameterization enables our rather coarse resolution model to mimic the calving of the small-scale outlet glaciers (i.e.removal in ice into the ocean by ice discharge) without an overestimation of contact regions of the ice sheet with the ocean, which leads to a smaller ice sheet.
Additionally, two extreme and unrealistic simulations, depicted by the red lines, were set up in order to demonstrate, what happens when a shape-only tuning applies in a coarseresolution model that disregards fast sub-grid processes of small outlet glaciers.Technically, we restrict the parameter space by setting c d = 0 (discharge parameterization off) and minimize the error measure err(H ) and, alternatively, the weaker error measure err(V ) to get the right present-day shape  2011) using observed present-day topography as well as outside the valid c m values in MBP space of this work (Fig. 5b) using simulated present-day topography to determine MBP.Because we consider the present-day ice-sheet shape as the only constraint (for demonstration), the model without the discharge parameterization (c d = 0) requires a rather high value of the melt parameter c m to minimized shape errors (Fig. 4).As one can see in Fig. 8a, around present-day the red line corresponding to minimal err(H ) is very close to the upper value of the range of the simulations with sub-grid discharge parameterization (blue shading), while the other red line (minimal err(V )) even merges with that valid range.Sim- ulated present-day elevations at NEEMup lie rather close to each other in different model versions.However, during the Eemian interglacial, the runs from the shape-only constraints show strong downward excursions for ice volume as well as for the NEEMup elevation (Fig. 8b).Whether such a small Eemian ice volume is still realistic might be disputable but in any case the simulated Eemian reduction in NEEMup elevation is by far larger than that estimated from the ice core (NEEM community members, 2013).The NEEMup position was even ice free during the Eemian in the simulation without the ice discharge parameterization and minimized err(V ), which certainly contradicts observational data.
Moreover, the strong drop in Eemian sea-level and NEEMup elevation hints at very different stability properties of the model version without ice discharge parameterization and shape-only tuning compared to all our valid model versions which contain the sub-grid scale discharge parameterization.Even more, the models with shape-only tuning are much less stable with respect to applied positive temperature anomalies than all the model versions that are constrained using the MBP and palaeo data, whether they include discharge parameterization or not.In other words, the models with shape-only tuning of the melt parameter are less stable than both the valid model versions of our former approach without the discharge parameterization (Robinson et al., 2011(Robinson et al., , 2012) ) and our present ones with the discharge parameterization.
To achieve more detailed information about the stability of the GrIS, we performed an analysis based on many steady state runs as in Calov and Ganopolski (2005), but in temperature space instead of insolation space.Namely, we performed a suite of steady state simulations each 300 kyr long imposing different spatially uniform and temporarily constant surface air temperature anomalies to the lateral boundary conditions of the REMBO model for each single simulation.We use the simulated present-day GrIS as an initial condition.We sample with a temperature increment of 0.25 • C, i.e. we add a T = 0.25 • C, 0.5 • C, . . ., 2.75 • C to the present-day temperature at the lateral boundaries of the REMBO model domain and run the model for each individual T into a steady state.If the GrIS decays for a certain T but covers most of the Greenland land surface (see Robinson et al., 2012) for T − 0.25 • C, we define T as the temperature threshold of GrIS decay.
We applied the procedure to three representative valid simulations with the discharge parameterization.From these simulations, we obtain thresholds of decay of the GrIS between 1.25 and 2.5 • C, depending on values of model parameters (see also Table 1).We also applied our stability analysis with the steady state method to a model version without the discharge parameterization using the shape-only constraint.The threshold estimated with this shape-only setting with err(V ) minimisation (c d = 0, c m = −40 W m −2 ) is much lower -only 0.25 • C. The higher values for GrIS decay between 1.25 and 2.5 • C from our new REMBO-SICOPOLIS version with discharge parameterization are in the range of those that we found previously without using the discharge parameterization (Robinson et al., 2012).This similarity clearly is an implication of the use of the MBP as one common constraint in both approaches.In this work, the shape and palaeo constraints are important as well, because they cover different regions in parameters space as discussed in Sect.7. It should be noted that a complete uncertainty analysis with our new model version is planned in the future, which will likely widen the range of our estimates.
In summary, if we optimize the melt parameter in the coarse resolution model without the sub-grid scale ice discharge parameterization for only err(H ) or err(V ), the re- sulting models all overestimate surface melt and violate the MBP criterion.This leads to a strong drop of the NEEM elevation far below the one reconstructed from palaeo data, and the Greenland ice sheet becomes too sensitive to temperature anomalies.This is why Robinson et al. (2011) and Robinson et al. (2012) used the MBP criterion together with a palaeo constraint for calibration of their coarse resolution model, ensuring the correct long-term stability properties reported in this work.In our improved model with sub-grid scale discharge parameterization, we found a stability behaviour similar to that found by Robinson et al. (2012) when using the MBP and NEEMup constraints, but now -still in the coarse resolution ice-sheet model -we can additionally fulfil a strong present-day shape constraint and achieve err(H ) < 20 %.We expect that all of our constraints would play a similar role in a model of the Greenland glacial system, which explicitly describes small-scale fast processes.Development of such a model is in our future plans.
One major advantage of our simple parameterization is that it applies easily for climates far away from present day -a fully explicit modelling of present-day outlet glaciers could fail for the Eemian, because many present-day outlet glaciers just vanish in the Eemian.Figure 3a and b compares simulated divergence of fast flow during present-day and the Eemian.
As stated in Sect.3, our discharge parameterization is not intended to resolve every small individual outlet glacier; it is rather designed to capture their draining effect on the spatial and temporal average in a sub-grid scale statistical approach.Consequently, the simulated ice discharge of a certain region is the integral of the divergence of fast flow (Eq. 3) over this region.In general, our simulated present-day ice discharge (Fig. 3a) is high over regions with many observed outlet glaciers and low over regions with fewer observed outlet glaciers, see Fig. 3 in Moon and Joughin (2008) for comparison.However, in the south-western part of Greenland, the model overestimates the (sectoral/regional) ice discharge, in part due to too high simulated accumulation over that region (see Sect. 9).
Furthermore, we demonstrate that the regions of fast flow can be reduced drastically for the Eemian time period compared to the present-day state.For the Eemian, there is practically no ice discharge from regions far away from the coast.In particular, the land bridge between the large ice sheet in the north and the smaller ice cap in the south of Greenland shows diminishing ice discharge.In general, our model results suggest that during the Eemian more ice calves into the ocean from the eastern coast of Greenland than from its western coast.In particular, the Kangerlussuaq Glacier region delivers ice into the ocean during the Eemian in all our valid model versions.

Comparison with present-day observations and findings by others
A direct comparison of our simulated Greenland surface elevation with the observed elevation by Bamber et al. (2001) and the former approach of Robinson et al. (2011) is shown in Fig. 9. Overall, we improved agreement with observations significantly.In particular, in the simulation with the discharge parameterization several regions are now ice-free, which look very similar to reality.The remaining deficiencies are partly due to the simple discharge parameterization and the limitations in the REMBO climatology, e.g.biases in representation of precipitation as discussed earlier.In this context, we would like to stress, that our model is a fully interactive one, where no observational data over Greenland are prescribed.This ensures that REMBO is applicable to climates far away from the present state, what is vital because the Eemian climate can deliver additional constraints for the model.
Figure 10 compares our simulated present-day sectoral and total ice discharge with findings by others.The sectors (see inset in the upper right of Fig. 10) correspond with those of Reeh (1994) and sub-divide the GrIS into northern (sector N), north-western (sector NW), north-eastern (sector NE), south-western (sector SW), and south-eastern (sector SE) parts.This subdivision is also adequate to the degree of complexity (Claussen et al., 2002) of our model in its current stage of development (a refinement of the sectors is planned for our later work).
Over the sectors N and SW, our simulated range of ice discharge compares well with previous estimates.While our simulated ice discharge range is somewhat low over sec-tor NW, it is certainly too high over sector NE.The latter can be explained by the overestimation of our simulated present-day accumulation over sector NE by some 10 Gt yr −1 compared to the compilation by Bales et al. (2009).In sector SE, our results are consistent with Reeh (1994) and Rignot et al. (2008) but are significantly lower than those by Enderlin et al. (2014).The reason for the mismatch of our simulated ice discharge over sector SE with that by Enderlin et al. (2014) could be because Enderlin's data, from the year 2000, already reflect the response of the outlet glaciers to warming in the fjords due to subtropic water transported by the ocean towards Greenland (Straneo et al., 2012).The SE sector may be particularly sensitive to this change, because in this region several outlet glaciers have contact with the ocean via fjords.
Our range of valid model versions is 326-479 Gt yr −1 simulated total ice discharge.The lower estimate matches that by Reeh (1994), while the upper end of this range is somewhat larger than the ice discharge by Enderlin et al. (2014).The relative small total ice discharge by Reeh (1994) cor- The Cryosphere, 9, 179-196, 2015 responds to the rather small accumulation estimate of the Ohmura and Reeh (1991) compilation, which Reeh used together with the assumption that the GrIS is in equilibrium to derive his discharge values.The data by Reeh (1994) can be regarded as roughly similar to pre-industrial, because it is based on accumulation, which contains several old data points, certainly before the 1990s.Meanwhile, the data by Enderlin et al. (2014) belongs to a more recent time, which corresponds to the upper value of the time interval of the reanalysis data used to prescribed the lateral boundary conditions of REMBO (1958REMBO ( -2001)).Because by the year 2000, Greenland already started to respond to global warming, it is reasonable that the upper end of our simulated ice discharge range approximately agrees with the finding by Enderlin et al. (2014).

Discussion
In spite of significant improvements in the simulated GrIS topography with our discharge parameterization, for all of our simulations it was impossible to yield an error in ice thickness smaller than about 18 %.These rather large errors partly underline the limits of our ice discharge parameterization and modelling approach in general.The errors can be reduced by incorporating additional parameters, in particular such parameters which affect the interior of the ice sheet, like the basal sliding parameter (e.g.Price et al., 2011;Larour et al., 2012a).In this paper, we restricted our analysis to the parameters for surface melt and ice discharge, which rather affect the marginal regions of the ice sheet.Nonetheless, we use the error metric of the entire ice sheet here in order to keep our approach simple, and because we intend to use the metric in forthcoming work, which will include a more comprehensive approach.
We designed this parameterization as a workaround until a more comprehensive whole-Greenland glacial system model becomes available.Of course, additional improvements are possible, like introducing physically based models for individual outlet glaciers and fjords.Nevertheless, note that the relative high error in ice thickness (up to 20 %) also results from the fact that this is a stronger measure of the error in ice-sheet shape than the error in total ice area or in ice volume.
Although our model enables us to reduce the cumulative error in ice thickness from about 30 to 18 %, there is still room for further improvement.For example, higher-order models can play an important role in capturing the GrIS dynamics (Larour et al., 2012b), as inclusion of membrane stress can reduce ice volume change in short-term climate projections by up to 20 % (Fürst et al., 2013).In particular, the role of stress transfer might become important in short-term climate projections (Nick et al., 2009).Inclusion of membrane stress (e.g.Bueler and Brown, 2009) in our modelling approach together with better resolved topogra-phy (Morlighem et al., 2014) is planned for the future.In addition, a comparison of our discharge parameterization with the high-resolution and high-order models will help to improve the parameterization, because it will still be useful for large-ensemble simulations.
The term d in the thickness evolution equation (Eq.2) represents the divergence of the ice flow associated with nonresolved ice streams/outlet glaciers, see also Appendix A. In this sense, d reduces the ice thickness over areas with fast flow.This is why we only consider total and sectoral ice discharge, which is equivalent to the spatial integral of d over areas near the ice margin where fast flow appears.Of course, our parameterization cannot explicitly simulate ice calving and we remove some mass over inland regions where in reality no calving appears.The distance parameter accounts for this, because it appears in the denominator of our parameterization and, therefore, our method removes the highest amount of mass near the ocean.One advantages of our method is that it enables us to distinguish between the mass that is removed by surface melt and by calving.
We restricted our simulations to the spatial horizontal resolution of 20 km and have not inspected a possible dependence of our ice discharge parameterization on resolution.We cannot rule out that a recalibration of the parameters will be necessary for a different resolution.For a better understanding, we plan to investigate the potential of resolution dependency in the future.At the same time, we regard such a parameterization as an important tool for use in a fixed resolution model.
While the model agrees reasonably well with observations overall, there are some significant biases in simulated ice discharge at the regional scale.For example, we have too much ice discharge in the north-eastern and too little in the northwestern sector.The disagreements can be partly attributed to regional biases of simulated precipitation by REMBO and to difficulties in interpretation of the data used for comparison.When designing our constraints, we took the reduction in Eemian surface elevation upstream of the NEEM ice core from the NEEM community members (2013).In their statistics, the NEEM community members gave a one σ error for this value.In principle, one could have included the more uncertain values too by using the two σ range.Nonetheless, all of our simulations with valid parameter sets show a strong retreat of the ice in northern Greenland during Eemian times.Such a retreat strongly influences the local climate and might lead to an additional Eemian temperature rise over that region, although it is unlikely to be as strong as that reported by the NEEM community members (2013).These and other uncertainties in Eemian temperature and precipitation will be examined in future work.
We introduced a new sub-grid scale ice discharge parameterization aimed at mimicking Greenland's fast outlet glaciers in a coarse resolution ice-sheet model.Our simulated ice discharge compares reasonably well with observations and other model estimates.The ice discharge parameterization enables us to simulate an ice sheet, whose shape is in good agreement with observations and whose partition between total ice discharge and total surface melt is in good agreement with stateof-the-art regional climate models.
We used various constraints to reduce the range of valid melt and discharge parameters of the REMBO-SICOPOLIS model: a shape constraint, a constraint on the mass balance partition (Robinson et al., 2011), and a palaeo constraint on Greenland's surface elevation drop (upstream of the NEEM borehole) during the Eemian interglacial compared to present.We favoured a measure of ice thickness error at each grid point instead of just considering total Greenland area or volume, since it is a stronger measure of the quality of simulated ice-sheet shape.
The NEEM constraint proved to be a complementary constraint to the other two present-day constraints.It was the strongest constraint in controlling the upper end of the range of valid melt parameter values and thereby the highest Greenland's contribution to Eemian sea-level rise.Taken individually, this constraint was also comparable to the shape constraint in determining the range of simulated present-day GrIS margins.This demonstrates the importance of palaeoclimate information for determining the range of model parameters applicable for future prediction of the contribution of the GrIS to sea level.
We can satisfy all constraints if our sub-grid scale ice discharge parameterization is included in a coarse resolution ice sheet model in order to mimic small-scale fast processes.When using a shape-only constraints in a coarse resolution model without the parameterization of fast processes, we obtained a very unstable ice sheet -i.e. a regional temperature rise of as low as 0.25 • C was sufficient to melt the GrIS almost completely on longer time scales.When applying the MBP constraint in a coarse resolution model without the sub-grid scale ice discharge parameterization, the model has about the same stability properties as with the discharge parameterization.
The inclusion of our ice discharge parameterization along with the above-described constraints leads to similar results concerning long-term stability as Robinson et al. (2012), with a decay threshold between 1.25 and 2.5 • C. Note that although this range is consistent with previous work (Robinson et al., 2012), it does not result from an exhaustive uncertainty analysis.An updated range comparable with Robinson et al. (2012) will be the provided in future work.Finally, complying with all three constraints leads to a GrIS contribution to sea level rise during the Eemian compared to present day in the range of 0.6-2.5 m, with an average of 1.4 m.Again, this range could widen if further uncertainties were included.

Appendix A: The ice discharge parameterization as divergence of fast flow
We stated that our ice discharge parameterization locally represents the divergence of fast flow.Here we cannot give a complete proof of that, however we can present a plausibility argument.
We start with rewriting Eq. ( 3) If we assume that the basal topography has a small spatial gradient, then the surface gradient equals approximately the gradient of ice thickness.Furthermore, let us rename l, the distance from a point in the ice sheet to the nearest ocean grid: l = x.Because the ice thickness at an ocean grid point is zero in our coarse resolution model, we can write h = h.With that we can very crudely define the surface gradient h/ x at the ice margin.We can then write Making x very small and using Eq.(A1) yields The fraction in the latter equation is a Weertman-type sliding law As sliding contributes by far the most to the velocity of ice streams, we can thus generalize to two dimensions where q fast = h u b .This we interpret as the divergence of the fast flow of outlet glaciers.Of course, this only corresponds to our rather simple approach.Nonetheless, we believe that the similarity of our parameterization to the divergence of a Weertman-type sliding velocity substantiates our approach.
The Supplement related to this article is available online at doi:10.5194/tc-9-179-2015-supplement.

Figure 1 .
Figure1.Principle sketch of the discharge parameterization over a part of the horizontal computational domain.The gray shading shows the ice-covered cells, while the dark gray shaded area R indicates the region over the ice sheet where the discharge parameterization applies.The (half) circle illustrates how that band R with the width of about R is determined in our scheme.Namely, the centre of such a circle is applied to every land point (open circles over the brown area).The smallest distance to the ocean l ij is depicted here for one example ice grid point.It is determined for every grid point inside the band R. For the discharge parameterization, the ice thickness h ij and the smallest distance to the ocean l ij are evaluated at grid points (i, j ), see Eq. (3).

Figure 2 .
Figure 2. Distance field over the entire Greenland area in km.It is determined by the minimal distance of every land grid point (icefree and ice-covered ones) to the coast (first ocean grid point, see Fig. 1).It defines the length l in the discharge parameterization (Eq.3).The yellow line indicates the ice margin of the present-day Greenland ice sheet.

Figure 3 .
Figure 3. Simulated parameterized divergence of fast ice flow in m yr −1 (Eq. 3) at (a) present day and, (b) at Eemian from the model version with the most reduced Eemian ice sheet.

Figure 4 .
Figure 4. Error measures for a modelled ice sheet in the {c d } × {c m } parameter space.Relative errors in (a) ice area, (b) ice volume and (c) ice thickness, all in percent.

Figure 5 .
Figure 5.Estimated constraints on the parameters c d (ice discharge) and c m (surface melt) illustrated with our ensemble simulations together with our estimate of the GrIS contribution to Eemian sea-level rise.(a) Relative error in ice thickness (%).(b) Mass balance partition (%).(c) Maximum elevation reduction during the Eemian compared to present day, 200 km upstream from NEEM (m).Here, regions where no Eemian ice is simulated at the upstream position of NEEM are displayed in white.(d) Simulated contribution of the Greenland ice sheet to sea-level rise between the Eemian and present day under our constraints.The black lines in (a)-(c) indicate our constraints: error values < 20 % for ice thickness, a 45 to 65 % range for mass balance partition and an Eemian to present-day surface elevation reduction of < = 430 m at the upstream position of NEEM.

Figure 6 .
Figure 6.Simulated geographical position of present-day ice margins for simulations with the discharge parameterization (gray areas) compared to observations (red line).The gray shaded areas cover the range of simulated ice margins determined by different constraints.(a)-(c) Single constraint applies: (a) mass balance partition, (b) elevation reduction during the Eemian referenced to present day at the upstream position of NEEM, and (c) error in ice thickness.(d) All three constraints apply.
Our estimates showing a strong retreat of the GrIS during the Eemian rather correspond to the simulations by Otto-Bliesner et al. (2006), while the medium to modest retreat of the Eemian GrIS was found in simulations by Helsen et al. (2013) and Quiquet et al. (2013).

Figure 7 .
Figure 7. Present-day (a-c) and Eemian (d-f) surface topography for varying melt and discharge parameters.Simulations correspond to those giving high medium and low contributions to Eemian sea-level rise (2.5, 1.5 and 0.6 m), respectively.The Eemian snapshots correspond to times with the simulated minimum ice volume during the Eemian for the respective simulation.NEEM locations are marked in magenta (square for borehole and circle for upstream).
. The former belongs to the parameter setting c d = 0, c m = −42 W m −2 and the latter to c d = 0, c m = −40 W m −2 .Please, note that these melt parameter values are outside the valid range of MBP as determined by Robinson et al. (

Figure 8 .
Figure 8.Time series of the simulated Greenland ice sheet evolution during the last two glacial cycles.Blue shading represents the range of valid model versions including our discharge parameterization.Black and red lines show simulations without the discharge parameterization (c d = 0).Solid black lines indicate the central run of a set of optimized simulations by Robinson et al. (2011).The red lines are from simulations found via shape-only tuning of the melt parameter (see main text for explanation).In particular, a simulation with c m = −42 W m −2 , found by minimizing err(H ) (solid red line) and, alternatively, with c m = −40 W m −2 , determined by minimizing err(V ) (dashed red line).(a) Ice volume of the Greenland ice sheet.(b) Surface elevation at the NEEMup location.

Figure 10 .
Figure 10.Simulated ice discharge (open bars) versus observations and findings by others (horizontal lines) at present day (i.e.pre-industrial conditions).The heights of the open rectangles indicate the range of our simulated discharge.Acronyms are as follows: Re94: Reeh (1994), Ri08: Rignot et al. (2008) and En14: Enderlin et al. (2014).Further on, SIM indicates our simulations.(a)-(e) Sectoral Greenland ice discharge in Gt yr −1 for the northern (N), north-western (NW), north-eastern (NE), south-western (SW), and south-eastern (SE) parts (colours of the sectors are indicated in the inset).(f) Total ice discharge in Gt yr −1 .Note that the y axes have different scales.

Table 1 .
Temperature threshold of the stability of the Greenland ice sheet for a number of valid model parameters.