Search

At present, exploring the space of rheological parameters in models of glacial isostatic adjustment (GIA) and relative sea level (RSL) which incorporate laterally variable Earth structure is computationally expensive. A single simulation using the Seakon model (Latychev et al., 2005), using contemporary high-performance computing hardware, requires several wall-days & ≈ 1 core-year for one RSL simulation from late Marine Isotope Stage 3 to present day. However, it is well established that the impact from laterally variable mantle viscosity and lithospheric thickness on RSL and GIA is significant (Whitehouse, 2018). We present initial results from using the Tensorflow (Abadi et al.) framework to construct artificial neural networks that emulate the difference in the rate of change of relative sea level and relative radial displacement between model configurations using spherically symmetric (SS) and laterally variable (LV) Earth structures. Using this emulator we can accurately sample the parameter space (≈ 360 realisations of the background (SS) structure) for a given realization of lateral Earth structure (e.g. viscosity variations derived from shear-wave tomographic models) using ≈ 1/10th the amount of parameter vectors as a training set. Average misfits are O(0.1-1%) of the total RSL signal when using the emulator to adjust SS GIA model output to incorporate the impact from LV. We shall report on two case studies which allow us to examine the influence of lateral Earth structure on inferences of background (i.e. global-mean) viscosity. For these case studies, the emulator, in conjunction with a fast SS GIA/RSL model, is used to determine optimal Earth model parameters (elastic lithosphere thickness, upper and lower mantle viscosities) by calculating the model misfits across the parameter space. The first case study uses the regional RSL database of Vacchi et al. (2018) which spans the Canadian Arctic and East Coast with several hundred sea level index points and limiting points for the early to late Holocene. The second case study uses a global database of several thousand contemporary uplift rates derived from GPS data (Schumacher et al., 2018). For the first case study we find two main features from incorporating LV structures compared to the SS configuration: a decrease in the best scoring misfit and a shift of the misfit distribution in the parameter space to favour a reduced upper mantle viscosity and reduced sensitivity to the lower mantle viscosity.ReferencesAbadi, M., Agarwal, A., Barham, P., et al.: TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, https://www.tensorflow. org/.Latychev, K., Mitrovica, J. X., Tromp, J., et al.: Glacial isostatic adjustment on 3-D Earth models: a finite-volume formulation, GJI, 161, 421–444, https://doi.org/10.1111/j.1365-246x.2005.02536.x, 2005.Schumacher, M., King, M. A., Rougier, J., et al.: A new global GPS data set for testing and improving modelled GIA uplift rates, GJI, 214, 2164–2176, https://doi.org/10.1093/gji/ggy235, 2018.Vacchi, M., Engelhart, S. E., Nikitina, D., et al.: Postglacial relative sea-level histories along the eastern Canadian coastline, QSR, 201, 124–146, https://doi.org/10.1016/j.quascirev.2018.09.043, 2018.Whitehouse, P. L.: Glacial isostatic adjustment modelling: historical perspectives, recent advances, and future directions, Earth Surface Dynamics, 6, 401–429, https://doi.org/10.5194/esurf-6-401-2018, 2018.

As in any environmental system, modeling instabilities within the glacial system is a numerical challenge of potentially high real-world relevance. Differentiating between the impacts of physical system processes and numerical noise is not straightforward. Here we use an idealized North American geometry and climate representation (similar to the HEINO experiments, Calov et al., 2010) to examine the numerical sensitivity of ice stream surge cycling in glaciological models. Through sensitivity tests, we identify some numerical requirements for a robust model configuration for such contexts. To partly address model-specific dependencies, we use both the Glacial Systems Model (GSM) and Parallel Ice Sheet Model (PISM). We show that modeled surge characteristics are resolution-dependent though converging (decreasing differences between resolutions) at higher horizontal grid resolutions. Discrepancies between high and coarse horizontal grid resolutions can be reduced by incorporating a resolution-dependent basal temperature ramp for basal sliding thermal activation. Inclusion of a diffusive bed thermal model reduces the surge cycling ice volume change by ∼33 % as the additional heat storage dampens the change in basal temperature during surge events. The inclusion of basal hydrology, as well as a non-flat topography, leads to increased ice volume change during surge events (∼20 and 17 %, respectively). Therefore, these latter three components are essential if one is endeavoring to maximize physical fidelity in ice stream surge cycle modeling. An abrupt transition between hard bedrock and soft sediment, as in the HEINO experiments, leads to ice stream propagation along this boundary but is not the cause of the main surge events.

The timing of Laurentide Ice Sheet deglaciation along its southwestern margin controlled the evolution of large glacial lakes and has implications for human migration into the Americas. Accurate reconstruction of the ice sheet's retreat also constrains glacial isostatic adjustment models and is important for understanding ice-sheet sensitivity to climate forcing. Despite its significance, retreat of the southwestern Laurentide Ice Sheet (SWLIS) is poorly constrained by minimum-limiting 14C data. We present 26 new cosmogenic 10Be exposure ages spanning the western Interior Plains, Canada. Using a Bayesian framework, we combine these data with geomorphic mapping, 10Be, and high-quality minimum-limiting 14C ages to provide an updated chronology. This dataset presents an internally consistent retreat record and indicates that the initial detachment of the SWLIS from its convergence with the Cordilleran Ice Sheet began by ca. 15.0 ka, concurrent with or slightly prior to the onset of the Bølling-Allerød interval (14.7–12.9 ka) and retreated >1200 km to its Younger Dryas (YD) position in ~2500 yr. Ice-sheet stabilization at the Cree Lake Moraine facilitated a meltwater drainage route to the Arctic from glacial Lake Agassiz within the YD, but not necessarily at the beginning. Our record of deglaciation and new YD constraints demonstrate deglaciation of the Interior Plains was ~60% faster than suggested by minimum 14C constraints alone. Numerical modeling of this rapid retreat estimates a loss of ~3.7 m of sea-level equivalent from the SWLIS during the Bølling-Allerød interval.

While subglacial hydrology is known to play a role in glacial dynamics on sub-annual to decadal scales, it remains unclear whether subglacial hydrology plays a critical role in ice sheet evolution on centennial or longer time-scales. Furthermore, several drainage topologies have been inferred but it is unclear which drainage topology is most applicable at the continental/glacial scale. More fundamentally, it is even unclear if the structural choice of subglacial hydrology truly matters for this context. Here we compare three subglacial hydrology topologies as to their contribution to surge behaviour for an idealized Hudson Strait like ice stream. We use the newly updated model BraHms2.0 and provide model verification tests. BraHms2.0 incorporates each of these systems: two process-based forms dominant in the literature (linked-cavity and poro-elastic) and a non-mass conserving zero-dimensional form (herein termed leaky-bucket) coupled to an ice sheet systems model (the Glacial Systems Model, GSM). We also assess the likely bounds on poorly constrained subglacial hydrology parameters and adopt an ensemble approach to study their impact and interactions within those bounds. We find that subglacial hydrology is an important system inductance for realistic ice stream surging but that the three formulations all exhibit similar surge behaviour. Even a detail as fundamental as mass conserving transport of subglacial water is not necessary for simulating a full range of surge frequency and amplitude. However, one difference is apparent: the combined positive and negative feedbacks of the linked-cavity system yields longer duration surges and a broader range of effective pressures than its poro-elastic and leaky-bucket counterparts.

We present an ensemble of last glacial inception (LGI) simulations for the Northern Hemisphere that captures a significant fraction of inferred ice volume changes within proxy uncertainties. This ensemble was performed with LCice 1.0, a coupled ice sheet and climate model, varying parameters of both climate and ice sheet components, as well as the coupling between them. Certain characteristics of the spatiotemporal pattern of ice growth and subsequent retreat in both North America (NA) and Eurasia (EA) are sensitive to parameter changes while others are not. We find that the initial inception of ice over NA and EA is best characterized by the nucleation of ice at high-latitude and high-elevation sites. Subsequent spreading and merger along with large-scale conversion of snowfields dominate in different sectors. The latter plays an important role in the merging of eastern and western ice regions in NA. The inception peak ice volume in the ensemble occurs approximately at 111 ka and therefore lags the summer 60∘ N insolation minimum by more than 3 kyr. Ice volumes consistently peak earlier over EA than NA. The inception peak in North America is characterized by a merged Laurentide and Cordilleran ice sheet, with the Davis Strait covered in ice in ∼80 % of simulations. Ice also bridges Greenland and Iceland in all runs by 114 ka and therefore blocks the Denmark Strait. This latter feature would thereby divert the East Greenland Current and Denmark Strait overflow with a potentially significant impact on ocean circulation. The Eurasian ice sheet at its inception peak varies across ensemble runs between a continuous ice sheet and multiple smaller ice caps. In both continents, the colder high latitudes (i.e. Ellesmere and Svalbard) tend to grow ice through the entire simulation (to 102 ka), while lower latitudes lose ice after ∼110 ka. We find temperature decreases over the initial phases of the inception lead to the expansion of NA ice sheet area and that subsequent precipitation increases contribute to its thickening. EA ice sheet area also expands with decreasing temperatures, but sea ice limits any increases in precipitation, leading to an earlier retreat away from the EA maximum ice sheet volume. We also examine the extent to which the capture of both LGI ice growth and retreat constrains the coupled ice–climate model sensitivity to changing atmospheric pCO2. The 55-member sub-ensemble that meets our criteria for “acceptable” ice growth and retreat has an equilibrium climate sensitivity lower bound that is 0.3 ∘C higher than that of the full ensemble. This suggests some potential value of fully coupled ice–climate modelling of the last glacial inception to constrain future climate change.

Highlights: • Climate model sensitivity experiments are performed using state-of-the-art ice sheet and freshwater reconstructions • Declining Northern Hemisphere ice sheets increase the sensitivity of the AMOC to North Atlantic meltwater discharge • Deglacial rise in atmospheric CO2 concentration decreases the sensitivity of the AMOC to North Atlantic meltwater discharge • Both effects provide a complementary perspective to existing explanations for abrupt AMOC transitions Abstract: The last deglaciation was characterized by a sequence of abrupt climate events thought to be linked to rapid changes in Atlantic meridional overturning circulation (AMOC). The sequence includes a weakening of the AMOC after the Last Glacial Maximum (LGM) during Heinrich Stadial 1 (HS1), which ends with an abrupt AMOC amplification at the transition to the Bølling/Allerød (B/A). This transition occurs despite persistent deglacial meltwater fluxes that counteract vigorous North Atlantic deep-water formation. Using the Earth system model COSMOS with a range of deglacial boundary conditions and reconstructed deglacial meltwater fluxes, we show that deglacial CO2 rise and ice sheet decline modulate the sensitivity of the AMOC to these fluxes. While declining ice sheets increase the sensitivity, increasing atmospheric CO2 levels tend to counteract this effect. Therefore, the occurrence of a weaker HS1 AMOC and an abrupt AMOC increase in the presence of meltwater, might be explained by these effects, as an alternative to or in combination with changes in the magnitude or routing of meltwater discharge.

In the geosciences, complex computational models have become a common tool for making statements about past earth system evolution. However, the relationship between model output and the actual earth system (or component thereof) is generally poorly specified and even more poorly assessed. This is especially challenging for the paleo sciences for which data constraints are sparse and have large uncertainties. Bayesian inference offers, in principle, a self-consistent and rigorous framework for assessing this relationship as well as a coherent approach to combining data constraints with computational modelling. Though “Bayesian” is becoming more common in paleoclimate and paleo ice sheet publications, our impression is that most scientists in these fields have little understanding of what this actually means nor are they able to evaluate the quality of such inference. This is especially unfortunate given the correspondence between Bayesian inference and the classical concept of the scientific method. Herein, we examine the relationship between a complex model and a system of interest, or in equivalent words (from a statistical perspective), how uncertainties describing this relationship can be assessed and accounted for in a principled and coherent manner. By way of a simple example, we show how inference, whether Bayesian or not, can be severely broken if uncertainties are erroneously assessed. We explain and decompose Bayes Rule (more commonly known as Bayes Theorem), examine key components of Bayesian inference, offer some more robust and easier to attain stepping stones, and provide suggestions on implementation and how the community can move forward. This overview is intended for all interested in making and/or evaluating inferences about the past evolution of the Earth system (or any of its components), with a focus on past ice sheet and climate evolution during the Quaternary.

Abrupt climate shifts of large amplitudes were common features of the Earth's climate as it transitioned into and out of the last full glacial state approximately 20 000 years ago, but their causes are not yet established. Midlatitude atmospheric dynamics may have played an important role in these climate variations through their effects on heat and precipitation distributions, sea ice extent, and wind-driven ocean circulation patterns. This study characterizes deglacial winter wind changes over the North Atlantic (NAtl) in a suite of transient deglacial simulations using the PlaSim Earth system model (run at T42 resolution) and the TraCE-21ka (T31) simulation. Though driven with yearly updates in surface elevation, we detect multiple instances of NAtl jet transitions in the PlaSim simulations that occur within 10 simulation years and a sensitivity of the jet to background climate conditions. Thus, we suggest that changes to the NAtl jet may play an important role in abrupt glacial climate changes. We identify two types of simulated wind changes over the last deglaciation. Firstly, the latitude of the NAtl eddy-driven jet shifts northward over the deglaciation in a sequence of distinct steps. Secondly, the variability in the NAtl jet gradually shifts from a Last Glacial Maximum (LGM) state with a strongly preferred jet latitude and a restricted latitudinal range to one with no single preferred latitude and a range that is at least 11∘ broader. These changes can significantly affect ocean circulation. Changes to the position of the NAtl jet alter the location of the wind forcing driving oceanic surface gyres and the limits of sea ice extent, whereas a shift to a more variable jet reduces the effectiveness of the wind forcing at driving surface ocean transports. The processes controlling these two types of changes differ on the upstream and downstream ends of the NAtl eddy-driven jet. On the upstream side over eastern North America, the elevated ice sheet margin acts as a barrier to the winds in both the PlaSim simulations and the TraCE-21ka experiment. This constrains both the position and the latitudinal variability in the jet at LGM, so the jet shifts in sync with ice sheet margin changes. In contrast, the downstream side over the eastern NAtl is more sensitive to the thermal state of the background climate. Our results suggest that the presence of an elevated ice sheet margin in the south-eastern sector of the North American ice complex strongly constrains the deglacial position of the jet over eastern North America and the western North Atlantic as well as its variability.

We model the change of Coulomb Failure Stress (δCFS) during the Weichselian glaciation up until today at 12 locations in Latvia, Lithuania and Russia that are characterised by soft-sediment deformation structures (SSDS). If interpreted as seismites, these SSDS may point to glacially-induced fault reactivation. The δCFS suggests a high potential of such reactivation when it reaches the instability zone. We show that δCFS at all 12 locations reached this zone several times in the last 120,000 years. Most notably, all locations exhibit the possibility of reactivation after ca. 15 ka BP until today. Another time of possible activity likely happened after the Saalian glaciation until ca. 96 ka BP. In addition, some models suggest unstable states after 96 ka BP until ca. 28 ka BP at selected locations but with much lower positive δCFS values than during the other two periods. For the Valmiera and Rakuti seismites in Latvia, we can suggest a glacially-induced origin, whereas we cannot exactly match the timing at Rakuti. Given the (preliminary) dating of SSDS at some locations, at Dyburiai and Ryadino our modelling supports the interpretation of glacially-induced fault reactivation, while at Slinkis, Kumečiai and Liciškėnai they likely exclude such a source. Overall, the mutual benefit of geological and modelling investigations is demonstrated. This helps in identifying glacially-induced fault reactivation at the south-eastern edge of the Weichselian glaciation and in improving models of glacial isostatic adjustment.

We present the official land uplift model NKG2016LU of the Nordic Commission of Geodesy (NKG) for northern Europe. The model was released in 2016 and covers an area from 49° to 75° latitude and 0° to 50° longitude. It shows a maximum absolute uplift of 10.3 mm/a near the city of Umea in northern Sweden and a zero-line that follows the shores of Germany and Poland. The model replaces the NKG2005LU model from 2005. Since then, we have collected more data in the core areas of NKG2005LU, specifically in Norway, Sweden, Denmark and Finland, and included observations from the Baltic countries as well. Additionally, we have derived an underlying geophysical glacial isostatic adjustment (GIA) model within NKG as an integrated part of the NKG2016LU project. A major challenge is to estimate a realistic uncertainty grid for the model. We show how the errors in the observations and the underlying GIA model propagate through the calculations to the final uplift model. We find a standard error better than 0.25 mm/a for most of the area covered by precise levelling or uplift rates from Continuously Operating Reference Stations and up to 0.7 mm/a outside this area. As a check, we show that two different methods give approximately the same uncertainty estimates. We also estimate changes in the geoid and derive an alternative uplift model referring to this rising geoid. Using this latter model, the maximum uplift in Umea reduces from 10.3 to 9.6 mm/a and with a similar reduction ratio elsewhere. When we compare this new NKG2016LU with the former NKG2005LU, we find the largest differences where the GIA model has the strongest influence, i.e. outside the area of geodetic observation. Here, the new model gives from − 3 to 4 mm/a larger values. Within the observation area, similar differences reach − 1.5 mm/a at the northernmost part of Norway and − 1.0 mm/a at the north-western coast of Denmark, but generally within the range of − 0.5 to 0.5 mm/a.

Projections of the sea level contribution from the Greenland and Antarctic ice sheets (GrIS and AIS) rely on atmospheric and oceanic drivers obtained from climate models. The Earth System Models participating in the Coupled Model Intercomparison Project phase 6 (CMIP6) generally project greater future warming compared with the previous Coupled Model Intercomparison Project phase 5 (CMIP5) effort. Here we use four CMIP6 models and a selection of CMIP5 models to force multiple ice sheet models as part of the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6). We find that the projected sea level contribution at 2100 from the ice sheet model ensemble under the CMIP6 scenarios falls within the CMIP5 range for the Antarctic ice sheet but is significantly increased for Greenland. Warmer atmosphere in CMIP6 models results in higher Greenland mass loss due to surface melt. For Antarctica, CMIP6 forcing is similar to CMIP5 and mass gain from increased snowfall counteracts increased loss due to ocean warming.

Freshwater, in the form of glacial runoff, is hypothesized to play a critical role in centennial to millennial scale climate variability such as the Younger Dryas and Dansgaard-Oeschger Events. Indeed, freshwater injection/hosing experiments with climate models have long shown that freshwater has the capability of generating such abrupt climate transitions. However, the relationship between freshwater and abrupt climate transitions is not straightforward. Large-scale glacial runoff events, such as Meltwater Pulse 1A, are not always temporally proximal to subsequent large-scale cooling. As well, the typical design of hosing experiments tends to artificially amplify the climate response. This study explores the impact that limitations in the representation of runoff in conventional hosing simulations has on our understanding of this relationship and addresses the more fundamental question of where coastally released freshwater is transported when it reaches the ocean. We focus particularly on the prior use of excessive freshwater volumes (often by a factor of 5) and present-day (rather than paleo) ocean gateways, as well as the injection of freshwater directly over sites of deep-water formation (DWF) rather than at runoff locations. We track the routing of glaciologically-constrained freshwater volumes from four different plausible injection locations in a suite of eddy-permitting glacial ocean simulations using MITGCM under both open and closed Bering Strait conditions. Restricting freshwater forcing values to realistic ranges results in less spreading of freshwater across the North Atlantic and indicates that the response of DWF depends strongly on the geographical location of meltwater input. In particular, freshwater released into the Gulf of Mexico has little impact on DWF regions as a result of turbulent mixing by the Gulf Stream. In contrast, freshwater released from the Eurasian Ice sheet or initially into the Arctic is found to have the largest impact on DWF in the North Atlantic and GIN seas. Additional experiments show that when the Bering Strait is open, much like present-day, the Mackenzie River source exhibits twice as much freshening of the Labrador sea as a closed Bering Strait. Finally, our results illustrate that applying a freshwater hosing directly into the North Atlantic with even realistic freshwater amounts still over-estimates the effect of terrestrial runoff on ocean circulation.

Projections of the sea level contribution from the Greenland and Antarctic ice sheets (GrIS and AIS) rely on atmospheric and oceanic drivers obtained from climate models. The Earth System Models participating in the Coupled Model Intercomparison Project phase 6 (CMIP6) generally project greater future warming compared with the previous Coupled Model Intercomparison Project phase 5 (CMIP5) effort. Here we use four CMIP6 models and a selection of CMIP5 models to force multiple ice sheet models as part of the Ice Sheet Model Intercomparison Project for CMIP6 (ISMIP6). We find that the projected sea level contribution at 2100 from the ice sheet model ensemble under the CMIP6 scenarios falls within the CMIP5 range for the Antarctic ice sheet but is significantly increased for Greenland. Warmer atmosphere in CMIP6 models results in higher Greenland mass loss due to surface melt. For Antarctica, CMIP6 forcing is similar to CMIP5 and mass gain from increased snowfall counteracts increased loss due to ocean warming.

This archive contains a high variance subset of last glacial inception simulations from LCice (coupled ice and climate model). Details in contained tableOfRuns.org and associated journal submission under review: Bahadory, T., Tarasov, L., and Andres, H.: The phase space of last glacial inception for the Northern Hemisphere from coupled ice and climate modelling, Clim. Past, https://doi.org/10.5194/cp-2020-1, in review, 2020. Contact: lev@mun.ca, {"references":["Bahadory, T., Tarasov, L., and Andres, H.: The phase space of last glacial inception for the Northern Hemisphere from coupled ice and climate modelling, Clim. Past, https://doi.org/10.5194/cp-2020-1, in review, 2020"]}

This study deals with modelling the distribution of the subsurface pore pressure in space and the respective evolution in time in response to variations in hydromechanical surface loading during a full glacial-interglacial cycle. The aim here is to better understand (i) the feedback mechanisms between the atmosphere and solid earth components, and (ii) to which degree this coupling might be relevant for subsurface hydromechanical modelling studies. The study area is the Central European Basin System (CEBS) in northern and Central Europe and state-of-the-art ice reconstructions for the last glacial-interglacial period have been used to model the surface hydromechanical loading conditions. Thereby, investigations on how transient ice coverage influences the pore pressure distribution with depth and over time within a heterogeneous sedimentary cover were carried out. The subsurface beneath the CEBS consists of more than 10 km thick sediments, which have been heavily restructured by salt movements during the whole Mesozoic evolution. Our 3D geological model resolves all major sedimentary and crustal domains, and we relied on the GLAC1-D (1.0 degree longitude by 0.5 degree latitude spatial resolution) ice sheet chronology. Starting from ice-free initial conditions, transient simulation runs are performed (hydraulic vs hydromechanical) which cover the entire last glacial cycle, i.e. encompassing 122ka BP till present day conditions. Results are discussed in terms of pore pressure evolution over time and space. The focus will lie on quantifying subsurface conditions favourable to the establishment and maintenance of overpressure evolution and the related equilibration time within the sedimentary pile. We also investigate how these transient conditions influence the subsurface hydrodynamics, showcasing representative time steps during the evolution of the system. We will finally attempt to quantify the memory effect of such loading conditions on the basin-wide hydromechanics, a feedback mechanism that has been neglected so far in 3D subsurface studies.

At the centers of previously glaciated regions such as Hudson Bay in Canada and the Gulf of Bothnia in Fennoscandia, it has been observed that the sea level history follows an exponential form and that the associated decay time is relatively insensitive to uncertainty in the ice loading history. We revisit the issue of decay time sensitivity by computing relative sea level histories for Richmond Gulf and James Bay in Hudson Bay and Ångerman River in Sweden for a suite of reconstructions of the North American and Fennoscandian Ice Sheets and Earth viscosity profiles. We find that while some Earth viscosity models do indeed show insensitivity in computed decay times to the ice history, this is not true in all cases. Moreover, we find that the location of the study site relative to the geometry of the ice sheet is an important factor in determining ice sensitivity, and based on our set of ice sheet reconstructions, conclude that the location of James Bay is not well-suited to a decay time analysis. We describe novel corrections to the RSL data to remove the effects associated with the spatial distribution of sea level indicators as well as for other signals unrelated to regional ice loading (ocean loading, rotation and global mean sea-level changes) and demonstrate that they can significantly affect the inference of viscosity structure. We performed a forward modelling analysis based on a commonly adopted 2-layer, sub-lithosphere viscosity structure to determine how the solution space of viscosity models changes with the input ice history at the three study sites. While the solution spaces depend on ice history, for both Richmond Gulf and Ångerman River there are regions of parameter space where solutions are common across all or most ice histories, indicating low ice load sensitivity for these mantle viscosity parameters. For example, in Richmond Gulf, upper mantle viscosity values of (0.3–0.5)x1021 Pa s and lower mantle viscosity values of (5–50)x1021 Pa s tend to satisfy the data constraint consistently for most ice histories considered in this study. Similarly, the Ångerman River solution spaces contain a solution with an upper mantle viscosity of 0.3 × 1021 Pa s and lower mantle viscosity values of (5–50)x1021 Pa s common to 9 of the 10 ice histories considered there. However, the dependence of the viscosity solution space on ice history suggests that joint estimation of ice and Earth parameters is the optimal approach.

Marine isotope stage (MIS) 3 (∼58–28 ka) is a period of intermediate global ice volume between MIS 4 and the last glacial maximum of MIS 2. Here we report geologic evidence for southern Laurentide ice-sheet rapid growth to near its last glacial maximum extent after a period with limited ice in the southernmost Hudson Bay lowland. A14C age on wood in lacustrine sediments interbedded with glacial tills in central-eastern Wisconsin dates a Laurentide ice-sheet advance southwards to an extent at least equivalent to at least its ∼17 ka deglacial limit by 39.1 ± 0.4 ka. This advance ended before 30.4 ± 0.9 ka based on another 14C date on wood in lacustrine sediment overlying the till layers. This advance is consistent with 14C ages from Michigan and Iowa, and Gulf of Mexico runoff records that support a concurrent southern Laurentide ice-sheet advance. We infer changes in North American ice volume using ice-sheet model simulations from a large ensemble that are consistent with 14C-data and Gulf of Mexico-discharge constraints. The simulations show the Laurentide ice sheet growing from a volume equivalent to 25–30 m of global mean sea level (GMSL) before ∼40 ka to 40–45 m of GMSL at ∼40 ka, and reaching 65–70 m GMSL by ∼30 ka, consistent with glacial isostatic adjustment assessments of near-to intermediate-field sea-level data. We thus show from our terrestrial field data and ice-sheet model simulations that an individual ice sheet can grow rapidly, which has only been inferred previously for global ice volume from GSML records.

We present BrAHMs (BAsal Hydrology Model): a physically based basal hydrology model which represents water flow using Darcian flow in the distributed drainage regime and a fast down-gradient solver in the channelized regime. Switching from distributed to channelized drainage occurs when appropriate flow conditions are met. The model is designed for long-term integrations of continental ice sheets. The Darcian flow is simulated with a robust combination of the Heun and leapfrog–trapezoidal predictor–corrector schemes. These numerical schemes are applied to a set of flux-conserving equations cast over a staggered grid with water thickness at the centres and fluxes defined at the interface. Basal conditions (e.g., till thickness, hydraulic conductivity) are parameterized so the model is adaptable to a variety of ice sheets. Given the intended scales, basal water pressure is limited to ice overburden pressure, and dynamic time stepping is used to ensure that the Courant–Friedrichs–Lewy (CFL) condition is met for numerical stability. The model is validated with a synthetic ice sheet geometry and different bed topographies to test basic water flow properties and mass conservation. Synthetic ice sheet tests show that the model behaves as expected with water flowing down gradient, forming lakes in a potential well or reaching a terminus and exiting the ice sheet. Channel formation occurs periodically over different sections of the ice sheet and, when extensive, displays the arborescent configuration expected of Röthlisberger channels. The model is also shown to be stable under high-frequency oscillatory meltwater inputs.

Sea level rise presents a hazard for coastal populations, and the Mississippi Delta (MD) is a region particularly at risk due to the high rates of land subsidence. We apply a gravitationally self-consistent model of glacial and sediment isostatic adjustment (SIA) along with a realistic sediment load reconstruction in this region for the first time to determine isostatic contributions to relative sea level (RSL) and land motion. We determine optimal model parameters (Earth rheology and ice history) using a new high-quality compaction-free sea level indicator database. Using the optimal model parameters, we show that SIA can lower predicted RSL in the MD area by several meters over the Holocene and so should be taken into account when modeling these data. We compare modeled contemporary rates of vertical land motion with those inferred using GPS. This comparison indicates that isostatic processes can explain the majority of the observed vertical land motion north of latitude 30.7°N, where subsidence rates average about 1 mm/yr; however, subsidence south of this latitude shows large data-model discrepancies of greater than 3 mm/yr, indicating the importance of nonisostatic processes. This discrepancy extends to contemporary RSL change, where we find that the SIA contribution in the Delta is on the order of 10−1 mm/yr. We provide estimates of the isostatic contributions to 20th and 21st century sea level rates at Gulf Coast Permanent Service for Mean Sea Level tide gauge locations as well as vertical and horizontal land motion at GPS station locations near the MD.

The Antarctic Ice Sheet is an important indicator of climate change and driver of sea-level rise. Here we combine satellite observations of its changing volume, flow and gravitational attraction with modelling of its surface mass balance to show that it lost 2,720 ± 1,390 billion tonnes of ice between 1992 and 2017, which corresponds to an increase in mean sea level of 7.6 ± 3.9 millimetres (errors are one standard deviation). Over this period, ocean-driven melting has caused rates of ice loss from West Antarctica to increase from 53 ± 29 billion to 159 ± 26 billion tonnes per year; ice-shelf collapse has increased the rate of ice loss from the Antarctic Peninsula from 7 ± 13 billion to 33 ± 16 billion tonnes per year. We find large variations in and among model estimates of surface mass balance and glacial isostatic adjustment for East Antarctica, with its average rate of mass gain over the period 1992–2017 (5 ± 46 billion tonnes per year) being the least certain.

We present an ensemble of Last Glacial Inception (LGI) simulations for the Northern Hemisphere that largely captures inferred ice volume changes within proxy uncertainties. This ensemble was performed with LCice 1.0, a coupled ice sheet and climate model, varying parameters of both climate and ice sheet components, as well as the coupling between them. Certain characteristics of the spatio-temporal pattern of ice growth and subsequent retreat in both North America (NA) and Eurasia (EA) are sensitive to parameter changes, especially with respect to regional rates of ice growth and retreat. We find that the initial inception of ice over NA and EA is best characterized by the nucleation of ice at high latitude and high elevation sites. Subsequent spreading and merger along with large-scale conversion of snow fields dominate in different sectors. The latter plays an important role in the merging of eastern and western ice regions in NA. The inception peak ice volume in the ensemble occurs approximately at 111 ka and therefore lags the summer 60° N insolation minimum by more than 3 kyr. Ice volumes consistently peak earlier over EA than NA. The inception peak in North America is characterized by a merged Laurentide and Cordilleran ice sheet, with Davis Strait covered in ice in 80 % of simulations. Ice also bridges Greenland and Iceland in all runs by 114 ka and therefore blocks Denmark Strait. This latter feature would thereby divert the East Greenland Current and Denmark Strait overflow and thereby potentially have a significant impact on ocean circulation. The Eurasian ice sheet at its inception peak varies across ensemble runs between a continuous ice sheet to multiple smaller ice caps. In both continents, the colder high latitudes (Ellsmere and Svalbard) tend to grow ice through the entire simulation (to 102 ka), while lower latitudes lose ice after 110 ka. We find temperature decreases over the initial phases of the inception lead to the expansion of NA ice sheet area, and that subsequent precipitation increases contribute to its thickening. EA ice sheet area also expands with decreasing temperatures, but sea ice limits any increases in precipitation, leading to an earlier retreat away from the EA maximum ice sheet volume. We also examine the extent to which the capture of both LGI ice growth and retreat constrains the coupled ice/climate model sensitivity to changing atmospheric pCO2. For a standard transient climate response experiment (1 % increase in pCO2 until doubled), warming ranges between 0.6–2.0 °C for our initial set of 500 simulations without LGI constraint. The warming is reduced to 0.7–1.4 °C for the 55 member ensemble that captures both LGI ice growth and retreat. This therefore underlines the potential value of fully coupled ice/climate modelling of last glacial inception to constrain future climate change.

The North American Ice Sheet Complex (NAISC; consisting of the Laurentide, Cordilleran and Innuitian ice sheets) was the largest ice mass to repeatedly grow and decay in the Northern Hemisphere during the Quaternary. Understanding its pattern of retreat following the Last Glacial Maximum is critical for studying many facets of the Late Quaternary, including ice sheet behaviour, the evolution of Holocene landscapes, sea level, atmospheric circulation, and the peopling of the Americas. Currently, the most up-to-date and authoritative margin chronology for the entire ice sheet complex is featured in two publications (Geological Survey of Canada Open File 1574 [Dyke et al., 2003]; 'Quaternary Glaciations - Extent and Chronology, Part II' [Dyke, 2004]). These often-cited datasets track ice margin recession in 36 time slices spanning 18 ka to 1 ka (all ages in uncalibrated radiocarbon years) using a combination of geomorphology, stratigraphy and radiocarbon dating. However, by virtue of being over 15 years old, the ice margin chronology requires updating to reflect new work and important revisions. This paper updates the aforementioned 36 ice margin maps to reflect new data from regional studies. We also update the original radiocarbon dataset from the 2003/2004 papers with 1541 new ages to reflect work up to and including 2018. A major revision is made to the 18 ka ice margin, where Banks and Eglinton islands (once considered to be glacial refugia) are now shown to be fully glaciated. Our updated 18 ka ice sheet increased in areal extent from 17.81 to 18.37 million km(2), which is an increase of 3.1% in spatial coverage of the NAISC at that time. Elsewhere, we also summarize, region-by-region, significant changes to the deglaciation sequence. This paper integrates new information provided by regional experts and radiocarbon data into the deglaciation sequence while maintaining consistency with the original ice margin positions of Dyke et al. (2003) and Dyke (2004) where new information is lacking; this is a pragmatic solution to satisfy the needs of a Quaternary research community that requires up-to-date knowledge of the pattern of ice margin recession of what was once the world's largest ice mass. The 36 updated isochrones are available in PDF and shapefile format, together with a spreadsheet of the expanded radiocarbon dataset (n = 5195 ages) and estimates of uncertainty for each interval. (C) 2020 Elsevier Ltd. All rights reserved.

The original version of this Article contained an error in Fig. 5, in which the grey line representing ‘EPICA Dust’ was inadvertently omitted from Fig. 5h. This has been corrected in both the PDF and HTML versions of the Article.

The cycling of carbon on Earth exerts a fundamental influence upon the greenhouse gas content of the atmosphere, and hence global climate over millennia. Until recently, ice sheets were viewed as inert components of this cycle and largely disregarded in global models. Research in the past decade has transformed this view, demonstrating the existence of uniquely adapted microbial communities, high rates of biogeochemical/physical weathering in ice sheets and storage and cycling of organic carbon (>104 Pg C) and nutrients. Here we assess the active role of ice sheets in the global carbon cycle and potential ramifications of enhanced melt and ice discharge in a warming world., Ice sheets have long been overlooked as regulators of the global carbon cycle. In this Review, Wadham and colleagues show how the growth and retreat of ice sheets support the productivity of the oceans and variably store or release organic carbon–in effect, these frozen landscapes must be considered in future assessments of climate impacts on biogeochemical cycling.

The penultimate deglaciation (PDG, ~ 138–128 thousand years before present, hereafter ka) is the transition from the penultimate glacial maximum to the Last Interglacial (LIG, ~ 129–116 ka). The LIG stands out as one of the warmest interglacials of the last 800 ka, with high-latitude temperature warmer than today and global sea level likely higher by at least 6 meters. Considering the transient nature of the Earth system, the LIG climate and ice-sheets evolution were certainly influenced by the changes occurring during the penultimate deglaciation. It is thus important to investigate, with coupled Atmosphere-Ocean General Circulation Models (AOGCMs), the climate and environmental response to the large changes in boundary conditions (i.e. orbital configuration, atmospheric greenhouse gas concentrations, ice-sheet geometry, and associated meltwater fluxes) occurring during the penultimate deglaciation. A deglaciation working group has recently been set up as part of the Paleoclimate Modelling Intercomparison Project (PMIP) phase 4, with a protocol to perform transient simulations of the last deglaciation (19–11 ka; although the protocol covers 26–0 ka). Similar to the last deglaciation, the disintegration of continental ice-sheets during the penultimate deglaciation led to significant changes in the oceanic circulation during Heinrich Stadial 11 (~ 136–129 ka). However, the two deglaciations bear significant differences in magnitude and temporal evolution of climate and environmental changes. Here, as part of the PAGES-PMIP working group on Quaternary Interglacials, we propose a protocol to perform transient simulations of the penultimate deglaciation under the auspices of PMIP4. This design includes time-varying changes in orbital forcing, greenhouse gas concentrations, continental ice-sheets as well as freshwater input from the disintegration of continental ice-sheets. This experiment is designed for AOGCMs to assess the coupled response of the climate system to all forcings. Additional sensitivity experiments are proposed to evaluate the response to each forcing. Finally, a selection of paleo records representing different parts of the climate system is presented, providing an appropriate benchmark for upcoming model-data comparisons across the penultimate deglaciation.

We have coupled an Earth system model of intermediate complexity (LOVECLIM) to the Glacial Systems Model (GSM) using the LCice 1.0 coupler. The coupling scheme is flexible enough to enable asynchronous coupling between any glacial cycle ice sheet model and (with some code work) any Earth system model of intermediate complexity (EMIC). This coupling includes a number of interactions between ice sheets and climate that are often neglected: dynamic meltwater runoff routing, novel downscaling for precipitation that corrects orographic forcing to the higher resolution ice sheet grid (“advective precipitation”), dynamic vertical temperature gradient, and ocean temperatures for sub-shelf melt. The sensitivity of the coupled model with respect to the selected parameterizations and coupling schemes is investigated. Each new coupling feature is shown to have a significant impact on ice sheet evolution. An ensemble of runs is used to explore the behaviour of the coupled model over a set of 2000 parameter vectors using present-day (PD) initial and boundary conditions. The ensemble of coupled model runs is compared against PD reanalysis data for atmosphere (2 m temperature, precipitation, jet stream, and Rossby number of jet), ocean (sea ice and Atlantic Meridional Overturning Circulation – AMOC), and Northern Hemisphere ice sheet thickness and extent. The parameter vectors are then narrowed by rejecting model runs (1700 CE to present) with regional land ice volume changes beyond an acceptance range. The selected subset forms the basis for ongoing work to explore the spatial–temporal phase space of the last two glacial cycles.

In Fennoscandia, tectonics, Glacial Isostatic Adjustment (GIA), and climatic changes cause ongoing crustal deformation of some millimetres per year, both vertically and horizontally. These displacements of the Earth can be measured to a high degree of precision using a Global Navigation Satellite System (GNSS). Since about three decades, this is the major goal of the Baseline Inferences for Fennoscandian Rebound, Sea-level, and Tectonics (BIFROST) project.We present a new velocity field for an extended BIFROST GNSS network in the ITRF2008 reference frame making use of the GNSS processing package GPS Analysis Software of MIT (GAMIT). Compared to earlier publications, we have almost doubled the number of stations in our analysis and increased the observation time span, thereby avoiding the early years of the network with many instrument changes. We also provide modelled vertical deformation rates from contributing processes, i.e. elastic deformation due to global atmospheric and non-tidal ocean loading, ice mass and hydrological changes as well as GIA. These values for the vertical component can be used for removal of these contributions so that the residual uplift signal can be further analysed, e.g., in the context of local or regional deformation processes or large-scale but low-magnitude geodynamics.The velocity field has an uplift maximum of 10.3 mm/yr in northern Sweden west of the Gulf of Bothnia and subsidence exceeding 1 mm/yr in northern Central Europe. The horizontal velocity field is dominated by plate motion of more than 20.0 mm/yr from south-west to north-east. The elastic uplift signal sums up to 0.7 - 0.8 mm/yr for most stations in Northern Europe. Hence, the maximum uplift related to the past glaciation is ca. 9.6 mm/yr. The residual uplift signal after removal of the elastic and GIA contribution may point to possible improvements of the GIA model, but may also indicate regional tectonic and erosional processes as well as local deformation effects. We show an example of such residual signal discussing potential areas of interest for further studies.

Borehole temperature profiles provide robust estimates of past ground surface temperature changes, in agreement with meteorological data. Nevertheless, past climatic changes such as the Last Glacial Cycle (LGC) generated thermal effects in the subsurface that affect estimates of recent climatic change from geothermal data. We use an ensemble of ice sheet simulations spanning the last 120ka to assess the impact of the Laurentide ice sheet on recent ground surface temperature histories reconstructed from borehole temperature profiles over North America. When the thermal remnants of the LGC are removed, we find larger amounts of subsurface heat storage [2.8 times] and an increased warming of the ground surface over North America by 0.75K, both relative to uncorrected borehole estimates.

We determine the contribution of glacial isostatic adjustment (GIA) to future relative sea-level change for the North American coastline between Newfoundland and Texas. We infer GIA model parameters using recently compiled and quality-assessed databases of past sea-level changes, including new databases for the United States Gulf Coast and Atlantic Canada. At 13 cities along this coastline, we estimate the GIA contribution to range from a few centimeters (e.g., 3 [−1 to 9] cm Miami) to a few decimeters (e.g., 18 [12–22] cm, Halifax) for the period 2085–2100 relative to 2006–2015 (1−σ ranges given). We provide estimates of uncertainty in the GIA component using two different methods; the more conservative approach produces total ranges (1−σ confidence) that vary from 3 to 16 cm for the cities considered. Contributions from ocean steric and dynamic changes as well as those from changes in land ice are also estimated to provide context for the GIA projections. When summing the contributions from all three processes at the 13 cities considered along this coastline, using median or best-estimate values, the GIA signal comprises 5–38% of the total depending on the adopted climate forcing and location. The contributions from ocean dynamic/steric changes and ice mass loss are similar in amplitude but with spatial variation that approximately cancels, resulting in GIA dominating the net spatial variability north of 35°N.

The contribution of the Greenland and West Antarctic ice sheets to sea level has increased in recent decades, largely owing to the thinning and retreat of outlet glaciers and ice streams1, 2, 3, 4. This dynamic loss is a serious concern, with some modelling studies suggesting that the collapse of a major ice sheet could be imminent5, 6 or potentially underway7 in West Antarctica, but others predicting a more limited response8. A major problem is that observations used to initialize and calibrate models typically span only a few decades, and, at the ice-sheet scale, it is unclear how the entire drainage network of ice streams evolves over longer timescales. This represents one of the largest sources of uncertainty when predicting the contributions of ice sheets to sea-level rise8, 9, 10. A key question is whether ice streams might increase and sustain rates of mass loss over centuries or millennia, beyond those expected for a given ocean–climate forcing5, 6, 7, 8, 9, 10. Here we reconstruct the activity of 117 ice streams that operated at various times during deglaciation of the Laurentide Ice Sheet (from about 22,000 to 7,000 years ago) and show that as they activated and deactivated in different locations, their overall number decreased, they occupied a progressively smaller percentage of the ice sheet perimeter and their total discharge decreased. The underlying geology and topography clearly influenced ice stream activity, but—at the ice-sheet scale—their drainage network adjusted and was linked to changes in ice sheet volume. It is unclear whether these findings can be directly translated to modern ice sheets. However, contrary to the view that sees ice streams as unstable entities that can accelerate ice-sheet deglaciation, we conclude that ice streams exerted progressively less influence on ice sheet mass balance during the retreat of the Laurentide Ice Sheet.

Understanding changes in the climate of the late Pleistocene and the Holocene has long been a research topic. Studies rely on different sources of information, ranging from terrestrial and marine archives to a hierarchy of climate modeling activities. In contrast to the climate of the last millennium, novel approaches are necessary to bridge the different temporal and spatial representations of the various archives and the climate models, and to achieve a robust understanding of climate variability and climate processes on centennial-to-millennial time scales.

The penultimate deglaciation (PDG, ∼138–128 thousand years before present, hereafter ka) is the transition from the penultimate glacial maximum (PGM) to the Last Interglacial (LIG, ∼129–116 ka). The LIG stands out as one of the warmest interglacials of the last 800 000 years (hereafter kyr), with high-latitude temperature warmer than today and global sea level likely higher by at least 6 m. Considering the transient nature of the Earth system, the LIG climate and ice-sheet evolution were certainly influenced by the changes occurring during the penultimate deglaciation. It is thus important to investigate, with coupled atmosphere–ocean general circulation models (AOGCMs), the climate and environmental response to the large changes in boundary conditions (i.e. orbital configuration, atmospheric greenhouse gas concentrations, ice-sheet geometry and associated meltwater fluxes) occurring during the penultimate deglaciation. A deglaciation working group has recently been set up as part of the Paleoclimate Modelling Intercomparison Project (PMIP) phase 4, with a protocol to perform transient simulations of the last deglaciation (19–11 ka; although the protocol covers 26–0 ka). Similar to the last deglaciation, the disintegration of continental ice sheets during the penultimate deglaciation led to significant changes in the oceanic circulation during Heinrich Stadial 11 (∼136–129 ka). However, the two deglaciations bear significant differences in magnitude and temporal evolution of climate and environmental changes. Here, as part of the Past Global Changes (PAGES)-PMIP working group on Quaternary interglacials (QUIGS), we propose a protocol to perform transient simulations of the penultimate deglaciation under the auspices of PMIP4. This design includes time-varying changes in orbital forcing, greenhouse gas concentrations, continental ice sheets as well as freshwater input from the disintegration of continental ice sheets. This experiment is designed for AOGCMs to assess the coupled response of the climate system to all forcings. Additional sensitivity experiments are proposed to evaluate the response to each forcing. Finally, a selection of paleo-records representing different parts of the climate system is presented, providing an appropriate benchmark for upcoming model–data comparisons across the penultimate deglaciation.

Abrupt climate shifts of large amplitude were common features of the Earth's climate as it transitioned into and out of the last full glacial state approximately twenty thousand years ago, but their causes are not yet established. Mid-latitudinal atmospheric dynamics may have played an important role in these oscillations through their effects on heat and precipitation distributions, sea ice extent, and wind-driven ocean circulation patterns. This study characterises deglacial winter wind changes over the North Atlantic (NAtl) in a suite of transient deglacial simulations we performed using the PlaSim earth system model, as well as in the TraCE-21ka simulation. We detect multiple instances of NAtl jet transitions that occur within 10 simulation years and a sensitivity of the jet to background climate conditions. Thus, we suggest that changes to the NAtl jet may play a critical role in abrupt glacial climate oscillations. We identify two types of simulated wind changes over the last deglaciation. Firstly, the latitude of the NAtl eddy-driven jet shifts northward over the deglaciation in a sequence of distinct steps. Secondly, the variability of the NAtl jet gradually shifts from a Last Glacial Maximum (LGM) state with a strongly preferred jet latitude and a restricted latitudinal range to one with no single preferred latitude and a range that is at least 11° broader. Changes to the position of the NAtl jet alter the location of the wind forcing driving oceanic surface gyres and the limits of sea ice extent, whereas a shift to a more variable jet reduces the effectiveness of the wind forcing at driving surface ocean transports. The processes controlling these two types of changes differ on the upstream and downstream ends of the NAtl eddy-driven jet. On the upstream side over eastern North America, the elevated ice sheet margin acts as a physical barrier to the winds in both the PlaSim simulations and the TraCE-21ka experiment. This constrains both the position and the latitudinal variability of the jet at LGM, so the jet shifts in sync with ice sheet margin changes. In contrast, the downstream side over the eastern NAtl is more sensitive to the thermal state of the background climate. Our results suggest that knowing the position of the south-eastern margin of the North American ice complex strongly constrains the deglacial position of the jet over eastern North America and the western North Atlantic as well as its variability.

The penultimate deglaciation (~ 138–128 thousand years before present, hereafter ka) is the transition from the penultimate glacial maximum to the Last Interglacial (LIG, ~ 129–116 ka). The LIG stands out as one of the warmest interglacials of the last 800 ka, with high-latitude temperature warmer than today and global sea level likely higher by at least 6 meters. The LIG therefore receives ever-growing attention, in particular to identify mechanisms and feedbacks responsible for such regional warmth that is comparable to that expected before 2100. Considering the transient nature of the Earth system, the LIG climate and ice-sheets evolution were certainly influenced by the changes occurring during the penultimate deglaciation. It is thus important to investigate the climate and environmental response to the large changes in boundary conditions (i.e. orbital configuration, atmospheric greenhouse gas concentrations, ice sheet geometry) occurring during this time interval. A deglaciation working group has recently been set up as part of the Paleoclimate Modelling Intercomparison Project (PMIP) phase 4, with a protocol to perform transient simulations of the last deglaciation (19–11 ka). Similar to the last deglaciation, the disintegration of continental ice-sheets during the penultimate deglaciation led to significant changes in the oceanic circulation during Heinrich Stadial 11 (~ 136–129 ka). However, the two deglaciations bear significant differences in magnitude and temporal evolution of climate and environmental changes. Here, as part of the PAGES-PMIP working group on Quaternary Interglacials, we propose a protocol to perform transient simulations of the penultimate deglaciation to complement the PMIP4 effort. This design includes time-varying changes in orbital forcing, greenhouse gas concentrations, continental ice-sheets as well as freshwater input from the disintegration of continental ice-sheets. This experiment is designed to assess the coupled response of the climate system to all forcings. Additional sensitivity experiments are proposed to evaluate the response to each forcing. Finally, a selection of paleo records representing different parts of the climate system is presented, providing an appropriate benchmark for upcoming model-data comparisons across the penultimate deglaciation.