2.1 Ice flow enhancement factors

Enhancement factors account for unresolved effects of grain size, fabric and impurities and have often been used as tuning parameters in ice sheet modeling entering the constitutive flow law (Eq. 1) that balances strain rates and stresses within the ice sheet. A value of 1.0 means “no enhancement”, but generally enhancement factors for the SSA tend to be smaller than 1.0 and larger than 1.0 for the SIA. Anisotropic ice-flow modeling suggests E_SSA values between 0.5 and 0.7 for ice shelves and between 0.6 and 1.0 for ice streams, while for SIA enhancement factors they should lie between 5.0 and 6.0 (Ma et al., 2010). Previous model ensembles that consider isotropic ice flow use values down to E_SSA=0.3 (for E_SIA=1.0) in Pollard and DeConto (PSU-ISM; 2012a) or up to E_SIA=9.0 (for E_SSA=0.8) in Maris et al. (ANICE model; 2014), both for 20 km resolution. PISM seems to favor enhancement factors closer to 1.0 (e.g. E_SSA=0.5–0.6 and E_SIA=1.2–1.5 in Golledge et al., 2015, for 10–20 km resolution).

Figure 3 shows the effects of ice flow enhancement factors on the simulated Antarctic Ice Sheet history. SIA enhancement generally produces thinner grounded ice. Compared to the reference simulations with E_SIA=2.0 the model simulates ice sheets for E_SIA=1.0 that are more than 5 m SLE thicker and ice sheets for E_SIA=5.0 that are about 6 m SLE thinner, both at glacial and interglacial states (dark and light blue). The effect of E_SSA variation is most predominant for ice sheet growth and for deglaciation, when ice flow across the grounding line influences its migration and stabilization (Schoof, 2007b). For E_SSA=1.0 (orange) we find earlier retreat and hence 1–2 m SLE thinner modern ice sheets than for E_SSA=0.6 (grey reference). For lower values of E_SSA=0.3 (light orange), in contrast, deglaciation is limited and modern ice volumes are more than 5 m SLE thicker than observations. Small values for the SSA enhancement factor produce slower and thicker ice streams and ice shelves. However, as the SIA enhancement factor, ESIA, has a larger influence on the ice sheet's volume on glacial scales, we find it a more suitable ice-internal parameter candidate for the ensemble analysis in Albrecht et al. (2020).

2.2 Ice rheological flow-law exponent

Generations of ice sheet modelers used variants of Glen's flow law (Eq. 1) as the constitutive relationship between stress and internal flow with the rheological exponent n=3. According to an analysis of data in Greenland by Bons et al. (2018) the ice rheological exponent for the SIA stress balance should instead be n=4. For the same strain rates at a reference pressure of 50 kPa, we would need to adjust the ice softness factor A accordingly.

We tested the effect of the flow-law exponent n=4 (and n=2 to cover the range according to Goldsby and Kohlstedt, 2001) in comparison to the reference simulation and found only small differences in the ice volume time series during glacials and moderate differences in interglacial periods, with on average less than 0.9 m SLE difference (Fig. 4). However, the flow-law exponent for the SSA has much stronger effects on ice volume, with 2–7 m SLE less Antarctic ice volume for n=4 and significantly earlier deglaciation, while for n=2 deglaciation is effectively damped.

2.3 Model grid resolution

Resolution is a key parameter which determines the misfit between the model results based on discretized model equations and the associated analytical solutions. Analytical solutions for the coupled system of ice sheet model equations, however, can only be found for simplified configurations. For the more realistic case we can get some impression of resolution requirements if we run regression tests to show that grid refinement leads to a convergence of the model solution (Cornford et al., 2016). The horizontal resolution of the boundary data, for instance, can control key parameters of ice-stream flow, such as basal roughness (Falcini et al., 2018). Here we also test for the vertical resolution of the three-dimensional enthalpy field. In order to deal with the changing geometry of the ice, especially in the case of a non-flat moving bed, the independent variable z relative to the geoid is replaced in PISM by the vertical coordinate relative to the base of the ice. The vertical resolution is highest at the ice sheet's base and lowest at the top of the computational domain using quadratic spacing. The enthalpy formulation allows the transition from cold to temperate ice (Aschwanden et al., 2012), which can form temperate ice layers of up to a few hundred meters.

We find that for our reference parameters and 16 km resolution a similar ice volume history can be reconstructed as for 8 km resolution (see Fig. 5), while computation costs are higher by about a factor of 10. For much coarser grids of the order of 30 km or more (Briggs et al., 2014, used 40 km) we find that relevant ice-stream dynamics cannot be resolved anymore in an adequate way (Aschwanden et al., 2016). In Fig. 5 we find how resolution can effect glacial–interglacial ice volume history, resulting in very different modern ice sheet configurations (see blue lines in Fig. 5). In fact, also for coarse resolutions we may find solutions that are closer to the reference simulation, e.g. by choosing different enhancement factors.

Figure 3Ice volume above flotation for simulations over the last two glacial cycles, with varied enhancement factors for SIA and SSA stress balance. Lower E_SIA values lead to thicker grounded ice sheets, while larger E_SSA values cause faster deglaciation and slightly thinner modern ice sheets.

Download

Figure 4Simulations with different flow-law exponents in the SIA and SSA stress balance, varying from 2 to 4. For variation in the SIA flow-law exponent we find only small impact on ice volume (blue, grey and light blue), while for the SSA it strongly affects ice sheet growth and deglaciation (orange and light orange).

Download

For the vertical grid, we define, in the reference simulation, the narrowest grid spacing at the base with around 20 m. Coarser vertical resolution (doubled spacing) does not change the simulation result much (orange line in Fig. 5). For finer resolution, in contrast, shear heating and the formation of temperate ice is expected to be better resolved, and the simulation results should converge. However, the simulated ice volume seems to increase by 3–5 m SLE for doubling vertical resolution (see light orange line in Fig. 5), as less temperate ice is formed in the lowest layers of the ice sheet. Benchmark experiments with respect to an analytical enthalpy solution (Kleiner et al., 2015) suggest adequate convergence for vertical resolution finer than 1 m at the base (see violet line in Fig. 5), although this comes with much higher computational costs and memory requirements. For the used set of model parameters, such a high vertical resolution yields a present-day state close to LGM, which is 14 m SLE thicker than observed. This finding emphasizes the fact that resolution is a relevant model parameter that should be taken into account in model tuning.

2.4 Iceberg calving

Currently, calving constitutes almost half of the observed mass loss from the Antarctic Ice Sheet (Depoorter et al., 2013). PISM provides different schemes for iceberg calving. For floating ice shelves we use the strain-rate-based eigencalving parameterization, which accounts for the average tabular iceberg calving flux, depending on ice shelf flow and confinement geometry (Levermann et al., 2012; Albrecht and Levermann, 2014). The minor eigenvalue of the horizontal strain rate basically determines where calving can occur, e.g. in the expansive flow regions beyond a critical arch between outer pinning points of ice rises or mountain ridges (Doake et al., 1998; Fürst et al., 2016). The average eigencalving rate is the product of the minor and major eigenvalue of the horizontal strain rate scaled with a constant K.

In our simulations we use a reference parameter value of $K=\mathrm{1}\times {\mathrm{10}}^{\mathrm{17}}$m s. For a much larger value of 1×10¹⁸ m s the calving front tends to retreat but is limited by the location of the compressive arch. As smaller ice shelves exert less buttressing on the ice flow, we hence find slower ice sheet growth for glacial climate conditions (see Fig. 6) but negligible effects on deglaciation or interglacial ice volume. For a smaller value of 1×10¹⁶ m s, in contrast, estimated calving rates tend to be smaller than the terminal ice shelf flow, and thus the calving front expands up to the edge of the continental shelf. The additional buttressing supports a slightly larger present-day ice volume, while in turn the more extended LGM ice shelves are more exposed to increasing sea level and ocean temperatures, leading to slightly earlier deglaciation (as already discussed in Kingslake et al., 2018).

In our simulations we define a maximum extent for ice shelves where the present ocean floor drops below 2 km depth, assuming that ice shelves can only exist on the shallow continental shelf (“deep-ocean-calv”). Additionally we avoid very thin ice shelves below 75 m, as enthalpy field evolution and hence the ice flow cannot adequately be represented for only a few vertical grid layers. Hence, ice at the calving front thinner than 75 m is removed. Both calving conditions are applied mainly for numerical reasons (adaptive time stepping) to avoid thin ice tongues, but they have negligible influence on the simulated ice volume history.

For higher lower bounds of terminal ice thickness of 150 m or even 225 m, as often used in other model studies, we find slower ice sheet growth but a negligible effect on deglaciation and interglacial ice volume (see Fig. 6). As eigencalving supports ice shelves within confinements we find that the effects of ice shelf extent beyond its compressive arch are of relatively low relevance to the glacial–interglacial ice sheet history. This is in contrast to previous model ensembles, which diagnosed high sensitivity of simulation results to varied calving parameters using different calving parameterizations (Briggs et al., 2013; Pollard et al., 2016).

Figure 5Ice volume above flotation for simulations with horizontal reference resolution of 16 km (grey), fine grid resolution of 8 km (light blue) and relatively coarse resolution of 32 km (blue), resulting in different Antarctic Ice Sheet histories. Refined vertical resolution of 10 m (light orange) and 1 m (violet) at the base yields larger ice volumes than the reference simulation with 20 m resolution (grey) or for coarse resolution of 40 m (orange).

Download

Figure 6Ice volume histories over two glacial cycles for different parameter values of eigencalving constant (blue and light blue) and calving thickness (orange and light orange) compared to the reference simulation (grey). Defining a maximum extent of the ice along the edge of the continental shelf (“deep-ocean-calv”) has only negligible effect on the sea-level-relevant ice volume (violet).

Download

2.5 Mantle viscosity and flexural rigidity

The topography of the bed has a decisive role in the stability of the marine ice, and it can change significantly on paleo-timescales. PISM incorporates an Earth model that reflects the deformation of an elastic plate overlying a viscous half-space, based on Lingle and Clark (1985). A key advantage of this approach over traditional elastic-lithosphere–relaxing-asthenosphere (ELRA) models is that the response time of the bed topography is not considered a constant but depends on the wavelength of the load perturbation for a given asthenosphere viscosity (Bueler et al., 2007). Calculations are carried out using the computationally efficient fast Fourier transform to solve the biharmonic differential equation for vertical bed displacement in response to (ice) load changes σ_zz (Bueler et al., 2007, Eq. 1). The Earth model can be initialized with a present-day uplift map (Whitehouse et al., 2012) and reproduces plausible uplift pattern and magnitudes for a given load history (Kingslake et al., 2018). However, it is still a simplification of the approach used within many glacial isostatic adjustment (GIA) models (Whitehouse, 2018), which are defined to account for the response of the solid Earth and the global gravity field to changes in the ice and water distribution on the Earth's surface (Whitehouse et al., 2019). With our modification of PISM's Earth model approach we can also account for vertical bed displacement in response to spatially varying water-load changes. However, the model is not able to account for self-consistent gravitational effects associated with local sea-level variations or the rotational state of the Earth (Gomez et al., 2013; Pollard et al., 2017).

We investigate two relevant parameters of the Earth model with regard to both the viscous and the elastic part.

Mantle viscosity affects the model behavior because it defines the rate and pattern of the deformation of the ice sheet's bed and the sea floor. GIA modeling indicates a range of viscosities of 10²⁰–10²¹ Pa s for the upper mantle, while the lower mantle is less constrained, with 10²¹–10²³ Pa s (Whitehouse, 2018). Our PISM simulations distinguish neither between the upper and lower mantle nor between the East or West Antarctic Ice Sheet (EAIS and WAIS, respectively). We choose a reference value of 5×10²⁰ Pa s, although for the relatively weak mantle beneath the WAIS even lower values are likely (Gomez et al., 2015), and for some regions mantle viscosity values <10²⁰ Pa s have been suggested (Hay et al., 2017; Barletta et al., 2018).

In our glacial-cycle simulations a lower mantle viscosity of 1×10²⁰ Pa s results in delayed deglaciation and thicker present-day ice sheets than for the reference value of 5×10²⁰ Pa s (cf. blue and grey in Fig. 7a). As the bed at the grounding line responds faster to unloading for lower viscosities, grounding-line retreat is hampered accordingly. In contrast, more viscous mantles of 25–100×10²⁰ Pa s result in slower ice sheet growth but in faster deglaciation and hence lower present-day ice volumes above flotation (light blue and orange in Fig. 7a). Within the range of 10²⁰–10²¹ Pa s the effect of mantle viscosity on grounding-line retreat and re-advance since last deglaciation has been discussed in a PISM study by Kingslake et al. (2018). In Appendix A2, we highlight the influence of further model decisions on grounding-line sensitivity and thus on the onset of deglacial retreat, although this is not the focus of this study.

The bed deformation model in PISM¹ up to version 1.0² considered all changes in ice thickness H to be loads, including changes in ice shelf thickness, although this does not make physical sense. Here we present simulations that consider changes in the load of the grounded ice sheet and of the ocean layer within the computational domain, with load per unit area defined as

with h_f being the flotation height as defined in Eq. (7). Be aware that since PISM v1.1 the bed deformation model has been fixed, and only the grounded ice sheet changes are considered to be a load.

Our simulations that additionally consider changes in ocean loads yield up to 3 m SLE higher ice volumes above flotation at glacial maximum and delayed deglaciation by a few thousand years, while interglacial ice volumes are comparable (cf. grey and light orange lines in Fig. 7a). The effect of changing ocean loads seems to be larger than the effect of accidentally added ice shelf loads (cf. violet line in Fig. 7b, which excludes both ocean and ice shelf loads as in PISM v1.1.

Flexural rigidity is associated with the thickness of the elastic lithosphere and has an influence on the horizontal extent to which bed deformation responds to changes in load. We have deactivated the elastic part of the Earth model in our reference simulation, as the numerical implementation was flawed. In order to evaluate the ice sheet volume's sensitivity to changes in the flexural rigidity parameter value, we used PISM v1.1 instead, with an additional fixed elastic part (https://github.com/pism/pism/pull/435, last access: 27 January 2020). However, PISM v1.1 considers only grounded ice changes to be loads and not changes in the ocean layer thickness, as in the reference. As most dynamic changes on glacial cycles occur in West Antarctica, previous studies based on gravity modeling suggest appropriate values lying within the range of 5×10²³ N m and 5×10²⁴ N m (Chen et al., 2018). The PISM default value marks the upper end of this range, assuming a thickness of 88 km for the elastic plate lithosphere (Bueler et al., 2007).

Extending this range to 1×10²⁵ N m, which is more than an order of magnitude, we find differences in ice volume above flotation of up to 4 m SLE, part of which might be due to increased temporal variability (see Fig. 7b). Compared to the reference simulation without the elastic part, we find earlier deglacial retreat but similar present-day ice volumes.

Figure 7Simulations over two last glacial cycles, with varying Earth model parameters for the viscous and the elastic part. (a) Mantle viscosities with a range over 2 orders of magnitude cause slower ice sheet growth but faster decay for increasing viscosity. (b) The tested flexural rigidity over a range of more than 1 order of magnitude yields smaller difference in simulated ice volume response but increased variability. Be aware that for (b) a different PISM v1.1 was used, with a fixed elastic model, but only accounting for changes in grounded ice loads, explaining generally lower glacial volumes than in the reference in (a).

Download

3.1 Air temperature

Air temperature is an important surface boundary condition for the enthalpy evolution which is thermodynamically coupled to the ice flow. The annual and summer mean temperature distribution are required in the positive-degree-day (PDD) scheme to estimate surface melt and runoff rates, assuming a sinusoidal yearly temperature cycle (Huybrechts and de Wolde, 1999, Eqs. C1–3). Hitherto, surface temperatures for Antarctica have been often parameterized based on a multiple regression fit to reanalysis data, e.g., as a function of latitude “lat” and surface elevation h. This provides a temperature field that adjusts to a changing geometry with a prescribed lapse rate and is hence convenient for paleo-timescale simulations. Using ERA-Interim C20 data (Dee et al., 2011) a multiple regression fit of summer mean temperatures (January) provides a temperature distribution with a RMSE (root-mean-square error) of 2.2 K over the entire ice sheet, while for annual mean temperatures the RMSE is 4.1 K. Temperatures are considerably overestimated by up to 11 K over the large Ronne–Filchner and Ross ice shelves and in large parts of inner East Antarctica, while temperatures in dynamically relevant regions along the WAIS Divide are underestimated by up to 5 K.

Regarding the comparably shallow ice shelves with less than 100 m surface height, this surface-height-dependent parameterization estimates temperatures close to those on the sea surface, although observed climatic conditions on the large ice shelves are much colder than those on the ocean, where sea ice can be absent in summer. Typically, the Antarctic Ice Sheet and shelves inhibit a negative radiation budget that efficiently cools the surface (van Wessem et al., 2014).

As a correction we assume that the ice shelves' surface (and all other icy regions below 1000 m altitude for consistency) is as cold as at 1000 m altitude to achieve a better fit to ERA-Interim surface temperature data. Furthermore, the climate of the much larger East Antarctic Ice Sheet is more isolated than the climate in West Antarctica, which can be accounted for with a longitudinal dependence “lon” and symmetry axis through the WAIS Divide (110^∘ W), as

Summer temperatures T_sml are well represented by Eq. (10), with a RMSE of 2.1 K over the entire ice sheet and a particularly good match in the large ice shelves and East Antarctica. Annual mean temperatures T_aml parameterized by Eq. (9) are overestimated by less than 5 K both in the inner East Antarctic Ice Sheet and the large ice shelves of Ross and Ronne–Filchner (see Fig. 8), while temperatures in the West Antarctic Ice Sheet are underestimated by less than 1 K (RMSE 3.1 K).

Figure 8Comparison of ERA-Interim annual mean temperatures (a) with multiple regression fit (b) according to Eq. (9), with 110^∘ W longitude indicated as blue dashed line, the 1000 m surface height isoline in green and the black dotted transects through the large ice shelf regions. Inner parts of the West Antarctic Ice Sheet are less than 1 K too cold (d), while temperatures in the shallow ice shelves of Ross and Ronne–Filchner are overestimated by up to 5 K (cf. c and f). (e) Root-mean-square errors of temperatures in grounded and floating ice sheet regions are 3.0 and 4.1 K, respectively.

Figure 9Comparison of ERA-Interim summer (January) mean temperatures with parameterization of Eq. (10). Root-mean-square errors of temperatures in grounded and floating ice sheet regions are 2.0 and 3.1 K, respectively, with temperatures in the large ice shelves close to observations.

The annual mean and summer air temperatures enter the PDD scheme to calculate melt rates and runoff at the surface. We assume melt rates of snow and ice of 3.0 and 8.8 mm for each day and per degree above freezing, assuming a daily temperature variability represented by a normally distributed white noise signal with a 5 K standard deviation.

In order to evaluate the transient effects of the choice of the modern climate boundary conditions we run two-glacial-cycle paleo-simulations with a climatic forcing that is introduced in the following sections of the paper. Here we compare the simulated histories of the ice sheet's volume above flotation with respect to different modern surface temperatures and PDD settings. A lower temperature standard deviation in the PDD scheme of 2 K instead of 5 K has no influence on glacial volume, but it can delay deglaciation slightly (see Fig. 10). Generally, the effect of the PDD scheme for Antarctic paleo-simulations seems of minor relevance. Even if we do not parameterize the surface mass balance via the PDD scheme and directly apply the annual mean temperature pattern from the parameterization, RACMO2.3p2 (van Wessem et al., 2018) or ERA-Interim (Dee et al., 2011), to the ice surface boundary, simulated ice volumes differ at most by up to 2 m SLE.

Figure 10Ice volume above flotation for variation in PDD standard deviation and for different modern mean temperature distributions. The sensitivity is low for different PDD parameters of 5 K (reference in grey) and 2 K and when no surface model with yearly cycle is used (blue and light blue respectively), also for interglacials. For more realistic annual mean surface temperatures from RACMO2.3p2 (van Wessem et al., 2018) or ERA-Interim (Dee et al., 2011), but without PDD and yearly cycle applied, the simulated volumes show most difference for glacial ice sheet growth and at interglacials (orange and light orange, respectively).

Download

3.2 Precipitation

The distribution of mean precipitation over the Antarctic continent is related to temperatures, but its pattern is strongly determined by the moisture transport over the ice and mountain surface such that parameterizations as for the temperature (see Sect. 3.1) are rather unrealistic. We use a more realistic precipitation field from the regional atmospheric climate model RACMO2.3p2 (van Wessem et al., 2018, as mean over CMIP5 reference period 1986–2005), which was forced with ERA-Interim reanalysis data for the recent past. RACMO2.3p2 incorporates all relevant physical processes at the ice surface simulated with a resolution of 27 km.

Analogous to the temperature parameterization, we apply a lapse correction to the RACMO precipitation field for changing surface elevation Δh to account for the elevation desertification effect (DeConto and Pollard, 2016), which we define as

with f_p=7 % K⁻¹ being a precipitation change factor according to the Clausius–Clapeyron relationship and with γ_T=7.9 K km⁻¹ being the temperature lapse rate. This correction ensures that topographical changes have an influence on local precipitation through their effect on local surface temperature. However, as reanalysis and regional climate models tend to underestimate present-day precipitation in the interior of the East Antarctic Ice Sheet, simulated ice volume may be biased towards being lower than observed values (Van de Berg et al., 2005).

3.3 Basal heat flux

Geothermal heat flux is one of the most poorly known boundary conditions that controls ice flow (Pittard et al., 2016). It can keep basal ice relatively warm and thus less viscous than the colder ice above. PISM uses a bedrock thermal model (1-D heat equation, similar to Ritz et al., 1996), with storage in an upper-lithosphere thickness of 2 km discretized in 20 equidistant layers. Geothermal heat flux is applied as a constant boundary condition to calculate the heat flux into the ice at the ice–bedrock interface depending on the ice base temperature. In combination with enhanced supply of meltwater at the ice sheet base, it supports rapid ice flow by sliding over the bed and deformation of the subglacial sediments (see Sect. 3.4.2). Various maps with substantially different patterns derived from satellite magnetic and seismological data have been made available for the whole Antarctic continent (https://blogs.egu.eu/divisions/cr/2018/03/23/image-of-the-week-geothermal-heat-flux-in-antarctica, last access: 27 January 2020) and have been used in ice sheet model simulations for longer then a decade now (Shapiro and Ritzwoller, 2004; Fox Maule et al., 2005; Purucker, 2013; An et al., 2015). Due to their sparse data coverage and significant methodological uncertainty, some modelers decided to use a simplified two-valued heat-flux pattern that distinguishes the geology of the East and West Antarctic plates (Pollard and DeConto, 2012a). In our reference simulation, we use the latest high-resolution heat-flux map by Martos et al. (2017), which is derived from spectral analysis of the most advanced continental compilation of airborne magnetic data.

For the different Antarctic basal heat-flux datasets we compare PISM-simulated quasi-equilibria after 50 kyr with constant boundary conditions and find only few differences, with 40–45 m RMSE between the resulting ice thickness distributions (Fig. 11), except for the lowest estimate of basal heat flux by Purucker (2013), with about 100 m RMSE.

Figure 11PISM present-day equilibrium results after 50 kyr for different basal heat-flux fields (upper row) affecting the near-ground ice temperatures and hence the ice thickness, shown here for the first equilibrium result as anomaly to Bedmap2 (first panel in lower row) or relative to the first result. First three columns show the results for the magnetic reconstructions by Martos et al. (2017) (as in the reference), Purucker (2013) and Fox Maule et al. (2005). Columns 4 to 5 show the seismic reconstructions by An et al. (2015) and Shapiro and Ritzwoller (2004), and last column shows, for comparison, a two-valued field as used in Pollard and DeConto (2012a), separating East and West Antarctica. The overall effect of the choice of geothermal heat flux on present-day equilibrium ice thickness measures up to 100 m RMSE for Purucker (2013), while for the other tested set-ups it differs only in individual spots, with a RMSE of 40–45 m.

Here we compare the simulated histories of the ice sheet's volume above flotation with respect to the different basal heat-flux boundary conditions. The overall effect on ice volume history seems rather small, with a variation of about 3 m SLE for glacial climates. The lowest mean heat flux of 54 mW m⁻² (Purucker, 2013) yields generally larger ice volumes and vise versa for the highest mean heat flux of 67 mW m⁻² (see Fig. 12; Martos et al., 2017). The transient sea-level-relevant volume shows most variance during interglacials, when in some relevant regions of West Antarctica, such as the Siple Coast, grounding-line migration is delayed, and the local anomaly in ice thickness reaches up to 500 m. In particular for the present-day state the largest range of simulated ice volumes is found with about 8 m SLE. Compared to other uncertainties discussed in this study, e.g. with regard to friction-related parameters (see Sect. 3.4), we evaluate the choice of the basal heat-flux distribution of low relevance to the total ice volume history. However, in another PISM ensemble study that covers the last 2 million years with a focus on ice domes and deep ice cores during interglacials, geothermal heat flux becomes the most relevant uncertain boundary condition (Sutter et al., 2019). Also, the thickness of the bed thermal unit seems to be of relatively low relevance, even when compared to the extreme case of basal heat flux directly applied to the ice sheet base, without a thermal unit (see Fig. 12b).

3.4 Basal friction

Subsurface boundary conditions and their influence on basal sliding are key to understanding Antarctic ice flow, in particular subglacial topography, basal morphology (e.g., presence of sediments) and subglacial hydrology (Siegert et al., 2018). In the following section we will investigate the model's sensitivity to the pseudo-plastic exponent of the basal sliding relationship, the till friction angle under grounded ice and underneath the modern ice shelves, and parameters related to subglacial hydrology, namely the till water decay rate and the critical fraction of overburden pressure.

3.4.1 Pseudo-plastic exponent

In PISM, basal sliding is parameterized according to Eq. (4) with sliding exponent q. A value of q=0.0 represents purely plastic (Coulomb) deformation of the till where ice flows over a rigid bed with filled cavities (Schoof, 2005, 2006). Many studies use $q=\mathrm{1}/\mathrm{3}$ (e.g., Schoof, 2007a; Pattyn et al., 2013; Gillet-Chaulet et al., 2016) or a linear sliding relationship between basal velocity and basal shear stress for q=1.0, as most commonly adopted for inversion methods (e.g., Larour et al., 2012; Gladstone et al., 2014; Yu et al., 2017). Note that u_b results from solving the non-local SSA stress balance (Bueler and Brown, 2009, Eq. 17) in which τ_b appears as one of the terms that balances the driving stress. The PISM default sliding velocity threshold is u₀=100 m yr⁻¹, where basal shear stress is independent of q.

Figure 12(a) Simulation results for different basal heat-flux distributions (Martos et al., 2017; Purucker, 2013; Fox Maule et al., 2005; An et al., 2015; Shapiro and Ritzwoller, 2004; Pollard and DeConto, 2012a). The range of the resulting ice volumes above flotation is rather small, with less than 3 m SLE on average. In interglacial periods the divergence in ice volume can reach more than 8 m SLE. (b) Sensitivity of transient ice volume to varied thickness of the bed thermal layer (model option “Lbz” in km, with “Mbz” the number of equidistant vertical grid subdivision of the layer; see legend). Reference thickness is 2 km (grey), which yields similar results to a 5 km thick bed thermal layer (light blue). Not taking any bed thermal unit into account leads to slightly lower ice volumes, in particular during LIG (blue).

Over the valid range of q we find a spread of reconstructed ice volumes of up to 12 m SLE (see Fig. 13) in transient simulations. Larger values lead to thinner ice sheets, in particular at interglacial states. Due to its large impact on ice volume the pseudo-plastic exponent q, or PPQ, is a key parameter in the ensemble analysis presented in a companion paper (Albrecht et al., 2020).

Figure 13Antarctic ice volume history over the last two glacial cycles for different values of the basal sliding exponent q. Between plastic (blue) and linear sliding (light orange), the ice volume above flotation can vary by 6 m SLE on average and by up to 12 m SLE in interglacial periods. Generally, smaller values of q lead to slower ice flow and hence to larger ice volumes.

Download

3.4.2 Till properties

The till friction angle ϕ is a shear strength parameter for the till in Eq. (5), associated with the geology of the bed. It can be parameterized in PISM as a piecewise-linear function of bed elevation (Martin et al., 2011), assuming that marine basins and ice-stream fjords have a rather loose sediment material, while being denser in the rocky regions above the sea level. The till friction angle is weakly confined; the observational range in ice streams is 18–40^∘ in Cuffey and Paterson (2010). For now, we assume that ϕ_min=2^∘ in the marine basins lower than −500 m, that ϕ_max=45^∘ above 500 m and that there is a linear gradient between those two levels.

We run 50 kyr equilibrium simulations, first with a spatially constant value of the till friction angle with ϕ=30^∘ (Cuffey and Paterson, 2010, see Fig. 14). The simulation shows a general overestimation of ice thickness, with anomalies of more than 800 m in West Antarctica and an overall RMSE of 372 m for the constant till friction angle. When using the piecewise-linear parameterization, the anomalies are negative in large parts of the East Antarctic Ice Sheet and at the WAIS Divide, with a lower RMSE of about 296 m. This suggests that estimates of till friction angle in parts of the submarine basins are too low, while along the Siple Coast and in the Transantarctic Mountains values seem too high. Nevertheless, the dependence on fjord topography supports narrowly confined structures of simulated ice streams, similar to those observed in many Antarctic regions.

In the following we want to demonstrate how till friction angle estimates can be optimized to fit the observed grounded surface elevation (or ice thickness) from Bedmap2 (Fretwell et al., 2013). We here follow the simple inversion method by Pollard and DeConto (2012b), but we invert for the till friction angle ϕ rather than for the basal sliding coefficient τ_c. In PISM the effective pressure N_till in Eq. (5) is physically determined by the subglacial hydrology model, while Pollard and DeConto (2012b) use basal temperature as a surrogate. As we run the simulation forward in time for constant boundary conditions the till friction angle is adjusted in every grounded grid cell every 500 years in incremental steps of Δϕ. These steps are proportional to the misfit to observed surface elevation (divided by 200 m) and bounded by $-\mathrm{0.5}{}^{\circ }<=\mathrm{\Delta }\mathit{\varphi }<={\mathrm{1}}^{\circ }$. Since the surface elevation is underestimated in the inner parts of the East Antarctic Ice Sheet by a couple hundred meters (due to underestimated precipitation forcing), the retrieved till friction angles reach a maximum value, set here as ${\mathit{\varphi }}_{i_{\max }}=\mathrm{70}$^∘, which enhances yield stress by about an order of magnitude (see middle row in Fig. 14). In contrast, on the Siple Coast the minimum values of ${\mathit{\varphi }}_{i_{\min }}=\mathrm{2}$^∘ compensate for the overestimated ice thickness (see bottom row in Fig. 14). Thus, the RMSE of ice thickness can be significantly reduced to 141 m (or even 123 m for ${\mathit{\varphi }}_{i_{\min }}=\mathrm{0.5}$^∘), and the modeled ice volume is only 0.5 % below observation. The retrieved distributions of till friction angles are rather independent of initial conditions and iteration parameters (not shown here). In fact, this method may overcompensate for inconsistent model boundary conditions or processes that are not accurately represented.

Figure 14PISM equilibrium simulations for different basal friction fields (τ_c in middle row) based on different till friction angle distributions (ϕ in first row). First column is with spatially uniform till friction angle of 30^∘, second column shows till friction angle according to PISM parameterization as function of bed topography and last columns use a simple inversion technique with different minimal ${\mathit{\varphi }}_{i_{\min }}$. Last row shows the resulting ice thickness anomaly to Bedmap2 observations (Fretwell et al., 2013), with improvements particularly in East Antarctica and Siple Coast region and overall RMSE reduced to 123 m (with parameters ESIA = 1 and C_d=5 mm yr⁻¹).

For the transient response of the ice sheet's volume to the different distributions of till friction angle, we find similar glacial ice volumes for the depth-dependent parameterization and the optimization, while in contrast interglacial ice volumes differ considerably (cf. light blue and grey in Fig. 15a). For the spatially constant till friction angle of 30^∘ underneath the grounded ice sheet and 2^∘ on the ocean floor, we simulate glacial-cycle volume histories that are a few meters larger than for the piecewise-linear parameterization (blue). In fact, we find a much larger effect on the ice sheet volume for variations in the till friction angle on the ocean floor.

Till properties under modern ice shelves

Friction (and also bed topography) at the ocean bed underneath the modern ice shelves is poorly constrained, as the optimization algorithms only apply to modern grounded regions (e.g., Pollard and DeConto, 2012b; Morlighem et al., 2017). The friction coefficient on the continental shelf has thus been chosen as one of the ensemble parameters in Pollard et al. (2016, 2017). In PISM the till friction angle accounts for the flow properties of the substrate and enters the yield stress definition as tan (ϕ) (see Eq. 5). As sandy sediments are prevalent in the ice shelf basins, low values of ϕ are likely in these regions (Halberstadt et al., 2018). Additionally, the till friction angle in the ice shelf basins is a crucial parameter which determines the thickness of the extended ice sheet for LGM conditions and hence the potential contribution to the global sea-level change.

Modeled LGM ice volumes increase by up to 2–4 m SLE per degree of change in the minimal till friction angle (see Fig. 15b). Compared to observations, we detect much higher volumes above flotation at present for ϕ_min≥3^∘. At the same time, relative volume changes between the LGM and modeled modern state become slightly smaller for rougher basins. This effect may be related to the effect of friction on the rate of grounding-line retreat. The till friction angle is an important uncertain parameter for possible WAIS collapse. As no (partial) WAIS collapse is induced in the simulations, we find very similar ice volumes for the Last Interglacial (LIG) and present-day periods. The spread of ice volumes among the four experiments with ϕ_min=1–5^∘ is on average 13 m SLE. We choose, as a reference, a lower bound for the till friction angle of ϕ_min=2^∘ in ocean regions, as simulated deglaciation shows a good match to modern ice volume (Fig. 15b; grey).

Figure 15Ice volume histories over two glacial cycles for different parameterizations of till friction angle. (a) For a spatially uniform ϕ (with 30^∘ in grounded ice and 2^∘ elsewhere; blue) and a till friction angle that is parameterized as function of the bed topography (light blue), the simulated ice volumes are generally larger than for the reference (grey), where ϕ is a result of the simple inversion technique. In those three cases the same minimal ϕ_min was used. For a smaller ϕ_min=1^∘ we find generally smaller ice volumes. (b) The choice of ϕ_min has strong influence on the reconstructed ice volume, with high till friction angles leading to more friction and hence thicker ice sheets. Within the plausible range from 1.0–5.0^∘ we find up to 18 m SLE difference in simulated ice volume above flotation.

Download

Till water distribution

The effective till water content in PISM's non-conserving hydrology scheme is a result of the balance between basal melting and a constant drainage rate (see Eq. 3). PISM's default drainage rate of 1 mm yr⁻¹ is smaller than the basal melting in most of the grounded Antarctic Ice Sheet regions such that till saturates over time. Higher decay rates can effectively drain the till water in the inner ice sheet regions, which generally cause less-extended and more-confined ice streams, less ice discharge, and hence thicker ice sheets.

In transient glacial-cycle simulations, this relationship applies for both present-day climate conditions (see Fig. 16) and colder-than-present climates. A till water drainage rate of 10 mm yr⁻¹ can cause up to 11 m SLE additional ice volume (light blue line in Fig. 17).

Figure 16Present-day result of glacial-cycle simulation, showing ice thickness anomaly in Bedmap2 (lower left) or in reference for different till water decay rates (from left to right: 1, 5 and 10 mm yr⁻¹), causing different till water distributions underneath the ice sheet (upper row). For a reference decay rate of 1 mm yr⁻¹ about 90 % of ice sheet's bed is saturated, while for 10 mm yr⁻¹ saturated till is only found in the coastal regions and underneath the fast-flowing ice stream.

Another relevant aspect is the initial till water fraction on ocean beds that become grounded. PISM assumes that the grounding line advances into dry till area W_till=0, where a till water layer can form over the following decades or centuries. If we assume a rigid till layer at the ocean floor instead, with $W_{\mathrm{till}}=W_{\mathrm{till}}^{\max }$, this affects grounding-line migration. We hence find slower growth of glacial ice sheet volume and much earlier deglaciation, while ice volumes are comparable to the reference case for present-day climate conditions (compare orange and grey line in Fig. 17).

Figure 17Simulation over two glacial cycles for different till water decay rates and for different fractions of the overburden pressure underneath the grounded ice sheet. (a) Higher decay rates than the reference value of C_d=1 mm yr⁻¹ (grey) cause less-extended ice streams, less ice discharge and hence thicker ice sheets (blue and light blue). When assuming maximum till water fraction across ocean beds, grounding-line advance of marine glaciers is decelerated. Accordingly, we find a less-extended and thinner ice sheet at glacial periods, earlier retreat, but similar present-day results (compare grey and orange lines). (b) The critical fraction of the effective overburden pressure δ has strong influence on the reconstructed ice volume, particularly at glacial maximum, with high parameter values leading to more friction and hence thicker ice sheets. For the evaluated range of 2 %–8 % we find up to 12 m SLE difference in the ice volume above flotation.

Download

4.1 Sea-level forcing

The Antarctic Ice Sheet, particularly in its western part, rests on a bed below sea level with floating ice shelves attached. The location of the grounding line in PISM is solely determined by the flotation criterion (cf. H≤h_f in Eq. 7) and therefore also by the current sea-level height z_sl. Marine ice sheet dynamics are hence sensitive to changes in sea level, which has been 120–140 m lower than today at the Last Glacial Maximum. We neglect regional sea-level effects due to changes of the rotation of the Earth or due to self-gravitation, which can have locally a stabilizing effect on the ice sheet (Konrad et al., 2015). Instead, we consider global mean sea-surface-height reconstructions prescribed by the GIA model ICE-6G_C, which accounts for a changing surface area of the oceans (Stuhne and Peltier, 2015, 2017, courtesy Dick Peltier). To analyze the sensitivity of the model's response to the choice of the sea-level forcing, we compare our results to other reconstructions by Lambeck et al. (2014), Bintanja and Van de Wal (2008), Imbrie and McIntyre (2006), and Spratt and Lisiecki (2016). Focusing on the last deglacial period, the timing of sea-level rise onset in response to the melting of the northern hemispheric land ice masses varies by a couple of thousand years among the different sea-level curves (Fig. 18). For instance, the reconstruction by Stuhne and Peltier (2015) peaks already before 25 kyr BP, with around 130 m below present (grey), while the much smoother SPECMAP sea-level curve by Imbrie and McIntyre (2006) has a minimum around 18 kyr BP and a comparably late relaxation to the present-day sea level (blue).

The modeled Antarctic Ice Sheet response at the Last Glacial Maximum is rather unaffected by the choice of the sea-level forcing, while in contrast the ice volume at present varies by up to 2.5 m SLE (Fig. 18d). In particular, the meltwater pulse 1a (MWP1a; Liu et al., 2016) at around 14.35 kyr BP, with a global sea-level rise of 9–15 m or more within a few hundred years (see grey vertical band in Fig. 18), is well represented as a step in the sea-level curve in the reference forcing time series of the ICE-6G_C model. About 4 kyr after MWP1a, this triggers a comparably early and quick grounding-line retreat in the Ross and Weddell Sea embayment, where the large modern ice shelves begin to float.

If self-gravitational effects were accounted for within a sea-level model, the local sea-level anomaly at the grounding line would be reduced compared with the global mean. A scaling of the sea-level forcing by 80 %–90 % would mimic the first-order feedback of self-gravitation on grounding-line motion. Interestingly, neither the ice volume at LGM nor that at present is significantly affected, while the onset of deglaciation occurs earlier and with a lower rate (Fig. 18; violet line).

4.2 Surface temperature forcing

Varying surface temperatures drive ice flow changes on glacial timescales. In PISM we model the nonlinear thermo-coupling via a Glen-type flow law with a generalized form of the Arrhenius equation. More specifically, PISM's default flow law is the polythermal Glen–Paterson–Budd–Lliboutry–Duval law (Lliboutry and Duval, 1985; Aschwanden et al., 2012), where the ice softness depends on both the temperature and the liquid water fraction. Hence it parameterizes the (observed) softening of ice at pressure-melting temperature as its liquid fraction increases. Since vertical diffusion processes in the ice are rather slow, the corresponding response of the ice sheet to surface temperature anomalies occurs with a considerable delay of up to a few thousand years, involving a long-lasting memory for past events (cf. Sect. 1.3). On paleo-timescales PISM can use temperature reconstructions from ice cores based on deuterium isotopes, such as the EPICA Dome C ice core (EDC; Jouzel et al., 2007). This is a well-resolved time series over the last 803 kyr (EDC3 age scale), representative of the inner East Antarctic Ice Sheet. However, most of the ice-dynamical changes on glacial timescales occur in the marine regions of the West Antarctic Ice Sheet. The much closer WAIS Divide ice core (WDC) provides a highly resolved temperature reconstruction (Cuffey et al., 2016), but it spans only the last 67.7 kyr. As surface temperature forcing ΔT(t) we consider the WDC temperature anomaly with respect to the year 1950 CE added to the parameterized temperature field of Sect. 3.1. For model periods before 67.7 kyr the temperature anomaly of EDC applies (with respect to the last 1000-year average), with a small jump of 0.2 K, which is within the variability (see Fig. 19a and b). The EDC temperature reconstruction shows generally higher variability than the WDC reconstruction.

As expected, the Antarctic Ice Sheet volume responds in our simulations with a delay of several thousand years to the surface temperature forcing. The LGM minimum in surface temperature reconstructions happened around 22 kyr BP in the WDC data, while the largest ice volume is simulated at 14 kyr BP. The main temperature rise at WDC occurred between 18 and 12 kyr BP, contributing to initial deglaciation at around 12 kyr BP, with major retreat between 8 and 4 kyr BP in our reference simulation. At the EDC location the reconstructed temperature rise happened about 1 kyr later and with more variability, leading to a more gradual deglaciation between 12 and 2 kyr BP (Fig. 19; blue). Comparisons with other ice core temperature reconstructions, however, suggest a superimposed lapse-rate effect due to surface height change during deglaciation at the WDC location (Werner et al., 2018). This means that Antarctic temperature anomalies with up to −10 K at glacial maxima may be systematically overestimated (personal communication with Eric Steig). We thus also tested for a scaled temperature reconstruction from WDC with a LGM temperature of only 6 K below present (see light blue lines in Fig. 19). Interestingly, this weaker temperature forcing results in slightly thicker glacial ice volume (probably an effect of temperature-coupled surface mass balance; see Sect. 4.4) and delayed deglaciation.

We also tested for the influence of temperature variability on the simulated ice volume (Mikkelsen et al., 2018) and found slightly earlier initial retreat for the 500-year moving-average WDC temperature forcing. In contrast, for added white noise of 1.5 K variance, the present-day ice volume is more than 5 m SLE larger than in the reference (compare orange, light orange and grey lines in Fig. 19).

Figure 19Time series of PISM-simulated ice volume above flotation relative to observations (Fretwell et al., 2013) over last 135 kyr (c) and last 31 kyr (d) forced with different surface temperature reconstructions (a, b) at WDC (Cuffey et al., 2016, grey) and EDC (Jouzel et al., 2007, blue), leading to different ice volume histories, particularly during deglaciation period. The WDC temperature reconstruction scaled by 60 % causes 2 m SLE larger glacial ice volume and slower deglaciation (light blue). When WDC temperature reconstruction is smoothed, we find a slightly earlier initial retreat (orange), while added white noise leads to higher ice volume at present (light orange).

Download

4.3 Ocean temperature forcing

Sub-shelf melting in PISM is calculated via PICO (Reese et al., 2018) from salinity and temperature in the lower ocean layers on the continental shelf (Schmidtko et al., 2014) averaged over 18 separate basins adjacent to the ice shelves around the Antarctic continent. While salinity change over time in the deeper layers is neglected in this study, the ocean temperature responds with some delay to changes in the global mean temperature. We analyzed simulations with the coupled climate model ECHAM5–MPIOM over more than 6000 years following a 4-fold increase in CO₂ forcing (courtesy Li et al., 2012) and identified the anomaly in global mean temperature, in Antarctic temperature (south of 66^∘ S) and in Antarctic ocean temperature at depth levels between 500 and 2,500 m. After a response time of about 3000 years the ocean temperature stabilizes at about f_o=0.75 of the global mean anomaly, while Antarctic surface temperature anomaly is amplified by a factor of 1.8 (see Fig. 20a). As we intend to estimate ocean temperature change from ice surface temperature change reconstructed from ice cores, we could have fitted the response function directly from the Antarctic mean surface temperature in the climate model data. But we prefer to use a more general relationship, defined as response to global mean temperature change, which makes it easier to compare to other approaches. We found that, for the timescales considered here, Antarctic surface temperatures respond almost linearly with global mean temperature change. We therefore assume in the analysis that both time series are interchangeable except for the inferred amplification factor.

Figure 20Climate model data and response function analysis – (a) global mean temperature (green); Antarctic mean temperature south of 66^∘ S (grey); and ocean temperature averaged over upper 500 m (light blue), intermediate (500–2500 m; blue) and deeper layers (below 2500 m; violet) of the MPIOM model coupled to ECHAM5 forced by a CO₂ quadrupling within 140 years (Li et al., 2012). (b) Reconstructed response function R(t) (blue) with fitted function R^*(t) (orange dashed) as in Eq. (14). For comparison the fit function with α=1 (green dotted) is shown.

Download

Using linear response theory (Winkelmann and Levermann, 2013) and assuming the global mean temperature anomaly ΔT_GM(t) via a convolution integral related to the ocean temperature as

we reconstruct a response function R(t) and a corresponding fit function with $R^{*}\left(t\right)\sim t^{-\mathit{\alpha }}$ (see Fig. 20b), which vanishes beyond the typical response time τ_r. For α=2 and with an integration constant in the numerator, such that

is valid, this yields

The inferred response fit function convoluted with a given time series of global mean surface temperature anomaly (or here with the scaled ice core temperature reconstruction) hence provides an estimate for the corresponding change in ocean temperature at intermediate depth. Figure 21 shows the estimated ocean temperature anomaly curve (blue) with some delay with respect to the WDC surface temperature reconstructions (grey; Cuffey et al., 2016), here for easier comparison also shown scaled by $\mathrm{0.75}/\mathrm{1.8}=\mathrm{5}/\mathrm{12}=\mathrm{42}$ % (purple). The WDC likely better reflects the ocean conditions in the widely marine West Antarctic Ice Sheet than the EDC, which is located in central East Antarctica (Jouzel et al., 2007).

The inferred ocean temperature time series is comparably smooth with a resolution of 500 years and serves PICO as ocean temperature forcing. A comparison to reconstructions with a general circulation model (GCM) in the TraCE-21ka project (Liu et al., 2009)³, suggests that short warming periods above the present level could have occurred at intermediate depth, e.g. during Antarctic Cold Reversal (ACR) at around 14 kyr BP, which cannot be adequately resolved with our approach. The GCM ocean data are bounded below by the pressure-melting point. As negative ocean temperature anomalies can result in unphysical values below the pressure-melting point, we leave it to the PICO module to assert this lower bound such that melting vanishes and overturning circulations halts accordingly. As a consequence, ocean forcing becomes irrelevant for much-colder-than-present glacial climates. The parameterization presented here assumes that ocean water masses at depth below 500 m can access ice shelf cavities and induce melting, which is certainly very simplified regarding the complex bathymetry and flow patterns around Antarctica. Also, we used data from a simplified sensitivity experiment with ECHAM5–MPIOM, for a warmer-than-present climate, which implies various model uncertainties. We also had to make assumptions about a suitable response function, which is fitted to model data that are averaged over certain regions and ocean depths, implying further uncertainties.

Figure 21Response of the intermediate ocean temperature (blue) to surface temperature anomaly, as reconstructed from WDC (Cuffey et al., 2016, grey line), assuming a polar amplification factor of 1.8 and an equilibrium ocean scale factor of 0.75, as identified in climate model analysis for a warming scenario. Convolution yields a delayed response (blue) with respect to the forcing, shown here as scaled curve (violet), with dots indicating bins of 500-year means. The time series includes the Antarctic Cold Reversal (ACR) of 14.6–12.7 kyr BP (Fogwill et al., 2017). For comparison, GCM ocean temperature anomaly is shown (in olive), as reconstructed by TraCE-21ka (Liu et al., 2009), for the mid-ocean depth 600–2,800 m and south of 66^∘ S.

Download

Transient PISM simulations reveal a delayed response of the ice volume to the ocean temperature forcing, with the main deglaciation occurring when temperature anomalies reach almost present-day levels (Fig. 22; cf. upper and lower panels). When we directly apply the smoothed and scaled surface temperature forcing as PICO forcing, we still get very similar results but a slightly earlier retreat (orange). Simulations reveal that the power of the response function (cf. Fig. 20b) is of minor relevance to the ice sheet's response (compare grey and blue lines). Also the amplitude of cooling at glacial stage, here scaled by 60 % (light orange), shows only little change in the simulation results. If the ocean forcing is related to the EDC temperature reconstruction (see Sect. 4.2) we find a later warming and hence delayed deglaciation (light blue). Ocean forcing likely plays a key role in warmer-than-present climates. However, we do not see this effect in our simulations during the Last Interglacial. Although ocean temperatures rise by more than 1 K above the present we find similar ice volumes to the present in all our simulations. Precipitation scaling and the till properties seem to play an important role in stabilizing WAIS and preventing collapse. However, a thorough investigation of necessary model settings for WAIS collapse during the LIG would fill a separate study.

As we employ an anomaly forcing, which becomes near zero at present, the modeled modern ice sheet configurations are rather independent of the applied ocean temperature forcing (except for EDC, which is generally colder through the Holocene). Also, there is literature suggesting periods of decoupled ocean and surface temperature evolution (e.g., ACR) with strong potential effects on the ice sheet deglaciation. Appendix A1 provides an estimate of the effect of intermediate ocean warming events from a simple perturbation analysis with PICO.

4.4 Precipitation forcing

Continental-scale precipitation change is closely related to temperature change. While colder temperatures lead to higher effective ice viscosity and hence to slower ice flow and larger ice masses, they also lead to dryer conditions and hence to less ice mass accumulation. This effect is based on the Clausius–Clapeyron relationship, which suggests higher atmospheric moisture capacity and hence more accumulation in a warmer atmosphere. On average over the Antarctic continent, the analysis of ice core and modeling data suggest a linear scaling relationship between precipitation and temperature change of $f_{\mathrm{p},\mathrm{l}}=\mathrm{5}\pm \mathrm{1}$ % K⁻¹ (Frieler et al., 2015). In PISM simulations, precipitation forcing P(t) is coupled directly to the temperature forcing ΔT(t) (Sect. 4.2) using an exponential relationship (Ritz et al., 1996; Quiquet et al., 2012, Eq. 2), which scales the present-day mean precipitation field P₀ as

The exponential function is hence compatible with the precipitation lapse correction (Eq. 11), and it allows for easier comparisons with other studies that use power law relationships, e.g., with $f_{\mathrm{p}}=\ln \left(\mathrm{2}\right)/\mathrm{10}$ for Pollard and DeConto (2009). In fact, precipitation change depends very much on the Antarctic region. In our simulations we use a combination of temperature reconstructions from the EDC and WDC for ΔT(t) (see Sect. 4.2). Compared to the Antarctic mean value, the study by Frieler et al. (2015) suggests slightly higher values at core locations of $f_{\mathrm{p},\mathrm{l}}=\mathrm{5.9}\pm \mathrm{2.2}$ % K⁻¹ and $f_{\mathrm{p},\mathrm{l}}=\mathrm{5.5}\pm \mathrm{1.2}$ % K⁻¹, respectively. In East Antarctica estimates of the linear precipitation scaling factor tend to be slightly higher than in West Antarctica.

Figure 22Sensitivity of transient ice volume above flotation to varied ocean temperature forcing. Reference simulation (grey) is based on EDC + WDC surface temperature with fit response function of power α=2. The resulting ice volume change is similar for a fit response function with power α=1 (blue). Even if ocean temperature forcing responds directly to the scaled surface temperature anomaly (orange), the resulting effects on ice volume are comparably small. If a different surface temperature reconstruction is used, here from EDC with α=2 (light blue), we find, in response to a later warming signal, delayed deglaciation. If the amplitude of WDC surface temperature anomalies were 40 % lower (light orange), this would have only negligible effect on the modeled ice volume.

Download

In glacial periods with much colder temperatures of $\mathrm{\Delta }T=-\mathrm{10}$ K, an exponential precipitation change with f_p=7–9 % K⁻¹ yields 50 %–60 % less precipitation when compared to modern times. Without any temperature-scaled precipitation change, our simulations suggest up to 7 m SLE thicker ice sheets at the Last Glacial Maximum (cf. grey and light blue line in Fig. 23). The reference value of 7 % K⁻¹ precipitation change corresponds to more than 50 % dryer conditions than the present and about 3–4 m SLE less ice volume than for 5 % K⁻¹ (cf. grey and blue line). Reconstructions of precipitation at ice core sites, e.g., at WDC (Buizert et al., 2015; Fudge et al., 2016), reveal precipitation changes relative to present of up to −60 % during LGM and up to 25 % more precipitation through the Holocene. Figure 23 shows the corresponding transient effects of reconstructed precipitation forcings on the ice sheet's volume above flotation (orange line). The additional accumulation during the Holocene at WDC prohibits deglaciation and causes a modern ice sheet configuration that is more than 10 m SLE larger. However, the reconstructed signal at WDC may be biased to some extent from a lapse-rate effect due to surface elevation changes during deglaciation (personal communication with Eric Steig). Interestingly, we find a similar ice volume response when applying no lapse-rate correction for changing surface elevation (as described in Eq. 11), which also diminishes deglaciation (light orange in Fig. 23; for the same precipitation forcing as the reference). The simulations hence suggest that the precipitation scaling parameter f_p is highly relevant to the ice sheet's extent at glacial maximum and will be considered to be ensemble parameter PREC in Albrecht et al. (2020).

4.5 Combined effects of climatic forcings in glacial-cycle simulations

In our simulations with PISM the climatic forcings described above (Sect. 4.1–4.4) have different effects on the Antarctic Ice Sheet evolution at different periods. Both temperature and sea-level forcing reveal a periodic signature not only of the dominant 100 kyr orbital cycle but also with smaller amplitudes of the higher-frequency cycles, e.g., Marine Isotope Stage 3 (MIS 3) around 57 kyr BP (Weichselian High Glacial) or the MIS 2 at the Last Glacial Maximum around 29 kyr BP (Lisiecki and Raymo, 2005). For glacial-cycle simulations, sea-level forcing has the strongest effect on the ice volume, as it alone can trigger a larger glacial ice volume than in the reference simulation (see Fig. 24b; orange and grey). However, the rising sea level during deglaciation alone induces only little ice sheet retreat. When sea-level forcing is turned off, the other forcings balance each other out such that ice volume remains approximately at modern level through glacial periods (Fig. 24c). Also surface temperature forcing alone (Fig. 24b; blue) can produce glacial volumes of similar extent to the sea-level forcing. But if surface temperatures remain at modern levels, the other forcings can still produce a glacial maximum volume that is only 3 m SLE lower than the reference. In our simulations, surface temperature anomalies also drive changes in ocean temperature and precipitation. Ocean temperature forcing has only little effect on glacial extent, but it influences the onset of deglaciation (light orange). Without ocean forcing the interglacial (and modern) ice volume is hence up to 7 m SLE larger than in the reference (Fig. 24c). While sea-level and temperature forcings cause a growth of ice sheets at glacial climates, precipitation forcing has the opposite effect (Fig. 24; light blue). Without precipitation forcing the Antarctic Ice Sheet can reach glacial extents of up to 7 m SLE above the reference (coupled to the surface temperature forcing with f_p=7 % K⁻¹). This effect of the precipitation forcing for glacial climates also explains why the individual responses to the sea-level and the surface temperature forcing exceed the reference ice volume by about 3–4 m SLE, in which all four forcings are superimposed.

Figure 23Sensitivity of transient ice volume above flotation to varied precipitation forcing. Grey curve is the reference, with an exponential scaling factor of 7 % K⁻¹ with respect to WDC temperature reconstruction (Cuffey et al., 2016). Ice volume is slightly larger for 5 % K⁻¹ (blue) and much larger without precipitation forcing (light blue). Deglaciation is prohibited when using WDC accumulation reconstruction (Fudge et al., 2016) with additional accumulation during Holocene (orange). A similar effect is seen when no lapse-rate correction is applied for the precipitation (light orange).

Download

Figure 24Antarctic Ice Sheet volume above flotation for combinations of surface temperature, ocean temperature, sea-level and precipitation forcings (a), as defined in previous sections. (b) The ice sheet's response to individual forcings (or constant equilibrium conditions in purple) or (c) in the absence of one of the forcings, compared to the reference simulation (with all forcings combined; in grey). Sea-level forcing (orange) and surface temperature forcing (blue) have the strongest effects: each alone can cause larger glacial ice volume than in the reference simulation. However, each alone cannot initiate effective glacial retreat during modern climate. In contrast, ocean temperature forcing (light orange) has no effect on glacial volume, but it amplifies deglaciation. Precipitation forcing (light blue) counteracts sea-level and temperature forcings in glacial climates.

Download

Our results are in line with another recent study using the PSU-ISM ice sheet model. In this study, they also find a dominance of atmospheric and sea-level forcing on the Antarctic ice volume over the last four glacial cycles (Tigchelaar et al., 2019), which together drive glacial–interglacial ice volume changes of 12–14 m SLE, while ocean temperature forcing is almost negligible, also during interglacials. Here, we do not want to go into the details of this study, which uses a comparably coarse output of an Earth system model of intermediate complexity for the atmospheric and ocean forcing instead of scaling with ice core temperature reconstructions. As a key result, Tigchelaar et al. (2019) find much smaller individual ice volume changes, which amount to less than half of the full ice volume response. In our simulations, however, the individual response to sea-level forcing (and surface temperature forcing) and the sum of all individual forcings exceed the combined response. This is partly due to the fact that precipitation forcing (up to 50 % less during glacial climate) provides a strong negative effect on the ice volume in the full-forcing case, which seems to be weaker in the “atm” forcing in Tigchelaar et al. (ca. 15 % in Fig. 8; 2019). If we consider the LGM and present-day state to be rather stable states, a certain perturbation threshold must be crossed to initiate the (nonlinear) transition (Termination I) into the other state. In our simulations this threshold can be reached with individual forcings, while in the other study the combined superposition is required.

Table 2Sensitivity (mean and standard deviation) of simulated ice volume at Last Interglacial (LIG), Last Glacial Maximum (LGM) and present day (PD) for varied input datasets and model parameters (some of which are already indicated in Table 1), and the used reference value. Selected ensemble parameters for each category are bold: VISC and PPQ are more relevant to deglacial and present WAIS dynamics, and ESIA and PREC are relevant rather to EAIS dynamics. Asterisk indicates a different PISM version used.

Download Print Version | Download XLSX