Search

Paleo records indicate significant variation in sea level and temperature proxies between different glacial cycles. What is unclear is the extent to which these differences are due to noise in the physical system versus a robust response to external forcings. When one considers what is happening with each individual ice sheet, variations between glacial cycles are largely unknown, given the few relevant records available to constrain ice sheet extent before the Eemian. To explore both the controls on past ice sheet and climate evolution and explore bounds on what the evolution might actually have looked like, we are running ensemble simulations of the last two glacial cycles with the fully coupled ice/climate model LCice. LCice is a coupled version of the Loveclim EMIC and GSM glacial systems model with hybrid shallow shelf and shallow ice flow and global visco-elastic glacio-isostatic adjustment. The current configuration includes all 4 ice sheet complexes and is subject to only orbital and greenhouse gas forcing.To answer the above questions, we present ensemble results for the last two glacial inceptions, focusing on what key ice sheet and climate characteristics are robust across the ensemble and what are not. The role of key forcings and feedbacks are also isolated through a set of sensitivity experiments.

The mid-Pleistocene Transition (MPT) from 41 kyr to 100 kyr glacial cycles was one of the largest changes in the Earth system over the past 2 million years. The transition happened in the absence of a relevant change in orbital forcing. As such, it presents a challenge for the Milankovitch theory of glacial cycles. A change from a low to high friction bed under the North American Ice Complex through the removal of pre-glacial regolith has been hypothesized to play a critical role in the transition. For testing, this hypothesis requires constraint on pre-glacial regolith cover and topography as well as mechanistic constraint on whether the appropriate amount of regolith can be removed from the required regions to enable MPT occurrence at the right time. To date, however, Pleistocene regolith removal has not been simulated for a realistic, 3D North American ice sheet fully resolving relevant basal processes. A further challenge is very limited constraints on pre-glacial bed elevation and sediment thickness.Herein, we address these challenges with an appropriate computational model and ensemble-based analysis addressing parametric and initial mean sediment cover uncertainties. We use the 3D Glacial Systems Model that incorporates relevant glacial processes. Specifically, it includes: 3D thermomechanically coupled hybrid SIA/SSA ice physics, fully coupled sediment production and transport, subglacial linked-cavity and tunnel hydrology, isostatic adjustment from dynamic loading and erosion, and climate from a 2D non-linear energy balance model and glacial index. The sediment model includes quarrying and abrasion for sediment production with both englacial and subglacial transport. The coupled system is driven only by atmospheric CO2 and insolation.We show that the ice, climate, and sediment processes encapsulated in this fully coupled glacial systems model enables capture of the evolution of the Pleistocene North American glacial system. Specifically and within observational uncertainty, our model captures: the shift from 41 to 100 kyr glacial cycles, early Pleistocene extent, LGM ice volume, deglacial ice extent, and the broad present-day sediment distribution. We also find that pre-glacial sediment thickness and topography have a strong influence on the strength and duration of early Pleistocene glaciations.

The interglacial period spanning ca. 423 to 398 ka and known as Marine Isotope Stage (MIS) 11c has been the subject of much study, due largely to the unique evolution of global temperatures, greenhouse gas levels, and sea levels relative to other interglacials of the late Pleistocene. Particularly concerning is some geological evidence and prior modeling studies which have suggested that a large majority of the Greenland ice sheet (GrIS) disappeared during this period, despite global mean air temperatures only modestly higher than those of the preindustrial period. However, uncertainty is high as to the extent and spatiotemporal evolution of this melt due to a dearth of direct geological constraints. Our study therefore endeavors to better constrain these large uncertainties by using spatiotemporally interpolated climate forcing from CESM v1.2 time slice simulations and an ensemble of ice sheet model parameter vectors derived from a GrIS history matching over the most recent glacial cycle from the Glacial Systems Model (GSM). The use of different ice sheet initialization states from simulations of the previous glacial-interglacial transition helps to capture the large initial condition uncertainty. Two different regional present-day climate modeling datasets are utilized for anomaly correction of CESM precipitation and temperature fields. Preliminary analysis indicates that the most robust retreat across most ensemble members happens in the northern, western, and central portions of the ice sheet, while the higher terrain of the south and east retain substantial amounts of ice. This is broadly consistent with indications that ice may have survived the MIS-11c interglacial at the Summit ice core location, but not at DYE-3. Simulations indicate a maximum MIS-11c sea level contribution from the GrIS centered between 408 and 403 ka, with minimum GrIS volumes reaching between 25% and 70% of modern-day values. In part due to the prior constraint of ice-sheet model ensemble parameters from history matching, ensemble parameters controlling downscaling and climate forcing bias correction are the largest parametric sources of output variance in our simulations. Though CESM uncertainties are unassessed in this study, it is likely they dominate given that the choice of present-day reference temperature climatology for anomaly correction of the climate model output has the largest effect on the GrIS melt response in our simulations.

Some ice sheets and glaciers experience long quiescent periods interspersed with short periods of rapid ice advance, such as the binge-purge-type cycling hypothesized to be associated with Heinrich Events. Modeling ice stream activation/de-activation, however, is numerically challenging given the relatively abrupt changes at surge onset and the high ice velocities. In spite of this, a number of high-profile modeling papers have explored Heinrich events and ice surges, but generally with very limited consideration of numerical aspects. Here we test the ability of the 3D Glacial Systems Model (GSM) and Parallel Ice Sheet Model (PISM) to simulate binge-purge-type surges and explore the stability of the simulations with respect to relevant numerical and discretization uncertainties. We find surge characteristics exhibit a resolution dependency but converge at higher horizontal grid resolutions (order 5 km). In accordance with theoretical and experimental work, our model results suggest that the thermal activation of basal sliding should start below the pressure melting point. A resolution-dependent basal temperature ramp for the thermal activation of basal sliding as well as a subglacial hydrology model can reduce the discrepancies between high and coarse horizontal grid resolutions. Furthermore, incorporating a bed thermal and at least a minimal complexity subglacial hydrology model significantly affects surge characteristics and is, therefore, essential for modeling large-scale ice stream cycling.

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.

A leading contender for explaining the mid-Pleistocene transition (MPT) from small 40 kyr glaciations to large, abruptly terminating 100 kyr ones is a shift to high friction bed under the Northern hemisphere ice sheets – the North American ice sheet in particular. The regolith hypothesis posits that this occurred with the removal of deformable regolith – laying bare higher-friction bedrock under ice sheet core domains. Is the regolith hypothesis consistent with the physics of glacial removal of mechanically weak surface material? Self-consistency of the regolith hypothesis has not been tested for a realistic, 3D North American ice sheet, capturing the transition from soft to hard bedded and 40 to 100 kyr cycles, fully considering basal processes and sediment production. To test self-consistency, we simulate the pace and distribution of regolith removal in a numerical ice sheet model incorporating the relevant glacial processes and their uncertainties. Specifically, the Glacial Systems Model includes: fully coupled sediment production and transport, subglacial hydrology, glacial isostatic adjustment, 3D thermomechanically coupled hybrid ice physics, and internal climate solution from a 2D non-linear energy balance model. The sediment model produces sediment via quarrying and abrasion while transporting material englacially and subglacially. The subglacial hydrology model employs a linked-cavity system with a flux based switch to tunnel drainage, giving dynamic effective pressure needed for realistic sediment and sliding processes. Deflection and rebound of the Earth's surface are calculated for a range of solid Earth visco-elastic rheologies. The coupled system is driven only by prescribed atmospheric CO2 and orbitally derived insolation. Starting from a range of initial sediment distributions and simulating an ensemble of model parameter values, we model the rate and spatial distribution of regolith dispersal and compare this against the inferred range of Pliocene regolith thickness, the present day sediment distribution, and the timing of the MPT. A first order fully coupled representation of ice, climate and sediment interactions captures the transition within parametric and observational uncertainty. The system gives the shift from 40 to 100 kyr glacial cycles while broadly reproducing the present day sediment distribution, inferred early Pleistocene extent, LGM ice volume and deglacial margin locations.

New research aims to improve relative sea-level (RSL) projections for the Norwegian coast. The main objectives are to: i) collect observations of past RSL changes, ranging from the end of the last ice age to the last century, ii) develop a high-quality database of post-glacial sea-level index points (SLIPs) for the Norwegian coast, and to iii) improve our understanding of past and future vertical land motion using glacial isostatic adjustment (GIA) modelling. To now, our collection of new empirical data has focussed on three significant, but enigmatic RSL histories that are not adequately reproduced in existing GIA models: very recent stillstands and transgressions documented by historical tide gauge records, rapid transgressions during the early- to mid-Holocene Tapes period, and abrupt transgressions during the latest Pleistocene Younger Dryas chronozone. Ongoing field sampling is focussed on developing high-resolution RSL trends from salt marshes, isolation basins, and raised beaches, using multiple biostratigraphic and geochemical proxies (i.e. micropaleontology, macrofossils, x-ray fluorescence, C/N) and dating techniques (i.e. Pb-210, Cs-137, C-14, tephrochronology, geochemical markers). Results from various localities spanning the Norwegian coast provide robust constraints for the timing and rate of RSL change during the Younger Dryas and Tapes chronozones. Additional results providing new estimates of very recent RSL trends in southwest Norway are presented by Holthuis et al. (Late Holocene sea-level change and storms in southwestern Norway based on new data from intertidal basins and salt marshes; Session CL5.2.2). These new and emerging constraints are being integrated into a post-glacial RSL database that incorporates high-quality data from the entire Norwegian coastline. Over 1000 SLIPs have been assembled from published studies. These existing data were updated using current radiocarbon calibration curves, high-resolution digital elevation models, new field observations, and new quantitative estimates of relevant uncertainties. Ongoing GIA modelling is utilizing the new RSL database, a glaciological model that freely simulates ice sheet changes, as well as geodetic and ice margin chronology constraints, to develop rigorous uncertainty estimates for present and future GIA along the Norwegian coast.

The lithospheric thickness beneath and around Greenland varies from a few tens of kilometres in offshore regions to several tens of kilometres (up to 200 – 250 km) in land areas. But, due to different datasets and techniques applied in geophysical studies, there are large differences between the different geophysical lithosphere models. As an example, lithosphere models from seismological datasets show generally larger values (above 100 km), while models using gravity or thermal datasets tend to be thinner (values mostly below 100 km). To model the deformation associated with the melting of the Greenland Ice Sheet a detailed lithosphere model is required. Nevertheless, seismologically obtained lithosphere models are the ones usually applied in these so-called glacial isostatic adjustment (GIA) models. Besides, GIA models can be used to provide additional constraints on the lithospheric thickness.Results from most 3D GIA models are compared to observed vertical velocities only, while horizontal velocities are known to be sensitive to the lateral variations of the Earth (e.g., lithospheric thickness). But, horizontal velocities from incompressible GIA models, which are commonly used, are not suitable due to the neglect of material parameter changes related to the dilatation. Compressible GIA models in turn can provide more accurate estimates of the horizontal and vertical viscoelastic deformations induced by ice-mass changes. Here, we use a variety of lithospheric thickness models, obtained from gravity, thermal, and seismological datasets, in a three-dimensional compressible GIA Earth model. The GIA model will be constructed using the finite-element software ABAQUS (Huang et al., under review in GJI) and applying recent ice history models Huy3 and GLAC-GR2a for Greenland in combination with the Little Ice Age deglaciation model by Kjeldsen et al. (2015). We will compare various lithosphere models, including their impact on the modelled 3D velocity field, and compare these against independent GNSS (Global Navigation Satellite System) observations.References:Huang, P., Steffen, R., Steffen, H., Klemann, V., van der Wal, W., Reusen, J., Wu, P., Tanaka, Y., Martinec, Z., Thomas, M. (under review in GJI): A finite element approach to modelling Glacial Isostatic Adjustment on three-dimensional compressible earth models. Geophysical Journal International. Under review.Kjeldsen, K., Korsgaard, N., Bjørk, A. et al. (2015): Spatial and temporal distribution of mass loss from the Greenland Ice Sheet since AD 1900. Nature 528, 396–400, https://doi.org/10.1038/nature16183.

Many coastal cities are an early casualty in climate-related coastal flooding because of processes resulting in land subsidence and thus enhanced relative sea-level (RSL) rise. Much of the Atlantic coast of North America has been sinking for thousands of years, at a maximum rate of ~20 cm per century as a consequence of solid Earth deformation in response to deglaciation of the Laurentide ice sheet (between ~18,000 and ~7,000 years ago) [e.g. Love et al., Earth's Future, 4(10), 2016]. Karegar et al. [Geophysical Research Letters, 43(7), 2016] have shown that vertical land motion along the Atlantic coast of the USA is an important control on nuisance flooding. A key finding in this study is that while glacial isostatic adjustment (GIA) is the dominant process driving land subsidence in most areas, there can be large deviations from this signal due to the influence of anthropogenic activity impacting hydrological processes. For example, between Maine (45°N) and New Hampshire (43°N), the GPS data show uplift while geological data show long-term subsidence. The cause of this discrepancy is not clear, but one hypothesis is increasing water mass associated with the James Bay Hydroelectric Project in Quebec [Karegar et al., Scientific Reports, 7, 2017].The primary aim of this study is to better constrain and understand the processes that contribute to contemporary and future vertical land motion in this region to produce improved projections of mean sea-level change and nuisance flooding. The first step towards achieving these aims is to determine a GIA model parameter set that is compatible with observations of past sea-level change for this region. We make use of two regional RSL data compilations: Engelhart and Horton [Quaternary Science Reviews, 54, 2012] for northern USA and Vacchi et al. [Quaternary Science Reviews, 201, 2018] for Eastern Canada, comprising a total of 1013 data points (i.e., sea level index points and limiting data points) over 38 regions distributed throughout our study region. These data are well suited to determine optimal GIA model parameters due to the magnitude of other signals being much smaller, particularly in near-field regions such as Eastern Canada. We consider a suite of 32 ice history models that is comprised mainly of a subset from Tarasov et al. [Earth and Planetary Science Letters, 315–316, 2012] as well as the ICE-6G and ANU models. We have computed RSL for these ice histories using a state-of-the-art sea-level calculator and 440 1-D Earth viscosity models per each ice history model to identify a set of Earth model parameters that is compatible with the observations.

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, eg. the Younger Dryas and Dansgaard Oeschger events. Freshwater injection, or hosing, model experiments demonstrate that freshwater has the capability to generate abrupt climate transitions. However, in an attempt to mitigate the inability of most models to resolve the smaller-scale features relevant to freshwater transport (such as boundary currents and mesoscale eddies), these hosing experiments commonly apply the entirety of the freshwater directly to the regions of deepwater formation (DWF). Our results indicate that this can inflate the freshwater signal in those regions by as much as four times. We propose a novel method of freshwater injection for such low-resolution models that spatially distributes the freshwater in accord with the results of eddy-permitting modelling. Furthermore, this “freshwater fingerprint” method not only impacts the timing of simulated climate transitions but also can allow us to evaluate how much we are overestimating the effects of freshwater when injected directly into sites of DWF. The freshwater fingerprints we develop are based on a suite of freshwater injection experiments performed using an eddy permitting Younger Dryas configuration of the MITGCM. Freshwater injection locations include the Mackenzie River, Gulf of St. Lawrence, Gulf of Mexico and a location off the coast of Norway, with flux amounts bounded by glacial reconstructions. These simulations indicate that freshwater from the Mackenzie River and Fennoscandia have the largest impact on salinity in most of the conventional sites of DWF (GIN and Labrador Seas, and in these simulations, predominantly the northern North Atlantic due to extensive sea ice), while freshwater from the Gulf of St. Lawrence is effective at freshening only the northern North Atlantic. The Gulf of Mexico has little impact on any DWF region we consider, mostly because the lower but continual flux in our simulations does not allow freshwater to penetrate northward past the Gulf Stream. The dilution of the freshwater signal as it is transported from the site of injection to the DWF zones leads to a reduction in the effective freshwater forcing, making hosing directly over DWF zones even with realistic freshwater amounts unrealistic. Thus, we construct freshwater fingerprints from these simulations by extracting the freshwater anomaly spatial pattern averaged over the last 5 simulation years, vertically integrating the field and normalizing it.The freshwater fingerprint is then implemented in the COSMOS Earth Systems Model, which is run at resolutions typical for paleoclimate simulations (non-eddy permitting). Initial results show that freshwater from the Mackenzie River using our fingerprint method leads to a more gradual cooling than if the meltwater is released directly over the hosing region (50-70N). We conclude that hosing over DWF zones, even with realistic freshwater amounts, produces an unrealistically large climate response. Additional results for the remaining injection locations and with the fingerprint implemented in a simpler climate model will be presented.

Is the regolith hypothesis consistent with the physics of glacial removal of mechanically weak surface material? The mid-Pleistocene transition (MPT) from small 40 kyr glacial cycles to large, abruptly terminating 100 kyr ones represents a major climate system reorganization for which a clear understanding is lacking. A leading mechanism for this transition is a stabilization of ice sheets due to a shift to higher friction substrate. The Pleistocene saw the removal of deformable regolith -- laying bare hard higher-friction bedrock that would help preserve regional ice during warm interstadials. This is the regolith hypothesis. The removal of regolith by Pleistocene ice sheets remains poorly constrained. To date, only models with a forced change in area of regolith cover or 1D flow line models with simplistic sediment transport have been used to probe the role of regolith in the MPT. It is therefore unclear if the appropriate amount of regolith removal can occur within the time-frame of the MPT. To properly test the hypothesis, at least three components are required: capable model, observational constraint, and a probe of uncertainties. A capable model must explicitly represent relevant processes in a fully coupled self-consistent manner. We have therefore configured a state of the art 3D glacial systems model (GSM). The GSM incorporates a state-of-the-art fully coupled sediment production/transport model, subglacial hydrology, visco-elastic glacial isostatic adjustment, 3D thermomechanically coupled hybrid shallow ice/shallow shelf ice dynamics, and internal climate solution from an energy balance model. The model generates sediment by quarrying and abrasion, and both subglacial and englacial sediment transport. The subglacial hydrology model employs a linked-cavity system with a flux based switch to tunnel drainage, giving dynamic effective pressure needed for realistic sediment and sliding processes. The coupled model is driven only by prescribed atmospheric CO2 and orbitally derived insolation. The required observational constraints include present-day regolith distribution and inferred Pleistocene ice volume from proxy records. The final component is a large ensemble of full Pleistocene simulations that probe both initial regolith distribution uncertainties and model parametric uncertainties. We present the results of such an ensemble, examining both rates of computed regolith removal and changes in ice volume cycling across the MPT interval.

One of the main contributors to palaeoclimate variability on millennial timescales is understood to be Dansgaard-Oeschger (D-O) cycles. Our awareness of these phenomena arises primarily from quasi-periodic, abrupt transitions of large magnitude detected in δ18O records from Greenland ice cores (e.g. Dansgaard et al, 1982; Johnsen et al, 1992), although there is evidence of similar variability in other archives and regions. D-O cycles have plenty to capture the imagination: the strength and rapidity of climate changes over Greenland, their regularity throughout MIS3 (~60 to 30 thousand years before present) and occurrence during the last deglaciation contrasting with their relative absence during the Last Glacial Maximum and Holocene, their opposed characteristics in Greenland and Antarctica, and that different models require different boundary conditions to reproduce this phenomena, if they can reproduce it at all. This talk characterises D-Olike cycles in two different models: Planet Simulator (PlaSim, an Earth System Model with simplified atmospheric physics, thermodynamic sea ice, and simplified ocean dynamics), and COSMOS (a CMIP3-era ESM). We identify four phases to D-O cycles and commonalities and differences in their representations in these models. Finally, we examine which phases of this type of variability continue to contribute to climate variability today and what that looks like.

The Pacific Coast of Central North America is a geodynamically complex region subject to various geophysical processes with different patterns of vertical land motion. It includes two distinct tectonic regimes: the Cascadia subduction zone and the strike-slip San Andreas fault system. The vertical land motion in this region reflects not only tectonic activity of these plate boundaries, but also isostatic signals associated with different loading effects such as the (de)glaciation of North American ice sheets and the more contemporary, anthropogenically-related groundwater extraction and mountain glacier mass loss. These processes occur over a broad range of timescales and are observed by a variety of measuring techniques.Here we combine geological measurements of relative sea level (RSL) change with contemporary observations of vertical land motion inferred from geodetic data to decipher and thus better understand the contribution from various individual processes. Our results suggest that contemporary vertical land motion is dominated by Cascadia interseismic deformation and the isostatic response to the retreat of the North American ice sheets but is also influenced by other contemporary processes. We present some model results that illustrate the contributions of the above-mentioned processes to RSL projections along this coastline.

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"]}

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.This work has been published in Steffen et al. (2019).Reference:Steffen, H., Steffen R., Tarasov L. 2019. Modelling of glacially-induced stress changes in Latvia, Lithuania and the Kaliningrad District of Russia. Baltica, 32 (1), 78–90.

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.

What determines the character of glacial inceptions? Does the spatio-temporal pattern of ice nucleation and expansion vary much between Late Pleistocene glacial inceptions? According to various benthic del18O stacks, the MIS 7 interglacial was the most anomalous in character of the last 4 interglacials. Key differences include a weaker interglacial state and an initial fast inception interrupted by a return to a similar and extended interglacial state. These anomalies of MIS 7 along with temporal proximity arguably make the last two glacial inceptions the best test case for addressing our opening questions. As part of a larger project to generate and analyze a data-constrained ensemble of fully coupled ice/climate transient simulations for the last two complete glacial cycles, herein we present initial results comparing the last two glacial inceptions (MIS 7 and 5d). We are using a new version of the fully coupled ice/climate model LCice. LCice now simulates all 4 paleo ice sheet complexes with hybrid shallow-shelf and shallow-ice physics. It has already been shown to capture northern hemispheric ice sheet growth and subsequent retreat consistent with inferences from global mean sea level proxies (Bahadory et al, 2019). Orbital and greenhouse gas changes are the only external forcings applied to the model. A 300 member ensemble probes parametric uncertainties in both the 3D Glacial Systems Model and LoveClim (Atmosphere/Ocean/Vegetation) components of LCice. Our presentation will compare the evolution and relative phasing of all 4 paleo ice sheets, and associated changes in the rest of the modelled climate system.

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.

We infer the GIA signal and its uncertainty along the central Pacific coast of North America using 680 sea-level index points and over 20,000 model runs sampling >700 (1-D) Earth viscosity models and 29 ice sheet reconstructions. Due to the large spatial extent and different tectonic settings of the study area, we divided it into three sub-regions (northern, central and southern) for which model parameters were inferred separately. Also, given that this region is tectonically active, the influence of this process (as well as sediment isostatic adjustment) was accounted for where possible by removing it from the data using published estimates. Our results indicate that it is not possible to produce an acceptable to for all of the RSL data with a single set of model parameters, suggesting significant lateral variability in viscous structure. Specifically, low viscosities (1019–1020 Pas) are inferred in the upper mantle within the northern region (southwestern British Columbia and northwest Washington) compared to those inferred (2–5 × 1020 Pas) for the central and southern regions (extending from southern Washington to southern California). High quality model fits were obtained for all data except those from the northern region where no single parameter set was able to capture both the rapid and large RSL fall during the late glacial and the monotonic rise during the mid-to-late Holocene at all localities. This suggests the need for an Earth model that incorporates departures from a linear Maxwell rheology (as applied here) and/or lateral variations in viscosity structure. Using our optimal model parameter sets, we show that GIA is a significant contributor to both contemporary vertical land motion and relative sea level change in our study region and so should be considered when interpreting observations of these signals and for making future relative sea level projections. Model output of present-day vertical land motion at 483 GPS stations and sea-level change at 56 tide gauge stations is provided (with estimated uncertainty) so that these data can be used to study non-GIA processes more accurately.

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.