Exploring the impact of atmospheric forcing and basal boundary conditions on the simulation of the Antarctic ice sheet at the Last Glacial Maximum

Little is known about the distribution of ice in the Antarctic ice sheet (AIS) during the Last Glacial Maximum (LGM). Whereas marine and terrestrial geological data indicate that the grounded ice advanced to a position close to the continentalshelf break, the total ice volume is unclear. Glacial boundary conditions are potentially important sources of uncertainty, in particular basal friction and climatic boundary conditions. Basal friction exerts a strong control on the large-scale dynamics of the ice sheet and thus affects its size, and is not well constrained. Glacial climatic boundary conditions determine the net 5 accumulation and ice temperature, and are also poorly known. Here we explore the effect of the uncertainty in both features on the total simulated ice storage of the AIS at the LGM. For this purpose we use a hybrid ice-sheet-shelf model that is forced with different basal-drag choices and glacial background climatic conditions obtained from the LGM ensemble climate simulations of the third phase of the Paleoclimate Modelling Intercomparison Project (PMIP3). For a wide range of plausible basal friction configurations, the simulated ice dynamics vary widely but all simulations produce fully extended ice sheets 10 towards the continental-shelf break. More dynamically active ice sheets correspond to lower ice volumes, while they remain consistent with the available constraints on ice extent. Thus, this work points to the possibility of an AIS with very active ice streams during the LGM. In addition, we find that the surface boundary temperature field plays a crucial role in determining the ice extent through its effect on viscosity. For ice sheets of a similar extent and comparable dynamics, we find that the precipitation field determines the total AIS volume. However, precipitation is highly uncertain. Climatic fields simulated by 15 climate models show more precipitation in coastal regions than a spatially uniform anomaly, which can lead to larger ice volumes. We strongly support using these paleoclimatic fields to simulate and study the LGM and potentially other time periods like the Last Interglacial. However, their accuracy must be assessed as well, as differences between climate model forcing lead to a range in the simulated ice volume and extension of about 6 m sea-level equivalent and one million km. 1 https://doi.org/10.5194/tc-2020-28 Preprint. Discussion started: 10 March 2020 c © Author(s) 2020. CC BY 4.0 License.


Introduction
Sea-level variations on long timescales are driven by the waxing and waning of large continental ice sheets. The characterisation of the sensitivity of ice sheets to past climate changes is fundamental to gaining insight into their underlying dynamics as well as their response to future climate change. In addition, understanding past sea-level changes is important for quantifying sealevel rise (Nicholls and Cazenave, 2010;Defrance et al., 2017;King and Harrington, 2018;Golledge et al., 2019;Robel et al., 5 2019) and for assessing the risk of crossing tipping points within the Earth System, such as the collapse of the West Antarctic Ice Sheet (Kopp et al., 2009;Sutter et al., 2016;Pattyn et al., 2018).
The Antarctic Ice Sheet (AIS), in particular, plays a fundamental role as it is the largest ice sheet on Earth and stores ca. 58 meters of sea-level equivalent (msle; Fretwell et al. (2013)). Due to its size it is potentially the largest contributor to future sea-level projections, but it is also the most uncertain (Collins et al., 2013). Assessing the AIS contribution to the total sea-level budget at different time periods has proven to be challenging. The Last Glacial Maximum (LGM, 21 ka BP) represents an ideal benchmark period since there is a large availability and variety of proxy data that, furthermore, indicate important AIS changes relative to present day (PD). Both, marine and terrestrial geological data, indicate that at the LGM, the AIS extended up to the continental-shelf break (Anderson et al., 2002(Anderson et al., , 2014Hillenbrand et al., 2012Hillenbrand et al., , 2014The RAISED Consortium, 2014;Mackintosh et al., 2014). However, its exact extent is not well constrained everywhere. Whereas its advance in the Amundsen 15 region, the Bellingshausen Sea and the Antarctic Peninsula is well established, in the Ross Sea and the East Antarctic region it remains controversial (Stolldorf et al., 2012;The RAISED Consortium, 2014). Furthermore, the total AIS ice volume is even less well constrained (Simms et al. (2019) and references therein). Geological data furthermore do not provide direct information on past thickness and volume of ice sheets, which must hence be inferred. There have been several approaches to infer past ice-volume change of an individual ice sheet as the AIS. One approach is to use direct ice-sheet modelling to simulate 20 the volume of the AIS at the LGM (e.g Huybrechts (2002); Whitehouse et al. (2012a); Golledge et al. (2012); Gomez et al. (2013); Maris et al. (2014); Briggs et al. (2014); Quiquet et al. (2018)). An alternative is to use Glacial Isostatic Adjustment (GIA) modelling, which describes the viscous response of the solid Earth to past changes in surface loading by ice and water (e.g. Ivins and James (2005); Bassett et al. (2007)). This approach has also been used in combination with direct ice-sheet modelling (e.g. Whitehouse et al. (2012b)) and/or by making use of constraints on ice-thickness from reconstructions based 25 on exposure age dating, as well as satellite observations of current uplift (Whitehouse et al., 2012b;Ivins et al., 2013;Argus et al., 2014b). Whereas older studies estimated large sea-level contributions generally above 15 m (e.g. Nakada et al. (2000); Huybrechts (2002); Peltier and Fairbanks (2006); Philippon et al. (2006); Bassett et al. (2007)), more recent modelling studies and reconstructions have lowered these estimates to 7.5-13.5 m (Mackintosh et al., 2011;Whitehouse et al., 2012a;Golledge et al., 2012Golledge et al., , 2014Gomez et al., 2013;Argus et al., 2014b;Briggs et al., 2014;Maris et al., 2014;Sutter et al., 2019). Several 30 factors have contributed to a decrease in the estimate of the LGM AIS volume. On one hand, the state of the art of ice-sheet modelling has considerably advanced in the last years, for example through the inclusion of more complex physics, increased spatial resolution and sub-grid scale grounding-line treatment (e.g. Goelzer et al. (2017); Pattyn (2018)). On the other hand, external processes, like the ice-ocean interaction or the GIA, are now treated with more accurate parameterisations and models 2 https://doi.org/10.5194/tc-2020-28 Preprint. Discussion started: 10 March 2020 c Author(s) 2020. CC BY 4.0 License.
Given that ablation and basal melting were probably negligible at the LGM in the AIS, ice-sheet dynamics and accumulation must have been the two main factors controlling ice-mass gain during this period. The representation of ice dynamics in icesheet models is a key feature that can potentially lead to important discrepancies. Most ice-sheet models simulating the past 5 long-term evolution of large-scale ice sheets are hybrid models that rely on the Shallow Ice Approximation (SIA) and the Shallow Shelf Approximation (SSA). Moreover, there is no universally accepted friction law, and basal friction is treated in different manners in ice-sheet models. Ritz et al. (2015) emphasize the importance of the basal friction, as it can favour the occurrence of the marine instability in future AIS projections. Generally, basal stress follows either a power-law formulation on the basal ice velocity (a special case being the Weertman (1957) friction law) or a Coulomb friction law (Schoof, 2005) with 10 different power-law coefficients, a friction coefficient and potentially a regularization term. Ice-sheet models thus use friction formulations that can range from linear viscous and regularized Coulomb friction laws, typical of hard bedrock sliding (Larour et al., 2012;Pattyn et al., 2013;Joughin et al., 2019) to Coulomb-plastic deformation, characteristic of ice flow over a soft bedrock with filled cavities (Schoof, 2005(Schoof, , 2006Nowicki et al., 2013). In the simplest cases a constant friction coefficient is prescribed over the whole domain (Golledge et al., 2012), but generally this parameter incorporates the dependency of 15 basal friction on the effective pressure exerted by the ice, as well as on bedrock characteristics by making use of assumed till properties Albrecht et al., 2019;Sutter et al., 2019) or basal temperature conditions (Pattyn, 2017;Quiquet et al., 2018). The sensitivity of the simulated ice volume to these features is substantial. For instance, Briggs et al. (2013) obtained differences of more than 5 msle for an Antarctic LGM state depending only on the friction coefficients used for hard and soft beds. Some studies have attempted to overcome the uncertainty in basal friction by optimising the friction 20 coefficient through inversion methods in order to obtain an accurate PD ice-sheet state Le clec'h et al., 2019). However, these optimizations are based on a particular configuration of the PD state, and it is unclear whether they remain valid for glacial conditions. All in all, basal friction is poorly characterised, and the potential consequences of the associated uncertainty should be considered in ice-sheet modeling.
Glacial atmospheric boundary conditions over Antarctica are also far from being well constrained. It is clear from ice-core 25 records and marine deep-sea sediment data that, at the continental scale, temperatures were lower than today and that the climate was drier (Frieler et al., 2015;Fudge et al., 2016). Typically, ice-sheet models use two approaches for simulating the atmospheric conditions at the LGM. On one hand, some studies prescribe a spatially-uniform temperature anomaly (generally between 8 K and 10 K below PD) and a uniform reduction in precipitation (generally by 40-50% compared to PD), as inferred from individual ice-core records (Huybrechts, 2002;Golledge et al., 2012;Whitehouse et al., 2012a;Gomez et al., 2013;30 Quiquet et al., 2018). However, this approach provides only a crude representation of glacial climate anomalies. In reality, even if ice cores show a similar temperature decrease, estimated precipitation changes are less homogeneous. Thus imposing a constant change over the whole domain will potentially misrepresent climatologies in localized areas (Frieler et al., 2015;Fudge et al., 2016). In addition, ice cores are extracted from domes and the recorded changes are not necessarily representative of coastal regions. Because the LGM is a cold state, with presumably no (or negligible) ablation and oceanic basal melt, the 35 reduction of precipitation with respect to the PD should have an important impact on the size of the simulated ice sheet.
In addition, because the temperature and/or precipitation anomalies are uniform, the PD pattern is imprinted on the LGM atmospheric forcing fields, and changes in atmospheric patterns are thus neglected.
Another commonly used method is to prescribe the LGM temperature and precipitation fields for the whole Antarctic domain from climate simulations (Briggs et al., 2013;Maris et al., 2014;Sutter et al., 2019). Output from simulations using 5 a hierarchy of climate models has been used in the literature, from global general circulation models (GCMs) (Sutter et al., 2019), sometimes downscaled with regional models (Maris et al., 2014), to Earth System Models of Intermediate Complexity (EMICs) (Blasco et al., 2019). Briggs et al. (2013) went a step forward to investigate the effect of uncertainty in the climate forcing fields by assessing the effect of the inter-model variance through an empirical orthogonal function (EOF) analysis. However, some model outputs do not simulate the temperature anomalies correctly at specific sites where proxies are available, 10 such as Vostok or Dome C. This may lead to an unrealistic configuration and thus it is necessary to evaluate the accuracy of model outputs (Cauquoin et al., 2015).
In this work we aim to assess the effects of the uncertainty in basal friction and climatic (in particular atmospheric) boundary conditions on the simulated LGM AIS. We focus on basal-drag choices which can lead to realistic LGM states. For these we then investigate the effect of different temperature and precipitation fields. To this end, we use a thermomechanical ice-sheet-

Methods and experimental setup
For this study we use the three-dimensional, hybrid, thermomechanical ice-sheet-shelf model Yelmo . where the computed basal velocity is corrected with the corresponding basal friction. Ice shelves are solved within the SSA solution without basal drag. The initial topographic conditions (ice thickness, surface and bedrock elevation) are obtained from the RTopo-2 dataset (Schaffer et al., 2016). The internal ice temperature is calculated via the advection-diffusion equation.
Yelmo computes the total mass balance (MB) as a sum of the surface mass balance (SMB), the basal mass balance at the ice 30 base and calving at the ice front. The SMB is obtained as a difference between the ice accumulation through precipitation and surface melting using the positive degree-day method (PDD; Reeh (1989)). Although there are more comprehensive methods that account for short-wave radiation for instance (Robinson et al., 2011), the PDD scheme is commonly used in ice models 4 https://doi.org/10.5194/tc-2020-28 Preprint. Discussion started: 10 March 2020 c Author(s) 2020. CC BY 4.0 License. in the Antarctic domain, because ablation at these latitudes is limited Pollard and DeConto, 2012;Pattyn, 2017). Furthermore, in this particular study, the transient character of the AIS evolution is not simulated, as we focus on the LGM period. Thus, there is no need for explicitly accounting for the effects of changes in insolation on melting. Calving occurs when the ice-front thickness decreases below an imposed threshold (200 m in this study) and the upstream ice flux is not large enough to provide the necessary ice for maintaining the previous thickness (Peyaud et al., 2007). Present-day basal 5 melting rates at the ice-shelf base and at the grounding line are obtained from Rignot et al. (2013) and extrapolated over all 27 basins identified by Zwally et al. (2012). Below grounded ice, the basal mass balance is determined through the heat equation as in Greve and Blatter (2009), where the geothermal heat flux field is obtained from Shapiro and Ritzwoller (2004). The glacial isostatic adjustment (GIA) is computed with the elastic lithosphere-relaxed asthenosphere (ELRA) method (Le Meur and Huybrechts, 1996), where the relaxation time of the asthenosphere is set to 3000 years. 10 Yelmo does not explicitly model the impact of ice anisotropy on the ice flow, so the classical "enhancement factor" are used as a tuning parameter (Ma et al., 2010;Pollard and DeConto, 2012;Maris et al., 2014;Albrecht et al., 2019). For this study we found realistic PD states for E grounded =1.0 and for ice shelves E floating =0.7.

Basal-drag law
As mentioned above basal sliding is calculated within the SSA solution, which is a function of the basal stress. Yelmo computes 15 the basal stress at the ice base (τ b ) through a linear viscous friction law. It depends on the basal ice velocity (u b ), the effective ice pressure (N eff ) and a tunable friction coefficient (c b ): and , is a coefficient that reflects the bedrock characteristics, and N eff is the effective ice pressure, given in [kPa]. Here we have parameterized c b as a function of the bedrock elevation, z b (positive above sea level), analogous to previous work (e.g., Martin et al. (2011)): Here, z 0 is an internal parameter that determines the bedrock e-folding depth over which the friction coefficient c b decreases 25 from a maximum value of c max reached for bedrock elevations above sea level (z b ≥ 0) and a minimum threshold value c min .
For higher values of z 0 (i.e., lower absolute values of z 0 ), c b falls more rapidly with depth. This parameterisation captures the phenomenon by which the occurrence of sliding (and its intensity) is favoured at low bedrock elevations and specifically within the marine sectors of ice sheets. It follows a similar approach as in Albrecht et al. (2019) and Martin et al. (2011), where the bedrock friction (in their case the "till friction angle") depends on the bedrock elevation. The effective pressure is represented by the Leguy et al. (2014) formulation, under the assumption that the subglacial drainage system is hydrologically well connected to the ocean so that there is full support from the ocean wherever the icesheet base is below sea level. We thus assume that the exerted basal pressure at the land-ice interface depends on the difference between the overburden pressure and the basal water pressure (i.e. the distance from flotation as measured in ice thickness), hence: where ρ i is the density of ice, g is gravity, H is the ice thickness and H f is the flotation thickness, given by H f = max 0, − ρw ρi z b , where ρ w is the seawater density, respectively, and z b is the bedrock elevation (positive above sea-level). In this way, far from the grounding line, H f = 0 and N eff = ρ i gH, while at the grounding line, where H = H f , N eff = 0. This ensures continuity of τ b at the grounding line.

Climate forcing
To simulate the AIS at the LGM, Yelmo is run over 80 kyr with constant LGM conditions. The atmospheric forcing field is given by the following equation: where T atm 0 is the PD temperature field at sea level obtained from RACMO2.3 forced by the ERA-Interim reanalysis data (Van Wessem et al., 2014) and ∆T atm LGM-PD is the LGM surface temperature anomaly relative to the PD. The monthly-mean temperature fields are obtained from each of the the eleven PMIP3 models, as well as by the ensemble mean (Fig. 1a). We apply a lapse rate correction that accounts for changes in elevations (0.008 K m −1 for annual temperatures and 0.0065 K m −1 for summer temperatures).
The LGM precipitation is calculated as where P 0 is the PD monthly-mean precipitation obtained in the same way as the PD temperature and δP LGM/PD is the relative anomaly between the LGM and PD obtained from the PMIP3 ensemble. Figure 1b shows the resulting precipitation field, P LGM , for the PMIP3 ensemble mean. Precipitation is corrected with local temperature anomalies through Clausius-Clapeyron scaling, which assumes more accumulation for warmer temperatures and therefore lower elevations. Note that precipitation is 25 given in water equivalent and transformed into accumulation via changes in density (i.e. 1 m yr −1 water equivalent ca. 1.09 m ice). Basal-melting rates for floating ice shelves are set to zero in the LGM state.

Basal friction
To investigate the impact of changes in basal friction on the LGM AIS we assess the sensitivity to the friction in marine zones via the minimum friction allowed (c min ) and the elevation parameter (z 0 ) in Eq. 3 that controls how quickly friction decreases with depth. For this purpose we force Yelmo with a single reference climatic state obtained from the average anomaly of the 5 PMIP3 ensemble for the LGM climate (Fig. 1) and a range of friction parameters. This range was determined in two steps.  Figure S1).

Climatic fields
To understand the impact of changes in climatic forcing on the ice sheet, we fix the friction parameter values to a single, reference set of values (z 0 = -175 m and c min = 1·10 −5 yr m −1 ) and analyze the AIS simulated at the LGM for the climatic 15 forcing derived from each of the 11 models in the PMIP3 ensemble, using the aforementioned forcings for temperature (Eq. 5) and precipitation (Eq. 6). We focus on how the temperature and precipitation fields control the size and extent of the ice sheet.
In all experiments the sea-level change estimates are computed with respect to the simulated PD state for the reference friction parameter values.

Impact of basal friction
Here we present our LGM simulated AIS for different basal friction parameters. Ice volume is converted into a sea-level contribution by subtracting the floating portion and taking isostatic depression of the bedrock into account .  (Lambeck and Johnston, 1998;Lambeck and Chappell, 2001;Lambeck et al., 2002Lambeck et al., , 2003. Note that in order to avoid biases due to Yelmo's coarse spatial resolution, these extensions were computed using the ice-sheet 5 margins of each of the reconstructions at Yelmo's spatial resolution. The three lowest bound simulations correspond to cases for which the corresponding PD AIS ice volume deviates from PD observations by more than 3.5 msle (see SI, Fig. S1, S2).
The simulated surface velocity pattern shows a distribution with low values near the summit and increasing values towards the margins (Fig. 3). Our friction parameterisation reproduces the fact that ice streams become faster on topographic lows with the Amery, Wilkes and Victorias Land showing active ice streams of more than 50 m yr −1 (Fig. 3a,b). The WAIS, due to its 10 marine character, is also a very active sector. Ice volume differences between a slowly decreasing friction (z 0 =-200 m) and a more rapidly decreasing friction (z 0 =-150 m) primarily originate in the WAIS and the coastal marine regions of the EAIS and its surroundings (Fig. 3c), and are the result of higher basal velocities with lower friction values (Fig. 3d) leading to thinner ice.
Subtle differences are found when comparing the extension of grounded ice in our simulated AIS with previous reconstruc- Our simulated extension stands between the ICE-6G model (green line in Fig. 4) and the RAISED Consortium (red line) and 20 the ANU (blue line) model. The largest discrepancies between models occur on the Ross shelf. Whereas ANU and RAISED estimate an advance close to the continental-shelf break, ICE-6G is more retreated, while our results support a nearly complete advance.

Impact of climatic forcing
Here we present the simulated LGM AIS of each individual PMIP3 model for the reference friction parameters (Fig. 5). The 25 simulated ice-volume anomaly ranges from 7.8 msle to 14.0 msle (Fig. 6), a difference of 6.2 msle. The total ice extension ranges from 14.6 million km 2 to 15.8 million km 2 , a difference of 1.2 million km 2 . Thus, while the spread in ice volume is somewhat smaller than found when investigating the sensitivity to friction, the spread in extension is significantly larger.
Because the underlying dynamics in Yelmo are the same in all cases, the differences in size and extension can only be explained by differences in the climatic fields. To determine the causes underlying these differences, we investigate the sen-30 sitivity of the ice thickness and extension to the climatic fields used to force the ice-sheet model (Fig. 7). We find that higher accumulation results in a thicker ice sheet (Fig. 7a), but has no appreciable effect on the ice extension (Fig. 7b). For model climatologies for which the LGM ice sheet extends close to the continental-shelf break (an extension of around 15.5 million km 2 , see Fig 7d), the AIS ice volume increases with increasing accumulation (Fig. 7c). However, there are four climate models 8 https://doi.org/10.5194/tc-2020-28 Preprint. Discussion started: 10 March 2020 c Author(s) 2020. CC BY 4.0 License.
(CNRM-CM5, GISS-E2-R-150, GISS-E2-R-151, FGOALS-g2) that despite having higher accumulation on average than the ensemble mean, do not allow the ice sheet to advance as much as the other models, leading in all cases to extensions below 15 million km 2 (Fig. 7b). Therefore, the simulated AIS volume is smaller for these less advanced ice sheets, despite the relatively high accumulation rates imposed. For all the others, for which extension is around 15.5 million km 2 , the AIS ice volume clearly increases with increasing accumulation (Fig 7c).

5
Further inspection allows us to identify the atmospheric temperature close to the grounding line (Fig 7d) as a critical factor in determining how far the AIS advances. Whereas low temperatures present similar ice extension, as it becomes warmer the ice sheet is more retreated. Given the low temperature values, ablation can be generally discarded as the source of this behaviour (SI Fig. S3; there is, however, one exception, as discussed below), so we turn our attention to ice viscosity. A necessary condition for marine-based ice sheets to advance is that the ice thickness at the grounding line overcomes the flotation criterion 10 as sustained through accumulation and/or by inland ice flow. This condition is fulfilled when the ocean depth (z b ) is shallower than ∼90% of the ice thickness. Warmer ice temperatures lower the ice viscosity (Fig. 7e) and prevent the grounding-line to thicken, as a consequence of enhanced ice flow, and advance towards more depressed bedrock zones. Therefore, simulations with lower ice viscosity as GISS-E2-R-150, GISS-E2-R-151 and FGOALS-g2 do not fully advance in the Ross shelf, Pine Island or the Amery (Fig. 5,6). 15 Finally, CNRM-CM5 is a particular case which does not fulfil any of our proposed hypotheses. Viscosity describes the fluidity of a material, therefore warmer temperatures enhance ice flow. Thus, following the same reasoning as before, one would expect a low viscosity as a consequence of a warmer ice column for CNRM-CM5, which is not the case (Fig. 7e). This model expands fully at the Ross shelf and Antarctic Peninsula zone, but the Ronne shelf is far from the grounding-line and the Amery shelf is even more retreated than PD (Fig. 5). The ice sheet does not advance in these regions due to the presence of abnormal 20 ablation, which impedes the ice expansion (see SI, Fig. S3). We argue that the unexpected large viscosity is a consequence of two competing effects. The fully advanced regions, as the Ross basin, contribute to a rather low ice temperature and hence a high viscosity. On the other hand, the ablation zones such as the Ronne and Amery basin, have warmer ice temperatures which conclude into low viscosity. Therefore Fig. 7e shows that CNRM-CM5 has on average a warm ice column and a high viscosity.
A similar reasoning can be applied to Figure 7a where the mean ice thickness is low despite its high accumulation.

25
In summary, we find that the choice of the boundary climate is crucial for the simulated LGM ice sheet. On one hand, the atmospheric temperatures near the coastal regions control the ice extension through viscosity. If the viscosity is too low, then the ice flows too fast, preventing the necessary thickening. Particularly, if the bedrock is too deep, the ice sheet's expansion will be hampered. Secondly, if the ice sheet extends close to the continental-shelf break, then the accumulation pattern will determine the total amount of ice volume. We find that for similarly extended ice sheets (IPSL-CM5A-LR and MRI-CGCM3), 30 the sea-level difference due to accumulation differences is about 3.5 msle.

Spatially homogeneous approach
Applying a simple scheme that lowers the ice accumulation and surface temperature homogeneously over the whole domain is a common and valid approach at first order, because during the LGM, at continental scale, a colder and drier climate is 9 https://doi.org/10.5194/tc-2020-28 Preprint. Discussion started: 10 March 2020 c Author(s) 2020. CC BY 4.0 License.
expected (Huybrechts, 2002;Golledge et al., 2012;Whitehouse et al., 2012a;Gomez et al., 2013;Quiquet et al., 2018). We thus tested a spatially homogeneous scaling (hereafter, the homogeneous method) for comparison. All simulations simulated realistic SLE and ice extensions during the LGM for the same friction coefficients. In overall, consistently lower ice volumes as well as reduced ice extensions are simulated with the homogeneous method (Fig. S4). Again, because the ice dynamics are the same, this difference can only be explained by the climatic forcing. Moreover, because temperatures are not sufficiently high 5 to produce ablation it points to ice accumulation differences. Fig. 8c illustrates the ice thickness difference between the two methods for a similar ice extension (Fig. 8a,b). The anomaly shows that the main source of this difference in ice volume and extension comes from the WAIS. The Antarctic Peninsula in particular shows a high positive thickness anomaly for the average PMIP3 climatic fields relative to the homogeneous case because the grounding-line does not advance there in the latter case. In the EAIS, the anomalies are not so pronounced, however inland ice is slightly thinner, whereas closer to the coast it is thicker.

10
This anomaly pattern can be explained by the difference between the accumulation fields (Fig. 8d). The spatially homogeneous method accumulates more ice inland and and leads to a reduced accumulation towards the continental-shelf break, especially at the Ross shelf, Pine Island and the Antarctic Peninsula. Because ice cores are generally extracted from dome regions with colder conditions, it is expected that precipitation and air temperatures near the coast are underestimated by the homogeneous approach.

Basal dragging law
Even at present-day it is difficult to estimate bed properties like basal temperature or ice velocities, which could improve our understanding of basal friction. Therefore, estimating bed properties at the LGM, where the total ice volume and extension is not fully constrained, adds a degree of difficulty. We covered a range of friction values which lead to realistic LGM and PD 20 configurations. The simulated sea-level differences were about 7 msle between the extreme cases (Fig. 2). We found that the choice of different bedrock frictions has an impact on ice-stream activity in marine-based regions. For example, an AIS that extends up to the continental-shelf break, but with a relatively low volume increase, can be achieved through a very dynamically active ice sheet. In that case, marine regions, and more specifically the WAIS, have the potential to maintain fast ice streams at the LGM and still agree with PD observations.

25
The choice of a given and unique friction law for the whole AIS is still somewhat arbitrary and unconstrained. We focused on a linear viscous friction law commonly used in other studies Quiquet et al., 2018;Alvarez-Solas et al., 2019). We are aware that other types of friction laws could have been tested, such as a regularized Coulomb law (Joughin et al., 2019) or a Coulomb-plastic behaviour (Nowicki et al., 2013), typically for ice flowing over a bedrock filled with cavities.
However, the importance of saturated tills is specially determinant for transient simulations with a retreating grounding line.

30
Given the large uncertainty we quantified for only one friction formulation, we expect that this range would increase further considering additional formulations.

Sea-level and ice extent uncertainty
For our reference friction parameters we used the individual climate simulations of the participating PMIP3 groups as surface boundary forcing. The sea-level difference between the models was about 6.2 msle. The lowest sea-level contribution was 7.8 msle (CNRM-CM5) and the largest 14.0 msle (IPSL-CM5A-LR). These sea-level estimates were inside the range of other studies and reconstructions. From this point of view, we were not able to discard any specific model field. Nonetheless, it seems 5 unrealistic that air temperatures were high enough to produce ablation during the LGM as seen in CNRM-CM5.
The simulated ice extension is determined through air temperatures. Warmer temperatures lower the ice viscosity. Due to the marine character of the AIS, a lower viscosity enhances ice flow leading to thin ice in regions where the bedrock is too deep, which prevents a complete advance towards the continental-shelf break. Forcing from the models GISS-E2-R-150 and GISS-E2-R-151 for instance do not allow a full advance in the Ross shelf, resembling the ICE-6G reconstruction (Fig. 4). 10 Similarly, with FGOALS-g2 the advance into the Pine Island region or the Amery shelf advance is impeded (Fig. 5). On the other hand, if temperatures are sufficiently cold, less than -20ºC or so, then the ice fully advances as in the ANU reconstruction ( Fig. 4). The RAISED Consortium has a similar extension, but presents two large ice shelves at the margins of the Ronne shelf, which we are not able to simulate. Again, the simulated ice extensions were inside the range of the reconstructions, and we could not exclude any case. But we found that in addition to the precipitation field, temperature fields play a crucial role as 15 they have the potential to accelerate the ice by lowering the viscosity and determine the total grounded ice area, which in turn affects the grounded ice volume.
Of course there are several sources which could impact AIS volume estimates, aside from the climatology and basal friction.
A change in bedrock depth, for instance, has profound implications on the simulated AIS, as it does not only change the local sea level, but it can also facilitate (or impede) the ice advance and retreat (Philippon et al., 2006). Here we used a simple 20 parameterization that accounts for the elasticity of the lithosphere and a non-local response caused by lateral shift (Le Meur and Huybrechts, 1996). This formulation does not capture differences in the mantle viscosity as it applies the same spatially homogeneous time response. Nonetheless, the Antarctic bedrock is a complex component with different rheological properties.
The WAIS for instance is a low-viscosity region where the bedrock deformation happens on a shorter timescale (Whitehouse, 2018;Whitehouse et al., 2019). The next generation of ice-sheet models coupled to GIA models may produce more realistic 25 bedrock responses and hence help to improve the sea-level budget at the LGM. This can be helpful for instance to constrain the phase space of friction parameters.

Forcing methods
Overall, a homogeneous anomaly relative to present day simulates a lower ice volume as a consequence of low accumulation near the ice-sheet margins (Fig. 8b). This indicates that the AIS could have stored more ice at the LGM than estimated by 30 studies applying such a scheme. As opposed to a spatially homogeneous method, GCM outputs are capable of representing local atmospheric effects, such as atmospheric circulation changes or localized precipitation structures. Thus, the latest icesheet models have begun to be forced by more detailed and arguably more realistic climatic fields (Briggs et al., 2013;Maris et al., 2014;Sutter et al., 2019). Nevertheless we have shown here that the spread of the simulated ice volume and ice extension for different climatic outputs can be equal to or larger than that resulting from different basal-dragging choices. The PMIP3 LGM climatologies are built with a prescribed ice extension and surface elevation (Abe-Ouchi et al., 2015). It is clear then, that by construction, ice models should be driven towards these particular configurations. Nonetheless GCM models may exhibit biases in the temperatures and precipitation in localized regions. A way to potentially test the plausibility of the employed 5 climatic fields is to compare with ice proxies. We strongly recommend that paleo ice-sheet simulations should be performed with GCM outputs, as they capture more complex processes than a spatially homogeneous method, but the choice of the climatic fields has to be consistent with reconstructions. In the future with PMIP4 results, more accurate climatic fields are expected.

10
The ice dynamics and the boundary climatology are two essential building blocks for the simulation of an Antarctic LGM state. Here we studied the uncertainty in LGM ice volume associated with these two factors, by investigating the effect of the representation of basal friction and of the atmospheric forcing, respectively, in simulations. First, we tested a range of potential basal friction values of marine zones which simulated plausible LGM states. We found that for a simple linear friction law lower (larger) friction values enhance (diminish) the ice dynamics of marine zones and result in ice sheet configurations with 15 less (more) ice volume, but still similar grounded ice extension. This led to several potential configurations of the AIS with a sea-level difference with respect to today in the range of 11.2 msle and 17.5 msle and with a total ice extension in the range of 15 to 16 million km 2 . Then, for a particular friction configuration within the estimates of ice volume and extension, we studied the individual sea-level contribution from simulations driven by LGM climates provided by the eleven PMIP3 participating groups. We found ice volume anomalies ranging from 7.8-14.0 msle and extensions of 14.6 to 15.8 million km 2 .

20
Imposing the PMIP3 fields, whose climate simulations include dynamic adjustment to the LGM boundary conditions, translate into higher precipitation rates along the Antarctic coast, hence leading to a larger simulated ice volume compared to using a homogeneous anomaly method. The grounding-line advance is strongly determined by the atmospheric temperatures as well.
Higher temperatures enhance ice flow reducing the ice viscosity. Because of the marine character of the AIS, relatively high temperatures near the coast can prevent ice expansion. Thus, along with improved knowledge of basal conditions, constraining 25 broader possible climatic changes during the LGM is imperative to be able to reduce uncertainty in the AIS volume estimates for this time period.

Code and data availability
Yelmo is maintained as a git repository hosted at https://github.com/palma-ice/yelmo under the licence GPL-3.0. Model documentation can be found at https://palma-ice.github.io/yelmo-docs/. The results used in this paper are archived on Zenodo