Performance analysis of high-resolution ice-sheet simulations

Abstract

Numerical glacier and ice-sheet models compute evolving ice geometry and velocity fields using various stress-balance approximations and boundary conditions. At high spatial resolution, with horizontal mesh/grid resolutions of a few kilometers or smaller, these models usually require time steps shorter than climate-coupling time scales because they update ice thickness after each velocity solution. High-resolution performance is degraded by the stability restrictions of such explicit time-stepping. This short note, which considers the shallow ice approximation and Stokes models as stress-balance end members, clarifies the scaling of numerical model performance by quantifying simulation cost per model year in terms of mesh resolution and the number of degrees of freedom. The performance of current-generation explicit time-stepping models is assessed, and then compared to the prospective performance of implicit schemes. The main results highlight the key roles played by the algorithmic scaling of stress-balance solvers and coupled, implicit-step solvers.

1 Introduction

Numerical ice-sheet and glacier models with evolving ice geometry are now routinely used for addressing scientific questions such as quantification of future sea-level rise from changes in the Antarctic (Seroussi and others, 2020) and Greenlandic (Goelzer and others, 2020) ice sheets, interpretation of the paleoglacial record (Weber and others, 2021) and evaluation of long-term glacial erosion rates (Seguinot and Delaney, 2021), among other applications. In order to resolve ice streams as fluid features it is now generally accepted that valid results need horizontal mesh (grid) cells smaller than ~10 km, but narrow outlet glacier flows and alpine topography need yet finer resolution. Resolving the physics of marine ice sheets also benefits from fine resolution, although careful choice of basal parameterizations will reduce resolution dependence (Gladstone and others, 2017). Thus, whether using local mesh refinement (Hoffman and others, 2018; Fischler and others, 2022, e.g.) or not, resolutions of 2 or 1 km (Seguinot and Delaney, 2021) or less (Clarke and others, 2015; Aschwanden and others, 2019) are increasingly used for science at ice-sheet scale.

Current-generation ice-sheet models apply a variety of stress balances, from the simplest shallow ice approximation (SIA), through ‘hybrid’ (Winkelmann and others, 2011; Robinson and others, 2022) and higher-order balances, up to the non-shallow and non-hydrostatic Stokes approximation. With very few exceptions, however, current time-stepping models alternate between solving the stress balance for velocity, using the geometry determined by the previous time step, and then updating the geometry using the just-computed velocity field. Therefore ice thickness and surface elevation, which are equivalent geometry variables for most modeling purposes, are updated after velocity is fixed. The interaction between the ice sheet and the surrounding climate occurs during this geometry-update operation, via the mass continuity or surface kinematical equations (Greve and Blatter, 2009). Such climatic coupling occurs through surface mass balance, sub-shelf (basal) mass balance and calving processes, in particular.

In other words, these current-generation schemes implement explicit time-stepping for the coupled mass and momentum system describing the dynamical evolution of ice sheets.¹ However, at least for simpler, textbook partial differential equation problems, the limited and conditional stability of such explicit time-stepping schemes is well-understood (LeVeque, 2007). Stability conditions of explicit SIA models appeared some time ago (e.g. Hindmarsh and Payne, 1996), but recent studies have focussed on the stability conditions of explicit hybrid, higher-order and Stokes dynamics models (Cheng and others, 2017; Robinson and others, 2022), or on lengthening their steps (Löfgren and others, 2022).

However, actual implicit time-stepping (LeVeque, 2007) should also be considered. Here the velocity and geometry are updated simultaneously by solving coupled mass and momentum conservation equations. While such implicit time-stepping requires the solution of systems of equations at each step, for glaciers and ice sheets an implicit step must simultaneously compute the velocity field and the domain on which the velocity is defined, namely the 3D extent of the ice once the coupled solution has converged. The problem is of free-boundary type, especially in map-plane (horizontal) directions (Schoof and Hewitt, 2013), but also in the more easily resolved vertical direction.

An implicit strategy has been demonstrated at high (900 m) resolution for the Greenlandic ice sheet in the simplest frozen-base, isothermal SIA case (Bueler, 2016). That work also shows how the steady-state SIA problem (c.f. Jouvet and Bueler, 2012) can be solved, demonstrating unconditional stability. (Observe that steady-state equations correspond to an implicit time step of infinite duration.)

The corresponding problem for the Stokes equations has not, to the author's knowledge, yet been attempted. However, important early work applying a semi-implicit time step using Stokes dynamics (Wirbel and Jarosch, 2020) to solve a free boundary problem illuminates some of the techniques and difficulties needed to make such a strategy work for a membrane-stress-resolving balance.

This is the context in which the current note relates time-stepping and stress-balance solver choices to computational effort. The simplified performance analysis here exposes the most important considerations and trade-offs. While this author expects that implicit time-stepping ice-sheet models will eventually be the fastest, and that advanced solver techniques like multigrid (Briggs and others, 2000) will lead to better modeling, such beliefs should be assessed quantitatively to the extent possible.

2 Coupled geometry-velocity modeling

Reader familiarity is assumed with the standard SIA and Stokes stress-balance equations (Greve and Blatter, 2009; Schoof and Hewitt, 2013). These (continuum) models are regarded here as end members of current-usage stress-balance approximations. Familiarity with the mass continuity and surface kinematical equations (Greve and Blatter, 2009) is also assumed. However, before analyzing the performance consequences of numerical modeling choices, the form of the mass continuity equation is examined, and then the terms ‘explicit’ and ‘implicit’ are carefully defined.

For an incompressible ice sheet with thickness H(t,  x), vertically averaged horizontal velocity u(t,  x) and climatic-basal mass balance a(t,  x), the mass continuity equation says

(1)$${\partial H\over \partial t} + \nabla_{\bf x} \cdot \left({\bf u} H\right) = a.$$

(Here x = (x,  y) denotes the horizontal coordinates.) Equation (1) suggests that ice sheets change geometry in an essentially advective manner, but this appearance is deceiving, or at least over-simplified, especially regarding the growth of numerical instabilities. This is because ice flows dominantly downhill. Indeed, ice-sheet flow has no characteristic curves, as would Eqn (1) if it were a true advection, because the velocity u actually depends on the gradient of thickness through the stress balance.

Thus, as can be addressed by linearized analysis (Robinson and others, 2022), when thickness perturbations grow unstably under explicit time-stepping, i.e. with too large a step, they do so by a mix of advective and diffusive mechanisms. Let s(t,  x) = H(t,  x) + b(x) denote the surface elevation, for bed elevation b(x). A numerical thickness perturbation will often cause the velocity u to respond by increasing in a direction close to downhill ($-\nabla _{\bf x} s$), a direction correlated to $-\nabla _{\bf x} H$ over large areas of an ice sheet. In membrane-stress-resolving models like Stokes this happens through the non-local solution of the stress balance, in which the gravitational source term effectively acts along the surface gradient $\nabla _{\bf x} s$ because of the ice geometry. One might write that Eqn (1) has velocity ${\bf u}( H,\; \, \nabla _{\bf x} s)$, a non-local function of geometry, but a stress-balance solution is required to evaluate this function. In any case, a numerical instability of Eqn (1) occurs when the ice thickness under/over-shoots its correct value because the numerical velocity from evaluating this non-local function is too strong for the given time step (LeVeque, 2007).

A diffusive description of the mass continuity equation is also valid in the small-aspect-ratio limit which generates the SIA (Schoof and Hewitt, 2013):

(2)$${\partial H\over \partial t} = \nabla_{\bf x} \cdot \left(d\, \nabla_{\bf x} s \right) + a.$$

Here $d = C H^{\text {n} + 2} \vert \nabla _{\bf x} s\vert ^{\text {n}-1}$ is the nonlinear diffusivity.² While Eqn (2) does not hold directly for Stokes or other membrane-stress-resolving dynamics, the same diffusivity d, an essentially geometric quantity, can be computed regardless. Generically across stress-balance choices, for grounded glaciers and ice sheets one will observe that large values of d indicate locations of unstable mode growth if explicit time steps are chosen too large.

However, either statement (1) or (2) of the mass continuity equation is fundamentally incomplete without modeling the evolving (map-plane) boundary position of the glacier or ice sheet. The problem of determining this updated margin position is only well-posed via an inequality constraint, namely that ice thickness is non-negative (H ≥ 0), equivalently that the ice surface elevation equals or exceeds the bed elevation (s ≥ b). The mathematical role of this constraint has been understood for some time in the SIA model (Calvo and others, 2003; Jouvet and Bueler, 2012; Schoof and Hewitt, 2013), but its role as the determinant of margin location is universal across fluid-layer problems with signed source terms (Bueler, 2021a).

In this context one can distinguish time-stepping types in coupled geometry–velocity, i.e. mass and momentum conserving, glacier and ice-sheet models. A coupled model computes (fully) implicit time steps if no significant aspects of the geometry or velocity are held fixed during the (coupled) solution. That is, in an implicit scheme a coupled and well-posed model for the ice surface elevation (or thickness) and velocity field updates is solved. A scheme which is implicit in the sense here can be unconditionally stable; it can stably compute time steps of arbitrary duration causing non-trivial changes in margin position. By contrast, a scheme which does not allow the map-plane ice-covered region to change during the time step cannot respond physically (i.e. conservatively with respect to mass and momentum conservation) to a change in climate inputs. We say that a coupled scheme is semi-implicit if, in particular, aspects of the ice geometry are held fixed during the geometry–velocity solution, or, for instance, if the velocity is held fixed during a time step as the geometry is updated (though perhaps ‘implicitly’ on its own). Thus a scheme which computes margin advance or retreat only after the velocity field update is accepted is only semi-implicit, for example. Finally, a scheme is (fully) explicit for the coupled problem if the ice geometry is held fixed during the velocity solution, and then this (accepted and now fixed) velocity solution is used to update the geometry through an explicit step for Eqn (1) or (2).

Thus a key point about time-stepping schemes for the coupled geometry and velocity problem for glaciers and ice sheets is that implicit schemes are not solving fixed-domain systems of coupled PDEs. The mathematical problem for the coupled and implicit time step includes the condition which controls the moving and free boundary, namely the inequality constraint of non-negative thickness. Implicit schemes which solve these free-boundary problems at each time step can, at high spatial resolution, transcend stability limitations on time-step duration. They, and only they, can perform time-stepping at a rate controlled only by ice–climate interaction time scales.

3 Performance analysis

Table 1 lists the parameters used in the numerical model performance analysis which follows. The primary parameters are Δx, which is a representative value for the horizontal mesh (grid) cell diameter, and m, the number of nodes (vertices) in the horizontal mesh. High resolution refers to the equivalent Δx → 0 and m → ∞ limits; asymptotic and big-O notation is used only in this limit.

Table 1. Parameters for performance analysis

Note α, β, γ and m are pure, unit-less numbers.

Assuming the map-plane extent of the model domain is 2D³ and of width L, these primary parameters are related by

(3)$$\Delta x \sim {L\over m^{1/2}} \quad \text{and} \quad {\it m} \sim {{\it L}^2\over \Delta {\it x}^2}.$$

Specifically, for fixed domain width L there are O(m ^1/2) mesh cells in each horizontal dimension.

It is assumed that a numerical glacier or ice-sheet model will use m ice thickness or surface elevation variables, i.e. one degree of freedom per mesh node, including at ice-free nodes where the thickness has value zero. Storing these model state variables, plus the thermodynamical state, requires O(m) memory if the mesh/grid has a priori bounded resolution in the vertical direction. The amount of computer memory needed by the simulation is also O(m); this assumes that prior states are discarded.

Such models also have O(m) velocity variables, but note that these are not state variables, essentially because the Stokes model lacks time derivatives. That is, a very-viscous stress balance computes velocity as a function of the true state variables, such as ice thickness and temperature or enthalpy.

The above assumption of fixed vertical resolution reflects common usage (Winkelmann and others, 2011; Leng and others, 2012; Brinkerhoff and Johnson, 2015; Hoffman and others, 2018; Aschwanden and others, 2019, e.g.), and it permits a rational comparison of asymptotics, but it is not the only possibility. Some solvers use 3D refinement (Brown and others, 2013; Isaac and others, 2015; Tuminaro and others, 2016), with various distinctions between how horizontal and vertical meshing is handled, but this and many other details cannot be pursued here.

Glacier and ice-sheet models resolve and integrate ice–climate interactions, especially via surface mass balance, on time scales which are dominated by an annual cycle, and on longer scales. Let q be the number of ice-dynamical time steps per year needed to capture this coupling. Typical values q = 0.1 a⁻¹,  1 a⁻¹,  12 a⁻¹ correspond to decadal, yearly and monthly frequency, respectively. Note that energy balance and degree-day schemes for computing surface mass balance (Hock, 2005) generally have much shorter time scales, but here q describes the coupling frequency on which ice geometry needs to be updated via solution of the mass continuity or surface kinematical equation.

Current-technology glacier and ice-sheet models use explicit and semi-implicit time-stepping which is only conditionally stable. For the spatial resolutions used in present-day scientific applications, maintenance of time-stepping stability requires steps substantially shorter than 1/q model years. That is, the performance of current-generation models is limited by numerical analysis choices, and not by scientific needs.

For an explicit SIA model the well-known stability restriction is Δt < O(D ⁻¹Δx ²) (Hindmarsh and Payne, 1996; Bueler and others, 2005), where D is a representative diffusivity value, i.e. of d in Eqn (2). For Stokes models, or other membrane-stress-resolving dynamics, the stability of explicit time-stepping is largely unexplored in any precise sense, but an advective restriction Δt < O(U ⁻¹Δx), for some representative horizontal velocity scale U, might be said to represent the optimistic paradigm. The corresponding pessimistic paradigm for Stokes models says Δt < O(D ⁻¹Δx ²), using a representative diffusivity value D computed in the SIA manner.

Explicit time-stepping with hybrid and higher-order schemes is somewhat better studied than for Stokes dynamics, especially over horizontal resolutions relevant to whole ice sheets. Some hybrid schemes apply the pessimistic paradigm as an adaptive restriction (Winkelmann and others, 2011) across all spatial scales. Other models apparently require the user to choose a fixed time step through trial and error in some circumstances (Fischler and others, 2022; Robinson and others, 2022, e.g.). The optimistic paradigm has theoretical support for a certain higher-order DIVA scheme (Robinson and others, 2022, Equation (52)), but practical Greenland simulations in the same work actually suggest an intermediate power Δt ≈ O(Δx ^1.6) (Robinson and others, 2022, Figure 3(a)). An intermediate power, a restriction $\Delta t = O( \Delta x^\omega )$ for some fuzzy value 1.5 < ω < 2, also aligns with this author's experience.

Unconditionally stable implicit schemes also have a maximum time-step restriction, namely Δt < O(q ⁻¹), but this restriction reflects the simulation purpose, not maintenance of stability. For an implicit scheme the desired frequency of ice–climate interaction will, however, determine the total simulation cost according to the (large) solution cost, at each time step, of solving a coupled mass and momentum-free boundary problem, including the surface kinematical or mass continuity equation.

For each explicit time step of a model which applies a membrane-stress-resolving balance, the computational cost of a velocity (or velocity/pressure) solution of the stress-balance equations is determined by solver design. Let us assume that one such solution requires O(m ^1+α) floating point operations (flops), with the power α ≥ 0 depending on the solver implementation. For example, a Stokes solver using direct linear algebra for each Newton step might yield α ≈ 1 if sparsity is exploited or α ≈ 2 if not (Bueler, 2021b). However, a multigrid method (Trottenberg and others, 2001) can greatly reduce α, and the ideal value α = 0 describes an optimal solver in the language of algorithmic scaling or solver complexity (Bueler, 2021b). For example, Antarctic ice-sheet results by Isaac and others (2015), for a Stokes solver implemented using algebraic multigrid, show that the total number of preconditioned Krylov iterations, over the non-linear solve, grows slowly under mesh refinement, suggesting perhaps α ≈ 0.2 (Isaac and others, 2015, Table 8.1). The Leng and others (2012) algebraic-multigrid Stokes solver may have similar scaling, but reported results do not constrain α. For a higher-order stress balance on the Greenlandic ice sheet, Tuminaro and others (2016) report total algebraic-multigrid-preconditioned iterations suggesting α ≈ 0.05,⁴ and a geometric multigrid method by (Brown and others, 2013) suggests α is close to zero for simplified geometries. Observe that the SIA velocity computation, a trivialization of the Stokes problem in which velocity is computed by a pointwise formula, requires optimal O(m) flops.

When analyzing solver scaling in our simplified form, one must be aware that the constant in ‘O(m ^1+α)’ above can be very large, and it depends strongly on solver design. Furthermore, many considerations are suppressed in any flops-based analysis of algorithmic scaling, as actual run time is also determined by memory latency, memory bandwidth and process/thread/GPU parallelism, among other factors. Indeed, many models (e.g. Leng and others, 2012; Brown and others, 2013; Isaac and others, 2015; Tuminaro and others, 2016; Fischler and others, 2022) also show good parallel scaling, something not addressed here.

Regardless of the stress balance, an explicit time-stepping scheme will first compute velocity from current geometry and then apply the mass continuity equation to update the ice thickness using O(m) work. That is, once the velocity is computed from the last-known geometry, let us assume that an explicit scheme replaces old thickness values by new ones using an optimal pointwise formula.⁵

Now, how many flops are needed to simulate one model year? Suppose a numerical model takes time steps of Δt model years duration, equivalently Δt ⁻¹ steps per model year. In these terms stability restrictions will give required numbers of steps per model year. That is, stability requires that Δt ⁻¹ be bounded below by a function of the horizontal resolution Δx or the degrees of freedom m. Recalling the scaling in Eqn (3), the explicit SIA and pessimistic-Stokes cases have

(4)$${1\over \Delta t} > \displaystyle O\left({D\over \Delta x^2}\right) = \displaystyle O\left({D m\over L^2}\right).$$

The explicit, optimistic-Stokes estimate becomes

(5)$${1\over \Delta t} > \displaystyle O\left({U\over \Delta x}\right) = \displaystyle O\left({U m^{1/2}\over L}\right).$$

The number of time steps per model year is then multiplied by the per-step computational cost, namely O(m ^1+α) for Stokes models and O(m) for SIA models, to give a work estimate for each model year in a simulation. The results so far are shown in the ‘explicit’ rows of Table 2.

Table 2. Asymptotic estimates of algorithmic scaling, measured by floating point operations per model year, for map-plane (2D) time-stepping numerical ice-sheet simulations, in the high-resolution limit where Δx → 0 and m → ∞.

See Table 1 for notation.

As already explained, unconditionally stable implicit methods may have a fixed time step Δt = 1/q, independent of Δx and determined only by the need to resolve ice–climate interactions. On the other hand, the per-step expense is much greater because non-trivial coupled equations, and a free-boundary problem, must be solved simultaneously for the updated velocity and geometry.

For SIA models, let us assume that the flops of such coupled solutions scale as O(m ^1+β) with β ≥ 0; a large constant is assumed. The only currently implemented, unconditionally stable, fully implicit time-step solvers use the SIA stress balance. For simplified (dome) geometry the scheme in Bueler (2016), based on Newton steps solved via direct linear algebra and single-grid ice margin determination, directly computes steady states with scaling β = 0.8. (An earlier steady, Picard iteration implementation by Jouvet and Bueler (2012) scales worse.) Convergence is robust for implicit time steps of years to centuries on kilometer-scale grids for the Greenland ice sheet, using realistic and irregular bed topographies.

Certain recent numerical models are semi-implicit in innovative ways. The SIA portions of hybrid time-stepping schemes by Jouvet and Gräser (2013) and Brinkerhoff and Johnson (2015) are solved implicitly, with non-negativity of thickness in their SIA portions enforced as part of a free-boundary solution. However, the time-stepping in these models is only semi-implicit because the membrane-stress-resolving portion of the velocity is held fixed as the ice thickness is advected. The SIA time-step solvers in these schemes have apparently only been tested for time steps satisfying an advective condition Δt < O(U ⁻¹Δx), for U derived from sliding speeds. (Also, Jouvet and Gräser (2013) use multigrid methods, but their published results do not constrain the power β.) The semi-implicit scheme in the open-source Úa numerical marine ice-sheet model (Gudmundsson, 2013, see also G. H. Gudmundsson, Úa Compendium, github.com/GHilmarG/UaSource) implicitly solves a mass continuity free-boundary problem for Eqn (1) using an active set method, but also while holding the (membrane-stress-resolved) velocity fixed.

For implicit Stokes time-stepping, a coupled and free-boundary velocity and geometry-update solve is assumed to be, in the absence of constraining research, O(m ^1+γ) for some γ ≥ α to be determined. One might also suppose γ ≥ β, but there are no implemented cases to measure. These comments complete Table 2. Note that all implicit scheme estimates involve an especially large scheme-dependent constant.

4 Discussion and conclusion

From Table 2 one first observes a well-known property of explicit time-stepping for 2D (map-plane) diffusion equations such as SIA equation (2), namely that computational effort, here flops per model year, scales as O(Δx ⁻⁴). This follows because Δt < O(Δx ²), and because the expense of one geometry-update operation is O(m) = O(Δx ⁻²). Spatial mesh refinement by a factor of two therefore imposes an impressive 16-times increase in effort.

The scaling of the Bueler (2016) implicit SIA solver is not enough better than such an explicit scheme, however. Although long time steps can be taken by this implicit solver, the β = 0.8 scaling gives O(Δx ^−3.6) effort. However, an improvement to β < 0.5, presumably by application of a multigrid method, would transform such an implicit solver into a tool with notably superior O(Δx ⁻³) performance or better.

Now bypassing the over-simplified SIA stress balance, we see from Table 2 how the algorithmic scaling of membrane-stress-resolving solvers dominates performance concerns. In particular, an α ≈ 1 explicit Stokes model, e.g. one using direct, sparsity-exploiting linear algebra on each Newton step system, will do work which scales at the horrific rate O(Δx ⁻⁶) under a pessimistic stability condition. Optimistic stability yields a still-bad O(Δx ⁻⁵) rate. These simple observations emphasize the key role which will be played by multigrid-based stress-balance solvers. That is, even retaining explicit time-stepping, nearly optimal (α ≈ 0) scaling of computation effort in Stokes and higher-order solvers will be necessary for routine application on high-resolution meshes.

On the other hand, suppose resolution (Δx) is fixed. The Table also shows why algorithmic scaling remains important in the large-domain L → ∞ limit when applying Stokes or higher-order dynamics. The computational work of a nearly optimal α ≈ 0 solver will be proportional to ice-sheet area L ². However, an α ≈ 1 method which might suffice for a smaller L = 100 km ice cap will struggle for a L = 1000 km ice sheet because the effort scales as the fourth (2 + 2α ≈ 4) power of L.

The promise of nearly optimal solvers is truly revealed, however, when implicit, coupled geometry–velocity updates are considered. If they become possible, methods which simultaneously satisfy the mass and momentum equations at each time step, doing work (essentially) proportional to the number of degrees of freedom, have great promise. A γ ≈ 0 implicit Stokes method would be capable of many new tasks. More achievably, a γ < 0.5 implicit, essentially unconditionally stable, Stokes time-stepping method, presumably based on multigrid solution of the free-boundary problem for the coupled mass and momentum equations, is an appropriate goal for coming decades of research on numerical ice-sheet models. The computation cost would scale at O(Δx ⁻³), better than explicit SIA models, and the method would only update ice geometry when needed by ice–climate coupling, while avoiding shallow approximations. (The same goal makes sense for any membrane-stress-resolving solver.) Little technical progress has yet been made on such a coupled, fully implicit and scalable Stokes design (but see Wirbel and Jarosch, 2020), so these aspirations are decidedly long term. However, the above discussion suggests why measured values for α,  β,  γ, or equivalent algorithmic scaling measures, are important performance metrics to report when describing new ice-sheet solvers.