Glacial cycles simulation of the Antarctic Ice Sheet with PISM-Part 2: Parameter ensemble analysis

Glacial cycles simulation of the Antarctic Ice Sheet with PISM Part 2: Parameter ensemble analysis Torsten Albrecht 1, Ricarda Winkelmann 1,2, and Anders Levermann 1,2,3 1Potsdam Institute for Climate Impact Research (PIK), Member of the Leibniz Association, Potsdam, Germany 2Institute of Physics and Astronomy, University of Potsdam, Potsdam, Germany 3Lamont-Doherty Earth Observatory, Columbia University, New York, USA Correspondence to: T. Albrecht (albrecht@pik-potsdam.de)


Introduction
Sea-level estimates involve high uncertainty in particular with regard to the potential instability of marine-based parts of the Antarctic Ice Sheet (e.g., Weertman, 1974;Mercer, 1978;Slangen et al., 2017). Processed-based models provide the tools to evaluate the currently observed ice sheet changes (Shepherd et al., 2018a, b), to better distinguish between natural drift, 20 variability or anthropogenic drivers (Jenkins et al., 2018) and to estimate future changes for possible climatic boundary conditions (Oppenheimer and Alley, 2016;Shepherd and Nowicki, 2017;Pattyn, 2018). Regarding the involved variety of uncertain parameters and boundary conditions, confidence of future projections from such models is strengthened by systematic validation against modern observations and past reconstructions. We can build on experience gained in several preceding Antarctic modeling studies (Briggs et al., , 2014Whitehouse et al., 2012a;Golledge et al., 2014;Maris et al., 2014Maris et al., , 2015Pollard 25 et al., 2016Pollard 25 et al., , 2017Quiquet et al., 2018), providing paleo dataset compilations or advanced scoring schemes. Modern datasets encompass ice thickness, grounding line and calving front position (Bedmap2, Fretwell et al., 2013), surface velocity (Rignot et al., 2011) as well as uplift rates from GPS measurements (Whitehouse et al., 2012b). Reconstructions of grounding line location at the Last Glacial Maximum (LGM) as provided by the RAISED Consortium (Bentley et al., 2014) are used as paleo constraints as well as grounding line locations and cosmogenic elevation-age data from the AntICEdat database (Briggs and 30 Tarasov, 2013) at specific sites during the deglaciation period.
In this study we run simulations of the entire Antarctic Ice Sheet with the Parallel Ice Sheet Model (PISM, Winkelmann et al., 2011;The PISM authors, 2017). The hybrid of two shallow approximations of the stress balance and the comparably coarse resolution of 16 km allow for running an ensemble of simulations of ice sheet dynamics over the last two (dominant) In all ensemble runs we used for the SSA stress balance an enhancement factor of 0.6 (see Sect. 2.1 in companion paper), which is relevant for ice stream and ice shelf regions.
PPQ: Exponent q used in "pseudo plastic" sliding law which relates bed-parallel basal shear stress τ b to sliding velocity u b in the form 95 as calculated from the Shallow Shelf Approximation (SSA) of the stress balance (Bueler and Brown, 2009), for threshold speed u 0 and yields stress τ c . The sliding exponent hence covers uncertainties in basal friction. Values are 0.25, 0.5, 0.75 and 1.0 (non PREC: Precipitation scaling factor f p according to temperature forcing ∆T motivated by Clausius-Clapeyron-relationship and data analysis (Frieler et al., 2015), which can be formulated as exponential function (Ritz et al., 1996;Quiquet et al., 2012) as 100 P (t) = P 0 exp (f p ∆T (t)) ≈ P 0 (1.0 + f p ∆T (t)) .
For given present-day mean precipitation field P 0 , the factor f p captures uncertainty in the climatic mass balance, in particular for glacial periods. Values are 2, 5, 7 and 10 %/K. VISC: Mantle viscosity determines the characteristic response time of the linearly viscous half-space of the Earth to changing ice and adjacent ocean loads (Bueler et al., 2007, Eq. 1). It covers uncertainties within the Earth model for values of 2. TOTI: Mean-square error mismatch of present-day floating ice shelf areas to observations (Fretwell et al., 2013) assuming uncertainty in grounding line and calving front location of 30 km, according to Pollard et al. (2016, Appendix B2).
3. TOTDH: Mean-square-error model misfit of present-day state to observed ice thickness (Fretwell et al., 2013) with respect to an assumed observational uncertainty of 10 m and evaluated over the contemporary grounded region, close to Pollard et al. (2016, Appendix B3). 125 4. TOTGL: Mean-square-error misfit to observed grounding line location for the modeled Antarctic grounded mask (ice rises excluded) using a two-dimensional distance field approximation 1 . This method is different to the GL2D constraint used in Pollard et al. (2016, Appendix B5), and is only applied to the present-day grounding line around the whole Antarctic Ice Sheet according to Fretwell et al. (Bedmap2;, while considering observational uncertainty of 30 km as in TOTI and TOTE above.  , Table S2) including individual observational uncertainty. Misfit is evaluated for the closest model grid point as in Pollard et al. (2016, Appendix B8), including intra-data type weighting (Briggs and Tarasov, 2013, Sect. 4.3.1).
6. TOTVEL: Mean-square error misfit in (grounded) surface ice speed compared to a remapped version of observational 135 data by Rignot et al. (2011) including their provided grid-cell wise standard deviation, bounded below by 1.5 m/yr.
Paleo-data type constraints are partly based on the AntICEdat compilation by Briggs and Tarasov (2013, Sect. 4.2), following closely their model-data misfit computation. Their compilation also includes records of regional sea-level change above present-day elevation (RSL), which has not been considered in this study, as PISM lacks a self-consistent sea-level model to account for regional self-gravitational effect of the order of up to several meters, which can be similar to the magnitude been taken into account here, as they would favor simulations that reveal a rather slow and progressive grounding-line retreat through the Holocene in both Ross and Ronne Ice Shelf, which has not necessarily been the case (Kingslake et al.,8. ELEV: Mean of squared misfit of past (cosmogenic) surface elevation vs. age in the last 120 kyr based on model-data differences at 106 individual sites (distributed over 26 regions, weighted by inverse areal density, see Sect. 4.3.1 in Briggs and Tarasov (2013)). For each data-point the smallest misfit to observations is computed for all past ice surface elevations (sampled every 1 kyr) of the 16 km model grid interpolated to the core location and datum as part of a thinning trend 155 (Briggs and Tarasov, 2013, Sect. 4.2). A subset of these data has been also used in Maris et al. (2015) and Pollard et al. (2016, Appendix B7).
9. EXT: Mean of squared misfit of observed ice extent at 27 locations around the AIS in the last 28 kyr with dates for the onset of open marine conditions (OMC) or grounding-line retreat (GLR). The modeled age is computed as the most recent transition from grounded to floating ice conditions considering the sea-level anomaly. The model output every 160 1 kyr is interpolated down to core location and linearly interpolated to 100 yr temporal resolution, while weighting is not necessary here, as described in Briggs and Tarasov (2013, Sect. 4.2). A subset of these data has been also used in Maris et al. (2015) 2.3 Score aggregation Each of the misfits above are first transformed into a normalized individual score for each data type i and each run j using 165 the median over all misfits M i,j for the 256 simulation. The procedure closely follows Approach (A) in Pollard et al. (2016, Sect. 2.4.1). Then the individual score S i,j is normalized according to the median to  The score values are computed versus geologic and modern data sets, normalized by the best score in the ensemble, and range from <0.01 (bright yellow, no skill) to 1.0 (dark red, best score), on a logarithmic color scale (cfs. Pollard et al., 2016, Figs. 2 + C1). The four parameters are the SIA enhancement factor ESIA (outer y-axis), the temperature-dependent precipitation scaling PREC (outer x-axis), the mantle viscosity VISC (inner y-axis) and the power-law sliding exponent PPQ (inner x-axis).
The parameter ESIA enhances the shear-dominated ice flow and hence yields ice thinning particularly in the interior of the ice sheet and therewith a decrease in the total ice volume. ESIA values of 4.0 or 7.0 have been used in other models (e.g., Maris et al., 2015) to compensate for overestimated ice thickness in the interior of East Antarctic Ice Sheet under present-day climate conditions. In our ensemble, we find a trend towards higher scores for small ESIA values of 1.0 or 2.0 (in the upper 185 half section of Fig. 1). This becomes more prominent when considering ensemble-mean score shares for individual parameter values as in Fig. 3, with a normalized mean score of 46% for ESIA = 1.0 as compared to a mean score of 6% for ESIA = 7.0.
Most of this trend is a result of the individual data-type score TOTDH (see Fig. 5, column 4, row 3) as it measures the overall misfit of modern ice thickness (and volume distribution). Partly this trend can be also attributed to the TROUGH data-type scores (Fig. 5, column 8, row 3), as for higher ESIA values grounding line motion tends to slow down, such that the time span 190 between LGM and present is not sufficient for a complete retreat back to the observed present-day location, at least in some ice shelf basins. The best-score ensemble members for small ESIA values are found in combination with both high values of mantle viscosity VISC and high values of friction exponent PPQ (center column panels in Fig. 4).

195
Regarding the choice of the precipitation scaling PREC the best-score members are found at the upper sampling range with values of 7 %/K or 10 %/K (see right half section in Fig. 1). Considering normalized ensemble-mean score for individual parameter values over the full range of 2-10 %/K, we can find a trend from 13% to 42% (see lowest panel in Fig. 3). Regarding parameter combinations with PREC (left-hand column in Fig. 4), we detect a weak trend towards lower ESIA and higher PPQ, while individual data-type scores (lower row in Fig. 5) show a rather uniform pattern, in particular regarding the misfit to 200 present-day observations. As the PREC parameter is linked to the temperature anomaly forcing, it affects the ice volume and hence the grounding line location particularly for temperature conditions different from present day. This suggests a stronger signal of PREC parameter variation in the paleo data-types scores.

205
A more complex pattern is found for PPQ in each of the sub-panels of Fig Regarding mantle viscosity VISC, scores are generally low with 9% for the smallest sampled value of the parameter 215 VISC =0.1 × 10 21 Pa s, while best scores are found in the ensemble for the five times larger viscosity of VISC =0.5 × 10 21 Pa s, with 44%. In our model, the mantle viscosity parameter has been applied to the whole Antarctic continent, although observations in some localized regions as in the Amundsen Sea suggest that upper mantle viscosities could be considerably smaller than the tested range, up to the order of 10 19 Pa s (Barletta et al., 2018). For the upper range of tested mantle viscosities, up to VISC =10.0 × 10 21 Pa s, we find a normalized ensemble mean of 27% and 20%, respectively. Note that VISC parameter 220 values have been sampled non-linearly over a range of two orders of magnitude. For the lowest value there is a clear trend towards smaller scores in the grounding-line and ice-thickness related data-types, such as TOTE, TOTGL, TROUGH and TOTDH respectively. As mantle viscosity determines the rate of response of the bed to changes in ice thickness a low viscosity corresponds to a rather quick uplift after grounding-line retreat and hence to a retarded retreat, which corresponds to a rather extended present-day state. This implies smaller ice shelves with slower flow and less velocity misfit, such that also TOTVEL 225 favors small VISC values. In contrast, a trend to rather high mantle viscosities in the aggregated score stems mainly from the misfit of present-day uplift rates expressed as data-type score TOTUPL, probably due to reduced sensitivity to fluctuations in grounding line location. High mantle viscosities involve a slow bed uplift and grounding-line retreat can occur faster. More specifically, in the partially over-deepened ice shelf basins, which have been additionally depressed at the Last Glacial Maximum by a couple of hundred meters as compared to present, grounding-line retreat can amplify itself in terms of a regional 230 Marine Ice Sheet Instability (Mercer, 1978;Schoof, 2007;Bart et al., 2016). In fact, the best score ensemble members are found for intermediate mantle viscosities of VISC =0.5 × 10 21 Pa s, and VISC =2.5 × 10 21 Pa s. This could be a result of the product formulation of the aggregated score, in which individual data types scores favor opposing extreme values.

235
The five best-score ensemble members and associated parameter combinations are listed in Table 1. With the best-fit simulation parameters we have participated in the initMIP-Antarctica model intercomparison (Seroussi et al., 2019, PISMPAL3).
The individual scores with respect to the nine data-types are visualized for the 20 best ensemble members in Fig. 2. The scores associated with the paleo data-types ELEV and EXT show only comparably little variation among the ensemble (both around 0.07 standard deviation). This also applies for the present-day ice shelf area mismatch TOTI (0.04), as no calving parame-240 ter has been varied. In contrast, present-day data types associated with velocity (TOTVEL) and uplift rates (TOTUPL) show strong variations among the 20 best ensemble members, with a standard deviation in score across the entire ensemble of 0.18 and 0.30 respectively. For data types that are related to grounding line position (TOTGL, TOTE, TROUGH) and ice volume (TOTDH) we find a similar order as for the TOTAL aggregated score ( Table 1. Five best-score ensemble parameter combinations with parameter values and total scores. The best-fit simulation parameters were used in the initMIP-Antarctica model intercomparison (Seroussi et al., 2019, PISMPAL3) and for the reference simulation in the companion paper (Albrecht et al., 2019).
across the basins in time, with distances of up to 1000 km, leaving behind the large floating ice shelves (Fig. 7). In about 10% of the score-weighted simulations grounding line remains at the extended position without significant retreat, linked to 255 an efficient negative feedback on grounding line motion, related to a fast responding bed (low VISC). In contrast, for rather low friction and high mantle viscosities, we find fast grounding-line retreat, with a stabilization of grounding line position at or even inland of the observed location in 50% or 75% of the score-weighted simulations in the Ross and Weddell Sea sector, respectively ( Fig. 8, upper panels). Due to the grounded ice retreat and the consequent unloading across the large ice shelf basins, the marine bed lifts up by up to a few hundred meters, which can lead to grounding line re-advance supported by the 260 formation of ice rises (Kingslake et al., 2018). The ensemble mean re-advance is up to 100 km, while some of the best-score  Red dots indicate the best-score member, green and blue the second and third best ensemble members (see Table 1). Grey-dashed line indicates mean score tendency over sampled parameter range.    Grounding line in the Ross Sea and Weddell Sea sector is found inland of its present location (vertically dotted) within the last 10 kyr simulation time in about 50% and 75% of score-weighted simulations, respectively. The score-weighted mean curve (black) reveals readvance of the grounding line of up to 100 km in about 20% of the score-weighted simulations, both in the Ross and Weddell Sea sector, as discussed in (Kingslake et al., 2018  depending on Greenland's contribution (Cuffey and Marshall, 2000;Tarasov and Peltier, 2003;Kopp et al., 2009). This value has thus been used as lower bound in terms of a "sieve" criterion in a previous Antarctic model ensemble analysis (Briggs et al., 2014).

285
Assuming, that the memory of the previous spin-up has vanished at the Last Glacial Maximum (in our simulations at around 15 kyr BP), the model ensemble yields a range of (grounded) Antarctic Ice Sheet volume of 9.4 ± 4.1 m above present-day observations, or 6.5 ± 2.0 ×10 6 km 3 . The histogram of score-weighted sea-level anomalies of all simulations at Last Glacial Maximum actually reveals four distinct maxima at around 4.5, 8.1, 9.0 and 13.0 m SLE (Fig. 10 a), which can be attributed to the five best-score simulations in Table 1. The ensemble spread is hence relatively wide, but still quite symmetric, as compar- Most of the deglacial retreat from LGM extent and hence most of Antarctica's sea-level rise contribution occurs in our simulations after 10 kyr BP (cf. Fig. 10 b, c). In particular, for higher mantle viscosities we find episodic self-amplified retreat with change rates of more than 0.5 cm SLE per year in West Antarctic basins (in the best-fit simulation at 7.5 kyr BP, see 300 below in Sect 3.3). This leads in some cases to grounding-line retreat beyond its present location and subsequent re-advance during Holocene, due to the uplift of the bed (discussed in Kingslake et al., 2018). However, these rapid episodes of retreat occur in our simulations consistently after Meltwater pulse 1A (MWP1a, around 14.5 kyr BP, see dashed line in Fig. 9). This delay supports the idea, that Antarctic Ice Sheet retreat has not been a source but rather a consequence of the relatively quick rise in global mean sea-level by about 15 m within 350 yr or ≈ 4 cm/yr at MWP1a (Liu et al., 2016), while core analysis of 305 iceberg-rafted debris suggest earlier and stronger recession of the Antarctic Ice Sheet at the time of MWP1a (Weber et al., 2014 hence to stronger sub-shelf melt (Golledge et al., 2014;Fogwill et al., 2017). As our sub-shelf melting module is forced with a modified surface temperature anomaly, PICO responds with less melt during the ACR period and hence prohibits significant ice sheet retreat. But even if the intermediate ocean temperature would rise by 1 or 2 K during ACR, the induced additional melt would correspond to less than -1 mm/yr SLE and hence far less than the ice volume change rate of -6 mm/yr SLE found by

320
The timing of deglaciation and possible rebound effects can explain a natural drift in certain regions that lasts through the Holocene until present-day. In the score-weighted average the ensemble simulations suggest a sea-level contributions over the last 3,000 model years of about 0.25 mm/yr, while for the reference simulation the Antarctic ice above flotation is on average even slightly growing (cf. Fig. 9 c), partly explained by net uplift in grounded areas (Fig. 12).   2016)). Violet lines indicate simulations over the penultimate glacial cycle with four different mantle viscosities, from which the full ensemble branches at 125 kyr BP. During deglaciation the score-weighted ensemble mean (green) shows most of the sea-level change rates (lower panels) between 9 kyr BP and 5 kyr BP with mean rates above 1 mm yr −1 , while the best-score simulation (red) reveals rates of sea-level rise of up to 5 mm yr −1 (100 yr bins) in the same period (cf. Golledge et al., 2014, Fig. 3 d). In contrast to the ensemble mean, the best-score member (red line) shows a minimum ice volume in the mid-Holocene (around 4 kyr BP) and subsequent regrowth. The simulations are based on the Bedmap2 dataset (Fretwell et al., 2013), remapped to 16 km resolution, which corresponds 325 to a total grounded modern Antarctic Ice Sheet volume of 56.85 m SLE (or 26.29 ×10 6 km 3 ). The ensemble mean at the end of the simulations (in the year 2000 or -0.05 kyr BP) underestimates the observed ice volume slightly by 0.6 ± 3.5 m SLE, or in terms of grounded ice volume by 0.7 ± 1.7 ×10 6 km 3 (see Fig. 9). The histogram of score-weighted sea-level anomalies at the end of all simulations can be well approximated by a normal distribution (Fig. 10 d). As for the LGM ice volume the ESIA parameter is also responsible for most of the present-day ice volume range of the ensemble with more than 10 m SLE for 330 the covered parameter range, while PREC has almost no effect with less than 1 m SLE on average, in contrast to the LGM, as expected. VISC and PPQ reveal on average a range for the present-day ice volume of about 6 m SLE and 5 m SLE, respectively.

Comparison of LGM sea-level estimates in previous studies
For the maximum Antarctic ice volume at the Last Glacial Maximum, the inferred ensemble range of 5.3 -13.5 m SLE excess relative to observations (or 4.5 -8.5 ×10 6 km 3 ) is at the upper range found in the recent literature (Fig. 11), except for the 335 "GRISLI" model results (Quiquet et al., 2018). The other previous model reconstructions are based on four different models:   Golledge et al. (2014) and this study provided anomalies based on the volume-above-flotation calculation (VAF), the corresponding values are smaller than the directly converted values, that still include the marine part below flotation (Fig. 11b). For a conversion factor of c = 2.5 our study would yield 11.3 -21.3 m SLE instead. For the LGM ice volume excess relative to the modeled present day our study yields 5.9 -14.1 m SLE (or 3.8 -7.8 ×10 6 km 3 ), both indicated in Fig. 11.

Best-fit ensemble simulation
The best-fit ensemble member simulation (no. 165, see Table 1) provides an Antarctic Ice Sheet configuration for the present day, which is comparably close to observations. Yet, the present-day ice volume of the West Antarctic Ice Sheet is overestimated (by around 25%), while the much larger East Antarctic Ice Sheet (EAIS) volume is rather underestimated (by around 5%), which is also valid for the ensemble mean (Fig. 6). Part of the overestimation can be explained by the relatively coarsely 375 resolved topography of the Antarctic Peninsula and weakly constrained basal friction in the Siple Coast and Transantarctic Mountain area. This results in a root-mean-square error (RMSE) of ice thickness of 266 m (see Fig. 12 a), a RMSE of grounding line distance of 67 km (see Fig. 13) and a RMSE for surface velocities of 66 m/yr (see Fig. 14). The best-fit simulation also   reproduces the general pattern of observed modern isostatic adjustment rates (see Fig. 12 b) with highest uplift rates of more than 10 mm yr −1 in the Weddell and Amundsen Sea Region in agreement with GIA model reconstructions (cf. Argus et al.,380 2014, Fig. 6). In contrast to these GIA reconstructions, our best-fit simulation shows depression rather than uplift in the Siple Coast regions as grounded ice is still re-advancing and hence adding load.  At the Last Glacial Maximum around 15 kyr BP the sea-level relevant volume history of the best-score simulation is close to the ensemble mean ( Fig. 9), which agrees well with reconstructions by the RAISED Consortium (Bentley et al., 2014, cf. Fig. 7 a). The LGM state is characterized by extended ice sheet flow towards the outer Antarctic continental shelf edges, with 385 more than 2,000 m thicker ice than today in the basins of the largest modern ice shelves (Ross, Weddell, Amery and Amundsen), while the inner East Antarctic Ice was a few hundred meters thinner than today (see Fig. 15).
Even though this is not the primary focus of this parameter ensemble study, it is worthwhile to have a closer look into the deglacial period. The last glacial termination (also known as Termination I, which is the end of Marine isotope stage 2), and 390 hence the period of major ice sheet retreat, initiates in our best-score simulation in the Ross and Amundsen sector at around 9 kyr BP, in the Amery sector at around 8 kyr BP and in the Weddell Sea Sector at around 7 kyr BP. Maximum ice volume change rates are found accordingly in the period between 10 and 8 kyr BP with on average -1.4 mm/yr SLE (or -660 Gt/yr) and in the period between 8 and 6 kyr BP with on average -2.4 mm/yr SLE (or -1,300 Gt/yr, Fig. 16). In the 100 yr running mean of the ice volume change rate we find a peak of around -5 mm/yr SLE at 7.5 kyr BP (or -3,300 Gt/yr, compare black and 395 khaki line in Fig. 17). This rate of change is significantly larger than in the ensemble mean with up to -2 mm/yr SLE, as the mean retreat becomes smoothed over a longer deglacial period (see Fig. 9 c). The total ice volume change during the period 10-5 kyr BP in the best-fit simulation amounts to -9.7 m SLE. Most of this change can be attributed to increased discharge by around 1,000 Gt/yr and increased sub-shelf melting by around 450 Gt/yr (partly due to increased floating ice shelf area), while surface mass balance increased only by around 300 Gt/yr (Fig. 17). Recent proxy-data reconstructions from the eastern 2017), which is consistent with our model simulations. In the reconstructions by the RAISED Consortium most of the retreat in the Ross Sea Sector (almost up to present-day grounding line location) occurred between 10 kyr BP and 5 kyr BP, while major retreat in the Weddell Sea Sector likely happened before 10 kyr BP in scenario A and after 5 kyr BP in scenario B (Bentley et al., 2014, cf. Fig. 7 b,c).

405
A Holocene minimum ice volume is reached in our simulations around 3 kyr BP with slight re-advance and thickening in the Siple coast and Bungenstock ice rise until present-day (see Fig. 16). This regrowth signal cannot be inferred from RAISED reconstructions with snapshots only every 5 kyr (Bentley et al., 2014). The corresponding mass change is rather small with 60 Gt/yr (or 0.07 mm/yr SLE) in the last 3,000 years, see Fig. 17. During this late Holocene period, surface mass balance of 410 around 3,700 Gt/yr is balanced by approximately 2,600 Gt/yr discharge, while sub-shelf melt plays a minor role with around 1,000 Gt/yr.  Figure 15. Snapshots of grounded ice thickness anomaly to present-day observations (Fretwell et al., 2013) over last 15 kyr in best-fit simulation, analogous to Fig. 2 in Golledge et al. (2014). At LGM state grounded ice extends towards the edge of the continental shelf with much thicker ice than present, mainly in West Antarctica. Retreat of the ice sheet occurs first in the Ross basin between 9 and 8 kyr BP, followed by the Amery basin around 1 kyr later and the Amundsen and Weddel Sea basin between 7 and 5 kyr BP. East Antarctic Ice Sheet thickness is underestimated throughout the deglaciation period (light blue shaded area).

Discussion of individual ensemble parameters
In this section we want to discuss the effects of individual ensemble parameters in more detail, also in comparison to previous model studies. We performed our analysis for an ensemble of 256 simulations of the entire Antarctic Ice Sheet over the last two  0). Interestingly, we find best scores for rather high sliding exponents of PPQ with values of 0.75 or 1.0 (rather linear relationship of sliding velocity and till strength).  used Coulomb friction and varied instead three parameters that control the basal sliding over soft and hard beds, based on an erosion parameterization. In our study, the till weakness is associated with the till friction angle, which is optimized for the present-day grounded Antarctic Ice Sheet (Pollard and DeConto, 2012b). In , basal 435 sliding additionally accounts for basal roughness and pinnings points (three parameters), which is otherwise underestimated as a result of the coarse model resolution.
Another sliding-related key parameter is the friction coefficient underneath the modern ice shelves as varied in Pollard et al. (2016Pollard et al. ( , 2017, who found it to be the most dominant ensemble parameter. As discussed in the companion paper (Albrecht et al., 2019), we also find till properties in the ice shelf regions highly relevant, in particular during deglaciation. As a consequence, 440 we have run an additional ensemble analysis for four parameters associated with basal sliding and hydrology, including friction underneath modern ice shelves and discussed the results in the Supplementray Material A. In the best fit simulations of this "basal ensemble", we find main deglacial retreat occurring a few thousand years earlier (closer to MWP-1A) than in the base ensemble. Hence, the corresponding scores are even better than for the best fit simulation of the base ensemble, for same sliding exponent but smaller minimal till friction angle of ≤ 1 • .

445
For a representation of the ice dynamical uncertainty we chose the ESIA enhancement factor as most relevant ensemble parameter, which mainly affects the grounded ice volume. We find best fits for rather hancement and a reference value for SSA enhancement close to 0.6 (see Table 1 in , which we have used in all our ensemble simulations. Maris et al. (2014) determined in their sensitivity study for the SIA enhancement an even larger reference value of 9.0 and for the SSA enhancement 0.8. In Quiquet et al. (2018) best fits to present-day thickness are found 455 for SIA enhancement between 1.5 and 4.0, for SSA enhancement between 0.2 and 0.5 respectively.
As the climate-related uncertain ensemble parameter we chose a parameter associated with the change of precipitation with temperatures, PREC. The best-fit parameter values of PREC = 7-10 %/K yield for 10 K colder glacial temperatures about 50-65 % less precipitation. This parameter is similar to the insolation scaling parameter in , where the 460 best fit value would result in about 60% less precipitation at insolation minimum. In total, Briggs et al. (2014) varied seven precipitation-related parameters based on three different precipitation forcings (one of which is similar to the one we used). Maris et al. (2014) used instead a linear temperature-based scaling between LGM and present-day surface mass balance (with about 58% anomaly) with a fixed parameter.

465
Beyond the four parameters varied in our ensemble, previous ensemble studies found for instance a high sensitivity in at least one of the five temperature related parameters . In contrast, we found only little effect of temperature related parameter variation on the sea-level relevant ice sheet volume, as discussed in the companion paper (Albrecht et al., 2019). Concerning iceberg calving, Pollard et al. (2016) included one related parameter in their ensemble analysis, while  varied three parameters for ice shelf calving and one parameter for tidewater calving. Our 'eigen-470 calving' parameterization also uses a strain-rate based calving estimate, combined with a minimal terminal ice thickness and provides a representation of calving front dynamics, which, to the first order, yields calving front locations close to present observations (Levermann et al., 2012). As this parameterization is rather independent of the climate conditions, variations of the 'eigencalving' parameter show only little effect on sea-level relevant ice volume (see Albrecht et al. (2019), Sect. 2.4).
The four selected ensemble parameters, representing uncertainties in interacting ice-Earth dynamics, basal sliding as well climate conditions, imply a large range of uncertainty for the total Antarctic ice volume change. They have been chosen, 485 such that the model yields a present-day ice volume close to observations, while the LGM ice volume differs significantly for parameter change. The probed parameter range has been chosen to be rather broad, which implies a low sampling density of the parameter space. With the knowledge gained in this ensemble analysis, this range could be further constrained in a (larger) sub-ensemble. Also, other parameters may be more relevant for certain regions of the Antarctic Ice Sheet or for the onset and rate of the last deglaciation, which in our ensemble occurs generally later than suggested by many paleo records. A closer look 490 into the details of the deglacial period and relevant parameters will be discussed in a separate follow-up study.
One deficiency of our model settings is the general underestimation of ice thickness in the inner ice sheet sections up to -500 m mainly in the EAIS, which could be a result of the underestimated RACMO precipitation. In contrast, ice thickness is overestimated in the outer terminal regions and at Siple Coast by up to 500 m, where the complex topography is not sufficiently resolved in the model, with implications for inferred basal conditions and temperature conditions. Accordingly, we find a 495 considerable misfit to most paleo elevation data (ELEV), which are located mainly in the marginal mountain regions. This could be improved, e.g. by parameterized basal roughness or erosion, as proposed in , or by higher model resolution and updated bed topography data sets 2 .
The score aggregation scheme according to Pollard et al. (2016) implies that the paleo data types have equal influence as the present-day constraints, although they cover only relatively small regions and periods of the modeled ice sheet history (Briggs 500 and Tarasov, 2013). However, as the variability in paleo misfits is comparably low among the ensemble, these data types have only relatively small imprint on the aggregated score (see more details in Supplementary Material B). This is also valid for a data-type weighted score , which applied to our ensemble results yields a similar set of best score runs.
Further work will consist in the determination of more realistic climate reconstructions using general circulation model re-505 sults and in the explicit computation of the local relative sea-level, which could potentially have an strong impact on grounding line migration for glacial cycles (Gomez et al., 2013).

Conclusions
We have run an ensemble of 256 simulations over the last two glacial cycles and have applied a simple averaging method with full factorial sampling similar to Pollard et al. (2016). Although this kind of ensemble method is limited to a comparably small 510 number of values for each parameter and hence the retrieved scores are somewhat blocky (as compared to advanced techniques that can interpolate in parameter space) we still recognize a general pattern of parameter combinations that provide best model fits to both present-day observations and paleo records. However, the selected ensemble parameters certainly cannot cover the full range of possible model response, in particular with regard to the self-amplifying effects during deglaciation.
For the four sampled parameters, best fits are found for mantle viscosity around VISC = 0.5-2.5×10 21 Pa s, rather lin- The grounding line extends at last glacial maximum to the edge of the continental shelf for nearly all simulations. The onset and rate of deglaciation, however is very sensitive to the choice of parameters and boundary conditions, in particular 525 those related to basal sliding. Due to the comparably coarse resolution and the high uncertainty (sensitivity) that comes with the strong non-linearity of the system, we here discuss rather general patterns of reconstructed ice sheet histories than exact numbers, which would require a much larger ensemble with an extended number of varied parameters.
The score-weighted likely range (one standard deviation) of our reconstructed ice volume histories suggest a contribution of the Antarctic Ice Sheet to the global mean sea level since the Last Glacial Maximum at around 15 kyr BP of 9.4 ± 4.1 m than in most recently published studies. The choice of basal sliding parameterization in the different models seems to have most impact on the corresponding estimates. The ensemble reproduces the observed present-day grounded ice volume with an score-weighted anomaly of 0.6 ± 3.5 m SLE (0.7 ± 1.7 ×10 6 km 3 ) and hence serves as a suitable inital state for future projections.

535
The reconstructed score-weighted ensemble range (1σ) is comparably large with up to 4.3 m SLE (or 2.0 ×10 6 km 3 ), which can be explained by a high model sensitivity (Albrecht et al., 2019), by a comparably large range of the sampled parameter values and of course due to the choice of the aggregated score scheme. By using "sieve" criteria the ensemble range could be reduced. For the much larger ensemble study by covering 31 parameters Briggs et al. (2014) a narrowed ensemble range of 4.4 mESL (different definition of sea-level equivalent volume change) or 1.8 ×10 6 km 3 was found for the best 5% of the 540 ensemble simulations, which is close to the range of our study.
The onset of deglaciation and hence major grounding-line retreat occurs in our model simulations after 12 kyr BP and hence well after MWP1a (≈14.3 kyr BP). A previous PISM study simulated much earlier and larger sea-level contributions from Antarctica for oceanic forcing at intermediate levels that is anticorrelated to the surface temperature forcing (Golledge et al., 2014), as likely happened during the two millennia of Antarctic Cold Reversal following the MWP1a.

545
The PISM model results in Kingslake et al. (2018) are based on this ensemble study, but have been published before with a slightly older model version (see data and model code availability therein). Meanwhile, we have improved the Earth model, which accounts for changes in the ocean water column induced by variations in bed topography or sea-level changes. In contrast to Kingslake et al. (2018), we used the remapped topography without local adjustment in the region of Bungenstock ice rise in the Weddell Sea sector, and found in about 20% of the score-weighted simulations an extensive retreat of the grounding line 550 and subsequent re-advance in both the Ross and Weddell Sea sector.
The paleo simulation ensemble analysis presented here provides a set of data-constrained parameter combinations for PISM simulations, that can be used as a reference for further sensitivity studies investigating specific episodes in the history of Antarctica, such as the the last deglaciation or the Last Interglacial, as well as for projections of Antarctic sea-level contributions.
Code and data availability. The PISM code used in this study can be obtained from https://github.com/talbrecht/pism_pik/tree/pism_pik_1.0 Acknowledgements. Development of PISM is supported by NASA grant NNX17AG65G and NSF grants PLR-1603799 and PLR-1644277.
Open-source software was used at all stages, in particular NumPy (www.numpy.org), CDO (https://code.mpimet.mpg.de/projects/cdo/), NCO (http://nco.sourceforge.net/) and matplotlib (https://matplotlib.org/). The authors gratefully acknowledge the European Regional Development Fund (ERDF), the German Federal Ministry of Education and Research and the Land Brandenburg for supporting this project by providing resources on the high performance computer system at the Potsdam Institute for Climate Impact Research. Computer resources for 565 this project have been also provided by the Gauss Centre for Supercomputing/Leibniz Supercomputing Centre (www.lrz.de) under Project-ID pr94ga and pn69ru. T.A. is supported by the Deutsche Forschungsgemeinschaft (DFG) in the framework of the priority program "Antarctic Research with comparative investigations in Arctic ice areas" by grant LE1448/6-1 and LE1448/7-1. We thank Dave Pollard for sharing ensemble analysis scripts and for valueable discussions and Stewart Jamieson for sharing gridded RAISED datasets. Finally, we appreciate the assistence by the editor Alexander Robinson and the many helpful suggestions and comments by Lev Tarasov and an anonymous reviewer, 570 which led us to considerable improvements of the manuscript.