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.

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.

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.