Research article 25 Feb 2022
Research article  25 Feb 2022
A comparison of the stability and performance of depthintegrated icedynamics solvers
 ^{1}Department of Earth Science and Astrophysics, Complutense University of Madrid, Madrid, Spain
 ^{2}Geosciences Institute CSICUCM, Madrid, Spain
 ^{3}Potsdam Institute for Climate Impact Research, Potsdam, Germany
 ^{4}Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder, CO 80305, USA
 ^{5}School of GeoSciences, University of Edinburgh, Edinburgh, UK
 ^{1}Department of Earth Science and Astrophysics, Complutense University of Madrid, Madrid, Spain
 ^{2}Geosciences Institute CSICUCM, Madrid, Spain
 ^{3}Potsdam Institute for Climate Impact Research, Potsdam, Germany
 ^{4}Climate and Global Dynamics Laboratory, National Center for Atmospheric Research, Boulder, CO 80305, USA
 ^{5}School of GeoSciences, University of Edinburgh, Edinburgh, UK
Correspondence: Alexander Robinson (robinson@ucm.es)
Hide author detailsCorrespondence: Alexander Robinson (robinson@ucm.es)
In the last decade, the number of icesheet 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 socalled Hybrid, L1L2SIA and DIVA solvers, as well as the wellknown SIA and SSA solvers as boundary cases. We investigate the numerical stability of these solvers as a function of grid resolution and the state of the ice sheet for an explicit time discretization scheme of the mass conservation step. Under simplified conditions with constant viscosity, the maximum stable time step of the Hybrid solver, like the SIA solver, has a quadratic dependence on grid resolution. In contrast, the DIVA solver has a maximum time step that is independent of resolution as the grid becomes increasingly refined, like the SSA solver. A simple 1D implementation of the L1L2SIA solver indicates that it should behave similarly, but in practice, the complexity of its implementation appears to restrict its stability. In realistic simulations of the Greenland Ice Sheet with a nonlinear rheology, the DIVA and SSA solvers maintain superior numerical stability, while the SIA, Hybrid and L1L2SIA solvers show markedly poorer performance. At a grid resolution of Δx=4 km, the DIVA solver runs approximately 20 times faster than the Hybrid and L1L2SIA solvers as a result of a larger stable time step. Our analysis shows that as resolution increases, the icedynamics solver can act as a bottleneck to model performance. The DIVA solver emerges as a clear outlier in terms of both model performance and its representation of the iceflow physics itself.
 Article
(2565 KB) 
Supplement
(97 KB)  BibTeX
 EndNote
Modeling ice sheets at the continental scale requires compromise. The Greenland and Antarctic ice sheets span hundreds to thousands of kilometers horizontally and interact with the rest of the Earth system on timescales of up to thousands of years. Solving the full Stokes stress balance at high resolution to calculate the ice dynamics over such a large domain for a long timescale is still computationally challenging. Therefore, it is important to reduce the complexity of simulations of ice dynamics.
Several factors reduce the need to simulate the whole system using the full Stokes equations. First, ice sheets typically have low aspect ratios, growing only a few kilometers thick near the domes but extending thousands of kilometers horizontally (i.e., ice sheets are shallow). This condition implies that over most length scales of interest, several terms in the stress balance can be ignored. Furthermore, some boundary conditions, particularly the basal friction at the bed, are poorly known, which means that the added accuracy in simulating velocities with a full Stokes approach may not translate into more robust estimates of past or future icesheet evolution.
Many approximations to the full Stokes stress balance exist, each with tradeoffs. The simplest and historically most widely used shallow approximation is the shallowice approximation (SIA) (Hutter, 1983; Morland, 1984). The SIA is valid for slowflowing grounded ice, frozen to the bedrock, which equates to shear flow induced by gravitational driving stress balanced by basal drag. It is a local solution, in that the velocity diagnosed at a given location is fully determined by the local driving stress. A complementary, but nonlocal, approximation is the shallowshelf approximation (SSA) (Morland, 1987; MacAyeal, 1989), which represents fastflowing ice that is either floating or sliding rapidly over the bed, resulting in plug flow (i.e., with a constant vertical profile of horizontal velocity).
An ad hoc approach to gaining the benefits of computational speed and validity in multiple flow regimes is to combine the SIA and SSA by summing their contributions to obtain the “hybrid” horizontal velocity field (hereafter called the Hybrid solver). This approach was proposed by Bueler (2009), who used a weighting function to transition between the two regimes. Later, Winkelmann et al. (2011) recognized that, in light of the approximations involved and the uncertainty in basal friction, it is more straightforward simply to sum the contributions. The fundamental assumption behind this approach is that the SSA represents a sliding regime, in which the depthaveraged velocity is equal to the basal velocity, while the SIA represents a frozen regime, in which the basal velocity is zero. In the former, it can be expected that SIA velocities will be negligible, and in the latter, SSA velocities will be zero. Only in the narrow transition between the two will the velocity solutions be 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).
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 depthintegrated approximations are especially appealing, since they are comparable in speed per time step to the SSA and Hybrid solvers. The depthintegrated 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 L1L2SIA solver. The stress balance is solved at one 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 shearstress components via the SIA, which facilitates its numerical solution. The L1L2SIA solver is used by the BISICLES model (Cornford et al., 2013).
The socalled depthintegratedviscosity approximation (DIVA) derived by Goldberg (2011) is also part of the L1L2 family of solvers. Here, the stress balance is formulated in terms of the depthaveraged velocity, but like the L1L2SIA solver, the viscosity accounts for both shear and stretching. The shear stress, however, is not approximated, but the longitudinal and lateral stresses are treated as depthindependent, which facilitates a 2D solution of the system of PDEs. The viscosity is vertically averaged, and basal friction is cast in terms of the depthaveraged velocity. This solver has been used by continentalscale models to investigate dynamics and interactions with climate for the Greenland and Antarctic ice sheets (Arthern et al., 2015; Arthern and Williams, 2017; Lipscomb et al., 2019, 2021).
All the approximations described above have been used effectively in a variety of glaciological contexts, but there has been no rigorous, comparative study of their numerical performance. Until recently, the primary concern at the regional/continental scale has been to represent the physics of ice flow, i.e., ensuring that both sliding and shearing regimes are well represented – a motivation that contributed to the development of the Hybrid, L1L2SIA 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 modeling 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 time stepping. 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 and Hybrid solvers, with no analytical solution found for the L1L2SIA solver. We investigate the underlying factors that can affect the maximum stable time step for each solver. The analytical results are also validated in idealized numerical tests. We then compare the performance of all five solvers in terms of stable timestep size and model computational speed in a realistic test case of quasisteadystate simulations of the Greenland Ice Sheet. This is followed by a discussion of the results and conclusions.
We describe the assumptions and equations behind the five solvers considered here, namely the SIA, SSA, DIVA, Hybrid and L1L2SIA solvers. These approximations can be obtained by considering the various terms in the firstorder Blatter–Pattyn (BP) approximation (Blatter, 1995; Pattyn, 2003). To give context to the depthintegrated solvers, we first write the basic equations of the 3D BP stressbalance 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 lefthand side (LHS) describe gradients of longitudinal stress, lateral shear stress, and vertical shear stress, respectively, and the righthand 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
where A is the temperaturedependent rate factor in Glen's flow law (Glen, 1955), n is the Glen exponent (many models set n=3) and ${\dot{\mathit{\epsilon}}}_{\mathrm{e}}$ is the effective strain rate, given in the BP approximation by
Here, u_{x} denotes the partial derivative $\partial u/\partial x$, and similarly for other derivatives.
2.1 SIA solver
The shallowice approximation (SIA) solver is a zeroorder 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 shearstressdriven flow within the ice column. In other words, the membrane stress gradients are ignored, which leads to a stressbalance 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. (1), we obtain the SIA stressbalance equations:
where the effective viscosity μ in Eq. (2) is obtained using an effective strain rate that includes only the shear terms:
Integrating Eq. (4) once, and by assuming a low aspect ratio and a freestress surface boundary condition, we obtain
which can be integrated from the bed b to the height z to provide a 3D solution for horizontal ice flow:
The 2D vertically averaged velocity components $\overline{u}$ and $\overline{v}$ can then be obtained as
Equation (8) relies implicitly on the velocity solution via μ. In practice, an explicit formulation can also be found, which is more typical in model codes (e.g., Greve and Blatter, 2009). Furthermore, we note that the simplified case of uniform μ, Eqs. (7) and (8), can be integrated to give
It is possible to account for sliding in an SIA model, but Eqs. (7)–(9) apply to the regime where ice is frozen to the bedrock.
2.2 SSA solver
The shallowshelf approximation (SSA) solver is complementary to the SIA solver in that the membrane stress gradients are retained, while the verticalshearstress 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 verticalshearstress gradient terms in Eq. (1) and integrating vertically while assuming a vertically uniform velocity, we obtain the 2D SSA stressbalance equations:
where τ_{b,i} indicates the basal stress in the direction i. For the SSA and subsequent approximations, we assume a basalfriction law of the form
where τ_{b} is the basal shear stress, ${\mathit{u}}_{\mathrm{b}}=({u}_{\mathrm{b}},{v}_{\mathrm{b}})$ is the basal velocity and β is a nonnegative friction coefficient that can depend on u_{b}. In the case of the SSA solver, ${\mathit{u}}_{\mathrm{b}}=\overline{\mathit{u}}$. 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:
2.3 DIVA solver
Goldberg (2011) derived a higherorder stress approximation that is similar in accuracy to the firstorder BP stress balance but like other depthintegrated solvers is computationally much cheaper. Since the stressbalance equations use a depthintegrated effective viscosity in place of a vertically varying viscosity, we refer to this scheme as a depthintegratedviscosity approximation, or DIVA.
In the DIVA solver, the horizontal velocity gradients and effective viscosity in the membrane stress terms are replaced by vertical averages $\overline{u}$, $\overline{v}$ and $\overline{\mathit{\mu}}$ (Eq. 8). 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 (14) has the same form as the SSA stress balance, but the dependence of basal stress and viscosity on velocity is more complex.
To solve Eq. (14) for the mean horizontal velocity components, the basal stress terms must be written in terms of $\overline{u}$ and $\overline{v}$. The following discussion refers to the u component; the v component is analogous. Following Arthern et al. (2015), we define some generalized integrals ℱ_{m} for clarity:
Note that when n>1, these integrals depend implicitly on velocity through μ. Goldberg (2011) showed that the vertical profile of u is related to u_{b} by
The surface velocity is thus related to the bed velocity as
Integrating u(z) from the bed to the surface gives the depthaveraged mean velocity $\overline{u}$:
which allows βu_{b} in Eq. (11) to be replaced with ${\mathit{\beta}}_{\mathrm{eff}}\overline{u}$, where
For a frozen bed (u_{b}=0, with nonzero basal stress τ_{b,x}), β will tend to infinity, whereby the above expression can be reduced to
To calculate the 3D viscosity field using the effective strain rate, Eq. (13), the verticalshearstrain terms u_{z} and v_{z} must first be calculated, with u_{z} (and similarly for v_{z}) given by the expression (Lipscomb et al., 2019)
This expression depends on both τ_{b,x} and η, which can be obtained from the previous iteration.
The velocity is found in two steps. First, Eq. (14) is solved iteratively for the mean velocity, using Eqs. (11), (19) and (20) to write the basal stress terms in terms of $\overline{u}$ and $\overline{v}$ and Eqs. (13) and (21) to calculate the effective viscosity. Then Eq. (16) is integrated vertically to find the 3D velocity.
2.4 Hybrid solver
The Hybrid solver, as defined here, follows the approach of Winkelmann et al. (2011). The horizontal velocity is defined as the sum of the depthdependent 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.
2.5 L1L2SIA solver
The L1L2SIA solver as defined here follows from Schoof and Hindmarsh (2010) and Perego et al. (2012). Like the DIVA and Hybrid balances, it is a twostep approach that begins by solving a depthaveraged stress balance. The balance is similar in form to the SSA, Eq. (10):
Aside from the fact that the equations are solved for basal, not depthuniform, velocity, the balance differs from the SSA in the form of the effective viscosity. $\overline{{\mathit{\mu}}^{\prime}}$ is a vertical average of μ^{′}, the L1L2 effective viscosity, which is given by
where
with τ_{b,ij}, $i,j=\mathit{\{}x,y\mathit{\}}$ being the longitudinal stress induced by the basal strainrate tensor. It should be noted that the form of μ^{′} is obtained by assuming that vertical shear stress can be approximated by the SIA solution (see Perego et al., 2012, for details). τ_{b,e} depends implicitly on ${\dot{\mathit{\epsilon}}}_{\mathrm{b},\mathrm{e}}$, the effective horizontal (2D) strain rate induced by the basal strainrate tensor (cf. Eq. 12):
For n=3, taking ${\dot{\mathit{\epsilon}}}_{\mathrm{b},\mathrm{e}}$ as a function of u_{b} and v_{b} from the previous iteration, Eq. (26) is a cubic equation for τ_{b,e} and can be solved analytically. Given $\overline{{\mathit{\mu}}^{\prime}}$ from Eqs. (24) and (26), we can then solve Eq. (23) for u_{b} and v_{b}. Note that the trial viscosity μ^{′} makes the approximation that vertical shear stress is equal to the shallowice shear stress, −ρ_{i}gH∇s. In the second step, threedimensional 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 L1L2SIA is the difference between μ^{′} and μ. The former is defined using the SIA shear stress as a proxy for the true vertical shear stress, while μ is formulated with an updated vertical shear stress based on the solution to the stress balance. This is in contrast to DIVA, where the viscosity is consistent between the two steps.
In this section, we assess the numerical stability of various schemes to solve the coupled stressbalance and continuity equations. Our analysis is closely related to von Neumann stability analysis (Isaacson and Keller, 2012), which has been developed for finite differences, but the nonlocal nature of 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 setting $\overline{v}=\mathrm{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. (2) and the rate factor is constant, set to $A={\left(\mathrm{2}\mathit{\mu}\right)}^{\mathrm{1}}$. This also implies that ${\mathcal{F}}_{\mathrm{2}}=H/\left(\mathrm{3}\mathit{\mu}\right)$ following Eq. (15).

We consider perturbations to an idealized geometry: an xperiodic domain with uniform thickness, H_{0}, and surface slope, $\mathrm{d}s/\mathrm{d}x=\mathit{\alpha}$, where α>0.
This problem is described by Dukowicz (2012), who derived exact solutions for the velocity profile as a function of α and the nondimensional parameter $\mathit{\eta}=\mathit{\beta}H/\mathit{\mu}$.
We specifically analyze the DIVA, Hybrid and L1L2SIA solvers. We begin with DIVA, which leads to the simplest expressions for stability in this framework; results for the Hybrid balance follow. Finally, we attempt to analyze the L1L2SIA balance but find that the complexity of the solver does not lend itself to an analytical result. 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 found previously (Hindmarsh, 2001; Cheng et al., 2017), showing that the maximum stable time step under the simplifications above is proportional to the square of the grid resolution:
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 finitedifference scheme to different stressbalance equations. As in Cheng et al. (2017) our analysis is applied to a linearized form of the coupled stress and continuity equations under small perturbations. Our analytical results provide context for the empirical timestep limits found in more realistic settings in Sect. 4 and give a theoretical basis for those results. Readers interested primarily in the performance of various depthintegrated stress balances in continentalscale icesheet models may wish to skip this section.
3.1 DIVA solver stability
We start with the x component of the DIVA stress balance, Eq. (14). By assuming a constant slope α and spatially uniform effective viscosity μ and by setting $\overline{v}=\mathrm{0}$, this equation can be written as a secondorder ordinary differential equation (ODE) for $\overline{u}$:
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 quantity in parentheses is β_{eff}, the effective basalfriction coefficient. This coefficient includes a nondimensional parameter, $\mathit{\eta}\equiv \mathit{\beta}H/\mathit{\mu}$, 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 ${\mathit{\beta}}_{\mathrm{eff}}\approx \mathrm{3}\mathit{\mu}/H$, and most of the flow is internal shearing. Note that in the context of this analysis, β is considered uniform.
To analyze the stability of a simple numerical scheme for the DIVA equations, we consider a reference state with uniform thickness H_{0} and surface slope α as described above. This reference state has uniform $\overline{u}$ and hence zero longitudinal stress, with depthmean velocity
The reference mean velocity can also be written as
where ${u}_{\mathrm{b},\mathrm{0}}={\mathit{\rho}}_{\mathrm{i}}g{H}_{\mathrm{0}}\mathit{\alpha}/\mathit{\beta}$ is the reference basal velocity. Following mass conservation, the thickness evolution equation can be written as a balance of horizontal flux divergence:
We will consider the evolution of thickness under a sinusoidal perturbation of the initial ice thickness:
Following Cheng et al. (2017), we substitute $H={H}_{\mathrm{0}}+\mathit{\delta}H$ and $\overline{u}={u}_{\mathrm{0}}+\mathit{\delta}\overline{u}$ in Eqs. (32) and (35), where $\mathit{\delta}\overline{u}\left(x\right)$ is a perturbation in the mean velocity. Expanding terms and subtracting the zeroorder balance, we can approximate the DIVA momentum–mass balance to first order in ε_{h}:
Next, we define a simple finitedifference scheme to solve this system of equations. While the theoretical domain has infinite extent, the computational domain must be finite, and we therefore consider a periodic domain of length L. The discretized perturbation thickness, h_{n}, is located at points x=nΔx, for n=0, … ,N (where $\mathrm{\Delta}x=L/N$), and the discretized perturbation velocity, u_{n}, lies at points $\left(n+\frac{\mathrm{1}}{\mathrm{2}}\right)\mathrm{\Delta}x$, for n=0, … ,N−1. We assume that L is much larger than the parameter L_{m}, defined as
As shown in Eq. (A1), L_{m} is equivalent to a membrane stress length scale, over which longitudinal stresses are important (Hindmarsh, 2006).
Equation (37) is discretized 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 firstorder finite difference of a discrete function on the numerical grid, i.e., ${\mathrm{\Delta}}_{n}\mathit{h}=({h}_{n+\mathrm{1}}{h}_{n}$). C is a symmetric, tridiagonal matrix with the entries
where
To proceed, an expression for $\mathit{\delta}\overline{u}$ in terms of δH and its derivatives 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. 40) 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 closedform solution can be derived. The mathematical details can be found in Appendices A1 and A2, but here we state the result:
where
The linearized mass balance equation, Eq. (38), is discretized with an explicit upwind scheme for the advective term and a centered difference for the divergence term:
With ${h}_{\mathrm{m}}^{j}={\mathit{\epsilon}}_{h}{e}^{ikm\mathrm{\Delta}x}$, and using Eq. (43) 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. (50) (i.e., the divergence term) can be written as
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. (51), which is purely imaginary. This is justified as we are primarily concerned with cases where Δx→0, where this term becomes negligible. In Sect. 3.4, however, we test the validity of the assumption over a wide range of Δx values. Using Eq. (45) to write $\mathit{\varphi}(\mathrm{\Delta}x,k,r)$ in full, the remaining term is
This expression is nonpositive, is continuously differentiable in θ, and has extremal points at θ=0 and θ=π. The maximum value occurs when $\mathrm{cos}\left(\mathit{\theta}\right)=\mathrm{cos}\left(\mathit{\pi}\right)=\mathrm{1}$:
We refer to the second term in brackets in Eq. (50) as γ. This term is related to advection and has a maximal value (when ${e}^{ik\mathrm{\Delta}x}=\mathrm{1}$) of
Together, γ and Ξ determine the stability of the timestepping 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, projects onto the mode e^{ikx} for which Eqs. (53) and/or (54) are realized. The scheme is stable when the real part of the bracketed expression in Eq. (50) has an absolute value less than 1 or, equivalently, when $\left\mathit{\gamma}\right+\left\mathrm{\Xi}\right<\mathrm{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 leadingorder behavior in the relationship between resolution and maximal stable time step. When γ≫Ξ, the equation system is essentially advective, and stability is governed by Eq. (54):
When Ξ≫γ, the maximum stable time step follows from Eq. (53):
The dependence of Δt_{dyn} on Δx is more complex than that of Δt_{adv}. It is useful, however, to consider two endmember 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 a limiting factor). The case q≫1 corresponds to high basal friction and/or coarse resolution. In this limit it can be shown (see Eq. A1) that ${a}_{\mathrm{0}}\approx r\approx \mathrm{1}/q$, so that ${a}_{\mathrm{0}}^{\mathrm{1}}(\mathrm{1}+r)/(\mathrm{1}r)\approx q$, and Eq. (56) becomes
If η≫1, as in the case of high basal friction, this simplifies to
In this case, the stability criterion is identical to that of the SIA solver (Eq. 31), with the maximum stable time step, Δ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.
In contrast, for q≪1 (i.e., low basal friction and/or fine resolution), it can be shown that ${a}_{\mathrm{0}}\approx \mathrm{1}/\left(\mathrm{2}\sqrt{q}\right)$ and $r\approx \mathrm{1}\sqrt{q}$. The terms in Eq. (56) containing a_{0} and r reduce to $\mathrm{2}(\mathrm{2}\sqrt{q})\approx \mathrm{4}$, resulting in
Thus, in the limit of small q (or equivalently, small Δx), the timestep limit arising from iceflux divergence depends on μ and H_{0} but not Δx. This is not to say that the limiting time step is independent of resolution – but such dependence arises from the advective term, and thus the maximum stable time step varies linearly, not quadratically, with Δx. For a typical ice thickness of 10^{3} m and effective viscosity of 10^{6} Pa yr, we would have Δt_{dyn}∼1 yr. With Δx=10^{3} m and u_{0}=10^{3} m yr^{−1}, the advective limit from Eq. (55) would also be ∼1 yr. At resolutions finer than ∼10^{3} m, the maximum time step would be determined by the advective limit.
These two regimes of limiting time step correspond to differing behavior of the DIVA equations at different scales. In the largeq limit, Eq. (57), 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 time step on resolution. In the smallq limit, Eq. (59), smallscale thickness oscillations are damped by dynamic thinning due to velocity divergence induced by the oscillations. The velocity gradients are scaleindependent due to averaging of associated driving stresses over the membrane scale, resulting in a scaleindependent damping rate and hence a scaleindependent timestep limit, provided the advective timestep limit is large enough.
3.2 Hybrid solver stability
We consider now the Hybrid stress balance, in which sliding velocity is determined by the SSA momentum balance, Eq. (10). Thickness transport is due to a vertically averaged velocity ${\overline{u}}_{\mathrm{hyb}}$, the sum of sliding velocity and the vertical average of the SIA velocity from Eq. (7), ${\overline{u}}_{\mathrm{sia}}={\mathit{\rho}}_{\mathrm{i}}g{H}^{\mathrm{2}}\mathit{\alpha}/\left(\mathrm{3}\mathit{\mu}\right)$.
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 lowfriction limit of the DIVA balance, Eq. (32), with u replacing $\overline{u}$ since velocity is depthuniform. The SSA mass balance equation is given by Eq. (35), again replacing $\overline{u}$ with u. The RHS of the matrix equation for velocity, Eq. (40), is identical, and the matrix C is modified with B=β. Thus, the dynamically limited time step Δt_{dyn} is given by Eq. (56) but with r and a_{0} modified by the new definition of B. At high resolution, this leads to the same dynamically limited behavior as for DIVA, Eq. (59). At coarse resolution, the dynamically limited time step is similar in form but offset from DIVA due to the modified frictional term.
The Hybrid equations for momentum balance and mass conservation in 1D with constant viscosity μ are
where the mean shallowice velocity ${\overline{u}}_{\mathrm{sia}}$ follows from Eq. (7). We consider the same reference state as in the DIVA analysis. The reference sliding velocity is
and the Hybrid depthaveraged velocity in this state is
Thus, ${\overline{u}}_{\mathrm{hyb},\mathrm{0}}={u}_{\mathrm{0}}$, the DIVA depthaveraged reference velocity from Eq. (34). As before (cf. Eqs. 37 and 38), we consider perturbations in H of ε_{h}e^{ikx} and write the equations to first order in ε_{h}:
where
Since ${u}_{\mathrm{hyb}}^{\mathrm{eff}}>{u}_{\mathrm{0}}$, the Hybrid advective stability constraint is slightly more restrictive than for DIVA.
Consider the third term on the RHS of Eq. (66), which follows from the dependence of u_{SIA} on perturbations in s. If this term is discretized using an explicit secondorder finite difference, it can be expressed as
Note that this term is missing from the analogous equation for the DIVA solver, Eq. (38). With this result and those from the previous section (cf. Eq. 50 for DIVA), it can be shown that the limiting time step is determined by the constraint that
In the above expression, r and a_{0} are as defined as in Eqs. (46) and (47) but with the substitution B=β, giving $q=\mathit{\beta}\mathrm{\Delta}{x}^{\mathrm{2}}/\left(\mathrm{4}\mathit{\mu}{H}_{\mathrm{0}}\right)$. The third term of Eq. (69) is equal to $\mathrm{\Xi}{}_{\mathrm{max}}$ as defined in Eq. (53) from the DIVA analysis. In the limit q≫1 (high basal friction and/or coarse resolution), this term reduces to
Recalling that $\mathit{\eta}\equiv \mathit{\beta}{H}_{\mathrm{0}}/\mathit{\mu}$, the third and fourth terms can be combined to give
Provided the advective timestep limit $\mathrm{\Delta}{t}_{\mathrm{adv}}=\mathrm{\Delta}x/{u}_{\mathrm{hyb}}^{\mathrm{eff}}$ is large, an approximate timestep limit in this regime is given by
which is identical to the largeq limit for DIVA, Eq. (57), and reduces to Eq. (58) when η≫1. When q≪1 (low basal friction and/or fine resolution), the third term in Eq. (69) is bounded independently of Δx (cf. the smallq limit Eq. 59 in the DIVA analysis). For small Δx, the timestep limit for small q is therefore governed by the fourth term:
Since Δt is proportional to Δx^{2}, Eq. (71) becomes very restrictive at high resolution, like the identical equation for the SIA, Eq. (31).
In summary (see Table 1), both DIVA and the Hybrid solver have SIAlike stability (Eqs. 57 and 70) for large q (i.e., high friction and/or coarse resolution). However, in the smallq limit (low friction and/or fine resolution), the Hybrid scheme is governed by the SIA limit, Eq. (71), whereas DIVA is usually governed by a less restrictive advective limit, Eq. (55).
3.3 L1L2SIA solver stability
In this section we attempt to analyze stability limits of the L1L2SIA coupled stress–continuity equations. The complexity of the solver prevents us from deriving analytical expressions for timestep limits as with the other solvers, but we nonetheless comment on possible implications of the linearized expressions.
We first derive the expression for ${\overline{u}}_{\mathrm{L}\mathrm{1}\mathrm{L}\mathrm{2}}$, the depthintegrated L1L2SIA velocity, in the context of our simplified analysis (i.e., 1D flow with uniform viscosity). With these simplifications, Eq. (27) of Perego et al. (2012) for shear stress becomes
Integrating twice in z yields the depthaveraged velocity:
The stress balance which is solved for u_{b} is identical to that of the Hybrid balance, and the reference velocity, ${\overline{u}}_{\mathrm{L}\mathrm{1}\mathrm{L}\mathrm{2},\mathrm{0}}$, is identical to that of the Hybrid and DIVA balances. Linearizing the L1L2SIA scheme around this state leads to a perturbed momentum balance identical to that of the Hybrid balance, and the perturbed mass balance is
The L1L2SIA perturbed mass balance above has a similar form to those found for the DIVA and Hybrid cases but contains higherorder derivatives of δu_{b}. Treating these terms analytically as in Sects. 3.1 and 3.2 is beyond the scope of our study, and so we cannot give a rigorous stability bound in terms of resolution and governing parameters.
It is worth pointing out, however, that the second term in the parentheses in Eq. (74) gives rise to a secondderivative term (a diffusive term) in δH, which leads to the quadratic stability limit for small Δx in the Hybrid balance. In contrast to the Hybrid balance, this term is balanced by a term involving the second derivative of δu_{b}. Rearranging the linearized momentum balance, Eq. (65), yields
in which the LHS is equivalent to the terms in parentheses in Eq. (74). In other words, there is a degree of cancellation in Eq. (74), and this cancellation might lead to stabilization. While this argument does not provide a bound for stability, it suggests that the stability might be tied to the degree of cancellation of higherorder derivative terms in discretized form. Note also that, with a nonlinear rheology, the viscosity in Eq. (74) differs from that of the momentum balance (cf. Sect. 2.5), so that such cancellation will not occur for n>1, regardless of the numerical scheme.
3.4 Numerical validation
To confirm the validity of the stability bounds found in the above sections, we run numerical simulations with 1D DIVA, Hybrid and L1L2SIA models derived using the assumptions above. Although we do not have an analytical expression to compare to for the L1L2SIA solver, its numerical behavior can be compared to the other solvers. For the DIVA solver we use a firstorder upwind finitevolume scheme for Eq. (35), and we use Eq. (32) for depthaveraged velocity. For the Hybrid solver, we use Eqs. (61) and (62), and for the L1L2SIA solver, we use Eqs. (61), (73) and (35). The 1D code therefore implements the full coupled mass and momentum equations for DIVA, Hybrid and L1L2SIA, rather than their linear perturbed forms. 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 icesheet models, Yelmo (Robinson et al., 2020) and CISM (Lipscomb et al., 2019). The three solvers have been implemented in both models, and all test cases are simulated.
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 $\mathit{\alpha}=\mathrm{d}s/\mathrm{d}x={\mathrm{10}}^{\mathrm{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
where 𝒳_{n} represents independent Gaussian random variables with zero mean and standard deviation of 0.1 m. We then run the model forward for 100 time steps. 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 realized).
For a given set of physics parameters ($\mathit{\mu},{H}_{\mathrm{0}},\mathit{\beta}$), Δx is varied over a range from 10 m to 40 km (12 values are chosen, equidistant in log space – see the MATLAB code in the Supplement to reproduce the exact values). At each resolution we run multiple tests with different values of Δt to determine the maximum stable time step. For each run, we calculate the ratio $\mathit{\sigma}/{\mathit{\sigma}}_{\mathrm{0}}$, where σ is the standard deviation of h at the final time step, 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 $\mathit{\sigma}/{\mathit{\sigma}}_{\mathrm{0}}\le \mathrm{1}$ to indicate stability and $\mathit{\sigma}/{\mathit{\sigma}}_{\mathrm{0}}>\mathrm{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 iceflow regimes. The first parameter set ($\mathit{\mu}=\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}$ Pa yr, H_{0}=1000 m, β=1000 Pa yr m^{−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 ($\mathit{\mu}=\mathrm{4}\times {\mathrm{10}}^{\mathrm{5}}$ Pa yr, H_{0}=500 m, β=30 Pa yr m^{−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 this case, the maximum stable time step Δt_{dyn} is generally larger (as a function of Δx) than in the thick, shearing case. Also, u_{0} is larger (150 m yr^{−1}), 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.
For the DIVA solver, the maximum stable time steps determined by the 1D model, Yelmo and CISM are in excellent agreement with the analytical estimate of Δt_{dyn} as given by Eq. (53) for both tests. In the shearing case, the advective timestep limit is large enough that it is not relevant, so the maximum stable time step 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 time step follows a similar dependence until Δx<10^{2} m, at which point the advective limit becomes the limiting factor.
For the Hybrid solver, we again find excellent agreement with the analytical estimate of Δt_{dyn}, as given by Eq. (69), for the simplified 1D model, Yelmo and CISM. In both tests, the maximum stable time step depends quadratically 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 time step drops quickly as Δx decreases, such that it is already as small as $\mathrm{\Delta}t\sim {\mathrm{10}}^{\mathrm{2}}$ yr when $\mathrm{\Delta}x\sim {L}_{\mathrm{m}}\sim \mathrm{1000}$ m. For the Hybrid solver, the very small timestep limits ensure that the advective timestep limit never becomes relevant in the two tests.
In the case of the L1L2SIA solver, we cannot compare to an analytical solution. It is notable, however, that each of the three numerical solvers tested (1D solver, Yelmo and CISM) gives different estimates of the maximum stable time step. The 1D solver is the most stable. In the shearing test at high resolution, the time step approaches a constant value like the DIVA solver, although offset to a lower value. For extremely high resolution (Δx<10 m), the 1D model becomes unstable in the shearing case – 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 show a quadratic dependence on grid resolution akin to that of the Hybrid solver for this test. At higher resolutions, their estimates differ, with CISM being somewhat more stable. In the sliding test, the 1D model again resembles the behavior of the DIVA solver, while Yelmo and CISM show stability more like that of the Hybrid solver.
It is not clear why both Yelmo and CISM have stability behavior for L1L2SIA similar to that for the Hybrid balance (yet slightly different from each other), whereas the 1D model has resolutionindependent stability (at least for high shearing). One possible explanation is suggested in Sect. 3.3: that the stability of the scheme is tied to the degree to which the terms in parentheses in Eq. (74), the linearized L1L2SIA balance, cancel each other (or, equivalently, the degree to which the terms in parentheses in Eq. 73, the depthaverage velocity, cancel). In the 1D model, these terms are discretized similarly to the way they are discretized in the momentum balance. That is, the second derivative of u_{b} in Eq. (73) is discretized at a velocity point m as
where ${u}_{\mathrm{b},m,m\pm \mathrm{1}}$ represent the numerical solution of sliding velocity at points $m,m\pm \mathrm{1}$. To study the impacts of discretization, we consider a 1D model for L1L2SIA which is identical to that described above but with Eq. (77) replaced by
Equation (78) is consistent with a second derivative but is distinct from that used in the momentum solve. The stability results for this model are very similar to those for Yelmo and CISM (not shown), which indicates that the L1L2SIA solver is quite sensitive to discretization choices.
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 L1L2SIA solvers. As shown in the analytical derivation, the DIVA solver has the same stability limits as the SSA solver at high resolution; thus, it can be expected to be more stable than the Hybrid solver. In contrast, previous work has shown that SIA stability depends quadratically on Δx at higher resolutions (Cheng et al., 2017), consistent with the results found for the Hybrid solver that incorporates the SIA solution. It appears that the L1L2SIA solver stability limits may be highly dependent on discretization choices. While a simple 1D model is able to reproduce DIVAlike stability limits at high resolution, the two comprehensive models indicate that the L1L2SIA solver may be less stable in practice.
The above analysis sheds light on the mathematical basis for stability differences between approximations in idealized 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 quasisteadystate simulations for the Greenland Ice Sheet (GrIS) at different resolutions. We use the icesheet model Yelmo (Robinson et al., 2020), which supports all five solvers (SIA, SSA, Hybrid, L1L2SIA and DIVA). By using Yelmo for all experiments, we ensure a clean comparison among solvers. We run experiments at four resolutions for each solver: 32, 16, 8 and 4 km. The simulations are designed to represent the complexity in typical simulations of real ice sheets. Each simulation starts from the observed GrIS geometry and uses presentday boundary forcing as in, e.g., the initMIP experiments (Goelzer et al., 2018; Seroussi et al., 2019), with a nonlinear rheology following Glen's flow law with n=3. The basal hydrology and thermodynamics are interactive, with a linear basalfriction law that depends on the basal water layer thickness. The exact configuration is not critical to the analysis, as similar behavior is observed under a variety of conditions and domains.
To save computation time, the simulations are not run fully to steady state. Rather, we run the model for a total of 2 kyr. For the first 1 kyr, we briefly spin up the thermodynamics and basal hydrology and let the model reach a selfconsistent state with the fixed, presentday topography. We then evaluate the timestep stability for a fully prognostic simulation during an additional 1 kyr.
Yelmo uses an adaptive timestepping scheme that determines the optimal time step 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 time step (in Yelmo, this is found by comparing the predicted to corrected ice thickness in a predictor–corrector method), 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 time step is reduced when η_{t}>ε_{t} and increased when η_{t}<ε_{t}, leading to a value that is as large as possible while maintaining stability. Over time, the time step oscillates around this optimal value. In the following analysis, the mean adaptive time step 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 Sect. 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, L1L2SIA and DIVA), the 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 fastsliding 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 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 presentday 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, L1L2SIA 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 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 time step 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 Sect. 3, since the simulations here correspond more closely to the lowresolution regime highlighted in Table 1. The three solvers that rely on the SIA in some form (SIA, Hybrid and L1L2SIA) have similarly high values of p, ≈2.6–2.8, and thus a large decrease in time step with finer resolution. In contrast, the SSA and DIVA solvers have a much weaker dependence on grid resolution (p≈1.7–1.8), resulting in larger stable time steps. At the coarsest resolution, Δx=32 km, all solvers allow a time step of a similar order of magnitude, between 1–5 yrs, although the SIA family requires time steps more than 2 times smaller than the SSA family. As the resolution increases to Δx=4 km, the SIA family requires time steps of less than Δt=0.01 yr, while the SSA and DIVA solvers remain stable for time steps at least an order of magnitude larger.
We also tested solver performance for a constant time step to ensure that the adaptive timestepping scheme provides good estimates of the maximum stable time step. We ran simulations for all solvers for the GrIS domain at 16 km resolution with fixed time steps of Δt=0.2, 0.5, 1.0, 2.0 and 5.0 yr. From the simulations with adaptive time stepping, we expect that all five solvers would be stable at Δt=0.5 yr but that the SIA, Hybrid and L1L2SIA solvers would become unstable for larger time steps. Additionally, we expect that the DIVA solver would become unstable for time steps larger than Δt=1.0 yr and that the SSA solver would be unstable for the largest time step tested here, Δt=5.0 yr. Figure 4 confirms that this is the case. This supports the notion that the adaptive timestepping scheme successfully diagnoses the maximum stable time step allowed in each simulation.
The difference in maximum stable time step 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 L1L2SIA solvers, while the SIA solver is on par with the fastest solvers because it is computationally much less intensive. As the resolution increases, the time step dependency on resolution becomes the limiting factor, and the solvers again separate into the faster SSA family and the slower SIA family.
With the exception of the SIA solver, all solvers are comparable in terms of computational time per time step. Ordered from fastest to slowest per time step, the solver ranking is SIA, SSA, DIVA, Hybrid and L1L2SIA. The SIA solver does not require any matrix solutions, so it is much less expensive. The SSA solver requires a computationally intensive matrix solution but does not require additional calculations for the 3D horizontal velocity field due to the assumption of no vertical shear. The Hybrid, L1L2SIA and DIVA solvers use the same 2D matrixsolution method as the SSA solver, and all require some vertical integration as well. The L1L2SIA solver requires more intermediate calculations than the Hybrid or DIVA solvers, but for a Glen's flow law exponent of n=3, as is the case here, an additional iterative step to determine the effective viscosity for the L1L2SIA solver can be avoided via an analytical solution. Further analysis (not shown) indicates that all solvers require a similar number of Picard iterations to arrive at the converged matrix solution. This shows that the model time step is the primary determinant of overall model speeds.
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 time step and fewer operations per time step. The performance advantage of DIVA over the remaining solvers increases with resolution. DIVA runs about 2–4 times faster than the Hybrid and L1L2SIA solvers with Δx=32 km and about 20 times faster at Δx=4 km. Interestingly, DIVA also runs faster than the simpler SIA solver at high resolution, reaching speeds up to 5 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 time step does not offset a time step nearly 2 orders of magnitude smaller than the DIVA time step. Furthermore, the SIA solver does not resolve the same complexity of iceflow physics. This is also true of the SSA solver, which does perform marginally better than the DIVA solver but cannot represent all largescale iceflow regimes.
We have shown that different stressbalance 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 diffusivelike quadratic dependence on grid resolution. Strictly speaking, the stability limit depends on q, which is a function of $\mathit{\eta}=\mathit{\beta}{H}_{\mathrm{0}}/\mathit{\mu}$ and H_{0} as well as Δx. However, resolution is the dominant factor in the tests performed here. 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 L1L2SIA solver is a notable exception, in that a simplified 1D model with constant viscosity indicates stability akin to the DIVA solver, albeit offset to a lower maximum timestep limit at high resolution. However, when using comprehensive icesheet models for either the simplified experiments or the more realistic Greenland experiments, the stability of the L1L2SIA solver is more consistent with the Hybrid solver. Cornford et al. (2013) noted that using the reconstructed depthaveraged L1L2SIA velocity in their 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 L1L2SIA solver (which is indeed more complex 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 grid is refined. The SSA solver maintains a systematically larger time step than DIVA, most likely because vertical shear stress does not contribute to a reduction in stability. In contrast, the SIA, Hybrid and L1L2SIA solvers show a stronger resolution dependence, with p values in the relation Δt∼Δx^{p} ranging from 2.6–2.8. The SIA solver requires a systematically lower time step than the Hybrid and L1L2SIA 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, fastflowing ice shelves, give similar results as shown here for the GrIS.
The maximum stable time step depends on resolution for all solvers in these realistic tests, in contrast to the 1D and analytical results. At least three factors may contribute to this result. First, the Glen flow law in the realistic tests is nonlinear with n=3. Second, viscosity is no longer constant and therefore makes a highly nonlinear contribution to the solution. Third, icefront boundary conditions must be considered, which were avoided in the idealized framework. For this reason, the analytical results provide a general guide of the relative stability of the solvers but not a realistic estimate of the stable time step under all conditions.
Continentalscale 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 time step. Also, at a high enough resolution, the advective limit will further restrict the time step 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 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 iceflow regimes that are predominantly driven by vertical shear stress. This is consistent with results of the ISMIPHOM experiments (Pattyn et al., 2008), which compare the ability of different solvers to reproduce expected features of ice dynamics in an idealized setting. Previous work has shown that, when membrane stresses dominate (experiment C of Pattyn et al., 2008), the DIVA, L1L2SIA 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 L1L2SIA solvers perform quite poorly, as expected, since without sliding, velocities are only represented by the SIA solution. 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 iceflow physics is least robust.
In general, it should be noted that the simple slab test presented here (Sect. 3) serves as an excellent benchmark for testing the maximum stable time step of an icesheet model formulation. It is computationally cheap, avoids complications related to lateral boundaries (e.g., calving) and is simple to implement. Most importantly, our results show that the relationship of the maximum stable time step versus resolution determined via this test should be representative of model stability in more realistic experiments.
We have investigated the numerical performance of several commonly used icedynamics solvers. We focused on three fast, depthintegrated solvers that permit continentalscale simulations of ice sheets, namely the Hybrid, L1L2SIA and DIVA solvers. These solvers treat both shear and membrane stresses, so they are appropriate for simulating largescale 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 time steps for the DIVA, Hybrid and L1L2SIA solvers in an idealized 1D configuration. This analysis showed that with coarse resolution, the maximum stable time step in all three solvers 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 time step. In contrast, this term drops out for the DIVA solver, whose maximum stable time step (assuming that stability is not limited by the advective time step) depends only on ice thickness and viscosity. The result is that the DIVA solver may use a larger time step than the Hybrid solver, especially at higher resolutions. The stability analysis shows that only the DIVA solver has similar stability characteristics to the SSA solver. The analysis was inconclusive for the L1L2SIA solver. A 1D model suggests that its stability should be similar to DIVA, but tests using Yelmo and CISM show that its stability is closer to that of the Hybrid solver.
We next performed quasisteadystate simulations of the GrIS using the five solvers for grid resolutions ranging from 4–32 km. Two families of solvers emerge, largely consistent with the analytical results. The SSA and DIVA solvers are stable for larger time steps, with a lessthanquadratic dependence on grid resolution. In contrast, the SIA, Hybrid and L1L2SIA solvers show reduced stability at high resolution, with an overall morethanquadratic dependence on Δx.
Overall, our analysis shows the DIVA solver to be superior to the Hybrid and L1L2SIA solvers for reasons of greater 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 length scales. As continentalscale simulations are performed at higher resolutions, the Hybrid and L1L2SIA solvers may become bottlenecks for model performance, while the DIVA solver remains computationally efficient.
In this appendix we derive Eq. (43) for u, the solution to Eq. (40). We rescale this linear system as
where the elements of A are
where q is given by
Thus, q is nonnegative and is determined by two nondimensional parameters, $\mathit{\eta}=\mathit{\beta}H/\mathit{\mu}$ and $\mathrm{\Delta}x/{H}_{\mathrm{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,\mathrm{1}}={a}_{\mathrm{1},N}=\mathrm{1}$,
The inverse of A, denoted by A^{−1}, has an approximate analytic expression G (derived in Eq. A2), where the approximation improves exponentially the larger the size of the domain L relative to L_{m} is, as defined in Eq. (39). The element on row m and column n of G can be written as
where
and
That is, ℒ 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}_{\mathrm{0}}\approx r\approx \mathrm{1}/q$, and for small q, we have ${a}_{\mathrm{0}}\approx \mathrm{1}/\left(\mathrm{2}\sqrt{q}\right)$ and $r\approx \mathrm{1}\sqrt{q}$. For any q>0, we have $\mathrm{0}<r<\mathrm{1}$. Thus, for matrix elements far from the main diagonal (i.e., with sufficiently large $M\equiv mn$), 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\mathrm{ln}r<p$. For small q, where $r\approx \mathrm{1}\sqrt{q}$, the condition for negligible matrix elements (recalling that $\mathrm{ln}(\mathrm{1}x)\approx x$ for small x) becomes
Comparing Eqs. (39) and (A3), we find $\sqrt{q}=\mathrm{\Delta}x/{L}_{\mathrm{m}}$, and therefore Eq. (A8) can be written as MΔx≳pL_{m}. Thus, L_{m} 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. (A1), is a localized weighted average of driving stress (i.e., the righthand side of Eq. A1), with L_{m} as a measure of the averaging kernel.
We refer to the righthand side of Eq. (A1) as R, with elements R_{n}. Recall that thickness and velocity points are staggered, with h_{un} located between h_{n} and h_{n+1}, so that ${\mathrm{\Delta}}_{n}h={h}_{n+\mathrm{1}}{h}_{n}$ and ${h}_{un}=({h}_{n+\mathrm{1}}+{h}_{n})/\mathrm{2}$. With δH as in Eq. (36), R_{n} is given by
where for compactness we have written the large bracketed expression as κ(k,Δx). Premultiplying Eq. (A1) by G, the velocity u_{m} can be computed as
where the sum is taken over all columns of G with nonnegligible entries in row m. The series in Eq. (A10) 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 $\mathrm{0}\ll m\ll N$ and replace $\mathcal{L}(m,n.N)$ by $mn$. Further, we assume that N is large enough that the terms in the series are negligible for large $mn$. Thus, the two series (n>m and n<m) can be viewed as infinite geometric sums of the form $z(\mathrm{1}+z+{z}^{\mathrm{2}}+\mathrm{\dots})$, where z=re^{ikΔx} for the first series and $z=r{e}^{ik\mathrm{\Delta}x}$ for the second series. These infinite series will converge to $z/(\mathrm{1}z)$ provided that $\leftz\right<\mathrm{1}$, which follows from $\mathrm{0}<r<\mathrm{1}$ as shown above. This results in the expression given in Sect. 3.1, Eq. (43):
where we have defined a function ϕ that can be simplified and shown to be real:
The matrix A in Eq. (A1) is a circulant matrix with diagonal terms 2+q and first offdiagonal terms −1, and zero elsewhere. Define by a_{m,n} the entry at row m and column n of A and by ${a}_{m,n}^{inv}$ the entry at row m and column n of A^{−1}. Since ${\mathbf{AA}}^{\mathrm{1}}=\mathbf{I}$, the identity matrix, we require that for each row m, the following must hold:
We also require 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. (B1) and (B2) 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 $\mathcal{L}(m,n.N)=\mathrm{min}\left(\rightmn,m+Nn,mN+n\left\right)$. That is, G is a circulant matrix and ℒ is the distance in columns between m and n but accounting for the circulant property of G. For m≠n, Eq. (B2) yields
yielding the solution
We take only the negative branch of the solution, leading to Eq. (A5). Although this choice is not immediately apparent, it should become clear below. In particular, Eq. (B2) is not satisfied everywhere by G. When $\mathcal{L}(m,n,N)=\lfloor N/\mathrm{2}\rfloor $ (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
if N is even, and it is
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 q goes to zero, we have $r\approx \mathrm{1}\sqrt{q}$, and therefore (using $\mathrm{log}(\mathrm{1}x)\approx x$ for small x):
from which it can be shown that ${r}^{N/\mathrm{2}}$ is asymptotic to e^{−p}, where $p=L/\left(\mathrm{2}{L}_{\mathrm{m}}\right)$. (Note that if the positive branch of Eq. B5 were taken, the nonzero offdiagonal entries of GA would instead grow without bound.) Thus, as long as the numerical domain is sufficiently large compared to L_{m}, the offdiagonal terms of the matrix product GA (and equivalently AG) are 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}. Equation (B1) becomes
yielding Eq. (A6). It can be seen numerically (Fig. B1) that the rows of G are a good approximation to those of A^{−1}.
Yelmo is maintained as a Git repository hosted at https://github.com/palmaice/yelmo (last access: 2 November 2021) under the license GPL3.0. Model documentation can be found at https://palmaice.github.io/yelmodocs/ (last access: 19 December 2021). The exact version of the model used to produce the results used in this paper is archived on Zenodo at https://doi.org/10.5281/zenodo.5791864 (Robinson, 2021) and has been tagged in the repository as solverstabilityv1.0. CISM is an opensource code developed on the Earth System Community Model Portal (ESCOMP) Git repository at https://github.com/ESCOMP/CISM (last access: 17 November 2021). The version used for these runs has been archived on Zenodo at https://doi.org/10.5281/zenodo.5889016 (Lipscomb et al., 2022) and tagged in the repository as cism_slab_tests.
The supplement related to this article is available online at: https://doi.org/10.5194/tc166892022supplement.
AR conceived of the study. DG derived and performed the analytical stability analysis with input from WHL. AR performed the experiments with Yelmo, and WHL ran the experiments with CISM. All authors contributed to the analysis of the results and writing of the manuscript.
At least one of the (co)authors is a member of the editorial board of The Cryosphere. The peerreview process was guided by an independent editor, and the authors also have no other competing interests to declare.
Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Alexander Robinson was supported by the Ramón y Cajal Programme of the Spanish Ministry for Science, Innovation and Universities (grant no. RYC201620587), as well as the Spanish Ministry of Science and Innovation project ICEAGE (grant no. PID2019110714RA100). Daniel Goldberg was supported by the Natural Environment Research Council (grant nos. 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.
This research has been supported by the Ministerio de Ciencia, Innovación y Universidades (grant no. RYC201620587), the Ministerio de Ciencia e Innovación (grant no. PID2019110714RA100), and the Natural Environment Research Council (grant nos. NE/S006796/1 and NE/T001607/1).
This paper was edited by Olivier Gagliardini and reviewed by Mauro Perego and one anonymous referee.
Arthern, R. J. and Williams, C. R.: The sensitivity of West Antarctica to the submarine melting feedback, Geophys. Res. Lett., 44, 2352–2359, https://doi.org/10.1002/2017GL072514, 2017. a
Arthern, R. J., Hindmarsh, R. C. A., and Williams, C. R.: Flow speed within the Antarctic ice sheet and its controls inferred from satellite observations, J. Geophys. Res.Earth, 120, 1171–1188, https://doi.org/10.1002/2014JF003239, 2015. a, b
Blatter, H.: Velocity and stress fields in grounded glaciers – a simple algorithm for including deviatoric stress gradients, J. Glaciol., 41, 333–344, 1995. a
Bueler, E.: Lectures at Karthaus: Numerical modelling of ice sheets and ice shelves, https://glaciers.gi.alaska.edu/sites/default/files/Notes_icesheetmod_Bueler2014.pdf (last access: 23 February 2022), 2009. a
Bueler, E. and Brown, J.: Shallow shelf approximation as a “sliding law” in a thermodynamically coupled ice sheet model, J. Geophys. Res., 114, F03008, https://doi.org/10.1029/2008JF001179, 2009. a
Cheng, G., Lötstedt, P., and von Sydow, L.: Accurate and stable time stepping in ice sheet modeling, J. Comput. Phys., 329, 29–47, https://doi.org/10.1016/j.jcp.2016.10.060, 2017. a, b, c, d, e, f, g
Cornford, S. L., Martin, D. F., Graves, D. T., Ranken, D. R., Le Brocq, A. M., Gladstone, R. M., Payne, A. J., Ng, E. G., and Lipscomb, W. H.: Adaptive mesh, finite volume modeling of marine ice sheets, J. Comput. Phys., 232, 529–549, 2013. a, b
Dukowicz, J. K.: Reformulating the fullStokes ice sheet model for a more efficient computational solution, The Cryosphere, 6, 21–34, https://doi.org/10.5194/tc6212012, 2012. a
Glen, J. W.: The creep of polycrystalline ice, Proc. R. Soc. London A, 228, 519–538, 1955. a
Goelzer, H., Nowicki, S., Edwards, T., Beckley, M., AbeOuchi, A., Aschwanden, A., Calov, R., Gagliardini, O., GilletChaulet, F., Golledge, N. R., Gregory, J., Greve, R., Humbert, A., Huybrechts, P., Kennedy, J. H., Larour, E., Lipscomb, W. H., Le clec'h, S., Lee, V., Morlighem, M., Pattyn, F., Payne, A. J., Rodehacke, C., Rückamp, M., Saito, F., Schlegel, N., Seroussi, H., Shepherd, A., Sun, S., van de Wal, R., and Ziemen, F. A.: Design and results of the ice sheet model initialisation experiments initMIPGreenland: an ISMIP6 intercomparison, The Cryosphere, 12, 1433–1460, https://doi.org/10.5194/tc1214332018, 2018. a, b, c
Goelzer, H., Nowicki, S., Payne, A., Larour, E., Seroussi, H., Lipscomb, W. H., Gregory, J., AbeOuchi, A., Shepherd, A., Simon, E., Agosta, C., Alexander, P., Aschwanden, A., Barthel, A., Calov, R., Chambers, C., Choi, Y., Cuzzone, J., Dumas, C., Edwards, T., Felikson, D., Fettweis, X., Golledge, N. R., Greve, R., Humbert, A., Huybrechts, P., Le clec'h, S., Lee, V., Leguy, G., Little, C., Lowry, D. P., Morlighem, M., Nias, I., Quiquet, A., Rückamp, M., Schlegel, N.J., Slater, D. A., Smith, R. S., Straneo, F., Tarasov, L., van de Wal, R., and van den Broeke, M.: The future sealevel contribution of the Greenland ice sheet: a multimodel ensemble study of ISMIP6, The Cryosphere, 14, 3071–3096, https://doi.org/10.5194/tc1430712020, 2020. a
Goldberg, D. N.: A variationally derived, depthintegrated approximation to a higherorder glaciological flow model, J. Glaciol., 57, 157–170, https://doi.org/10.3189/002214311795306763, 2011. a, b, c, d, e, f
Greve, R. and Blatter, H.: Dynamics of Ice Sheets and Glaciers, 1st edn., SpringerVerlag, Berlin, https://doi.org/10.1007/9783642034152, 2009. a, b
Hindmarsh, R.: A numerical comparison of approximations to the Stokes equations used in ice sheet and glacier modeling, J. Geophys. Res., 109, F01012, https://doi.org/10.1029/2003JF000065, 2004. a, b
Hindmarsh, R.: The role of membranelike stresses in determining the stability and sensitivity of the Antarctic ice sheets: back pressure and grounding line motion, Philos. T. R. Soc. A, 364, 1733–1767, https://doi.org/10.1098/rsta.2006.1797, 2006. a, b
Hindmarsh, R. C.: Notes on basic glaciological computational methods and algorithms, in: Continuum mechanics and applications in geophysics and the environment, edited by Straughan, B., Greve, R., Ehrentraut, H., and Wang, Y., 222–249, Springer, Berling, http://nora.nerc.ac.uk/id/eprint/19612/ (last access: 23 February 2022), 2001. a, b
Hutter, K.: Theoretical Glaciology, Mathematical Approaches to Geophysics, 1st edn., edited by: Rikitaki, T. and Hazenwinkel, M., D. Reidel Publishing Company, Dordrecht, Boston, Lancaster, https://doi.org/10.1007/9789401511674, 1983. a
Isaacson, E. and Keller, H. B.: Analysis of numerical methods, 1st edn., Courier Corporation, ISBN 9780486680293, 2012. a
Lipscomb, W. H., Price, S. F., Hoffman, M. J., Leguy, G. R., Bennett, A. R., Bradley, S. L., Evans, K. J., Fyke, J. G., Kennedy, J. H., Perego, M., Ranken, D. M., Sacks, W. J., Salinger, A. G., Vargo, L. J., and Worley, P. H.: Description and evaluation of the Community Ice Sheet Model (CISM) v2.1, Geosci. Model Dev., 12, 387–424, https://doi.org/10.5194/gmd123872019, 2019. a, b, c, d
Lipscomb, W. H., Leguy, G. R., Jourdain, N. C., AsayDavis, X., Seroussi, H., and Nowicki, S.: ISMIP6based projections of oceanforced Antarctic Ice Sheet evolution using the Community Ice Sheet Model, The Cryosphere, 15, 633–661, https://doi.org/10.5194/tc156332021, 2021. a
Lipscomb, W., Sacks, B., Price, S., Hoffman, M., Martin, D. F., Ranken, D., Hagdorn, M., Salinger, A., Evans, K., Norman, M., Shannon, S., ThayerCalder, K., Leguy, G., Tezaur, I. K., Kennedy, J. H., Goelzer, H., Kluzek, E., Rutt, I. C., and Iesulauro Barker, E.: ESCOMP/CISM: CISM Slab Tests (cism_slab_tests), Zenodo [code], https://doi.org/10.5281/zenodo.5889016, 2022. a
MacAyeal, D. R.: Largescale ice flow over a viscous basal sediment – Theory and application to Ice Stream B, Antarctica, J. Geophys. Res., 94, 4071–4087, 1989. a
Morland, L.: Unconfined iceshelf flow, in: Dynamics of the West Antarctic Ice Sheet, edited by: Van der Veen, C. J. and Oerlemans, J., Glac. Quat. G., 4, 99–116, Springer, https://doi.org/10.1007/9789400937451_6, 1987. a
Morland, L. W.: Thermomechanical balances of ice sheet flows, Geophys. Astro. Fluid, 29, 237–266, https://doi.org/10.1080/03091928408248191, 1984. a
Pattyn, F.: A new threedimensional higherorder thermomechanical ice sheet model: Basic sensitivity, ice stream development, and ice flow across subglacial lakes, J. Geophys. Res., 108, 2382, https://doi.org/10.1029/2002JB002329, 2003. a
Pattyn, F.: Sealevel response to melting of Antarctic ice shelves on multicentennial timescales with the fast Elementary Thermomechanical Ice Sheet model (f.ETISh v1.0), The Cryosphere, 11, 1851–1878, https://doi.org/10.5194/tc1118512017, 2017. a
Pattyn, F., Perichon, L., Aschwanden, A., Breuer, B., de Smedt, B., Gagliardini, O., Gudmundsson, G. H., Hindmarsh, R. C. A., Hubbard, A., Johnson, J. V., Kleiner, T., Konovalov, Y., Martin, C., Payne, A. J., Pollard, D., Price, S., Rückamp, M., Saito, F., Souček, O., Sugiyama, S., and Zwinger, T.: Benchmark experiments for higherorder and fullStokes ice sheet models (ISMIP–HOM), The Cryosphere, 2, 95–108, https://doi.org/10.5194/tc2952008, 2008. a, b, c, d
Perego, M., Gunzburger, M., and Burkardt, J.: Parallel finiteelement implementation for higherorder ice sheet models, J. Glaciol., 58, 76–88, https://doi.org/10.3189/2012JoG11J063, 2012. a, b, c, d, e
Pollard, D. and DeConto, R. M.: Description of a hybrid ice sheetshelf model, and application to Antarctica, Geosci. Model Dev., 5, 1273–1295, https://doi.org/10.5194/gmd512732012, 2012. a
Quiquet, A., Dumas, C., Ritz, C., Peyaud, V., and Roche, D. M.: The GRISLI ice sheet model (version 2.0): calibration and validation for multimillennial changes of the Antarctic ice sheet, Geosci. Model Dev., 11, 5003–5025, https://doi.org/10.5194/gmd1150032018, 2018. a, b
Robinson, A.: palmaice/yelmo: Solver Stability v1.0, Zenodo [code], https://doi.org/10.5281/zenodo.5791864, 2021. a
Robinson, A., AlvarezSolas, J., Montoya, M., Goelzer, H., Greve, R., and Ritz, C.: Description and validation of the icesheet model Yelmo (version 1.0), Geosci. Model Dev., 13, 2805–2823, https://doi.org/10.5194/gmd1328052020, 2020. a, b, c, d
Schoof, C. and Hindmarsh, R. C. A.: Thinfilm flows with wall slip: an asymptotic analysis of higher order glacier flow models, Q. J. Mech. Appl. Math., 63, 73–114, 2010. a, b, c
Seroussi, H., Nowicki, S., Simon, E., AbeOuchi, A., Albrecht, T., Brondex, J., Cornford, S., Dumas, C., GilletChaulet, F., Goelzer, H., Golledge, N. R., Gregory, J. M., Greve, R., Hoffman, M. J., Humbert, A., Huybrechts, P., Kleiner, T., Larour, E., Leguy, G., Lipscomb, W. H., Lowry, D., Mengel, M., Morlighem, M., Pattyn, F., Payne, A. J., Pollard, D., Price, S. F., Quiquet, A., Reerink, T. J., Reese, R., Rodehacke, C. B., Schlegel, N.J., Shepherd, A., Sun, S., Sutter, J., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., and Zhang, T.: initMIPAntarctica: an ice sheet model initialization experiment of ISMIP6, The Cryosphere, 13, 1441–1471, https://doi.org/10.5194/tc1314412019, 2019. a, b
Seroussi, H., Nowicki, S., Payne, A. J., Goelzer, H., Lipscomb, W. H., AbeOuchi, A., Agosta, C., Albrecht, T., AsayDavis, X., Barthel, A., Calov, R., Cullather, R., Dumas, C., GaltonFenzi, B. K., Gladstone, R., Golledge, N. R., Gregory, J. M., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huybrechts, P., Jourdain, N. C., Kleiner, T., Larour, E., Leguy, G. R., Lowry, D. P., Little, C. M., Morlighem, M., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Reese, R., Schlegel, N.J., Shepherd, A., Simon, E., Smith, R. S., Straneo, F., Sun, S., Trusel, L. D., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., Zhao, C., Zhang, T., and Zwinger, T.: ISMIP6 Antarctica: a multimodel ensemble of the Antarctic ice sheet evolution over the 21st century, The Cryosphere, 14, 3033–3070, https://doi.org/10.5194/tc1430332020, 2020. a
Winkelmann, R., Martin, M. A., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISMPIK) – Part 1: Model description, The Cryosphere, 5, 715–726, https://doi.org/10.5194/tc57152011, 2011. a, b, c, d, e
 Abstract
 Introduction
 Icedynamics solvers
 Analytical stability limits
 Greenland experiments
 Discussion
 Conclusions
 Appendix A: Velocity derivation in stability analysis
 Appendix B: Approximate matrix inverse G
 Code availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Icedynamics solvers
 Analytical stability limits
 Greenland experiments
 Discussion
 Conclusions
 Appendix A: Velocity derivation in stability analysis
 Appendix B: Approximate matrix inverse G
 Code availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement