Ice dynamics will remain a primary driver of Greenland ice sheet mass loss over the next century

The mass loss of the Greenland Ice Sheet is nearly equally partitioned between a decrease in surface mass balance from enhanced surface melt and an increase in ice dynamics from the acceleration and retreat of its marine-terminating glaciers. Much uncertainty remains in the future mass loss of the Greenland Ice Sheet due to the challenges of capturing the ice dynamic response to climate change in numerical models. Here, we estimate the sea level contribution of the Greenland Ice Sheet over the 21st century using an ice-sheet wide, high-resolution, ice-ocean numerical model that includes surface mass balance forcing, thermal forcing from the ocean, and iceberg calving dynamics. The model is calibrated with ice front observations from the past eleven years to capture the recent evolution of marine-terminating glaciers. Under a business as usual scenario, we find that northwest and central west Greenland glaciers will contribute more mass loss than other regions due to ice front retreat and ice flow acceleration. By the end of century, ice discharge from marine-terminating glaciers will contribute 50 ± 20% of the total mass loss, or twice as much as previously estimated although the contribution from the surface mass balance increases towards the end of the century. Changes in ice dynamics in the Greenland Ice Sheet are projected to contribute to significantly more mass loss than previously thought, suggesting that proper calibration and accurate representation of ice-ocean interactions are critical.

T he Greenland ice sheet has been losing mass over the last several decades 1 and contributing to sea level at a rate of0 .8 mm/yr since 2002 [2][3][4] . Although the increase in melting at the ice sheet surface has been accounting for a large share of this mass loss in recent years, glacier dynamics has ultimately been responsible for 66% of the total mass loss over the last 46 years 1 . The physical processes responsible for the increase in glacier dynamics along the ice margins include grounding line retreat due to ice thinning, forced-grounding line retreat due to melt by the ocean, and enhanced calving. A leading cause of the changes in ice dynamics is thought to be an increase in subsurface ocean temperature due to an enhanced transfer of ocean heat into the fjords from the North Atlantic Ocean [5][6][7] . As marineterminating glaciers became exposed to warmer subsurface waters and enhanced subglacial meltwater discharge, ocean-induced undercutting at glacier termini and melting under floating ice shelves have considerably increased, which presumably lead to widespread glacier retreat and ice flow acceleration 5,[8][9][10] . Glacier dynamics will continue playing a major role in glacier mass loss if ocean and air temperatures stay high or keep increasing in the future. To project future changes in Greenland realistically, it is essential to employ numerical models validated by decades of precise observations of glacier evolution that include a realistic description of ice-ocean interaction mechanics leading to ice shelf melt and terminus undercutting as well as rapid iceberg calving mechanisms.
Several studies have been dedicated to projecting the future sea level contribution of Greenland with a focus on glacier dynamics [11][12][13][14] . While these studies include modeled changes in surface mass balance (SMB), they often use a coarse grid resolution or simplified physics, and do not calibrate their models at the glacier scale and with decades of observations. Many studies do not have the required spatial and temporal resolution to represent individual marine-terminating glaciers of Greenland with fidelity, which is a major limitation if ice dynamics is a major control on mass loss. Models operating at a coarse resolution without calibrating the retreat of individual glaciers might produce overly conservative estimates of mass loss 11,13 because the glacier response to changes in ice-ocean interactions is incompletely investigated. These estimates are currently the only projections mentioned in the Intergovernmental Panel on Climate Change (IPCC) Assessment Report 15 . According to these projections, the future contribution of Greenland to sea level will be dominated by changes in SMB and the effect of ocean thermal forcing will remain small 14 , which is not consistent with observations of accelerated retreat since the 1990s 1,16 . Although some studies relied on flowband models 11 or have employed higher horizontal resolution models (e.g.,~400 m) 14 to include individual outlet glaciers, they remained limited to a few outlet glaciers without resolving critical details in the bed topography geometry, or the models did not capture recent changes in these glaciers. Uncalibrated models yield large uncertainties in ice discharge for both hindcast and forecast simulations.
Here, we model the response of the Greenland ice sheet to oceanic and atmospheric forcings to investigate its evolution with a spatial resolution as high as 200 m along the coast (Fig. 1). We initialize the high-resolution model with data collected going back to year 2007. After the model initialization, we run the model forward from 2007 to 2017 and calibrate the calving parameterization for each glacier drainage basin so that the model is consistent with observed velocities, changes in ice front positions, and ice thinning rates (see "Methods"). This calibration process captures not only the current state of the ice sheet but also the trend in mass loss for the first decade, which improves its reliability for short-term projections 17 . This approach has been used in a recent modeling of Northwest Greenland 18 but not at the scale of the entire ice sheet. We use a 3D higher-order model (HO) 19,20 to account for both membrane stresses and vertical shear, which is adequate to model both fast outlet glaciers and the slower moving regions inland of the ice sheet. We simulate the future sea-level contribution of Greenland until year 2100 using output products from general circulation models (GCMs) from the Coupled Model Intercomparison Project Phase 5 (CMIP5) and Phase 6 (CMIP6), which were selected by the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6) 21 . Contrary to the data provided by ISMIP6, we use a temporal resolution of one month for these climatic forcings, instead of one year, to capture the seasonal cycle in SMB and ocean temperatures, which is important in the context of glacier de-stabilization 22,23 .

Results and discussion
Marine-terminating glaciers. The model includes 215 marineterminating glaciers, including six glaciers from northern Greenland which have an extensive floating ice shelf. These 215 marine-terminating glaciers control more than 90% of the ice sheet discharge and 85% of the observed mass loss from Greenland (ice discharge and SMB) 1 . By calibrating the stress threshold for calving, we are able to match the observed retreat distance of most glaciers (see Supplementary Tables 1-6 and Supplementary  Fig. 1). The ice fronts of 115 glaciers are calibrated within 1 km from the observed retreat. For 33 additional glaciers, we overestimate the retreat of 8 glaciers by up to 3 km and underestimate the retreat of 25 glaciers by up to 24 km from year 2007 to 2017. Finally, the dynamics of the ice front is not captured for 67 glaciers because of poorly constrained bed topography. For these glaciers, we keep the ice front fixed in time to avoid nonphysical behavior during the simulations. Interestingly, the 115 wellsimulated glaciers controlled 79% of the ice discharge in year 2007, versus 3% for the glaciers with overestimated retreat, 13% for the 25 glaciers with underestimated retreat, and 5% for the glaciers with poorly-known bathymetry. We note here that if we use a earlier bathymetry/topography map of Greenland (e.g., bed topography maps from the previous studies 24,25 or RTopo-2 26 ) to repeat the same exercise, many glaciers are impossible to calibrate because some glaciers were land terminating glaciers in the old map instead of marine terminating glaciers (not shown here). Additionally, in previous datasets, the bed topography is generally too shallow and the fjord bathymetry is usually non-existent 27 .
For each marine-terminating glacier, we project the changes in ice mass, cumulative SMB and ice discharge from year 2007 to 2100 ( Fig. 2 Fig. 2 indicates the rate of mass changes, SMB and ice discharge per year. Our simulations suggest that many outlet glaciers continue to experience ice front retreat under all future climate change scenarios. The pattern of modeled retreat varies significantly from one glacier to the next, which is consistent with observations and prior modeling studies 14,28 . For example, in central Greenland, Eqip Sermia glacier, and Jakobshavn Isbrae respond differently to the same climate forcings ( Fig. 2a and b). Under MIROC5 RCP8.5 scenario, the mass change of Eqip Sermia glacier is mainly controlled by SMB because ice discharge decreases once the ice front stabilizes 5 km upstream of its current position. For Jakobshavn Isbrae, the model suggests that ice discharge dominates the mass loss over the entire century as the ice front continues to retreat along an overdeepening in bed topography upstream of the current ice front position until year 2070, after which the grounding line will temporarily stabilize about 56 km upstream of its current position at a location where the bed elevation rises again. Our model does not capture the recent slow down of Jakobshavn Isbrae because our modeled ocean thermal forcing does not capture a cooling of the ocean waters over the past two years, which has not been seen for the past nearly 20 years 29 . The main factor responsible for these contrasting behaviors between glaciers is the bed geometry but the ocean thermal forcing also could affect the retreat of glaciers as shown in the recent deceleration of Jakobshavn Isbrae.
It is well known that bed topography plays a critical role in determining stable positions of calving fronts and in turn on the glacier mass loss 9,30 . For two glaciers in different regions under the RCP8.5 scenario ( Fig. 2c and d), Narsap Sermia in southwest (SW) and Kakivfaat Sermiat in northwest (NW), ice discharge initially dominates the changes in mass balance during the first 20-30 years, until the ice front stabilizes upstream, which leads to a glacier slow down. After 2060, SMB decreases significantly to dominate the glacier mass loss in the last 30 years of the simulation. The primary factors of mass loss, therefore, changes with time, depending on the retreat and the dynamics of individual glaciers. While the extent of the retreat for these two glaciers differs significantly, both retreats are controlled by bed topography.
A deep bed topography does not imply a fast glacier retreat as illustrated by Kangerlussuaq Gletscher in southeast Greenland under MIROC5 RCP8.5 scenario (Fig. 2e). This glacier loses mass continuously due to ice discharge, but the SMB is projected to remain largely positive over the region, which compensates for the mass loss from ice discharge. This positive SMB effect continues until the end of the simulation, keeping the ice front at its current position instead of retreating in a deep trough. The shape of the bed topography provides clues about the potential stability and vulnerability of a given glacier, but several factors such as lateral shear or buttressing may delay or slow the retreat, even within regions of deep bed topography or retrograde bed slopes 31 . A numerical model is therefore necessary to determine whether a glacier will retreat or not for a given scenario and it is important to evaluate the uncertainties of the model associated with uncertainties in climate forcing and geometry.
Regional changes. We divide Greenland into six large regions ( Fig. 1 and calculate the cumulative mass balance, SMB, and ice discharge for each region over the simulation period (Fig. 3). We find a better agreement between modeled and observed trends in mass loss for NW (−71.1 Gt/yr of observation versus −58.5 Gt/yr of modeled mass loss) and N (−25.5 Gt/yr of observation versus −21.6 Gt/yr of modeled mass loss) than the east regions of Greenland (SE and NE). Both regions gain mass in the model during the calibration period while they actually lost about 952 Gt of ice between 2007 and 2017 1 . This region remains affected by residual uncertainties in bed topography of its many glaciers, which yields a weak agreement between model and observations, as in previous studies 14,32 . We keep the ice front positions of many glaciers in this region fixed to stabilize the model, with the consequence that the regional mass loss is underestimated. If we have more accurate bed topography data with deeper fjords to match those glaciers, we could have more mass loss from this region by ice discharge with retreating ice fronts. Figure 3 shows the cumulative change in ice mass along with SMB and ice discharge from 2007 to 2100 for the six different regions forced by MIROC5 under the RCP8.5 scenario. SW loses 9581 Gt of ice by 2100, or 26.5 mm of sea level rise. The primary driver for the mass loss is a large decrease in SMB starting in 2060. The ice discharge from 15 marine-terminating glaciers decreases after 2060 in response to the decrease in SMB, which reduces the ice flux at the terminus. In CW, the simulation shows that the glaciers will raise sea level by 26.1 mm (i.e., 9462 Gt of ice). The overall SMB is positive until 2060 and then decreases until year 2100. Most of the mass loss (65%) is due to the extensive retreat of Jakobshavn Isbrae glacier. The remaining glaciers lose 3360 Gt of ice by the end of the century. The NW region has a similar pattern and loses 10,809 Gt of ice (29.8 mm sea level equivalent, SLE) by 2100 with ice discharge from over 60 glaciers contributing 66% of the total signal. The N region will contribute 10.9 mm to sea level by 2100, with 33% from Humboldt glacier alone, principally from SMB. NE is projected to lose 2237 Gt of ice mass (or 6.2 mm SLE) by 2100, principally through enhanced discharge from Nioghalvfjerdsfjorden (79North) and Zachariae Isstrøm glaciers. SE shows a similar pattern to NE: a positive SMB of 8594 Gt until 2100 versus 141 Gt/year mass loss via ice discharge. This result may be viewed as an underestimate of the mass loss given that SE has been one of the largest contributor to sea level rise in the past 46 years, including in the recent decades.
Overall, in the simulations forced by CMIP5 RCP8.5 and CMIP6 SSP585, the model shows that all sectors of Greenland will continue to lose mass and the rate of mass loss increases during this century, including in the north. This mass loss is due to a decrease in SMB caused by an increase in surface temperature and a decrease in surface albedo, which leads to higher melt rates at the surface. Additionally, ice discharge will either be sustained at high levels or increase for all marine-terminating glaciers. The increase in ice discharge is a major driver of the mass loss over the entire time period, except in SW, where about 60% of glaciers are land terminating. In the recent multi-model study of the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6) 21 , SW was highlighted as the largest contributor to sea level. In contrast, our model reveals that sectors further north, i.e., CW and NW, where over 90% of glaciers are marine terminating, will contribute almost equally to sea level rise under all scenarios. Models participating in ISMIP6 use a parameterized ice front retreat or a coarse grid resolution, leading to a weak dynamic response of CW and NW sectors under RCP8.5. Our calibrated model, however, projects that these two sectors will remain major contributors to the ice sheet mass change over this century, consistent with the past decades of observations. The Northern sectors of Greenland (N and NE) will continue to lose mass due to a sustained anomaly in ice discharge combined with a decrease in SMB (Fig. 3). The simulated sea level contribution from these regions is still underestimated because of the weaker calibration of the model in this region, the limitations of our ice shelf melt parameterization, and residual uncertainties in bed topography. For instance, the model does not capture the observed ice front retreat for 21 out of 32 glaciers in N and NE, 4 of which have a floating ice shelf e.g., Petermann and Storstrømmen. The ice fronts of 13 glaciers was fixed during the simulations to prevent the ice front from behaving unrealistically. Better observations of bed topography and ocean thermal forcing will be needed in these regions to estimate future changes in mass balance of N and NE more realistically.
Sea level contribution from Greenland. Our simulations project that, overall, the Greenland ice sheet contribution to sea level by the end of the century will range from 79 to 147 mm under RCP8.5 climate forcing scenarios (Fig. 4), which is on the high end or higher than ISMIP6 estimates 21 . Based on the latest CMIP6 SSP585 climate forcings, the simulated Greenland sea level contribution until 2100 ranges from 94 to 167 mm, which is at or above the range from CMIP5 RCP8.5 simulations. This is due to the increased future warming in many CMIP6 models compared to CMIP5 models 33 . Based on CMIP5 RCP4.5, the Greenland ice sheet will raise sea level between 54 and 79 mm by the end of this century. Overall, we find that the rate of mass loss will continue to increase but the rate of increase depends on the climate forcing applied. Interestingly, the simulations forced by CMIP6 SSP585 forcings (CESM2, CNRM-CM6, CNRM-ESM2, and UKESM1-CM6) lose mass at a higher rate after about 2080 compared to the ones forced by CMIP5 RCP8.5 data (MIROC5, CanESM2, and NorESM1), which we attribute to a larger decrease in SMB (Fig. 4). Figure 4 shows the partitioning of the mass loss between ice discharge and SMB under CMIP5 RCP8.5 (Fig. 4b) and CMIP6 SSP585 (Fig. 4c) simulations. The ice loss from changes in ice discharge continues until 2100 at a similar rate for both scenarios. Changes in SMB are small until 2050 and then increase significantly around year 2060. Although changes in SMB increase rapidly, ice discharge remains a major contributor to the total mass loss over this century, accounting for 38-70% for CMIP5 RCP8.5 and 22-56% for CMIP6 SSP585, respectively, in 2100. These estimates are significantly larger than from the previous studies 14,21,34 and suggests that it is essential to include changes in ice dynamics in projections.
The mass loss from changes in ice discharge is similar under different CMIP5 and CMIP6 models, while the mass loss from SMB varies largely between CMIP5 and CMIP6 models. We attribute this behavior to the fact that many marine-terminating glaciers engage in a state of ice front retreat triggered by ocean warming, and then proceed with their ice retreat independent of SMB as the retreat is primarily driven by bed topography and secondarily by the initial increase in ocean thermal forcing. For example, under both RCP4.5 and RCP8.5 scenarios, Jakobshavn Isbrae glacier continues to retreat along its deep trough over the same distance (56 km), at a similar rate, once the ice front is Conclusions. We projected the future sea-level contribution of Greenland with a high-resolution model that accounts for changes in ice front of individual marine-terminating glaciers, iceberg calving dynamics, and thermal forcing from the ocean. Under CMIP5 RCP8.5 and CMIP6 SSP585 emission scenarios, Greenland will contribute 79-167 mm to sea level by 2100, which is higher than recent Greenland-wide modeling estimates 21 . We find that the ice discharge accounts for 22-70% of the total mass loss until 2100, or 1.5 to 3 times larger than in previous studies 14,21,34 although the contribution from the surface mass balance to mass loss increases towards the end of the century. This higher contribution to sea level is caused by a large mass loss from NW and CW Greenland, where extensive retreat of marineterminating glaciers are projected by our model due to changes in ocean temperature. We note that the replication of glacier changes in NW and CW was only possible due to recent advances in bathymetry and bed topography mapping. In SE, N, and NE Greenland, our projections still underestimate the glacier loss because of incomplete bathymetry and bed topographic constraints. This limitation underlines the fundamental role of deriving proper boundary conditions and climate forcing, including the ocean, in order to replicate the last decades of observations and in turn improve the projections of sea level rise from Greenland and minimize the risk of significantly underestimating these changes.

Methods
Model setup. We use the Ice-sheet and Sea-level System Model (ISSM) 35 to model Greenland tidewater glaciers. Our model is based on 3-D higher-order approximation 19,36 , with subelement grounding line parameterizations 37 and level set-based moving boundaries 38 . The horizontal mesh resolution varies from 200 m to 20 km depending on observed surface velocity (Fig. 1). The mesh is vertically extruded into 4 layers which are sufficient for the model without a transient thermal field 39 . The ice temperature is kept constant as it has been shown that changes in surface temperature do not affect simulations over the time scales considered here 40 . We divide the Greenland ice sheet into six domains (Fig. 1) to reduce computational cost compared to running simulations with one large domain. For each domain, we perform the same model initialization approach. The surface and bed geometry are from BedMachine Greenland version 3 41 . To initialize the model, we use the inversions to infer the basal friction coefficients under grounded ice and ice viscosity parameters on floating ice 42 based on 2007-2008 surface velocities 43 . The ice viscosity, μ, is defined by Glen's law 44 : where B is ice viscosity parameter, _ ε e is the effective strain rate, and n is the Glen's law exponent set to 3. The friction in this study follows Budd friction law 45 with a linear relationship with the effective pressure, N, defined as the difference between ice overburden pressure and basal water pressure: where k is the friction coefficient. The exponent s represents the relation between velocity (v) and basal friction (τ b ) and is chosen to be 1 in this study. The inferred friction coefficient is kept constant in time during all simulations while the effective pressure changes in time. The modeled velocities for each glacier at initial state are shown in the Supplementary Data 1.
Model calibration. To calibrate our model, we try to match observed ice front retreat from 2007 to 2017 following previous studies 18,46 . To track the dynamic motion of the ice front, we rely on the level set method 38 , where the velocity of the calving front is defined as follows: where v is the ice velocity vector, c is the calving rate, _ M is the undercutting rate at the calving face and n is a unit normal vector that points outward from the ice domain. c and _ M are assumed to be independent. We use the undercutting parameterization 47 to estimate _ M and specifically optimize c to calibrate modeled ice front. We use von Mises tensile stress calving law 48 that provides the calving rate at each time step 46 . The calving rate is assumed to be proportional to the tensile von Mises stress,σ, which only accounts for the tensile component of the stress in the horizontal plane: where σ max is a stress threshold that is calibrated, B is the ice viscosity parameter, n = 3 is Glen's exponent, and_ ε e is the effective tensile strain rate mentioned in previous studies 46,48 . The stress threshold to be calibrated does not depend on the mesh resolution but boundary condition that varies from glacier to glacier. We determine the stress threshold value for each glacier by simulating the ice front changes during calibration period, between 2007 and 2017, and comparing the modeled retreat distance to observed retreat [49][50][51] . We manually adjust the stress threshold for each basin to best capture the observed variations in ice front positions. We limit the range of the stress threshold between 150 kPa and 3100 kPa following the ice tensile strength measurements 52 . The calibrated stress thresholds are assumed to be constant through time during all simulations because we assume that the ice tensile strength, which may also depend on other properties of ice (e.g., the existence of crevasses, presence of englacial water, ice temperature), would not change over the time scales considered here, and could not be well constrained by observations. We run the simulations with one σ max value for each glacier that gives the best possible matched distance. As shown in the previous study 46 , some glaciers might be sensitive to σ max value, which may add uncertainty to model results. We also did not calibrate σ max over a stable period of ice front position because, as stated in a previous study 18 , the stable glaciers are generally stable due to the bed topography under their terminus and have a wide range of σ max for which the model is also stable 46 . Further investigation is needed to validate the stress threshold value and improve this calving parameterization.
Supplementary Tables 1-6 show the observed retreat distance from 2007 to 2017 and the best possible matched distance from the model. We highlight the retreat distance within 1 km from the observations in blue. The calving law used in this study may not capture all modes of calving, which likely explains the poor fit for glaciers highlighted in red. In the calibration phase, after matching the observed retreat from 2007 to 2017 along a central flow line for each glacier, we compare our model to: the observed velocity, the changes in calving front positions and thinning rates (see Supplementary Data 1), which show a qualitatively good agreement for many glaciers. We also compare the modeled ice mass changes to observed mass changes from 1 to show the calibration results from our model (Fig. 3).
Atmospheric forcings. The ice sheet is forced by the surface mass balance (SMB) of Regional Atmospheric Climate Model (RACMO) 2.3 monthly data for the hindcast simulations 53 . For the forecast simulations, we use the anomalies of SMB from MAR (Modéle Atmosphérique Régional) 54 simulations forced by global climate models (GCMs). Here we use MIROC5, CanESM2 and NorESM1 outputs from CMIP5 based on RCP4.5 and RCP8.5 scenarios and CESM2, CNRM-CM6, CNRM-ESM2, and UKESM1-CM6 outputs from CMIP6 based on SSP585, that are selected for the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6). The RCPs represents 21st century pathways of greenhouse gas (GHG) emissions 34 . We choose one intermediate scenario (RCP4.5) and one with high GHG emissions (RCP8.5), with radiative forcing leading to 4.5 and 8.5 W/m 2 , respectively. The SSPs are new emissions scenarios combined with socio-economic development in the future 55 . We use one extreme scenario (SSP585) which represents fossil-fueled development with 8.5 W/m 2 radiative forcing until 2100. We use a monthly forcing instead of yearly averaged to account for seasonal cycle of change in SMB. Since the MAR simulations are performed for fixed ice surface elevation, we correct the surface mass balance following the gradient method 56 in order to consider changes in ice surface elevation under future warming scenarios. We apply the gradient method only for the ablation regime since this method is not well defined for the accumulation regime 13,56 . Here we do not consider the impact of surface runoff on ice dynamics through damage or hydro-fracture or through subglacial hydrology system because they are not well constrained by observations.
Oceanic forcings. For ocean forcing, we employ parameterizations of the rate of ice shelf melt beneath floating extensions 57,58 and the rate of ice undercutting for marine-terminating glaciers with no floating sections 47 . For the melting under floating ice, we use the parameterization 57 as described: where ρ M is the ocean layer density, C pM is the specific heat capacity of the mixed layer, γ T is the thermal exchange velocity andT is the thermal forcing. We also follow the simple basal melting-water depth parameterization 58 to consider changes in basal melt with ice shelf bottom depth. We assume that basal melting increases linearly with depth between the top-water and the deep-water elevation. We find the different top-water and deep-water elevations for each model domain.
We use the undercutting parameterization 47 to estimate the melt rate ( _ M) at the nearly-vertical ice front of tidewater glaciers: where h is the water depth, q sg is subglacial discharge, A = 3 × 10 −4 m −α day α−1 ∘ C −β , α = 0.39, b = 0.15 ∘ C −β , and β = 1.18. To estimate subglacial discharge for calibration, we integrate the RACMO2.3 runoff field over the drainage basin assuming that surface runoff is the dominant source of subglacial fresh water in summer 47 . The rates of ice shelf melt and ice undercutting define the ocean forcing on the glacier retreat, i.e. an increase in these rates above a steady state value translates into a retreat of the glacier grounding line. For the calibration period, we use a converted thermal forcing on the shelf and into the fjords using bathymetric information and the transfer function developed by the recent study 51 after combining discrete conductivity, temperature and depth (CTD) data and output products from an ocean model operating at 4 km resolution 59 . The recent study 51 provide a detailed description of thermal forcing. For the future simulations, we apply the anomalies of thermal forcing calculated from ocean temperature and salinity of GCM simulations outputs to the average of thermal forcing for the calibration period (Eq. (8)). We use the same GCM outputs as the SMB forcing. We also apply anomalies of runoff from MAR simulations forced by the corresponding GCM models 54 to the average subglacial discharge for the calibration period (Eq. (9)).

Code availability
The ISSM is open source and is available at http://issm.jpl.nasa.gov (Version 4.17, released on April 1st 2020).