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

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 depthintegrated 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 5 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. Under simplified conditions with constant viscosity, the maximum stable timestep of the Hybrid solver, like the SIA solver, has a quadratic dependence on grid resolution. In contrast, the DIVA solver has a maximum timestep that is independent of resolution as the grid become increasingly refined, like the SSA solver. Analysis indicates that the L1L2-SIA solver should behave similarly, but in practice, the complexity of its implementation can make it difficult to 10 maintain stability. In realistic simulations of the Greenland ice sheet with a non-linear 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 15 times faster than the Hybrid and L1L2-SIA solvers as a result of a larger stable timestep. 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 15 representation of the ice-flow physics itself.


Introduction
Modeling ice sheets at the continental scale requires compromise. The Greenland and Antarctic ice sheets span thousands of kilometers, and interact with the rest of the Earth system on timescales of up to thousands of years. Solving the full Stokes stress balance to calculate the ice dynamics over such a large domain for a long timescale is still computationally infeasible. 20 Therefore, it is important to reduce the complexity of simulations of ice dynamics. mixed, and the Hybrid approach should provide a smooth transition. The Hybrid solver is used by several models today (e.g., Winkelmann et al., 2011;Pattyn, 2017;Quiquet et al., 2018;Robinson et al., 2020). In fact, roughly half the models participating in ISMIP6-Greenland and ISMIP6-Antarctica used a Hybrid solver (Goelzer et al., 2018;Seroussi et al., 2019;Goelzer et al., 2020;Seroussi et al., 2020). 25 Others have rigorously derived approximations based on variational principles that intrinsically account for both shear and stretching (Hindmarsh, 2004;Schoof and Hindmarsh, 2010;Goldberg, 2011). The depth-integrated approximations are especially appealing, since they are comparable in speed per timestep to the SSA and Hybrid solvers. The depth-integrated solver derived by Schoof and Hindmarsh (2010) and formalized by Perego et al. (2012) is part of the L1L2 family of solvers, following the terminology of Hindmarsh (2004). We refer to it here as the L1L2-SIA solver. The stress balance is solved at one 30 layer in the ice sheet -in this case, at the bed -while the effective viscosity accounts for stress arising from both shearing and stretching. Thus the solver enables a 2D solution of the system of partial differential equations (PDEs) and fast performance compared to 3D solvers. This solver naturally incorporates both shearing and sliding regimes while approximating the shear stress components via the SIA, which facilitates its numerical solution. The L1L2-SIA solver is used by the BISICLES model (Cornford et al., 2013). a motivation that contributed to the development of the Hybrid, L1L2-SIA, and DIVA solvers. Today, however, the community is focused on running such simulations at high resolutions (i.e., ∆x < 5 km), as the tools and computational resources become available.
Given the importance of computational efficiency in modelling glacial dynamics at the continental scale, the goal of this paper is to investigate the numerical performance of the solvers described above, when coupled to mass-conservation timestep- 15 ping. We focus on the latter three solvers, which account for the dominant stress terms over most of an ice sheet, while the SIA and SSA solvers are useful boundary cases. Below, we first outline the numerical approach used by each of these solvers.
Next, we derive stable timestep limits of ice-thickness advection in an analytical test case for the DIVA, Hybrid and L1L2-SIA solvers, and we investigate the underlying factors that can affect the maximum stable timestep given each solver. The analytical results are validated in idealized numerical tests. We then compare the performance of all five solvers in terms of stable timestep 20 size and model computational speed in a realistic test case of quasi steady-state simulations of the Greenland ice sheet. This is followed by a discussion of the results and conclusions.

Ice dynamics solvers
We describe the assumptions and equations behind the five solvers considered here, namely: the SIA, SSA, DIVA, Hybrid and L1L2-SIA solvers. These approximations can be obtained by considering the various terms in the first-order Blatter-Pattyn (BP) 25 approximation (Blatter, 1995;Pattyn, 2003). To give context to the depth-integrated solvers, we first write the basic equations of the 3D BP stress-balance approximation: where u and v are the components of horizontal velocity, µ is the effective viscosity, ρ i is the density of ice (assumed constant), g is gravitational acceleration, s is the surface elevation, and x, y, z are 3D Cartesian coordinates. In each equation the three terms on the left-hand side (LHS) describe gradients of longitudinal stress, lateral shear stress, and vertical shear stress, respectively, and the right-hand side (RHS) is the gravitational driving force. The longitudinal and lateral shear stresses together are also referred to as membrane stresses (Hindmarsh, 2006).
The effective viscosity µ is defined as 5 where A is the temperature-dependent rate factor in Glen's flow law (Glen, 1955), n is the Glen exponent (many models set n = 3), andε e is the effective strain rate, given in the BP approximation bẏ Here, u x denotes the partial derivative ∂u/∂x, and similarly for other derivatives.

SIA solver 10
The shallow-ice approximation (SIA) solver is a zero-order approximation to ice flow that assumes a balance between the basal drag and the gravitational driving stress (Greve and Blatter, 2009). This leads to purely shear-stress-driven flow within the ice column. In other words, the membrane stress gradients are ignored, which leads to a stress balance equation that can be solved locally. This is a typical flow regime for grounded ice that is not sliding.
By setting membrane stress gradients to zero in Eq.
(2) is obtained using an effective strain rate that includes only the shear terms: By assuming a low aspect ratio and a free-stress surface boundary condition, these equations can readily be integrated from the bed b to the surface s to provide a 2D solution for depth-average ice flow: whereū andv are 2D vertically averaged velocity components, andμ is the vertically averaged effective viscosity: 2.2 SSA solver The shallow-shelf approximation (SSA) solver is complementary to the SIA solver in that the membrane stress gradients are retained, while the vertical shear stress gradients are ignored. This is a typical flow regime for floating ice or rapidly sliding ice streams, where the driving stress is low and the ice column has a uniform velocity profile (i.e., plug flow). By dropping the vertical shear-stress gradient terms in Eq.
(1) and integrating vertically while assuming a vertically uniform velocity, we obtain 5 the 2D SSA stress balance equations: where τ b,i indicates the basal stress in the direction i. We assume a basal friction law of the form where τ b is the basal shear stress, is the basal velocity, and β is a non-negative friction coefficient that can 10 depend on u b . The effective viscosity is calculated following Eq. (2) for the BP solver, using a vertically averaged rate factor and computing the effective strain rate in 2D without the shear terms:

DIVA solver
Goldberg (2011) derived a higher-order stress approximation that is similar in accuracy to the first-order BP stress balance, but 15 like other depth-integrated solvers is computationally much cheaper. Since the stress-balance equations use a depth-integrated effective viscosity in place of a vertically varying viscosity, we refer to this scheme as a depth-integrated-viscosity approximation, or DIVA.
In the DIVA solver, the horizontal velocity gradients and effective viscosity in the membrane stress terms are replaced by vertical averagesū,v, andμ (Eq. (7)). The 3D effective viscosity is given by Eq.
(2), but with the DIVA effective strain rate: The vertically integrated stress balance in DIVA can be written as where the boundary conditions at b and s have been used to evaluate the vertical stress terms, and the basal stress is defined as for the SSA solver. Equation (12) has the same form as the SSA stress balance, but the dependence of basal stress and viscosity 25 on velocity is more complex. To solve Eq. (12) for the mean horizontal velocity components, the basal stress terms must be written in terms ofū andv.
The following discussion refers to the u component; the v component is analogous. Goldberg (2011) showed that the vertical profile of u is related to u b by Following Arthern et al. (2015), we define some generalized integrals F m for clarity: The surface velocity is thus related to the bed velocity as Integrating u(z) from the bed to the surface gives the depth-averaged mean velocityū: 10 which allows βu b in Eq. (9) to be replaced with β effū , where For a frozen bed (u b = 0, with nonzero basal stress τ b,x ), the βu b term on the RHS of Eq. (13) is replaced by τ b,x , leading tō Then the basal stress term βu b in Eq. (9) can be replaced by β frz effū , where The velocity is found in two steps. First, Eq. (12) is solved for the mean velocity, using Eqs. (9), (17), and (19) to write the basal stress terms in terms ofū andv. Then Eq. (13) is integrated vertically to find the 3D velocity.

Hybrid solver
The Hybrid solver, as defined here, follows the approach of Winkelmann et al. (2011). The horizontal velocity is defined as the 20 sum of the depth-dependent internal shear velocity (u i , v i ) and the basal sliding velocity (u b , v b ): The sliding velocity is calculated via the SSA and the internal shear velocity is calculated via the SIA, as defined above. This approach generally works well because over most regions of an ice sheet, either plug flow (sliding) or shear flow is dominant, so one approximation is sufficient to represent the flow, and the other term goes to zero (Winkelmann et al., 2011;Pollard and DeConto, 2012). Although this approach is used widely with success at reproducing observed ice flow (e.g. Winkelmann et al., 2011;Quiquet et al., 2018), it is less clear how well it resolves the transition between the two regimes.

L1L2-SIA solver
The L1L2-SIA solver as defined here follows from Schoof and Hindmarsh (2010) and Perego et al. (2012). Like the DIVA and Hybrid balances, it is a two-step approach that begins by solving a depth-averaged stress balance. The balance is similar in in form to the SSA (Eq. (8)): 10 Aside from the fact that the equations are solved for basal, not depth-uniform, velocity, the balance differs from the SSA in the form of the effective viscosity.μ is a vertical average of µ , the L1L2 effective viscosity, which is given by where 15 with τ b,ij , i, j = {x, y} the longitudinal stress induced by the basal strain-rate tensor. τ b,e depends implicitly onε b,e , the effective horizontal (2D) strain rate induced by the basal strain-rate tensor (cf Eq. (10)): For n = 3, takingε b,e as a function of u b and v b from the previous iteration, Eq. (24) is a cubic equation for τ b,e and can be solved analytically. Givenμ from Eqs. (22) and (24), we can then solve Eq. (21) for u b and v b . Note that the "trial" viscosity 20 µ makes the approximation that vertical shear stress is equal to the shallow-ice shear stress, −ρ i gH∇s. In the second step, three-dimensional vertical shear stress is diagnosed from the solution: where τ yz is diagnosed with a similar expression. Finally, we can define an updated viscosity: Horizontal velocities are then given by An important point regarding L1L2-SIA is the difference between µ and µ. The former is defined using the SIA shear stress 5 as a proxy for the true vertical shear, while µ is formulated with an updated vertical shear based on the solution to the stress balance. This is in contrast to DIVA, where the viscosity is consistent between the two steps.

Analytical stability limits
In this section, we assess the numerical stability of various schemes to solve the coupled stress-balance and continuity equations. Our analysis is closely related to von Neumann stability analysis (Isaacson and Keller, 2012), but the non-local nature of 10 the equations requires us to consider global solutions to the stress balance. For this analysis, we impose several simplifications: -We reduce the problem to one horizontal dimension by settingv = 0 and ignoring the y derivatives.
-We assume that the effective viscosity µ is spatially uniform. This is equivalent to assuming that the ice is isothermal with n = 1 in Eq.

15
-We consider perturbations to an idealised geometry: an x-periodic domain with uniform thickness, H 0 , and surface slope, This problem is described by Dukowicz (2012), who derived exact solutions for the velocity profile as a function of α and the non-dimensional constant η = βH/µ.
We specifically analyze the DIVA, Hybrid and L1L2-SIA solvers. We begin with DIVA, which leads to the simplest expres-20 sions for stability in this framework; results for the Hybrid and L1L2-SIA balances follow. The SSA analysis arises naturally, since the SSA equations have the same form as the DIVA solver. The numerical stability of an SIA solver under the above assumptions has been done previously (Hindmarsh, 2001;Cheng et al., 2017), showing that the maximum stable timestep under the simplifications above is proportional to the square of the grid resolution: ∆t < 3µ 2ρ i gH 3 ∆x 2 (29) 25 We caution that the results depend on details of the numerical schemes and may not apply to all situations, such as when the rheology is nonlinear. Our aim is not to consider all possible schemes, or to produce a numerical scheme of optimal stability as in Cheng et al. (2017), but rather to examine stability properties when applying a representative finite-difference scheme to different stress-balance equations. Our analytical results provide context for the empirical timestep limits found in more realistic settings in Section 4, and give a theoretical basis for those results. Readers interested primarily in the performance of various depth-integrated stress balances in continental-scale ice sheet models may wish to skip this section.
3.1 DIVA solver stability 5 We start with the x component of the DIVA stress balance, Eq. (12). By assuming a constant slope α and spatially uniform effective viscosity µ and by settingv = 0, this equation can be written as a second-order ordinary differential equation (ODE) forū: The first term on the LHS is the longitudinal stress, and the RHS is the driving stress. In the second term on the LHS, the 10 quantity in parentheses is β eff , the effective basal friction coefficient. This coefficient includes a non-dimensional parameter, η ≡ βH/µ, that relates basal friction to ice viscosity. If η is small, we have β eff ≈ β, and the flow is dominated by basal sliding.
If η is large, then β eff ≈ 3µ/H, and most of the flow is internal shearing. For the remainder of the analysis, β is considered uniform.
To analyze the stability of a simple numerical scheme for the DIVA equations, we consider a reference state with uniform 15 thickness and surface slope as described above. This reference state has uniformū and hence zero longitudinal stress, with depth-mean velocity The thickness evolution equation is written as 20 We will consider the evolution of thickness under an initial sinusoidal perturbation: Following Cheng et al. (2017), we approximate the DIVA momentum-mass balance, Eqs. (30) and (32), to first order in ε h : where δū(x) is a perturbation in the mean velocity.
Next, we define a simple finite-difference scheme to solve this system of equations. While the theoretical domain has infinite extent, the computational domain must be finite, and we consider a periodic domain of length L. The discretised perturbation thickness, h n , is located at points x = n∆x, for n = 0,..,N (where ∆x = L/N ), and the discretised perturbation velocity, u n , lies at points (n + 1 2 )∆x, for n = 0,..,N − 1. We assume that L is much larger than the parameter L m , defined as As shown in Appendix A1, L m is equivalent to a membrane stress length scale, over which longitudinal stresses are important (Hindmarsh, 2006).
Equation (34) is discretised as where u is the vector of u n ; h is the vector of h n ; h u is the vector of thickness h n averaged to u-points; and ∆ n (·) indicates the 10 first-order finite difference of a discrete function on the numerical grid, i.e.
where 15 To proceed, an expression for u is needed. For a linear stability analysis of SIA schemes (Hindmarsh, 2001;Cheng et al., 2017), this is straightforward, since velocity depends locally on thickness and surface slope. For the DIVA balance, however, a linear system (Eq. 37) must be solved. Analytical solutions to matrix equations are not always available, but for our assumptions (constant viscosity and flow in one horizontal direction) and geometry (periodic with constant surface slope and thickness), an asymptotically accurate closed-form solution can be derived. The mathematical details can be found in Appendices A1 and 20 A2, but here we state the result: where The linearised mass-balance equation, Eq. (35), is discretised with an explicit upwind scheme for the advective term and a centered difference for the divergence term: With h j m = ε h e ikm∆x , and using Eq. (40) for u m , this becomes If the real part of the bracketed expression has an absolute value greater than 1 for an initial condition ε h e ikx , this scheme is unstable.
Writing κ in full, the third term in brackets in Eq. (47) (i.e., the divergence term) can be written 15 Defining θ ≡ k∆x, we seek the value of θ which maximizes the magnitude of the real part of Ξ. In doing so, we effectively ignore the second term in brackets in Eq. (48), which is purely imaginary. This is justified as we are primarily concerned with cases where ∆x → 0, where this term becomes negligible. In Section 3.4, however, we test the validity of the assumption over a wide range of ∆x values. Using Eq. (42) to write φ(∆x, k, r) in full, the remaining term is

20
This expression is non-positive, continuously differentiable in θ, and has extremal points at θ = 0 and θ = π. The maximum value occurs when cos(θ) = cos(π) = −1: We refer to the second term in brackets in Eq. (47) as γ. This term is related to advection and has a maximal value (when Together, γ and Ξ determine the stability of the time-stepping scheme. γ is related to advection, while Ξ arises from dynamic thinning and thickening under divergence. The presumption is that an arbitrary initial condition, unless carefully constructed, 5 projects onto the mode e ikx for which Eqs. (50) and/or (51) are realised. The scheme is stable when the real part of the bracketed expression in Eq. (47) has an absolute value less than 1, or equivalently, when |γ| + |Ξ| < 2. When γ and Ξ are of similar magnitude, the associated limit on ∆t does not have a simple expression. But when one or the other is dominant, there is differing leading-order behaviour in the relationship between resolution and maximal stable timestep.
When γ Ξ, the equation system is essentially advective, and stability is governed by Eq. (51): When Ξ γ, the maximum stable timestep follows from Eq. (50): The dependence of ∆t dyn on ∆x is more complex than that of ∆t adv . It is useful, however, to consider two end-member regimes: large q and small q (in the following, it is assumed unless stated otherwise that ∆t dyn ∆t adv , i.e., advection is not 15 a limiting factor). The case q 1 corresponds to high basal friction and/or coarse resolution. In this limit it can be shown (see If η 1, as in the case of high basal friction, this simplifies to

20
In this case, the stability criterion is identical to that of the SIA solver (Eq. (29)), with the maximum stable timestep, ∆t max , proportional to ∆x 2 . This condition usually is not very restrictive for the DIVA solver, since large q typically implies coarse resolution if the ice is not too thin.
The terms in Eq. (53) containing a 0 and r reduce to 2(2 − √ q) ≈ 4, resulting in 25 ∆t < 8µ Thus, in the limit of small q (or equivalently, small ∆x), the time-step limit arising from ice-flux divergence depends on µ and H 0 , but not ∆x. This is not to say that the limiting timestep is independent of resolution -but such dependence arises from the advective term, and thus the maximum stable timestep varies linearly, not quadratically, with ∆x. For a typical ice thickness of 10 3 m and effective viscosity of 10 7 Pa y, we would have ∆t dyn ∼ 10 y, which is not very restrictive. The maximum timestep would be determined by the advective limit, Eq. (52). With ∆x = 10 3 m and u 0 = 10 3 m yr −1 , the advective limit would be ∆t adv ∼ 1 yr.
These two regimes of limiting timestep correspond to differing behaviour of the DIVA equations at different scales. In 5 the large-q limit, Eq. (54), variations in velocity correlate closely with driving stress, and the dynamic response behaves like a diffusion process (similar to the SIA) in which the flux of thickness is proportional to thickness gradients, leading to a quadratic dependence of the limiting timestep on resolution. In the small-q limit, Eq. (56), small-scale thickness oscillations are damped by dynamic thinning due to velocity divergence induced by the oscillations. The velocity gradients are scale-independent due to averaging of associated driving stresses over the membrane scale, resulting in a scale-independent damping rate, and hence 10 a scale-independent timestep limit (provided the advective timestep limit is large enough).

Hybrid solver stability
We consider now the Hybrid stress balance, in which sliding velocity is determined by the SSA momentum balance, Eq. (8).
Thickness transport is due to a vertically averaged velocityū hyb , the sum of sliding velocity and the vertical average of the SIA velocity from Eq. (6),ū sia = ρ i gH 2 α/(3µ). 15 To investigate the Hybrid solver stability, we will first modify the DIVA analysis (Sect. 3.1) to treat the SSA case. Thus, we treat the SSA balance as a first step in the analysis of the Hybrid scheme. The SSA momentum balance (in one dimension, with constant viscosity) is This is the low-friction limit of the DIVA balance, Eq. 30, with u replacingū since velocity is depth-uniform. The SSA mass 20 balance equation is given by Eq. (32), again replacingū with u. The RHS of the matrix equation for velocity, Eq. (37), is identical, and the matrix C is modified with B = β. Thus, the dynamically-limited timestep ∆t dyn is given by Eq. (53), but with r and a 0 modified by the new definition of B. At high resolution, this leads to the same dynamically-limited behaviour as for DIVA, Eq. (56). At coarse resolution, the dynamically-limited timestep differs from DIVA due to the modified frictional term; we do not give details here.

25
The Hybrid equations for momentum balance and mass conservation in 1D with constant viscosity µ are We consider the same reference state as in the DIVA analysis. The reference sliding velocity is Thus,ū hyb = u 0 , the DIVA depth-averaged reference velocity from Eq. (31). As before (cf. Eqs. (34) and (35)), we consider perturbations in H of ε h e ikx and write the equations to first order in ε h : where Since u eff hyb > u 0 (effectively because there is no '3' in the denominator of the second term in parentheses), the Hybrid advective stability constraint is slightly more restrictive than for DIVA.
Consider the third term on the RHS of Eq. (63), which follows from the dependence of u SIA on perturbations in s. If this term is discretised using an explicit second-order finite difference, it can be expressed as 15 Note that this term is missing from the analogous equation for the DIVA solver (Eq. (35)). With this result and those from the previous section (cf. Eq. (47) for DIVA), it can be shown that the limiting timestep is determined by the constraint that In the above expression, r and a 0 are as defined as in Eqs. (43) and (44), but with the substitution B = β, giving q = β∆x 2 /(4µH 0 ). In the limit q 1 (high basal friction and/or coarse resolution), the third term in Eq. (66)  Recalling that η ≡ βH 0 /µ, the third and fourth terms can be combined to give Provided the advective timestep limit ∆t adv = ∆x/u eff hyb is large, an approximate time step limit in this regime is given by which is identical to the large-q limit for DIVA, Eq. (54), and reduces to Eq. (55) when η 1. When q 1 (low basal friction and/or fine resolution), the third term in (66) is bounded independently of ∆x. For small ∆x, the timestep limit for small q is therefore governed by the fourth term: Since ∆t is proportional to ∆x 2 , Eq. (68) becomes very restrictive at high resolution, like the identical equation for SIA 5 (Eq. (29)).

L1L2-SIA solver stability
As described in Section 2.5, the L1L2-SIA solver consists of a two-step process, the first being a solve for the sliding velocity u b . With a constant viscosity, the resulting equation is identical to that of the Hybrid model, and thus u b is found by solving the SSA equation. The depth-averaged velocity is found by depth-integrating the vertical shear stress. In other words, the 10 longitudinal strain rate is approximated by the gradient of u b , and the resulting x-momentum balance, given by is rearranged and integrated with depth (applying the free-surface boundary condition) to yield This yields the depth-averaged velocity Since the depth-averaged velocities in the Hybrid and L1L2-SIA approximations differ only by a term involving the second derivative of u b , the reference depth-averaged velocity,ū L1L2 , is again equal to u 0 . The linearised version of the perturbed mass balance is given by where the replacement of the terms in parentheses in the second equality comes from Eq. (58) and Functionally, Eq. (72) has the same form as (35). The differences are the prefactors on the gradients of perturbed thickness 25 and perturbed velocity, and the fact that the basal (sliding) velocity is considered, and not depth-averaged velocity. Since u b satisfies the same solution as in the Hybrid balance, the limiting time step is constrained according to Similarly, approximate bounds can be found for the regimes q 1 and q 1. When q 1 the limiting time step can be approximated by This expression is identical to the large-q bounds for DIVA and Hybrid, Eqs. (54), and (67). When q 1, the limiting time step asymptotes to ∆t < 8µ which is smaller than Eq. (56), the small-q limit for DIVA, by a factor related to η. For large η (i.e., high basal friction and low 10 sliding), ∆t is small. But like Eq. (56) for DIVA, Eq. (76) has no resolution dependence. At high resolutions, the maximum stable time step is either independent of ∆x, or, where the advective stability limitation is important, linear in ∆x. Thus the above analysis suggests that, unlike in the Hybrid scheme, the maximum time step does not depend quadratically on ∆x.
As will be shown, however, simulations carried out with realistic models and/or nonlinear rheology suggest that L1L2-SIA stability is similar to that of the Hybrid scheme. We explore this contradiction further in the Discussion section.
15 Table 1 summarizes the low-and high-resolution timestep limits for the DIVA, Hybrid and L1L2-SIA solvers. At low resolution, the limits are identical, since in all three cases the variations in velocity and divergence correlate with those in thickness.
At high resolution, the DIVA solver is not dependent on resolution. In the Hybrid balance, however, the SIA component of depth-averaged velocity correlates with thickness variations. As such, the bound on timestep is still quadratic in grid spacing, but with a larger bounding constant. The L1L2-SIA stability limit is similar to DIVA for high resolution, albeit offset to lower 20 values by the additional term 3 3+η , which is less than 1 and approaches 0 for large η.

Numerical validation
To confirm the validity of the above analysis, we run numerical simulations with 1D DIVA, Hybrid and L1L2-SIA models derived using the assumptions above. For the 1D DIVA solver, we discretise and solve the simplified equations (30) and (32), using a first-order upwind finite volume scheme for Eq. (32). We discretise the Hybrid solver analogously using Eqs. (58) and

25
(59), and L1L2-SIA using Eq. (32) with (71) for depth-averaged velocity. These simplified solvers have been implemented in Matlab and are provided along with code to run the tests cases as a Supplement to the paper. As further validation, we run the same simulations using two comprehensive ice sheet models, Yelmo (Robinson et al., 2020) and CISM (Lipscomb et al., 2019). The three solvers have been implemented in Yelmo (Robinson et al., 2020), and all test cases are simulated. CISM has no Hybrid solver, but is used to test the DIVA and L1L2-SIA solvers.
Solver Low resolution Eq. number High resolution Eq. number Table 1. Summary of asymptotic timestep limits derived for the DIVA, Hybrid, and L1L2-SIA solvers under simplified conditions (1D, uniform viscosity, infinitely long ice slab of uniform thickness H0, explicit time-stepping scheme). Limits are shown for when the advective limit is sufficiently large so as not to apply, and are defined loosely for "low resolution" and "high resolution" regimes. See the text for details.
The test problem is solved with a periodic uniform slab as described in the beginning of this section, with the constant bedrock slope set to α = ds/dx = −10 −3 . For each model and parameter set, we verified that the diagnosed velocity profile (not shown) is consistent with the exact solutions.
To test stability, we add a random Gaussian perturbation to the initial ice thickness, such that h n (t = 0) = H 0 + X n , n = 1, 2, ..N, 5 where X n are independent Gaussian random variables with zero mean and standard deviation of 0.1 m. We then run the model forward for 100 timesteps. While not guaranteed, it is likely that h(t = 0) will project onto the least stable numerical mode (i.e., that for which the expressions defining ∆t dyn and ∆t adv are realised).
For a given set of physics parameters (µ, H 0 , β), ∆x is varied over a range from 10 m to 40 km. At each resolution we run multiple tests with different values of ∆t to determine the maximum stable timestep. For each run, we calculate the ratio σ/σ 0 , 10 where σ is the standard deviation of h at the final timestep, and σ 0 is the standard deviation of h(t = 0). The variance in ice thickness should decrease for a stable scheme, so this ratio serves as a metric of stability. We consider σ/σ 0 ≤ 1 to indicate stability and σ/σ 0 > 1 to indicate instability. The numerical results are compared to the timestep limits determined analytically above.
We test two cases that are representative of different ice-flow regimes. The first parameter set (µ = 1 × 10 5 Pa yr, H 0 = 15 1000 m, β = 1000 Pa m yr −1 , η = 10) corresponds to thicker, less viscous ice with a strong bed, i.e., conditions that favor vertical shear over sliding. In this case, u 0 is equal to 39 m yr −1 , and L m is approximately 1.3 km. With such a low background velocity, the advective timestep limit ∆t adv is not restrictive except at extremely high resolution, and so stability is determined by dynamic divergence alone. The second parameter set (µ = 4 × 10 5 Pa yr, H 0 = 500 m, β = 30 Pa m yr −1 , η = 0.0375) corresponds to thinner, more viscous ice with a weak bed, i.e., conditions that favor fast sliding and vertically uniform flow. In 20 this case, the maximum stable timestep ∆t dyn is generally larger (as a function of ∆x) than in the thick, shearing case. Also, u 0 is larger (150 m/a), which means that ∆t adv may impose a stability limit as ∆x becomes small. Figure 1 shows results for the two parameter sets from the simple 1D models and from Yelmo and CISM, as compared to the analytical solutions.
In the case of the DIVA solver, the maximum stable timesteps determined by the 1D model, Yelmo, and CISM are in excellent agreement with the analytical estimate of ∆t dyn as given by Eq. (50) for both tests. In the shearing case, the advective timestep limit is large enough that it is not relevant, so the maximum stable timestep transitions from a quadratic dependence on grid resolution for low resolutions (∆x > L m ) to a constant value. In the sliding case, the maximum stable timestep follows a similar dependence until ∆x < 10 2 m, at which point the advective limit becomes the limiting factor.

5
For the Hybrid solver, we again find excellent agreement with the analytical estimate of ∆t dyn , as given by Eq. (66), for both the simplified 1D model and the ice sheet model Yelmo. In both tests, the maximum stable timestep follows a quadratic dependence on resolution, although a transition occurs to a slightly more stable regime below ∆x ∼ L m in the sliding case. In the shearing case in particular, the maximum stable timestep drops quickly as ∆x decreases, such that it is already as small as ∆t ∼ 10 −2 when ∆x ∼ L m ∼ 1000 m. For the Hybrid solver, the very small timestep limits ensure that the advective timestep 10 limit never becomes relevant in the two tests.
In the case of the L1L2-SIA solver, the simple 1D model agrees with the analytical results, but the comprehensive models do not. According to Eq. (74), at high resolution the L1L2-SIA solver's maximum stable timestep should become independent of grid resolution. The 1D model confirms this result. The maximum stable timestep in the shearing test at high resolution is a constant value, although offset to a lower value than that of the DIVA solver due to the scaling dependence of the equation on 15 η. Here as well, for extremely high resolution (∆x < 10 m), the model becomes unstable in the shearing case -though this may be due to our use of finite perturbations to validate the linear stability analysis, as smaller perturbations allow further filling of the shaded regions at high resolutions (not shown). In contrast to the 1D model, both Yelmo and CISM exhibit a quadratic dependence on grid resolution akin to that of the Hybrid solver for this test. In the sliding test, the 1D model again agrees with the analytical solution. Yelmo's results again resemble the Hybrid solver here, while CISM maintains a purely quadratic 20 dependence on resolution.
It is not clear why CISM's L1L2-SIA solver is less stable than Yelmo's for the thin sliding case, nor why the theoretical high-resolution stability limit for L1L2-SIA is recovered by the simple 1D model but not by Yelmo and CISM. We point out that the result given by Eq. (76) depends critically on a cancellation of terms in Eq. (72), and it is possible that subtle details of the finite-difference schemes of these more sophisticated schemes prevent such cancellation. Moreover, it is worth noting that 25 this theoretical result also depends on a constant viscosity, which is not the case in realistic models.
Overall, these experiments show that the stability limits derived under simplified assumptions are valid in numerical simulations. Moreover, they highlight the difference in stability regimes between the DIVA, Hybrid and L1L2-SIA solvers. As shown in the analytical derivation, the SSA solver has the same stability limits as the DIVA solver at high resolution, thus it can also be expected to be more stable than the Hybrid solver. In contrast, previous work has shown that SIA stability depends 30 quadratically on ∆x at higher resolutions (Cheng et al., 2017), consistent with the results found for the Hybrid solver that incorporates the SIA solution. While the L1L2-SIA solver appears to have stability limits similar to the DIVA solver in the simplified test setup, the more comprehensive models indicate that the solver may not be so stable in practice. for Hybrid, and Eq. (74) for L1L2-SIA. Light blue shading shows where a 1D model following the derivation in the text has a stable solution.
The green and magenta lines with solid points show the maximum stable timestep as determined using the ice sheet models Yelmo and CISM, respectively. The Hybrid solver was not implemented in CISM, so results are shown only for Yelmo. In some other panels, the CISM and Yelmo results overlap.

Greenland experiments
The above analysis sheds light on the mathematical basis for stability differences between approximations in idealised cases.
We are most interested, however, in comparing the solvers under more realistic conditions with, for example, a nonlinear rheology. To this end, we perform several quasi-steady-state simulations for the Greenland Ice Sheet (GrIS) at different resolutions.
We use the ice sheet model Yelmo (Robinson et al., 2020), which now supports all five solvers (SIA, SSA, Hybrid, L1L2-SIA To save computation time, the simulations are not run fully to steady state. Rather, we run the model for a total of 3 kyr. For the first 2 kyr, we briefly spin up the thermodynamics and basal hydrology and let the model reach a self-consistent state with 5 the topography. We then evaluate the timestep stability during an additional 1 kyr simulation. Yelmo uses an adaptive timestepping scheme that determines the optimal timestep for the current model conditions (Robinson et al., 2020;Cheng et al., 2017). Given a measure of the truncation error τ t in ice thickness at each timestep, the scheme ensures that the maximum value of the truncation error over the domain η t = max (τ t ) remains, on average, within a given tolerance level, ε t . The timestep is reduced when η t > ε t and increased when η t < ε t , leading to a value that is as large as pos-10 sible while maintaining stability. In the following analysis, the mean adaptive timestep over the final 1 kyr of the simulations is used as a metric of the model's computational performance. This metric can be compared to the analytically derived values in the simpler 1D formulation in Section 3.
First, it is instructive to compare the simulated GrIS using each solver, given the same boundary conditions and other model settings (Fig. 2). For the three solvers that include both shear and membrane stresses (Hybrid, L1L2-SIA, and DIVA), the 15 surface velocity fields are similar. All three schemes capture the balance between the gravitational driving stress and vertical shear stress in the interior, between driving stress and membrane stress in ice shelves, and among all three kinds of stresses in fast-sliding ice streams. The ice thickness and velocity distributions are nearly indistinguishable at the regional and large scale.
In contrast, the SIA solver simulates slow inland velocities well, but margin velocities are markedly reduced compared to the solvers with membrane stresses. Since the overall velocity field is dominated by rather slow flow (< 1000 m yr −1 ), the SIA 20 solver generates a reasonable solution for Greenland, but would clearly fail for, e.g., floating ice. Meanwhile, the SSA solver does not account for vertical shear and thus fails to reproduce inland velocities. The SSA solver can be tuned to reproduce present-day velocities well with the proper optimization of basal friction (e.g., Goelzer et al., 2018), but it is not designed for this regime a priori. For purposes of this analysis, the Hybrid, L1L2-SIA, and DIVA solvers produce the most realistic velocity fields, but it is useful to include the SIA and SSA solutions as representative of the two extreme flow regimes (pure shearing 25 and pure sliding).
An analysis of model stability, illustrated in Fig. 3, shows large, systematic differences among solvers. Two solver families emerge: the "SIA" and "SSA" families. For all solvers, the maximum stable timestep decreases nonlinearly with increasing grid resolution (Fig. 3a), as shown by the fitted exponent p in the relation ∆t ∼ ∆x p . Such a relationship could be expected from the stability analysis in Section 3, since the simulations here correspond more closely to the "low-resolution" regime 30 highlighted in Table 1. The three solvers that rely on the SIA in some form (SIA, Hybrid and L1L2-SIA) have similarly high values of p, ≈ 2.5-2.8, and thus a large decrease in timestep with finer resolution. In contrast, the SSA and DIVA solvers have a much weaker dependence on grid resolution (p ≈ 1.5-1.6), resulting in larger stable timesteps. At the coarsest resolution, ∆x = 32 km, all solvers allow a timestep of a similar order of magnitude, between 1-5 yrs, although the SIA family requires timesteps more than two times smaller than the SSA family. As the resolution increases to ∆x = 4 km, the SIA family requires We also tested solver performance for a constant timestep, to ensure that the adaptive timestepping scheme provides good estimates of the maximum stable timestep. We ran simulations using only the Hybrid, L1L2-SIA and DIVA solvers for the GrIS domain at 8 km resolution with fixed timesteps of ∆t = 0.05, 0.1 and 0.5 yr. From the simulations with adaptive timestepping, 5 we expect that the three solvers would be stable at ∆t = 0.05 yr, but that the Hybrid and L1L2-SIA solvers would become unstable for larger timesteps. A DIVA simulation with ∆t = 0.5 yr runs approximately 10 times faster than a simulation with ∆t = 0.05 yr.

10
The difference in maximum stable timestep between solvers results in a marked difference in computational speed ( Fig.   3b). At lower resolutions, the differences are modest. The SSA and DIVA solvers run about twice as fast as the Hybrid and L1L2-SIA solvers, while the SIA solver is on par with the fastest solvers because it is computationally much less intensive. As the resolution increases, the timestep dependency on resolution becomes the limiting factor, and the solvers again separate into the faster SSA family and the slower SIA family.

25
Taking DIVA as a reference, we can benchmark its performance against the other solvers (Fig. 3c). The DIVA solver is comparable in computational performance to the SSA solver at all resolutions, although SSA is systematically faster due to a somewhat larger stable timestep and fewer operations per timestep. The performance advantage of DIVA over the other solvers increases with resolution. DIVA runs about two times faster than the Hybrid and L1L2-SIA solvers with ∆x = 32 km, and about 15 times faster at ∆x = 4 km. Interestingly, DIVA also runs faster than the simpler SIA solver at high resolution, 30 reaching speeds up to five times faster at ∆x = 4 km. It may be that the numerical implementation of Yelmo's SIA solver can be improved, but in its current form, the lower computational demand per timestep is not sufficient to compensate for a timestep nearly two orders of magnitude smaller than the DIVA timestep. Furthermore, the SIA solver does not resolve the same complexity of ice-flow physics. This is also true of the SSA solver, which does perform marginally better than the DIVA solver, but cannot represent all large-scale ice-flow regimes.

Discussion
We have shown that different stress-balance approximations can be subject to different stability constraints that are not immediately apparent. The analytical stability analysis in Sect. 3 showed that at low resolution (i.e., ∆x > 5 − 10 km), all solvers are subject to a diffusive-like quadratic dependence on grid resolution. Strictly speaking, the stability limit depends on q, which is a function of η = βH 0 /µ and H 0 as well as ∆x. However, resolution is the dominant factor in the tests performed here.

5
At high resolution, two solver families emerge: those whose dynamic timestep limit becomes independent of grid resolution (as shown for the SSA and DIVA solvers) and those whose timestep limit maintains the quadratic dependence. These results were derived for simplified conditions (1D flow with uniform viscosity and constant slope), but the analysis is consistent with Greenland simulations performed without these simplifications. The L1L2-SIA solver is a notable exception, in that a simplified analysis with constant viscosity indicated stability akin to the DIVA solver, albeit offset to a lower maximum timestep limit at high resolution. Moreover, a simplified numerical model based on these assumptions was consistent with the analysis. However, when using comprehensive ice sheet models for either the simplified experiments or the more realistic Greenland experiments, the stability of the L1L2-SIA solver was more consistent with the Hybrid solver. Cornford et al. (2013) noted that using the reconstructed depth-averaged L1L2-SIA velocity in their 5 mass continuity scheme leads to instability as well, and they resort to using only the basal sliding velocity to evolve thickness.
Thus it appears to be challenging to implement the L1L2-SIA solver (which is indeed more complex to implement than the other solvers) in a way that retains the more stable timestepping behavior.
The results of the GrIS experiments confirm the emergence of the two solver families, with the SIA and SSA solvers serving as extreme bounds. The SSA and DIVA solvers show a reduced (less than quadratic) dependence on grid resolution as the 10 grid is refined. The SSA solver maintains a systematically larger timestep than DIVA, most likely because vertical shear stress does not contribute to a reduction in stability. In contrast, the SIA, Hybrid and L1L2-SIA solvers show a strong resolution dependence, with p-values in the relation ∆t ∼ ∆x p ranging from 2.5-2.8. The SIA solver requires a systematically lower timestep than the Hybrid and L1L2-SIA solvers, likely because the driving stress is not dissipated in any way via membrane stresses as in the other two solvers. Analogous simulations for the Antarctic Ice Sheet (not shown), which has large, fast-flowing 15 ice shelves, give similar results as shown here for the GrIS.
Continental-scale simulations using Yelmo at higher resolutions than ∆x = 4 km were not attempted, but the analytical results indicate that a further performance advantage of DIVA could be expected. At such high resolution, factors other than the dynamics, such as basal hydrology and thermodynamics, may play limiting roles in the maximum stable timestep. Also, at a high enough resolution, the advective limit will further restrict the timestep for the DIVA solver.
All solvers were tested in Yelmo with the same numerical treatment of mass conservation. This served to compare the solvers on an equal basis. However, some of the timestep limitations presented here might be alleviated by the use of mass conservation schemes tailored to the solver in question. For example, applying an implicit scheme to the SIA contribution to velocity can 5 improve stability (Bueler and Brown, 2009). Nonetheless, specific schemes tailored to the choice of dynamics solver, including fully implicit approaches, come with their own tradeoffs and limit the flexibility of the model.
Based on our analysis, the key difference in performance between the two families of solvers emerges in ice-flow regimes that are predominantly driven by vertical shear stress. This is consistent with results of the ISMIP-HOM experiments (Pattyn et al., 2008), which compare the ability of different solvers to reproduce expected features of ice dynamics in an idealized 10 setting. Previous work has shown that, when membrane stresses dominate (experiment C of Pattyn et al., 2008), the DIVA, L1L2-SIA, and SSA solvers compare well to the benchmark Stokes solutions for a range of length scales (Pattyn et al., 2008;Goldberg, 2011;Perego et al., 2012). However, when sliding is not permitted and shear stress plays a stronger role (experiment A), the results from all three solvers deviate from those of Stokes solvers as the length scale decreases (Fig. 5). The Hybrid and L1L2-SIA solvers perform quite poorly, as expected, since without sliding, velocities are only represented by the SIA solution.

15
The DIVA solver gives results closer to the Stokes solutions, but with decreasing fidelity at smaller length scales (Goldberg, 2011;Lipscomb et al., 2019). From this, we can understand that the solvers that are less stable numerically as resolution increases (those that reduce to SIA) are also those whose representation of ice-flow physics is least robust.
In general, it should be noted that the simple slab test presented here (Sec. 3) serves as an excellent benchmark for testing the maximum stable timestep of an ice-sheet model formulation. It is computationally cheap, avoids complications related to 20 lateral boundaries (e.g., calving) and is simple to implement. Most importantly, our results show that the relationship of the maximum stable timestep versus resolution determined via this test should be representative of model stability in more realistic experiments.

Conclusions
We have investigated the numerical performance of several commonly used ice-dynamics solvers. We focused on three fast, 25 depth-integrated solvers that permit continental-scale simulations of ice-sheets, namely the Hybrid, L1L2-SIA, and DIVA solvers. These solvers treat both shear and membrane stresses, so they are appropriate for simulating large-scale ice sheets. We included the SIA and SSA solvers as useful boundary cases that treat only shear or membrane stresses, respectively.
As a first step, we derived expressions for maximum stable timesteps for the DIVA, Hybrid, and L1L2-SIA solvers in an idealized 1D configuration. This analysis showed that with coarse resolution, the maximum stable timestep in all three solvers 30 is proportional to the grid resolution squared (i.e., ∆t ∼ ∆x 2 ). For high resolution, the timestep limit in the Hybrid solver maintains this resolution dependence, greatly restricting the timestep. In contrast, this term drops out for the DIVA solver, and the maximum stable timestep (assuming that stability is not limited by the advective timestep) depends on ice thickness and viscosity. The result is that the DIVA solver may use a larger timestep than the Hybrid solver, especially at higher resolutions.
The stability analysis shows that the SSA solver alone has similar stability characteristics to the DIVA solver. The analysis suggests that L1L2-SIA stability should be similar to DIVA, but tests using Yelmo and CISM showed that its numerical stability is closer to that of the Hybrid solver.
We next performed quasi-steady-state simulations of the GrIS using the five solvers for grid resolutions ranging from 4-5 32 km. We found that two families of solvers emerge, largely consistent with the analytical results. The SSA and DIVA solvers are stable for larger timesteps, with a less than quadratic dependence on grid resolution. In contrast, the SIA, Hybrid, and L1L2-SIA solvers show reduced stability at high resolution, with an overall more than quadratic dependence on ∆x. Again, the L1L2-SIA solver falls into the latter family.
Overall, our analysis shows the DIVA solver to be superior to the Hybrid and L1L2-SIA solvers for reasons of greater 10 numerical stability as resolution increases, and preferable to the SIA and SSA solvers because of greater physical fidelity in different parts of the ice sheet. Its representation of the stress balance is consistent with full Stokes solutions over a range of

A1: Velocity derivation in stability analysis
In this appendix we derive Eq. (40) for u, the solution to Eq. (37). We rescale this linear system as where the elements of A are where q is given by Thus, q is non-negative and is determined by two non-dimensional parameters, η = βH/µ and ∆x/H 0 . Here, it is understood that A is circulant. That is, each row is displaced by one element to the right compared to the row above, with a N,1 = a 1,N = −1, 15 The inverse of A, denoted by A −1 , has an approximate analytic expression G (derived in appendix A2) where the approximation improves exponentially the larger the size of the domain L relative to L m , as defined in Eq. (36). The element on row m and column n of G can be written as and L (m, n.N ) = min(|m − n|, |m + N − n|, |m − N + n|).
That is, L is the closest "distance" between m and n accounting for periodicity. Note that q uniquely determines the elements of G. For large q, it can be shown that a 0 ≈ r ≈ 1/q, and for small q, we have a 0 ≈ 1/(2 √ q) and r ≈ 1 − √ q.
For any q > 0, we have 0 < r < 1. Thus, for matrix elements far from the main diagonal (i.e., with sufficiently large M ≡ 5 |m−n|), g m,n becomes negligible compared to a 0 . Suppose we define "negligible" as smaller than a 0 e −p for some p > 0. Then the condition for element g m,n to be negligible is r M < e −p , or equivalently, M ln r < −p. For small q, where r ≈ 1 − √ q, the condition for negligible matrix elements (recalling that ln(1 − x) ≈ −x for small x) becomes M √ q p.
Comparing Eq. (36) and Eq. (79), we find √ q = ∆x/L m , and therefore Eq. (84) can be written as M ∆x pL m . Thus, L m 10 can be interpreted as a length scale over which the elements of the matrix inverse become small, and the condition L m L ensures that entries far from the diagonal approach zero. Moreover, the expression for G shows that velocity, as obtained by inverting Eq. (78), is a localised weighted average of driving stress (i.e., the right hand side of Eq. (78)), with L m as a measure of the "averaging kernel".
We refer to the right hand side of Eq. (78) as R, with elements R n . Recall that thickness and velocity points are staggered, 15 with h un located between h n and h n+1 , so that ∆ n h = h n+1 −h n , and h un = (h n+1 +h n )/2. With δH as in (33), R n is given by R n = −ε h e ikn∆x ρ i g 4µ ∆x(e ik∆x − 1) − ρ i gα 8µH 0 ∆x 2 (e ik∆x + 1) where for compactness we have written the large bracketed expression as κ(k, ∆x). Pre-multiplying Eq. (78) where the sum is taken over all columns of G with non-negligible entries in row m. The series in Eq. (86) can be written as a sum of a constant term (i.e., 1) and two series -one corresponding to n > m and the other to n < m. Without loss of generality, we assume 0 m N and replace L (m, n.N ) by |m − n|. Further, we assume that N is large enough that the terms in the series are negligible for large |m − n|. Thus, the two series (n > m and n < m) can be viewed as infinite geometric sums of 25 the form z(1 + z + z 2 + ...), where z = re ik∆x for the first series and z = re −ik∆x for the second series. These infinite series will converge to z/(1 − z) provided that |z| < 1, which follows from 0 < r < 1 as shown above. This results in the expression given in Sect. 3.1, Eq. (40): where we have defined a function φ that can be simplified and shown to be real: φ(∆x, k, r) = 1 + re ik∆x 1 − re ik∆x + re −ik∆x 1 − re −ik∆x = 1 − r 2 1 + r 2 − 2r cos(k∆x) .

A2: Approximate matrix inverse G
The matrix A in Eq. (78) is a circulant matrix with diagonal terms 2 + q and first off-diagonal terms −1, and zero elsewhere.
Define by a m,n the entry at row m and column n of A, and a inv m,n the entry at row m and column n of A −1 . Since AA −1 = I, the identity matrix, we require that for each row m, the following must hold: −a inv m,m−1 − a inv m,m+1 + (2 + q) a inv m,m = 1, 10 and also that the inner product of row m of A −1 with any column n of A, m = n, is zero, i.e.
We do not derive an analytical expression for a matrix that satisfies Eqs. (89) and (89) for all m and n (and we do not know that one exists). Rather, we derive here an analytically defined matrix G which is close to A −1 in the sense that (GA − I) is small. We choose an ansatz G as follows: where L (m, n.N ) = min(|m − n|, |m + N − n|, |m − N + n|). That is, G is a circulant matrix and L is the "distance" in columns between m and n, but accounting for the circulant property of G. For m = n, Eq. (90) yields a 0 r −1 r L (m,n,N ) − r 2 + (2 + q)r − 1 = 0, yielding the solution 20 r = 1 + q 2 ± 1 + q 2 2 − 1.
We take only the negative branch of the solution, leading to Eq. (81). Although this choice is not immediately apparent, it should become clear below. In particular, Eq. (90) is not satisfied everywhere by G. When L (m, n, N ) = N/2 (where · indicates the floor function), the inner product of G m , the mth row of G, and A n , the nth column of A, is a 0 −r N/2−1 + (2 + q)r N/2 − r N/2−1 = a 0 r N/2 2 + q − 2/r if N is even, and a 0 −r N/2 −1 + (2 + q)r N/2 − r N/2 = a 0 r N/2 1 + q − 1/r if N is odd. We are interested in the limit of high resolution (i.e., small q) and its influence on stability. As discussed in A1, as 5 q goes to zero, we have r ≈ 1 − √ q, and therefore (using log(1 − x) ≈ −x for small x): −log r ∼ √ q = ∆x L m , from which it can be shown that r N/2 is asymptotic to e −p , where p = L/(2L m ). (Note that if the positive branch of Eq. (93) were taken, the non-zero off-diagonal entries of GA would instead grow without bound.) Thus, as long as the numerical domain is sufficiently large compared to L m , the off-diagonal terms of the matrix product GA (and equivalently AG) are 10 bounded by e −p , where the correct choice of a 0 (given below) ensures that the diagonal entries are 1. It remains to find a 0 . Eq.
(89) becomes a 0 2 + q − 2r = 1, yielding Eq. (82). It can be seen numerically (Fig. 6) that the rows of G are a good approximation to those of A −1 . Competing interests. The authors declare that they have no competing interests.
no. PID2019-110714RA-I00). Daniel Goldberg was supported by the Natural Environment Research Council (grants NE/S006796/1 and NE/T001607/1). This material is based upon work supported by the National Center for Atmospheric Research, which is a major facility sponsored by the National Science Foundation under Cooperative Agreement No. 1852977. Figure 6. Difference between coefficients of row N/2 of matrix G and that of A −1 for q = 0.001.