Benchmarking the vertically integrated ice-sheet model IMAU-ICE (version 2.0)

Ice-dynamical processes constitute a large uncertainty in future projections of sea-level rise caused by anthropogenic climate change. Improving our understanding of these processes requires ice-sheet models that perform well at simulating both past and future ice-sheet evolution. Here, we present version 2.0 of the ice-sheet model IMAU-ICE, which uses the depth-integrated viscosity approximation (DIVA) to solve the stress balance. We evaluate its performance in a range of benchmark experiments, including simple analytical solutions and both schematic and realistic model intercomparison exercises. IMAU-ICE has adopted recent developments in the numerical treatment of englacial stress and sub-shelf melt near the grounding line, which result in good performance in experiments concerning grounding-line migration (MISMIP, MISMIP+) and buttressing (ABUMIP). This makes it a model that is robust, versatile, and user-friendly, which will provide a firm basis for (palaeo-)glaciological research in the coming years.publishedVersio

Miocene Antarctic Ice Sheet area adapts significantly faster than volume to CO2-induced climate change

The strongly varying benthic δ18O levels of the early and mid-Miocene (23 to 14 Myr ago) are primarily caused by a combination of changes in Antarctic Ice Sheet (AIS) volume and deep-ocean temperatures. These factors are coupled since AIS changes affect deep-ocean temperatures. It has recently been argued that this is due to changes in ice sheet area rather than volume because area changes affect the surface albedo. This finding would be important when the transient AIS grows relatively faster in extent than in thickness, which we test here. We analyse simulations of Miocene AIS variability carried out using the three-dimensional ice sheet model IMAU-ICE forced by warm (high CO2, no ice) and cold (low CO2, large East AIS) climate snapshots. These simulations comprise equilibrium and idealized quasiorbital transient runs with strongly varying CO2 levels (280 to 840 ppm). Our simulations show a limited direct effect of East AIS changes on Miocene orbital-timescale benthic δ18O variability because of the slow build-up of volume. However, we find that relative to the equilibrium ice sheet size, the AIS area adapts significantly faster and more strongly than volume to the applied forcing variability. Consequently, during certain intervals the ice sheet is receding at the margins, while ice is still building up in the interior. That means the AIS does not adapt to a changing equilibrium size at the same rate or with the same sign everywhere. Our results indicate that the Miocene Antarctic Ice Sheet affects deep-ocean temperatures more than its volume suggests

Modelling feedbacks between the Northern Hemisphere ice sheets and climate during the last glacial cycle

During the last glacial cycle (LGC), ice sheets covered large parts of Eurasia and North America, which resulted in ∼120 m of sea level change. Ice sheet-climate interactions have considerable influence on temperature and precipitation patterns and therefore need to be included when simulating this time period. Ideally, ice sheet-climate interactions are simulated by a high-resolution Earth system model. While these models are capable of simulating climates at a certain point in time, such as the pre-industrial (PI) or the Last Glacial Maximum (LGM; 21 000 years ago), a full transient glacial cycle is currently computationally unfeasible as it requires a too-large amount of computation time. Nevertheless, ice sheet models require forcing that captures the gradual change in climate over time to calculate the accumulation and melt of ice and its effect on ice sheet extent and volume changes. Here we simulate the LGC using an ice sheet model forced by LGM and PI climates. The gradual change in climate is modelled by transiently interpolating between pre-calculated results from a climate model for the LGM and the PI. To assess the influence of ice sheet-climate interactions, we use two different interpolation methods: the climate matrix method, which includes a temperature-albedo and precipitation-topography feedback, and the glacial index method, which does not. To investigate the sensitivity of the results to the prescribed climate forcing, we use the output of several models that are part of the Paleoclimate Modelling Intercomparison Project Phase III (PMIP3). In these simulations, ice volume is prescribed, and the climate is reconstructed with a general circulation model (GCM). Here we test those models by using their climate to drive an ice sheet model over the LGC. We find that the ice volume differences caused by the climate forcing exceed the differences caused by the interpolation method. Some GCMs produced unrealistic LGM volumes, and only four resulted in reasonable ice sheets, with LGM Northern Hemisphere sea level contribution ranging between 74-113 m with respect to the present day. The glacial index and climate matrix methods result in similar ice volumes at the LGM but yield a different ice evolution with different ice domes during the inception phase of the glacial cycle and different sea level rates during the deglaciation phase. The temperature-albedo feedback is the main cause of differences between the glacial index and climate matrix methods

Compensating errors in inversions for subglacial bed roughness: same steady state, different dynamic response

Subglacial bed roughness is one of the main factors controlling the rate of future Antarctic ice-sheet retreat and also one of the most uncertain. A common technique to constrain the bed roughness using ice-sheet models is basal inversion, tuning the roughness to reproduce the observed present-day ice-sheet geometry and/or surface velocity. However, many other factors affecting ice-sheet evolution, such as the englacial temperature and viscosity, the surface and basal mass balance, and the subglacial topography, also contain substantial uncertainties. Using a basal inversion technique intrinsically causes any errors in these other quantities to lead to compensating errors in the inverted bed roughness. Using a set of idealised-geometry experiments, we quantify these compensating errors and investigate their effect on the dynamic response of the ice sheet to a prescribed forcing. We find that relatively small errors in ice viscosity and subglacial topography require substantial compensating errors in the bed roughness in order to produce the same steady-state ice sheet, obscuring the realistic spatial variability in the bed roughness. When subjected to a retreat-inducing forcing, we find that these different parameter combinations, which per definition of the inversion procedure result in the same steady-state geometry, lead to a rate of ice volume loss that can differ by as much as a factor of 2. This implies that ice-sheet models that use basal inversion to initialise their model state can still display a substantial model bias despite having an initial state which is close to the observations

Benchmarking the vertically integrated ice-sheet model IMAU-ICE (version 2.0)

Ice-dynamical processes constitute a large uncertainty in future projections of sea-level rise caused by anthropogenic climate change. Improving our understanding of these processes requires ice-sheet models that perform well at simulating both past and future ice-sheet evolution. Here, we present version 2.0 of the ice-sheet model IMAU-ICE, which uses the depth-integrated viscosity approximation (DIVA) to solve the stress balance. We evaluate its performance in a range of benchmark experiments, including simple analytical solutions and both schematic and realistic model intercomparison exercises. IMAU-ICE has adopted recent developments in the numerical treatment of englacial stress and sub-shelf melt near the grounding line, which result in good performance in experiments concerning grounding-line migration (MISMIP, MISMIP+) and buttressing (ABUMIP). This makes it a model that is robust, versatile, and user-friendly, which will provide a firm basis for (palaeo-)glaciological research in the coming years

A computationally efficient depression-filling algorithm for digital elevation models, applied to proglacial lake drainage

Many processes govern the deglaciation of ice sheets. One of the processes that is usually ignored is the calving of ice in lakes that temporarily surround the ice sheet. In order to capture this process a "flood-fill algorithm" is needed. Here we present and evaluate several optimizations to a standard flood-fill algorithm in terms of computational efficiency. As an example, we determine the land-ocean mask for a 1 km resolution digital elevation model (DEM) of North America and Greenland, a geographical area of roughly 7000 by 5000 km (roughly 35 million elements), about half of which is covered by ocean. Determining the land-ocean mask with our improved flood-fill algorithm reduces computation time by 90 % relative to using a standard stack-based flood-fill algorithm. This implies that it is now feasible to include the calving of ice in lakes as a dynamical process inside an ice-sheet model. We demonstrate this by using bedrock elevation, ice thickness and geoid perturbation fields from the output of a coupled ice-sheet-sea-level equation model at 30 000 years before present and determine the extent of Lake Agassiz, using both the standard and improved versions of the flood-fill algorithm. We show that several optimizations to the flood-fill algorithm used for filling a depression up to a water level, which is not defined beforehand, decrease the computation time by up to 99 %. The resulting reduction in computation time allows determination of the extent and volume of depressions in a DEM over large geographical grids or repeatedly over long periods of time, where computation time might otherwise be a limiting factor. The algorithm can be used for all glaciological and hydrological models, which need to trace the evolution over time of lakes or drainage basins in general

The Utrecht Finite Volume Ice-Sheet Model: UFEMISM (version 1.0)

Improving our confidence in future projections of sea-level rise requires models that can simulate ice-sheet evolution both in the future and in the geological past. A physically accurate treatment of large changes in ice-sheet geometry requires a proper treatment of processes near the margin, like grounding line dynamics, which in turn requires a high spatial resolution in that specific region, so that small-scale topographical features are resolved. This leads to a demand for computationally efficient models, where such a high resolution can be feasibly applied in simulations of 105–107 years in duration. Here, we present and evaluate a new ice-sheet model that solves the hybrid SIA–SSA approximation of the stress balance, including a heuristic rule for the grounding-line flux. This is done on a dynamic adaptive mesh which is adapted to the modelled ice-sheet geometry during a simulation. Mesh resolution can be configured to be fine only at specified areas, such as the calving front or the grounding line, as well as specified point locations such as ice-core drill sites. This strongly reduces the number of grid points where the equations need to be solved, increasing the computational efficiency. A high resolution allows the model to resolve small geometrical features, such as outlet glaciers and sub-shelf pinning points, which can significantly affect large-scale ice-sheet dynamics. We show that the model reproduces the analytical solutions or model intercomparison benchmarks for a number of schematic ice-sheet configurations, indicating that the numerical approach is valid. Because of the unstructured triangular mesh, the number of vertices increases less rapidly with resolution than in a square-grid model, greatly reducing the required computation time for high resolutions. A simulation of all four continental ice sheets during an entire 120 kyr glacial cycle, with a 4 km resolution near the grounding line, is expected to take 100–200 wall clock hours on a 16-core system (1600–3200 core hours), implying that this model can be feasibly used for high-resolution palaeo-ice-sheet simulations

Reconstructing the evolution of ice sheets, sea level, and atmospheric CO2during the past 3.6 million years

Understanding the evolution of, and the interactions between, ice sheets and the global climate over geological timescales is important for being able to project their future evolution. However, direct observational evidence of past CO2 concentrations, and the implied radiative forcing, only exists for the past 800 000 years. Records of benthic δ18O date back millions of years but contain signals from both land ice volume and ocean temperature. In recent years, inverse forward modelling has been developed as a method to disentangle these two signals, resulting in mutually consistent reconstructions of ice volume, temperature, and CO2. We use this approach to force a hybrid ice-sheet climate model with a benthic δ18O stack, reconstructing the evolution of the ice sheets, global mean sea level, and atmospheric CO2 during the late Pliocene and the Pleistocene, from 3.6 million years (Myr) ago to the present day. During the warmerthan-present climates of the late Pliocene, reconstructed CO2 varies widely, from 320 440 ppmv for warm periods to 235 250 ppmv for the early glacial excursion ∼ 3:3 million years ago. Sea level is relatively stable during this period, with maxima of 6 14 m and minima of 12 26 m during glacial episodes. Both CO2 and sea level are within the wide ranges of values covered by available proxy data for this period. Our results for the Pleistocene agree well with the ice-core CO2 record, as well as with different available sea-level proxy data. For the Early Pleistocene, 2.6 1.2 Myr ago, we simulate 40 kyr glacial cycles, with interglacial CO2 decreasing from 280 300 ppmv at the beginning of the Pleistocene to 250 280 ppmv just before the Mid-Pleistocene Transition (MPT). Peak glacial CO2 decreases from 220 250 to 205 225 ppmv during this period. After the MPT, when the glacial cycles change from 40 to 80 120 kyr cyclicity, the glacial interglacial contrast increases, with interglacial CO2 varying between 250 320 ppmv and peak glacial values decreasing to 170 210 ppmv

Reconstructing the evolution of ice sheets, sea level, and atmospheric CO2during the past 3.6 million years

Understanding the evolution of, and the interactions between, ice sheets and the global climate over geological timescales is important for being able to project their future evolution. However, direct observational evidence of past CO2 concentrations, and the implied radiative forcing, only exists for the past 800 000 years. Records of benthic δ18O date back millions of years but contain signals from both land ice volume and ocean temperature. In recent years, inverse forward modelling has been developed as a method to disentangle these two signals, resulting in mutually consistent reconstructions of ice volume, temperature, and CO2. We use this approach to force a hybrid ice-sheet climate model with a benthic δ18O stack, reconstructing the evolution of the ice sheets, global mean sea level, and atmospheric CO2 during the late Pliocene and the Pleistocene, from 3.6 million years (Myr) ago to the present day. During the warmerthan-present climates of the late Pliocene, reconstructed CO2 varies widely, from 320 440 ppmv for warm periods to 235 250 ppmv for the early glacial excursion ∼ 3:3 million years ago. Sea level is relatively stable during this period, with maxima of 6 14 m and minima of 12 26 m during glacial episodes. Both CO2 and sea level are within the wide ranges of values covered by available proxy data for this period. Our results for the Pleistocene agree well with the ice-core CO2 record, as well as with different available sea-level proxy data. For the Early Pleistocene, 2.6 1.2 Myr ago, we simulate 40 kyr glacial cycles, with interglacial CO2 decreasing from 280 300 ppmv at the beginning of the Pleistocene to 250 280 ppmv just before the Mid-Pleistocene Transition (MPT). Peak glacial CO2 decreases from 220 250 to 205 225 ppmv during this period. After the MPT, when the glacial cycles change from 40 to 80 120 kyr cyclicity, the glacial interglacial contrast increases, with interglacial CO2 varying between 250 320 ppmv and peak glacial values decreasing to 170 210 ppmv