# Numerical Model Explores Subglacial Cavity Formation in Ice Sheets

Ice sheets—also called continental glaciers—are large masses of glacial ice that comprise more than 50,000 square kilometers. Though they blanketed much of the Earth many years ago, only two ice sheets remain: the Antarctic ice sheet (covering roughly 98 percent of Antarctica, see Figure 1) and the Greenland ice sheet (covering nearly 80 percent of Greenland). These sheets, which have been losing mass at an increasing rate since the 1990s, collectively hold a volume of ice that would raise the sea level by 65 meters if it were to melt. “Understanding ice sheet dynamics requires a combined knowledge of many processes,” Gonzalo Gonzalez de Diego of the University of Oxford said. These processes include the way in which an ice sheet melts, flows, and slides over bedrock, as well as glacial meltwater’s percolation into the sheets and subsequent effect on the ice’s dynamics.

**Figure 1.**Satellite composite image of the Antarctic ice sheet, which covers 98 percent of the continent. Public domain image courtesy of Dave Pape.

During a minisymposium presentation at the 2023 SIAM Conference on Computational Science and Engineering, which is currently taking place in Amsterdam, the Netherlands, de Diego presented a numerical model for subglacial cavity formation that couples the Stokes equations with a free boundary equation. “We model the formation of subglacial cavities by solving the subglacial boundary layer problem in two dimensions,” de Diego said. He began with a Stokes equation that has three contact boundary conditions, which account for the ice’s detachment from bedrock. The first condition asserts that the ice cannot penetrate the bedrock; the second affirms that the ice can press against the bedrock but not pull away from it; and the third contends that the ice either detaches or compresses from the bedrock, but not simultaneously. “What makes this problem very different from its elastic counterpart is the fact that we’re working with a viscous medium, so we have velocity of displacement,” de Diego said. “Subglacial cavitation is one example of a viscous contact problem.”

To solve the problem numerically, de Diego decoupled the Stokes equation from the free boundary equation. “The free boundary equation is the novelty of this problem,” he said. He then acknowledged the need to account for certain numerical considerations, such as the formulation of the discrete variational inequality for the Stokes equation and the decision to solve the variational inequality with the penalty method or Lagrange multiplier method. A compatibility issue with the discrete variational inequality also exists. “We found that for some problems, solving the Stokes variational inequality may not deliver sufficiently accurate computations of the velocity field,” de Diego said. “This leads to a breakdown of the numerical scheme.”

De Diego’s proposed numerical algorithm efficiently handled these and other considerations. Specifically, he employed the Taylor-Hood finite element pair to find the velocity and pressure and utilized a Lagrange multiplier to enforce the contact boundary conditions. As a numerical example, de Diego computed a friction law over a sinusoidal bedrock by finding steady states for different pairs of sliding velocities. After obtaining a steady cavity, he repeatedly calculated the integral for basal shear, computed many steady states, and mapped the points—each of which is a steady state—in space. As one advances down the resulting curve, the size of the cavity continues to grow; a small change in contact area ultimately leads to a huge change in basal shear.

Given the significance of subglacial cavitation in the context of ice sheet dynamics, de Diego and his colleagues are currently expanding their research to understand the future evolution of marine ice sheets, which slide off continents and float in the ocean. They are also exploring subglacial lake problems—which involve the relaxation of ice sheets—and want to understand the mechanisms that transpire under the sheet when relaxation occurs. Finally, the team is investigating possible extensions of their numerical model to three dimensions, which would allow for the possibility of new formulations of the Stokes variational inequality.

Lina Sorg is the managing editor of SIAM News. |