1. Introduction
It is established scientific consensus that the Earth's climate system is warming, and that human influence has been the dominant cause of the observed warming since the mid-20th century (e.g. IPCC, 2021). A major consequence of global warming is sea-level rise, currently (for the period 2006–2018) occurring at a global mean rate of 3.69 ± 0.48 mm a−1. The main sources are melting/discharge of ice sheets, ice caps and glaciers ($\sim {}45\percnt$), thermal expansion of ocean water ($\sim {}38\percnt$), and changes in land water storage ($\sim {}17\percnt$) (Fox-Kemper and others, Reference Fox-Kemper2021). In the long term, the two ice sheets of Antarctica (AIS) and Greenland (GrIS) are the largest potential contributors to global sea-level rise because of their enormous volumes, together amounting to ~ 65 m SLE (sea-level equivalent) (Morlighem and others, Reference Morlighem2017, Reference Morlighem2020). The ice sheets have therefore been the focus of intensive observational as well as modelling efforts.
The Coupled Model Intercomparison Project Phase 6 (CMIP6) is a major international climate modelling initiative (Eyring and others, Reference Eyring2016). As a part of this project, the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6) brought together a consortium of ice-sheet modellers to explore the sea-level-rise contribution from the GrIS and AIS (Nowicki and others, Reference Nowicki2016, Reference Nowicki2020). ISMIP6 focussed on the CMIP6 period from 2015 until the end of 2100. The main findings for the GrIS, when forced by output from CMIP5 global climate models (GCMs), were contributions of 90 ± 50 and 32 ± 17 mm SLE for the unabated warming pathway RCP8.5 (RCP: Representative Concentration Pathway) and the reduced emissions pathway RCP2.6, respectively (Goelzer and others, Reference Goelzer2020). The CMIP6 GCMs tend to feature a warmer atmosphere, which results in higher mass loss due to increased surface melt (Payne and others, Reference Payne2021). For the AIS and CMIP5 climate forcings, ISMIP6 found a mass loss in the range of −7.8 to 30.0 cm SLE under RCP8.5 (Seroussi and others, Reference Seroussi2020). The limited number of results for RCP2.6 fall within this range, and so do the results obtained with CMIP6 climate forcings (Payne and others, Reference Payne2021). This rather unclear picture for the AIS is a consequence of the counteracting effects of mass loss due to ocean warming and mass gain from increased snowfall.
The full suite of ISMIP6 experiments with both CMIP5 and CMIP6 forcings was carried out with the ice-sheet model SICOPOLIS (SImulation COde for POLythermal Ice Sheets, www.sicopolis.net), as documented in detail by Greve and others (Reference Greve2020a, Reference Greve, Chambers and Calov2020b). Chambers and others (Reference Chambers, Greve, Obase, Saito and Abe-Ouchi2022) extend the ISMIP6 simulations for the AIS with SICOPOLIS until the year 3000, assuming a sustained late-21st-century climate beyond 2100. Compared to the uncertain response projected over the ISMIP6 period, a radically different picture emerges, demonstrating that the consequences of the high-emissions scenario RCP8.5/SSP5-8.5 (SSP: Shared Socioeconomic Pathway) are much greater in the long term even if no further climate trend is applied beyond 2100.
The response of the GrIS to longer-term climate change has also been investigated. In addition to an ensemble of projections for the 21st century with their higher-order ice-sheet model, the study by Fürst and others (Reference Fürst, Goelzer and Huybrechts2015) conducted projections until 2300 for the two low-emission scenarios RCP2.6 and RCP4.5, forced by selected CMIP5 GCMs. Vizcaino and others (Reference Vizcaino2015) carried out simulations until 2300 with a coupled ECHAM5.2/MPI-OM/SICOPOLIS model for the pathways RCP2.6, RCP4.5, and a modified RCP8.5 with a 4 × CO2 limit. Calov and others (Reference Calov2018) drove an extended version of SICOPOLIS (IGLOO) with RCP4.5 and RCP8.5 surface temperature and surface mass-balance anomalies created by the regional climate model MAR with boundary conditions from simulations with three CMIP5 GCMs. Similar to our approach, prolongation of the climatic forcing beyond 2100 was done by assuming no further warming trend. Aschwanden and others (Reference Aschwanden2019) used projections based on four CMIP5 GCMs until 2300 for RCP2.6, RCP4.5 and RCP8.5, extrapolated until 3000, to force the ice-sheet model PISM. Their climatic forcing was processed by the regional climate model HIRHAM5, delivering the surface-temperature anomaly as the main driver, from which precipitation changes were parameterized, and runoff was computed by a positive-degree-day (PDD) method. Van Breedam and others (Reference Van Breedam, Goelzer and Huybrechts2020) projected the response of the GrIS and AIS 10000 years into the future with the Earth system model of intermediate complexity LOVECLIMv1.3 (including the ice-sheet model AGISM), forced by the extended concentration pathways ECP2.6, 4.5, 6.0 and 8.5 until 2300 and zero emissions thereafter. In the present study, we transfer the approach by Chambers and others (Reference Chambers, Greve, Obase, Saito and Abe-Ouchi2022) to the GrIS. The objective is to assess its long-term response to late-21st-century climatic conditions for the full ensemble of ISMIP6 climate forcings, which consists of 14 scenarios from ten different CMIP5 and CMIP6 GCMs.
2. Methods
The main tool used for this study is the ice-sheet model SICOPOLIS. We apply it to the GrIS with hybrid shallow-ice–shelfy-stream dynamics (Bernales and others, Reference Bernales, Rogozhina, Greve and Thomas2017), a Weertman–Budd-type sliding law tuned separately for 20 different regions (Greve and others, Reference Greve, Chambers and Calov2020b), and ice thermodynamics treated by the one-layer melting-CTS enthalpy scheme (CTS: cold-temperate transition surface; Blatter and Greve, Reference Blatter and Greve2015; Greve and Blatter, Reference Greve and Blatter2016). The horizontal resolution is 5 km. In the vertical, we use terrain-following coordinates (sigma transformation) with 81 layers in the ice domain and 41 layers in the thermal lithosphere layer below. For details on the set-up, the initialization procedure by a paleoclimatic spin-up, comparisons between the simulated and observed ice thickness and surface velocity for our initialization year 1990, as well as the historical run (hist) that bridges the gap between 1990 and the start date of the projections in January 2015 by employing MIROC5/RCP8.5 surface mass balance (SMB) and surface temperature (ST) forcing, we refer to Greve and others (Reference Greve, Chambers and Calov2020b).
Following the ISMIP6 protocol, climate forcing from 2015 until the end of 2100 has an atmospheric and an oceanic component. The atmospheric forcing consists of a 1960–1989 reference climatology, plus space-time-dependent anomalies for SMB and ST. These were derived from a systematic sampling of CMIP5 GCMs that reflects their spread in future projections (Barthel and others, Reference Barthel2020), while CMIP6 GCMs were added on the basis of availability only (Payne and others, Reference Payne2021). All GCM forcings were downscaled to the GrIS surface with the regional climate model MAR v3.9.6 (Fettweis and others, Reference Fettweis2017; Delhasse and others, Reference Delhasse2020). Although MAR uses a static GrIS, it also provides vertical gradients of SMB and ST, thus allowing to include SMB–height and ST–height feedbacks in the ice-sheet simulations (Franco and others, Reference Franco, Fettweis, Lang and Erpicum2012; Nowicki and others, Reference Nowicki2020). The oceanic forcing is based on a retreat parameterization for tidewater glaciers, forced by MAR runoff and ocean temperature changes specified for seven ice–ocean sectors around Greenland (Slater and others, Reference Slater2019, Reference Slater2020).
For the period from 2101 until the end of 3000, we extend the simulations in a similar way than Chambers and others (Reference Chambers, Greve, Obase, Saito and Abe-Ouchi2022) do for the AIS. For every year of this extended period, the atmospheric forcing (SMB, ST, vertical gradients) for the 10-year interval 2091–2100 is randomly sampled such that no further trend is applied, but some inter-annual fluctuations remain (similar to Calov and others, Reference Calov2018). The oceanic forcing (prescribed retreat maps) does not show any notable year-to-year fluctuations, so we simply keep it fixed at 2100 conditions.
An overview of our extended ISMIP6 experiments is given in Table 1, and the magnitude of the atmospheric forcing is shown in Supplementary Table S1. Twelve experiments are for the 21st-century unabated warming pathway RCP8.5 (CMIP5) / SSP5-8.5 (CMIP6), and two are for the reduced emissions pathway RCP2.6 (CMIP5) / SSP1-2.6 (CMIP6) that is largely in line with the commitments of the Paris Agreement (maintaining the global mean temperature well below a 2°C increase above pre-industrial levels). In two experiments, the impact of different sensitivities of the retreat parameterization due to oceanic forcing (‘high’ and ‘low’ vs the normal, ‘medium’ sensitivity, thereby exploring the uncertainty of the parameterization; Slater and others, Reference Slater2019, Reference Slater2020) is tested. In addition, a projection control run for the period 2015–3000 (ctrl_proj) employs constant climate conditions based on a 1960–1989 climatology and no explicit oceanic forcing.
See Nowicki and others (Reference Nowicki2020) for references for the GCMs.
3. Results
The simulated mass change of the GrIS, expressed as a sea-level contribution, and ice area are shown in Figure 1. For the control run ctrl_proj, the ice sheet remains nearly stable, showing a slight mass gain of 6.4 mm SLE and area loss of 4.7 × 103 km2 during the 986 years model time, which is of the order of permilles of the present-day values. For all future projections, the ice sheet keeps losing both mass and extent over the entire period. (Values of the resulting sea-level contribution by 2100, 2300 and 3000 relative to ctrl_proj are provided in the Supplementary Table S2.) The largest rate of change occurs typically around the year 2100, beyond which it slows down to some extent; however, without reaching or coming close to a new steady state. This demonstrates that the committed mass loss due to 21st-century climate change extends way beyond the 21st century and impacts the ice sheet on a much longer time scale. Corroborating the findings for the 21st century (Goelzer and others, Reference Goelzer2020; Greve and others, Reference Greve, Chambers and Calov2020b), the GrIS responds much more strongly to the ensemble of RCP8.5/SSP5-8.5 simulations than to the two RCP2.6/SSP1-2.6 simulations. By the year 3000, the mass loss amounts to 1.79 ± 0.80 m SLE (mean ± 1-sigma range) for RCP8.5/SSP5-8.5, while it is limited to 0.28 ± 0.12 m SLE for RCP2.6/SSP1-2.6.
The percentage area loss (Fig. 1b) is similar to the mass loss. However, it is striking that the variability in area varies significantly between the different simulations. This results from the different variability of the SMB forcing and affects mainly the thin, near-margin parts of the ice sheet, which do not contribute much to the total ice mass. Therefore, the variability is not paralleled in Figure 1a.
The influence of the ice retreat due to oceanic forcing is explored by Exps. 5, 9, 10 (MIROC5/RCP8.5 with ‘medium’, ‘high’ and ‘low’ sensitivity, respectively). The results are shown by the olive lines and olive-shaded regions in Figure 1. By 3000, the simulated mass loss is $1.62^{ + 0.051}_{ - 0.031}\, {\rm m\, SLE}$. Thus, the uncertainty due to these three calibrations is very small in the long range. Relative to the uncertainty due to the different climate forcings, it is more pronounced for the 21st century (Greve and others, Reference Greve, Chambers and Calov2020b). This is because the continued retreat of the ice sheet decreases its contact with the ocean, so that the oceanic forcing plays a smaller role in the longer term.
As reported by Greve and others (Reference Greve, Chambers and Calov2020b) and Payne and others (Reference Payne2021), for both the 21st-century RCP8.5/SSP5-8.5 and RCP2.6/SSP1-2.6 pathways, the CMIP6 climate models produce a larger response of the ice sheet than the CMIP5 ones. While the significance of this statement is limited in the case of RCP2.6/SSP1-2.6 (only one experiment each), it is more robust for RCP8.5/SSP5-8.5, where the ensemble contains eight and four experiments forced by CMIP5 and CMIP6 models, respectively. By 3000, the mean mass loss for the four CMIP6 SSP5-8.5 experiments is 2.73 m SLE, and the maximum value from Exp. B4 (CESM2/SSP5-8.5) is as large as 3.54 m SLE, almost 50% of the entire present-day ice mass.
We now discuss in more detail the results of Exp. 5 (MIROC5/RCP8.5), which was already focused on in the original ISMIP6-Greenland study (Goelzer and others, Reference Goelzer2020). Mainly due to its large SMB forcing, it produces, along with Exp. A1 (IPSL-CM5A-MR/RCP8.5), the strongest mass loss among the CMIP5 forcings, while the mass loss is about average for our combined CMIP5/CMIP6 ensemble. Figure 2 shows the components of the global mass balance (integrated over the ice sheet, all counted as positive for mass gain): surface mass balance (SMB), basal mass balance (BMB), calving and ice volume change (dV/dt). On a mean-annual basis, the residual, Res = |SMB + BMB + Calving − dV/dt|, is always <$10^{6}\, {\rm m^{3}\, a^{-1}}$. This is five to six orders of magnitude smaller than the typical range of values in the figure, so that the model conserves mass very well (see also Calov and others, Reference Calov2018).
As already stated above, the ice sheet keeps losing volume (∝ mass) over the entire period, with maximal rates of change occurring shortly before the year 2100. The SMB is initially positive, but changes its sign in the second half of the 21st century and stays negative beyond that. Calving into the surrounding ocean peaks during 2080–2085, when it contributes approximately the same amount to ice volume loss than negative SMB. After that, calving decreases continuously due to ice-sheet retreat from the coast (loss of ocean contact) and becomes almost negligible towards the end of the 3rd millennium. The decrease is likely accelerated by a limitation of the oceanic forcing approach: the fixed retreat mask beyond 2100 does not follow the retreating ice margin further, making it ineffective as the ice sheet recedes beyond its reach. The negative SMB shows a slightly decreasing trend after 2100 because, as the ice sheet shrinks, less area is available for further melting. Both effects together cause the magnitude of dV/dt (loss rate of ice volume) to decrease gradually, which results in the slight flattening of the curves in Figure 1 (the same mechanisms are effective for the other experiments). BMB is relatively small over the entire model time. The inter-annual variability of the volume change is due to that of the SMB, which reflects the variability of the atmospheric forcing.
Figure 3 shows snapshots of the ice thickness and surface velocity for Exp. 5 for the initial year 2015 and the final year 3000. Comparing the thickness distributions demonstrates nicely that the ice sheet retreats from the coast almost all around its perimeter, and contact to the ocean is very limited by the end of the simulation, which entails the low calving rates mentioned above. By contrast, the ice sheet does not suffer much change in its interior parts north of ~ 68°N. The large-scale pattern of the ice flow and the organization of the ice sheet into major drainage basins remain largely intact. However, on the regional scale, more pronounced changes occur. The fast-flowing outlet glaciers in south-western Greenland disappear entirely due to the extreme retreat in this area. The north-western outlet glaciers, including Petermann Glacier, also slow down substantially. The central-western Jakobshavn Ice Stream loses its clear delimitation to the surrounding glaciers, but remains an area of fast-flowing ice. The major features in East Greenland, e.g. the North-East Greenland Ice Stream, Kangerdlugssuaq and Helheim glaciers, are less affected and remain well identifiable.
4. Discussion and conclusion
The future climate simulations carried out in this study for the GrIS over the 3rd millennium confirm and continue the trends that were reported by ISMIP6-Greenland for the 21st century (Goelzer and others, Reference Goelzer2020; Greve and others, Reference Greve, Chambers and Calov2020b; Payne and others, Reference Payne2021). The response of the ice sheet is mainly governed by a negative SMB due to increased surface melting near the ice margin. Marine-terminating glacier retreat, caused by increased oceanic thermal forcing and increased meltwater runoff, constitutes a further negative contribution to the total mass balance, but becomes less important in the longer term. Under the unabated warming pathway RCP8.5/SSP5-8.5, this leads to a severe mass loss during the 3rd millennium, while the loss is much smaller under the reduced emissions pathway RCP2.6/SSP1-2.6. Results obtained with forcings from the newer CMIP6 climate models consistently produce larger mass losses than those obtained with the older CMIP5 models, for SSP5-8.5 in the range of a ~ 25–50% loss of the present-day ice mass (and area) by 3000. For comparison, Aschwanden and others (Reference Aschwanden2019) reported a mass loss of ~ 75–100% by 3000 for their ensemble of RCP8.5 simulations, for which a warming trend is assumed to continue until 2500. Efficient climate change mitigation during the next decades is therefore crucial for limiting the contribution of the GrIS to long-term sea-level rise.
As for interpreting the stronger response of the GrIS to CMIP6-derived forcings compared to CMIP5-derived ones, the different strategies of sampling the respective GCM ensembles for ISMIP6-Greenland must be considered. As we mentioned in Section 2, the six CMIP5 GCMs were chosen systematically, whereas the four CMIP6 GCMs were the only ones available for downscaling at the time. Subsequent analysis of the entire CMIP6 model ensemble revealed that the results cluster around two groups of climate sensitivities (global mean temperature response to doubled atmospheric CO2), and the four models that were available for ISMIP6-Greenland all fall in the high-sensitivity group (Meehl and 7 others, Reference Meehl2020; Payne and others, Reference Payne2021). Future work should therefore aim at conducting additional simulations with a more representative sampling of the CMIP6 GCMs.
Our study is limited to investigating the impact of a sustained late-21st-century climate (without imposing a further trend beyond 2100) on the GrIS. However, climate change is projected to continue beyond 2100 (e.g. Bakker and others, Reference Bakker2016; Lyon and others, Reference Lyon2022), with potentially even more devastating effects on the GrIS than reported here, plus a significant decay of the AIS in the long term (Van Breedam and others, Reference Van Breedam, Goelzer and Huybrechts2020). Further, the unidirectional coupling approach (climate model → ice-sheet model) employed by ISMIP6, and thus here, lacks a detailed accounting of feedbacks of the changing ice sheet on the climate. As we explained in Section 2, the climate forcing for Greenland includes vertical gradients of the surface mass balance and surface temperature. Therefore, the changing ice-sheet geometry acts back on these climatic forcing fields. However, the linearized approach was derived for small perturbations of the present-day state only, and it cannot be validated for large changes of the ice sheet. This shortcoming becomes more severe in our simulations over almost a millennium compared to the 86-year scope of ISMIP6, and it is not possible to judge a priori whether it rather leads to an over- or underestimation of the simulated mass loss. Using the PDD method for modelling runoff, like in the study by Aschwanden and others (Reference Aschwanden2019), allows in principle to handle arbitrarily large changes of the ice sheet. PDD is also easy to implement in an ice-sheet model, which adds to its appeal. However, it is a grossly simplifying parameterization of the complex energy balance at the surface of an ice sheet and difficult to calibrate for climatic conditions other than present-day ones (e.g. Bauer and Ganopolski, Reference Bauer and Ganopolski2017), so that we do not consider it superior to the ISMIP6-type approach of our study.
Future work in the direction of long-term simulations of ice-sheet response to climate change should aim at employing GCM projections beyond 2100 and improving the representation of feedback processes. The ultimate solution would be to carry out such simulations in a fully coupled way, with the ice-sheet model integrated in the GCM. This approach has been pursued (e.g. Vizcaino and others, Reference Vizcaino2015; Sellevold and others, Reference Sellevold2019; Gregory and others, Reference Gregory, George and Smith2020; Muntjewerf andand others, Reference Muntjewerf2020; Muntjewerf and others, Reference Muntjewerf2020); however, fully coupled simulations are demanding and computationally expensive, which makes it difficult to run large ensembles, involving many different climate and ice-sheet models, over long time scales and at adequate resolution. Intermediary, more manageable solutions may consist of involving snapshots of climate-model results combined with more refined parameterizations for the climatic forcing, similar to the approach by Abe-Ouchi and others (Reference Abe-Ouchi2013) for the paleo-glaciation of the Northern Hemisphere.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2022.9.
Code and data availability
SICOPOLIS is free and open-source software, available through a persistent Git repository hosted by the Alfred Wegener Institute for Polar and Marine Research (AWI) in Bremerhaven, Germany (Greve and SICOPOLIS Developer Team, Reference Greve2021). Detailed instructions for obtaining and compiling the code are at http://www.sicopolis.net (last access: 2 November 2021). The output data produced for this study are available at Zenodo, https://doi.org/10.5281/zenodo.5880517.
Acknowledgements
We thank the two anonymous reviewers, the Scientific Editor Frank Pattyn and the Chief Editor Hester Jiskoot for constructive remarks and suggestions that helped to improve the manuscript. We thank Jorge Bernales (MARUM Bremen), Reinhard Calov (PIK Potsdam), Takashi Obase (University of Tokyo) and Fuyuki Saito (JAMSTEC Yokohama) for their recent contributions to the development of the SICOPOLIS model, and Ayako Abe-Ouchi for fruitful discussions about ice-sheet and climate modelling. Some colour schemes of our figures were taken from Paul Tol's (SRON Netherlands Institute for Space Research) online resource at https://personal.sron.nl/̃pault. We thank the Climate and Cryosphere (CliC) effort, which provided support for ISMIP6 through sponsoring of workshops, hosting the ISMIP6 website and wiki, and promoting ISMIP6. We acknowledge the World Climate Research Programme, which, through its Working Group on Coupled Modelling, coordinated and promoted CMIP5 and CMIP6. We thank the climate modelling groups for producing their model output and making it available; the Earth System Grid Federation (ESGF) for archiving the CMIP data and providing access to it; the University at Buffalo for ISMIP6 data distribution and upload; and the multiple funding agencies who support CMIP5, CMIP6 and ESGF. We thank the ISMIP6 steering committee, the ISMIP6 model selection group and ISMIP6 dataset preparation group for their continuous engagement in defining ISMIP6. This is ISMIP6 contribution No. 26.
Financial support
Ralf Greve and Christopher Chambers were supported by Japan Society for the Promotion of Science (JSPS) KAKENHI Grant No. JP17H06323, and by a Leadership Research Grant (Category 2) of Hokkaido University's Institute of Low Temperature Science. Ralf Greve was supported by JSPS KAKENHI Grant No. JP17H06104, and by the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) through the Arctic Challenge for Sustainability project ArCS II (program grant number JPMXD1420318865).