Ice losses from the Greenland and Antarctic ice sheets have accelerated since the 1990s, accounting for a significant increase in the global mean sea level. Here, we present a new 29-year record of ice sheet mass balance from 1992 to 2020 from the Ice Sheet Mass Balance Inter-comparison Exercise (IMBIE). We compare and combine 50 independent estimates of ice sheet mass balance derived from satellite observations of temporal changes in ice sheet flow, in ice sheet volume, and in Earth's gravity field. Between 1992 and 2020, the ice sheets contributed 21.0±1.9g€¯mm to global mean sea level, with the rate of mass loss rising from 105g€¯Gtg€¯yr-1 between 1992 and 1996 to 372g€¯Gtg€¯yr-1 between 2016 and 2020. In Greenland, the rate of mass loss is 169±9g€¯Gtg€¯yr-1 between 1992 and 2020, but there are large inter-annual variations in mass balance, with mass loss ranging from 86g€¯Gtg€¯yr-1 in 2017 to 444g€¯Gtg€¯yr-1 in 2019 due to large variability in surface mass balance. In Antarctica, ice losses continue to be dominated by mass loss from West Antarctica (82±9g€¯Gtg€¯yr-1) and, to a lesser extent, from the Antarctic Peninsula (13±5g€¯Gtg€¯yr-1). East Antarctica remains close to a state of balance, with a small gain of 3±15g€¯Gtg€¯yr-1, but is the most uncertain component of Antarctica's mass balance. The dataset is publicly available at 10.5285/77B64C55-7166-4A06-9DEF-2E400398E452 (IMBIE Team, 2021).

In this study, we examine the effect of transient mantle creep on the prediction of glacial isostatic adjustment (GIA) signals. Specifically, we compare predictions of relative sea level (RSL) change from GIA from a set of Earth models in which transient creep parameters are varied in a simple Burgers model to a reference case with a Maxwell viscoelastic rheology. The model predictions are evaluated in two ways: first, relative to each other to quantify the effect of parameter variation, and second, for their ability to reproduce well-constrained sea level records from selected locations. Both the resolution and geographic location of the RSL observations determine whether the data can distinguish between model cases. Model predictions are most sensitive to the inclusion of transient mantle deformation in regions that are near-field and peripheral relative to former ice sheets. This sensitivity appears particularly true along the North American west coast in the region of the former Cordilleran Ice Sheet, which experienced rapid sea-level fall following deglaciation between 14 and 12 kyr BP. Relative to the Maxwell case, Burgers models better reproduce this rapid phase of regional postglacial sea-level fall. As well, computed goodness-of-fit values in this region show a clear preference for models where transient deformation is present in the whole or lower mantle, and for models where the rigidity of the Kelvin element is weakened relative to the rigidity of the Maxwell element. In contrast, model predictions of relative sea-level change in the far-field show weak sensitivity to the inclusion of transient deformation.

In this study, we focus on improved constraint of the glacial isostatic adjustment (GIA) signal at present-day, and its role as a contributor to present-day sea level budgets. The main study area extends from the coastal regions of northwestern Europe to northern Europe. Both Holocene relative sea level (RSL) data as well as vertical land motion (VLM) data are incorporated as constraints in a semi-empirical GIA model. 71 geological rates of GIA-driven RSL change are inferred from Holocene proxy data and 108 rates of vertical land motion from GNSS provide an additional measure of regional GIA deformation. Within the study area, the geological RSL data complement the spatial gaps of the VLM data and vice versa. Both data sets are inverted in a semi-empirical GIA model to yield updated estimates of regional present-day GIA deformations. A regional validation using tide gauges is presented for the North Sea, where the GIA signal may be complicated by lateral variations in Earth structure and existing predictions of regional and global GIA models show discrepancies. The model validation in the North Sea region suggests that geological data are needed to fit independent estimates of GIA-related RSL change inferred from tide gauge rates, indicating that geological rates from Holocene data do provide an important additional constraint for data-driven approaches to GIA estimation.

This work provides a comparison of four approaches that can be used to describe uncertainty in models of the long-term glacial isostatic adjustment (GIA) process. The four methods range from pessimistic to optimistic representations of GIA uncertainty. Each estimation method is applied to selected one dimensional GIA model predictions and compared with vertical land motion data from Global Positioning System (GPS) measurements across Fennoscandia and North America. The methods are evaluated relative to two main properties: (1) their expected ability to separate non-GIA from GIA signals and (2) their estimated statistical appropriateness given a specific GIA model and data set. For the first point, non-GIA signals are considered isolated from the long-term (millennial time scale) GIA signal at sites where measurement and model uncertainties do not overlap. Across methods, the frequency and accuracy with which non-GIA signals are separated from GIA signals in GPS data display both consistent similarities and disparities. For the second point, we compare model predictions with rates of vertical land motion and relative sea level change that have been cleaned of non-GIA signals to determine the most appropriate value of model uncertainty and relate the findings to the four approaches. Best fit inferences suggest that within deglaciation centers, GIA model uncertainty is up to ~2 mm/yr (vertical land motion). Likewise, away from the former ice sheet centers, GIA uncertainty for relative sea level change is inferred to be ~0.3–0.5 mm/yr along the U.S. East Coast and ~0.6–0.8 mm/yr in the North Sea.

The glacial isostatic adjustment (GIA) signal at present day is constrained via the joint inversion of geodetic observations and GIA models for a region encompassing northern Europe, the British Isles, and the Barents Sea. The constraining data are Global Positioning System (GPS) vertical crustal velocities and GRACE (Gravity Recovery and Climate Experiment) gravity data. When the data are inverted with a set of GIA models, the best-fit model for the vertical motion signal has a χ 2 value of approximately 1 and a maximum a posteriori uncertainty of 0.3-0.4mm yr-1. An elastic correction is applied to the vertical land motion rates that accounts for present-day changes to terrestrial hydrology as well as recent mass changes of ice sheets and glaciered regions. Throughout the study area, mass losses from Greenland dominate the elastic vertical signal and combine to give an elastic correction of up to +0.5mm yr-1 in central Scandinavia. Neglecting to use an elastic correction may thus introduce a small but persistent bias in model predictions of GIA vertical motion even in central Scandinavia where vertical motion is dominated by GIA due to past glaciations. The predicted gravity signal is generally less well-constrained than the vertical signal, in part due to uncertainties associated with the correction for contemporary ice mass loss in Svalbard and the Russian Arctic. The GRACE-derived gravity trend is corrected for present-day ice mass loss using estimates derived from the ICESat and CryoSat missions, although a difference in magnitude between GRACE-inferred and altimetry-inferred regional mass loss rates suggests the possibility of a non-negligible GIA response here either from millennial-scale or Little Ice Age GIA.

The ocean load in glacial isostatic adjustment (GIA) modelling is represented by the so-called sea level equation (SLE). The SLE describes the mass redistribution of water between ice sheets and oceans on a deforming Earth. Despite various teams independently investigating GIA, there has been no systematic intercomparison among the numerical solvers of the SLE through which the methods may be validated. The goal of this paper is to present a series of synthetic examples designed for testing and comparing the numerical implementations of the SLE in GIA modelling. The 10 numerical codes tested combine various temporal and spatial parametrizations. The time-domain or Laplace-domain discretizations are used to solve the SLE through time, while spherical harmonics, finite differences or finite elements parametrize the GIA-related field variables spatially. The surface ice-water load and solid Earth's topography are represented spatially either on an equiangular grid, a Gauss-Legendre or an equiarea grid with icosahedron-shaped spherical pixels. Comparisons aremade in a series of five benchmark examples with an increasing degree of complexity. Due to the complexity of the SLE, there is no analytical solution to it. The accuracy of the numerical implementations is therefore assessed by the differences of the individual solutions with respect to a reference solution. While the benchmark study does not result in GIA predictions for a realistic loading scenario, we establish a set of agreed-upon results that can be extended in the future by including more complex case studies, such as solutions with realistic loading scenarios, the rotational feedback in the linear-momentum equation, and by considering a 3-D viscosity structure of the Earth's mantle. The test computations performed so far show very good agreement between the individual results and their ability to capture the main features of sea-surface variation and the surface vertical displacement. The differences found can often be attributed to the different approximations inherent in the various algorithms. This shows the accuracy that can be expected from different implementations of the SLE, which helps to assess differences noted in the literature between predictions for realistic loading cases.

Sea-level rise and decadal variability along the northwestern coast of the North Atlantic Ocean are studied in a self-consistent framework that takes into account the effects of solid-earth deformation and geoid changes due to large-scale mass redistribution processes. Observations of sea and land level changes from tide gauges and GPS are compared to the cumulative effect of GIA, present-day mass redistribution, and ocean dynamics over a 50 year period (1965–2014). GIA explains the majority of the observed sea-level and land motion trends, as well as almost all interstation variability. Present-day mass redistribution resulting from ice melt and land hydrology causes both land uplift and sea-level rise in the region. We find a strong correlation between decadal steric variability in the Subpolar Gyre and coastal sea level, which is likely caused by variability in the Labrador Sea that is propagated southward. The steric signal explains the majority of the observed decadal sea-level variability and shows an upward trend and a significant acceleration, which are also found along the coast. The sum of all contributors explains the observed trends in both sea-level rise and vertical land motion in the region, as well as the decadal variability. The sum of contributors also explains the observed acceleration within confidence intervals. The sea-level acceleration coincides with an accelerating density decrease at high latitudes.

Geodetic measurements of vertical land motion and gravity change are incorporated into an a priori model of present-day glacial isostatic adjustment (GIA) in North America via least-squares adjustment. The result is an updated GIA model wherein the final predicted signal is informed by both observational data, and prior knowledge (or intuition) of GIA inferred from models. The data-driven method allows calculation of the uncertainties of predicted GIA fields, and thus offers a significant advantage over predictions from purely forward GIA models. In order to assess the influence each dataset has on the final GIA prediction, the vertical land motion and GRACE-measured gravity data are incorporated into the model first independently (i.e., one dataset only), then simultaneously. The relative weighting of the datasets and the prior input is iteratively determined by variance component estimation in order to achieve the most statistically appropriate fit to the data. The best-fit model is obtained when both datasets are inverted and gives respective RMS misfits to the GPS and GRACE data of 1.3 mm/yr and 0.8 mm/yr equivalent water layer change. Non-GIA signals (e.g., hydrology) are removed from the datasets prior to inversion. The post-fit residuals between the model predictions and the vertical motion and gravity datasets, however, suggest particular regions where significant non-GIA signals may still be present in the data, including unmodeled hydrological changes in the central Prairies west of Lake Winnipeg. Outside of these regions of misfit, the posterior uncertainty of the predicted model provides a measure of the formal uncertainty associated with the GIA process; results indicate that this quantity is sensitive to the uncertainty and spatial distribution of the input data as well as that of the prior model information. In the study area, the predicted uncertainty of the present-day GIA signal ranges from ∼0.2-1.2 mm/yr for rates of vertical land motion, and from ∼3-4 mm/yr of equivalent water layer change for gravity variations.

The thickness and equivalent global sea level contribution of an improved model of the central and northern Laurentide Ice Sheet is constrained by 24 relative sea level histories and 18 present-day GPS-measured vertical land motion rates. The final model, termed Laur16, is derived from the ICE-5G model by holding the timing history constant and iteratively adjusting the thickness history, in four regions of northern Canada. In the final model, the last glacial maximum (LGM) thickness of the Laurentide Ice Sheet west of Hudson Bay was ~3.4-3.6 km. Conversely, east of Hudson Bay, peak ice thicknesses reached ~4 km. The ice model thicknesses inferred for these two regions represent, respectively, a ~30 per cent decrease and an average ~20-25 per cent increase to the load thickness relative to the ICE-5G reconstruction, which is generally consistent with other recent studies that have focussed on Laurentide Ice Sheet history. The final model also features peak ice thicknesses of 1.2-1.3 km in the Baffin Island region, a modest reduction relative to ICE-5G and unchanged thicknesses for a region in the central Canadian Arctic Archipelago west of Baffin Island. Vertical land motion predictions of the final model fit observed crustal uplift rates well, after an adjustment is made for the elastic crustal response to present-day ice mass changes of regional ice cover. The new Laur16 model provides more than a factor of two improvement of the fit to the RSL data (?2 measure of misfit) and a factor of nine improvement to the fit of the GPS data (mean squared error measure of fit), compared to the ICE-5G starting model. Laur16 also fits the regional RSL data better by a factor of two and gives a slightly better fit to GPS uplift rates than the recent ICE-6G model. The volume history of the Laur16 reconstruction corresponds to an up to 8 m reduction in global sea level equivalent compared to ICE-5G at LGM.