A rapidly converging spin-up method for the present-day Greenland ice sheet using the GRISLI icesheet model

Providing reliable projections of the ice-sheet contribution to future sea-level rise has become one of the main challenges of the ice-sheet modelling community. To increase confidence in future projections, a good knowledge of the present-day state of the ice flow dynamics, which is critically dependent on basal conditions, is strongly needed. The main difficulty is tied to the scarcity of observations at the ice-bed interface at the scale of the whole ice sheet, resulting in poorly constrained parameterisations in ice-sheet models. To circumvent this drawback, inverse modelling approaches can be developed and validated 5 against available data to infer reliable initial conditions of the ice sheet. Here, we present a spin-up method for the Greenland ice sheet using the thermo-mechanical hybrid GRISLI ice-sheet model. Our approach is based on the adjustment of the basal drag coefficient that relates the sliding velocities at the ice-bed interface to basal shear stress in unfrozen bed areas. This method relies on an iterative process in which the basal drag is periodically adjusted in such as way that the simulated ice thickness matches the observed one. The process depends on three parameters controlling the duration and the number of iterations. The 10 best spin-up parameters are chosen according to two criteria to minimize errors in sea-level projections: the final difference between the simulated and the observed Greenland ice volume as well as the final ice volume trend which must both be as low as possible. To increase confidence in the inferred parameters, we also make sure that the final ice thickness root mean square error from the observations is not greater than a few tens of meters. Our best results are obtained after only 420 years of simulation, highlighting a rapid convergence and demonstrating that our method can be used for computationally expensive 15 ice-sheet models.


Introduction
Recent observations provide evidence that the rate of mass loss of the Greenland ice sheet (GrIS) is continuously increasing (Rignot et al., 2015;Mouginot et al., 2015).Simulating the GrIS response under future warm periods is therefore crucial to establish reliable projections of future sea-level rise at decade to century time scales (Bindschadler et al., 2013;Edwards et al., 2014), but also to investigate the effects of ice-sheet changes on the climate system Hansen et al. (2016); Defrance et al. (2017); Swingedouw et al. (2013); Böning et al. (2016).As a result, better constraining the GrIS evolution has become a key objective of the climate and ice-sheet modelling communities.Reliable simulations of the GrIS require a proper initialisation (i.e.spin-up) procedure to avoid an unrealistic evolution of the ice sheet caused by inconsistencies between the ice-sheet model initial conditions and the external forcing fields.For short-term projections (next decades to next centuries) starting from the present-day ice-sheet configuration, recent observations, such as surface and bedrock topographies (Bamber et al., 2013) and horizontal surface velocity (Joughin et al., 2010) offer only a partial description of the GrIS characteristics and the major source of uncertainty lies in the poor knowledge of the basal properties, such as the water content in the sediment and basal sliding, and of the vertical temperature profile.Indeed, the basal conditions have a strong impact on the ice motion (Boulton and Hindmarsh, 1987;Weertman, 1957;Kulessa et al., 2017).Optimizing the initial conditions of ice-sheet models is therefore an active area of research and a multidisciplinary effort.The initMIP project (Goelzer et al., 2017) gives a recent example of this effort.Its goal is to compare different initialisation techniques and to assess their impact on the dynamic responses of the models.Three main classes of initialisation techniques have been developed: 1.The free spin-up method allows the ice-sheet topography to evolve freely over a long enough time.This approach has long been the most commonly used technique to initialise ice-sheet models (Charbit et al., 2007;Huybrechts and de Wolde, 1999;Huybrechts et al., 2002, and other refecerence in (Rogozhina et al., 2011)).It consists in simulating the ice sheets during one or more glacial-interglacial cycles to account for the long-term ice-sheet history and thereby to obtain internal consistency between the simulated ice sheet and the climate forcing derived from ice core records.Since the ice-sheet topography evolves freely during the entire spin-up experiment, this may lead to a significant mismatch between modelled and observed presentday ice-sheet topography.Such spin-up methods can only be used with low computational cost models, which are often unable to properly capture fast ice flow processes.
2. The fixed topography spin-up method is similar to the free spin-up method except that during all the simulation the ice-sheet topography is kept constant and equal to its present-day observed value, while vertical temperature is allowed to freely evolve (e.g.Seddik et al., 2012;Sato and Greve, 2012).The disadvantage of this method is that an artificial drift may arise when free evolving topography is restored due to inconsistencies between internal and surface ice sheet fields.

3.
The third kind of spin-up technique is based on an inverse method of the poorly known basal conditions in such a way that simulated surface velocities match the observed surface velocities (e.g.Gillet-Chaulet et al., 2012;Arthern and Gudmundsson, 2010;Morlighem et al., 2010;Gudmundsson and Raymond, 2008).However this approach may lead to internal inconsistencies between the simulated internal conditions (temperature and velocities) and the actual ones.The inconsistencies within the different observational datasets (surface and bedrock topography, velocities) can also have an impact on the results.An alternative approach, which avoids the previously mentioned shortcomings, consists in considering the observed ice-sheet geometry as the final target by finding appropriate basal conditions that minimise the differences between observed and simulated ice thickness (Pollard and DeConto, 2012).
Here, we present a spin-up approach that relies on the same basic principles as those developed by Pollard and DeConto (2012) (referred to as PDC12 in the following) for the Antarctic ice sheet.Similarly to PDC12, we compute the basal drag coefficient that minimises the error in the simulated ice thickness and relates basal stresses to basal velocities.However, while PDC12 requires long (multi-millennial) integrations for the method to converge, we suggest instead an iterative method of short (decadal to centennial) integrations starting from the observed ice thickness.Our method ensures a more rapid convergence and is thus suitable for computationally expensive models.The paper is organised as follows.In section 2 we present the main characteristics of the GRISLI ice-sheet model used in this study.Section 3 describes the spin-up method in detail.The main results are presented in section 4 and sensitivity experiments in section 5.These sections are followed by a discussion and the main conclusions of the present study.

The ice-sheet model GRISLI
The GRISLI ice-sheet model was first designed to describe the Antarctic ice sheet (Ritz et al., 2001) and further adapted to the northern hemisphere ice sheets (Peyaud et al., 2007).The version used in this study has been specifically developed for Greenland (Quiquet, 2012) with an horizontal resolution of 5 km x 5 km (301 x 561 grid points) and 21 vertical unevenly spaced levels, with the smallest grid spacing near the ice-bedrock interface to better resolve the basal motion.GRISLI is a hybrid model accounting for the coupled behaviour of temperature and velocity fields.It relies on basic principles of mass, heat and momentum conservation.The evolution of ice-sheet geometry is a function of surface mass balance, velocity fields and bedrock altitude.Since this study only deals with present-day steady-state simulations, the module describing the isostatic adjustment is not activated here.The evolution of the ice thickness is governed by the mass balance equation: where H is the ice thickness, U is the depth-averaged velocity, M is the surface mass balance and b melt is the basal melting.
The ice flow velocity is derived from a simplified formulation of the Stokes equations (i.e. the stress balance) using the shallow-ice (Hutter, 1983) and shallow-shelf (MacAyeal, 1989) approximations.The shallow-ice approximation (SIA) assumes that, owing to the small ratio of vertical to horizontal dimensions of the ice sheet, longitudinal stresses can be neglected with respect to vertical shearing along the steepest slope.Conversely, in the shallow-shelf approximation (SSA), the horizontal strain rates become dominant and the horizontal velocities do not vary with depth.In the model, the velocities are computed as the sum of the SSA and the SIA components, as in Bueler and Brown (2009), with the SSA velocity used as the sliding velocity.
We assume no-slip conditions for a frozen bed (i.e.basal temperature below the melting point), and in these conditions, the SSA velocity is set to 0. In the model version used in this study, we assume a linear till, in which the basal shear stress (τ b ) and basal velocity (u b ) are related via the following expression: where β is the basal drag coefficient and varies with space.
To describe the effect of ice rheology, the deformation rate and stresses are related via the Glen's flow law (Glen et al., 1957).As in other large scale ice-sheet models, GRISLI uses a flow enhancement factor in the Glen's flow law to artificially account for the impact of ice anisotropy on the deformation rate.This enhancement factor (Ef) typically ranges from 1 to 5. The grounding line position is defined according to a flotation criterion and floating points are treated following the SSA assumption only.Calving is not explicitly computed, but the ice-shelf front position is determined for a cut-off criterion of 250 m (Peyaud et al., 2007).The amount of ice obeying this criterion (ice thickness > 250 m) is computed as the calving flux.Numerous studies are based on fitting the modelled ice velocities (e.g.Gillet-Chaulet et al., 2012;Arthern and Gudmundsson, 2010;Morlighem et al., 2010;Gudmundsson and Raymond, 2008), while Pollard and DeConto (2012) opted for fitting ice surface elevation.Here, we decided to adjust the basal sliding velocities via the adjustment of the β coefficient to fit the GrIS The GRISLI climate forcing is provided by the surface mass balance and the surface air temperature simulated by the state of the art regional atmospheric model MAR (Fettweis et al., 2013) forced at its boundary by the ERA-Interim reanalyses (Berrisford et al., 2011).Both fields are averaged over the 1979-2014 period (Fig. 1).They are interpolated on the GRISLI grid (5 km x 5 km) and corrected for surface elevation differences between MAR and GRISLI by applying the method developed by Franco et al. (2012).We use the reconstruction from Maule et al. (2005) for the geothermal heat flux.Using these boundary conditions, GRISLI is run forward starting from the present-day observed ice thickness (Fig. 2a), from which the ice volume is inferred Bamber et al. (2013), and from the bedrock elevation.Initial vertical temperature and velocity profiles as well as the initial map of the basal sliding coefficient (Fig. 2a) are derived from previous GRISLI simulations carried out with  the different datasets used as boundary and initial conditions, GRISLI is first run for 5 years.After this relaxation period, we start the spin-up procedure.This procedure is based on an iterative process set up to adjust the basal drag coefficient in such a way that the mismatch between observed and simulated ice thickness is reduced.At the end of the iterative process, we allow GRISLI to evolve freely for 200 years in order to assess the model performance in terms of trend and error in simulated ice volume compared to observations.The iterative process itself is divided in two main steps (Fig. 4).
U corr can be seen as a the vertically-averaged velocity field corrected by a factor representing the difference between the observed and the simulated ice thicknesses.As seen before (section 2), the mean velocity field U G in GRISLI is the sum of two velocity components: the sliding velocity U sli and the velocity U def due to vertical ice deformation: Considering that the differences of velocity between U G and U corr are only due changes of the sliding velocity U sli , we can also write: Following Eqs. ( 4) and ( 5), we can deduce the corrected sliding velocity (U corr sli ) needed to reduce the difference between H G and H obs : The new value of the basal drag coefficient allowing to reduce the gap between H G and H obs is deduced from the β old value, inferred from the previous iteration and from the ratio between uncorrected and corrected sliding velocities: with β new calculated at each GRISLI grid point.H G , U G , U sli corr andβ new are updated during Nb iter time steps.2 nd step: With this new basal drag coefficient we let the model to freely evolve.After Nb year of the free-evolving simulation, we obtain a new GrIS topography and new corrected velocity fields computed from the mismatch between the simulated ice thickness after Nb year and the observations.With this, we can start a new cycle in which the 1 st and 2 nd steps are repeated.
This new cycle uses the same set of spin-up parameters (Nb iter and Nb year ) and an initial guess of β coming from the previous iteration.All the iterations use the same initial conditions presented previously.The number of cycles carried out in this way is noted Nb cycle .For all the experiments presented in the following, we performed a maximum of nine cycles.
To assess the spin-up performance and the quality of the inferred β coefficient, we perform, at the end of each cycle, a free-evolving simulation of 200 years with a β coefficient fixed to the values computed during the last cycle.The best Nb cycle for a given set of Nb iter and Nb year will be the one that provides a final volume as close as possible to the observation and a minimal trend over the last ten years of this free-evolving simulation.In addition to these two criteria, we also take into account the ice thickness root mean square error from the observations.Once the best fit is obtained, steady-state or transient GrIS simulations can be performed with reliable initial conditions, as done in Le clec 'h et al. (2017) and Goelzer et al. (2017).The annual mean climatological SMB for the 1979-2014 period integrated over the whole GrIS is 381 Gt yr −1 (Fig. 1a) with strongly positive values in southeastern Greenland (up to 0.04 Gt yr −1 ), and largely negative ones over the ablation zone at the edges of the ice sheet, with values reaching -0.10 Gt yr −1 in the western area (Fig. 1a).The annual mean surface temperature is negative over the whole Greenland ice sheet, ranging from -29 • C in the highest altitude regions to -0.5 • C near the coast (Fig. 1b).To illustrate the need for a spin-up procedure, we performed a 200-year-long free-evolving simulation with this mean climatic forcing and with the initial conditions obtained after the 5-year-long relaxation (see Sect. 3), without any spin-up procedure.The simulated GrIS ice volume obtained in this experiment is smaller, by 20000 Gt, than the one estimated by Bamber et al. (2013) from observations (i.e.2.71 10 6 Gt) with negative ice thickness anomalies in the interior of the ice sheet, which are even stronger in the northwestern, northeastern and central eastern parts (Fig. 2b).Moreover, the ice volume decrease contributes by 0.6 mm yr −1 to the global sea level rise.Thus, despite an overall positive SMB, the model drift and the lack of spin-up procedure result in a decrease of the GrIS volume as large as the present melting due to the global warming Church et al. (2013).Therefore the use of a spin-up method to minimise the model drift is not avoidable if the goal is to produce reliable sea-level projections.

Spin-up performance
To assess the spin-up performance, we first examined the ice thickness root mean square error (RMSE), the ice volume anomaly

Root mean square error
The RMSE behavior is approximately the same for all the Nb year values (Fig. 5).A strong decrease between the first two cycles is obtained meaning that the departure between simulated and observed ice thickness is rapidly reduced.This decrease is then followed by a stabilisation occurring for Nb cycle between 4 and 6 depending on the Nb year value.Increasing Nb year results in lower RMSE values.For example, for Nb cycle = 6, the RMSE decreases from 84.8 m (Nb year = 50) to 57.4 m (Nb year = 400).
This can be explained by the fact that for longer free-evolving simulations, the basal velocity (computed through the previously determined β coefficient) exerts a longer influence on the vertically averaged velocity, which in turn impacts the simulated ice thickness.This results in larger differences between simulated and observed ice thickness.This implies that the corrections of the β coefficient are more significant for the following cycle, and finally, that the method is more efficient to correct for the differences of the ice thickness with respect to the observed one.

Volume
The ice volume anomalies (∆Vol) obtained for the same set of spin-up parameters are displayed in Figure 6a.A strong ∆Vol decrease is obtained, starting from a highly positive value (∆V ∼28 000 Gt for Nb cycle = 1) and reaching negative values when Nb cycle increases.This illustrates that our method tends to underestimate the ice volume with respect to observations Bamber et al. (2013).This underestimation is also more pronounced for higher Nb year values.As a consequence, the combination of spin-up parameters providing the lowest RMSE values are also those for which the ice volume anomalies are the largest ones.

10
For example, for Nb cycle = 6, ∆Vol equals ∼10 000 Gt for Nb year = 400, while it is two order of magnitude below (∆Vol = 107 Gt) with Nb year = 50.
This behaviour can be explained when examining the ice thickness anomaly (∆H, Fig. 7a).This anomaly tends to be negative

Ice volume trend
This analysis demonstrates that satisfying the ice volume anomaly criterion is in direct competition with the RMSE minimi-  Even if our approach is based on minimising the volume trend and fitting ice volume and ice thickness to the observed ones, the reliability of the method also depends of its capability to simulate ice velocities in good agreement with observations.
We therefore compare our results to the surface ice velocity dataset provided by Joughin et al. (2010) (Fig. 8).The simulated results are slightly different from the observations especially in the central plateau where the region of low ice velocities is less extended in the simulations than in the observations.However, the overall patterns are in a good agreement, especially in regions of fast ice flow, providing confidence in our method.

Sensitivity to Nb iter values
Results obtained for other Nb iter values (40, 80 and 160) are reported in Table 1.None of these experiments fulfill both the ∆Vol and IVT criteria.The Nb 40 iter -Nb 100 year -Nb 3 cycle provides the lowest ice volume trend (5.4 10-3 mm yr −1 ) but the simulated ice volume anomaly is more than 20 times larger than for Nb 20 iter -Nb 50 year -Nb 6 cycle .Conversely, Nb 160 iter -Nb 200 year -Nb 2 cycle simulates a GrIS ice volume anomaly only 2.5 as large as for Nb 20 iter -Nb 50 year -Nb 6 cycle but the ice volume trend (0.041 mm yr −1 ) is four times larger.Moreover, the duration of the spin-up procedure for Nb 160 iter -Nb 200 year -Nb 2 cycle is 1080 model years while it is only 420 years for Nb 20 iter -Nb 50 year -Nb 6 cycle .This confirms that Nb 20 iter -Nb 50 year -Nb 6 cycle appears as the best protocol in terms ∆Vol and IVT criteria and that it is also designed to ensure a more rapid convergence.
5 Sensitivity to initial conditions and model parameters

Temperature equilibrium
In the work presented in section 4, the vertical temperature and ice velocities profiles taken as initial conditions came from previous experiments carried out with GRISLI in the framework of the Ice2Sea project (Edwards et al., 2014).These profiles have been chosen because they were assumed to reflect the present-day conditions.However, they are not necessarily in equilibrium with the climatic forcing taken from the MAR simulations (Fig. 1).Indeed, at the ice-sheet surface, the temperature obtained   from MAR is about 5 • C warmer than the Ice2Sea one (Fig. 9a).Therefore, we performed a 30 000-year-long simulation to make the vertical temperature and velocity profiles consistent with the surface climate forcing.This new experiment has been carried out with the ice topography fixed to the observed one (Bamber et al., 2013)  These results illustrate the limitations of our spin-up method, as also shown in Figure 9b which displays the ice thickness anomaly for Nb 20 iter -Nb 50 year -Nb 6 cycle (∆Vol = -56102 Gt, IVT = -1.19mm yr −1 , RMSE = 118.3m) after a 200-year freeevolving simulation.Actually, the warmer temperatures obtained after the 30 000 year-long simulation induce higher ice 5 velocities due to the thermo-mechanical coupling.In the ice-sheet interior (SIA areas), these velocities are not corrected by our spin-up approach as shown by the β coefficient which reaches maximum values in these regions (Fig. 3b).Thus, an increased ice flux takes place from the central part to the peripheral regions leading to amplified negative ice thickness anomalies (Fig. 9b).

Sensitivity to the enhancement factor
In the ice-sheet interior, the ice flow is mainly due to internal ice deformation which is controlled by the temperature and thus by the viscosity.A possibility to reduce errors in ice surface elevation in these locations is to adjust the enhancement factor of the Glen's flow law, which relates viscosity to deformation rates.Lowering the Ef value allows to decrease the deformation and thus to slow down the ice flow velocities.Therefore we performed new sensitivity tests with Ef = 1 (instead of Ef = 3, as in the experiments presented in section 4) with the same spin-up parameters used in the previous section.
As expected, for given Nb iter and Nb year values, the ice thickness RMSE is improved when the number of iterations outlet glaciers and for their interactions with floating ice that strongly influence the ice-sheet mass balance (e.g.Aschwanden et al., 2016).While this effect is less crucial for Greenland than for Antarctica, recent observations have highlighted increasing thinning rates in most coastal regions (Thomas et al., 2009;Rignot et al., 2015) causing grounding line retreat and significant destabilization of grounded glaciers.
The reliability of the method also depends on the quality of observations data.Errors in observed surface or bedrock topography would give rise to errors in the present-day estimated ice thickness and thus to an erroneous choice of the best spin-up parameters.In the same way, large uncertainties remain in the reconstruction of the geothermal heat flux that strongly impacts the basal temperature.Finally, we would like to stress that in our simulations, the spatial distribution of the basal drag coefficient does not change through time.However, changes in basal hydrologic conditions along with changes in ice surface elevation and ice extent are likely to occur in a changing climate.While a constant spatial distribution of the β coefficient may seem reasonable for short-term projections, it is more questionable at the century time scale, and future modelling efforts should therefore be undertaken to compute interactively the basal drag coefficient as a function of changes in basal conditions.

Code and data availability
The developments of the GRISLI source code are hosted at https://forge.ipsl.jussieu.fr/mailman/listinfo.cgi/grisli.Access to those who conduct research in collaboration with the GRISLI users group can be granted upon request to C. Ritz (catherine.ritz@univgrenoble-alpes.fr).The model outputs from the simulations described in this paper are freely available from the authors upon request.

Figure 1 .
Figure 1.Climate forcing averaged over the 1979-2014 period simulated by the atmospheric regional model MAR (Fettweis et al., 2013) and interpolated on the GRISLI ice-sheet model grid (5 km x 5 km).a/ Mean surface mass balance (in Gt yr −1 ).The black line represents the equilibrium line indicating the frontier between accumulation and ablation areas.b/ Mean annual surface temperature (in • C).The white dashed lines represent the 5 • C isocontours.
Dev. Discuss., https://doi.org/10.5194/gmd-2017-322Manuscript under review for journal Geosci.Model Dev. Discussion started: 8 February 2018 c Author(s) 2018.CC BY 4.0 License.ice thickness to the observed one.Our choice is motivated by the need to refine the estimates of GrIS contribution to future sea-level rise.
10boundary conditions close to those of the present study, and performed within the Ice2Sea project, which aimed at reducing the uncertainties on future sea-level rise projectionsEdwards et al. (2014).In order to avoid large inconsistencies between

Figure 2 .
Figure 2. a/ Observed Greenland ice thickness (in m) from Bamber et al. (2013) interpolated on the GRISLI grid.Grey areas represent non ice-covered areas.b/ Difference between the simulated and the observed ice thickness (in m) obtained at the end of a 200-year-long simulation without spin-up procedure.The simulation has been carried out using the Ice2Sea initial conditions (see main text and Edwards et al. (2014)) and the climate forcing simulated by MAR.

Figure 3 .
Figure 3. Spatial distribution of the basal drag coefficient (in log10 Pa m −1 an) a/ for the initial condition as used in the GRISLI ice2Sea simulations, b/ obtained for the best fit Nb 20iter -Nb 50 year -Nb 6 cycle and c/ obtained at the end of a spin-up procedure using the same spin-up parameters as those inferred from the best fit but starting from a uniform spatial distribution of the basal drag coefficient (β = 1).
Figure 4. Schematic representation of the spin-up method.
Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-322Manuscript under review for journal Geosci.Model Dev. Discussion started: 8 February 2018 c Author(s) 2018.CC BY 4.0 License.Four values of Nb iter (20, 40, 80, 160 years) and of Nb year (50, 100, 200 and 400 model years) have been tested with Nb cycle ranging from 1 to 9, giving a total of 144 combinations of the spin-up parameters.The corresponding simulations are referred to as Nb X iter -Nb Y year -Nb Z cycle where X, Y and Z stand respectively for the Nb iter , Nb year and Nb cycle values.the spin-up needed?

(
computed -observed) and ice volume trend for each combination of the spin-up parameters.An illustrative example is given here for Nb iter = 20 and Nb year ranging from 50 to 400 years, with Nb cycle varying from 1 to 9 (Figs5 and 6).For Nb cycle = 1, all the tests corresponding to different Nb year values start from the same initial conditions and the basal sliding velocity has not yet been updated with the new β coefficient.

Figure 5 .
Figure 5. Ice thickness root mean square error w.r.t.observations from Bamber et al. (2013), in meters for Nbiter = 20 and the four Nbyear values (50, 100, 200, 400) as a function of the number of iterations (Nb cycle ).

Figure 6 .Figure 7 .
Figure 6.Same as figure 5 for the GrIS ice volume anomaly (a) and the ice volume trend (b) represented in Gt and mm yr −1 respectively.
10 sation, and therefore that a compromise needs to be found.Since our main objective is to obtain a simulated ice volume as 11 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-322Manuscript under review for journal Geosci.Model Dev. Discussion started: 8 February 2018 c Author(s) 2018.CC BY 4.0 License.close as possible to the observed one (Bamber et al., 2013), the criterion related to the ice volume trend (IVT) must be also examined.Figure 6b shows that this trend follows a behavior similar to the ice volume anomaly with increasing values of Nb year and Nb cycle .The key feature appearing in this figure is the strong decrease towards negative values (down to -0.5 mm yr −1 ) for most of the spin-up parameter combinations.However small IVT values are obtained for three set of parameters: Nb year = 100 and Nb cycle = 3 (IVT = -0.04mm yr −1 ), Nb year = 200 and Nb cycle = 2 (IVT = +0.03mm yr −1 ) and Nb year = 50 and Nb cycle = 6 (IVT = 8.5 10 −3 mm yr −1 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-322Manuscript under review for journal Geosci.Model Dev. Discussion started: 8 February 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 8
Figure 8. a/ Observed surface ice velocities coming from a compilation of interferometric synthetic aperture radar measurements obtained from RADARSAT data at different periods of the 2000s.b/ Simulated surface ice velocities obtained for our best fit (Nb 20 iter -Nb 50 year -Nb 6 cycle ).Values are given in log10 m s −1 .

Figure 9
Figure 9. a/ Vertical profile of temperature (in • C) in a central region of Greenland (73-74.5 • N, 40-43 • W) taken as initial condition from the Ice2Sea project (Edwards et al., 2014, see black dashed line) at the end of the spin-up procedure (Nb 20 iter -Nb 50 year -Nb 6 cycle ) and at different periods of the temperature equilibrium experiment.b/ Difference between simulated and observed Ice thickness (in m) obtained after the temperature equilibrium and a new spin-up procedure performed with the spin-up parameters providing the best fit fit in terms of ice volume and ice volume trend.

(
Figure 10.Difference between simulated and observed Ice thickness (in m) obtained for the spin-up parameters (Nb 20 iter -Nb 400 year -Nb 8 cycle ) providing the best fit when Ef = 1.
).While these two first combinations provide reasonable RMSE values (see Table 1), the ice volume anomalies are respectively 1079 and 4343 Gt.This must be compared to ∆Vol = 107 Gt obtained with Nb year = 50 and Nb cycle = 6 which provide the smallest ∆Vol and IVT values.Despite the corresponding RMSE being the highest one among all the tests which have been performed, it is less than 20 m higher than the lowest RMSE value.Thus the experiment Nb 20 iter -Nb 50 year -Nb 6 cycle matches our two main criteria (i.e.minimum ∆Vol and IVT values) and the corresponding spin-up parameters appear as good candidates for the overall spin-up procedure.

Table 1 .
Ice volume trend (in mm yr −1 ) and ice volume anomalies (simulated -observed) obtained for the 16 combinations of the NBiter and NByear parameters.The values of NB cycle correspond to the number of iterations providing the lowest IVT and ∆Vol values.Correponding ice thikness RMSE (w.r.t observations, Bamber et al. (2013) are also indicated.Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2017-322Manuscript under review for journal Geosci.Model Dev. Discussion started: 8 February 2018 c Author(s) 2018.CC BY 4.0 License.