2 Coupled geometry-velocity modeling

Reader familiarity is assumed with the standard SIA and Stokes stress-balance equations (Greve and Blatter, Reference Greve and Blatter2009; Schoof and Hewitt, Reference Schoof and Hewitt2013). 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, Reference Greve and Blatter2009) 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, Reference Robinson, Goldberg and Lipscomb2022), 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, Reference LeVeque2007).

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

(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.Footnote ² 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, Reference Calvo, Diaz, Durany, Schiavi and Vázquez2003; Jouvet and Bueler, Reference Jouvet and Bueler2012; Schoof and Hewitt, Reference Schoof and Hewitt2013), but its role as the determinant of margin location is universal across fluid-layer problems with signed source terms (Bueler, Reference Bueler2021a).

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 2DFootnote ³ 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, Reference Winkelmann2011; Leng and others, Reference Leng, Ju, Gunzburger, Price and Ringler2012; Brinkerhoff and Johnson, Reference Brinkerhoff and Johnson2015; Hoffman and others, Reference Hoffman2018; Aschwanden and others, Reference Aschwanden2019, 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, Reference Brown, Smith and Ahmadia2013; Isaac and others, Reference Isaac, Stadler and Ghattas2015; Tuminaro and others, Reference Tuminaro, Perego, Tezaur, Salinger and Price2016), 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, Reference Hock2005) 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, Reference Hindmarsh and Payne1996; Bueler and others, Reference Bueler, Lingle, Kallen-Brown, Covey and Bowman2005), 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, Reference Winkelmann2011) 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, Reference Fischler2022; Robinson and others, Reference Robinson, Goldberg and Lipscomb2022, e.g.). The optimistic paradigm has theoretical support for a certain higher-order DIVA scheme (Robinson and others, Reference Robinson, Goldberg and Lipscomb2022, Equation (52)), but practical Greenland simulations in the same work actually suggest an intermediate power Δt ≈ O(Δx ^1.6) (Robinson and others, Reference Robinson, Goldberg and Lipscomb2022, 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, Reference Bueler2021b). However, a multigrid method (Trottenberg and others, Reference Trottenberg, Oosterlee and Schuller2001) can greatly reduce α, and the ideal value α = 0 describes an optimal solver in the language of algorithmic scaling or solver complexity (Bueler, Reference Bueler2021b). For example, Antarctic ice-sheet results by Isaac and others (Reference Isaac, Stadler and Ghattas2015), 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, Reference Isaac, Stadler and Ghattas2015, Table 8.1). The Leng and others (Reference Leng, Ju, Gunzburger, Price and Ringler2012) 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 (Reference Tuminaro, Perego, Tezaur, Salinger and Price2016) report total algebraic-multigrid-preconditioned iterations suggesting α ≈ 0.05,Footnote ⁴ and a geometric multigrid method by (Brown and others, Reference Brown, Smith and Ahmadia2013) 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, Reference Leng, Ju, Gunzburger, Price and Ringler2012; Brown and others, Reference Brown, Smith and Ahmadia2013; Isaac and others, Reference Isaac, Stadler and Ghattas2015; Tuminaro and others, Reference Tuminaro, Perego, Tezaur, Salinger and Price2016; Fischler and others, Reference Fischler2022) 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.Footnote ⁵

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 (Reference Bueler2016), 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 (Reference Jouvet and Bueler2012) 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 (Reference Jouvet and Gräser2013) and Brinkerhoff and Johnson (Reference Brinkerhoff and Johnson2015) 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 (Reference Jouvet and Gräser2013) 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, Reference Gudmundsson2013, 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.