Improving Ice Sheet Models Through Algebraic Multigrid Methods

By Karthika Swamy Cohen

Improved modeling of ice-sheet dynamics is essential to better estimate and reduce uncertainties in the projections of sea-level rise resulting from ice sheet mass loss, which is challenging due to the complexity of modeling ice sheet behavior.  

“An overarching goal in modeling the Greenland and Antarctic ice sheets is to help make improved sea-level rise projections and also to understand the impact that ice sheet mass loss will have on other parts of the climate system e.g., changes in ocean circulation,” says Ray Tuminaro, whose group at Sandia National Laboratories is working on improved ice sheet modeling, along with researchers at Los Alamos National Lab. “This requires reliable and efficient ice sheet models coupled with models for other earth system components such as the atmosphere, ocean, and sea ice.” 

In 2007, the Intergovernmental Panel on Climate Change (IPCC) called for better ice sheet models, due to the inability of ice sheet models at the time to predict or explain dynamic changes observed on ice sheets and ice caps (e.g., acceleration and thinning).  Several models have been developed in recent years.  Although current generation models are able to better represent ice sheet behavior than models from a decade ago, they are computationally very expensive. 

“The simulation of ice sheet dynamics frequently includes numerical techniques—nonlinear solvers, steady-state calculations, and implicit time advancement methods—that rely heavily on iterative linear solvers,” Tuminaro explains. “These linear solvers are generally the most time-consuming component within the entire simulation, and the convergence of the solver is often poor due to long-range effects associated with ice shelves associated with the non-local nature of ice flow and the highly-stretched meshes used to discretize ice sheets.

Sample mesh hierarchy with finest level (top row, left), two semicoarsening created levels (top row, right), and two SA-AMG created levels (bottom row). Figure credit: Ray Tuminaro, Mauro Perego, Irina Tezaur, Andrew Salinger and Stephen Price.

Specifically, ice sheet models require a (global) elliptic solve for the solution of the momentum balance equations and this generally makes up 90% or more of the total cost of the ice sheet model.”    Even though this may change as models mature and include more complex physical processes, it is likely that the elliptic solve will continue to be the costliest part of the model. 

In a paper published in the SIAM Journal on Scientific Computing, Tuminaro and co-authors Mauro Perego, Irina Tezaur, Andrew Salinger, and Stephen Price present a new multigrid solver that addresses this time-intensive component of large-scale, high-resolution ice sheet simulations.  The authors present a new multigrid solver that combines matrix-dependent multigrid ideas, semi-coarsening, and algebraic multigrid. Matrix-dependent multigrid is a family of algorithms that takes advantage of structure and can be viewed as ideal grid transfers for one-dimensional simplifications of higher-dimensional problems.

Multigrid linear solvers often converge rapidly by employing inexpensive, low-resolution models to accelerate convergence. However, these low-resolution models must then accurately capture long-range interactions in the higher-resolution linear system. 

“The newly developed algebraic multigrid solver rectifies significant deficiencies in the low resolution approximations and in so doing gives rise to rapid convergence rates even within very large-scale simulations,” Tuminaro says. “Requiring only a small number of iterations, say 20, as opposed to thousands, implies that the new solver is much less time intensive.”

The challenge with using a standard algebraic multigrid method in tackling ice sheet linear systems is that it does not sufficiently delineate between horizontal and vertical directions. Specifically, the low-resolution models tend to poorly capture subtle yet important horizontal couplings as the solver naturally focuses on larger matrix coefficients, which in this case correspond to vertical coupling.  

“Ice sheets are much broader in horizontal directions than in the vertical direction - they are often represented as ‘thin’ flows. For example, Antarctica covers an area larger than the United States while being only a couple of kilometers thick,” explains Tuminaro. 

“Yet the horizontal flow coupling is critical for accurately modeling important parts of the system, like ice streams and ice shelves, and this coupling must be accurately represented within the multigrid solver,” adds co-author Stephen Price.

 “Related multigrid solver work has leveraged properties of the vertical and horizontal directions in ice sheets,” Tuminaro says. “However, our algorithm is the first algebraic approach that we are aware of that truly leverages the extruded nature of the meshes used in ice sheet modeling. In particular, ice sheet meshes are generated by taking a two dimensional unstructured grid (map plane view) and extruding it in the vertical direction.”

Figure caption: Coastlines: Greenland triangle mesh (left) and Antarctica quadrilateral mesh (right). Two hinges are highlighted with yellow ovals. Figure credit: Ray Tuminaro, Mauro Perego, Irina Tezaur, Andrew Salinger and Stephen Price.

Extruded meshes are fairly routine in the case of thin structures, such as ice sheets, and often give rise to anisotropic problems since the mesh spacing in the thin direction is much smaller than the spacing in the broad direction. 

“The new solver explicitly takes advantage of the special nature of the extruded direction: low resolution approximations are constructed by first coarsening only in the extruded direction and then only coarsening in the horizontal direction,” says Tuminaro.

The new solver builds on several important multigrid ideas for structured meshes and for addressing anisotropic problems, developed over many decades. The group’s goal was to generalize structured mesh ideas to the situation where there is only structure in one dimension.

The proposed algorithm is fully algebraic, which enables it to integrate in a variety of different application systems.

“Perhaps the biggest and most important improvement in ice sheet models in the last decade or so is that they now routinely include some ‘higher order’ approximation of the governing nonlinear Stokes-flow equations,” says Tuminaro. “The first generation of ice sheet models that were explicitly pointed out in the IPCC fourth assessment report as being deficient at helping to explain contemporary observations of ice sheet change did not account for horizontal stress gradients in the momentum balance, only vertical shear-stress gradients. Almost all models—including this one—do now use a more sophisticated approximation of the governing equations, which has led to a large improvement in their ability to be used as tools to help understand and interpret contemporary changes on ice sheets. By inference, they should also be better tools to use in predicting future changes.”

A number of related ongoing research activities are taking place under the Department of Energy’s Predicting Ice Sheet and Climate Evolution at Extreme Scales (PISCEES) project. This includes the coupling of the ice sheet model to other earth system components, improvements to various aspects of model physics (e.g., subglacial processes that control the rate of basal motion), leveraging high performance computing characteristics within key computational kernels, higher fidelity simulations, extended and more stable time integration, verification, validation, model calibration, the quantification of uncertainties, and a better understanding of sensitivities to changes in climate forcing.

Acknowledgments: The author thanks the paper's authors Ray Tuminaro, Stephen Price, Mauro Perego, Irina Tezaur, and Andrew Salinger for their review, quotes, and detailed feedback on the article. 

Karthika Swamy Cohen is the managing editor of SIAM News.