Viscoelastic rebound of the solid Earth upon the removal of ice loads has the potential to inhibit marine ice sheet instability, thereby forestalling ice-sheet retreat and global mean sea-level rise. The timescale over which the solid Earth - ice sheet system responds to changes in ice thickness and bedrock topography places a strong control on the spatiotemporal influence of this negative feedback mechanism. In this study, we assess the impact of solid-earth rheological structure on model projections of the retreat of Thwaites Glacier, West Antarctica, and the concomitant sea-level rise by coupling the dynamic ice sheet model MALI to a regional glacial isostatic adjustment (GIA) model. We test the sensitivity of model projections of ice-sheet retreat and associated sea-level rise across a range of four different solid-earth rheologies, forced by standard ISMIP6 ocean and atmospheric datasets for the RCP8.5 climate scenario. These model parameters are applied to 500-year, coupled ice-sheet - GIA simulations. For the mantle viscosity best supported by observations, the negative GIA feedback leads to a reduction in mass loss that remains above 20% after about a hundred years. Mass-loss reduction peaks at 50% around 2300, which is when a control simulation without GIA experiences its maximum rate of retreat. For a weaker solid-earth rheology that is unlikely but compatible with observational uncertainty, mass loss reduction remains above 50% after 2150. At 2100, mass loss reduction is 10% for the best-fit rheology and 25% for the weakest rheology. At the same time, we estimate that water expulsion from the rebounding solid Earth beneath the ocean near Thwaites Glacier may increase sea-level rise by up to 20% at five centuries. Additionally, the reduction in ice-sheet retreat caused by GIA is substantially reduced under stronger climate forcings, suggesting that the stabilizing feedback of GIA will also be an indirect function of emissions scenario. We hypothesize that feedbacks between the solid Earth - ice sheet system are controlled by a competition between the spatial extent and timescale of bedrock uplift relative to the rate of grounded ice retreat away from the region of most rapid unloading. Although uncertainty in solid-earth rheology leads to large uncertainty in future sea-level rise contribution from Thwaites Glacier, under all plausible parameters the GIA effects are too large to be ignored for future projections of Thwaites Glacier of more than a century.

Inaccurate representations of iceberg calving from ice shelves are a large source of uncertainty in mass-loss projections from the Antarctic ice sheet. Here, we address this limitation by implementing and testing a continuum damage-mechanics model in a continental scale ice-sheet model. The damage-mechanics formulation, based on a linear stability analysis and subsequent long-wavelength approximation of crevasses that evolve in a viscous medium, links damage evolution to climate forcing and the large-scale stresses within an ice shelf. We incorporate this model into the BISICLES ice-sheet model and test it by applying it to idealized (1) ice tongues, for which we present analytical solutions and (2) buttressed ice-shelf geometries. Our simulations show that the model reproduces the large disparity in lengths of ice shelves with geometries and melt rates broadly similar to those of four Antarctic ice shelves: Erebus Glacier Tongue (length ~ 13 km), the unembayed portion of Drygalski Ice Tongue (~ 65 km), the Amery Ice Shelf (~ 350 km) and the Ross Ice Shelf (~ 500 km). These results demonstrate that our simple continuum model holds promise for constraining realistic ice-shelf extents in large-scale ice-sheet models in a computationally tractable manner.

The ice sheets of the Amundsen Sea Embayment (ASE) are vulnerable to the marine ice sheet instability (MISI), which could cause irreversible collapse and raise sea levels by over a meter. The uncertain timing and scale of this collapse depend on the complex interaction between ice, ocean, and bedrock dynamics. The mantle beneath the ASE is likely less viscous (~1018 Pa s) than the Earth’s average mantle (~1021 Pa s). Here we show that an effective equilibrium between Pine Island Glacier’s retreat and the response of a weak viscoelastic mantle can reduce ice mass lost by almost 30 per cent over 150 years. Other components of solid Earth response—purely elastic deformations and geoid perturbations—provide less stability than the viscoelastic response alone. Uncertainties in mantle rheology, topography, and basal melt affect how much stability we expect, if any. Our study indicates the importance of considering viscoelastic uplift during the rapid retreat associated with MISI.

The viscoelastic load Love numbers encapsulate the Earth’s rheology in a remarkably efficient fashion. When multiplied by a sudden increment of spherical harmonic load change, they give the horizontal and vertical surface displacements and gravity change at all later times. Incremental glacial load changes thus need only be harmonically decomposed, multiplied by the Love numbers and summed to predict the Earth’s response to glacial load redistributions. The computation of viscoelastic Love numbers from the elastic, viscous and adiabatic profiles of the Earth is thus the foundation upon which many glacial isostatic adjustment models are based. Usually, viscoelastic Love numbers are computed using the Laplace transform method, employing the correspondence principle to convert the viscoelastic equations of motion into the elastic equations with complex material parameters. This method works well for a fully non-adiabatic Earth, but can accommodate realistic partially adiabatic and fully adiabatic conditions only by changing the Earth’s density profile. An alternative method of Love number computation developed by Cathles (1975) avoids this dilemma by separating the elastic and viscous equations of motion. The separation neglects a small solid-elastic/fluid-elastic transition for compressible deformation, but allows freely defining adiabatic, partially adiabatic or fully non-adiabatic profiles in the mantle without changing the Earth’s density profile. Here, we update and fully describe this method and show that it produces Love numbers closely similar to those computed for fully non-adiabatic earth models computed by the correspondence principle, finite element and other methods. The time-domain method produces Love numbers as good as those produced by other methods and can also realistically accommodate any degree of mantle adiabaticity. All method implementations are available open source.

Elastic deformation of the solid Earth in response to ice mass loss offers a promising constraint on the density of glacial material lost. Further, the elastic response to modern deglaciation is important to constrain for studies of glacial isostatic adjustment to determine the mantle’s structure and rheology. Models of this elastic uplift are commonly based on the 1-D, seismically derived global average Preliminary Reference Earth Model and typically neglect uncertainties that can arise from regional differences in elastic structure from that of the global average, lateral heterogeneities within the region, and inelastic behavior of the crust. We quantify these uncertainties using an ensemble of 1-D local elastic structure models and empirical relations for the effects of inelasticity in the upper 10 km of the crust. In Southeast Alaska, modeling elastic uplift rates with local elastic structures results in up to a 20–40 percent difference from those modeled with the Preliminary Reference Earth Model. Although these differences are limited to regions near to ice-covered areas, they are comparable to the differences in uplift rates expected from the loss of firn versus loss of ice. Far from ice-covered areas, where most of the region’s GPS observations were made, these differences become insignificant and do not affect previous glacial isostatic adjustment studies in the region. The methods presented here are based on the globally available LITHO1.0 seismic model and open source software, and the approach of using an ensemble of 1-D elastic structures can be easily adapted to other regions around the world.

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 are made 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.

The ice load configuration of the Barents Sea Ice Sheet (BSIS) over the last glacial cycle is in dispute. The traditional reconstruction, motivated by the observation that paleo-shoreline emergence increases toward the center of the Barents Sea, places a single dome in the center of the Barents Sea at the last glacial maximum (LGM) that collapses to island centered loads during deglaciation. Observations that suggest that ice flowed from the islands into the Barents even at the LGM motivate another reconstruction that places the ice loads over the islands with minimal marine ice. We analyze an ensemble of ice loads that are consistent with the geophysical observations using relatively new statistical methods and show that current relative sea level, GPS and gravity measurement do not and cannot distinguish a central dome from an island-centered BSIS. What is needed are constraints in the central Barents Sea. Improving the gravity data sufficiently will be difficult. However, obtaining even a single GPS uplift rate measurement in the central Barents would resolve the central dome versus island centered BSIS geometry question. The Barents Sea ice load geometry uncertainty provides a good illustration of statistical methods that may be useful in other areas of glaciology.

A freely cooling granular gas in a gravitational field undergoes a collapse to a multicontact state in a finite time. Previous theoretical [D. Volfson et al., Phys. Rev. E 73, 061305 (2006)] and experimental work [R. Son et al., Phys. Rev. E 78, 041302 (2008)] have obtained contradictory results about the rate of energy loss before the gravitational collapse. Here we use a molecular dynamics simulation in an attempt to recreate the experimental and theoretical results to resolve the discrepancy. We are able to nearly match the experimental results, and find that to reproduce the power law predicted in the theory we need a nearly elastic system with a constant coefficient of restitution greater than 0.993. For the more realistic velocity-dependent coefficient of restitution, there does not appear to be a power-law decay and the transition from granular gas to granular solid is smooth, making it difficult to define a time of collapse.

We study shock wave formation and propagation in an experimental vertically driven quasi-two-dimensional granular gas. We measure the moments of the single particle velocity distribution as a function of space and time. The space-time fields of the velocity moments show acoustic waves with a serrated substructure on the scale of a particle diameter. We show that this substructure is the result of collisional transport in which sequential collisions each transport momentum and energy by one particle diameter.

The ice sheets of the Amundsen Sea Embayment (ASE) are vulnerable to the marine ice sheet instability (MISI), which could cause irreversible collapse and raise sea levels by over a meter. The uncertain timing and scale of this collapse depend on the complex interaction between ice, ocean, and bedrock dynamics. The mantle beneath the ASE is likely less viscous (~1018 Pa s) than the Earth’s average mantle (~1021 Pa s). Here we show that an effective equilibrium between Pine Island Glacier’s retreat and the response of a weak viscoelastic mantle can reduce ice mass lost by almost 30 per cent over 150 years. Other components of solid Earth response—purely elastic deformations and geoid perturbations—provide less stability than the viscoelastic response alone. Uncertainties in mantle rheology, topography, and basal melt affect how much stability we expect, if any. Our study indicates the importance of considering viscoelastic uplift during the rapid retreat associated with MISI.

The viscoelastic load Love numbers encapsulate the Earth’s rheology in a remarkably efficient fashion. When multiplied by a sudden increment of spherical harmonic load change, they give the horizontal and vertical surface displacements and gravity change at all later times. Incremental glacial load changes thus need only be harmonically decomposed, multiplied by the Love numbers and summed to predict the Earth’s response to glacial load redistributions. The computation of viscoelastic Love numbers from the elastic, viscous and adiabatic profiles of the Earth is thus the foundation upon which many glacial isostatic adjustment models are based. Usually, viscoelastic Love numbers are computed using the Laplace transform method, employing the correspondence principle to convert the viscoelastic equations of motion into the elastic equations with complex material parameters. This method works well for a fully non-adiabatic Earth, but can accommodate realistic partially adiabatic and fully adiabatic conditions only by changing the Earth’s density profile. An alternative method of Love number computation developed by Cathles (1975) avoids this dilemma by separating the elastic and viscous equations of motion. The separation neglects a small solid-elastic/fluid-elastic transition for compressible deformation, but allows freely defining adiabatic, partially adiabatic or fully non-adiabatic profiles in the mantle without changing the Earth’s density profile. Here, we update and fully describe this method and show that it produces Love numbers closely similar to those computed for fully non-adiabatic earth models computed by the correspondence principle, finite element and other methods. The time-domain method produces Love numbers as good as those produced by other methods and can also realistically accommodate any degree of mantle adiabaticity. All method implementations are available open source.

Elastic deformation of the solid Earth in response to ice mass loss offers a promising constraint on the density of glacial material lost. Further, the elastic response to modern deglaciation is important to constrain for studies of glacial isostatic adjustment to determine the mantle’s structure and rheology. Models of this elastic uplift are commonly based on the 1-D, seismically derived global average Preliminary Reference Earth Model and typically neglect uncertainties that can arise from regional differences in elastic structure from that of the global average, lateral heterogeneities within the region, and inelastic behavior of the crust. We quantify these uncertainties using an ensemble of 1-D local elastic structure models and empirical relations for the effects of inelasticity in the upper 10 km of the crust. In Southeast Alaska, modeling elastic uplift rates with local elastic structures results in up to a 20–40 percent difference from those modeled with the Preliminary Reference Earth Model. Although these differences are limited to regions near to ice-covered areas, they are comparable to the differences in uplift rates expected from the loss of firn versus loss of ice. Far from ice-covered areas, where most of the region’s GPS observations were made, these differences become insignificant and do not affect previous glacial isostatic adjustment studies in the region. The methods presented here are based on the globally available LITHO1.0 seismic model and open source software, and the approach of using an ensemble of 1-D elastic structures can be easily adapted to other regions around the world.

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 are made 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.

The ice load configuration of the Barents Sea Ice Sheet (BSIS) over the last glacial cycle is in dispute. The traditional reconstruction, motivated by the observation that paleo-shoreline emergence increases toward the center of the Barents Sea, places a single dome in the center of the Barents Sea at the last glacial maximum (LGM) that collapses to island centered loads during deglaciation. Observations that suggest that ice flowed from the islands into the Barents even at the LGM motivate another reconstruction that places the ice loads over the islands with minimal marine ice. We analyze an ensemble of ice loads that are consistent with the geophysical observations using relatively new statistical methods and show that current relative sea level, GPS and gravity measurement do not and cannot distinguish a central dome from an island-centered BSIS. What is needed are constraints in the central Barents Sea. Improving the gravity data sufficiently will be difficult. However, obtaining even a single GPS uplift rate measurement in the central Barents would resolve the central dome versus island centered BSIS geometry question. The Barents Sea ice load geometry uncertainty provides a good illustration of statistical methods that may be useful in other areas of glaciology.

A freely cooling granular gas in a gravitational field undergoes a collapse to a multicontact state in a finite time. Previous theoretical [D. Volfson et al., Phys. Rev. E 73, 061305 (2006)] and experimental work [R. Son et al., Phys. Rev. E 78, 041302 (2008)] have obtained contradictory results about the rate of energy loss before the gravitational collapse. Here we use a molecular dynamics simulation in an attempt to recreate the experimental and theoretical results to resolve the discrepancy. We are able to nearly match the experimental results, and find that to reproduce the power law predicted in the theory we need a nearly elastic system with a constant coefficient of restitution greater than 0.993. For the more realistic velocity-dependent coefficient of restitution, there does not appear to be a power-law decay and the transition from granular gas to granular solid is smooth, making it difficult to define a time of collapse.

We study shock wave formation and propagation in an experimental vertically driven quasi-two-dimensional granular gas. We measure the moments of the single particle velocity distribution as a function of space and time. The space-time fields of the velocity moments show acoustic waves with a serrated substructure on the scale of a particle diameter. We show that this substructure is the result of collisional transport in which sequential collisions each transport momentum and energy by one particle diameter.