A comparison of the stability and performance of depth-integrated ice-dynamics solvers

Abstract

Abstract. In the last decade, the number of ice-sheet models has increased substantially, in line with the growth of the glaciological community. These models use solvers based on different approximations of ice dynamics. In particular, several depth-integrated dynamics solvers have emerged as fast solvers capable of resolving the relevant physics of ice sheets at the continental scale. However, the numerical stability of these schemes has not been studied systematically to evaluate their effectiveness in practice. Here we focus on three such solvers, the so-called Hybrid, L1L2-SIA and DIVA solvers, as well as the well-known SIA and SSA solvers as boundary cases. We investigate the numerical stability of these solvers as a function of grid resolution and the state of the ice sheet for an explicit time discretization scheme of the mass conservation step. Under simplified conditions with constant viscosity, the maximum stable time step of the Hybrid solver, like the SIA solver, has a quadratic dependence on grid resolution. In contrast, the DIVA solver has a maximum time step that is independent of resolution as the grid becomes increasingly refined, like the SSA solver. A simple 1D implementation of the L1L2-SIA solver indicates that it should behave similarly, but in practice, the complexity of its implementation appears to restrict its stability. In realistic simulations of the Greenland Ice Sheet with a nonlinear rheology, the DIVA and SSA solvers maintain superior numerical stability, while the SIA, Hybrid and L1L2-SIA solvers show markedly poorer performance. At a grid resolution of Δx=4 km, the DIVA solver runs approximately 20 times faster than the Hybrid and L1L2-SIA solvers as a result of a larger stable time step. Our analysis shows that as resolution increases, the ice-dynamics solver can act as a bottleneck to model performance. The DIVA solver emerges as a clear outlier in terms of both model performance and its representation of the ice-flow physics itself.

Simulating the Laurentide ice sheet of the Last Glacial Maximum

“…Numerical experiments are conducted with a higher-order three-dimensional ice-sheet model Yelmo (Robinson et al, 2020(Robinson et al, , 2022.…”

Benchmarking the vertically integrated ice-sheet model IMAU-ICE (version 2.0)

“…A comprehensive derivation of these equations is provided by Lipscomb et al (2019). The way the DIVA is solved numerically in IMAU-ICE v2.0 is mostly adopted from Yelmo (Robinson et al, 2022), and is very similar to typical SSA solvers, using an "outer loop" where the effective viscosity, the effective basal friction, and the velocities are iteratively updated until the solution converges. The stress-balance discretization on the staggered Arakawa-C grid is described in Appendix A.…”

“…Goldberg (2011) applied the DIVA to the experiments from the ISMIP-HOM intercomparison (Pattyn et al, 2008) and showed that it produces results that agree with those from higher-order and full-Stokes models down to spatial scales as small as 10 km. Goldberg (2011) also showed that the DIVA was nearly as simple to implement as the hybrid SIA/SSA approximation, and Robinson et al (2022) showed that it was significantly faster than other vertically integrated schemes at high (< 5 km) resolutions due to the superior numerical stability and consequently larger time steps. Despite these advantages, the DIVA is, to our knowledge, currently only used in CISM and Yelmo (Robinson et al, 2022).…”

“…Goldberg (2011) also showed that the DIVA was nearly as simple to implement as the hybrid SIA/SSA approximation, and Robinson et al (2022) showed that it was significantly faster than other vertically integrated schemes at high (< 5 km) resolutions due to the superior numerical stability and consequently larger time steps. Despite these advantages, the DIVA is, to our knowledge, currently only used in CISM and Yelmo (Robinson et al, 2022).…”

Multifidelity Deep Operator Networks

“…The proposed multifidelity DeepOnet framework can be used to leverage the available hierarchy of ice-sheet models [e.g. 15,32,40] of different fidelity and cost, and significantly reduce the cost of generating model data. This hierarchy of models is based on different approximations of the Stokes equation that exploit the shallow nature of ice sheets (see Appendix D).…”