Adaptive Mesh Refinement Predicts Future Melting of the Antarctic Ice Sheet

By Lina Sorg

On January 21, 2002, researchers observing satellite images of the Antarctic Peninsula watched the beginning of an alarmingly unexpected phenomenon. The Larsen B Ice Shelf, located in the northwest part of the Weddell Sea, began to splinter and collapse. The shelf continued to disintegrate over the next month, ultimately losing 1,250 square miles of ice. Scientists suspect that the shelf’s collapse was spurred by a sequence of warm summers in Antarctica, culminating with the remarkably warm summer of 2002. The warmer air triggered melting on the ice’s surface, creating melt pools that deepened existing crevices. “Melt pools on the surface turned into knives and chopped the shelf up into little bits,” said Daniel Martin of Lawrence Berkeley National Laboratory. Warm ocean temperatures in the Weddell Sea during this time simultaneously eroded ice underneath the shelf. “Why do we care?” Martin asked. “The answer is sea level rise.”

Since 1993, the average global sea level has risen approximately 3.2 millimeters per year. Melting ice sheets—defined as masses of glacial land ice that are over 50,000 square kilometers in area—are responsible for roughly 33 percent of sea level rise. Two ice sheets exist on Earth today: the Greenland ice sheet and the Antarctic ice sheet (AIS). As the largest single mass of ice on the planet, the AIS covers 98 percent of Antarctica and contains 26.5 million cubic kilometers of ice; this translates to 61 percent of Earth’s freshwater supply and a sea level equivalent of 60 meters.

As ice sheets expand and move towards/over the ocean, they transition from grounded ice to floating ice shelves. The point at which this transition occurs is called the grounding line. In present-day conditions, the AIS remains one of the largest uncertainties in the projection of sea level rise. The section in West Antarctica, which is not situated on a land mass and extends over 2,500 meters below sea level, is particularly vulnerable to collapse due to the incursion of warm water under existing ice shelves at the grounding line. Such intrusion causes a loss of buttressing—leaving the underwater ice shelves unsupported—and a corresponding retreat of the grounding line. Satellite images suggest that the shrinking AIS is one of the biggest contributors to rising sea levels. If the AIS were to collapse, as it did tens of thousands of years ago under similar climate conditions, the resulting sea level rise could generate a global flood.

Until recently, researchers were unable to accurately model the physical processes causing the AIS to retreat. However, during a minisymposium presentation at the 2019 SIAM Conference on Computational Science and Engineering, currently taking place in Spokane, Wash., Daniel Martin of Lawrence Berkeley National Laboratory (LBNL) presented a scalable adaptive mesh refinement (AMR) ice sheet model that examines full-continent AIS response to climate forcing. The model—called BISICLES—is the product of a collaboration between computational scientists at LBNL and climate scientists at Los Alamos National Laboratory and the University of Bristol. It is in its second round of funding from the Scientific Discovery through Advanced Computing.

BISICLES was built using LBNL’s Chombo software framework for block-structured AMR, which supports AMR discretizations, includes scalable solvers, and improves modeling accuracy. “Ice sheets are funny because there are really localized regions which require high resolution to accurately resolve ice-sheet dynamics,” Martin said. “But there are also really large regions where you don’t need fine resolutions.” This makes ice sheet modeling particularly well-suited for AMR, which allows researchers to study regions and processes of interest—such as the AIS’s retreating edge, glacier surges, ice streams, and grounding-line migration—at very high resolutions while reserving computationally cheaper coarse resolutions for areas that change more slowly and thus do not demand extreme detail. The high resolutions allowed by AMR identify a significant grounding line retreat by 2100, a phenomenon that satellites have already observed from space but previous simulations (with coarser resolutions) failed to identify.

Martin then briefed the audience on the physics of ice. “Ice is a non-Newtonian viscous fluid,” he said. “As you heat it, it becomes less viscous.” Although these properties are seemingly well-described by a full Stokes system, Martin chooses an approximate Stokes system that uses scaling arguments to produce a simpler set of equations and approximate systems that are easier to manage. In short, he utilizes the asymptotic structure of the full Stokes system to construct a higher-order approximation. This technique is computationally cheaper but still allows the recovery of a three-dimensional flow model.

Martin’s baseline model consists of a two-dimensional (2D) equation for ice thickness and a logically-rectangular grid obtained via time-dependent uniform mapping. The vertically-integrated momentum balance yields a 2D nonlinear viscous tensor solve for velocity at the base of the ice. Next, Martin employs finite volume discretizations to simplify the implementation of local refinement. Since 90 percent of computational time is spent in an iterative nonlinear viscous tensor solve for the ice velocity field, good linear solver performance is essential. This is this point at which the Chombo native geometric multigrid solvers come into play. “Chombo is a very useful framework on which to build these applications,” Martin said. It follows naturally from AMR hierarchy, is matrix-free, and is relatively efficient. When it works, it works really well. A select few problems do cause Chombo solvers to stall and thus instead require algebraic multigrid solvers.

After setting up his model, Martin applies it to an experiment with 1,000-year Antarctic simulations. At the initial time, he subjects ice shelves to extreme (outlandish) depth-dependent melting in which the floating ice gets clobbered. The grounding lines in the resulting simulation rapidly evolve and carve out kilometers of area in every direction. To specifically evaluate Antarctic vulnerability, Martin restricts the extreme melting to regional forcing and divides the AIS into susceptible areas with the potential to contribute to West Antarctica’s ultimate breakdown. “It’s so vulnerable that you can attack it from any direction and it’ll eventually collapse,” he said. Additionally, because resource limitations often force scientists to evaluate areas based on individual regions (rather than in full), Martin broke the western region of the AIS into sections and evaluated one region’s effect on another. However, this type of regional analysis is not without limitations, as Antarctic drainage bins are dynamically independent for only one or two centuries. “If you are creating a regional model, you can have confidence that it will be reasonable for 100-200 years," he said. "Then you have to throw it away."

Ultimately, Martin presented a fully-resolved, AMR-based systematic study of millennial-scale ice sheet response to regional sheet collapse. His BISICLES model yields an accurate, high-resolution snapshot of these processes at a manageable computational cost, which is crucial in the study and accurate prediction of future ice mass loss.

Lina Sorg is the associate editor of SIAM News.