A rapidly converging initialisation method to simulate the present-day Greenland ice sheet using the GRISLI ice sheet model ( version 1 . 3 )

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 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 to infer initial conditions for ice sheet models that best reproduce available data. Most often such approaches allow for a good representation of the mean presentday state of the ice sheet but are accompanied with unphysical trends. Here, we present an initialisation method for the Greenland ice sheet using the thermo-mechanical hybrid GRISLI (GRenoble Ice Shelf and Land Ice) 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 a way that the simulated ice thickness matches the observed one. The quality of the method is assessed by computing the root mean square errors in ice thickness changes. Because the method is based on an adjustment of the sliding velocities only, the results are discussed in terms of varying ice flow enhancement factors that control the deformation rates. We show that this factor has a strong impact on the minimisation of ice thickness errors and has to be chosen as a function of the internal thermal state of the ice sheet (e.g. a low enhancement factor for a warm ice sheet). While the method performance slightly increases with the duration of the minimisation procedure, an ice thickness root mean square error (RMSE) of 50.3 m is obtained in only 1320 model years. This highlights a rapid convergence and demonstrates that the method can be used for computationally expensive ice sheet models.

velocities) and observations (e.g., Arthern and Gudmundsson, 2010;Morlighem et al., 2010). However, this approach may lead to internal inconsistencies between the simulated internal conditions (temperature and velocities) or between the simulated ice velocities and the observational datasets, such as surface and bedrock topography. These inconsistencies may have a strong impact on results of forward simulations. To circumvent this drawback, other authors (e.g., Perego et al., 2014;Mosbeux et al., 2016) developed a multi-parameter inversion technique to optimise both the sliding velocities and 5 the bedrock topography in such a way that the modelled surface ice velocities match with the observed ones. This allows the set of initial conditions to be self-consistent. However, if not constrained by observed ice thickness, these methods may lead to unrealistic simulated topography. An alternative approach, which avoids the previously mentioned shortcomings, consists in considering only the observed ice sheet geometry as the final target by finding appropriate basal conditions (generally the basal drag coefficient, see Sect. 2) that minimise the differences between observed and simulated ice thickness (Pollard and DeConto,10 2012; Pattyn, 2017). However, methods that choose to invert the basal drag coefficient only are not able to correct ice thickness errors in regions where there is no sliding (i.e. where bed is frozen). Moreover, while inverse methods are designed to produce an ice sheet state close to observations, the inferred basal drag coefficient may cancel errors coming from erroneous simulated basal temperatures and/or model physics shortcomings. Yet, as outlined by Pollard and DeConto (2012), the risk of cancelling errors is of lesser importance compared to those related to inconsistencies between internal conditions and surface properties 15 that will likely to be considerably reduced with expected future improvements in ice sheet models and better observations of basal conditions.
Here, we present a new iterative initialisation procedure that relies on the same basic principles as those developed by Pollard and DeConto (2012) (referred to as PDC12 in the following) and applied by Pattyn (2017) for the Antarctic ice sheet using 20 linear and non-linear sliding laws. 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 iterative method ensures a more rapid convergence and is thus suitable for computationally expensive models. 25 The paper is organised as follows. In Sect. 2 we present the main characteristics of the ice sheet model used in this study.
Section 3 describes the iterative minimisation procedure in detail. The main results are presented in Sect. 4 and sensitivity experiments in Sect. 5. These sections are followed by a discussion and the main conclusions of the present study (Sect. 6).
2 The ice sheet model GRISLI 30 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 . The version used in this study has been specifically developed for Greenland (Quiquet et al., 2012) with a horizontal resolution of 5 km x 5 km (301 x 561 grid points) and 21 evenly spaced vertical levels. GRISLI accounts 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, ice dynamics 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 G is the depth-averaged velocity (2D vector), SM B 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 10 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 heuristic sum of the SSA and the SIA components, as in Bueler and Brown (2009) but with no-weighting function (Winkelmann et al., 2011). In this case, the SSA velocity is used as the sliding velocity. We assume no-slip conditions for a 15 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 viscous till with a uniform thickness, 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. 20 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 (Ef) in the Glen's flow law to artificially account for the impact of ice anisotropy on the deformation rate. This enhancement factor depends on the stress regime (e.g., Huybrechts, 1990). Lower enhancement factors lead to lower deformation rates and as such to slower ice velocities. The 25 grounding line position is defined according to a flotation criterion and floating points are treated following the SSA only.
Calving physics is not explicitly computed, but if a grid point at the ice-shelf front fails at maintaining a thickness threshold, it is automatically calved . The ice thickness cut-off threshold is set to 250 m.
Since GRISLI is thermo-mechanically coupled, the ice temperature influences the ice velocity via the viscosity. The temper-30 ature is computed both in the ice and in the bedrock by solving a time-dependent heat equation. The temperature signal itself depends on ice deformation, surface temperature forcing and geothermal heat flux.
The basic principle of inverse modelling approaches for ice sheet initialisation procedure is to adjust the basal drag coefficient (β) which varies spatially, in order to reduce the mismatch between the simulated surface ice velocities and/or the ice sheet geometry and the observed ones.
While numerous studies are based on fitting the modelled ice velocities (e.g., Gudmundsson and Raymond, 2008;Arthern 5 and Gudmundsson, 2010;Morlighem et al., 2010;Perego et al., 2014), or both surface velocities and basal topography (Perego et al., 2014;Mosbeux et al., 2016), only few authors opted for fitting ice surface elevation (Pollard and DeConto, 2012;Pattyn, 2017). Here, we decided to adjust the basal sliding velocities via the adjustment of the β coefficient to fit the GrIS ice thickness to the observed one. Similarly to Perego et al. (2014), our choice is motivated by the need to refine the estimates of GrIS contribution to future sea-level rise without the sea-level rise signal being contaminated by 10 unphysical transients from the initial condition. However, while Perego et al. (2014) adopted a formal minimisation approach (i.e. adjoint-based model) we suggest instead an ad hoc method potentially applicable to any ice sheet model.
The GRISLI climate forcing, i.e. surface mass balance and surface air temperature ( Fig. 1), is provided by the regional atmospheric model MAR (Fettweis et al., 2013) forced at its boundary by the ERA-Interim reanalyses (Berrisford et al., 2011).

15
Both forcing fields are averaged over the 1979-2005 period (Figs. 1a and. 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). For the geothermal heat flux we use the data generated for the SEARISE project (Fox Maule et al., 2005). Initial geometry consists in the present-day observed ice thickness and bedrock elevation taken from Bamber et al. (2013). To compute initial conditions consistent with the boundary conditions, we run a 30 000 year-integration of the model imposing a fixed 20 topography. For this long experiment, similar to the fixed topography spin-up method, we assumed a perpetual present-day climate forcing (Figs. 1a and b) and we used a basal drag coefficient ( Fig. 2a) coming from a previous simulation carried out within the Ice2Sea project (Edwards et al., 2014). The resulting basal temperature after this long integration, presented as a difference with respect to the pressure melting point, is shown in Fig. 1c. It shows areas with temperature largely below the pressure melting point, associated with frozen bed, and areas with temperature at the pressure melting point (red colours), asso-25 ciated with thawed bed. Compared to the recent synthesis of GrIS basal temperatures (see Fig. 11 in MacGregor et al. (2016)), our initial basal temperature agrees generally well with the reconstructions in the north-western and north-eastern parts of the GrIS but are probably overestimated, with a too large thawed bed area, in the eastern and central parts of the GrIS (not shown).
The impact of ice temperature on the minimisation procedure is discussed in Sect. 5.1.

30
In order to avoid inconsistencies between the different datasets used as boundary and initial conditions, GRISLI is first run forward (free-evolving surface elevation and temperature) for 5 years (relaxation step, blue box in Fig. 3). After this short relaxation period, we start the iterative minimisation procedure (red box in Fig. 3). 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. Instead of optimizing the basal drag coefficient every 5 000 years as in PDC12, here the optimization is done at every time step (which is set to one year for the present study), using an ice thickness ratio to correct the simulated sliding velocity with the help of a modification of the basal drag coefficient.
The iterative minimisation procedure itself consists in repeated cycles, each cycle being divided in two main steps (red box 5 in Fig. 3): 1 st step: The first step consists in a free-evolving simulation (thickness and temperature) during which we adjust, at each model time step, the basal drag coefficient so that the ice thickness difference with respect to the observations becomes minimal. To this end, from the simulated vertically-averaged velocity (U G ) computed from the previous time step (or from the values obtained after the relaxation for the first iteration), we compute a corrected vertically-averaged velocity 10 field (U corr ) as a function of the computed (H G ) and observed ice thickness (H obs ): As seen before (in Sect. 2), the mean velocity field U G is the sum of two velocity components: the sliding velocity U sli and the vertically-averaged velocity U def due to vertical ice deformation: Assuming that the differences between U G , the simulated vertically-averaged velocity field, and U corr , the idealised vertically-averaged velocity field, are only due to changes in the sliding velocity U sli , we can write:

20
Following Eqs. 4 and 5, we can deduce the corrected sliding velocity (U sli corr ): U sli corr represents the corrected sliding velocity whose difference with U sli indicates how the simulated sliding velocity must change to reduce the mismatch between H G and H obs .

25
As such, we use the ratio between the simulated and the corrected sliding velocities ( U sli U sli corr ) to compute a new basal drag coefficient (β new ). This results in slowing down or speeding up the simulated sliding velocity and acts to reduce the gap between H G and H obs : Equation 7 is essentially identical to what is done in Price et al. (2011) except that they use observed and modelled velocities rather than observed and modelled ice thickness to adjust the basal drag coefficient. It should be noted that U sli corr can be lower or equal to 0, leading to infinite or negative basal drag coefficient. This can happen when the velocity due to vertical shearing U def is greater or equal to U corr . In this case we artificially impose a no-slip condition by assigning to the basal drag coefficient a maximum value set to 5 10 5 Pa yr m −1 . On the other hand, in case of too small U def 5 velocity, β may be as low as 1 Pa yr m −1 to facilitate ice sliding. Owing to its design, the method is only able to correct for the ice thickness mismatch where sliding occurs, i.e. where the base of the ice sheet is at the pressure melting point.
Throughout this step, the basal drag coefficient is updated at each time step for each model grid point. In the following, the duration of this step is referred to as Nb inv and has typical value of a few decades. 10 Note that, using Eq. 3 and Eq. 4, we can show that Eq. 7 can be rewritten as: As such, the adjustment of the basal drag coefficient is stronger in regions dominated by ice deformation.
2 nd step: The second step consists in running a new free-evolving simulation but this time using a time constant (but spatially 15 varying) basal drag coefficient, i.e. the last inferred basal drag coefficient of the first step. The duration of this second step, referred to as Nb f ree in the following, is generally longer than that of the first step, typically a few decades to a few centuries. This step aims at quantifying the model drift and the model mismatch with observations for the inferred basal drag coefficient. The simulated ice sheet velocity and topography at the end of this second step are used to compute a new U corr value in order to start a new cycle from the first step. The number of iterative cycles will be noted Nb cycle in 20 the following.
In summary, our iterative minimisation procedure consists in: (i) Adjustment of the basal drag coefficient at each time step (each year) year for Nb inv years (1 st step, Eqs 3 to 7).
(ii) Free-evolving simulation with the last inferred basal drag coefficient from (i) for Nb f ree years (2 nd step).
(iii) The steps (i) and (ii) are repeated Nb cycle times.

25
In addition, to assess the performance of the minimisation procedure (i.e. the quality of the inferred basal drag coefficient), we compute some quality metrics at the end of each cycle (green box in Fig. 3). The metrics are computed at the year 200 of the free-evolving simulation of the second step, independently from its duration (i.e. Nb f ree ). If Nb f ree is shorter than 200 years, we simply extend the simulation until 200 years. The quality metrics discussed in Sec. 4 include in particular the root 30 mean square error (RMSE) of the simulated ice thickness with respect to the observations and the drift in geometry (integrated ice thickness changes). These metrics help to decide whether an additional cycle is required or not. In the following, we also discuss the spatial patterns of ice thickness and ice velocity mismatches with respect to observations. Our method does not use the observed surface velocity as a constraint. However, at the end of the minimisation procedure (e.g., minimal thickness error and minimal drift), the simulated velocity tends nonetheless to approximate the balance velocity, that is the depth-averaged velocity required to maintain the steady-state of the ice sheet.
Once the optimal basal drag coefficient is found, it can be used to run prognostic forward simulations such as in Le clec'h The simulated GrIS ice volume obtained for this experiment is 1.4 % higher than the one estimated by Bamber et al. (2013) from observations (2.71 10 6 Gt). This overestimation is driven by large positive ice thickness differences (> 200 m) with re-15 spect to observations in the margin regions (Fig. 4a). There are also negative ice thickness differences in the interior of the ice sheet, in particular in the central eastern region. On top of this geometry mismatch, this experiment also presents a drift at the end of the 200 years with a negative contribution to global sea level of 0.7 mm yr −1 (i.e. 263 Gt yr −1 ice mass gain). Compared to observations (Joughin et al., 2018), the simulated ice velocity presents the same large-scale pattern but with important local differences (Fig. 4b). In particular, the main GrIS glaciers are generally too slow. 20 These results show the limitations of the simulated GrIS under constant climate forcing without appropriate initialisation procedure. In this specific case, the simulated model drift can potentially counterbalance the effect of climate warming expected in the future, leading to unrealistic projected Greenland melting contribution to global sea level rise. Therefore, the use of an initialisation procedure to minimise the model drift with a realistic simulated topography is not avoidable if the goal is to 25 produce reliable sea-level projections.

Iterative minimisation performance for a range of enhancement factor values
An increase (respectively a decrease) of the basal drag coefficient (β) slows down (resp. speeds up) the sliding velocity and thus the ice flow. Based on the adjustment of the sliding velocity, our iterative minimisation procedure, allows for a tuning of β only in regions where the basal temperature is at the pressure melting point, i.e. where the ice can slide over the bedrock.
no sliding occurs. In order to slow down or speed up the ice flow in such regions, the value of the enhancement factor, Ef (see Sec. 2), can be tuned. As explained in Sect. 2, this factor is used to increase (when > 1) or decrease (when < 1) the ice deformation velocity. The more the ice deformation is increased (respectively decreased), the more the ice flow in frozen base region speed-ups (resp. slow-downs) and thus decrease (resp. increase) the ice thickness.
The enhancement factor for the SIA regime (slow ice flow) is expected to have a large influence on shear-stress driven ve-5 locities (Quiquet et al., 2018). Generally set to 3 (Ritz et al., 2001), the Ef can be chosen within a large range of values between 1 to 10 (Ma et al., 2010). In the following, we assess our iterative minimisation procedure for a range of Ef values: 0.1, 0.5, 1, 1. 5, 2, 2.5, 4, 3, 3.5, 4, 4.5, 5 The ice thickness RMSE defined with respect to observations is displayed in Fig. 5 as a function of the number of cycles performed for the different enhancement factors. For a given Ef value, the RMSE quickly decreases during the first cycles and generally stabilises after Nb cycle ≈ 5-6. This means that the procedure is very effective in reducing the ice thickness error for the first iterations but does not entirely correct the mismatch with observations. Depending on the enhancement factor considered, the overall improvement represents a reduction of about 20 to 40 m in ice thickness RMSE with respect to the first 20 iterative cycle.
The RMSE is largely different for the different enhancement factors. For Ef 2, we have systematically a larger RMSE for a larger Ef value regardless of the number of iterative cycles performed. This is no longer the case for smaller Ef since the experiment (i.e. the Ef value) providing the lowest RMSE is different for the Nb cycle considered. Note that for Ef = 0.5 the RMSE value is often larger than that obtained with Ef = 2.5 even with increasing Nb cycle . Indeed, Ef = 0.5 implies a too small 25 deformation rate that leads to a too slow ice flow velocity. As such, the departure from the observations is mainly characterized by positive ice thickness anomalies at the edges and in the half southern part of the ice sheet. In addition to the RMSE criterion, which is an integrated metric, the maps of the difference between the simulated and the observed ice thickness bring valuable information to understand the model structural biases. In Fig. 6, we can distinguish two 5 main patterns. Except for Ef = 0.5, all the Ef experiments with a Nb cycle producing the minimum RMSE value (Fig. 6 and Table 1) are marked with an underestimation in ice thickness in the interior and an overestimation at the edges of the GrIS.
This overestimation can be slightly reduced using higher Ef values, the underestimation being nonetheless larger in this case.
As explained above, larger Ef values amplify ice deformation and therefore speed up the ice velocity, explaining the spread 10 of the regions where the ice thickness is underestimated (Fig. 6). Some of these regions, such as a in some cases, the velocity due to deformation is too fast and the basal drag coefficient is set to its maximal value (β max = 5 10 5 Pa yr m −1 ) so that the sliding velocity becomes virtually zero. This is visible in Fig. 7 where the area for which the basal drag coefficient is set to β max (dark red colour) is getting larger with increasing Ef. On the contrary, in Fig. 7, where the model overestimates the ice thickness (i.e. too slow ice flow) and where basal temperature is at the pressure melting point (i.e. sliding can occurs), the computed basal drag coefficient is weaker in order to increase basal sliding. Similar to the β max region, 20 our iterative initialisation method could reach a minimum basal drag coefficient value (set to β min = 1 Pa yr m −1 ) in regions where the sliding velocity must be as strong as allowed by the flow law equation (i.e. meaning no basal friction). Reducing the enhancement factor, and thus the ice deformation in these regions can locally increase the ice thickness overestimation.
Regions with β max or β min values are an indication of the limit of our iterative ice thickness error minimisation procedure.

25
The ice thickness errors shown in Fig. 6 correspond to a median value ranging from +15 m to -99 m from the lowest (Ef = 0.5) to the highest (Ef = 5) enhancement factor. The decrease in the median of the error with increasing Ef values is mostly driven by the underestimation of the ice thickness in the interior regions. Our results show that the Ef = 1 experiment produces the best ice thickness error pattern, ranging from +133 m (5 th quantile) to -39 m (95 th quantile) and reaching a median error equal to +3 m.

b) Total ice volume and compensating biases
Because most of the Ef experiments have both positive ice thickness biases at the margins and negative biases over the central part (Fig. 6), the global ice sheet volume is not a good metric for model performance due to compensating biases. Fig-ure 8 shows the total ice volume difference with respect to observations for varying enhancement factors as a function of the number of iterative cycle. Some specific experiments show a very small error in global ice volume with respect to observations for given Ef values even though they have a poor RMSE (Fig. 5). Also, for Nb cycle = 6, RMSE value for Ef = 0.5 and Ef = 2.5 are identical (68.8 m) but ice volume anomalies are drastically different with 57994 Gt and -38841 Gt respectively (Table 1). Thus, a small global ice sheet volume difference does not necessary mean a minimisation of the ice thickness difference. 5 For the same reasons, the trend in global ice volume is not a good metric for assessing the ice sheet drift because local changes in ice thickness can compensate for each other. As an illustration, using a range of enhancement factors, Fig. 9 shows for free-evolving simulations, the temporal evolution of the total ice volume difference with respect to observations, along with the evolution of the RMSE. This figure confirms that the GrIS volume equilibrium can be reached by biases compensation as 10 we have a near-zero error in volume with Ef = 1.5 while the RMSE is very similar to that obtained with Ef = 1 and Ef = 2. For Ef 2, the negative biases in ice thickness dominate, with a decrease in ice volume as Ef increases. For Ef < 2, the positive biases in ice thickness dominate, leading to an increase of the global ice volume.
To assess the simulated ice sheet drift and in order to avoid the biases compensation, we compute the geometry trend as the 15 root mean square ice thickness change (ξ t , expressed in cm yr −1 ): where < (H t − H t−1 ) 2 > represents the averaged squared ice thickness change over the whole GrIS.
Values of ξ computed from the last 5 years of the 200 yr free-evolving simulation in the 2 nd step (green box in Fig. 3) are 20 reported in Table 1 for a given iteration and varying enhancement factors. The lowest values are generally obtained with the experiments that provide the lowest RMSE, which means that these simulated ice sheets are the closest to equilibrium. The minimal trends are about 15 cm yr −1 and are obtained with enhancement factors between 1 and 2.

c) Ice dynamics 25
Our iterative minimisation procedure aims at simulating an ice thickness as close as possible to observations. Hence, the observed ice velocity is not used as a target by the model. However, because our procedure generates an ice sheet at quasiequilibrium (trend ξ close to 0), the simulated velocities are close to the balance velocities, which in turn are supposedly close to present-day observations. As a result, our method simulates an ice flow pattern similar to the observations (Fig. 10).

30
The simulated velocity field is particularly sensitive to the choice of the enhancement factor (Fig. 10). In particular, for the highest Ef values, the simulated velocity is overestimated for the major ice streams where deformation due to vertical shearing is expected to be of lesser importance compared to basal sliding. For Ef = 1.5, the ice flow pattern in the margin regions is well reproduced compared to observations. Only some glaciers ice velocities can be faster (e.g. Jakobshavn or Kangerlussuaq) or slower (e.g. Petermann or NEGIS). While the best GrIS geometry (lowest RMSE) is obtained with Ef = 1.5, the experiments with Ef = 1 or Ef = 0.5 best reproduce the observed surface velocities (RMSE about 150 m yr −1 , Fig. S1).
Interestingly the extent of the NEGIS is particularly well represented, in particular for lower enhancement factors (Fig.   5 S2). This can be a relic of the long temperature equilibrium performed with a time constant basal drag coefficient taken from Ice2Sea experiments (Edwards et al., 2014), in which the NEGIS is well delimited (Fig. 2a). However, because this feature is still present when starting the iterations from a spatially homogeneous basal drag coefficient (see Sec. 5.2), it can also suggests that there is some topographic control of this feature as the adjustment of our local basal drag coefficient is very effective in reproducing the observed velocity in this area. Having a good representation of the NEGIS is an encouraging sign for the 10 performance of our minimisation procedure, especially since most models fail to achieve this (Goelzer et al., 2018). where the bed is frozen our iterative mininisation procedure is unable to correct for the ice thickness mismatch. This leads to a predominant role of the enhancement factor. The aim of this section is to investigate the sensitivity of our procedure to the initial temperature profile. To this end, we followed the same methodology as in Sec  (Fig. 12). In the following, the temperature profile taken from Edwards et al. (2014) is referred to as the non-equilibrated temperature as opposed to the 30-kyr equilibrated temperature used in the rest of the manuscript.
25 Figure 13 shows the evolution of the RMSE for 9 iterative cycles for the experiment performed with the non-equilibrated temperature profile with Ef = 3 (dark blue dots). Similarly to what was shown in Sec. 4.2, the minimisation procedure reduces the RMSE from +76.0 m after Nb cycle = 1 to a minimum after Nb cycle = 9 around +47 m. Figure 13 shows also the evolution of the RMSE for two experiments with Ef = 1 and Ef = 3 but using the equilibrated temperature profile (cyan and orange dots 30 in Fig. 13). If the pattern is essentially the same between the different experiments, the RMSE is higher when using the equilibrated temperature. For the same Ef value, the RMSE is 11.4 m higher (Nb cycle = 8) when using the equilibrated temperature.
This is because the warmer equilibrated temperature with respect to the non-equilibrated one leads to higher velocities which ultimately favour the ice thickness underestimation in the central regions (shown in Fig. 6). Using a smaller enhancement factor with the equilibrated temperature reduces the gap (3.3 m for Nb cycle = 8) and provides a closer response to that obtained for Ef = 3 with the non-equilibrated temperature.
If the RMSE is lower when the non-equilibrated temperature profile is used, the trend ξ is nonetheless largely higher (24.7 5 cm yr −1 for Nb cycle = 6) compared to the experiments with an equilibrated temperature (16.5 cm yr-1 for Ef=3 and 16.3 cm yr −1 for Ef = 1, for Nb cycle = 6 and 8 respectively). This is expected as there is an important thermal adjustment when using a profile that is not consistent with the climatic forcing.
However, despite existing differences with results obtained with the equilibrated temperature profile, this shows that our 10 minimisation procedure is able to reduce the mismatch between simulated and observed ice thickness independently from the initial temperature profile.

Sensitivity to the initial basal drag coefficient
As explained in Sec. 3, the initial basal drag coefficient β for the first iteration of the minimisation procedure is the one used Using Nb inv = 20, Nb f ree = 200, and Nb cycle varying from 1 to 15 with Ef = 1, we obtain a minimum ice thickness RMSE of 49.9 m and a trend ξ of 15.1 cm yr −1 . While there are some minor spatial differences in terms of the inferred basal drag coefficient (Fig. 2c), the aggregated metric such as the RMSE and the trend are identical to the results presented in Table 1. In 20 the same way, the simulated ice thickness and surface velocities obtained with β = 1 present very small differences with those obtained when starting from the Ice2Sea basal drag coefficient (Figs S3 and S4). This illustrates the robustness of the method and shows that it does not depend on the chosen initial distribution of the basal drag coefficient.

Sensitivity to the duration (Nb inv and Nb f ree ) of the minimisation procedure
In this section we assess the sensitivity of the minimisation procedure to the coefficients Nb inv (i.e. the duration of the period 25 during which the basal drag coefficient is iteratively computed -1 st step) and Nb f ree (duration of the free-evolving simulations -2 nd step). While in Sec. 4.2 we used Nb inv = 20 and Nb f ree = 200, here we explore a range of combinations of these parameters exploring four values for Nb inv (20, 40, 80, 160 yr) and Nb f ree (50, 100, 200 and 400 yr). Using an enhancement factor of 1, we iterate on 15 cycles (Nb cycle from 1 to 15). The initial conditions are the same as in Sec. 4.2.
30 Figure 14 shows the evolution of the RMSE as a function of the number of cycles performed for a range of Nb f ree values.
As previously shown, there is a strong decrease in RMSE between the first two cycles and only a limited improvement when using more than 6 cycles. The response is very linear: using larger Nb f ree leads to a smaller RMSE. This can be explained by the fact that the correction computed at the end of the 2 nd step, after Nb f ree , is greater if the duration of the free-evolving simulation is longer. This means that the changes imposed to the new basal drag coefficient computation (Eq. 7) from one cycle to another is larger for longer Nb f ree .
In Fig. 15 we show the evolution of the RMSE as a function of the number of cycles performed for a range of Nb inv values 5 (20,40,80 and 160 yr). The RMSE difference for a given Nb cycle is generally less than 10 m, while this difference is sometimes larger than 20 m when Nb f ree varies from one value to the other. This suggests that Nb inv is of second importance relative to Nb f ree . The RMSE appears to be slightly smaller for longer Nb inv . For example, for Nb f ree = 200 yr, increasing Nb inv from 20 yr to 40, 80 or 160 yr slightly reduces the minimum RMSE by 0.1 m, 1.7 m or 3.5 m respectively, and decreases the trend ξ by 13.7%, 7.2% and 21.1% for Nb cycle equal to 12, 11 and 8 respectively. The minimum RMSE value (46.1 m) and trend ξ 10 (12.3 cm yr-1) are reached with Nb cycle = 10 and with Nb inv = 160 yr. Performing more cycles once the minimum RMSE is reached does not improve the results.
Overall, the combination of the highest Nb inv (160 yr) with the highest Nb f ree (400 yr) leads to the smallest RMSE (44.1 m) with a trend ξ of 9.9 cm yr −1 for Nb cycle = 11. However, this minimum represents a considerable amount of computing 15 time (6160 years) and does not represent the most efficient combination. As shown in Figures 5, 8 and 13, the minimum RMSE generally stabilises between Nb cycle equal to 4 and 6. This means that similar RMSE and trend ξ could be obtained using fewer computing resources. For each set of combination, the mean value of the best RMSE values is equal to 51.1 m and is associated to a mean trend ξ of 15.5 cm yr −1 . The experiment with Nb inv = 20 yr, Nb f ree = 200 yr and Nb cycle = 6 produces a RMSE 0.6 m lower than the mean and is more than three times faster than best of the RMSE (1320 years to be compared to 6160 years).

Summary and discussion
In order to improve the reliability of Greenland ice sheet simulations in a future transient climate, an accurate evaluation of the present-day trend of ice flow dynamics is required. One the major difficulties in addressing this need lies in the poorly constrained observational data of the basal conditions that strongly control the ice motion in the entire ice sheet. Here, we present an inverse method to infer the spatial distribution of the basal drag coefficient in such a way that the mismatch between 25 simulated and observed GrIS ice thickness is minimised. As such, our target criteria are defined for the sets of minimisation procedure parameters providing minimum values of ice thickness RMSE (with respect to observations) and ice thickness trend, which are respectively as low as ∼ 50 m and 15 cm yr −1 for our best fit. This remains in the range of PDC12 results. The great advantage of the method is its rapid convergence (i.e. 1320 years) making it suitable for more computationally expensive models. Moreover, we have also shown that it only poorly depends on the initial guess of the spatial distribution of the basal Based on the adjustment of the basal sliding, our method cannot be applied in regions of frozen bed and is only effective in thawed bed areas where basal sliding may occur. However, in case of a too large deformation rate in these regions, the basal drag coefficient is set to its maximum value to counteract the too fast ice flow. The limit of applicability of the method led us to investigate the impact of the enhancement factor which is expected to have a large influence on the deformation rate, and thus, on the ice flow and subsequently on the simulated ice thickness. We performed a series of simulations with a range of various 5 values of the enhancement factor (from 0.5 to 5) and showed that the mismatch between the simulated and the observed GrIS topography is reduced with an appropriate tuning of the enhancement factor. This highlights that the overall performance of the method is critically dependent on the basal thermal state and points out that the finding of appropriate initial conditions with a simple adjustment procedure remains an undetermined issue. Actually, multiple combinations of the enhancement factor and the basal drag coefficient can produce a simulated ice thickness close the observed one, but this cannot discard the possibility Finally, we have shown in this paper that the iterative adjustment of β produces modelled surface velocities that compare well with the observed ones. This suggests that future work could include an additional metric related to surface ice velocities 20 so as to further reduce the uncertainties associated with the choice of model parameters and variables.
Another limitation of the method may come from the model resolution. The succession of higher/lower ice thickness due to the succession of valleys/ridges in mountain areas may be poorly resolved. Owing to the insulation effect of the ice, this may lead to an erroneous representation of the basal temperature patterns, and SSA regions may be erroneously interpreted as 25 frozen bed regions and vice versa (Pattyn, 2010). This drawback is clearly illustrated in our study in Fig. 6 (Ef = 1). Indeed, the simulated ice thickness obtained with the inversion procedure is generally less than 50 m in most GrIS areas, but can be greater than several hundred meters in coastal mountain ranges such the central eastern margin area where ice flow occurs in deep valleys. An alternative solution consists in correcting the basal temperature to account for bedrock roughness, similarly to what was done in PDC12 to improve their inversion procedure in the Transantarctics. On the other hand, higher resolution 30 models can also better account for the dynamics of small-scale outlet glaciers and for their interactions with floating ice that strongly influence the ice sheet mass balance (e.g., Aschwanden et al., 2016). However, due to the elliptic character of the SSA equation (e.g., Quiquet et al., 2018), the local adjustment of the basal drag coefficient impact the ice velocity of neighbouring points. As a result increased resolution may increase the noise, unless introducing a smoothing function that filters the high frequency noise (Pattyn, 2017).
The reliability of the method also depends on the quality of observations data and of climate forcing. Errors in observed surface or bedrock topography or in SMB patterns different from those associated to the observed ice thickness 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 5 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 hydrological 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 shortterm projections, it is more questionable at the century time scale, and future modelling efforts should therefore be undertaken 10 to compute interactively the basal drag coefficient as a function of changes in basal conditions.     Ef=1.5 Ef=2 Ef=2.5

T H 5 yrs
Ef=3.5 Ef=4 Ef=4.5 Ef=3 Ef=5 Ef=0.5 Figure 6. Difference between the simulated and the observed (Bamber et al., 2013) ice thickness (in meters) for Ef ranging from 0.5 to 5 for the iterative cycle Nb cycle that produces the lowest RMSE (Table 1). Here, Nbinv and Nb f ree are set to 20 and 200 years respectively.

Ef=1
Ef=1.5 Ef=2 Ef=2.5 Ef=3.5 Ef=4 Ef=4.5 Ef=3 Ef=5 Ef=0.5 0 -1000 -500 -200 -100 -10 1000 500 200 100 10 Ice velocity di erence (Log 10 m yr -1 ) Figure 11. Simulated ice surface velocity difference (in m yr −1 ) with respect to observations (Joughin et al., 2018) using Ef ranging from 0.5 to 5 for Nbinv = 20 yr, Nb f ree = 200 yr and Nb cycle that corresponds to the one producing the lowest ice thickness RMSE (see Table 1).  yr as a function of the number of iterations (Nbinv). Dark blue dots are for the experiment that uses the non-equilibrated temperature profile as initial condition and Ef = 3. Cyan and orange dots are for the experiments using the equilibrated temperature and Ef = 3 and Ef = 1, respectively.