Ice flux evolution in fast flowing areas of the Greenland ice sheet over the 20th and 21st centuries

Abstract

This study investigates the evolution of Greenland ice sheet flux focusing on five of the main fast flowing regions (Petermann glacier, North East Greenland Ice Stream, Kangerdlugssuaq glacier, Helheim glacier and Jakobshavn glacier) in response to 20th and 21st century climate change. A hybrid (shallow ice and shallow shelf) ice-sheet model (ISM) is forced with the combined outputs of a set of seven CMIP5 models and the regional climate model MAR. The ISM simulates the present-day ice velocity pattern, topography and surface mass balance (SMB) in good agreement with observations. Except for the Kangerdlugssuaq glacier, over the 21st century all the fast-flowing areas have exhibited a decrease in ice flux as a result of a negative SMB rather than dynamical changes. Only the fronts of Kangerdlugssuaq and Helheim glaciers have shown an interannual variability driven by dynamical rather than climate changes. Finally, the results predict a substantial inland ice margin retreat by the end of the 21st century, especially along the northern coasts.

1. INTRODUCTION

Over the last few decades, the Greenland ice sheet (GIS) has increasingly lost mass from 0.09 mm a⁻¹ (1992–2001) to 0.59 mm a⁻¹ (2002 to 2011) (e.g. Shepherd and others, 2012; Vaughan and others, 2013; Velicogna and others, 2014). Numerical models suggest that the cumulative mass loss from Greenland at the end of the 21st century could reach 27 cm of sea level equivalent (SLE) for the worse climate scenario (e.g. Church and others, 2013; Fettweis and others, 2013; Yan and others, 2014; Fürst and others, 2015). About half of the mass lost from Greenland today is due to surface melt, through runoff, while the other half results from solid ice discharge due to calving at the marine margins of the ice sheet (e.g. Rignot and Kanagaratnam, 2006; van den Broeke and others, 2009). The relative importance of these two mechanisms varies among the different regions of Greenland. For example ice discharge is predominant in the Baffin Bay, Davis strait and Atlantic sectors, while the surface melt dominates in the Arctic, Labrador Sea and Greenland Sea sectors (e.g. van den Broeke and others, 2009; Sasgen and others, 2012; Vaughan and others, 2013).

The dynamical ice discharge is controlled by the acceleration of fast-flowing marine terminating outlet glaciers (e.g. Howat and others, 2007; Nick and others, 2013) triggered by different mechanisms, such as: decreased ice flow resistance due to glacier frontal retreat (e.g. Howat and others, 2005, 2007), reduced basal and lateral drag (e.g. van der Veen and others, 2011), increased melting at the ice/ocean interface (e.g. Holland and others, 2008; Murray and others, 2010; Straneo and others, 2010), changes in basal lubrication (e.g. Zwally and others, 2002; Parizek and Alley, 2004; Sundal and others, 2011) and response to fast re-equilibration of calving-front geometry after calving events (e.g. Howat and others, 2007). Because of the interplay among all these mechanisms (e.g. Vieli and Nick, 2011), these regions have been extensively monitored in order to improve present sea-level contribution estimates (e.g. Sohn and others, 1998; Howat and others, 2007; Holland and others, 2008; Joughin and others, 2008a, b, c; van der Veen and others, 2011; Moon and others, 2012; Khan and others, 2014; Münchow and others, 2014; Mouginot and others, 2015).

The Jakobshavn (JKB) glacier (West Greenland, Fig. 1d) is one of the major outlet glaciers since it drains ~7% of GIS ice mass (Csatho and others, 2008). It is also one of the outlet glaciers that has shown the largest retreat in recent decades. In fact, between 1995 and 2005 the ice velocities near the terminus of the glacier doubled and the glacier has been continuously thinning since the late 1990s (e.g. Sohn and others, 1998; Holland and others, 2008; Joughin and others, 2008c; van der Veen and others, 2011). Many mechanisms have been proposed to explain the behaviour of the JKB glacier, however ocean warming is considered as the most probable trigger for the recent frontal speed-ups.

Fig. 1. (a) Observed Greenland topography (Bamber and others, 2013). (b) Simulated Greenland topography at the end of the spin-up run. (c) Differences between simulated topography at the end of the 24 ka spin-up and the observed topography. (d) Observed Greenland velocities (Joughin and others, 2010). The red boxes show the locations of the five studied areas. (e) Simulated Greenland velocities at the end of the spin-up run. (f) Differences between simulated velocities at the end of the 24 ka spin-up and the observed velocities.

Similarly, on the eastern coast, the frontal ice velocities of Helheim (HLH) and Kangerdlugssuaq (KGL) glaciers increased from 2001 to 2006 (Joughin and others, 2008b). Ice velocities accelerated substantially in 2002 in relation to the HLH glacier and in 2005 for the KGL glacier, but slowed down after 2006 (e.g. Howat and others, 2007; Murray and others, 2010). Changes in sea surface temperature (Murray and others, 2010), and the seismicity created by many large iceberg-calving episodes might be behind the variations in ice velocities (Joughin and others, 2008b). The different response timings of these glaciers has been related to the differences in the bedrock topography of the glaciers (Joughin and others, 2008b).

In contrast, the Petermann glacier (PTM) (North Greenland, Fig. 1d) was more stable than the others (e.g. Joughin and others, 2010; Nick and others, 2013) up to the end of the first decade of the 21st century, when two large calving events, in 2010 and 2012, led to a 33 km reduction in the length of its floating ice shelf and to ~15–30% increase in velocity compared with velocities observed before 2010 (Münchow and others, 2014). Münchow and others (2014) relate these calving events to the formation of crevasses and to a thinning of the ice shelf, which weaken the structural integrity of the ice shelf.

In the north-eastern area of the GIS, three outlet glaciers, namely Nioghalvfjerdsfjorden Glacier (NG), Zachariae Isstrøm (ZI) and Storstrømmen Glacier (SG), drain the North East Greenland Ice Stream (NEGIS). The NEGIS is a unique glaciological feature in Greenland since this ice stream expands far inland, near the central dome of the ice sheet (Fahnestock and others, 2001). After a period of relative stability (1978–2003), the NEGIS has been affected by sustained dynamical changes probably caused by regional atmospheric and ocean warming (Khan and others, 2014). While the ZI has exhibited an accelerated retreat, since 2012 the NG has retreated much more slowly in comparison, which Mouginot and others (2015) propose could be linked to the differences in their bed geometry. On the other hand, the SG is a surge glacier (a surge occurred between 1978 and 1984, Reeh and others, 2003), which explains the observed thickening of this glacier in the period 2003–12 (Khan and others, 2014).

To investigate the mechanisms and processes behind present-day and future changes in the GIS and to reproduce the complex pattern of GIS ice velocities, numerical models of different spatial and physical complexities, spanning from single outlet glacier models (e.g. Nick and others, 2010, 2013) to large-scale ice-sheet models (ISMs) (e.g. Price and others, 2011; Seddik and others, 2012; Goelzer and others, 2013; Greve and Herzfeld, 2013), have been used to investigate the dynamics of the main fast-flowing areas. The studies based on these models suggest that the sea-level contribution of the GIS for the next century is dominated by surface mass balance (SMB) changes, while the outlet glacier dynamics only accounts for a small fraction of this contribution. These two processes are mutually competitive in removing mass from the GIS. For example, Goelzer and others (2013) show that the relative importance of outlet glacier dynamics in removing mass from the GIS decreases with decreasing SMB conditions. Using a three-dimensional (3-D) higher order model combined with a specific outlet glacier flow-line model, Goelzer and others (2013) estimate that outlet glacier dynamics will only account for 6–18% of the sea-level contribution by 2200 under a A1B SRES scenario (Nakićenović and others, 2000). Price and others (2011) use a 3-D higher-order ice flow model to which they apply a perturbation at the stress boundary condition at the marine termini of JKB, HLH and KGL outlet glaciers to evaluate their dynamical response. They estimate a total 21st century cumulative sea-level contribution from the GIS of ~46 mm SLE, with 13% coming from dynamical perturbations. Using an outlet glacier model, Nick and others (2013) calculate that ~80% of the sea-level contribution up to 2100 under the A1B SRES scenario (Nakićenović and others, 2000) from PTM, KGL, JKB and HLH is likely to be due to dynamical changes, such as margin retreat or increased calving.

The mass conservation equation implies that the ice fluxes (in particular their divergence) together with surface and basal melting of the ice sheet determine the rate of change in the ice thickness (e.g. Greve and Blatter, 2009; Seroussi and others, 2011). While SMB and outlet glacier basal melting can be computed independently from the ISMs, this is not the case for the dynamical changes resulting from ice thickness and calving variations. Therefore, an ISM is crucial to assess the variations in ice fluxes, especially for the fast-flowing areas, and, ultimately, ice volume variations (Seroussi and others, 2011).

To model the dynamical changes in the fast flowing areas, high-resolution and higher-order numerical models could provide more accurate estimates than shallow ice approximation (SIA) and/or shallow shelf approximation (SSA) models. This is because the latter are not able to properly capture all the dynamical changes occurring at the ice-sheet margins due to the assumptions on which SIA relies (e.g. Seroussi and others, 2011). However, Greve and Herzfeld (2013) show that a SIA-only model, namely SICOPOLIS, with a high horizontal resolution of 5 km, satisfactorily reproduces the observed velocity fields in areas characterized by a complex topography and by fast-flowing conditions. Also, the hybrid SIA/SSA model GRISLI has been shown to perform very well in a GIS projection of sea-level contribution compared with full-Stokes models (Edwards and others, 2014). Compared with more comprehensive ISMs, the SIA and SIA/SSA-models suffer from various limitations and assumptions (see Section 4), however they are the only ISMs that have been coupled to climate models to date (e.g. Ganopolski and others, 2010; Lipscomb and others, 2013). Further advantages of SIA and SIA/SSA models include: (1) the possibility to perform long-term simulations, which are required, for example, in determining the tipping point of the GIS. On the other hand, for more computationally demanding higher order models, this solution is still not feasible; (2) higher-order and full-Stokes models rely on SIA or hybrid models to spin-up the vertical ice temperature and ice velocities (e.g. Seddik and others, 2012) in order to benefit from a long-term initialization of the ice-sheet thermodynamics; (3) similarly to what is done with climate models, in order to assess the GIS response to ongoing and future climate change, it is important to keep a hierarchy of the physical complexities of ISMs (e.g. Blatter and others, 2011; Kirchner and others, 2011). This facilitates an objective perspective of the system.

The present work evaluates the link between the changes in SMB induced by the 20th and 21st century climate changes, the GIS geometrical evolution and the ice flux changes. We focus on the ice stream evolution of five of the main fast flowing areas in Greenland, namely PTM, NEGIS, KGL, HLH and JKB. We perform our simulations using the SIA/SSA ISM GRISLI (Ritz and others, 2001). Consequently, the present work also assesses the ability of a hybrid ISM to properly simulate observed and predicted changes in the GIS. We force GRISLI with the surface air temperatures and precipitations of the 20th and 21st centuries simulated by a set of seven atmosphere-ocean coupled general circulation models (AOGCMs) within the framework of the 5th phase of the Coupled Model Intercomparison Project (CMIP5, Taylor and others, 2012) and with 20th century surface air temperature and precipitation from the regional atmospheric model MAR. The 21st century projections are based on the IPCC's Representative Concentration Pathways (RCP) 4.5 and 8.5 emissions scenarios (Moss and others, 2010).

We now outline the main physical characteristics of the ISM GRISLI, and its initialization method as well as the climate forcing from the CMIP5 AOGCMs. We then analyze the results obtained with our simulations and discuss their implications with regard to the use of SIA/SSA ISMs and on the GIS ice flux evolution.

2. METHODS

2.1. The GRISLI ISM

The GRenoble Ice Shelf and Land Ice model (GRISLI, Ritz and others, 2001) is a 3-D thermo-mechanical ISM in which the dynamics of grounded ice is described by the SIA (Hutter, 1983) and fast ice streams and ice shelves are described by the SSA (MacAyeal, 1989). The grounding line is determined by a floatation criterion. The basal drag coefficient is inferred with an iterative inverse method based on the observed surface velocities (Edwards and others, 2014, and Supplementary Figure S1). The inverse method relies on the minimization of a cost function that measures the differences between observed and simulated velocities as a function of the basal drag β (e.g. Arthern and Gudmundsson, 2010; Morlighem and others, 2010; Gillet-Chaulet and others, 2012).

GRISLI is run on a rectangular regular Cartesian grid at a 5 km × 5 km horizontal resolution using the polar stereographic projection with the standard parallel at 71°N and the central meridian at 39°W. In the vertical, GRISLI accounts for 21 layers of cold ice. The initial present-day bedrock elevation and ice-sheet thickness maps are derived from Bamber and others (2013).

GRISLI is forced by monthly surface air temperature and total precipitation provided by the combination of MAR reference climatology and CMIP5 anomaly fields (further details in Section 2.3). To account for surface elevation changes during runtime, air surface temperature is corrected by an atmospheric uniform lapse rate (λ = 6.309 °C km⁻¹, Fausto and others, 2009). Changes in precipitation are related to changes in temperature by a precipitation correction factor (γ = 0.07, Quiquet and others, 2012, and further details in Table 1), assuming that the saturation pressure of water vapour depends exponentially on the temperature (Charbit and others, 2002).

Table 1. Main characteristics and model parameters of GRISLI

The accumulation corresponds to the solid fraction of total precipitation, which is calculated using the empirical relationship introduced by Marsiat (1994). Ablation is computed using the positive degree-day (PDD) semi-empirical method by Reeh (1991), which linearly relates the number of positive degree days to the amount of ice and snow melt through the degree-day factors (DDF) C _snow and C _ice . Both factors are derived from observations and vary spatially as a function of July temperature following Tarasov and Peltier (2002). The main features of GRISLI are summarized in Table 1, and further details on the model physics can be found in Ritz and others (2001).

2.2. Initialization of the ISM

Because of the slow diffusivity of surface temperature through the ice column and its influence on ice velocities, a spin-up is required to obtain realistic initial present-day vertical temperature conditions (Rogozhina and others, 2011; Yan and others, 2013). We thus performed a long-term transient spin-up based on the index method (e.g. Cuffey and Marshall, 2000; Lhomme and others, 2005; Greve and others, 2011; Applegate and others, 2012; Quiquet and others, 2013), according to the initialization protocol of the SeaRISE project (e.g. Greve and others, 2011; Bindschadler and others, 2013; Nowicki and others, 2013). To balance initialization-time and computational-time, the spin-up simulation covers the last 24 ka of the last glacial cycle, starting slightly before the Last Glacial Maximum (21 ka BP), similarly to Herzfeld and others (2012). We performed a 125 ka spin-up for comparison. The temperature profiles obtained both inland (at GRIP site) and along the margins (at DYE3 site) are very similar to those obtained after 24 ka. We believe that the discrepancy between the two spin-ups does not significantly affect the conclusions of this work (Supplementary Figure S2). Finally, the ice-sheet topography is tuned by adding a correction term (Supplementary Table S1) to the melting factors C _snow and C _ice following Greve and others (2011).

At the end of the spin-up, compared with the observed ice topography (Bamber and others, 2013), the simulated GIS topography exhibits a RMSE of ~183 m on average, while the simulated surface ice velocities show a RMSE of 192 m a⁻¹ on average over the GIS compared with the observations (Joughin and others, 2010). We attribute this good agreement with the observations above all to the inverse method used to retrieve the basal drag (β) at the base of the ice sheet (see Section 2.1 and Edwards and others, 2014). This means that the discrepancies result from the physics of the model and the use of a free-surface spin-up simulation (other studies, for example Gillet-Chaulet and others, 2012; Seddik and others, 2012; Edwards and others, 2014, maintain a fixed topography for the spin-up).

More specifically, the simulated GIS topography is ~150 m thinner than in Bamber and others (2013) in the northern and central parts of Greenland (Fig. 1c). Along the eastern coasts, the simulated ice is ~200 m thicker than in Bamber and others (2013), this implies steeper slopes and 80% faster surface velocities than observed (Fig. 1c, Joughin and others, 2010). However, along the western coast, the simulated ice is ~300 m thinner (up to 700 m in the north-western GIS region) compared with Bamber and others (2013), this implies smoother slopes and 60% slower surface velocities than in the observations (Fig. 1f).

These discrepancies between our simulated spun-up ice topography and velocities and observations are similar to simulated findings in previous studies (e.g. Greve and others, 2011; Seddik and others, 2012; Nowicki and others, 2013). Nowicki and others (2013) show that, in general, long-term transient spin-ups, as in the present study, simulate a thicker peripheral ice sheet and a thinner interior, as obtained at the end of our spin-up (Fig. 1c). In the following we use the spun-up ice-sheet topography, vertical ice temperature and velocities as initial conditions for all the subsequent 20th and 21st centuries GIS simulations. In doing so, we show that long-term glacial climate history thus impacts on the future dynamics of ice.

2.3. Climate forcing

To investigate the future evolution of GIS ice fluxes, we used the climate fields computed within the framework of the CMIP5 project (Taylor and others, 2012). More specifically, and in order to account for the full range of model uncertainty, we used a set of climate forcing computed by seven AOGCMs whose climates were simulated from pre-industrial times (1850) to the end of the 21st century (Table 2). This ensures a continuity of the climate forcing used to force our simulations. The projected climates (2006–2100) follow two future emission scenarios, the RCP 4.5 and 8.5 (Supplementary Figure S3, Moss and others, 2010).

Table 2. Simulated climate forcing from a set of seven CMIP5 AOGCMs (http://cmip-pcmdi.llnl.gov/cmip5/) used in this study. The last two columns report the number of realizations for each CMIP5 AOGCM accounted for in the RCP 4.5 and RCP 8.5 future scenarios

In general, Yan and others (2014) show that the spatial distribution of the simulated CMIP5 air surface temperature and precipitation are in good agreement with the ERA Interim reanalysis (Dee and others, 2011) for 1979 to 2005. Individually, CMIP5 AOGCMs do not properly capture the interannual variability and linear trends over the same period (Yan and others, 2014). However, Yan and others (2014) show that the multi-model ensemble mean (MME) provides better results than individual models in terms of spatial distribution, interannual variability and trend.

Due to the coarse horizontal resolution of the AOGCM horizontal grid compared with that of GRISLI and due to persisting climate biases between observed and simulated present-day GIS climates (e.g. Bromwich and others, 2007; Walsh and others, 2008), the use of these climate forcing could lead to misrepresentations of the SMB (e.g. Quiquet and others, 2012). Thus GRISLI is not directly forced with the CMIP5 AOGCMs simulated climates. In order to minimize the biases, the climate forcing are reconstructed following an anomaly method (e.g. Huybrechts and others, 2004; Yan and others, 2014), which consists of adding climate anomalies from a coarser climate model on top of a high resolution climate dataset. In the present work, the simulations of the regional climate model (RCM) MAR (version 3.5, 25 km × 25 km, Fettweis and others, 2013) averaged between 1980 and 1999 are considered as the high-resolution reference climatology on top of which the coarse CMIP5 AOGCM monthly climate anomalies varying from 1850 to 2100 are added. For each month of the climate time series, the climate forcing $(T_{forcing}^{1850 - 2100}, P_{forcing}^{1850 - 2100} )$ are calculated as follows:

$$T_{\,forcing}^{1850 - 2100} (t) = T_{ref,mean}^{1980 - 1999} + (T_{aogcm}^{1850 - 2100} (t) - T_{aogcm,mean}^{1980 - 1999} )$$

(1) $$P_{\,forcing}^{1850 - 2100} (t) = P_{ref,mean}^{1980 - 1999} \times (P_{aogcm}^{1850 - 2100} (t)/P_{aogcm,mean}^{1980 - 1999} )$$

where $T_{ref,mean}^{1980 - 1999} (P_{ref,mean}^{1980 - 1999} )$ is the MAR surface monthly air temperature (precipitation) averaged over the period 1980–1999, $T_{aogcm}^{1850 - 2100} \;(P_{aogcm}^{1850 - 2100} )$ is the monthly surface air temperature (precipitation) from the CMIP5 AOGCMs, and $T_{aogcm,mean}^{1980 - 1999} \;(P_{aogcm,mean}^{1980 - 1999} )$ is the averaged surface air temperature (precipitation) over the period 1980–1999 from the CMIP5 AOGCMs.

Note that all climate forcing fields (for an overview of the models, see Table 2) are downscaled on the ISM grid by means of a simple downscaling method, using an atmospheric lapse rate (6.309 °C km⁻¹) and a precipitation correction factor (7 %/°C) to correct temperature and precipitation for the difference in elevation (see Section 2.1).

3. RESULTS

3.1. Surface mass balance

Although GRISLI uses simple melting and accumulation parametrizations, the model is able to capture the present-day spatial features of the observed GIS SMB. SMB simulated values result in a good agreement with previous modelling and observational studies (212 Gt a⁻¹ on average, Supplementary Figure S4 and Tables S2 and S3). GRISLI simulates accumulation values of ~629 Gt a⁻¹ on average, which is within the range of values from the previous literature and close to the value simulated with MAR (Supplementary Table S3, Fettweis and others, 2013). The snow/precipitation fractioning method from Marsiat (1994) yields reasonable snow accumulation rates. Conversely, the simulated ablation values (417 Gt a⁻¹ on average) overestimate the literature range by ~40% (Supplementary Figure S3).

The PDD method generally overestimates the melting especially under warm climate conditions (e.g. van de Berg and others, 2011) and uses a spatially uniform amount of refreezing. Refreezing is the main cause of discrepancies between the PDD and the surface energy balance models (SEBMs) (e.g. van de Wal, 1996; Bougamont and others, 2007). Although GRISLI overestimates the melting and slightly underestimates the accumulation compared with MAR, the MME SMB evolution is in good agreement with the SMB simulated with MAR (Supplementary Figure S5). However, our simulations fail to capture the SMB interannual variability, probably because the CMIP5 AOGCMs fail to simulate the observed interannual climate variability (e.g. Yan and others, 2014).

At the regional scale, the simulated SMB agrees relatively well with RCM results (Supplementary Figure S6, Tedesco and Fettweis, 2012; Vernon and others, 2013), except in the northern and western regions, where the simulated melting is overestimated (Supplementary Figure S6). We simulated negative SMB values in the northern and north-eastern basins, while Vernon and others (2013) and Tedesco and Fettweis (2012) simulated SMB values close to zero, but positive. In fact, Vernon and others (2013) show a large decrease towards negative SMB values in these regions already in the period 1996–2008. As a result, the overestimation of the melting in the northern and north-eastern basin may slightly accelerate the ice margin retreat and, thus, affect the simulated ice flux evolution of the PTM glacier and NEGIS, which are located in these regions.

By the end of the 21st century, GRISLI predicts a total SMB decrease in the range 0–28.8 cm SLE under the RCP 4.5 scenario, and in the range 2.7–51.6 cm SLE under the RCP 8.5 scenario. When forcing GRISLI with the MME climate, our simulations yield a mean cumulative SMB decrease of ~10 cm SLE (~20 cm SLE) by the end of the 21st century under the RCP 4.5 (RCP 8.5) emission scenario, in line with recent studies also using the CMIP5 simulations as climate forcing (e.g. Yan and others, 2014; Fürst and others, 2015).

At the regional scale (Fig. 2), our ice simulation ensemble mean exhibits the largest reduction in the northern and western basins. During the 21st century, Fettweis and others (2013), for example, suggest that the northern half of Greenland will be influenced by the projected reduction in the Arctic sea-ice extent and the subsequent changes in local albedo. For the western basin, Belleflamme and others (2013) show that the GIS is more sensitive to the changes in atmospheric circulation. Conversely, the southern and south-eastern basins exhibit small variations during the 21st century under the RCP 4.5 scenario due to simulated balanced conditions between ablation and accumulation in these regions. Under the RCP 8.5 scenario, a general acceleration in the SMB decrease is predicted after 2030–40 over the entire GIS.

Fig. 2. Regional distribution of the time evolution of the SMB. The values are obtained from the mean of the multi-model ensemble of the complete set of simulations. The extremes of the ensemble spread are shown as dashed lines. The ‘light coloured’ lines correspond to the RCP 4.5 scenario, the ‘dark coloured’ lines correspond to the RCP 8.5 scenario. Note that each panel has a different vertical scale.

3.2. Ice stream evolution over the 20th and 21st centuries

3.2.1. Simulated ice velocities between 1980 and 1999

As shown in Section 2.2, the initialization method that we used resulted in modern surface ice velocities, which are broadly in good agreement with Joughin and others (2010). We now analyze the discrepancies in ice velocities between our 20th century simulations and Joughin and others (2010), focusing on the five largest fast-flowing GIS areas: the Petermann glacier (PTM), the NEGIS, the Kangerdlugssuaq glacier (KGL), the Helheim glacier (HLH) and the Jakobshavn glacier (JKB) (Fig. 1d).

The ice stream velocity patterns are generally well captured by GRISLI (Fig. 3). The HLH and KGL glacier regions are the areas where GRISLI shows the best performance (Supplementary Figure S7). GRISLI is able to simulate the multiple-tributary structure of the KGL glacier, which is fed by ice fluxes from three separate valleys that converge in the vicinity of the coast in a single tongue (Figs 3e and f). This good performance results from the spun-up initial topography, which is in very good agreement with Bamber and others (2013), particularly in these areas (Fig. 1).

Fig. 3. (a, c, e, g, i) Observed (Joughin and others, 2010) and (b, d, f, h, j) simulated velocities (m a⁻¹) in the five regions studied. In particular (a, b) for the PTM glacier region, (c, d) for the NEGIS region, (e, f) for the KGL glacier region, (g, h) for the HLH glacier region and (i, j) for the JKB glacier region. Note that the grounding lines are drawn for the PTM and NEGIS regions to highlight the location of the floating points. The red lines exhibit the section for the upstream ice flux calculation in each region, likewise the blue lines for the downstream ice flux.

In contrast, compared with Joughin and others (2010), GRISLI underestimates the ice velocities for the JKB glacier (Figs 3i and j), especially downstream where the simulated ice velocities reach ~1.2 km a⁻¹ against ~11 km a⁻¹ in the observations. This discrepancy arises from the horizontal resolution (5 km), which is too coarse to properly capture the narrow valley of the JKB glacier. Aschwanden and others (2016) show that at least a horizontal resolution of 2 km is required to properly capture the JKB glacier valley topography and reproduce its high surface velocities.

Finally, GRISLI is not able to maintain the floating tongue of the PTM glacier and of the NEGIS (Figs 3a–d). There are two reasons for this. Firstly, the floating tongues retreat at the end of the spin-up as a result of the calving criterion, which correspond to a prescribed uniform cutting ice front thickness. GRISLI calculates whether or not the ice fluxes through the grounding line can sustain the floating tongue, if they cannot, or if the front of the floating tongue reaches the cutting thickness value, the front of the floating tongue is calved. Secondly, the inaccuracy in the SMB during the spin-up phase and in the 20th century simulations, stops the floating tongues from growing again (Greve and Herzfeld, 2013). Consequently, GRISLI fails to simulate the velocities and the ice topography in the coastal part of these regions. In contrast, the ice velocities in upstream regions of the PTM glacier and the NEGIS are well simulated by GRISLI (Figs 3a–d, and Supplementary Figure S7).

3.2.2. Projected ice flux evolution up to 2100

To estimate the future evolution of the ice flux in the five main GIS ice streams, calculations were performed through different vertical sections in each ice stream, one (or two) located upstream (red lines in Fig. 3) and one located downstream near the coasts (blue lines in Fig. 3). For each section, ice fluxes were computed by multiplying the vertical section and the vertically integrated ice velocity perpendicular to the section. In the following analysis, we relate the evolution of the ice fluxes to the variability in ice dynamics and to climate changes, through the equation of continuity $dH/dt = - \nabla Q + SMB + b_{melt} $ (H = ice thickness; Q = ice fluxes; b _melt  = basal melting).

As a result of the simulated negative SMB in the northern half of Greenland during the 21st century, the PTM glacier and the NEGIS ice fluxes decrease upstream (Figs 4a and b), as well as downstream (Figs 5a and b). In fact, the negative SMB causes a retreat of the grounding line of several kilometres inland during the 21st century, which influences the ice flux evolution in these glaciers (Table 3). In the NEGIS region, the ice sheet retreats forever behind the location of the section chosen for the coastal flux calculation by the end of the 21st century (Fig. 5b). Both glaciers exhibit a decreasing trend in ice thickness and ice velocities during the 21st century (Figs 6 and 7). The ice velocity through the upstream section of the NEGIS slightly increases during the 21st century (Fig. 7), due to the ice-sheet retreats (Table 3) and the subsequent steepening of the ice-sheet slopes. This increase in ice velocity, however, is not large enough to counteract the strong reduction in ice thickness resulting from the negative SMB.

Fig. 4. Evolution of the fluxes through the upstream sections (red lines in Fig. 3) of the five regions studied: (a) PTM glacier region, (b) NEGIS region, (c) KGL glacier region, (d) HLH glacier region, and (e) JKB glacier region. The values, which are given in ${\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ , cover the period from 1970 to 2100 and are shown as anomalies with respect to the 1980–99 fluxes: ${\rm 7}{\rm. 6}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (PTM), ${\rm 15}{\rm. 1}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (NEGIS), ${\rm 11}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (KGL), ${\rm 10}{\rm. 8}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (HLH), ${\rm 20}{\rm. 4}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (JKB). The ensemble spreads are shown as dashed areas.

Fig. 5. As Figure 4, apart from the fluxes through the downstream sections (blue lines in Fig. 3). The reference 1980–99 ice fluxes are: ${\rm 4}{\rm. 1}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (PTM), ${\rm 3}{\rm. 3}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (NEGIS), ${\rm 13}{\rm. 3}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (KGL), ${\rm 20}{\rm. 9}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (HLH), ${\rm 13}{\rm. 3}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (JKB). Note that the black vertical dotted lines in panels (a) and (b) show the years when part of the ice sheet retreats behind the coastal gates.

Fig. 6. Evolution of ice fluxes ( ${\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ ), ice velocities (m a⁻¹) and ice thickness (m) in the upstream portion of the PTM glacier. The values cover the period 1970–2100 and are shown as anomalies with respect to the 1980–99 values: ${\rm 7}{\rm. 6}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (ice fluxes), 130 m a⁻¹ (ice velocities), and 744 m (ice thickness). The ensemble spreads are shown as dashed areas.

Fig. 7. As Figure 6, apart from the upstream portion of the NEGIS. The reference 1980–99 values for the MME are: ${\rm 7}{\rm. 6}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (ice fluxes), 106 m a⁻¹ (ice velocities), and 870 m (ice thickness).

Table 3. The 21st century grounding line retreat of the GIS (in km) of the front of the five studied ice stream regions

The retreat values correspond to the maximum ice retreat of the ice stream front for each glacier. Note that, given the GRISLI horizontal resolution, a minimum of 5 km retreat or advance is required to be simulated.

Along the eastern coast, the mean simulated SMB remains positive for most of the 21st century, except for the last decades under the RCP 8.5 scenario (Fig. 2). Consequently, both the KGL and the HLH glaciers, which are located in this area, remain stable or advance up to the last decade of the 21st century (Table 3). The KGL glacier exhibits an increase of ~ ${\rm 0}{\rm. 3}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} {\rm (\sim\!\! 0}{\rm. 1}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} )$ in ice flux along its upstream section during the 21st century under the RCP 4.5 (RCP 8.5) scenario (Fig. 4c). Note that this is the only section of the five fast flowing regions analyzed in our study where simulated ice fluxes increase during the 21st century as a result of an increase in ice velocity (Fig. 8). This acceleration is caused by steepening ice-sheet slopes resulting from positive SMB values and thus increasing ice thickness in the upstream part of this glacier. Conversely, negative SMB and decreasing ice thickness are simulated downstream of the same section (not shown).

Fig. 8. As Figure 6, apart from the upstream portion of the KGL glacier. The reference 1980–99 values for the MME are: ${\rm 5}{\rm. 5}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (ice fluxes), 165 m a⁻¹ (ice velocities), and 946 m (ice thickness).

The HLH glacier exhibits a different ice flux evolution compared with the KGL glacier, although both glaciers are located in the same region, and are thus forced by similar climate changes. This different response derives from the differences in the evolution of local SMB, and to the bedrock and fjord geometry (e.g. Joughin and others, 2008b). The upstream section of the HLH glacier is characterized by positive SMB values leading to a small increase in ice elevation (Fig. 9). Therefore, the ice fluxes remain steady through the upstream section (Fig. 4). Due to the buttressing effect created by the bedrock geometry of the fjord (not shown), the amount of ice that reaches the front of the glacier piles up, which increases the ice thickness of the glacier at the front (Fig. 10). This increase in ice thickness at the front, in turn, leads to a reduction in the slope, which results in a reduction in the ice velocities at both upstream and downstream sections during the second half of the 21st century (Figs 9 and 10). This decrease in ice velocities is the main cause for the negative trend in ice fluxes (Figs 4d and 5d). Only in the last decades of the 21st century, is the increase in ice thickness at the downstream section large enough to balance the impact of the decrease in ice velocities, inducing almost steady ice fluxes conditions (Fig. 5d).

Fig. 9. As Figure 6, apart from the upstream portion of the HLH glacier. The reference 1980–99 values for the MME are: ${\rm 10}{\rm. 8}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (ice fluxes), 194 m a⁻¹ (ice velocities), and 1408 m (ice thickness).

Fig. 10. As Figure 6, apart from the downstream portion of the HLH glacier. The reference 1980–99 values for the MME are: ${\rm 20}{\rm. 9}\,{\rm km}^{\rm 3} {\kern 1pt} {\rm a}^{ - {\rm 1}} $ (ice fluxes), 1371 m a⁻¹ (ice velocities), and 694 m (ice thickness).

Observations show that the KGL and the HLH glaciers undergo rapid changes in the coastal area due to processes occurring at the calving front (e.g. Joughin and others, 2008b; Nettles and others, 2008; Andresen and others, 2011). Currently, the HLH glacier front is characterized by many calving events (e.g. Joughin and others, 2008b; Nettles and others, 2008), making this glacier one of the most prolific iceberg exporters in Greenland (e.g. Andresen and others, 2011). Furthermore, Andresen and others (2011) show that the HLH glacier responds to atmosphere-ocean variability on short timescales (3–10 years). In our simulations, small calving events at the front of both glaciers are simulated at a very high frequency (1–2 years) leading to a strong interannual variability in their ice flux evolution at the downstream sections during the 20th and 21st centuries (Figs 5c, d). In fact, GRISLI simulates the fronts of the KGL and the HLH glaciers close to the grounding line, and the absence of a proper description of the grounding line dynamics could explain the high frequency of simulated calving events and the strong instability of these glacier fronts.

Finally, in our simulations, the JKB glacier exhibits a quasi-linear decrease over the 21st century in its upstream region (Fig. 4e), while the coastal section shows a ‘step-like’ evolution (Fig. 5e). In both regions a similar behaviour is simulated under a constant 1980–99 climate (not shown), which indicates that this behaviour is not related to climate changes over the 21st century. Consequently, we relate this behaviour to the horizontal resolution of our grid. In fact, the 5 km horizontal resolution of our simulations inhibits the proper simulation of sub-grid topographic features and, thus, of the JKB glacier dynamics (Fig. 3, and e.g. Aschwanden and others, 2016).

To summarize, the ice flux variations during the 21st century are a consequence of SMB changes, which lead to a reduction in ice thickness and the retreat of the GIS in the northern area of Greenland. However, along the eastern coast of Greenland, local topographic features, such as the buttressing resulting from the bedrock geometry in the HLH glacier region, combined with dynamical features, such as the calving events at the front of the HLH and KGL glaciers, determine how the ice flux evolves at the front. Along the western coast, instead, the JKB glacier exhibits changes that reflect the inaccuracy of the coarse 5 km horizontal grid in simulating the geometry of the narrow glacier valley where the JKB glacier lies.

4. DISCUSSION

4.1. Reference climatology

The experiments performed with GRISLI in this work were forced by means of a set of climate forcing from the CMIP5 AOGCM simulations (Taylor and others, 2012, and Table 2). The surface air temperature and precipitation used to force our ISM were computed using an anomaly method [Eqn (1)]. The climate simulated by the MAR regional climate model (Fettweis and others, 2013) is used as a reference for the 20th century. This choice was based on the ability of this model to simulate the Greenland climate compared with observations (e.g. Fettweis and others, 2011; Franco and others, 2012; Rae and others, 2012). However, other climates could be used as reference climatology, such as reanalysis or other regional climate models outputs, for example RACMO (Ettema and others, 2009). For example, when performing the same experiments with ERA Interim global reanalysis as the reference climate, the mean 20th century SMB is ~10 Gt a⁻¹ larger than when using MAR as the reference climate (not shown). This discrepancy is explained by the fact that MAR surface air temperatures are warmer along the GIS margins than in the ERA-Interim (Supplementary Figure S8). The SMB time evolution over the 20th and 21st centuries simulated with the ERA Interim reference climate is similar to that simulated with MAR reference climatology (not shown). In both cases, the MME simulates a reduction of ~7 Gt a⁻² (~19 Gt a⁻²) under the RCP 4.5 (RCP 8.5) scenario. The ice flux evolution up to 2100 in the five main GIS ice streams show similar trends to those using MAR reference climate, although ice fluxes are lower in some of those areas, such as for the HLH glacier, the JKB glacier and the NEGIS. The lower ice fluxes, obtained when using ERA-Interim climate, might result from the low horizontal resolution grid of these climate forcing. In fact, in these main glacier outlets, there is a systematic underestimation of the precipitation compared with MAR (Supplementary Figure S8a). This could be explained by the fact that the lower reanalysis grid resolution does not capture the complex fjord topography of these glaciers. In addition, MAR simulates a warmer climate along Greenland coasts compared with the ERA-Interim (Supplementary Figure S8b and S8c), leading to faster ice velocities and larger ice fluxes in these regions.

4.2. SMB calculation

Our model computes the ablation by a PDD method. Previous studies have shown that the PDD method is more sensitive to temperature variations in warm climates than more physically-based methods, such as the SEBM. For example, van de Wal (1996) shows that the sensitivity of a PDD model for a 1°C perturbation in air surface temperature is ~20% higher than the sensitivity of an SEBM. Furthermore, Bougamont and others (2007) found that the description of the snow-pack thermo-dynamics has a large impact on the ice-sheet SMB. These discrepancies with the PDD method are mainly due to the fact that the PDD method does not account for albedo changes, and thus does not update surface temperature in relation to the type of glaciated surface (wet snow, dry snow, ice). Therefore, for example, Robinson and others (2010) use the insolation-temperature melt method, based on explicit insolation and albedo changes, which provides more satisfying results of the melt distribution than the basic PDD method. The PDD method also accounts for a fixed and uniform fraction of refreezing of the meltwater, which is the source of the main discrepancy in SMB with the SEBM, together with the lack of radiative feedback, as shown by Bougamont and others (2007). However, in Bougamont and others (2007), the melting calculation from the PDD and from the SEBM are very similar, and most of the difference comes from the refreezing component.

4.3. Ice-sheet initialization

At the end of our spin-up simulation, our experiments exhibit some biases with respect to observed ice-sheet topography and velocities, especially along the coastal areas, where, for example, the floating tongue of the Petermann glacier and of the NEGIS retreat in our simulations. We attributed these biases to the lack of spatial resolution combined with the effect of the long-term spin-up. To reduce these discrepancies, for example, Edwards and others (2014) fix the GIS geometry to observed present-day topography (Bamber and others, 2013), then estimate the temperature structure and retrieve the basal drag coefficient from observed ice velocities (Joughin and others, 2010) using an iterative inverse method. Finally, Edwards and others (2014) relax the ice-sheet geometry to the simulated temperature and ice flow fields. Lee and others (2015) show that the use of an inverse method in ice-sheet initialization helps to obtain better model-data agreement. However, in the Edwards and others (2014) initialization method, the resulting vertical temperature profile and the basal temperature in particular are probably inaccurate since they do not derive from the climate evolution over the last glacial cycle (e.g. Seroussi and others, 2013). Here, we initialize GRISLI with a free surface simulation forced by an index method in which the ice-sheet basal drag is retrieved using an iterative inverse method (Edwards and others, 2014, and Supplementary Figure S1). This methodology slightly reduces the model-data agreement achievable when using fixed surface and inverse methods together (e.g. Edwards and others, 2014; Lee and others, 2015), but it provides a more realistic temperature profile. Consequently, an unique technique to properly initialize ISMs has not yet been found (e.g. Rogozhina and others, 2011).

The index method also takes into account the impact of a Holocene drift on the ice flux evolution. We evaluate this drift by performing a reference simulation forced by a constant present-day climate (i.e. MAR, 1980–1999, Fettweis and others, 2013) over the 21st century. The drift amounts to ~36% of the MME ice flux evolution, on average (not shown). For example, this drift drives the ice flux evolution (~80%) in the JKB glacier region (not shown). In contrast, the evolution of ice fluxes in the upstream area of the NEGIS is slightly (~9%) influenced by the spin-up simulation and, consequently, the climatic forcing is the main driver of ice flux evolution in this region (not shown).

4.4. Ocean/ice interactions

In addition to the impact of the spin-up on the outlet of marine glaciers, spatial resolution and various missing processes, such as the treatment of the grounding line dynamics (e.g. Schoof, 2010), hydro-fracturing of the glacier front (e.g. Pollard and others, 2015), or ice/ocean interaction (e.g. Holland and others, 2008), limit the ability of our model to simulate the dynamics features of the GIS more accurately. With a 5 km grid, only the largest outlet glaciers are captured properly. As such, in our modelling framework, oceanic warming has only a very limited role for GIS evolution. Note that GRISLI simulates the inland retreat of some of the GIS margins, which also means that probably, in the future, the ocean will not directly impact on the GIS dynamics. However, the complete melting of the Arctic sea ice will lead to regional warming as a result of the subsequent decease in albedo (e.g. Day and others, 2013; Fettweis and others, 2013). Thus, a finer horizontal grid resolution is preferable, as shown by Aschwanden and others (2016) to capture all the topography of the fjords. However, large uncertainties in ocean-induced basal melting parametrizations due to the lack of marine observations from the fjords, and in the buttressing process limit the modelling impact of the use of finer horizontal grid resolutions. In addition, the simulation of a climate forcing at a comparable resolution with that of the ISM used is necessary to properly reproduce the dynamics of the fjord glaciers.

4.5. GRISLI basal drag

Our simulations were developed by assuming a constant basal drag over the 20th and 21st centuries. This assumption derives from the lack of any meltwater infiltration scheme in GRISLI that is able to redirect surface meltwater from the surface to the ice-sheet base, where it could increase the basal lubrication and thus modify ice-sheet velocities (e.g. Zwally and others, 2002). Fürst and others (2015) show that the effect of enhanced basal lubrication on the ice volume evolution is negligible on centennial timescales. However, infiltration models have been developed and applied at the regional scale (e.g. Clason and others, 2012; Banwell and others, 2013; Clason and others, 2015), and their implementation in the ISMs could be of interest for future work.

4.6. Hybrid model limitations

GRISLI is a hybrid model (e.g. Kirchner and others, 2011), and is based on SIA and SSA, which restrict its ability to resolve the dynamics in complex topography areas, such as the fjords located along the coasts of the GIS, which are also poorly resolved with a 5 km resolution grid. Higher-order or full-Stokes ISMs could overcome some of these limitations. However, higher-order or full-Stokes models cannot be run in the long-term, for example for paleoclimate spin-up simulations (Gillet-Chaulet and others, 2012; Seddik and others, 2012) due to their high demand on computational resources.

Despite various dynamical limitations, a hybrid ISM can be run over multiple timescales under a large variety of climate forcing and physical parametrizations, thus enabling the sources of uncertainty to be evaluated. In addition, SIA/SSA ISMs satisfactorily simulate present-day and future GIS conditions compared with more complex ISMs (e.g. Edwards and others, 2014). For example, in the case of the main drainage areas, such as the NEGIS, SIA, hybrids, higher-order and full-Stokes all fail to accurately reproduce the observations (e.g. Goelzer and others, 2013; Greve and Herzfeld, 2013).

Finally, future ISM development, especially for fast flowing areas, should include upgraded model physics, increased resolution, and better floating outlet glacier termination representation.

5. CONCLUSION

In the present study, we analyzed the evolution of the ice fluxes of five ice streams and outlet glaciers. The evolution of the ice dynamics was computed using a hybrid ISM, i.e. GRISLI, which simulates a 20th century SMB amplitude and distribution, as well as ice velocity patterns in good agreement with observations and the literature.

In agreement with previous works (e.g. Yan and others, 2014; Fürst and others, 2015), GRISLI simulates a mean present-day SMB of ~212 Gt a⁻¹ (170–469 Gt a⁻¹, literature range), and a mean cumulative decrease of ~10 cm SLE (~20 cm SLE) by 2100 under the RCP 4.5 (RCP 8.5) scenario. However, GRISLI overestimates ablation due to the sensitivity of the PDD method to warm climate conditions.

In general, our results confirm that the SMB evolution is the main driver of the evolution of the GIS over the 21st century. Nevertheless, local scale dynamical features have a significant influence on the interannual variability and on the local ice flux evolution.

Along the eastern coast, the simulated ice flux trends are mostly dominated by climate change up to 2100. However, the calving events and oscillations in the sliding velocities lead to a strong interannual variability of Kangerdlugssuaq and Helheim glacier outlets. Along the western coast, the simulated ice flux evolution in the Jakobshavn glacier region reflects the inaccuracy of the coarse 5 km horizontal grid in simulating the geometry of the narrow glacier valley where the Jakobshavn glacier lies, and, thus, its dynamics.

In the northern area, the GIS is characterized by negative SMB values, which lead to a retreat of the ice sheet in these regions by the end of the 21st century, as shown by the decrease in the ice flux evolution simulated in the Petermann glacier and NEGIS regions.

Finally, in our simulations, the changes in the GIS topography, such as the ice margin retreat and changes in the ice-sheet slopes, triggered by climate changes during the 21st century, exhibit a significant impact on the future GIS evolution, especially in the fast flowing areas.

AUTHOR CONTRIBUTION

Daniele Peano designed and performed the experiments, analyzed the outputs and wrote the manuscript; Florence Colleoni discussed the design of the experiments, discussed the results and participated in the writing; Aurélien Quiquet discussed the design of the experiments, discussed the results and participated in the writing; Simona Masina discussed the results and participated in the writing.