2.2.1 Model description

The GRISLI (GRenoble Ice Shelf and Land Ice) ice sheet model was first developed to compute the dynamical evolution of the Antarctic ice sheet (Ritz et al., 2001; Philippon et al., 2006; Alvarez-Solas et al., 2011a; Quiquet et al., 2018). It has then been successfully applied to the Northern Hemisphere ice sheets (Peyaud et al., 2007; Alvarez-Solas et al., 2011b; Charbit et al., 2013) and the Greenland ice sheet (Quiquet et al., 2012, 2013). In the present work, we use a 5 km resolution grid covering Greenland (301 × 561 grid points) and 21 evenly spaced vertical levels in the ice and 4 levels in the bedrock. The regular grid is projected on a polar stereographic grid with a standard parallel at 71^∘ N and a central meridian at 39^∘ W (same as in Bamber et al., 2013). GRISLI is a 3-D thermomechanically coupled ISM computing the temporal evolution of the ice sheet, which is a function of surface mass balance, ice flow and basal melting:

where t is time, H the ice thickness, U^G the vertically averaged velocity, ∇(U^GH) the ice flux divergence, SMB the surface mass balance and b_melt the basal melting. Basal melting occurs when the basal temperature is at the pressure melting point. The ice temperature plays a crucial role in the dynamics of the ice sheet because it affects the viscosity, and thus the ice flow in the entire ice column (Ritz et al., 1996, 2001). In turn, heat released by internal ice deformation and basal dragging over the bedrock modifies the temperature. The temperature field is computed by solving a time-dependent heat equation both in the ice and in the bedrock, accounting for advection and vertical diffusion processes. At the surface, the boundary condition is provided by the prescribed surface temperature. At the base of the ice sheet, the boundary condition is given either by the geothermal heat flux or by the temperature melting point at the ice-bed interface.

The ice flow is computed using both the shallow ice (Hutter, 1983) and shallow shelf (MacAyeal, 1989) approximations to solve the Stokes equations (Ritz et al., 2001). The shallow ice approximation (SIA) assumes that ice flow is caused only by vertical shear stress, neglecting the longitudinal stresses. This assumption is only valid for slow-flowing ice. For fast-flowing regions, vertical shearing becomes smaller than longitudinal shearing and the shallow-shelf approximation (SSA), which neglects the vertical stresses, is used. The ice thickness and the ice sheet surface slopes control the SIA and the SSA velocity components, but the SSA is also governed by basal dragging. Using a hybrid model (i.e. based on both SIA and SSA approximations) allows for better representation of the different deformation regimes found in an ice sheet. In GRISLI, the SSA velocity is used as a sliding velocity (Bueler and Brown, 2009) when the basal temperature is at the pressure melting point. In this case, we assume here a power-law basal friction (Weertman, 1957) and the presence of sediments allowing for viscous deformation. The relationship between the basal shear stress (τ_b) and the basal velocity (u_b) is expressed as follows:

where β is a time constant but spatially variable basal drag coefficient. For cold base conditions, the sliding velocity is set to zero.

The resulting velocity for every model grid point is the addition of the SIA and SSA components. For floating ice points (ice shelves), we assume no basal drag. In addition, if the ice thickness of the floating ice shelves is below 250 m and if no neighbouring points are grounded, the point is removed and the corresponding ice mass loss is considered a calving flux. Determination of the grounding line position is based on a floatation criterion.

The isostatic adjustment in response to ice loading changes is governed by the relaxation of the asthenosphere with a characteristic time constant of 3000 years and by the deformation of an elastic lithosphere (Le Meur and Huybrechts, 1996).

The climatic forcing is given by the mean annual SMB and the mean annual ST. Because seasonal variations of surface temperature are rapidly dampened, ST is considered as a good approximation of the bottom snowpack temperature. The initial GrIS surface and bedrock topographies come from Bamber et al. (2013) and the geothermal heat flux is taken from Fox Maule et al. (2005).

2.2.2 Initialization procedure

Due to the long timescale response of the ice sheet to a given climate forcing, a proper initialization of the model is required before performing forward experiments. For future sea-level projections, the aim of the initialization is to start the simulations from a present-day equilibrated ice sheet geometry as close as possible to the observed one while ensuring consistency between internal properties of the ice sheet (e.g. basal sliding velocities and vertical profile of temperature) with the climate forcing. Here we use an inverse method of the basal sliding velocities in order to reduce the mismatch between the simulated GRISLI ice thickness and the observed one (Bamber et al., 2013). The method is fully described in Le clec'h et al. (2018). Below we only remind the reader of the basic principles of our approach.

The procedure starts from initial conditions (i.e. vertical temperature and velocity profiles, first guess of the β coefficient) coming from previous GRISLI simulations carried out within the framework of the Ice2Sea project (Edwards et al., 2014a) and from the present-day observed topography (Bamber et al., 2013). The SMB and ST climatological means computed by MAR-MIROC5 (Fettweis et al., 2013) are used as climate forcing.

Our initialization procedure is based on an iterative process divided into two main steps.

1. The first step consists in the iterative adjustment of the spatially varying basal drag coefficient related to sliding velocities through Eq. (3). At each model time step, the ratio of modelled to observed ice thickness (i.e. the data–model mismatch) is used to adjust the sliding velocities via the basal drag coefficient, while the shearing velocities (SIA velocities) remain unchanged. This adjustment is typically performed for short time periods (a few decades).

2. After this adjustment phase, the model is left to freely evolve (second step) for anywhere from a few decades to a few centuries with the last inferred basal drag coefficient.

At the end of this second step, we obtain a new GrIS topography and, thus, a new data–model mismatch, which is used to start a new cycle in which the first and the second steps are repeated. This new cycle starts with exactly the same initial and boundary conditions as those used for the first cycle, except for the basal drag coefficient and the data–model ice thickness mismatch, which are inferred from the values computed at the previous cycle. The overall process is stopped when the ice thickness root mean square error is not significantly improved from one cycle to the other. This ensures a good compromise between the reduction of the mismatch between observed and simulated ice thickness and the rapidity of the convergence of the initialization procedure. In the present paper, the best fit with observations (RMSE = 63 m) is obtained for a first step duration of 20 years, a duration of 200 years for the relaxation GRISLI simulation (second step) and eight iterative cycles.

This method is based on the same basic principles as that of Pollard and DeConto (2012) except that their basal drag coefficient is adjusted as a function of the difference between modelled and observed ice surface elevation while we use the ice thickness ratio instead. Moreover, while the method suggested by Pollard and DeConto (2012) requires long (multi-millennial) integrations for the method to converge, we use an iterative method of short (decadal to centennial) integrations starting from the observed ice thickness allowing a more rapid convergence.

To further reduce the model drift in terms of ice volume and before starting the forward GRISLI experiments, a 2000-year GRISLI relaxation run is performed as a continuation of the second step of the last iteration (i.e. after the end of the eighth cycle). As such, the value of the basal drag coefficient is that obtained at the end of the first step of the eighth cycle. Over the last 150 years of this free-evolving simulation, the model drift is only $\pm {\mathrm{10}}^{-\mathrm{5}}$ mm yr⁻¹ sea-level equivalent. At the end of this 2000-year simulation, the simulated GrIS topography is slightly different from the observations (RMSE = 132 m; see Fig. S1). It will be referred to hereafter as S_ctrl and will be used as initial topography for the transient GRISLI simulations described in the following.

For all the simulations presented in this study, including those carried out within the iterative initialization framework, we apply a strong negative SMB value outside the observed ice sheet extent (Bamber et al., 2013). This avoids ice growth where there is none in reality and allows correction for both the potential atmospheric model biases (e.g. positive SMB values over tundra areas) and the initialization procedure biases (i.e. too strong ice export towards the margins). However, for GrIS projections run the RCP8.5 forcing scenario, this condition has only a limited impact on GrIS contribution to sea-level rise since the ice extent will likely to keep on retreating over the next few centuries.

4.1.1 Changes in the forcing climate

The evolution of the SMB and of its different components simulated in the 2W experiment (Eq. 1) and integrated over the entire GrIS is displayed in Fig. 1. During the 2000–2040 period, the averaged SMB remains positive with a mean value equal to 280±95 Gt yr⁻¹ (where the notation ± represents the standard deviation computed from yearly values) but slightly decreases by 4 Gt yr⁻¹. This decrease becomes substantially stronger from 2040–2100 (−17 Gt yr⁻¹ on average), and the mean SMB reaches strong negative values ($-\mathrm{638}\pm \mathrm{271}$ Gt yr⁻¹ over the 2090–2100 period). As the same MIROC5 year 2095 is repeatedly used to force MAR after 2100, there is no longer inter-annual variability and the integrated SMB remains quite stable between 2100 and 2150 ($-\mathrm{812}\pm \mathrm{13}$ Gt yr⁻¹). Indeed, MAR generates its own surface boundary layer fields which are not impacted by the MIROC5 forcing, explaining a slight SMB increase of ∼1 Gt yr⁻¹ over the last 50 years. Throughout the simulation, the evolution of the SMB signal is dominated by surface runoff whose increase rate (in absolute value) ranges from 5 Gt yr⁻¹ (2000–2040) to 19 Gt yr⁻¹ (2040–2100). After 2100, it slightly decreases (∼2 Gt yr⁻¹), explaining the slight SMB increase.

Figure 1Evolution of the SMB (black line) and its components from 2005 to 2150 (in m Gt yr⁻¹) simulated by MAR in the 2W experiment and integrated over the ice sheet mask taken from Bamber et al. (2013). The SMB components are snowfall (SF, red line), rainfall (RF, orange line), sublimation (-SU, light blue line) and runoff (-RU, dark blue line). Dashed lines correspond to regression lines ranging from 2000 to 2039 and from 2041 to 2099. The solid grey line corresponds to the zero line.

Download

The SMB anomaly between the beginning and the end of the 2W experiment is displayed in Fig. 2a. A total of 65 % of the grid points having surface elevations higher than 2000 m are characterized by a positive SMB anomaly, ranging from 0.07 m yr⁻¹ (5th percentile) to 0.2 m yr⁻¹ (95th percentile) at the end of the simulation. This SMB increase is particularly pronounced in the eastern part of the GrIS between 67 and 70^∘ N and in the northern-central part. It is due to a strong increase in snowfall (>0.5 m yr⁻¹, Fig. 2b), which occurs mainly during the winter season in the east and during autumn in the north (Fig. S2). On the other hand, 87 % of the GrIS grid points with surface elevation lower than 2000 m are dominated by an increase in surface runoff (Figs. 2c, S3 in the Supplement) and by an increase in the fraction of rainfall over snowfall in summer and in autumn (Fig. S4). As a result, strong negative SMB anomalies are found in these regions ranging from −3.3 m yr⁻¹ (5th percentile) to −0.1 m yr⁻¹ (95th percentile) and reaching more than −6 m yr⁻¹ along the western and the southeastern margins (Fig. 2a).

Figure 2Anomalies of (a) mean annual surface mass balance in m yr⁻¹, (b) annual snowfall in m yr⁻¹, (c) annual runoff in m yr⁻¹ and (d) mean annual surface temperature (in ^∘C). These anomalies are given between the last (2140–2150) and the first (2000–2010) decade of the 2W experiment. Panels (a) and (d) are computed on the GRISLI grid and (b) and (c) are given on the MAR grid. The dashed lines correspond to the 500 m surface elevation iso-contours for the present-day observed topography. The grey shading represents the non-ice-covered areas. Note for (a), (b) and (c) that the colour scale is not symmetrical for positive and negative values.

Download

The equilibrium line altitude (ELA, i.e. altitude for which SMB = 0) increases significantly between the beginning and the end of the 2W experiment, as a consequence of increased runoff for areas below 2000 m. As an example, at around 73.5^∘ N, on the eastern side of the ice sheet, the ELA moves from ∼1000 to ∼2500 m (Fig. 3). In other regions, at the end of the 2W experiment, the ELA is generally situated between 1500 and 2000 m high, except in the northern part where it is between 1000 and 1500 m. This shift of ELA towards higher altitudes represents an increase of 24 % of the ablation area between the beginning and the end of the experiment.

Figure 3GrIS surface elevation in 2150 simulated in the 2W experiment (in m). The solid black and purple lines represent the equilibrium line altitudes (ELA, limit between the accumulation and the ablation zones) in the 2000–2010 and 2090–2100 mean periods, respectively, for the NF, PF and 2W experiments. The yellow, red and orange solid lines indicate the ELA position over 2140–2150 for the NF, PF and 2W experiments, respectively. The dashed lines correspond to the 500 m surface elevation iso-contours for the present-day observed topography. The grey shading represents the non ice-covered areas.

Download

The ST anomaly (Fig. 2d) ranges from 2.2 ^∘C (5th percentile) to 6.5 ^∘C (95th percentile) and is characterized by a south–north gradient with the highest values found in the northern part. Beyond 78^∘ N the ST anomaly locally reaches values greater than 11 ^∘C. This temperature increase from 2000 to 2150 contributes to the amplification of the ablation processes below the ELA. However, while the stronger temperature anomaly is found in the northeastern part of the ice sheet, this region is also marked by the increasing snowfall in 2150 compared to 2000 (Fig. 2b), which counteracts the ablation processes.

Figure 4Ice thickness anomaly (2150–2000) simulated in the 2W experiment (in m). The dashed lines correspond to the 500 m surface elevation iso-contours for the present-day observed topography. The grey shading represents the non-ice-covered areas. A non-linear colour scale is used for positive values.

Download

4.1.2 Changes in Greenland ice sheet geometry

Figure 4 displays the ice thickness anomaly between the beginning and the end of the 2W experiment. For surface elevations higher than 2000 m in the northern part, and higher than 2500 m in the central and southern parts of the ice sheet, the ice thickness increases by 5 m on average, with the increase ranging from 1.5 m (5th percentile) to 17 m (95th percentile). On the other hand, in regions whose surface elevation is lower than 2000 m, the ice thickness decreases from −248 m (5th percentile) to −3 m (95th percentile) with a mean value equal to −100 m. As a result of these GrIS ice thickness changes, the surface slope between the central part of the ice sheet and the margins increases. On top of that, the ice sheet mask (defined as the fraction of a MAR grid cell with permanent ice cover, e.g. Fig. S5a) decreases by 2.8±0.1 % (mean ± standard deviation computed from yearly values) over the 2140–2150 mean period compared to the 2000–2010 mean period, and some GrIS margin regions become ice-free.

Figure 5Cumulated SMB (a) and ice flux divergence (b) throughout the entire 2W experiment (from 2000 to 2150) given in metres and computed on the GRISLI grid. The dashed lines correspond to the 500 m surface elevation iso-contours for the present-day observed topography. The grey shading represents the non-ice-covered areas. A non-linear colour scale is used for positive values.

Download

4.1.3 Impact of ice flow changes

The ice thickness anomaly is due to the complex combination of changes in surface atmospheric conditions (SMB, Fig. 5a), ice flow (ice flux divergence, Fig. 5b) and basal melting (not shown), following the continuity equation (Eq. 2). To quantify the role of ice flow on the GrIS geometry (Fig. 4), we plotted the ice flux divergence integrated over 150 years (2000–2150; see Fig. 5b). In particular, over the central plateau, the cumulated SMB (Fig. 5a) reaches about 50 m, 40 m of which is transported away by the ice flow (Fig. 5b). As a result, the ice thickness anomaly is reduced to only ∼10 m in this region (Fig. 4). An opposite behaviour is found near the western coast, where the runoff is partly compensated by ice convergence, resulting in a less negative ice thickness anomaly than that related to the SMB forcing. This shows that ice flow acts to counteract ice loss from surface melting, as previously noticed by several authors (Huybrechts and de Wolde, 1999; Goelzer et al., 2013; Edwards et al., 2014a). As a consequence, it appears to be essential to account for ice dynamics to accurately estimate the mass balance of the whole ice sheet.

Figure 6(a) Mean surface velocity anomaly (in m yr⁻¹) between the 2140–2150 and the 2000–2010 mean periods for the 2W experiment. (b) Mean surface velocity difference between the NF and the 2W experiments for the 2140–2150 mean period. (c) Same as (b) for the PF and 2W experiments. These differences are computed on the GRISLI grid. The dashed lines correspond to the 500 m surface elevation iso-contours for the present-day observed topography. The grey shading represents the non-ice-covered areas. A non-linear colour scale is used for both positive and negative values.

Download

In turn ice flow is impacted by changes in ice sheet geometry as illustrated by the mean surface velocity anomaly (Fig. 6a). For regions with surface altitudes between 2000 and 2500 m, the anomaly of the ice flow increases from the inner GrIS areas towards the edges of the ice sheet. The increase in the mean ice flow for the 2140–2150 period compared to 2000–2010 period ranges from 0.08 m yr⁻¹ (5th percentile) to 17 m yr⁻¹ (95th percentile). These faster ice velocities at the end of the 2W experiment are mainly explained by a larger surface slope between the central and the margin regions of the ice sheet. This is consistent with information inferred from ice flux divergence as shown in Fig. 5b.

On the contrary, for the margin regions with altitudes lower than 1500 m, the anomalies of surface ice velocities strongly decrease (Fig. 6a). Compared to the 2000–2010 period, this decrease ranges from −213 m yr⁻¹ (5th percentile) to −0.2 m yr⁻¹ (95th percentile), and agrees with the decrease in ice thickness (Fig. 4).

The changes in local ice flow between the first and the last 10 years of the 2W experiment are also related to changes in surface slope and ice thickness, particularly at the margins. To investigate the ice flow changes at the local scale, we used the examples of the Jakobshavn (western coast) and the Kangerlussuaq (eastern coast) glaciers, for which the fine-scale structures of the ice velocity, obtained after the GRISLI initialization procedure, are relatively well reproduced compared to the observations (Fig. 7a–b).

Figure 7Regional zoom over the Jakobshavn (left panels) and the Kangerlussuaq (right panels) glaciers for the 2W experiment. Panels (a) and (b) are the observed velocities (in m yr⁻¹) from Joughin et al. (2018). Panels (c) and (d) are the simulated velocity (in m yr⁻¹) after the initialization procedure. Panels (e) to (j) represent the velocity anomalies (in m yr⁻¹) between the 2140–2150 and the 2000–2010 mean periods of the 2W experiment for the vertically averaged surface velocity (e and f), the SIA velocity component (g and h) and the SSA velocity component (i and j). Note that a logarithmic scale is used for the velocity anomalies. The dashed lines correspond to the 500 m surface elevation iso-contours for the present-day observed topography. The grey shading represents the non-ice-covered areas and the blue shading is the ocean mask.

Download

For the Jakobshavn glacier, the ice sheet areas located above 1500 m are mainly characterized by an increase of more than 15 m yr⁻¹ (i.e. 10 %) of the vertically averaged velocity as a result of increasing surface slopes (Fig. 7c). Conversely, areas below 1000 m are dominated by a slow-down of the ice flow of more than 200 m yr⁻¹ (i.e. 29 %) due to the decreasing ice thickness (Fig. 7e). For altitudes above 500 m, the vertically averaged velocity is mainly driven by the SIA velocity (Fig. 7e–g). On the contrary, below 500 m, basal sliding velocities are large due to low basal drag coefficient (see Fig. 3 in Le clec'h et al., 2018) and the SSA velocity component dominates the ice flow (Fig. 7e–i). However, while basal drag is lower in locations below 500 m, the ice flow is limited by the strongly reduced ice thickness (Fig. 4).

The Kangerlussuaq glacier is located in regions where the bedrock is characterized by a succession of valleys surrounded by mountains merging in a canyon where the deepest part is located 100 km away from the coast (Morlighem et al., 2017). The ice flow of the Kangerlussuaq is therefore divided in different branches with increasing ice velocities towards the ice sheet margin and becoming even larger when merging in the canyon (Fig. 7d). As for the Jakobshavn glacier, the ice flow accelerates at the end of the 2W experiment as a consequence of the increase in surface slope for high altitudes (∼2000–2500 m; see Fig. 4). Conversely, a strong decrease of the ice flow is found in most of margin regions (Fig. 7f) directly related to the ice thinning (Fig. 4). Contrary to the case of the Jakobshavn glacier that presents large basal sliding velocities only below 500 m, the Kangerlussuaq shows low basal drag coefficients in the entire glacier (see Fig. 3 in Le clec'h et al., 2018) and thus the ice flow is mainly governed by the SSA component (Fig. 7j).

These results are in line with Peano et al. (2017), who also found a decreasing ice flux at the end of the 21st century (w.r.t 1970) in downstream regions of the Jakobshavn and Kangerlussuaq glaciers as a consequence of ice sheet thinning at the margins.

4.2.1 Impact on SMB and ST

To assess the importance of the atmosphere–GrIS feedbacks, we now compare the NF and the 2W experiments. The main SMB differences between both experiments, averaged over the 2140–2150 period, highlight higher SMB values in NF compared to 2W for altitudes below 2000 m, with the exception of some margin locations in the eastern part (Fig. 8a). This SMB difference is driven by a snowfall increase in low-altitude areas (Fig. S6) and by the runoff decrease in NF with respect to 2W (Fig. S7). This decreased runoff results from colder temperatures over the whole GrIS (up to −1 ^∘C in the western and northern parts, Fig. 8b), except in the regions at the very edge of the GrIS, which sees a significant warming (up to 8 ^∘C, Fig. 8b) despite an increase in ice thickness (w.r.t 2W). The cooling can be explained by the absence of the temperature–altitude feedback in the NF experiment. Indeed, taking this feedback into account in 2W, results in lower altitudes as the ice thickness decreases (Sect. 4.1.2 and Fig. 4) and therefore in warmer 2W temperatures compared to NF. The warming simulated in NF (w.r.t 2W) over the very edge of the ice sheet can be explained by changes in atmospheric circulation.

Figure 8Mean differences (NF–2W) for the 2140–2150 mean period (a) annual surface mass balance (in m yr⁻¹), (b) annual surface temperature (in ^∘C), (c) ice thickness (in m) in 2150 between NF and 2W. Note that these differences are computed on the GRISLI grid. The dashed lines correspond to the 500 m surface elevation iso-contours of the present-day observed topography. The grey shading represents the non-ice-covered areas.

Download

Indeed, unlike 2W, NF allows for an explicit computation of changes in ice sheet surface slopes due to increased melting at the margin. This has important consequences on the atmospheric circulation and in particular on the katabatic winds (Fig. 9). Over the ice sheet, the surface slopes simulated in NF in 2150 are less steep compared to those simulated in 2W (discussed in Sect. 4.1.2) lead to a slight increase in katabatic winds (Fig. 9). However, at the ice sheet margin, i.e. where the ice mask in MAR is below 100 %, there is a substantial decrease in surface winds. This stems from the fact that the changes in surface elevation seen by the atmospheric model are computed from the aggregated changes in GRISLI at 5 km. As such, a non-zero fraction of tundra, which presents no change in surface elevation, results in smaller elevation changes compared to grid cells with permanent ice cover only. This induces in 2W lower surface slopes at the margin with respect to the interior and thus weakened surface winds in these regions. Altogether, the slightly stronger katabatic winds over the ice sheet and their weakening at the margin lead to a cold air convergence towards the ice sheet edge, absent in NF (Figs. 8b, 9 and S8–S9). Another consequence of the stronger katabatic winds in 2W (w.r.t NF) due to increased surface slopes in the GrIS interior is to enhance the atmospheric exchanges along the slope of the ice sheet. The area with lower atmospheric pressure generated by the stronger katabatic winds is filled in by the warmer air coming from higher atmospheric levels in the boundary layer. The warming of the upper part of the boundary layer in 2W combined with the lower surface elevation explains the ST increases in the interior of the GrIS.

Figure 9Surface wind speed difference (shaded) between the NF and the 2W experiments for the 2140–2150 mean period. Black arrows represent the wind direction in the 2140–2150 mean period of the 2W experiment. The length of the arrows indicates the magnitude of wind speed. The grey shaded area stands for the extent of the region for which the permanent ice fraction is 100 % (no tundra). Red solid line indicate the ice sheet extent. The dashed black lines correspond to the 500 m surface elevation iso-contours.

Download

4.2.2 Impact on ice thickness and ice dynamics

The most important ice thickness difference between the last decade of the NF and the 2W experiments is a higher ice thickness in NF compared to 2W. As mentioned in the previous section, this is mainly explained by the positive SMB–elevation feedback in 2W that results in increased surface temperatures compared to NF, and thus increased runoff, when surface elevation decreases. Areas with this type of behaviour cover most of the Greenland ice sheet slopes and reach the interior of the ice sheet from the western or the northeastern margins. The largest changes occur over the western edge of the GrIS, where the thickening between NF and 2W reaches more than 50 m (Fig. 8c). The ice thickness difference pattern is essentially mimicking the SMB differences between NF and 2W (Fig. 8a), suggesting that the two-way coupling induces only a relatively limited change in ice flow, as shown by the ice flux divergence anomaly (Fig. S10), although the surface velocities (Fig. 6b) are slightly higher in NF due to higher ice thickness (Fig. 8c).