2.1. Overview

We applied ensembles of PISM simulations to study the twenty-first century evolution of ice in the Filchner-Ronne region. The workflow is illustrated in Fig. 1. We started with a spun-up, thermodynamically evolved model of the Filchner-Ronne region of Antarctica. We then calibrated a number of parameters controlling ice flow with a Bayesian optimization approach using a comparison to real world observed ice speeds. We applied this model to several scenarios, including using surface and ocean forcing from RCP 2.6 and 8.5, and also to a set of scenarios involving an immediate step change in ocean temperature with constant climate forcing in order to test the response of this region to increased ocean warming.

Fig. 1. Description of modeling workflow. Names of topics in the workflow refer to section names.

2.2. Model description

PISM uses a hybrid stress balance that solves numerical approximations to the shallow ice and shallow shelf equations (Bueler and Brown, Reference Bueler and Brown2009). The effective viscosity of glacier ice, η, is given by:

(1)$$2 \eta = \left(E\, A\right)^{-1} \left(\tau_{\textrm e}^2 + \epsilon^2\right)^{{1-n\over 2n}},\; $$

where $\tau _{\textrm e}$ is the effective stress, E is the enhancement factor, A is the enthalpy-dependent rate factor, and n is the exponent of the power law. The small constant $\epsilon$ (units of stress) regularizes the flow law at low effective stress, avoiding the problem of infinite viscosity at zero deviatoric stress. The pseudo-plastic power law (Schoof, Reference Schoof2010) relates bed-parallel shear stress, τ_b, and the sliding velocity, u_b:

(2)$${\boldsymbol \tau}_{{\rm b}} = - \tau_{{\rm c}}\, {{\bf u}_{{\rm b}}\over \vert {\bf u}_{{\rm b}}\vert ^{( 1-q) } u_0^q},\; $$

where τ _c is the yield stress, q is the pseudo-plasticity exponent, and u ₀ = 100 m yr⁻¹ is a threshold speed. We assume that yield stress τ _c is proportional to effective pressure N (‘Mohr-Coulomb criterion’: Cuffey and Paterson, Reference Cuffey and Paterson2010),

(3)$$\tau_{{\rm c}} = ( \tan\phi) \, N,\; $$

where ϕ is the till friction angle given as a function of bed elevation z:

(4)$$\phi = \left\{\matrix{ \phi_{{\rm min}},\; \hfill & z \le b_{{\rm min}},\; \hfill \cr \phi_{{\rm min}} + ( z - b_{{\rm min}}) \, {\phi_{{\rm max}}-\phi_{{\rm min}}\over b_{{\rm range}}},\; \hfill & b_{{\rm min}} < z < b_{{\rm max}},\; \hfill \cr \phi_{{\rm max}},\; \hfill & b_{{\rm max}} \le z. \hfill }\right.$$

The parameter b _range (m) is defined as b _range = b _max − b _min and is used as a tuneable parameter in this study. The effective pressure N is given by Tulaczyk and others (Reference Tulaczyk, Kamb and Engelhardt2000) and Bueler and van Pelt (Reference Bueler and van Pelt2015):

(5)$$N = \delta P_o \, 10^{( e_0/C_{{\rm c}}) \left(1 - ( W/W^{{\rm max}}) \right)}.$$

Here δ is a lower limit of the effective pressure, expressed as a fraction of overburden pressure, e ₀ is the void ratio at a reference effective pressure N ₀, C _c is the coefficient of compressibility of the sediment, W is the effective thickness of water in the till, and W ^max is the maximum amount of basal water. We use a non-conserving hydrology model that connects W to the basal melt rate $\dot B_{{\rm b}}$ (kg m⁻² a⁻¹; Tulaczyk and others, Reference Tulaczyk, Kamb and Engelhardt2000):

(6)$${\partial W\over \partial t} = { \dot B_{{\rm b}}\over \rho_{{\rm w}}} - C_{{\rm d}},\; $$

where ρ _w is the density of water and C _d = 1 mm yr⁻¹ is a fixed drainage rate.

PISM has been used in many ice sheet model intercomparison projects (e.g. Seroussi and others, Reference Seroussi2019; Cornford and others, Reference Cornford2020; Levermann and others, Reference Levermann2020; Seroussi and others, Reference Seroussi2020; Sun and others, Reference Sun2020), and PISM modeling studies have helped projections of the trajectories of the Greenland and Antarctic ice sheets (e.g. Aschwanden and others, Reference Aschwanden2019; Golledge and others, Reference Golledge2019).

2.3. Model setup and initialization

The Filchner-Ronne drainage basin is represented in PISM using a 1799 km × 2159 km domain (Fig. 2). The maximum extent of the ice shelf is constrained to the extent in 2016 to prevent further ice shelf advance, although ice shelf calving fronts are allowed to retreat. At the domain boundaries, we forced the gradients of ice thickness and surface elevation to be zero, and set ice thicknesses to zero beyond the boundary. Inside that region ice is free to evolve.

Fig. 2. Map of modeled domain showing observed and modeled ice velocities. (a) Observed 2018 ice velocities from ITS_LIVE, with ice streams and ice shelves labeled. (b) Modeled ice velocities, showing median 2018 velocity for each gridcell. Grounding lines and ice shelf bounds are shown in black, from the 2020 outlines from SCAR, which are close to the grounding line locations of at the start of the model run.

As ice is a good insulator and the thermal regime needs multiple millennia to be in equilibrium with the climate conditions at the surface, we perform a step-wise initialization method. The ice sheet requires time to reach thermal equilibrium as well as dynamic equilibrium, as dynamics are a function of temperature. A simulation of temperature evolution at 10 km grid resolution keeping the geometry fixed was run for approximately 100 000 years, followed by a 10 000 year model run to reach dynamic equilibrium, and then 1000 years at 1 km grid resolution with fixed ice thicknesses. Finally, we updated bed topography and present day thicknesses from Bedmap2 (Fretwell and others, Reference Fretwell2013) to BedMachine (Version 2; Morlighem and others, Reference Morlighem2020) and performed a 500 year simulation to reach a new thermal equilibrium. The switch from Bedmap2 to BedMachine was done in order to utilize existing work at recreating modern the dynamic and thermodynamic conditions of Antarctica within PISM. Directly following the switch from Bedmap2 to BedMachine, ice was frozen to the bed along some ice streams due to energy conservation as some ice thickness increased, and the additional time was required to restart the ice stream activity.

2.4. Model forcing

PISM was forced by the 1950–2014 mean near-surface air temperature and climatic mass-balance fields from the regional climate model RACMO2.3p2 remapped to 1 km resolution (van Wessem and others, Reference van Wessem2018). RACMO was forced with the Community Earth System Model Version 2 (CESM2; Danabasoglu and others, Reference Danabasoglu2020). During the model setup and initialization, ocean temperature and ice shelf basal mass fluxes were generated with the Finite Element/columE Sea ice-Ocean Model (FESOM; Wang and others, Reference Wang2014), and applied using constant melt rates. In the Filchner-Ronne region FESOM is capable of reproducing the subshelf circulation to produce realistic shelf basal melt fields (Timmermann and Goeller, Reference Timmermann and Goeller2017).

For model calibration and future projections, we calculated ice shelf basal mass flux using PICO (Reese and others, Reference Reese, Albrecht, Mengel, Asay-Davis and Winkelmann2018b). PICO is an ocean module parametrization which assumes warm water reaches the grounding line of the ice shelf and then basal melt is modeled using the overturning as that water rises and cools with melt and sub-ice shelf geometry. The inputs for PICO include the mean ocean salinity and potential temperature at the depth of the continental shelf off the ice shelf front across the entire drainage basin and a few other parameters, most importantly an ocean overturning coefficient $\gamma _T^\ast$ and heat exchange coefficient C. Although PICO does not capture the complexities of ocean circulation underneath the ice shelf of a more advanced ocean model (e.g. FESOM), it reproduces typical melt patterns within the range of observed melt rates around Antarctica (Reese and others, Reference Reese, Albrecht, Mengel, Asay-Davis and Winkelmann2018b) and provides a realistic method to model melt using different ocean temperatures (Hill and others, Reference Hill, Rosier, Gudmundsson and Collins2021).

Calving is controlled by an eigen-calving scheme which has been found to be appropriate for Antarctic ice shelves (Levermann and others, Reference Levermann2012). Calving is proportional to the product of eigenvalues of strain rates for diverging flow, with a coefficient K. Overall ice mass has not been found to be sensitive to reasonable values of the parameter K. Albrecht and others (Reference Albrecht, Winkelmann and Levermann2020a) tested values of K between 10¹⁶ and 10¹⁸ m s and found this law caused present-day calving to occur largely in areas that provide little buttressing; calving occurred within the ‘passive ice shelf’ areas described by Fürst and others (Reference Fürst2016). As such we set the value of K here to be 2 × 10¹⁵ m s based on a manual calibration, but we do not tune it further. An additional minimum ice shelf thickness criteria of 100 m is applied, with thinner ice calved off.

2.5. Model calibration

Uncertainties in ice flow parameters can have a large impact on projections of ice mass changes in modeling studies (Albrecht and others, Reference Albrecht, Winkelmann and Levermann2020a; Zeitz and others, Reference Zeitz, Levermann and Winkelmann2020). To mitigate this issue, we calibrated six parameters controlling ice flow to better reproduce observed surface speeds These parameters are: the Shallow Ice and Shallow Shelf enhancement factors (SIAe and SSAe, respectively), exponent of the pseudo-plastic sliding law (q), lower limit of till effective pressure as a fraction of overburden (δ), till friction angle depth minimum (b _min) and till friction angle depth range (b _range). Naming conventions follow PISM documentation (The PISM Authors, 2018), except for δ which is called till effective fraction overburden in PISM, and b _range. The applications of these parameters are described by equations 1–5. The values of model parameters which were not varied come from PISM default values, Albrecht and others (Reference Albrecht, Winkelmann and Levermann2020a), or manual calibration.

The enhancement factors SIAe and SSAe are coefficients applied to the Glen-Paterson-Budd-Lliboutry-Duval flow law (Lliboutry and Duval, Reference Lliboutry and Duval1985), with enhancement factors for shelf/stream plug flow usually being expected to be below 1 (Ma and others, Reference Ma2010). Sliding is controlled by a power law relating velocities and basal shear stress (Eqn. (2)), with exponent q such that q = 0 is coulomb-plastic sliding and q = 1 is purely linear sliding. Sliding occurs when the driving stress is greater than a basal shear stress, which is a function of a number of parameters including the till friction angle. The till friction angle varies linearly with bed depth from a minimum angle ϕ _min at b _min to a maximum angle ϕ _max at b _max (Eqn. (4)). The ranges of parameter values used in this calibration are shown in Table 1. To determine the range of values for the calibrated parameters, the entire calibration process was done twice, the first time with a wider range of values, and the second time with a more narrow range. The parameter values for each ensemble were drawn using Sobol Sequences, which is a quasi-random sampling method intended to evenly search the parameter space.

Table 1. Model parameters used in PISM. Parameters that were tuned during the calibration process have the values listed as ‘varied’ range of values tested listed. The varied PICO parameters were not tuned in the calibration

Our calibration closely follows the methodology of Brinkerhoff and others (Reference Brinkerhoff, Aschwanden and Fahnestock2021) which estimated the joint parameter distributions given observations using Markov-chain Monte Carlo sampling (MCMC). We rendered the problem tractable by constructing an artificial neural network to act as a surrogate for the ice flow model, which provided a mapping from ice flow parameters to surface speeds.

To begin we created a training dataset with a 300 member model ensemble. For each ensemble member we then ran PISM forward for 50 years at 2 km resolution while keeping atmosphere and ocean forcing constant. Ocean forcing was provided by PICO using modern temperatures and salinities (Schmidtko and others, Reference Schmidtko, Heywood, Thompson and Aoki2014), and using the PICO default parameters for ocean overturning and heat exchange coefficients in PISM (2 × 10⁻⁵ m s⁻¹ and 1 × 10⁶ m⁶ s⁻¹ kg⁻¹, respectively). The values of the constant model parameters are given in Table 1. For each of the six calibrated parameters listed above, we drew values using Sobol Sequences from within their respective ranges.

We extracted the ice surface speeds from the training model runs and used principal component analysis to produce a set of ≈50 basis functions, enough to describe 99.95% of the velocity variability. Then we trained an artificial neural network-based surrogate model to map parameter combinations to the principal components and recover the velocity field by multiplying the principal components and the basis functions. Finally, we estimated the joint posterior parameter distributions by using MCMC to draw 10⁵ samples, which were compared to observed surface speeds taken from the ITS_LIVE map of 2018 velocities (Gardner and others, Reference Gardner, Fahnestock and Scambos2019). The MCMC sampling was done in log-space for ice speeds. The observed ice velocities were regridded to the 2 km model grid using bilinear interpolation. Velocities from 2018 were used because they were the most recent available when the comparison was carried out. We used one year of velocity data for this analysis. The annual variation surface speed is very low for locations with speed above 50 m a⁻¹ in the ITS_LIVE record (between 2014–2018 the average standard deviation of surface speed was 0.005 of its magnitude for these locations). Surfaces with slower speeds had higher annual variability, although the interannual variability appears randomly distributed in space so is unlikely to have a strong impact on the parameter calibration which optimizes parameters across the entire domain. From the posterior distributions (shown in the supplementary material) we drew parameter samples for our study scenarios.

We use these parameter sets for each of the climate scenarios in this study, and one of these sets is shown in Fig. 3. A comparison of the observed surface speeds in 2018 and the simulated median speed for each gridcell is given in Fig. 2. The modeled surface velocities are drawn from a scenario with the same climate forcing as the calibration model runs. Overall, speeds are reproduced to within 100 m a⁻¹ near the grounding lines of some of ice streams, including the Bailey, Slessor, and Support Force Glacier (Fig. 2). The most prominent differences between model and observations are higher observed than modeled ice speeds along the ice streams of the southwestern Ronne Ice Shelf, especially including the Evans and Rutford ice streams. Our model has ice draining into the southwestern Ronne Ice Shelf from Carlson Inlet and the Rutford Ice stream, yet presently all ice flows out through the Rutford Ice Stream instead (Fig. 2a). It has been suggested that Carlson Inlet was an active ice stream as recently as several hundred years ago (Vaughan and others, Reference Vaughan, Corr, Smith, Pritchard and Shepherd2008), although other field observations instead suggest Carlson Inlet is not a periodically activating ice stream at all (King, Reference King2011). In our model the Möller Ice Stream has much higher ice velocities and a much greater extent than observations. The goal of the ensemble parameter optimization was to best reproduce ice velocities across the entire region, and as such parameters are not tuned to correctly represent the behavior of individual ice streams.

Fig. 3. Frequency of parameter values in the ocean sensitivity analysis ensembles, shown here as the total occurrences for each range. Parameter values were sampled from the parameter posterior distributions, shown in Fig. S1.

2.6. RCP scenarios

In a first set of experiments, we assessed the region's response to the projected range of 21st century climate forcing. We used anomalies of mean annual ocean potential temperature, near-surface air temperature, and climatic mass balance for 2015–2100. For ocean temperatures, we applied a linearly increasing annual ocean potential temperature anomaly in PICO, with the annual rate of increase obtained from the mean value of ocean potential temperature change in the Weddell Sea region over 1995–2100 at 400–600 m depth from Little and Urban (2016) for RCP 2.6 and 8.5. These anomalies were added to the observed 2014 ocean potential temperatures from Schmidtko and others (Reference Schmidtko, Heywood, Thompson and Aoki2014), with salinity held constant.

Several climate parameters varied between ensemble members, and these were drawn in a different manner. Within PICO, the parameters for overturning coefficient $\gamma _T^\ast$ and heat exchange coefficient C have a strong impact on melt (Reese and others, Reference Reese, Levermann, Albrecht, Seroussi and Winkelmann2020). We drew values for these parameters from within the most plausible range found by Reese and others (Reference Reese, Gudmundsson, Levermann and Winkelmann2018a). Specifically, Fig. 4 of Reese and others (Reference Reese, Gudmundsson, Levermann and Winkelmann2018a) identified the 14 most realistic $\gamma _T^\ast$ and C value combinations, and for each ensemble member we randomly chose one of those combinations and then uniformly drew $\gamma _T^\ast$ and C values from within the range of values closest to that combination. This resulted an equal chance of values of C to be less and greater than 1 × 10⁶ m⁶ s⁻¹ kg⁻¹. We drew values in this manner to match our knowledge of the plausible parameter configurations, and not skew our values of C to be too high as to likely overestimate basal melt.

For the near-surface air temperatures and climatic mass balance, we obtained annual anomaly values from three Earth System Models (ESMs) for RCP 2.6 and 8.5. These three models, as prepared by Barthel and others (Reference Barthel2020) at 2 km resolution, are the Norwegian Climate Center's ESM (NorESM1-M), Community Climate System Model (CCSM4), and the atmospheric chemistry-coupled ESM from the Model for Interdisciplinary Research on Climate (MIROC-ESM-CHEM). The temperature and mass balance anomalies were applied to the mean 1950–2015 climatic mass balance and near-surface temperature fields (van Wessem and others, Reference van Wessem2018). We ran 200 simulations at 2 km resolution for both RCP 2.6 and 8.5, and we randomly assigned one of the three surface anomaly models to each ensemble member, with each of the three models being used between 66 and 68 times.

2.7. Ocean warming scenarios

For a second set of experiments, we tested the sensitivity of ice in this region to thinning of the ice shelf in response to step change increases in ocean temperatures. These are not based on direct climate projections but rather help to explore the characteristics of the response of this region to increased ice shelf melt. We produced a set of scenarios involving raising present-day ocean potential temperatures at depth by 1–6 $^\circ$C at 1 $^\circ$C intervals. These potential temperature changes were applied to present day ocean temperatures (Schmidtko and others, Reference Schmidtko, Heywood, Thompson and Aoki2014) as a step change. For surface temperature and climatic mass balance we used the RACMO 1950–2014 mean values. We will refer to these as the +1 $^\circ$C to +6 $^\circ$C ocean warming scenarios and an additional control scenario was run with constant climate and no ocean warming enforced. For each scenario we ran 100 ensemble members at 2 km resolution from 2015–2100, holding ocean and surface forcing constant, aside from the step change in ocean potential temperature at the beginning. The ocean overturning and heat exchange coefficient parameters in PICO were also held constant.

Only the lowest warming scenario presented here is consistent with projected ocean warming at depth, whereas the ocean temperature scenarios beyond +1 $^\circ$C are well outside the range of 21st century warming under the higher RCP 8.5 scenario (Little and Urban, Reference Little and Urban2016). However, changes in ocean circulation forced by changes in surface winds, potentially possible in the 21st century, can cause an order of magnitude increase in ice shelf melting in the region (Hazel and Stewart, Reference Hazel and Stewart2020; Naughten and others, Reference Naughten2021). The higher end of these warming scenarios are well outside the range found in century-scale projections, but the +6 °C warming scenario provides a forcing comparable to, and perhaps more realistic than, removing all floating ice instantaneously as in the Antarctic Buttressing Model Intercomparison Project (Sun and others, Reference Sun2020).