Interplay of grounding-line dynamics and sub-shelf melting during retreat of the Bjørnøyrenna Ice Stream

The Barents Sea Ice Sheet was a marine-based ice sheet, i.e., it rested on the Barents Sea floor during the Last Glacial Maximum (21 ky BP). The Bjørnøyrenna Ice Stream was the largest ice stream draining the Barents Sea Ice Sheet and is regarded as an analogue for contemporary ice streams in West Antarctica. Here, the retreat of the Bjørnøyrenna Ice Stream is simulated by means of two numerical ice sheet models and results assessed against geological data. We investigate the sensitivity of the ice stream to changes in ocean temperature and the impact of grounding-line physics on ice stream retreat. Our results suggest that the role played by sub-shelf melting depends on how the grounding-line physics is represented in the models. When an analytic constraint on the ice flux across the grounding line is applied, the retreat of Bjørnøyrenna Ice Stream is primarily driven by internal ice dynamics rather than by oceanic forcing. This suggests that implementations of grounding-line physics need to be carefully assessed when evaluating and predicting the response of contemporary marine-based ice sheets and individual ice streams to ongoing and future ocean warming.

• Description of the method employed to force the ocean basal melting parametrization from Pollard & DeConto (2012) in the ice sheet simulations; • Description of the statistical approach used in this study to set the parameters for GRISLI and PSU models.

Air temperature and precipitation forcing
TraCE-21ka transient climate simulation Liu et al. (2009) is used to construct different indexes for air temperature and precipitation representative of Fennoscandia, Svalbard/Barents Sea and Siberia/Kara Sea macro-regions. For each macro-region, seven different nodes in the TraCE-21ka grid are selected . In each grid node, annual mean air temperature T a and annual mean precipitation P a are computed every 100 years from the Last Glacial Maximum (around 21,000 years ago, LGM) to PI (1850 a.d., PI), as shown in Fig.1a-f. In each macro-region, the average annual mean near-surface air temperature T a c and precipitation P a c are considered for each macro-region and at each time step ( Fig.1a-f). For each macro-region and at each time step i, the temperature and precipitation indexes I T and I P are computed as follows, so that the all the indexes take value 1 at LGM and 0 at PI.

Surface mass balance
In both GRISLI and PSU ISMs the surface mass balance takes into account accumulation and ablation. The snow accumulation is computed from the annual mean total precipitation following Marsiat (1994). A linear transition between solid and liquid precipitation depending on the annual mean near-surface air temi ii perature is considered to compute the solid precipitation fraction P sf , The accumulation is then computed from the annual mean total precipitation as follows, The ablation is computed using the Positive-Degree-Days (PDD) semi-empirical method Reeh (1989). The PDD method relates the number of days with positive near-surface air temperature with snow and ice melting rates through observationbased melting coefficients. First of all, the near-surface air temperature is assumed to vary sinusoidally over time. Such an assumption allows to reconstruct the annual near-surface air temperature cycle from the annual and July mean nearsurface air temperature, where A = 1 year. Following Reeh (1989), the number of PDD is then computed from the normal probability distribution around the monthly mean temperatures during the year, i.e., where t is the time variable, T ( o C) is the near-surface air temperature variable, and σ is the standard deviation of the near-surface air temperature. In this study, the altitude-latitude dependent parametrization for the near-surface air temperature standard deviation proposed by Fausto et al. (2011) is used. Based on near-surface air temperature observations at the automatic weather stations on the Greenland Ice Sheet, annual and July mean standard deviations (σ a and σ j respectively) of the near-surface air temperature are parametrized as follows, where φ is the latitude expressed in o N. In such a way, in each grid point of the Eurasian Ice Sheet Complex domain mean annual and July standard deviation of near-surface air temperature are computed depending on the surface elevation and latitude. Then, similarly as for the near-surface air temperature (see Equation 5), a sinusoidally variation of the near-surface air temperature standard deviation over the year is assumed, iii Once the number of PDD is calculated, the amount of snow and ice melted is related to the number of PDD by means of snow and ice melt coefficients (C s and C i respectively) depending on the annual mean near-surface air temperature Tarasov & Peltier (2002), Once values for snow and ice melting coefficients are assigned, ablation takes place in three steps: • fresh snow, if present, is melted first at rate C s . Meltwater resulting from the melted snow is able to percolate into the snow layer and refreeze as superimposed ice. If the amount of superimposed ice exceeds 60% of the annual snowfall, runoff occurs; • when all the snow is melted, superimposed ice, if present, is melted at rate C i ; • when all the snow and the superimposed ice are melted, glacier ice is melted at rate C i .
Depending on the snow accumulation and on the available melting energy in each time step, not all the steps (or even none of the steps) will be carried out, Colleoni et al. (2014).

Ocean basal melting forcing and parameter tuning
In order to force the ocean basal melting parametrization from Pollard & De-Conto (2012), we derived from TraCE-21ka transient climate simulation Liu et al. (2009) four different vertical profiles of ocean temperature and salinity, respectively, which vary every 100 years from the LGM to PI and are representative of four macro-regions (Norwegian Sea, South-Western and North Western Barents Sea, Arctic Ocean, Fig.2, 8). For each macro-region, the ocean temperature and salinity profiles are included both in GRISLI and PSU numerical models in the following way: • first, the LGM and PI ocean temperature and salinity at five selected vertical levels (0 m, 200 m, 400 m, 600 m, 800 m) are considered. The vertical levels are selected depending on typical Barents Sea depths. At each vertical level, the LGM and PI ocean temperature and salinity are obtained from TraCE-21ka simulation and included in GRISLI and PSU. Given a vertical level j , iv LGM and PI ocean temperature and salinity are denoted with T j LGM , T j PI , S j LGM and S j PI respectively. At a given depth z 0 , the ocean temperature and salinity profiles at LGM and PI are computed by linearly interpolating between the two vertical levels including z 0 .
• at each vertical level j , an index for ocean temperature (i j T ) and salinity (i j S ), respectively, is computed every 100 years from TraCE-21ka data. At a given time t the temperature and salinity index, respectively, is given by In such a way all the indexes take value 1 at LGM and 0 at PI. All the ocean temperature and salinity indexes are included in GRISLI and PSU.
During a GRISLI or PSU simulation, given an ice shelf grid point the ISM identifies first the macro-region in which the point is contained depending on its longitude/latitude. The depth of the ice shelf base z b is then computed, and the vertical levels j and j + 1 such that j ≤ z b < j + 1 are identified, where j = 0, 200, 400, 600, 800. At a given time t, the ocean temperature and salinity at vertical levels j and j + 1 are computed as follows, By linearly interpolating between T j (t), S j (t) and T j+1 (t), S j+1 (t) the ocean temperature and salinity at depth z b are computed and used to force the ocean basal melting parametrization from Pollard & DeConto (2012).
Once an adequate oceanic forcing is set, a new strategy to tune the model parameter F m is needed. In fact, the model parameter F m has been hand-tuned in Martin et al. (2011) in order to match simulated and observed grounding-line position in the West Antarctic Ice Sheet. Therefore, the oceanic conditions used to force the ocean basal melting parametrization are drastically different from those derived from TraCE-21ka climate simulation and used in this study ( Fig.2  and 8). An a-priori analysis of the oceanic forcing used in this study is performed, in order to identify a range of values for F m . First of all, a typical grounding-line depth of 400 m in the north-western Barents Sea is considered. Looking at TraCE-21ka simulation, at such depth the warmer water temperature (around 7.5 o C, see Fig.2) is reached around 10 ka BP in the south-western and north-western Barents Sea. A water salinity of around 35.5 psu is reached at the same depth and time frame in the south-western and north-western Barents Sea. By forcing the basal melting parametrization from Pollard & DeConto (2012) with such "warm" ocean v conditions we obtain the maximum ocean basal melting rate between LGM and PI (b max ) as a function of the model parameter F m only, i.e., b max = 41.97 · F m · (9.73) 2 .
In such a way it is possible to identify a range of values for the model parameter F m so that b max matches the ocean basal melting values observed under the present-day Antarctica ice shelves. Estimates for the annual basal mass-loss rates of Antarctic ice shelves from Depoorter et al. (2013) and Rignot et al. (2013) range between 0.1 − 22 m/yr. However, the averaged annual basal mass-loss rate of ice shelves in the West Antarctica Peninsula are lower 0.8−1.0 m/yr, Depoorter et al. (2013). Therefore, we decided to set a range for b max between 0.1 − 6 m/yr. In such a way, the range identified for the model parameter F m is 0.02 · 10 −3 ≤F m ≤ 1.56 · 10 −3 .

Latin Hypercube Sampling of model parameters
This statistical approach was previously used by Stokes & Tarasov (2010) to simulate the Laurentide Ice Sheet and Applegate et al. (2015) to simulate the Greenland Ice Sheet. In this study five GRISLI model parameters are considered for the Latin Hypercube Sampling (see Table 3): • the topographic lapse-rate λ, which represents an approximation of how much the near-surface air temperature changes with elevation. This parameter is poorly constrained in large-scale ice sheet modeling studies, whereas topographic lapse-rate values computed from climate simulations range from 4 to 7 o C/km Abe-Ouchi et al. (2007). Following Stone et al. (2010) and Colleoni et al. (2016), in this study the range 4 to 8.2 o C/km is explored; • the precipitation-correction factor γ , which approximates the saturation pressure of water vapour Charbit et al. (2002). In large-scale ice sheet modeling studies this parameter ranges between 0.03 and 0.078 o C −1 (e.g., Charbit et al. (2002)), although climate modelling studies suggest that γ assumes higher values up to 0.11 o C −1 . Following Colleoni et al. (2016), in this study the range 0.03 − 0.1 o C −1 is explored; • the SIA-enhancement factor E SIA , accounting for the anisotropy of polycristalline ice under condition of simple-shear flow, Ma et al. (2010). In large-scale ice sheet modeling studies this parameter is set to values ranging from 1 to 5  and references therein). However, an higher value of 5.6 is suggested by Ma et al. (2010) in a study where an anistropic full-Stokes model is used, thus explicitly accounting for grain orientation (fabric). Therefore, in accord with Colleoni et al. (2016) the range 1 − 5.6 is explored; • the basal drag coefficient c f , which regulates the resistive force acting at the base of the ice sheet in ice stream regions treated with the shallow-shelf approximation. In Wekerle et al. (2016) the strong impact of this parameter on the GRISLI-simulated ice thickness is shown. This parameter is set in early works with GRISLI to 1 · 10 −5 in Peyaud et al. (2007), to 9 · 10 −5 in Dumas (2002) and between 10·10 −5 −100·10 −5 in Álvarez Solás et al. (2011). In this study the range 1 · 10 −5 -10 · 10 −5 is explored; • the ocean basal melting parameter F m the ocean basal melting parametrization from Pollard & DeConto (2012). In the previous section we show the approach adopted to set a range for this parameter. The range explored in this study is 0.02 · 10 −3 -1.56 · 10 −3 .
Following this statistical approach, a group of 101 simulations is performed with GRISLI. The Eurasian Ice Sheet Complex ice volume evolution during the deglaciation simulated with GRISLI is compared with those simulated with the global glacio-isostasy model ICE- 5G Peltier (2004), and a best fit simulation (referred to as GBvol) out of the statistical ensemble is identified. We did not use ICE-6G to constrain our ice sheet simulations since some flaws have been recently found in the more recent ICE- 6G Peltier et al. (2015) reconstruction Purcell et al. (2016). Due to large computational costs, a similar statistical approach was not feasible with PSU. However, GRISLI and PSU have in common four of the five model parameters that are explored with GRISLI LHS, namely λ, γ , E SIA and F m . As concerns the SIA-enhancement factor E SIA , in PSU an high value of E SIA = 10 is necessary to match the Hughes et al. (2016) reconstructed groundingline position at the LGM in Bjørnøyrenna Trough. Such high value is out of the range explored with the LHS with GRISLI. In contrast, the optimal values for λ, γ and F m employed in GBvol are employed to run a corresponding PSU best fit transient simulation, referred to as PBvol. Table 3         10.0 · 10 −5 7.1 · 10 −5 -Ocean basal melting parameter Fm -0.02 · 10 −3 1.56 · 10 −3 1.08 · 10 −3 1.08 · 10 −3