The role of history and strength of the oceanic forcing in sea-level projections from Antarctica with the Parallel Ice Sheet Model

Mass loss from the Antarctic Ice Sheet constitutes the largest uncertainty in projections of future sea-level rise. Ocean-driven melting underneath the floating ice shelves and subsequent acceleration of the inland ice streams is the major reason for currently observed mass loss from Antarctica and is expected to become more important in the future. Here we show that for projections of future mass loss from the Antarctic Ice Sheet, it is essential (1) to better constrain the sensitivity of sub-shelf melt rates to ocean warming, and (2) to include the historic trajectory of the ice sheet. In particular, we find that 5 while the ice-sheet response in simulations using the Parallel Ice Sheet Model is comparable to the median response of models in three Antarctic Ice Sheet Intercomparison projects – initMIP, LARMIP-2 and ISMIP6 – conducted with a range of ice-sheet models, the projected 21st century sea-level contribution differs significantly depending on these two factors. For the highest emission scenario RCP8.5, this leads to projected ice loss ranging from 1.4 to 4.0 cm of sea-level equivalent in the ISMIP6 simulations where the sub-shelf melt sensitivity is comparably low, opposed to a likely range of 9.2 to 35.9 cm using the 10 exact same initial setup, but emulated from the LARMIP-2 experiments with a higher melt sensitivity based on oceanographic studies. Furthermore, using two initial states, one with and one without a previous historic simulation from 1850 to 2014, we show that while differences between the ice-sheet configurations in 2015 are marginal, the historic simulation increases the susceptibility of the ice sheet to ocean warming, thereby increasing mass loss from 2015 to 2100 by about 50%. Our results emphasize that the uncertainty that arises from the forcing is of the same order of magnitude as the ice-dynamic response for 15 future sea-level projections. 1 https://doi.org/10.5194/tc-2019-330 Preprint. Discussion started: 21 January 2020 c © Author(s) 2020. CC BY 4.0 License.

The sub-surface ocean forcing informs parameterizations that provide melt rates underneath the ice shelves for ice-sheet mod-55 els. For the ISMIP6 experiments, a depth-dependent, non-local parameterization and a depth-dependent, local parameterization have been proposed (Jourdain et al., under review) that both mimic a quadratic dependency of melt rates on thermal forcing (Holland et al., 2008). As an alternative, more complex modules that capture the basic physical processes within ice-shelf cavities have been developed recently (Lazeroms et al., 2018;Reese et al., 2018a). We here analyse results as submitted to ISMIP6 that apply the Potsdam Ice-shelf Cavity mOdel (PICO; Reese et al., 2018a) which extends the ocean box model (Olbers and 60 Hellmer, 2010) for application in three-dimensional ice-sheet models. The model has been tested and compared to other parameterizations for an idealized geometry (Favier et al., 2019). In this case, the induced ice-sheet response matches the response driven by a three-dimensional ocean model. In contrast to ISMIP6, the LARMIP-2 experiments are forced by basal-melt rate changes directly. Scaling factors between global mean temperature changes and Antarctic sub-surface temperature changes are determined from CMIP5 models. These are used to generate ocean temperature forcing under different RCP scenarios 65 emulated from MAGICC-6.0 RCP realisations (Meinshausen et al., 2011). Sub-shelf melt rates are assumed to increase by 7 to 16 m a −1 per degree of sub-surface ocean warming, based on Jenkins (1991) and Payne et al. (2007).
Here we compare simulations with the Parallel Ice Sheet Model as submitted to ISMIP6 with results obtained following the LARMIP-2 protocol and analyse (1) the effect of the oceanic forcing and (2) the effect of a historic simulation preceding the projections. In Sect. 2 we describe the methods used and the initial configurations of PISM. This is followed by an analysis of 70 the experiments for ISMIP6 with only ocean forcing applied and the results obtained when following the LARMIP-2 protocol in Sect. 3. These are compared and discussed in Sect. 4 and 5.

Methods
We use the comprehensive, thermo-mechanically coupled Parallel Ice Sheet Model (PISM, Bueler and Brown, 2009;Winkelmann et al., 2011;the PISM authors, 2019) which employs a superposition of the Shallow-Ice and Shallow-Shelf Approxima-75 tions (Hutter, 1983;Morland, 1987;MacAyeal, 1989). We apply a power-law relationship between SSA basal sliding velocities and basal shear stress with a Mohr-Coulomb criterion relating the yield stress to parameterized till material properties and the effective pressure of the overlaying ice on the saturated till (Bueler and Pelt, 2015). Basal friction and sub-shelf melting are linearly interpolated on a sub-grid scale around the grounding line (Feldmann et al., 2014). In order to improve the approximation of driving stress across the grounding line, the surface gradient is calculated using centered differences of the ice thickness geometry to obtain a thermodynamic equilibrium with present-day climate on 16 km. Based on this, an ensemble of simulations with varying model parameters is run for several thousand years towards dynamic equilibrium on 8 km horizontal resolution.
The simulations employ 121 vertical layers with a quadratic spacing from 13 m at the ice shelf base to 100 m towards the surface. We vary parameters of PICO (heat exchange coefficient γ T and overturning coefficient C) as well as the minimum till 90 friction angle in the parameterized till material properties (Φ min ). The initial configuration is selected in two steps: after 5000 years of model simulation, 5 candidates that compare best to present-day observations of ice geometry and speed (Fretwell et al., 2013;Rignot et al., 2011) are selected and continued. After 12, 000 years the best fit equilibrium result was selected among them and used as initial configuration for the projections, see Fig. S.1. We assess the ensemble members at each step using a scoring method Albrecht et al., 2020) that tests for root-mean-square deviation to present-day ice 95 thickness, ice-stream velocities, as well as deviations in grounded and floating area, and the average distance to the observed grounding line position. We lay a specific focus on the Amundsen region, Filchner-Ronne and Ross ice shelves by additionally evaluating each indicator for these drainage basins individually.
The historic simulation is based on the same initial steady-state configuration and additionally applies atmospheric and oceanic forcing over the period from 1850 to 2014 as described below. The initial state for the experiments without historic simulation, 100 hereafter referred to as INIT*, and the initial configuration after the historic simulation, hereafter referred to as INIT, are shown in Fig S.2. The INIT* configuration is very close to a steady-state with ice volume change rates being 5 mm over 85 years while the INIT state is out of balance with ice volume change rates being −1.5 cm over 85 years, see Table 1

Experiments
We here present experiments based on the ISMIP6, LARMIP-2 and initMIP protocols that were done for both initial configurations. A list of all experiments is given in Table S.1. The initMIP experiments employ idealized forcing designed to test the model response to simplified forcing of the surface mass balance (experiment 'asmb'), and the basal mass balance (experiment 'abmb') which increase linearly for 50 years and are kept constant afterwards .

115
For LARMIP-2, constant step-forcing perturbations of the basal mass balance (4, 8 and 16 m a −1 ) are applied in five Antarctic regions (Antarctic Peninsula, East Antarctica, Ross Sea, Amundsen Sea, Weddell Sea). From the modeled sea-level response, linear response functions are derived that can be used to emulate the model's response to arbitrary melt forcing.
We run experiments for both initial configurations with * indicating simulations starting from the pseudo-steady state in 2015, INIT*. The control experiments for both initial configurations employ constant climate conditions as described in the following two subsections.

Atmosphere forcing
Surface mass balance and ice surface temperatures for the initial configuration without historic forcing are from RACMOv2.3p2 (1986to 2005averages, Van Wessem et al., 2018, remapped from 27 km resolution. The historic simulation is started from the same conditions with historic surface mass balance and surface temperature changes following the NorESM1-M simulation as suggested by ISMIP6 (Bentsen et al., 2013). The historic forcing from NorESM1-M is normalized to its initial period (1950- In contrast to ISMIP6, where surface mass balance and surface temperature changes are driven by GCM data, we here keep surface conditions -in line with LARMIP-2 -constant throughout the projections. Note that due to changes in the ice-sheet 135 extent, surface mass balance integrated over the entire ice sheet might change slightly, see Table 1. Surface mass balance and temperatures in the projections that start from the pseudo steady-state INIT* are given by the RACMO climatology. For the projections based on the historic simulation we created a new climatology to account for increases in surface mass balance and temperatures in the historic simulation. We avoid using exceptionally high or low values that arise from interanual variability at a specific snapshot in time by using the 1995 to 2014 average of the respective fields.

Ocean forcing
Sub-shelf melt rates are calculated by PICO which extends the ocean box model by Olbers and Hellmer (2010) for application in 3-dimensional ice sheet models (Reese et al., 2018a). It mimics the vertical overturning circulation in ice-shelf cavities and has two model parameters that apply for all Antarctic ice shelves simultaneously: C related to the strength of the overturning circulation and γ T related to the vertical heat exchange across the ice-ocean boundary layer. We here use parameters C = 1 × 10 6 m 6 s −1 kg −1 and γ T = 3 × 10 −5 m s −1 that were found to yield realistic melt rates in comparison to present-day estimates (Reese et al., 2018a;Rignot et al., 2013). The value of γ T is slightly higher than the reference value as an outcome of the ensemble study, see We initialize PICO with an ocean data compilation from the World Ocean Atlas 2018 pre-release Zweng et al., 2018) and Schmidtko et al. (2014). PICO is driven by ocean temperature and salinity averaged over the depth 150 of the continental shelf within each drainage basin. The data from the WOA2018 pre-release is processed by determining the relevant depth from bathymetric access to ice-shelf cavities. In Dronning Maud Land (PICO basins 2 to 5), where ocean temperatures have a warm bias due to the lack of data along narrow continental shelves, values from Schmidtko et al. (2014) were used. Using the currently observed 'warm' conditions in the Amundsen Sea, we found that region to collapse in the initial ensemble irrespective of basal sliding parameters. As this region is out of balance today due to oceanic forcing (e.g., Konrad 155 et al., 2018; Shepherd et al., 2018), it would be inconsistent to initialize our model by running it towards equilibrium over several thousand years applying constant present-day climate forcing. We hence reduced temperatures in the Amundsen Sea to 'cold' conditions (−1.25 • C, Jenkins et al., 2018).
Ocean temperature and salinity forcing is calculated from CMIP5 models using an anomaly approach as suggested for ISMIP6 (Barthel et al., 2019;Jourdain et al., under review). We average these values over 400 to 800 m depth to obtain suitable input 160 for PICO. The historic forcing is based on NorESM1-M (as suggested for ISMIP6) and anomalies are normalized to the initial period (here 1850-1900), similar to the atmosphere forcing. A new ocean climatology for the experiments starting from the historic simulation is obtained from the 1995 to 2014 average conditions.
For LARMIP-2, we add melt rate anomalies to the underlying PICO melt rates in different Antarctic regions as described in Levermann et al. (2020). Using linear response theory, the probability distribution of the sea-level contribution for RCP8.5 is 165 then calculated following the LARMIP-2 protocol.

Results
We present here (1) the results for the two initial configurations submitted to ISMIP6 and (2) the sea-level estimates for RCP8.5 obtained following the LARMIP-2 and ISMIP6 experiments based on the historic configuration.

Initial configurations and historic simulation 170
The two initial configurations for 2015, one based on a pseudo-equilibrium and one on a historic simulation from 1850 to 2014, do not differ much in terms of state variables such as ice thickness, volume or speed (see Sect. 2.1). However, the configurations have opposed change rates: INIT* has a small tendency to gain mass and INIT is clearly out of balance and loses mass (compare the control simulations in Table 1). Over the historic period, the ice sheet thins along its margins through increased sub-shelf melting and at the same time thickens in the interior due to more snowfall. These signals are smaller than 50 m over grounded

Comparison to initMIP Antarctica
Results from the idealized surface mass balance experiment 'asmb' as described in initMIP Antarctica  are very similar for both initial states with 119 mm SLE of mass gains for the 'historic' configuration INIT and 118 mm SLE for the 'cold-start' configuration INIT* after 85 years of simulation with respect to the control simulations, see Table 1. This is close to the response of the different models that participated in initMIP Antarctica which showed mass gains between 125 190 and 186 mm SLE after 100 years.
For the idealized basal melt rate experiment 'abmb' from initMIP Antarctica, both states are also quite similar with mass loss of 43 and 40 mm SLE after 85 years, respectively, see  (Bamber et al., 2019). Furthermore, we find that the historic simulation makes  the configuration more susceptible to ocean forcing, see Fig. 2. Ocean-driven mass loss in comparison to the control run is increased by about 50% (factor 1.5) when starting from the historic simulation in contrast to the 'cold-start'.

LARMIP-2 basal melt rate forcing experiments 205
In LARMIP-2, sea-level probability distributions from the Antarctic Ice Sheet are derived using linear response functions as described in Levermann et al. (2020). The response functions are derived from experiments in which constant basal melt rate forcing is applied for five different regions of Antarctica. We here perform the same experiments for both configurations described in Sect. 2.
We find that for all regions the ice sheet response compares with the responses found in LARMIP-2 as, for example, in the  Following the procedure in LARMIP-2, we derive response functions from the idealized experiments for the five Antarctic regions. We then convolve the response function with basal melt rate forcing, given in Fig. 4, to obtain a probability distribution of the future sea-level contribution for RCP8.5 which is given in Fig. 5. The ocean-driven mass loss from 2015 to 2100 has a very likely range of 3.5 to 55.6 cm SLE, a likely range of 9.1 to 35.8 cm SLE and a median of 18.3 cm SLE (5 to 95%, 16.6 to 83.3%, and 50% percentiles, respectively, see Table 2). Similar to the ISMIP6 simulations, these results obtained for 220 the historic initial configuration are larger than the results for the steady-state configuration, with increases between 5 and 7%.
In comparison, the PISM-PIK contribution of LARMIP-2 has a very likely range of 7 to 48 cm SLE, a likely range of 11 to 31 cm SLE and a median of 19 cm SLE for the 21st century. The resulting sea-level probability distribution is hence in line with the estimates presented in LARMIP-2.

225
In the following, we compare the results found in the ISMIP6 and LARMIP-2 experiments, discuss the role of the ocean forcing and of the historic simulation.

Comparison of LARMIP-2 and ISMIP6 sea-level projections
The projected mass loss in LARMIP-2 is an order of magnitude larger than the ocean-driven mass loss in our ISMIP6 experiments for RCP8.5, see Sect. 3. In order to understand this difference better, we here investigate the ocean forcing in more   vulnerability of the Totten region then causes higher overall mass loss.
In Figure 6 we assess for each region the mass loss by applying the response functions to the corresponding PICO melt rate changes driven by NorESM1-M ocean forcing, once averaged over the entire ice shelves and once close to the grounding lines. When comparing the respective mass loss with the ISMIP6 simulation, we find that indeed the changes at the grounding line provide an upper limit while the changes over the entire ice shelf provide a lower limit for the actual mass loss.
Overall, we find that mass losses in the ISMIP6 projections are generally lower than the likely range in LARMIP-2, and in the Weddell Sea losses are smaller than the very likely range, as the basal melt rate changes in the LARMIP experiments are an order of magnitude higher than those estimated with PICO and ISMIP6 forcing.

Role of ocean forcing and basal melt rate sensitivity
In order to gain a better understanding of the conversion of ocean forcing to basal melt rates in LARMIP-2 and in our ISMIP6 contribution, we further analyse the sensitivity to ocean warming.
We perform step-forcing experiments for both initial configurations and diagnose the effect on basal melt rates, see Fig. 7.
Ocean temperatures are increased by 0.5, 1, 2, 3 and 4 • C and the corresponding basal melt rates for constant ice-shelf geome-270 tries are diagnosed using PICO. We find that the sensitivity in the Amundsen Sea Region is comparatively high with about 10 ma −1 K −1 , while the sensitivity in the Weddell Sea is lower with about 1.5 ma −1 K −1 , which yields for the entire Antarctic ice shelves an overall sensitivity of about 2.2 ma −1 K −1 . The sensitivities for melting close to the grounding line are as expected a bit higher: 10.5 ma −1 K −1 for the Amundsen Sea region, 3.9 ma −1 K −1 for the Weddell Sea and 5.3 ma −1 K −1 on average for all Antarctic ice shelves. In both cases, the Antarctic-wide sensitivity is substantially lower than the sensitivity used 275 in LARMIP-2. In the latter study, a sensitivity between 7 and 16 ma −1 K −1 , based on Jenkins (1991)  is assumed to translate ocean forcing into sub-shelf melt rates. This is consistent with our findings in the previous section that in the ISMIP6 simulations mass loss from the Antarctic Ice Sheet, and especially from the regions that drain the large Filchner-Ronne and Ross ice shelves, is smaller than the likely range estimated following the LARMIP-2 protocol. Jourdain et al. (under review) report that a different tuning of the ISMIP6 basal melt parameterization to fit observations in the Amundsen Sea (from 280 Dutrieux et al., 2014;Jenkins et al., 2018) substantially increases the sensitivity to ocean changes and Seroussi et al. (under review) find that this enhances the sea-level contribution by a factor of six.
Since the sensitivity in PICO depends on the parameters used, with the overturning coefficient C affecting the sensitivity in large ice shelves and the heat exchange coefficient γ T affecting the sensitivity in small ice shelves, a different tuning could improve basal melt rate sensitivities and thereby lead to higher mass losses in the ISMIP6-experiments. A consistency of sub-shelf 285 melt rates with present-day observations could be achieved by introducing additional degrees of freedom through temperature corrections that reflect uncertainties in ocean properties, as for example used in Lazeroms et al. (2018) and Jourdain et al. (under review). In addition, tuning to realistic melt rates close to the grounding lines (in PICO's first box) is potentially more important than fitting shelf-wide melt rates (e.g., Goldberg et al., 2019;Reese et al., 2018b).
Few observations exist for targeted tuning of the sensitivity of basal melt rate parameterizations to ocean temperatures, hence 290 the use of dynamic modelling of the ocean circulation in ice-shelf cavities could be explored. We estimate that the sensitivity in Seroussi et al. (2017) varies between 6 and 16 ma −1 K −1 with an average of 9.4 ma −1 K −1 over the first twenty years of model simulation, which would be in line with the sensitivities used in LARMIP-2, see Fig. S.5.
Note that we provide linear estimates of the sensitivity of PICO in the discussion above, while Holland et al. (2008) report a quadratic dependency of melt rates on thermal forcing. They also discuss that the sensitivity depends on ice shelf cavity 295 properties such as the slope of the ice-shelf draft and shape of the ice shelf and that sensitivities are generally higher close to the grounding line. Further factors that influence ocean circulation, such as bathymetry, also affect the ocean sensitivity. While PICO takes into account, that not all heat content of the ocean water masses that enter the cavity might be used for melting, it does not capture three-dimensional circulations in ice-shelf cavities that play a role especially for larger ice shelves such as Filchner-Ronne.  could be linked to a local instability that is kicked-off for the simulations starting from the historic configuration but not for those starting from the pseudo steady state. This is less pronounced for CCSM4 and MIROC, maybe due to differences in the ocean forcing and regions contributing to sea level rise. In the idealized experiments for LARMIP-2 (Fig. 3), differences in simulations starting from the two initial states arise especially in East Antarctica, the Weddell Sea and the Amundsen Sea, less in the Ross Sea. The effect of the historic simulation is reduced for the stronger basal melt rate forcing applied in the 315 LARMIP-2 experiments, with mass loss increases in the projections between 5 and 7%.
Furthermore, the ice sheet's response might have changed after the historic simulation due to changes in boundary conditions. Moreover, changes in the ice-sheet state could result since, for instance, the underlying equation system depends non-linearly on the three-dimensional temperature field. The grounding line retreats during the historic simulation slightly into deeper regions, where the local freezing point at the ice shelf base near the grounding line is decreased due to its pressure dependence. Hence 320 more heat is available for melting the ice-shelf base, which shows also in an increased sensitivity to ocean changes, see Fig. 7.
In particular for lower temperatures, PICO shows a non-linear sensitivity of melt rates to ocean temperatures, as discussed in Reese et al. (2018a). Further investigations would be required to disentangle the reasons for the increased susceptibility to ocean warming after the historic simulation, also considering the strength of the forcing applied.
The sea-level contribution over the historic period from 1850 to 2014 is 3.6 mm in comparison to the control simulation. This 325 is smaller than the reported mass loss of the Antarctic Ice Sheet that amounts to 7.6 ± 3.0 mm SLE between 1992 and 2017 . An improved understanding of the basal melt rate sensitivity, potential biases in the atmospheric or oceanic forcing, as well as an extension of the scoring with observed patterns of thickness changes would allow for performing 'hindcasting' experiments that, in a next step, could inform future projections.

330
In this study we compare sea-level projections for RCP8.5 from the Antarctic Ice Sheet as submitted to ISMIP6, using the PICO basal melt rate parameterization and constant surface mass balance forcing, and projections derived following the LARMIP-2 protocol that scales global temperatures to subsurface temperatures and melt rates, both using the Parallel Ice Sheet Model.
Overall, we find that the sea-level contribution driven by ocean forcing in ISMIP6 is smaller than the likely range of the sealevel probability distribution in LARMIP-2. This difference can be explained by the comparably low sensitivity of melt rates to 335 ocean temperature changes for the parameter tuning in PICO in comparison to LARMIP-2 where a sensitivity of 9 to 16m/a/K is used that we found to be consistent with a coupled simulation of Thwaites glacier (Seroussi et al., 2017). Future sealevel projections should hence carefully consider the sensitivity of basal melt rates to ocean changes. Additional observations of ocean conditions and ocean-induced melt rates in combination with ocean modelling are needed to better constrain this sensitivity for the diverse ice-shelf cavities in Antarctica. Furthermore, we find that while the initial state resulting from a 340 historic simulation from 1850 to 2014 is virtually indistinguishable from a steady-state simulation, the historic simulation increases the projected mass loss in 2100 by up to 50%. This means that not only the currently committed sea-level contribution should be considered in projections, but also the effect of the historic forcing on the ice sheet's susceptibility to ocean changes.
'Hindcasting' experiments, that reproduce observed thinning rates and ice loss over the past decades, would be valuable to better constrain model parameters and improve confidence in projections. Hence, further investigations are needed to assess through its Working Group on Coupled Modelling, coordinated and promoted CMIP5. We thank the climate modeling groups for producing and making available their model output, the Earth System Grid Federation (ESGF) for archiving the CMIP data and providing access, the University at Buffalo for ISMIP6 data distribution and upload, and the multiple funding agencies who support CMIP5 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 X. Development of PISM is supported by NASA grant NNX17AG65G and