Description and evaluation of the Community Ice Sheet Model (CISM) v2.1

3.1.1 Blatter–Pattyn approximation

The basic equations of the Blatter–Pattyn approximation are

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 3-D Cartesian coordinates. These and other variables and parameters used in CISM are listed in Tables 1 and 2 for reference. 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) gives the gravitational driving force. The longitudinal and lateral shear stresses together are sometimes called membrane stresses (Hindmarsh, 2006). Neglecting membrane stress gradients leads to the much simpler SIA, and neglecting vertical shear stress gradients leads to the SSA.

Table 1Model variables defined in the text.

Download Print Version | Download XLSX

Table 2Model constant and parameter values used for simulations described in the text. The calving and sub-shelf melting parameters are used only for simulations with ice shelves; other parameters are used for all Greenland simulations.

Download Print Version | Download XLSX

The equations are discretized on a structured 3-D mesh. In the map plane, the mesh consists of rectangular cells, each with four vertices. The nz vertical levels of the mesh are based on a terrain-following sigma coordinate system, with $\mathit{\sigma }=(s-z)/H$, where H is the ice thickness. Each cell layer is treated as a 3-D hexahedral element with eight nodes. (In other words, cells and vertices are defined to lie in the map plane, whereas elements and nodes live in 3-D space.) Scalar 2-D fields such as H and s are defined at cell centers, and 3-D scalars such as the ice temperature T lie at the centers of elements. Gradients of 2-D scalars (e.g., the surface slope ∇s) live at vertices. The velocity components u and v are 3-D fields defined at nodes.

For problems on multiple processors, a cell or vertex (and the associated elements and nodes in its column) may be either locally owned or part of a computational halo. Each processor is responsible for computing u and v at its locally owned nodes. Any cell that contains one or more locally owned nodes and has H exceeding a threshold thickness (typically 1 m) is considered dynamically active, as are the elements in its column. Likewise, any vertex of an active cell is active, as are the nodes in its column.

The effective viscosity η is defined in each active element by

where A is the temperature-dependent rate factor in Glen's flow law (Glen, 1955), and ${\dot{\mathit{\epsilon }}}_{\mathrm{e}}$ is the effective strain rate, given in the BP approximation by

where the components of the symmetric strain rate tensor are

The rate factor A is given by an Arrhenius relationship:

where T^* is the absolute temperature corrected for the dependence of the melting point on pressure ($T^{*}=T+\mathrm{8.7}\times {\mathrm{10}}^{-\mathrm{4}}(s-z)$, with T in Kelvin), a is a temperature-independent material constant from Paterson and Budd (1982), Q is the activation energy for creep, and R is the universal gas constant.

The coupled partial differential equations (PDEs) (1) are discretized using the finite-element method (e.g., Hughes, 2000; Huebner et al., 2001). This method is more often applied to unstructured grids but was chosen for CISM's rectangular grid because of its robustness and natural treatment of boundary conditions (Dukowicz et al., 2010). The PDEs, with appropriate boundary conditions, are converted to a system of algebraic equations by dividing the full domain into subdomains (i.e., hexahedral elements), representing the velocity solution on each element, and integrating over elements. The solution is approximated as a sum over trilinear basis functions φ. Each active node is associated with a basis function whose value is φ=1 at that node, with φ=0 at all other nodes. The solution at a point within an element can be expanded in terms of basis functions and nodal values:

where the sum is over the nodes of the element; u_n and v_n are nodal values of the solution; and φ_n varies smoothly between 0 and 1 within the element.

Glissade's finite-element scheme is formally equivalent to that described by Perego et al. (2012). Equation (1) can be written as

where

Following Perego et al. (2012), these equations can be rewritten in weak form. This is done by multiplying Eq. (7) by the basis functions and integrating over the domain, using integration by parts to eliminate the second derivative:

where Ω represents the domain volume, Γ_B denotes the lower boundary, Γ_L denotes the lateral boundary (e.g., the calving front of an ice shelf), β is a basal traction parameter, p is the pressure at the lateral boundary, and n₁ and n₂ are components of the normal to Γ_L. These equations can also be obtained from a variational principle as described by Dukowicz et al. (2010). The four terms in Eq. (9) represent internal ice stresses, basal friction, lateral pressure, and the gravitational driving force, respectively.

At the basal boundary, we assume a friction 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 non-negative friction parameter that is defined at each vertex and can vary spatially. For some basal sliding laws (see Sect. 3.4), β depends on the basal velocity.

The lateral pressure p applies at marine-terminating boundaries. The net pressure is equal to the pressure directed outward from the ice toward the ocean by the ice, minus the (smaller) pressure directed inward from the ocean by the hydrostatic water pressure. The outward pressure is found by integrating ρ_ig(s−z)dz from (s−H) to s and then dividing by H; it is given by

The inward pressure is found by integrating ρ_ogzdz (where ρ_o is the seawater pressure) from (s−H) to 0 and then dividing by (H−s); it is given by

Assuming hydrostatic balance (i.e., a floating margin), we have $s-H=({\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{o}})H$, in which case Eqs. (11) and (12) can be combined to give the net pressure

directed from the ice to the ocean. However, Eqs. (11) and (12) are used in the code because they are valid at both grounded and floating marine margins. They are combined to give the net pressure $p=p_{\text{out}}-p_{\text{in}}$ that is integrated over vertical cliff faces Γ_L in Eq. (9).

The gravitational forcing terms require evaluating the gradients $\partial s/\partial x$ and $\partial s/\partial y$ at each active vertex. For vertex (i,j) (which lies at the upper right corner of cell (i,j)), these are computed using a second-order-accurate centered approximation:

and similarly for $\partial s/\partial y$. This is equivalent to computing $\partial s/\partial x$ at the edge midpoints adjacent to the vertex in the y direction, and then averaging the two edge gradients to the vertex. In some cases (e.g., near steep margins), Eq. (14) can lead to checkerboard noise (i.e., 2-D patterns of alternating positive and negative deviations in s at the grid scale), because centered averaging of s permits a computational mode. To damp this noise, Glissade also supports upstream gradient calculations that do not have this computational mode but are formally less accurate.

Equation (14) is ambiguous at the ice margin, where at least one of the four cells neighboring a vertex is ice-free. One option is to include all cells, including ice-free cells, in the gradient. This approach works reasonably well (albeit with numerical errors; see, e.g., Van den Berg et al., 2006) for land-based ice but can give large gradients and excessive ice speeds at floating shelf margins. A second option is to include only ice-covered cells in the gradient. For example, suppose cells (i,j) and $(i,j+\mathrm{1})$ have ice, but cells $(i+\mathrm{1},j)$ and $(i+\mathrm{1},j+\mathrm{1})$ are ice-free. Then, lacking the required information to compute x gradients at the adjacent edges, we would set $\partial s/\partial x=\mathrm{0}$. With a y gradient available at one adjacent edge, we would have $\partial s/\partial y=\left(s\right(i,j+\mathrm{1})-s(i,j\left)\right)/\left(\mathrm{2}\mathrm{\Delta }y\right)$. This option works well at shelf margins but tends to underestimate gradients at land margins. A third, hybrid option is to compute a gradient at edges where either (1) both adjacent cells are ice-covered or (2) one cell is ice-covered and lies higher in elevation than a neighboring ice-free land cell. Thus, gradients are set to zero at edges where an ice-covered cell (either grounded or floating) lies above ice-free ocean, or where an ice-covered land cell lies below ice-free land (i.e., a nunatak). Since this option works well for both land and shelf margins, it is the default.

Given T, s, H, and an initial guess for u and v, the problem is to solve Eq. (1) for u and v at each active node. This problem can be written as

or more fully,

In Glissade, A is always symmetric and positive definite.

Since A depends (through η and possibly β) on u and v, the problem is nonlinear and must be solved iteratively. For each nonlinear iteration, Glissade computes the 3-D η field based on the current guess for the velocity field and solves a linear problem of the form (16). Then, η is updated and the process is repeated until the solution converges to within a given tolerance. This procedure is known as Picard iteration.

Appendix A describes how the terms in Eq. (9) are summed over elements and assembled into the matrix A and vector b. Section 3.1.5 discusses solution methods for the inner linear problem and the outer nonlinear problem.

3.1.2 Shallow-ice approximation

The shallow-ice equations follow from the Blatter–Pattyn equations if membrane stresses are neglected. The SIA analogs of Eq. (1) are

The SIA could be considered a special case of the Blatter–Pattyn equations and solved using the same finite-element methods, ignoring the horizontal stresses. With these methods, however, each ice column cannot be solved independently, because each node is linked to its horizontal neighbors by terms that arise during element assembly. As a result, a finite-element-based SIA solver is not very efficient.

Instead, Glissade has an efficient local SIA solver. The solver is local in the sense that u and v in each column are found independently of u and v in other columns. It resembles Glide's SIA solver as described by Payne and Dongelmans (1997) and Rutt et al. (2009), except that Glide incorporates the velocity solution in a diffusion equation for ice thickness, whereas Glissade's SIA solver computes u and v only, with thickness evolution handled separately as described in Sect. 3.3.

For small problems that can run on one processor, there is no particular advantage to using Glissade's local SIA solver in place of Glide. Glide's implicit thickness solver permits a longer time step and thus is more efficient for problems run in serial. For whole-ice-sheet problems, however, Glissade's parallel solver can hold more data in memory and may have faster throughput, simply because it can run on tens to hundreds of processors.

3.1.3 Shallow-shelf approximation

The SSA equations can be derived by vertically integrating the BP equations, given the assumption of small basal shear stress and vertically uniform velocity. The shallow-shelf analog of Eq. (1) is

where $\overline{\mathit{\eta }}$ is the vertically averaged effective viscosity, and u and v are vertically averaged velocity components. The SSA equations in weak form resemble Eqs. (8) and (9) but without the vertical shear terms. Thus, the SSA matrix and RHS can be assembled using the methods described in Sect. 3.1.1 and Appendix A, using 2-D rectangular (instead of 3-D hexahedral) elements and omitting the vertical shear terms. The effective viscosity is defined as in Eq. (2) but with a vertically averaged flow factor and with Eq. (3) replaced by

The element matrices are 2-D analogs of Eqs. (A3)–(A6), with vertical derivatives excluded.

3.1.4 Depth-integrated-viscosity approximation

Goldberg (2011) derived a higher-order stress approximation that in most cases is similar in accuracy to BP but is much cheaper to solve, because (as with the SSA) the matrix system is assembled and solved in two dimensions. 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. Arthern and Williams (2017) used a similar model, also based on Goldberg (2011), to simulate ice flow in the Amundsen Sea sector of West Antarctica. To our knowledge, CISM is the first model to adapt this scheme for multimillennial whole-ice-sheet simulations. Here, we summarize the DIVA stress balance and describe Glissade's solution method. Section 4.2 compares DIVA results to BP results for the ISMIP-HOM experiments (Pattyn et al., 2008).

Dukowicz et al. (2010) showed that the Blatter–Pattyn equations (Eq. 1) and associated boundary conditions can be derived by taking the first variation of a functional and setting it to zero. The functional includes terms that depend on the effective strain rate (Eq. 3), which can be written as

where u_x denotes the partial derivative $\partial u/\partial x$, and similarly for other derivatives. Goldberg (2011) derived an alternate set of equations using a similar functional but with the horizontal velocity gradients in the effective strain rate replaced by their vertical averages:

where

The resulting equations of motion, analogous to Eq. (1), are

where the depth-integrated effective viscosity is given by

The vertical shear terms still contain η(z).

Since the horizontal stress terms are depth-independent, Eq. (23) can be integrated in the vertical to give

where the boundary conditions at b and s have been used to evaluate the vertical stress terms, with a sliding law of the form Eq. (10). Equation (25) is equivalent to the SSA stress balance (Eq. 18), except for the addition of basal stress terms. Thus, the methods used to assemble and solve the SSA equations can be applied to the DIVA equations to find the mean velocity components $\overline{u}$ and $\overline{v}$.

In order to solve Eq. (25), however, the basal stress terms must be rewritten in terms of $\overline{u}$ and $\overline{v}$. Following Goldberg (2011), we show how this is done for the x component of velocity; the y component is analogous. The x component of Eq. (25) can be rearranged and divided by H:

From Eq. (23), the RHS of Eq. (26) is just −(ηu_z)_z, giving

Integrating Eq. (27) from z to the surface s gives

Dividing by η and integrating from b to z gives u(z) in terms of u_b and η(z):

Following Arthern et al. (2015), we can define some useful integrals F_n as

In this notation, the surface velocity is related to the bed velocity by

Integrating u(z) from the bed to the surface gives the depth-averaged mean velocity $\overline{u}$:

We can think of F₂ as a depth-integrated inverse viscosity. The less viscous the ice, the greater the value of F₂, and the greater the difference between $\overline{u}$ and u_b.

Given Eq. (32), we can replace βu_b and βv_b in Eq. (25) with ${\mathit{\beta }}_{\mathrm{eff}}\overline{u}$ and ${\mathit{\beta }}_{\mathrm{eff}}\overline{v}$, respectively, where

For a frozen bed ($u_{\mathrm{b}}=v_{\mathrm{b}}=\mathrm{0}$, with nonzero basal stress τ_bx), the βu_b term on the RHS of Eq. (29) is replaced by τ_bx, leading to

Then, the basal stress terms in Eq. (25) can be replaced by ${\mathit{\beta }}_{\text{eff}}^{\text{frz}}\overline{u}$ and ${\mathit{\beta }}_{\text{eff}}^{\text{frz}}\overline{v}$, where

With these substitutions, Eq. (25) can be written in terms of mean velocities $\overline{u}$ and $\overline{v}$ and is fully analogous to the SSA.

In order to compute $\overline{\mathit{\eta }}$, the effective viscosity η(z) must be evaluated at each level. We use Eq. (2) with the effective strain rate (Eq. 21). The horizontal terms in ${\dot{\mathit{\epsilon }}}_{\text{DIVA}}^{\mathrm{2}}$ are found using the mean velocities from the previous iteration. The vertical shear terms are computed using Eq. (28):

where both η(z) and τ_b are from the previous iteration.

The DIVA solution procedure can be summarized as follows:

1. Starting with the current guess for the velocity field (u,v), assemble the DIVA solution matrix in analogy to the SSA solution matrix, with the following DIVA-specific computations:

   • Compute η(z) using Eqs. (2), (21), and (36), and integrate to find the depth-averaged $\overline{\mathit{\eta }}$.

   • Evaluate numerically the vertical integrals in Eqs. (29) and (30).

   • Interpolate F₂ to vertices and compute β_eff.

2. Solve the matrix system for $\overline{u}$ and $\overline{v}$.

3. Given $(\overline{u},\overline{v})$ at each vertex, use Eq. (32) to find (u_b,v_b), and then use Eq. (29) to find u(z) at each level.

4. Iterate to convergence.

The code is initialized with $(u,v)=\mathrm{0}$ everywhere, and thereafter each iterated solution starts with (u,v) from the previous time step.

DIVA assumes that horizontal gradients of membrane stresses are vertically uniform. It is less accurate than the BP approximation in regions of slow sliding and rough bed topography, where the horizontal gradients of membrane stresses vary strongly with height, and may also be less accurate where vertical temperature gradients are large. The accuracy of DIVA compared to BP is discussed in Sect. 4.2 with reference to the ISMIP-HOM benchmark experiments.

3.1.5 Solving the matrix system

After assembling the matrices and right-hand side vectors for the chosen approximation (SSA, DIVA, or BP), we solve the linear problem of Eq. (16). Glissade supports three kinds of solvers: (1) a native Fortran preconditioned conjugate gradient (PCG) solver, (2) links to Trilinos solver libraries, and (3) links to the Sparse Linear Algebra Package (SLAP).

The native PCG solver works directly with the assembled matrices (A_uu, A_uv, A_vu, and A_vv) and the right-hand-side vectors (b_u and b_v). Instead of converting to a sparse matrix format such as compressed row storage, Glissade maintains these matrices and vectors as structured rectangular arrays with i and j indices, so that they can be handled by CISM's halo update routines. Each node location (i,j) (in 2-D) or $(i,j,k)$ (in 3-D) corresponds to a row of the matrix, and an additional index m ranges over the neighboring nodes that can have nonzero terms in the columns of that row. The maximum number of nonzero terms per row is 9 in 2-D and 27 in 3-D.

The matrices and RHS vectors are passed to either a “standard” PCG solver (Shewchuk, 1994) or a Chronopoulos–Gear PCG solver (Chronopoulos, 1986; Chronopoulos and Gear, 1989). The PCG algorithm includes two dot products (each requiring a global sum), one matrix–vector product, and one preconditioning step per iteration. For small problems, the dominant computational cost is the matrix–vector product. For large problems, however, the global sums become increasingly expensive. In this case, the Chronopoulos–Gear solver is more efficient than the standard solver, because it rearranges operations such that both global sums are done with a single MPI call.

As described by Shewchuk (1994), the convergence rate of the PCG method depends on the effectiveness of the preconditioner (i.e., a matrix M which approximates A and can be inverted efficiently). Glissade's native PCG solver has two preconditioning options:

• diagonal preconditioning, in which M consists of the diagonal terms of A and is trivial to invert, and

• shallow-ice-based preconditioning (for the BP approximation only). In this case, M includes only the terms in A that link a given node to itself and its neighbors above and below. Thus, M is tridiagonal and can be inverted efficiently. This preconditioner works well when the physics is dominated by vertical shear but not as well when membrane stresses are important. Since the preconditioner is local (it consists of independent column solves), it scales well for large problems.

The Trilinos solver (Heroux et al., 2005) uses C++ packages developed at Sandia National Laboratories. As described by Evans et al. (2012), an earlier version of CISM's higher-order dycore, known as Glam, was parallelized and linked to Trilinos. The Trilinos links developed for Glam were then adapted for Glissade's SSA, DIVA, and BP solvers. The CISM documentation gives instructions for building and linking to Trilinos and choosing appropriate solver settings.

SLAP is a set of Fortran routines for solving sparse systems of linear equations. SLAP was part of Glimmer and is used to solve for thickness evolution and temperature advection in Glide. It can also be used to solve higher-order systems in Glissade, using either GMRES or a biconjugate gradient solver. The SLAP solver, however, is a serial code unsuited for large problems.

The default linear solver is native PCG using the Chronopoulos–Gear algorithm. Although it has been found to run efficiently on platforms ranging from laptops to supercomputers, its preconditioning options are limited, so convergence can be slow for problems with dominant membrane stresses (e.g., large marine ice sheets). SLAP solvers are generally robust and efficient but are limited to one processor. Trilinos contains many solver options, but in tests to date – using a generalized minimal residual (GMRES) solver with incomplete lower–upper (ILU) preconditioning – it has been found to be slower than the native PCG solver because of the extra cost of setting up Trilinos data structures, and possibly the slower performance of C++ compared to Fortran for the solver.

The linear solver is wrapped by a nonlinear solver that does Picard iteration. Following each iteration, a new global matrix is assembled, using the latest velocity solution to compute the effective viscosity. Glissade then computes the new residual, $\mathit{r}=\mathit{b}-\mathbf{A}\mathit{x}$. If the l₂ norm of the residual, defined as $\sqrt{(\mathit{r},\mathit{r})}$, is smaller than a prescribed tolerance threshold, the nonlinear system of equations is considered solved. Otherwise, the linear solver is called again, until the solution converges or the maximum number of nonlinear iterations (typically 100) is reached. The number of nonlinear iterations per solve – and optionally, the number of linear iterations per nonlinear iteration – is written to standard output. In the event of non-convergence, a warning message is written to output, but the run continues. In some cases, convergence may improve later in the simulation, or solutions may be deemed sufficiently accurate even when not fully converged.

3.2.1 Vertical diffusion

Glissade uses a vertical σ coordinate, with $\mathit{\sigma }\equiv (s-z)/H$. Thus, the vertical diffusion terms can be written as

The central difference formulas for first derivatives at the upper and lower interfaces of layer k are

where ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_k$ is the value of σ at the midpoint of layer k, halfway between σ_k and σ_k+1. The second partial derivative, defined at the midpoint of layer k, is

Inserting Eq. (40) into (41), we obtain the vertical diffusion term:

3.2.2 Heat dissipation

In higher-order models, the internal heating rate Φ in Eq. (38) is given by the tensor product of strain rate and stress:

where τ_ij is the deviatoric stress tensor. The effective stress (cf. Eq. 3) is defined by

The stress tensor is related to the strain rate and the effective viscosity by

Dividing each side of Eq. (45) by 2η, substituting in Eq. (43), and using Eq. (44) gives

Both terms on the RHS of Eq. (46) are available to the temperature solver, since the higher-order velocity solver computes η during matrix assembly and diagnoses τ_e from η and ${\dot{\mathit{\epsilon }}}_{ij}$ at the end of the calculation.

3.2.3 Boundary conditions

The temperature T₀ at the upper boundary is set to $min(T_{\mathrm{air}},\mathrm{0})$, where the surface air temperature T_air is usually specified from data or passed from a climate model. The diffusive heat flux at the upper boundary (defined as positive up) is

where the denominator contains just one term because σ₀=0.

The lower ice boundary is more complex. For grounded ice, there are three heat sources and sinks. First, the diffusive flux from the bottom surface to the ice interior (positive up) is

Second, there is a geothermal heat flux F_g which can be prescribed as a constant (∼0.05 W m⁻²) or read from an input file. Finally, there is a frictional heat flux associated with basal sliding, given by

where τ_b and u_b are 2-D vectors of basal shear stress and basal velocity, respectively. With a friction law given by Eq. (10), this becomes

If the basal temperature T_nz<T_pmp (where T_pmp is the pressure melting point temperature), then the fluxes at the lower boundary must balance:

In other words, the energy supplied by geothermal heating and sliding friction is equal to the energy removed by vertical diffusion. If, on the other hand, T_nz=T_pmp, then the net flux is nonzero and is used to melt or freeze ice at the boundary:

where M_b is the melt rate and L is the latent heat of melting. Melting generates basal water, which may either stay in place, flow downstream (possibly replaced by water from upstream), or simply disappear from the system, depending on the basal water parameterization. While basal water is present, T_nz is held at T_pmp.

For floating ice, the basal boundary condition is simpler; T_nz is set to the freezing temperature T_f of seawater, taken as a linear function of depth based on the pressure-dependent melting point of seawater. Optionally, a melt rate can be prescribed at the lower surface (Sect. 3.6).

3.2.4 Vertical temperature solution

Equation (38) can be discretized for layer k as

where the coefficients a_k and b_k are inferred from Eq. (42), and n is the current time level. The vertical diffusion terms are evaluated at the new time level, making the discretization backward Euler (i.e., fully implicit) in time. Equation (53) can be rewritten as

where

At the lower surface, we have either a temperature boundary condition (T_nz=T_pmp for grounded ice, or T_nz=T_f for floating ice) or a flux boundary condition:

which can be rearranged to give

In each ice column, the above equations form a tridiagonal system that is solved for T_k in each layer.

Occasionally, the solution T_k can exceed T_pmp for a given layer. If so, we set T_k=T_pmp and use the extra energy to melt ice internally. This melt is assumed to drain immediately to the bed. If Eq. (52) applies, we compute M_b and adjust the basal water depth. When the basal water goes to zero, T_nz is set to the temperature of the lowest layer (less than T_pmp at the bed), so that the flux boundary condition will apply during the next time step.

3.3.1 Incremental remapping

The incremental remapping scheme has several desirable features:

• It conserves the quantities being transported (including mass and internal energy).

• It is non-oscillatory; that is, it does not create spurious ripples in the transported fields.

• It preserves tracer monotonicity. That is, it does not create new extrema in tracers such as temperature; the values at time n+1 are bounded by the values at time n. Thus, T never rises above the melting point under advection.

• It is second-order accurate in space and therefore is less diffusive than first-order schemes. The accuracy may be reduced locally to first order to preserve monotonicity.

• It is efficient for large numbers of tracers. Much of the work is geometrical and is performed only once per cell edge instead of being repeated for each field being transported.

The model's upwind scheme, like IR, is conservative, non-oscillatory, and monotonicity preserving, but because it is first order it is highly diffusive.

The IR time step is limited by the requirement that trajectories projected backward from grid cell corners are confined to the four surrounding cells; this is what is meant by incremental as opposed to general remapping. This requirement leads to an advective Courant–Friedrichs–Lewy (CFL) condition,

The remapping algorithm can be summarized as follows:

1. Given mean values of the ice thickness and tracer fields in each grid cell, construct linear approximations of these fields that preserve the mean. Limit the field gradients to preserve monotonicity.

2. Given ice velocities at grid cell vertices, identify departure regions for the fluxes across each cell edge. Divide these departure regions into triangles and compute the coordinates of the triangle vertices.

3. Integrate the thickness and tracer fields over the departure triangles to obtain the mass and energy transported across each cell edge.

4. Given this transport, update the state variables.

These steps are carried out for each of nz−1 ice layers, where nz is the number of velocity levels.

3.3.2 CFL checks

As mentioned above, the time step for explicit mass transport is limited by the advective CFL condition (Eq. 60). For ice flow parallel to ∇s, ice thickness evolution is diffusive, giving rise to a diffusive CFL condition (Bueler, 2009):

where D is the ice flow diffusivity. Flow governed by the shallow-ice approximation is subject to this diffusive CFL. The stability of Glissade's BP, DIVA, and SSA solvers, however, is limited by Eq. (60) in practice; Eq. (61) is too restrictive. For this reason, Glissade checks for advective CFL violations at each time step. Optionally, the transport equation can be adaptively subcycled within a time step to satisfy advective CFL. This can prevent the model from crashing, though possibly with a loss of accuracy.