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

We describe and evaluate version 2.1 of the Community Ice Sheet Model (CISM). CISM is a parallel, 3-D thermomechanical model, written mainly in Fortran, that solves equations for the momentum balance and the thickness and temperature evolution of ice sheets. CISM's velocity solver incorporates a hierarchy of Stokes flow approximations, including shallow-shelf, depth-integrated higher order, and 3-D higher order. CISM also includes a suite of test cases, links to third-party solver libraries, and parameterizations of physical processes such as basal sliding, iceberg calving, and sub-ice-shelf melting. The model has been verified for standard test problems, including the Ice Sheet Model Intercomparison Project for Higher-Order Models (ISMIP-HOM) experiments, and has participated in the initMIP-Greenland initialization experiment. In multimillennial simulations with modern climate forcing on a 4 km grid, CISM reaches a steady state that is broadly consistent with observed flow patterns of the Greenland ice sheet. CISM has been integrated into version 2.0 of the Community Earth System Model, where it is being used for Greenland simulations under past, present, and future climates. The code is open-source with extensive documentation and remains under active development.

1 Introduction

As mass loss from the Greenland and Antarctic ice sheets has accelerated (Shepherd et al., 2012; Church et al., 2013; Hanna et al., 2013; Shepherd et al., 2017), climate modelers have recognized the importance of dynamic ice sheet models (ISMs) for predicting future mass loss and sea level rise (Vizcaino, 2014). Meanwhile, ISMs have become more accurate and complex in their representation of ice flow dynamics. Early ISMs used either the shallow-ice approximation (SIA; Hutter, 1983) or the shallow-shelf approximation (SSA; MacAyeal, 1989). The SIA, which assumes that vertical shear stresses are dominant, is valid for slow-moving ice sheet interiors, whereas the SSA, which assumes that flow is dominated by lateral and longitudinal stresses in the horizontal plane, is valid for floating ice shelves. Neither approximation is valid for ice streams and outlet glaciers where both vertical shear and horizontal plane stresses are important. Advanced ISMs developed in recent years solve the Stokes equations or various higher-order approximations (Pattyn et al., 2008; Schoof and Hindmarsh, 2010). Among the models that solve Stokes or higher-order equations, or otherwise combine features of the SIA and SSA, are the Parallel Ice Sheet Model (PISM; Bueler and Brown, 2009; Winkelmann et al., 2011), the Ice Sheet System Model (ISSM; Larour et al., 2012), the Penn State Model (Pollard and DeConto, 2012), BISICLES (Cornford et al., 2013), Elmer-Ice (Gagliardini et al., 2013), and MPAS-Albany Land Ice (MALI; Tezaur et al., 2015; Hoffman et al., 2018). Higher-order models have performed well for standard test cases (Pattyn et al., 2008, 2012) and have been applied to many scientific problems.

Here, we describe and evaluate the Community Ice Sheet Model (CISM) version 2.1, a higher-order model that evolved from the Glimmer model (Rutt et al., 2009). The current name reflects the model's evolution as a component of the Community Earth System Model (CESM; Hurrell et al., 2013). Like Glimmer, CISM is written mainly in Fortran 90 and its extensions, to maximize efficiency and to simplify coupling to climate models. Glimmer, however, is a serial SIA model, whereas CISM is a parallel model that solves not only the SIA but also higher-order Stokes approximations.

CISM development was guided by the following goals:

• The model should be well documented and easy to install and run on a variety of platforms, ranging from laptops to local clusters to high-performance supercomputers.

• It should solve a range of Stokes approximations, including the SIA, SSA, and higher-order approximations. Velocity solvers for these approximations are included in a new dynamical core called Glissade. (The dynamical core, or “dycore”, is the part of the model that solves equations for conservation of mass, energy, and momentum.)

• It should remain backward-compatible with Glimmer, allowing continued use of the older Glide SIA dycore.

• It should be well verified for standard test cases such as the Ice Sheet Model Intercomparison Project for Higher-Order Models (ISMIP-HOM) (Pattyn et al., 2008), with a user-friendly verification framework.

• It should run efficiently – supporting whole-ice-sheet applications even on small platforms – and should scale to hundreds of processor cores, enabling century- to millennial-scale simulations with higher-order solvers at grid resolutions of ∼5 km or finer.

• It should support not only stand-alone ice sheet simulations but also coupled applications in which fields are exchanged with a global climate or Earth system model (CESM in particular).

• It should support simulations of the Greenland ice sheet, a scientific focus of CESM. Support for Antarctic applications is deferred to future model releases and publications.

• The code should be open-source, with periodic public releases.

Many of these features were present in CISM v.2.0, which was released in 2014 (Price et al., 2015).

Changes between versions 2.0 and 2.1 have been made primarily to support robust, accurate, and efficient Greenland ice sheet simulations, both as a stand-alone model and in CESM. These changes include a depth-integrated higher-order velocity solver (Sect. 3.1.4), new parameterizations of basal sliding (Sect. 3.4), iceberg calving (Sect. 3.5), and sub-ice-shelf melting (Sect. 3.6), a new build-and-test structure (Sect. 4.6), and many improvements in model numerics.

We begin with an overview of CISM, including the core model and its testing and coupling infrastructure (Sect. 2). We then describe the model dynamics and physics, focusing on the new Glissade dycore (Sect. 3). To verify the model, we present results from standard test cases (Sect. 4). We then present results from long spin-ups of the Greenland ice sheet, with and without floating ice shelves (Sect. 5). Finally, we summarize the model results and suggest directions for future work (Sect. 6).

2 Model overview

CISM is a numerical model – a collection of software libraries, utilities, and drivers – used to simulate ice sheet evolution. It is modular in design and is coded mainly in standard-compliant Fortran 90/95. CISM consists of several components:

cism_driver

is the high-level driver (i.e., the executable) that is used to run the stand-alone model in all configurations, including idealized test cases with simplified climate forcing, as well as model runs with realistic geometry and climate forcing data.

Glide

is the serial dycore based on shallow-ice dynamics. Glide solves the governing conservation equations and computes ice velocities, internal ice temperature, and ice geometry evolution. Apart from minor changes, this is the same dycore used in Glimmer and described by Rutt et al. (2009). It will not be discussed further here.

Glissade

is the dycore that solves higher-order approximations of the Stokes equations for ice flow. Glissade, unlike Glide, is fully parallel in order to take advantage of multiprocessor, high-performance architectures. It is described in detail in Sect. 3.

Glint

is the original climate model interface for Glimmer. Glint allows the core ice sheet model to be coupled to a global climate model or any other source of time-varying climate data on a lat–long grid. Glint computes the surface mass balance (SMB) on the ice sheet grid using a positive-degree-day (PDD) scheme.

Glad

is a lightweight climate model interface that has replaced Glint in CESM. The CESM coupler supports remapping and downscaling between general land-surface grids and ice sheet grids, and thus is able to send CISM an SMB that is already downscaled. The Glad interface simply sends and receives fields on the ice sheet grid, accumulating and averaging as needed based on the ice dynamic time step and the coupling interval.

Test cases

are provided for the Glide and Glissade dynamical cores. These are used to confirm that the model is working as expected and to provide a range of simple model configurations from which new users can learn about model options and create their own configurations. CISM test cases are described in Sect. 4.

Shared code

consists of modules shared by different parts of the code. Examples include modules for defining derived types, physical constants, and model parameters, and modules that parse CISM configuration files and handle data input/output (I/O).

In order to reduce development effort, CISM runs on a structured rectangular grid and thus lacks the flexibility of models that run on variable-resolution or adaptive meshes (e.g., ISSM and BISICLES). Although CISM includes several common Stokes approximations, it does not solve the more complex “full Stokes” equations (Pattyn et al., 2008).

CISM is distributed as source code and therefore requires a reasonably complete build environment to compile the model. For UNIX- and LINUX-based systems, the CMake (http://www.cmake.org/, last access: 20 January 2019) system is used to build the model. Sample build scripts for a number of standard architectures are included, as are working build scripts for a number of high-performance-computing architectures including Cheyenne (NCAR Computational and Information Systems Laboratory), Titan (Oak Ridge Leadership Computing Facility), Edison (National Energy Research Scientific Computing Center), and Cartesius (the Dutch national supercomputer).

The source code can be obtained by downloading a released version from the CISM website or by cloning the code from a public git repository; see Sect. . In either case, a Fortran 90 compiler is required. Other software dependencies include the NetCDF library (https://www.unidata.ucar.edu/software/netcdf/, last access: 20 January 2019) (used for data I/O) and a Python (http://www.python.org, last access: 20 January 2019) distribution (used to analyze dependencies and to automatically generate parts of the code) with several specific Python modules. Parallel builds require Message Passing Interface (MPI), and users desiring access to Trilinos packages (Heroux et al., 2005) will need to build Trilinos and then link it to CISM. Finally, CMake and Gnu Make are needed to compile the code and link to the various third-party libraries. The CISM documentation contains detailed instructions for downloading and building the code.

CISM is run by specifying the names of the executable (usually cism_driver) and a configuration file. Typically, the configuration file includes the input and output filenames, the grid dimensions, the time step and length of the run, the dycore (Glide or Glissade), and various options and parameter values appropriate for a given application. If not set in the config file, each option or parameter takes a default value. Supported options are described in the model documentation.

Multiple input, forcing, and output files can be specified, containing any subset of a large number of global scalars and fields (1-D, 2-D, and 3-D). If a given field is “loadable” and is present in the input file, it is read automatically at startup; otherwise, it is set to a default value. Loadable fields include the initial ice thickness and temperature, bedrock topography, surface mass balance, surface air temperature, and geothermal heat flux. Forcing files are input files that are read at every time step (not just at initialization) so that time-dependent forcing can be applied during a simulation.

Each output file includes a user-chosen set of variables (listed in the configuration file) and can contain multiple time slices, written at any frequency. A special kind of I/O file is the restart file, which includes all the fields needed to restart the model exactly. Whatever configuration options are chosen, model results are exactly reproducible (i.e., bit for bit) for a given computer platform and processor count, regardless of how many times a simulation is stopped and restarted.

Some basic information is sent to standard output during the run, and more verbose output is written at regular intervals to a log file. The log file lists the options and parameter values chosen for the run and also notes the simulation time when CISM reads or writes I/O. In addition, the log file includes diagnostic information about the global state of the model (e.g., the total ice area and volume, total surface and basal mass balance, and maximum surface and basal ice speeds), along with vertical profiles of ice speed and temperature for a user-chosen grid point.

CISM2.1 has been implemented in CESM version 2.0, released in June 2018. Earlier versions of CESM supported one-way forcing of the Greenland ice sheet (using SIA dynamics) by the SMB computed in CESM's land model (Lipscomb et al., 2013). CESM2.0 extends this capability by supporting conservative, interactive coupling between ice sheets and the land and atmosphere. Coupling of ice sheets with the ocean is not yet supported but is under development. CESM2.0 does not have an interactive Antarctic ice sheet, in part because of the many scientific and technical issues associated with ice sheet–ocean coupling. However, CISM2.1 includes many of the features needed to simulate marine ice sheets, and a developmental model version (including a grounding-line parameterization to be included in a near-future code release) has been used for Antarctic simulations.

3 Model dynamics and physics

CISM includes a parallel, higher-order dynamical core called Glissade, which solves equations for conservation of momentum (i.e., an appropriate approximation of Stokes flow), mass, and internal energy. Glissade numerics differ substantially from Glide numerics:

• Velocity: Glide solves the SIA only, but Glissade can solve several Stokes approximations, including the SIA, SSA, a depth-integrated higher-order approximation based on Goldberg (2011), and a 3-D higher-order approximation based on Blatter (1995) and Pattyn (2003). Glide uses finite differences, whereas the Glissade velocity solvers use finite-element methods.

• Temperature: To evolve the ice temperature, Glide solves a prognostic equation that incorporates horizontal advection as well as vertical heat diffusion and internal dissipation. In Glissade, temperature advection is handled by the transport scheme, and a separate module solves for vertical diffusion and internal dissipation in each column.

• Mass and tracer transport: Glide solves an implicit diffusion equation for mass transport, incorporating shallow-ice velocities. Glissade solves explicit equations for horizontal transport of mass (i.e., ice thickness) and tracers (e.g., ice temperature) using either an incremental remapping scheme (Lipscomb and Hunke, 2004) or a simpler first-order upwind scheme. Horizontal transport is followed by vertical remapping to terrain-following sigma coordinates.

Glissade numerics are described in detail below.

3.1 Velocity solvers

Glissade computes the ice velocity by solving approximate Stokes equations, given the surface elevation, ice thickness, ice temperature, and relevant boundary conditions. Section 3.1.1 describes the solution method for the Blatter–Pattyn (BP; Blatter, 1995; Pattyn, 2003) approximation, which is the most sophisticated and accurate solver in CISM. Subsequent sections discuss simpler approximations.

3.1.1 Blatter–Pattyn approximation

The basic equations of the Blatter–Pattyn approximation are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial x}}\left[\mathrm{2}\mathit{\eta }\left(\mathrm{2}{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}\right)\right]+{\displaystyle \frac{\partial }{\partial y}}\left[\mathit{\eta }\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\right]+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\eta }{\displaystyle \frac{\partial u}{\partial z}}\right)\\ {\displaystyle }& {\displaystyle }={\mathit{\rho }}_{\mathrm{i}}g{\displaystyle \frac{\partial s}{\partial x}},\\ {\displaystyle }& {\displaystyle \frac{\partial }{\partial y}}\left[\mathrm{2}\mathit{\eta }\left(\mathrm{2}{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial u}{\partial x}}\right)\right]+{\displaystyle \frac{\partial }{\partial x}}\left[\mathit{\eta }\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\right]+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\eta }{\displaystyle \frac{\partial v}{\partial z}}\right)\\ \text{(1)}& {\displaystyle }& {\displaystyle }={\mathit{\rho }}_{\mathrm{i}}g{\displaystyle \frac{\partial s}{\partial y}},\end{array}$$

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 1 Model variables defined in the text.

Table 2 Model 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.

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

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathit{\eta }\equiv {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}^{\frac{-\mathrm{1}}{n}}{\dot{\mathit{\epsilon }}}_{\mathrm{e}}^{\frac{\mathrm{1}-n}{n}},\end{array}$$

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

$$\begin{array}{}\text{(3)}& {\displaystyle }{\dot{\mathit{\epsilon }}}_{\mathrm{e}}^{\mathrm{2}}={\dot{\mathit{\epsilon }}}_{xx}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{yy}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{xx}{\dot{\mathit{\epsilon }}}_{yy}+{\dot{\mathit{\epsilon }}}_{xy}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{xz}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{yz}^{\mathrm{2}},\end{array}$$

where the components of the symmetric strain rate tensor are

$$\begin{array}{}\text{(4)}& {\displaystyle }{\dot{\mathit{\epsilon }}}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_{\mathrm{i}}}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_{\mathrm{i}}}}\right).\end{array}$$

The rate factor A is given by an Arrhenius relationship:

$$\begin{array}{}\text{(5)}& {\displaystyle }A\left(T^{*}\right)=a{e}^{-Q/\mathrm{R}{T}^{*}},\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }u(x,y,z)=\sum _n{\mathit{\phi }}_{\mathrm{n}}(x,y,z)u_{\mathrm{n}},\\ \text{(6)}& {\displaystyle }& {\displaystyle }v(x,y,z)=\sum _n{\mathit{\phi }}_{\mathrm{n}}(x,y,z)v_{\mathrm{n}},\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }-\mathbf{\nabla }\cdot \left(\mathrm{2}\mathit{\eta }{\dot{\mathit{\epsilon }}}_{\mathrm{1}}\right)=-{\mathit{\rho }}_{\mathrm{i}}g{\displaystyle \frac{\partial s}{\partial x}},\\ \text{(7)}& {\displaystyle }& {\displaystyle }-\mathbf{\nabla }\cdot \left(\mathrm{2}\mathit{\eta }{\dot{\mathit{\epsilon }}}_{\mathrm{2}}\right)=-{\mathit{\rho }}_{\mathrm{i}}g{\displaystyle \frac{\partial s}{\partial y}},\end{array}$$

where

$$\begin{array}{}\text{(8)}& {\displaystyle }{\dot{\mathit{\epsilon }}}_{\mathrm{1}}=\left[\begin{array}{c}\mathrm{2}{\dot{\mathit{\epsilon }}}_{xx}+{\dot{\mathit{\epsilon }}}_{yy}\\ {\dot{\mathit{\epsilon }}}_{xy}\\ {\dot{\mathit{\epsilon }}}_{xz}\end{array}\right],{\dot{\mathit{\epsilon }}}_{\mathrm{2}}=\left[\begin{array}{c}{\dot{\mathit{\epsilon }}}_{xy}\\ {\dot{\mathit{\epsilon }}}_{xx}+\mathrm{2}{\dot{\mathit{\epsilon }}}_{yy}\\ {\dot{\mathit{\epsilon }}}_{yz}\end{array}\right].\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathrm{\Omega }}{\int }\mathrm{2}\mathit{\eta }{\dot{\mathit{\epsilon }}}_{\mathrm{1}}(u,v)\cdot \mathrm{\nabla }{\mathit{\phi }}_{\mathrm{1}}\text{ d}\mathrm{\Omega }+\underset{{\mathrm{\Gamma }}_{\mathrm{B}}}{\int }\mathit{\beta }u{\mathit{\phi }}_{\mathrm{1}}\text{ d}\mathrm{\Gamma }+\underset{{\mathrm{\Gamma }}_{\mathrm{L}}}{\int }p{n}_{\mathrm{1}}{\mathit{\phi }}_{\mathrm{1}}\text{ d}\mathrm{\Gamma }\\ {\displaystyle }& {\displaystyle }+\underset{\mathrm{\Omega }}{\int }{\mathit{\rho }}_{\mathrm{i}}g{\displaystyle \frac{\partial s}{\partial x}}{\mathit{\phi }}_{\mathrm{1}}\text{ d}\mathrm{\Omega }=\text{0},\\ {\displaystyle }& {\displaystyle }\underset{\mathrm{\Omega }}{\int }\mathrm{2}\mathit{\eta }{\dot{\mathit{\epsilon }}}_{\mathrm{2}}(u,v)\cdot \mathrm{\nabla }{\mathit{\phi }}_{\mathrm{2}}\text{ d}\mathrm{\Omega }+\underset{{\mathrm{\Gamma }}_{\mathrm{B}}}{\int }\mathit{\beta }u{\mathit{\phi }}_{\mathrm{2}}\text{ d}\mathrm{\Gamma }+\underset{{\mathrm{\Gamma }}_{\mathrm{L}}}{\int }p{n}_{\mathrm{2}}{\mathit{\phi }}_{\mathrm{2}}\text{ d}\mathrm{\Gamma }\\ \text{(9)}& {\displaystyle }& {\displaystyle }+\underset{\mathrm{\Omega }}{\int }{\mathit{\rho }}_{\mathrm{i}}g{\displaystyle \frac{\partial s}{\partial y}}{\mathit{\phi }}_{\mathrm{2}}\text{ d}\mathrm{\Omega }=\text{0},\end{array}$$

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

$$\begin{array}{}\text{(10)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{b}}=\mathit{\beta }{\mathit{u}}_{\mathrm{b}},\end{array}$$

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

$$\begin{array}{}\text{(11)}& {\displaystyle }{p}_{\text{out}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}gH}{\mathrm{2}}}.\end{array}$$

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

$$\begin{array}{}\text{(12)}& {\displaystyle }{p}_{\text{in}}={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{o}}g(s-H{)}^{\mathrm{2}}}{\mathrm{2}H}}.\end{array}$$

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

$$\begin{array}{}\text{(13)}& {\displaystyle }p={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}gH}{\mathrm{2}}}\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{w}}}}\right),\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial s}{\partial x}}(i,j)=\\ \text{(14)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{s(i+\mathrm{1},j+\mathrm{1})+s(i+\mathrm{1},j)-s(i,j+\mathrm{1})-s(i,j)}{\mathrm{2}\mathrm{\Delta }x}},\end{array}$$

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

$$\begin{array}{}\text{(15)}& {\displaystyle }\mathbf{A}\mathit{x}=\mathit{b},\end{array}$$

or more fully,

$$\begin{array}{}\text{(16)}& {\displaystyle }\begin{array}{c}\left[\begin{array}{cc}{\mathbf{A}}_{uu}& {\mathbf{A}}_{uv}\\ {\mathbf{A}}_{vu}& {\mathbf{A}}_{vv}\end{array}\right]\left[\begin{array}{c}\mathit{u}\\ \mathit{v}\end{array}\right]=\left[\begin{array}{c}{\mathit{b}}_u\\ {\mathit{b}}_v\end{array}\right].\end{array}\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\eta }{\displaystyle \frac{\partial u}{\partial z}}\right)={\mathit{\rho }}_{\mathrm{i}}g{\displaystyle \frac{\partial s}{\partial x}},\\ \text{(17)}& {\displaystyle }& {\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\eta }{\displaystyle \frac{\partial v}{\partial z}}\right)={\mathit{\rho }}_{\mathrm{i}}g{\displaystyle \frac{\partial s}{\partial y}}.\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial x}}\left[\mathrm{2}\overline{\mathit{\eta }}H\left(\mathrm{2}{\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}\right)\right]+{\displaystyle \frac{\partial }{\partial y}}\left[\overline{\mathit{\eta }}H\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\right]\\ {\displaystyle }& {\displaystyle }={\mathit{\rho }}_{\mathrm{i}}gH{\displaystyle \frac{\partial s}{\partial x}},\\ {\displaystyle }& {\displaystyle \frac{\partial }{\partial y}}\left[\mathrm{2}\overline{\mathit{\eta }}H\left(\mathrm{2}{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\partial u}{\partial x}}\right)\right]+{\displaystyle \frac{\partial }{\partial x}}\left[\overline{\mathit{\eta }}H\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)\right]\\ \text{(18)}& {\displaystyle }& {\displaystyle }={\mathit{\rho }}_{\mathrm{i}}gH{\displaystyle \frac{\partial s}{\partial y}},\end{array}$$

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

$$\begin{array}{}\text{(19)}& {\displaystyle }{\dot{\mathit{\epsilon }}}_{\mathrm{e}}^{\mathrm{2}}={\dot{\mathit{\epsilon }}}_{xx}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{yy}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{xx}{\dot{\mathit{\epsilon }}}_{yy}+{\dot{\mathit{\epsilon }}}_{xy}^{\mathrm{2}}.\end{array}$$

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

$$\begin{array}{}\text{(20)}& {\displaystyle }{\dot{\mathit{\epsilon }}}_{\mathrm{BP}}^{\mathrm{2}}=u_x^{\mathrm{2}}+v_y^{\mathrm{2}}+u_x{v}_y+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left(u_y+v_x\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{u}_z^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{v}_z^{\mathrm{2}},\end{array}$$

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:

$$\begin{array}{}\text{(21)}& {\displaystyle }{\dot{\mathit{\epsilon }}}_{\mathrm{DIVA}}^{\mathrm{2}}={\overline{u}}_x^{\mathrm{2}}+{\overline{v}}_y^{\mathrm{2}}+{\overline{u}}_x{\overline{v}}_y+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\overline{u}}_y+{\overline{v}}_x\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{u}_z^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{v}_z^{\mathrm{2}},\end{array}$$

where

$$\begin{array}{}\text{(22)}& {\displaystyle }\overline{u}={\displaystyle \frac{\mathrm{1}}{H}}\underset{\mathrm{b}}{\overset{s}{\int }}u\left(z\right)\mathrm{d}z,\overline{v}={\displaystyle \frac{\mathrm{1}}{H}}\underset{\mathrm{b}}{\overset{s}{\int }}v\left(z\right)\mathrm{d}z.\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{1}}{H}}{\displaystyle \frac{\partial }{\partial x}}\left[\mathrm{2}\overline{\mathit{\eta }}H\left(\mathrm{2}{\displaystyle \frac{\partial \overline{u}}{\partial x}}+{\displaystyle \frac{\partial \overline{v}}{\partial y}}\right)\right]+{\displaystyle \frac{\mathrm{1}}{H}}{\displaystyle \frac{\partial }{\partial y}}\left[\overline{\mathit{\eta }}H\left({\displaystyle \frac{\partial \overline{u}}{\partial y}}+{\displaystyle \frac{\partial \overline{v}}{\partial x}}\right)\right]\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\eta }{\displaystyle \frac{\partial u}{\partial z}}\right)={\mathit{\rho }}_{\mathrm{i}}g{\displaystyle \frac{\partial s}{\partial x}},\\ {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{H}}{\displaystyle \frac{\partial }{\partial y}}\left[\mathrm{2}\overline{\mathit{\eta }}H\left(\mathrm{2}{\displaystyle \frac{\partial \overline{v}}{\partial y}}+{\displaystyle \frac{\partial \overline{u}}{\partial x}}\right)\right]+{\displaystyle \frac{\mathrm{1}}{H}}{\displaystyle \frac{\partial }{\partial x}}\left[\overline{\mathit{\eta }}H\left({\displaystyle \frac{\partial \overline{u}}{\partial y}}+{\displaystyle \frac{\partial \overline{v}}{\partial x}}\right)\right]\\ \text{(23)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\eta }{\displaystyle \frac{\partial v}{\partial z}}\right)={\mathit{\rho }}_{\mathrm{i}}g{\displaystyle \frac{\partial s}{\partial y}},\end{array}$$

where the depth-integrated effective viscosity is given by

$$\begin{array}{}\text{(24)}& {\displaystyle }\overline{\mathit{\eta }}={\displaystyle \frac{\mathrm{1}}{H}}\underset{\mathrm{b}}{\overset{s}{\int }}\mathit{\eta }\left(z\right)\mathrm{d}z.\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\partial }{\partial x}}\left[\mathrm{2}\overline{\mathit{\eta }}H\left(\mathrm{2}{\displaystyle \frac{\partial \overline{u}}{\partial x}}+{\displaystyle \frac{\partial \overline{v}}{\partial y}}\right)\right]+{\displaystyle \frac{\partial }{\partial y}}\left[\overline{\mathit{\eta }}H\left({\displaystyle \frac{\partial \overline{u}}{\partial y}}+{\displaystyle \frac{\partial \overline{v}}{\partial x}}\right)\right]-\mathit{\beta }{u}_{\mathrm{b}}\\ {\displaystyle }& {\displaystyle }={\mathit{\rho }}_{\mathrm{i}}gH{\displaystyle \frac{\partial s}{\partial x}},\\ {\displaystyle }& {\displaystyle \frac{\partial }{\partial y}}\left[\mathrm{2}\overline{\mathit{\eta }}H\left(\mathrm{2}{\displaystyle \frac{\partial \overline{v}}{\partial y}}+{\displaystyle \frac{\partial \overline{u}}{\partial x}}\right)\right]+{\displaystyle \frac{\partial }{\partial x}}\left[\overline{\mathit{\eta }}H\left({\displaystyle \frac{\partial \overline{u}}{\partial y}}+{\displaystyle \frac{\partial \overline{v}}{\partial x}}\right)\right]-\mathit{\beta }{v}_{\mathrm{b}}\\ \text{(25)}& {\displaystyle }& {\displaystyle }={\mathit{\rho }}_{\mathrm{i}}gH{\displaystyle \frac{\partial s}{\partial y}},\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathit{\beta }{u}_{\mathrm{b}}}{H}}={\displaystyle \frac{\mathrm{1}}{H}}{\displaystyle \frac{\partial }{\partial x}}\left(\mathrm{2}H\overline{\mathit{\eta }}\left(\mathrm{2}{\displaystyle \frac{\partial \overline{u}}{\partial x}}+{\displaystyle \frac{\partial \overline{v}}{\partial y}}\right)\right)\\ \text{(26)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{H}}{\displaystyle \frac{\partial }{\partial y}}\left(H\overline{\mathit{\eta }}\left({\displaystyle \frac{\partial \overline{u}}{\partial y}}+{\displaystyle \frac{\partial \overline{v}}{\partial x}}\right)\right)-{\mathit{\rho }}_{\mathrm{i}}g{\displaystyle \frac{\partial s}{\partial x}}.\end{array}$$

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

$$\begin{array}{}\text{(27)}& {\displaystyle \frac{\mathit{\beta }{u}_{\mathrm{b}}}{H}}=-{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\eta }{\displaystyle \frac{\partial u}{\partial z}}\right).\end{array}$$

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

$$\begin{array}{}\text{(28)}& {\displaystyle }\mathit{\eta }{\displaystyle \frac{\partial u}{\partial z}}={\displaystyle \frac{\mathit{\beta }{u}_{\mathrm{b}}\left(s-z\right)}{H}}.\end{array}$$

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

$$\begin{array}{}\text{(29)}& {\displaystyle }u\left(z\right)=u_{\mathrm{b}}+{\displaystyle \frac{\mathit{\beta }{u}_{\mathrm{b}}}{H}}\underset{\mathrm{b}}{\overset{z}{\int }}{\displaystyle \frac{\left(s-z^{\prime }\right)}{\mathit{\eta }\left(z^{\prime }\right)}}\mathrm{d}{z}^{\prime }.\end{array}$$

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

$$\begin{array}{}\text{(30)}& {\displaystyle }{F}_{\mathrm{n}}\equiv \underset{\mathrm{b}}{\overset{s}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\eta }}}{\left({\displaystyle \frac{s-z}{H}}\right)}^n\mathrm{d}z.\end{array}$$

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

$$\begin{array}{}\text{(31)}& {\displaystyle }{u}_{\mathrm{s}}=u_{\mathrm{b}}\left(\mathrm{1}+\mathit{\beta }{F}_{\mathrm{1}}\right).\end{array}$$

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

$$\begin{array}{}\text{(32)}& {\displaystyle }\overline{u}=u_{\mathrm{b}}\left(\mathrm{1}+\mathit{\beta }{F}_{\mathrm{2}}\right).\end{array}$$

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

$$\begin{array}{}\text{(33)}& {\displaystyle }{\mathit{\beta }}_{\text{eff}}={\displaystyle \frac{\mathit{\beta }}{\mathrm{1}+\mathit{\beta }{F}_{\mathrm{2}}}}.\end{array}$$

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

$$\begin{array}{}\text{(34)}& {\displaystyle }\overline{u}={\mathit{\tau }}_{\mathrm{b}x}{F}_{\mathrm{2}}.\end{array}$$

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

$$\begin{array}{}\text{(35)}& {\displaystyle }{\mathit{\beta }}_{\text{eff}}^{\text{frz}}={\displaystyle \frac{\mathrm{1}}{F_{\mathrm{2}}}}.\end{array}$$

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):

$$\begin{array}{}\text{(36)}& {\displaystyle }{u}_z\left(z\right)={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{b}x}(s-z)}{\mathit{\eta }\left(z\right)H}},v_z\left(z\right)={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{b}y}(s-z)}{\mathit{\eta }\left(z\right)H}},\end{array}$$

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 Temperature solver

The thermal evolution of the ice sheet is given by

$$\begin{array}{}\text{(37)}& {\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{k_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}}}{\mathrm{\nabla }}^{\mathrm{2}}T-\mathbf{u}\cdot \mathbf{\nabla }T-w{\displaystyle \frac{\partial T}{\partial z}}+{\displaystyle \frac{\mathrm{\Phi }}{{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}}},\end{array}$$

where T is the temperature in ^∘C, $\mathit{u}=(u,v)$ is the horizontal velocity, w is the vertical velocity, k_i is the thermal conductivity of ice, c_i is the specific heat of ice, ρ_i is the ice density, and Φ is the rate of heating due to internal deformation and dissipation. This equation describes the conservation of internal energy under horizontal and vertical diffusion (the first term on the RHS), horizontal and vertical advection (the second and third terms, respectively), and internal heat dissipation (the last term). Unlike Glide, which solves this equation in a single calculation, Glissade divides the temperature evolution into separate advection and diffusion/dissipation components. The temperature module uses finite-difference methods to solve for vertical diffusion and internal dissipation in each ice column:

$$\begin{array}{}\text{(38)}& {\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{k_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{i}}c}}{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial z^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{\Phi }}{{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}}},\end{array}$$

as described in this section. The advective part of Eq. (37) is described in Sect. 3.3. Horizontal diffusion is assumed to be negligible compared to vertical diffusion, giving ${\mathrm{\nabla }}^{\mathrm{2}}T\simeq \frac{{\partial }^{\mathrm{2}}T}{\partial z^{\mathrm{2}}}$.

In Glissade's vertical discretization of temperature, T is located at the midpoints of the (nz−1) layers, staggered relative to the velocity. This staggering makes it easier to conserve internal energy under transport, where the internal energy in an ice column is equal to the sum over layers of ρ_ic_iTΔz, with layer thickness Δz. The upper surface skin temperature is denoted by T₀ and the lower surface skin temperature by T_nz, giving a total of nz+1 temperature values in each column.

The following sections describe how the terms in Eq. (38) are computed, how the boundary conditions are specified, and how the equation is solved.

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

$$\begin{array}{}\text{(39)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial z^{\mathrm{2}}}}={\displaystyle \frac{\mathrm{1}}{H^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial {\mathit{\sigma }}^{\mathrm{2}}}}.\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\left.{\displaystyle \frac{\partial T}{\partial \mathit{\sigma }}}\right|}_{{\mathit{\sigma }}_k}={\displaystyle \frac{T_k-T_{k-\mathrm{1}}}{{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_k-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{k-\mathrm{1}}}},\\ \text{(40)}& {\displaystyle }& {\displaystyle }{\left.{\displaystyle \frac{\partial T}{\partial \mathit{\sigma }}}\right|}_{{\mathit{\sigma }}_{k+\mathrm{1}}}={\displaystyle \frac{T_{k+\mathrm{1}}-T_k}{{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{k+\mathrm{1}}-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_k}},\end{array}$$

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

$$\begin{array}{}\text{(41)}& {\displaystyle }{\left.{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial {\mathit{\sigma }}^{\mathrm{2}}}}\right|}_{{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_k}={\displaystyle \frac{{\left.\frac{\partial T}{\partial \mathit{\sigma }}\right|}_{{\mathit{\sigma }}_{k+\mathrm{1}}}-{\left.\frac{\partial T}{\partial \mathit{\sigma }}\right|}_{{\mathit{\sigma }}_k}}{{\mathit{\sigma }}_{k+\mathrm{1}}-{\mathit{\sigma }}_k}}.\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\left.{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial {\mathit{\sigma }}^{\mathrm{2}}}}\right|}_{{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_k}={\displaystyle \frac{T_{k-\mathrm{1}}}{\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_k-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{k-\mathrm{1}}\right)\left({\mathit{\sigma }}_{k+\mathrm{1}}-{\mathit{\sigma }}_k\right)}}\\ {\displaystyle }& {\displaystyle }-T_k\left({\displaystyle \frac{\mathrm{1}}{\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_k-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{k-\mathrm{1}}\right)\left({\mathit{\sigma }}_{k+\mathrm{1}}-{\mathit{\sigma }}_k\right)}}+{\displaystyle \frac{\mathrm{1}}{\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{k+\mathrm{1}}-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_k\right)\left({\mathit{\sigma }}_{k+\mathrm{1}}-{\mathit{\sigma }}_k\right)}}\right)\\ \text{(42)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{T_{k+\mathrm{1}}}{\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{k+\mathrm{1}}-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_k\right)\left({\mathit{\sigma }}_{k+\mathrm{1}}-{\mathit{\sigma }}_k\right)}}.\end{array}$$

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:

$$\begin{array}{}\text{(43)}& {\displaystyle }\mathrm{\Phi }={\dot{\mathit{\epsilon }}}_{ij}{\mathit{\tau }}_{ij},\end{array}$$

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

$$\begin{array}{}\text{(44)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{e}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\tau }}_{ij}{\mathit{\tau }}_{ij}.\end{array}$$

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

$$\begin{array}{}\text{(45)}& {\displaystyle }{\mathit{\tau }}_{ij}=\mathrm{2}\mathit{\eta }{\dot{\mathit{\epsilon }}}_{ij}.\end{array}$$

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

$$\begin{array}{}\text{(46)}& {\displaystyle }\mathrm{\Phi }={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{e}}^{\mathrm{2}}}{\mathit{\eta }}}.\end{array}$$

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

$$\begin{array}{}\text{(47)}& {\displaystyle }{F}_{\mathrm{d}}^{\mathrm{top}}={\displaystyle \frac{k_{\mathrm{i}}}{H}}{\displaystyle \frac{T_{\mathrm{1}}-T_{\mathrm{0}}}{{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{\mathrm{1}}}},\end{array}$$

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

$$\begin{array}{}\text{(48)}& {\displaystyle }{F}_{\mathrm{d}}^{\mathrm{bot}}={\displaystyle \frac{k_{\mathrm{i}}}{H}}{\displaystyle \frac{T_{nz}-T_{nz-\mathrm{1}}}{\mathrm{1}-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{nz-\mathrm{1}}}}.\end{array}$$

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

$$\begin{array}{}\text{(49)}& {\displaystyle }{F}_{\mathrm{f}}={\mathit{\tau }}_{\mathrm{b}}\cdot {\mathit{u}}_{\mathrm{b}},\end{array}$$

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

$$\begin{array}{}\text{(50)}& {\displaystyle }{F}_{\mathrm{f}}=\mathit{\beta }(u_{\mathrm{b}}^{\mathrm{2}}+v_{\mathrm{b}}^{\mathrm{2}}).\end{array}$$

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:

$$\begin{array}{}\text{(51)}& {\displaystyle }{F}_{\mathrm{g}}+F_{\mathrm{f}}=F_{\mathrm{d}}^{\mathrm{bot}}.\end{array}$$

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:

$$\begin{array}{}\text{(52)}& {\displaystyle }{M}_{\mathrm{b}}={\displaystyle \frac{F_{\mathrm{g}}+F_{\mathrm{f}}-F_{\mathrm{d}}^{\mathrm{bot}}}{{\mathit{\rho }}_{\mathrm{i}}L}},\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{T_k^{n+\mathrm{1}}-T_k^n}{\mathrm{\Delta }t}}={\displaystyle \frac{k_{\mathrm{i}}}{{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}{H}^{\mathrm{2}}}}\left[a_k{T}_{k-\mathrm{1}}^{n+\mathrm{1}}-(a_k+b_k)T_k^{n+\mathrm{1}}+b_k{T}_{k+\mathrm{1}}^{n+\mathrm{1}}\right]\\ \text{(53)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{{\mathrm{\Phi }}_k}{{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}}},\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }-{\mathit{\alpha }}_k{T}_{k-\mathrm{1}}^{n+\mathrm{1}}+(\mathrm{1}+{\mathit{\alpha }}_k+{\mathit{\beta }}_k)T_k^{n+\mathrm{1}}-{\mathit{\beta }}_k{T}_{k+\mathrm{1}}^{n+\mathrm{1}}\\ \text{(54)}& {\displaystyle }& {\displaystyle }=T_k^n+{\displaystyle \frac{{\mathrm{\Phi }}_k\mathrm{\Delta }t}{{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}}},\end{array}$$

where

$$\begin{array}{}\text{(55)}& {\displaystyle }{\mathit{\alpha }}_k={\displaystyle \frac{k_{\mathrm{i}}{a}_k\mathrm{\Delta }t}{{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}{H}^{\mathrm{2}}}},{\mathit{\beta }}_k={\displaystyle \frac{k_{\mathrm{i}}{b}_k\mathrm{\Delta }t}{{\mathit{\rho }}_{\mathrm{i}}{c}_{\mathrm{i}}{H}^{\mathrm{2}}}}.\end{array}$$

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:

$$\begin{array}{}\text{(56)}& {\displaystyle }{F}_{\mathrm{f}}+F_{\mathrm{g}}-{\displaystyle \frac{k_{\mathrm{i}}}{H}}{\displaystyle \frac{T_{nz}^{n+\mathrm{1}}-T_{nz-\mathrm{1}}^{n+\mathrm{1}}}{\mathrm{1}-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{nz-\mathrm{1}}}}=\mathrm{0},\end{array}$$

which can be rearranged to give

$$\begin{array}{}\text{(57)}& {\displaystyle }-T_{nz-\mathrm{1}}^{n+\mathrm{1}}+T_{nz}^{n+\mathrm{1}}={\displaystyle \frac{\left(F_{\mathrm{f}}+F_{\mathrm{g}}\right)H\left(\mathrm{1}-{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_{nz-\mathrm{1}}\right)}{k_{\mathrm{i}}}}.\end{array}$$

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 Mass and tracer transport

The transport equation for ice thickness H can be written as

$$\begin{array}{}\text{(58)}& {\displaystyle \frac{\partial H}{\partial t}}+\mathbf{\nabla }\cdot \left(H\mathit{U}\right)=B,\end{array}$$

where U is the vertically averaged 2-D velocity and B is the total surface mass balance. This equation describes the conservation of ice volume under horizontal transport. With the assumption of uniform density, volume conservation is equivalent to mass conservation. There is a similar conservation equation for the internal energy in each ice layer:

$$\begin{array}{}\text{(59)}& {\displaystyle \frac{\partial \left(hT\right)}{\partial t}}+\mathbf{\nabla }\cdot \left(hT\mathit{u}\right)=\mathrm{0},\end{array}$$

where h is the layer thickness, T is the temperature (located at the layer midpoint), and u is the layer-average horizontal velocity. (Vertical diffusion and internal dissipation are handled separately as described in Sect. 3.2.) If other tracers are present, their transport is described by conservation equations of the same form as Eq. (59). Glissade solves Eqs. (58) and (59) in a coordinated way, one layer at a time. All tracers, including temperature, are advected using the same algorithm.

Glissade has two horizontal transport schemes: a first-order upwind scheme and a more accurate incremental remapping (IR) scheme (Dukowicz and Baumgardner, 2000; Lipscomb and Hunke, 2004). The IR scheme is the default. This scheme was first implemented in the Los Alamos sea ice model, CICE, and has been adapted for CISM. Since the scheme is fairly complex, we give a summary in Sect. 3.3.1 and refer the reader to the earlier publications for details.

After horizontal transport, the mass balance is applied at the top and bottom ice surfaces. The new vertical layers generally do not have the desired spacing in σ coordinates. For this reason, a vertical remapping scheme is applied to transfer ice volume, internal energy, and other conserved quantities between adjacent layers, thus restoring each column to the desired σ coordinates while conserving mass and energy. (This is a common feature of arbitrary Lagrangian–Eulerian (ALE) methods.) Internal energy divided by mass gives the new layer temperatures, and similarly for other tracers.

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,

$$\begin{array}{}\text{(60)}& {\displaystyle \frac{max\left|\mathit{u}\right|\mathrm{\Delta }t}{\mathrm{\Delta }x}}\le \mathrm{1}.\end{array}$$

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):

$$\begin{array}{}\text{(61)}& {\displaystyle }(maxD)\mathrm{\Delta }t\le \mathrm{0.5}\mathrm{\Delta }{x}^{\mathrm{2}},\end{array}$$

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.

3.4 Basal sliding

Glissade assumes a basal friction law of the form Eq. (10). If β were independent of velocity, then Eq. (10) would be a simple linear sliding law. Allowing β to depend on velocity allows more complex and physically realistic sliding laws. The following options are supported in CISM:

• Spatially uniform β, possibly a large value that makes sliding negligible.

• No sliding, enforced by a Dirichlet boundary condition (u_b=0) during finite-element assembly.

• A 2-D β field specified from an external file. Typically, β is chosen at each vertex to optimize the fit to observed surface velocity.

• β is set to a large value where the bed is frozen (T_b<T_pmp) and a lower value where the bed is thawed (T_b=T_pmp).

• Power-law sliding, based on the sliding relation of Weertman (1957). Basal velocity is given by

  $$\begin{array}{}\text{(62)}& {\displaystyle }{\mathit{u}}_{\mathrm{b}}=k{\mathit{\tau }}_{\mathrm{b}}^p{N}^{-q},\end{array}$$

  where N is the effective pressure (the difference between the ice overburden pressure P₀=ρ_igH and the basal water pressure p_w) and k is an empirical constant. This can be rearranged to give

  $$\begin{array}{}\text{(63)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{b}}=k^{-\mathrm{1}/p}{N}^{q/p}{\mathit{u}}_{\mathrm{b}}^{\mathrm{1}/p}.\end{array}$$

• Plastic sliding, in which the bed deforms with a specified till yield stress.

• Coulomb friction as described by Schoof (2005), with notation and default values from Pimentel et al. (2010). The form of the sliding law is

  $$\begin{array}{}\text{(64)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{b}}=C{\left|\mathit{u}\right|}^{\frac{\mathrm{1}}{n}-\mathrm{1}}\mathit{u}{\left({\displaystyle \frac{N^n}{\mathit{\kappa }\left|\mathit{u}\right|+N^n}}\right)}^{\frac{\mathrm{1}}{n}},\end{array}$$

  where C is a constant, $\mathit{\kappa }=m_{max}/\left({\mathit{\lambda }}_{max}{A}_{\mathrm{b}}\right)$, m_max is the maximum bed obstacle slope, λ_max is the wavelength of bedrock bumps, and A_b is a basal flow-law parameter. Equation (64) has two asymptotic behaviors. In the interior, where the ice is thick and slow-moving, κ|u|≪Nⁿ and the basal friction is independent of N:

  $$\begin{array}{}\text{(65)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{b}}\approx C{\left|\mathit{u}\right|}^{\frac{\mathrm{1}}{n}-\mathrm{1}}\mathit{u}.\end{array}$$

  The Coulomb-friction limit, where the ice is thin and fast-moving, we have κ|u|≫Nⁿ, giving

  $$\begin{array}{}\text{(66)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{b}}\approx {\displaystyle \frac{C}{{\mathit{\kappa }}^{\frac{\mathrm{1}}{n}}}}N{\displaystyle \frac{\mathit{u}}{\left|\mathit{u}\right|}}.\end{array}$$

• Pseudo-plastic sliding, as described by Schoof and Hindmarsh (2010) and Aschwanden et al. (2013) and implemented in PISM. The basal friction law is

  $$\begin{array}{}\text{(67)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{b}}=-{\mathit{\tau }}_{\mathrm{c}}{\displaystyle \frac{{\mathit{u}}_{\mathrm{b}}}{{\left|{\mathit{u}}_{\mathrm{b}}\right|}^{\mathrm{1}-q}{u}_{\mathrm{0}}^q}},\end{array}$$

  where τ_c is the yield stress, q is a pseudo-plastic exponent, and u₀ is a threshold speed. This law incorporates linear (q=1), plastic (q=0), and power-law ($\mathrm{0}<q<\mathrm{1}$) behavior in a single expression. The yield stress is given by τ_c=tan (ϕ)N, where ϕ is a friction angle that can vary with bed elevation, resulting in a lower yield stress at lower elevations. This scheme has been shown to give a realistic velocity field for most of the Greenland ice sheet with tuning of just a few parameters, instead of adjusting a basal friction parameter in every grid cell (Aschwanden et al., 2016).

The power-law, Coulomb-friction and pseudo-plastic sliding laws give lower basal friction and faster sliding as the effective pressure N decreases from the overburden pressure. CISM offers several options for computing N at the bed:

• N=ρ_igH, the full overburden pressure. That is, the water pressure p_w=0 at the bed.

• N is reduced as the basal water depth increases, reaching a small fraction of overburden pressure (typically δ=0.02) when the water depth reaches a prescribed threshold.

• Following Leguy et al. (2014), N is reduced where the bed is below sea level, to account for partial or full connectivity of the basal water system to the ocean. The effective pressure is given by

  $$\begin{array}{}\text{(68)}& {\displaystyle }N={\mathit{\rho }}_{\mathrm{i}}gH{\left(\mathrm{1}-{\displaystyle \frac{H_{\mathrm{f}}}{H}}\right)}^p,\end{array}$$

  where $H_{\mathrm{f}}=max(\mathrm{0},-{\mathit{\rho }}_{\mathrm{o}}b/{\mathit{\rho }}_{\mathrm{i}})$ is the flotation thickness, ρ_o is seawater density, b is the bed elevation (negative below sea level), and $\mathrm{0}\le p\le \mathrm{1}$. For p=0, there is no connectivity and N is the full overburden pressure. For p=1, there is full connectivity, and the basal water pressure is equal to the ocean pressure at the same depth.

Glissade does not yet support sophisticated subglacial hydrology. (The model of Hoffman and Price (2014) was implemented in serial in an earlier version of CISM but is not currently supported.) The following (non-conserving) options are available for handling basal meltwater:

• All basal water immediately drains.

• Basal water depth is set to a constant everywhere, to force T=T_pmp.

• Basal water depth is computed using a local till model, as described by Bueler and van Pelt (2015). In this model, water depth W evolves according to

  $$\begin{array}{}\text{(69)}& {\displaystyle \frac{\partial W}{\partial t}}={\displaystyle \frac{-B_{\mathrm{b}}}{{\mathit{\rho }}_{\mathrm{w}}}}-C_{\mathrm{d}},\end{array}$$

  where W is the water depth (capped at W_max), B_b is the basal mass balance (negative for melting), ρ_w is the density of water, and C_d is a fixed drainage rate. The effective pressure N is related to water depth by

  $$\begin{array}{}\text{(70)}& {\displaystyle }N=min\left[P_{\mathrm{0},}{N}_{\mathrm{0}}{\left({\displaystyle \frac{\mathit{\delta }{P}_{\mathrm{0}}}{N_{\mathrm{0}}}}\right)}^{W/W_{\max }}{\mathrm{10}}^{(e_{\mathrm{0}}/C_{\mathrm{c}})(\mathrm{1}-W/W_{\max })}\right],\end{array}$$

  where e₀ is the void ratio at reference effective pressure N₀, C_c is the till compressibility, and δ is a scalar between 0 and 1. Default values of these terms are taken from Bueler and van Pelt (2015). The effect of Eq. (70) is to drive N from P₀ down to δP₀ as W increases from 0 to W_max.

3.5 Iceberg calving

Glissade supports several schemes for calving ice at marine margins. Two of these are very simple: (1) calve all floating ice and (2) calve ice where the bedrock (either the current bedrock or the bedrock toward which a viscous asthenosphere is relaxing) lies below a prescribed elevation. In addition, Glissade includes several new calving schemes: mask-based calving, thickness-based calving, and eigencalving.

With mask-based calving, any floating ice that forms in cells defined by a calving mask is assumed to calve instantly. By default, the mask is based on the initial ice extent. Floating ice calves in all cells that are initially ice-free ocean, and thus the calving front cannot advance (but it can retreat). Alternatively, users can define a calving mask that is read at initialization.

Thickness-based calving is designed to remove floating ice that is thinner than a user-defined thickness, $H_{\mathrm{c}}^{\text{min}}$. As discussed by Albrecht et al. (2011), accurate thickness-based calving requires a subgrid-scale parameterization of the calving front. Suppose, for example, that $H_{\mathrm{c}}^{\text{min}}=\mathrm{100}$ m, and an ice shelf with a calving front thicker than 100 m is advancing. During the next time step, ice is transported to grid cells just downstream of the initial calving front. Typically, these cells have $H<H_{\mathrm{c}}^{\min }$. If this thin ice immediately calves, the calving front cannot advance.

Glissade avoids this problem using an approach similar to that of Albrecht et al. (2011). Floating cells adjacent to ice-free ocean are identified as calving-front (CF) cells. For each CF cell, an effective thickness H_eff is determined as the minimum thickness of adjacent ice-filled cells not at the CF. CF cells with H<H_eff are deemed to be partially filled. For example, a cell with H=50 m and H_eff=100 m is considered to be half-filled with 100 m ice. It is dynamically inactive and thus cannot transport ice downstream. It can thicken, however, as ice is transported from active cells upstream. Once the ice in this cell thickens to H≥H_eff, it becomes dynamically active. The downstream cells, previously ice-free, then become inactive CF cells and can thicken in turn. In this way, the calving front can advance. Similarly, the calving front can retreat when an inactive CF cell becomes ice-free; its upstream neighbor, formerly an active cell in the shelf interior, becomes an inactive CF cell.

Thickness-based calving is applied not to CF cells with $H<H_{\mathrm{c}}^{\min }$ but rather to cells with $H_{\text{eff}}<H_{\mathrm{c}}^{\min }$. In other words, the CF thickness is inferred from the active cells just interior to the CF. Where $H_{\text{eff}}<H_{\mathrm{c}}^{\min }$, the effective rate of thickness change is given by

$$\begin{array}{}\text{(71)}& {\displaystyle \frac{\mathrm{d}{H}_{\text{eff}}}{\mathrm{d}t}}=-{\displaystyle \frac{\left(H_{\mathrm{c}}^{min}-H_{\text{eff}}\right)}{t_{\mathrm{c}}}},\end{array}$$

where t_c is a calving timescale. This effective rate is converted to dH∕dt (i.e., the rate of ice volume change per unit cell area) as

$$\begin{array}{}\text{(72)}& {\displaystyle \frac{\mathrm{d}H}{\mathrm{d}t}}=min\left(H/H_{\text{eff}},\mathrm{1}\right){\displaystyle \frac{\mathrm{d}{H}_{\text{eff}}}{\mathrm{d}t}}.\end{array}$$

Glissade's eigencalving scheme is related to the method described by Levermann et al. (2012) but differs in the details. In the method of Levermann et al. (2012), the lateral calving rate is proportional to the product of the two eigenvalues of the horizontal strain rate tensor, provided that both eigenvalues are positive. This scheme was tested in CISM but found to give noisy, erratic results. Instead, we compute the lateral calving rate c as

$$\begin{array}{}\text{(73)}& {\displaystyle }c=k_{\mathit{\tau }}{\mathit{\tau }}_{\mathrm{ec}},\end{array}$$

where k_τ is an empirical constant with units of m yr⁻¹ Pa⁻¹, and τ_ec is an effective calving stress defined by

$$\begin{array}{}\text{(74)}& {\displaystyle }{\mathit{\tau }}_{\mathrm{ec}}^{\mathrm{2}}=max({\mathit{\tau }}_{\mathrm{1}},\mathrm{0}{)}^{\mathrm{2}}+w_{\mathrm{2}}max({\mathit{\tau }}_{\mathrm{2}},\mathrm{0}{)}^{\mathrm{2}}.\end{array}$$

Here, τ₁ and τ₂ are the eigenvalues of the 2-D horizontal deviatoric stress tensor, and w₂ is an empirical constant. The stresses τ₁ and τ₂ are positive when the ice is in tension; τ₁ corresponds to along-flow tension and τ₂ to across-flow tension. These stresses are computed for dynamically active cells and then are extrapolated to inactive neighbors. When w₂>1, across-flow tension drives calving more effectively than does along-flow tension, as is the case for the calving law of Levermann et al. (2012). Equation (74) is analogous to Eq. (6) in Morlighem et al. (2016) (in a calving law for Greenland's Store Gletscher) but with stress replacing strain rate and with a weighting term added.

Using Eqs. (73) and (74), we can compute a lateral calving rate c for each CF cell. The lateral calving rate is converted to a rate of volume change using

$$\begin{array}{}\text{(75)}& {\displaystyle \frac{\mathrm{d}H}{\mathrm{d}t}}=-{\displaystyle \frac{H_{\text{eff}}c}{\sqrt{\mathrm{\Delta }x\mathrm{\Delta }y}}},\end{array}$$

where Δx and Δy are the grid cell dimensions. Typically, CISM is run with Δx=Δy, but this is not a model requirement.

Eigencalving, when applied on its own, can allow very thin ice to persist at the calving front where stresses are small. For this reason, the eigencalving algorithm is automatically followed by thickness-based calving as described above, with $H_{\mathrm{c}}^{\text{min}}$ set to a value large enough to remove thin ice, where present, but not so large as to be the dominant driver of calving.

An additional option related to calving is iceberg removal. An iceberg is defined as a connected region where every cell is floating and has no connection to grounded ice. Icebergs can be present in initial data sets and also can arise during simulations, depending on the calving scheme. Finding the velocity for such a region is an ill-posed problem, so it is best to remove icebergs as soon as they occur. This is done by first marking all the grounded ice cells, and then using a parallel flood-fill algorithm to mark every floating cell that is connected to a grounded cell. (That is, one could travel from a grounded cell to a given floating cell along a path at least one cell wide, without leaving the ice sheet.) Any unmarked floating cells then calve immediately.

Another option is to limit the thickness of ice cliffs, defined as grounded marine-based cells adjacent to ice-free ocean. Bassis and Walker (2012) pointed out that there is an upper thickness limit for marine-terminating cliffs of outlet glaciers. If the ice cliff sits more than ∼100 m out of the water, the longitudinal stresses exceed the yield strength, and the ice will fail. The maximum stable thickness at the terminus is given by

$$\begin{array}{}\text{(76)}& {\displaystyle }{H}_{max}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{c}}}{{\mathit{\rho }}_{\mathrm{i}}g}}+\sqrt{{\left({\displaystyle \frac{{\mathit{\tau }}_{\mathrm{c}}}{{\mathit{\rho }}_{\mathrm{i}}g}}\right)}^{\mathrm{2}}+{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{o}}}{{\mathit{\rho }}_{\mathrm{i}}}}{d}^{\mathrm{2}}},\end{array}$$

where τ_c∼1 MPa is the depth-averaged shear strength, d is the ocean depth, and ρ_i and ρ_o are densities of ice and seawater. When cliff limiting is enabled in CISM, marine-grounded cells adjacent to ice-free ocean are limited to H≤H_max, with any excess thickness added to the calving flux. This thinning mechanism does not trigger the rapid ice sheet retreat seen by Pollard et al. (2015) in Antarctic simulations that combined cliff failure with hydrofracture. CISM2.1 does not simulate hydrofracture or lateral cliff retreat.

3.6 Sub-shelf melting

By default, melting beneath ice shelves is set to zero. The following simple options may be enabled for simulations with marine ice:

• Sub-shelf melting is set to a constant value for all floating ice.

• CISM reads in a 2-D field of basal melt rates and applies the prescribed rates to floating ice.

• Sub-shelf melting is prescribed as for the Marine Ice Sheet Model Intercomparison Project (MISMIP+) Ice1 experiments described by Asay-Davis et al. (2016). In this parameterization, the melt rate is set to zero above a reference elevation z₀, then increases linearly as a function of depth. In shallow sub-shelf cavities, the melt rate is reduced to zero over a characteristic length scale H₀. The depth dependence is motivated by observations of temperature increasing with depth in polar regions but is not necessarily a realistic treatment of melting near grounding lines.

• The sub-shelf melt rate is a piecewise linear function of depth. Above a reference elevation z₀ there is freezing, and below z₀ there is melting. Above z₀, the freezing rate increases linearly from 0 to $B_{\text{frz}}^{\text{max}}$ at $z_{\text{frz}}^{\text{max}}$; above $z_{\text{frz}}^{\text{max}}$, it is capped at $B_{\text{frz}}^{\text{max}}$. Similarly, below z₀, the melt rate increases linearly from 0 to $B_{\text{mlt}}^{\text{max}}$ at $z_{\text{mlt}}^{\text{max}}$; below $z_{\text{mlt}}^{\text{max}}$, it is capped at $B_{\text{mlt}}^{\text{max}}$. In shallow cavities, the melt rate is reduced to zero over a scale H₀ as for MISMIP+. Like the MISMIP+ parameterization, this is an ad hoc scheme that should be used with caution. More realistic treatments of sub-shelf melting are being actively developed in the ocean and ice modeling communities (Asay-Davis et al., 2017).

Sub-shelf melting is applied only to cells that are floating based on the criterion $b<-({\mathit{\rho }}_{\mathrm{i}}/{\mathit{\rho }}_{\mathrm{o}})H$, where b is bed elevation (negative below sea level) and H is ice thickness. In partly filled cells at the calving front, basal melt is applied to the effective thickness H_eff rather than the grid cell mean thickness H (see Sect. 3.5). For example, a melt rate of 10 m yr⁻¹ applied to a cell that is 10 % full would reduce H at a rate of only 1 m yr⁻¹. Basal melting for grounded ice, including marine-based ice, is described in Sect. 3.2.3.

3.7 Isostasy

Isostatic adjustment is handled as in Rutt et al. (2009), with several options for the lithosphere and underlying asthenosphere. The lithosphere can be described as either local (i.e., floating on the asthenosphere) or elastic (taking flexural rigidity into account). The asthenosphere is either fluid (reaching isostatic equilibrium instantaneously) or relaxing (responding to the load on an exponential timescale of several thousand years). The default, when isostasy is active, is an elastic lithosphere and relaxing asthenosphere.

In the elastic lithosphere calculation, data from many grid cells must be summed to compute the load in each location. For parallel runs, this is done by gathering data to one processor to compute the load and then distributing the result to the local processors. This calculation scales poorly, but it is done infrequently (every ∼100 model years) and therefore has a minimal computational cost at grid resolutions of ∼5 km.

4 Model results: standard test cases

Test cases for CISM include (1) problems with analytic solutions, (2) standard experiments without analytic solutions but for which community benchmarks are available, and (3) some CISM-specific experiments used for code development and testing. These tests are organized in directories, each of which includes a README file with instructions on running the test. Most tests have Python scripts that are used to set up the initial conditions and run the model, and some tests have an additional Python script to analyze and plot the output. Each test has a default configuration file, whose settings can be adjusted manually to make changes (e.g., swapping velocity solvers). Several tests are described below, with plots generated by Python scripts included in the code release.

4.1 Halfar dome test

The Halfar test case describes the time evolution of a parabolic dome of ice (Halfar, 1983; Bueler et al., 2005). For a flat-bedded SIA problem without accumulation, this case has an analytic solution for the time-varying ice thickness. The SIA ice evolution equation can be written as

$$\begin{array}{}\text{(77)}& {\displaystyle \frac{\partial H}{\partial t}}=\mathbf{\nabla }\cdot \left(\mathrm{\Gamma }{H}^{n+\mathrm{2}}\right|\mathbf{\nabla }H{|}^{n-\mathrm{1}}\mathbf{\nabla }H),\end{array}$$

where n is the exponent in the Glen flow law, commonly taken as 3, and $\mathrm{\Gamma }=\mathrm{2}A(\mathit{\rho }g{)}^n/(n+\mathrm{2})$ is a positive constant. For n=3, the time-dependent solution is

$$\begin{array}{}\text{(78)}& {\displaystyle }H(t,r)=H_{\mathrm{0}}{\left({\displaystyle \frac{t_{\mathrm{0}}}{t}}\right)}^{\frac{\mathrm{1}}{\mathrm{9}}}{\left[\mathrm{1}-{\left({\left({\displaystyle \frac{t_{\mathrm{0}}}{t}}\right)}^{\frac{\mathrm{1}}{\mathrm{18}}}{\displaystyle \frac{r}{R_{\mathrm{0}}}}\right)}^{\frac{\mathrm{4}}{\mathrm{3}}}\right]}^{\frac{\mathrm{3}}{\mathrm{7}}},\end{array}$$

where

$$\begin{array}{}\text{(79)}& {\displaystyle }{t}_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{18}\mathrm{\Gamma }}}{\left({\displaystyle \frac{\mathrm{7}}{\mathrm{4}}}\right)}^{\mathrm{3}}{\displaystyle \frac{R_{\mathrm{0}}^{\mathrm{4}}}{H_{\mathrm{0}}^{\mathrm{7}}}},\end{array}$$

and H₀ and R₀ are the central dome height and dome radius, respectively, at time t=t₀. As the dome evolves, the ice margin advances and the thickness decreases. The test can be run using either Glide or Glissade solvers.

Figure 1a–c show Halfar dome results for the Glide solver with $H_{\mathrm{0}}=\mathrm{500}\sqrt{\mathrm{2}}\approx \mathrm{700}$ m, $R_{\mathrm{0}}=\mathrm{15}\sqrt{\mathrm{2}}\approx \mathrm{21}$ km, a grid resolution of 2 km, and a time step of 5 years. The three panels show the modeled thickness, analytic thickness, and thickness difference, respectively, at t=200 years. Differences are largest near the ice margin. The rms error in the solution is 6.43 m, with a maximum absolute error of 32.8 m.

Figure 1 Results from the Halfar dome test using SIA solvers from Glide (a, b, c) and Glissade (d, e, f). (a, d) Modeled ice thickness. (b, e) Analytic ice thickness. (c, f) Difference between modeled and analytic thickness.

The bottom panels of Fig. 1 show Halfar results for the Glissade SIA solver described in Sect. 3.1.2, with the same grid and initial conditions and a time step of 1 year (Glissade's explicit transport solver requires a shorter time step than Glide for numerical stability). The errors are larger than for Glide; the rms error is 9.06 m and the maximum error is 34.2 m. Since Glissade's SIA velocity solutions are very similar to Glide's for a given ice sheet geometry, the larger errors can be attributed to the explicit transport solver (Sect. 3.3), which is less accurate for this problem than Glide's implicit diffusion-based solver.

4.2 ISMIP-HOM tests

ISMIP-HOM, described by Pattyn et al. (2008), consists of six experiments, labeled A through F, for higher-order velocity solvers. CISM supports all six experiments, and here we show results for experiments A and C. These two tests are particularly useful for benchmarking higher-order models, since they gauge the accuracy of simulated 3-D flow over a bed with large- and small-scale variations in basal topography and friction.

Experiment A is a test for ice flow over a bumpy bed. The domain is doubly periodic, and the bed topography is sinusoidal in both the x and y directions, with an amplitude of 500 m and wavelength L=5, 10, 20, 40, 80, or 160 km. The mean ice thickness is 1000 m, and the surface slopes smoothly downward from left to right. The velocity field is diagnosed from this geometry given a no-slip basal boundary condition. Figure 2 shows the surface ice speed as a function of x along the bump at $y=L/\mathrm{4}$, computed using Glissade's BP solver. For L=10 km or more, the solution agrees well with the mean from the Stokes models in Pattyn et al. (2008). At L=5 km, the amplitude is too large (as is typical of BP models), but the magnitude is close to the Stokes solution.

Figure 2 Results for ISMIP-HOM experiment A (ice flow over a bumpy bed). Each plot shows the surface ice speed across the bump at $y=L/\mathrm{4}$ for a different length scale L. The solid black line shows output from Glissade's BP velocity solver; the dotted red and blue lines show output from full Stokes and first-order (i.e., higher-order) models, respectively, in Pattyn et al. (2008), and the red and blue shaded regions show the corresponding standard deviations.

Figure 3 shows test A results using Glissade's DIVA solver (cf. Fig. 1 in Goldberg, 2011, which shows similar results from a depth-integrated flowline solver). The DIVA solution is much less accurate than the BP solution; the modeled ice speed is greater than for the Stokes solution, with the relative error increasing as L decreases. These errors are expected, since a bumpy, frozen bed implies that the horizontal velocity gradients are not well approximated (as DIVA assumes) by their vertical averages. The DIVA solution is more accurate, however, than the L1L2 scheme implemented by Perego et al. (2012). The latter scheme does a 2-D matrix solve for the basal velocity and then integrates upward, reducing to the very inaccurate SIA when the basal velocity is zero. DIVA, by solving for the 2-D mean velocity instead of the basal velocity, is more accurate than the SIA.

Figure 3 Same as Fig. 2, except that the solid black line shows output from Glissade's DIVA velocity solver. Vertical scales are adjusted relative to Fig. 2 to show the full range of the model solution.

Figure 4 Results for ISMIP-HOM experiment C (ice flow over a bed with variable basal friction). Each plot shows the surface ice speed across the bump at $y=L/\mathrm{4}$ for a different length scale L. The solid black line shows output from Glissade's BP velocity solver; the dotted red and blue lines show output from Stokes and first-order models, respectively, in Pattyn et al. (2008), and the red and blue shaded regions show the corresponding standard deviations.

Experiment C is a test for flow over a bed with variable friction. The ice surface and mean thickness are as in experiment A, but the bed is flat, with sinusoidal variations in the basal friction parameter β. Figure 4 shows the surface speed at $y=L/\mathrm{4}$ for each of six wavelengths using the BP solver, and Fig. 5 shows the surface speed using the DIVA solver (cf. Fig. 2 from the depth-integrated flowline solver in Goldberg, 2011). For both BP and DIVA, the CISM results are similar to the Stokes results at all wavelengths. There is a modest difference between BP and DIVA at L=5 km, where the DIVA velocity has a small variation across the domain, whereas the BP and Stokes velocities are nearly uniform.

Figure 5 Same as Fig. 4, except that the solid black line shows output from Glissade's DIVA velocity solver.

For both experiments A and C, the Glissade BP results are nearly identical to those shown by Tezaur et al. (2015), who used similar finite-element methods to assemble and solve the BP equations.

4.3 Stream tests

CISM's stream tests simulate flow over an idealized ice stream underlain by subglacial till with a specified yield stress distribution. For the two yield stress distributions specified in this test case, analytical solutions are available from Raymond (2000) and Schoof (2006). For the Raymond test case, the yield stress is given within the ice stream by a uniform value of 5.2 kPa (below the gravitational driving stress of about 10 kPa), and outside the ice stream by a uniform value of 70 kPa (much larger than the driving stress). That is, the yield stress distribution is approximated by a step function. For the Schoof test case, the yield stress across the ice stream varies continuously between about 45 kPa and 0. In both cases, the basal properties and resulting velocity solutions vary in the across-flow direction only and are symmetric about the ice stream centerline.

Figures 6 and 7 compare model results to the analytic Raymond and Schoof solutions, respectively. CISM is run using the BP velocity solver at a grid resolution of 2 km. (DIVA results, not shown, are nearly identical.) The top panels show across-flow velocity profiles, and the bottom panels show the prescribed yield stress and driving stress. In both cases, there is excellent agreement between the modeled and analytic solutions. For the Raymond test, this excellent agreement requires a modified assembly procedure for basal friction terms, as described in Sect. A2, to resolve the step change in yield stress. With the standard finite-element procedure, there is some smoothing of the friction parameter β over neighboring vertices, and the results at a given grid resolution are less accurate. For the Schoof test with a smooth transition in yield stress, either assembly procedure works well.

Figure 6 Results from the Raymond stream test using Glissade's BP velocity solver. (a) Across-flow surface and basal velocity (m yr⁻¹) at x=15 km compared to the analytic solution. At most points, the analytic and simulated solutions are indistinguishable. (b) Prescribed yield stress and gravitational driving stress (Pa).

Figure 7 Same as Fig. 6 but for the Schoof stream test with a smooth transition in yield stress. The analytic and simulated solutions are virtually indistinguishable.

4.4 Shelf tests

CISM includes three shelf tests:

• The first is a confined-shelf test based on tests 3 and 4 from the Ice Shelf Model Intercomparison exercise (Rommelaere, 1996). This test simulates the flow of a 500 m thick ice shelf within a rectangular embayment that is confined on three sides. Figure 8 shows the 2-D ice speed using Glissade's SSA solver (Sect. 3.1.3) on a 5 km mesh. As expected, the solution is nearly identical to those obtained with the DIVA and BP solvers (not shown). No benchmark is available, but the Glissade solutions are similar to the first-order (BP) solutions found by Perego et al. (2012) and Tezaur et al. (2015) for the same test.

• The second is a circular-shelf test that simulates the radially symmetric flow of a 1000 m thick ice shelf that is pinned at one point in the center. The SSA solution (not shown) is similar but not identical to the DIVA and BP solutions, since the latter allow vertical shear above the pinning point.

• The third test simulates the flow of the Ross Ice Shelf in Antarctica under idealized conditions (e.g., a spatially uniform flow-rate factor), as described by MacAyeal et al. (1996). Figure 9, which is based on Figs. 1 and 2 of MacAyeal et al. (1996), shows the ice speed computed by Glissade's SSA solver on the 6.8 km finite-difference grid used in that paper, compared to observed ice speeds. The model velocities agree well with the published velocities from the finite-element model in Fig. 3a of MacAyeal et al. (1996) and, like the published values, tend to be biased fast compared to observations. As noted by MacAyeal et al. (1996), the model–data agreement could be improved by allowing the flow-rate factor to vary spatially. Results from the DIVA and BP solvers (not shown) are nearly identical.

Figure 8 Ice speed (m yr⁻¹) for the confined-shelf test (based on tests 3 and 4 from the Ice Shelf Model Intercomparison exercise of Rommelaere, 1996) using Glissade's SSA velocity solver. The ice shelf is 500 m thick. The rectangular domain is open to the south and confined on the other three sides.

4.5 Dome test

CISM's dome test has a simple configuration with a parabolic, radially symmetric dome on a flat bed. It is a general-purpose time-dependent problem that can be used to test not just the velocity solver but also the transport and temperature solvers and various physics options. There is no analytic solution, but the test has proven useful for day-to-day testing.

4.6 Build-and-test structure

To facilitate testing, CISM includes a build-and-test structure (BATS) that can automatically build the model and then run a set of regression tests, including the tests discussed above. BATS allows users and developers to quickly generate a set of regression tests for use with the Land Ice Verification and Validation toolkit (LIVVkit; Kennedy et al., 2017). LIVVkit is designed to provide robust and automated numerical verification, software verification, performance validation, and physical validation analyses for ice sheet models. Instructions on using BATS and LIVV with CISM can be found on the LIVV website (https://livvkit.github.io/Docs/, last access: 20 January 2019).

5 Model results: Greenland ice sheet simulations

The CISM v2.0 release included the tests described in Sect. 4 but did not support robust multimillennial simulations of whole ice sheets. Recent work, included in CISM v2.1, has made the model more efficient, reliable, and realistic for Greenland ice sheet simulations. These improvements support the use of CISM in CESM2.0 for century- to millennial-scale Greenland simulations under paleoclimate (e.g., Pliocene and Eemian), present-day, and future climate conditions.

Coupled ice sheet–climate simulations can be problematic, because ice sheets require 10⁴ years or longer to equilibrate with a given climate (Vizcaino, 2014), and because climate is never truly constant on these timescales. If the initial ice sheet conditions in a coupled simulation are not consistent with the climate, then the transient adjustment can swamp any climate change signal. One way to minimize the adjustment is to adjust selected model parameters (e.g., basal traction coefficients or bed topography at each grid point) to satisfy an optimization problem, reducing the mismatch between model variables and observations. In this way, one can generate ice sheet thickness and velocity fields consistent with the SMB from a climate model (Perego et al., 2014). The risk of this approach is that the model can be overtuned to present-day conditions that might not hold on long timescales or in different climates. An alternate approach is to spin up the ice sheet to equilibrium under forcing from a climate model: e.g., steady pre-industrial forcing or forcing spliced together from two or more climate time slices (Fyke et al., 2014). Because of climate model biases, however, the spun-up ice thickness and velocity can differ substantially from modern observations.

Here, we take the second approach, spinning up the Greenland ice sheet with surface forcing from the regional climate model RACMO2.3p1 statistically downscaled to 1 km resolution (Noël et al., 2016). Since RACMO2 is run at high resolution, is constrained by reanalysis at model boundaries, and is well validated against observations, its SMB is more realistic than that of a global climate model. RACMO2 provides the SMB only for the region included within its ice sheet mask; outside this region, we prescribe a negative SMB. The ice sheet is initialized with present-day extent and thickness, and then is spun up for 50 000 years (long enough to reach quasi-equilibrium) on a 4 km grid (the standard grid for CESM ice sheet simulations).

We analyze two experiments. The first experiment uses a “no-float” calving scheme, in which floating ice is assumed to calve immediately. A similar configuration was used to initialize CISM for the initMIP-Greenland experiment (Goelzer et al., 2018). The no-float assumption simplifies model setup and analysis, while still simulating ice flow for the vast majority of the Greenland ice sheet. The second experiment uses the eigencalving scheme described in Sect. 3.5, along with the depth-dependent sub-shelf melting scheme of Sect. 3.6. This simulation tests the model's ability to generate robust, stable ice shelves that at least roughly resemble Greenland's present-day ice shelves, including the floating shelf of Petermann Glacier in northern Greenland and two floating termini of the Northeast Greenland Ice Stream (NEGIS): Nioghalvfjerdsbræ (79 North Glacier) and Zachariae Isstrom.

Although the model has been configured to give results that are reasonably realistic, these spin-ups should be viewed primarily as a demonstration of CISM capabilities, not as scientifically validated simulations. For example, model tuning has likely compensated for biases in the forcing data. Generally speaking, model parameters should always be tested and reviewed depending on the science application.

5.1 Simulation without ice shelves

The simulation without floating ice shelves is configured as follows:

• The model is initialized with present-day ice sheet geometry. The ice thickness and bed topography are based on the mass-conserving-bed method of Morlighem et al. (2011, 2014).

• The SMB forcing over the Greenland ice sheet is a 1958–2015 climatology from RACMO2.3p1 (Noël et al., 2016). For ice-free regions where RACMO2 does not compute an SMB, the SMB is arbitrarily set to −4 m yr⁻¹, thus limiting deviations from present-day ice extent. The surface temperature is from the 20th century RACMO2 climatology of Ettema et al. (2009).

• The basal geothermal heat flux is prescribed from Shapiro and Ritzwoller (2004).

• Where the SMB is negative, the initial temperature profile in each column is linear, with $T=\min (T_{\mathrm{air}},\mathrm{0})$ at the surface and $T=T_{\mathrm{pmp}}-\mathrm{5}{}^{\circ }$ at the bed. Where the SMB is positive, the temperature is initialized to an analytic profile based on a balance between vertical conduction and cold advection (Cuffey and Paterson, 2010, Sect. 9.5.1). This profile is relatively cold in the upper part of the ice sheet.

• The 3-D velocity field is computed at 11 vertical levels using the DIVA solver described in Sect. 3.1.4. Thickness and tracer evolution are given by the incremental remapping scheme of Sect. 3.3, with a time step dt=0.2 years.

• The basal friction coefficient β is computed at each vertex using a pseudo-plastic sliding law as described in Sect. 3.4. This scheme has several parameters that can be tuned to reduce effective pressure and increase sliding in warm, low-elevation regions. Parameter values, mostly following Aschwanden et al. (2016), are q=0.5 and u₀=100 m yr⁻¹, with ϕ decreasing from 40 to 5^∘ as the bed elevation decreases from 700 to −700 m.

• Effective pressure is a function of basal water depth as computed by the local till model described in Sect. 3.4 with till parameters from Bueler and van Pelt (2015). In particular, W_max=2 m, C_d=1 mm yr⁻¹, and δ=0.02.

• Floating ice calves immediately, and marine cliff heights are limited as described in Sect. 3.5.

• Some of the model physics is constrained to make the simulation more robust. The surface gradients $\partial s/\partial x$ and $\partial s/\partial y$ in the gravitational driving stress are limited to a magnitude of 0.10 to prevent unrealistically large ice speeds in coastal regions with steep topography, and β is held to a minimum of 100 Pa m⁻¹ yr for grounded ice to prevent very fast sliding.

After 50 000 model years, the ice sheet extent and volume have equilibrated to 1.63×10⁶ km² and 2.97×10⁶ km³, respectively, in close agreement with values of 1.67×10⁶ km² and 2.95×10⁶ km³ in the observational data set (Morlighem et al., 2014). The close agreement between observed and modeled extent and volume depends, in part, on the negative SMB applied outside the RACMO2 ice sheet mask to inhibit ice advance. We checked for equilibrium by comparing thickness snapshots at 500-year intervals near the end of the run. Thickness variations over the last 1000 years (not shown) are less than 50 m everywhere on the ice sheet.

Figure 9 Results of the Ross Ice Shelf test described by MacAyeal et al. (1996), computed using Glissade's SSA solver. Panel (a), based on Fig. 1 of MacAyeal et al. (1996), shows the vertically uniform ice speed (m yr⁻¹). The circles indicate locations where ice speed was measured by the Ross Ice Shelf Geophysical and Glaciological Survey (RIGGS); measured speeds are shaded with the same color scale as the observations. Axis coordinates match the coordinates in MacAyeal et al. (1996). Panel (b), based on Fig. 2 of MacAyeal et al. (1996), compares simulated to measured ice speeds at the RIGGS locations.

Figure 10 shows the surface ice speed at the end of the run, compared to observations. In most of the ice sheet, the simulated ice flow is in good agreement with observations. The model captures slow-flowing regions in the interior where the bed is frozen, as well as fast-flowing outlet glaciers along the margins. A major shortcoming of the simulation is its failure to capture NEGIS, which in reality extends far into the northeast interior but is simulated to be slower, shorter, and more diffuse. Aschwanden et al. (2016), who ran PISM using the pseudo-plastic basal sliding scheme adopted by CISM, obtained similar results. While capturing complex outlet glacier flow, both models have a poor representation of NEGIS, possibly due to a missing geothermal heat source or simplified subglacial hydrology.

Figure 10 Surface ice speed (m yr⁻¹, log scale) for the Greenland ice sheet. (a) Observed speed (Joughin et al., 2010). (b) Simulated speed at the end of a 50 kyr Greenland spin-up without floating ice shelves. The gray background shows ice-free regions. Figures 10–17 were created using the NCAR Command Language (2017).

Figure 11 Difference (m) between (1) simulated ice thickness at the end of a 50 kyr Greenland spin-up without ice shelves and (2) observed thickness from Morlighem et al. (2014). The observed thickness is the model initial condition. The gray background shows regions that are ice-free at both the beginning and end of the simulation.

Figure 11 shows the difference between the final ice thickness and the observed initial thickness. In most of the interior, the thickness errors are ∼100 m or less, with the main exception being the upper part of NEGIS, where slow flow correlates with a positive thickness bias of ∼200 m. Errors are larger along the margins, with positive biases of up to ∼500 m in the southeast, where accumulation is large, and negative biases of up to ∼500 m in the north and northwest, where accumulation is lower and there is significant ablation at the margins. Low thicknesses along the coast could result from excessive marginal ablation in the SMB data set, insufficient basal friction, or a combination. High thickness in the southeast could be attributed to positive accumulation biases in the data set or too-slow flow in outlet glaciers, among other factors.

Figure 12 illustrates the basal water state. The ice sheet is shaded to indicate three regions: a frozen central region (extending to parts of the northern and eastern margins) with basal water depth W=0; thawed coastal regions, especially along the western margin and in major outlet glacier basins where W is capped at 2 m; and intermediate regions where $\mathrm{0}<W<\mathrm{2}$ m. The color scheme in Fig. 12 was chosen for comparison to Fig. 11 from MacGregor et al. (2015), who presented a synthesis of Greenland's basal state from observations and models. The agreement is generally excellent (though the two results are not fully independent, since models similar to CISM were part of the synthesis). The largest area of discrepancy is the upstream part of NEGIS, which is thawed or uncertain in the synthesis but frozen in the model.

Figure 12 Basal water state at the end of a 50 kyr Greenland spin-up. Red: fully saturated bed with basal water depth W capped at the upper limit of 2 m. Light blue: partly saturated bed with $\mathrm{0}<W<\mathrm{2}$ m. Dark blue: frozen bed with W=0. The color scale is based on Fig. 11 of MacGregor et al. (2015). The gray background shows ice-free regions.

The intermediate regions in Fig. 12 have striped patterns with century-scale temporal variability, consistent with the band-like patterns in driving and basal shear stresses found by Sergienko et al. (2014) in regions with rapid basal sliding. Sergienko et al. (2014) attributed these patterns to instabilities related to subglacial water. A more detailed investigation would be needed to determine how robust these features are in the model, and to what extent they are real features as opposed to numerical artifacts.

As shown by Schoof and Hindmarsh (2010), depth-integrated approximations like DIVA are computationally much cheaper than BP while incurring only a small loss of accuracy, mainly in slow-sliding regions. However, their analysis did not consider the effects of vertically varying temperature. In CISM simulations of real ice sheets, the depth-integrated viscosity depends on the vertical temperature profile (see Eqs. 2 and 24), with dynamic effects that are not included in idealized problems with a uniform temperature or flow factor. To assess possible errors associated with the DIVA solver (taking BP to be “truth”), we ran CISM for 10 000 years using the BP solver, with settings and forcing otherwise identical to the DIVA run. After 10 kyr, the ice sheet area and volume are 1.635×10⁶ km² and 2.991×10⁶ km³, respectively, for the BP run, compared to 1.634×10⁶ km² and 2.993×10⁶ km³ for DIVA. Figures 13 and 14 show differences in thickness and surface ice speed, respectively, between the two runs. In the vast majority of the ice sheet, thickness differences are less than 50 m, and velocity differences are less than 20 m yr⁻¹. In several outlet glaciers, however – notably Humboldt Glacier in the northwest – thickness differences exceed 100 m, and velocity differences are greater than 100 m yr⁻¹. These differences might merit further study. Some, although not all, of the larger differences can be attributed to transient behavior. For example, Humboldt Glacier is not fully equilibrated at 10 kyr, and the two simulations could capture the glacier in different phases.

Figure 13 Difference (m) in simulated ice thickness after 10 kyr between (1) a run using the BP velocity solver and (2) a run using the DIVA velocity solver and otherwise configured identically. The gray background shows regions that are ice-free in both simulations.

Figure 14 Difference (m yr⁻¹) in simulated ice speed after 10 kyr between (1) a run using the BP velocity solver and (2) a run using the DIVA velocity solver and otherwise configured identically. The gray background shows regions that are ice-free in both simulations.

Table 3 Computational performance on NCAR's Cheyenne supercomputer for various simulations: DIVA vs. BP solvers, with and without floating ice shelves, and 144 vs. 288 processor cores.

Table 3 shows the cost of these runs on 144 and 288 cores of NCAR's Cheyenne supercomputer. On 144 cores, the DIVA run without ice shelves completes 1350 model years per wall-clock hour or 9.41 model years per core hour. Running on 288 cores, the throughput increases to 2480 model years per wall-clock hour, at slightly reduced efficiency of 8.61 model years per core hour. With the BP solver and an otherwise identical configuration (including 11 vertical velocity levels), the throughput decreases by nearly a factor of 50. Table 3 also shows the cost of DIVA runs with ice shelves, which are described next. Although we have found that solver convergence can be slower when floating ice is present, the runs presented here have a similar cost with or without ice shelves.

5.2 Simulation with ice shelves

The simulation with ice shelves is configured like the no-shelf simulation but with these differences:

• No-float calving is replaced by the eigencalving scheme of Sect. 3.5, with k_τ=0.0025 m yr⁻¹ Pa⁻¹ and w₂=25. These parameters give a low calving rate for along-flow tensile stresses and a much higher rate for across-flow stresses.

• Beneath floating ice, we apply the depth-dependent melt rate of Sect. 3.6, with a basal mass balance b=0 at $z=-\mathrm{200}$ m. The freezing rate is 3 m yr⁻¹ above $z=-\mathrm{100}$ m, with a linear ramp between −200 and −100 m. The melt rate is 100 m yr⁻¹ below $z=-\mathrm{500}$ m, with a linear ramp between −200 m and −500 m. The scale for reducing melting in shallow cavities is H₀=20 m.

Choi et al. (2017) used a similar combination of stress-based calving and depth-dependent sub-shelf melting in ISSM to simulate glacier evolution in the NEGIS region. For the CISM simulations, calving and melting parameters were tuned to improve agreement with Greenland's observed ice shelves (or lack thereof). For example, higher values of k_τ were found to remove existing ice shelves in northern Greenland. With the chosen value of k_τ, lower values of w₂ are insufficient to calve long, unrealistic ice tongues along the southeast coast. Similarly, sub-shelf freezing is needed near the surface to prevent substantial calving-front retreat for the Petermann and 79 North glaciers. While near-surface freezing is physically plausible, basal freezing in the simulation might be compensating for a negative SMB bias (Noël et al., 2018) or excessive calving. A maximum melt rate of 100 m yr⁻¹ exceeds inferred melt rates of ∼30 m yr⁻¹ near the grounding line (GL) of Petermann Glacier (Cai et al., 2017), but smaller rates in the model permit unrealistic GL advance. Simulated GL locations could differ, however, if the model were run at different resolution or with a grounding-line parameterization (Gladstone et al., 2010; Leguy et al., 2014).

By the end of a 50 000-year simulation with ice shelves, the total ice area and volume reach equilibrium values of 1.67×10⁶ and 3.08×10⁶ km², respectively. The area agrees closely with the observational data set, but the volume is high by about 4 %. The floating area more than doubles compared to observations, from 3.7×10³ to 9.2×10³ km². The additional floating ice is contained in many small shelves along the coasts.

Figure 15 Difference (m) between (1) simulated ice thickness at the end of a 50 kyr Greenland spin-up with ice shelves and (2) observed thickness from Morlighem et al. (2014). The observed thickness is the model initial condition. The thickness difference is shown for grounded ice only. Regions shaded light green are floating at the end of the simulation, and black contours enclose regions that are floating in observations. The three boxes indicate floating regions highlighted in Fig. 17: Petermann Glacier (red), 79 North Glacier and Zachariae Isstrom (dark green), and Kangerlussuaq Glacier (dark blue). The gray background shows regions that are ice-free at both the beginning and end of the simulation.

Figure 15 shows (for grounded ice) the thickness difference between the simulation and observations. Simulated floating ice is shaded light green, and observed floating ice is enclosed by black contours. The pattern of thickness errors is generally similar to the run without ice shelves (Fig. 11), with positive biases in the northeast interior and along the southeast coast, and negative biases along the northern and northwest margins. The run with shelves, however, has additional positive biases in east central and west central Greenland, upstream of large outlet glaciers such as Jakobshavn in the west and Helheim and Kangerlussuaq in the southeast. The basins that have thickened relative to the no-shelf simulation generally have ice shelves at their termini, where no shelves exist in reality. Although these shelves are small, they appear to provide buttressing that leads to thickening upstream. Figure 16 shows the difference in surface ice speed between the simulations with and without ice shelves. Differences are largest in outlet glaciers with unrealistic floating termini. The general pattern is of faster speeds near glacier termini in the run with shelves (where ice can be much thicker than in the no-shelf simulation), with slower speeds upstream.

Figure 16 Difference (m yr⁻¹) in simulated surface ice speed at the end of 50 kyr Greenland spin-up runs with and without ice shelves. Black contours enclose regions that are floating in the run with ice shelves. The gray background shows regions that are ice-free in both simulations.

Figure 17 Close-ups of three regions with floating ice at the end of a 50 kyr Greenland spin-up with ice shelves: Petermann Glacier (a), 79 North Glacier and Zachariae Isstrom (b), and Kangerlussuaq Glacier (c). These regions are enclosed by boxes in Fig. 15. The shaded color shows the bed topography (m). Black and green contours enclose regions that are floating at the start and end of the simulation, respectively. Calving fronts are located at the north, northeast, and southeast shelf boundaries, respectively, from top to bottom; grounding lines are located at the south, southwest, and northwest shelf boundaries, respectively.

Figure 17 shows closeups of three regions enclosed by boxes in Fig. 15. These regions illustrate floating ice simulated in the vicinity of Petermann Glacier in the northwest, 79 North Glacier/Zachariae Isstrom in the northeast, and Kangerlussuaq Glacier in the southeast. The Petermann simulation (Fig. 17a) is the most realistic, with modest retreat of both the grounding line and calving front, and a total area similar to observations. As mentioned above, however, both the CF and GL locations are sensitive to model parameters. The 79 North Glacier ice shelf (upper left of Fig. 17b) shows modest GL retreat with somewhat larger CF retreat. The simulated Zachariae shelf (lower right of Fig. 17b) substantially retreats; we found no combination of calving and basal melting/freezing parameters that maintains Zachariae in its observed location. Figure 17c, showing Kangerlussuaq Glacier, illustrates the model's tendency to grow unobserved shelves. While there is little floating ice in this region in reality, there is a moderate-sized simulated shelf at the fjord mouth, buttressing the glacier upstream.

These results demonstrate a modeling challenge for Greenland: how to simulate existing ice shelves in northern Greenland without permitting unrealistic shelves in other regions. An eigencalving parameterization, suitably tuned, can prevent shelves from extending into the open ocean, but in enclosed bays there are unrealistic shelves with non-negligible dynamic effects. We were unable to eliminate the unobserved shelves in CISM without also removing observed shelves. This result suggests that important calving or melting mechanisms could be missing in the model. For example, Bassis and Ma (2015) suggested that basal melting, which is spatially variable, is a major driver of crevassing and calving. Lacking ocean forcing, the model does not simulate spatially dependent basal melt rates.

Despite these biases, the results show that CISM is capable of simulating ice velocities broadly consistent with observations for both slow-flowing and fast-flowing parts of the Greenland ice sheet, with reasonable computational costs for multimillennial simulations with a higher-order Stokes approximation.

6 Conclusions

We have described CISM v2.1, which includes many innovations to support robust, accurate, and efficient ice sheet simulations in both idealized and real-world applications. The model incorporates a hierarchy of Stokes approximations, including SIA, SSA, depth-integrated higher-order, and 3-D higher-order. To solve the large elliptic systems associated with membrane stresses, CISM has an efficient native Fortran preconditioned conjugate gradient solver, along with links to third-party solver libraries (Trilinos and SLAP). CISM also adds test cases for higher-order models, including ISMIP-HOM and various shelf and stream cases. For each velocity solver and test case, we have verified that CISM solutions are consistent with community benchmarks. That is, the solutions are as accurate as expected given the simplifications in the various approximations.

CISM's structured rectangular grid has limitations; in particular, the grid cannot be refined selectively near grounding lines and other regions requiring high resolution. On the other hand, the structured grid lends itself to straightforward algorithm development, debugging, creation of forcing data sets, and analysis of results. CISM runs efficiently on a 4 km Greenland grid, with throughput of ∼2000 model years per wall-clock hour on NCAR's Cheyenne supercomputer using the depth-integrated higher-order solver. Performance improvements might be needed, however, to support whole-ice-sheet simulations on grids as fine as 1 or 2 km.

CISM has participated in the initMIP-Greenland experiments (Goelzer et al., 2018) for model initialization. When spun up for 50 kyr with modern climate forcing and without floating ice, CISM gives a steady-state simulation in good agreement with the observed ice extent, volume, and velocity structure of the Greenland ice sheet. The quality of the simulation can be attributed, in part, to a higher-order solver that simulates a realistic stress state in fast outlet glaciers, along with a pseudo-plastic basal sliding law that enhances sliding on thawed beds at low elevation. The simulation also benefits from an accurate SMB (from the regional climate model RACMO2), adjusted to inhibit ice sheet advance beyond modern boundaries. When floating ice is allowed, CISM can maintain ice shelves that resemble observed shelves in northern Greenland, but the model also simulates many small shelves that do not exist in reality.

CISM can be used for both stand-alone and coupled ice sheet simulations. CISM v2.1 is the dynamic ice sheet component of CESM2.0, released in June 2018. Compared to earlier CESM versions, CESM2.0 has more sophisticated ice sheet dynamics and physics (as described here), along with interactive coupling capabilities. These improvements will support increasingly realistic simulations of coupled ice sheet evolution.

Appendix A: Matrix assembly for the Blatter–Pattyn approximation

Here, we describe the assembly of the terms in Eq. (9) into a global matrix A and vector b. The four terms in Eq. (9) describe internal ice stresses, basal friction, lateral pressure, and the gravitational driving force, respectively. We begin with the internal stress term, which is the most complex.

A1 Internal ice stresses

The first term on the LHS of each component of Eq. (9) can be rewritten in terms of velocity components:

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle }x:\underset{\mathrm{\Omega }}{\int }\mathrm{2}\mathit{\eta }\left[\begin{array}{ccc}\mathrm{2}\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}& \frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)& \frac{\mathrm{1}}{\mathrm{2}}\frac{\partial u}{\partial z}\end{array}\right]\left\{\begin{array}{c}\frac{\partial \mathit{\phi }}{\partial x}\\ \frac{\partial \mathit{\phi }}{\partial y}\\ \frac{\partial \mathit{\phi }}{\partial z}\end{array}\right\},\\ {\displaystyle }& {\displaystyle }y:\underset{\mathrm{\Omega }}{\int }\mathrm{2}\mathit{\eta }\left[\begin{array}{ccc}\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)& \mathrm{2}\frac{\partial v}{\partial y}+\frac{\partial u}{\partial x}& \frac{\mathrm{1}}{\mathrm{2}}\frac{\partial v}{\partial z}\end{array}\right]\left\{\begin{array}{c}\frac{\partial \mathit{\phi }}{\partial x}\\ \frac{\partial \mathit{\phi }}{\partial y}\\ \frac{\partial \mathit{\phi }}{\partial z}\end{array}\right\},\end{array}\end{array}$$

where brackets denote row vectors and braces denote column vectors. Glissade evaluates Eq. (A1) for each active element. Hexahedral elements have eight nodes, with u and v to be determined at each active node. Inserting the velocity expressions (Eq. 6) into Eq. (A1), we obtain

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }x:\underset{\mathrm{\Omega }}{\int }\mathrm{2}\mathit{\eta }\left({\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial x}}\left(\mathrm{2}\left[{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial x}}\right]\left\{u_j\right\}+\left[{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial y}}\right]\left\{v_j\right\}\right)\right.\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial y}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial y}}\right]\left\{u_j\right\}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial x}}\right]\left\{v_j\right\}\right)\\ {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial z}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial z}}\right]\left\{u_j\right\}\right)\right),\\ {\displaystyle }& {\displaystyle }y:\underset{\mathrm{\Omega }}{\int }\mathrm{2}\mathit{\eta }\left({\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial y}}\left(\mathrm{2}\left[{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial y}}\right]\left\{v_j\right\}+\left[{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial x}}\right]\left\{u_j\right\}\right)\right.\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial x}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial x}}\right]\left\{v_j\right\}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial y}}\right]\left\{u_j\right\}\right)\\ \text{(A2)}& {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial z}}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial z}}\right]\left\{v_j\right\}\right)\right).\end{array}$$

Each row or column vector has eight terms, one for each node of the element. These terms can be evaluated to form a set of four 8×8 element matrices, denoted as K_uu, K_uv, K_vu, and K_vv. Each row of an element matrix is associated with u or v at a given node. The columns in this row contain terms linking the node to u or v at the other nodes of the same element (with the diagonal term linking the node to itself).

In the x component of Eq. (A2), the terms that multiply u_j are given by

$$\begin{array}{}\text{(A3)}& {\displaystyle }\underset{\mathrm{\Omega }}{\int }\mathit{\eta }\left(\mathrm{4}{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial x}}{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial x}}+{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial y}}{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial y}}+{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial z}}{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial z}}\right)\mathrm{d}\mathrm{\Omega }.\end{array}$$

Letting i and j range from 1 to 8, Eq. (A3) gives the 64 terms of the 8×8 element matrix K_uu, which links the u value at each node to the u values at all eight nodes. Similarly, the 64 terms of element matrix K_uv, which links u at each node to v at each of the eight nodes, are given by

$$\begin{array}{}\text{(A4)}& {\displaystyle }\underset{\mathrm{\Omega }}{\int }\mathit{\eta }\left(\mathrm{2}{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial x}}{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial y}}+{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial y}}{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial x}}\right)\mathrm{d}\mathrm{\Omega }.\end{array}$$

Likewise, two 8×8 matrices are associated with the y component of Eq. (A2). The terms of K_vu are

$$\begin{array}{}\text{(A5)}& {\displaystyle }\underset{\mathrm{\Omega }}{\int }\mathit{\eta }\left(\mathrm{2}{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial y}}{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial x}}+{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial x}}{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial y}}\right)\mathrm{d}\mathrm{\Omega },\end{array}$$

and the terms of K_vv are

$$\begin{array}{}\text{(A6)}& {\displaystyle }\underset{\mathrm{\Omega }}{\int }\mathit{\eta }\left(\mathrm{4}{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial y}}{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial y}}+{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial x}}{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial x}}+{\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial z}}{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial z}}\right)\mathrm{d}\mathrm{\Omega }.\end{array}$$

Because of the symmetry of the underlying PDEs, K_uu and K_vv are symmetric, and ${\mathbf{K}}_{uv}={\mathbf{K}}_{vu}^T$. The terms containing z (i.e., the vertical shear stresses) appear only in K_uu and K_vv, whereas the terms containing x and y (i.e., the membrane stresses) appear in all four element matrices.

In the weak form of the equations, each of the 64 coefficients in each element matrix must be integrated over the element. (Since φ varies over the element, the integrands in Eqs. A3–A6 have a different value at each point.) The integrals can be computed exactly by evaluating the integrand at each of eight quadrature points and summing over quadrature points.

We now specify the form of the basis functions and transform these functions to the geometry of the element (which is irregular in the vertical direction because of the sigma coordinate). We can then evaluate the basis function derivatives at quadrature points. Glissade uses trilinear basis functions defined on a reference cube. This cube is centered at the origin $(\mathrm{0},\mathrm{0},\mathrm{0})$ in local reference coordinates $(\widehat{x},\widehat{y},\widehat{z}$). The eight nodes of the reference cube are located at $(\widehat{x},\widehat{y},\widehat{z})=(\pm \mathrm{1},\pm \mathrm{1},\pm \mathrm{1})$. By convention, nodes 1–4 are the nodes of the lower face, proceeding counterclockwise from the southwest corner $(\widehat{x},\widehat{y})=(-\mathrm{1},-\mathrm{1})$. Nodes 5–8 are the nodes of the upper face, also moving counterclockwise from the southwest corner. Thus, we have

$$\begin{array}{}\text{(A7)}& {\displaystyle }\begin{array}{c}{\mathit{\phi }}_{\mathrm{1}}=(\mathrm{1}-\widehat{x})(\mathrm{1}-\widehat{y})(\mathrm{1}-\widehat{z})/\mathrm{8},\\ {\mathit{\phi }}_{\mathrm{2}}=(\mathrm{1}+\widehat{x})(\mathrm{1}-\widehat{y})(\mathrm{1}-\widehat{z})/\mathrm{8},\\ {\mathit{\phi }}_{\mathrm{3}}=(\mathrm{1}+\widehat{x})(\mathrm{1}+\widehat{y})(\mathrm{1}-\widehat{z})/\mathrm{8},\\ {\mathit{\phi }}_{\mathrm{4}}=(\mathrm{1}-\widehat{x})(\mathrm{1}+\widehat{y})(\mathrm{1}-\widehat{z})/\mathrm{8},\\ {\mathit{\phi }}_{\mathrm{5}}=(\mathrm{1}-\widehat{x})(\mathrm{1}-\widehat{y})(\mathrm{1}+\widehat{z})/\mathrm{8},\\ {\mathit{\phi }}_{\mathrm{6}}=(\mathrm{1}+\widehat{x})(\mathrm{1}-\widehat{y})(\mathrm{1}+\widehat{z})/\mathrm{8},\\ {\mathit{\phi }}_{\mathrm{7}}=(\mathrm{1}+\widehat{x})(\mathrm{1}+\widehat{y})(\mathrm{1}+\widehat{z})/\mathrm{8},\\ {\mathit{\phi }}_{\mathrm{8}}=(\mathrm{1}-\widehat{x})(\mathrm{1}+\widehat{y})(\mathrm{1}+\widehat{z})/\mathrm{8}.\end{array}\end{array}$$

For each n, we have φ_n=1 at node n, with φ_n=0 at the other nodes.

The integrands in Eqs. (A3)–(A6) are written in terms of real Cartesian coordinates $(x,y,z)$ rather than reference coordinates $(\widehat{x},\widehat{y},\widehat{z})$. Spatial derivatives in real coordinates are related to derivatives in reference coordinates by

$$\begin{array}{}\text{(A8)}& {\displaystyle }\left\{\begin{array}{c}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{x}}\\ \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{y}}\\ \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{z}}\end{array}\right\}=\left[\begin{array}{ccc}\frac{\partial x}{\partial \widehat{x}}& \frac{\partial y}{\partial \widehat{x}}& \frac{\partial z}{\partial \widehat{x}}\\ \frac{\partial x}{\partial \widehat{y}}& \frac{\partial y}{\partial \widehat{y}}& \frac{\partial z}{\partial \widehat{y}}\\ \frac{\partial x}{\partial \widehat{z}}& \frac{\partial y}{\partial \widehat{z}}& \frac{\partial z}{\partial \widehat{z}}\end{array}\right]\left\{\begin{array}{c}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial x}\\ \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial y}\\ \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial z}\end{array}\right\}=\left[J\right]\left\{\begin{array}{c}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial x}\\ \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial y}\\ \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial z}\end{array}\right\},\end{array}$$

where [J] is the Jacobian of the transformation between coordinate systems. Given the finite-element expansion,

$$\begin{array}{}\text{(A9)}& {\displaystyle }x=\sum _n{\mathit{\phi }}_{\mathrm{n}}{x}_{\mathrm{n}},\end{array}$$

along with the spatial derivatives of φ at $(\widehat{x},\widehat{y},\widehat{z})$ (which are easily derived from Eq. A7), we can compute $\left[J\right(\widehat{x},\widehat{y},\widehat{z}\left)\right]$ as

$$\begin{array}{}\text{(A10)}& {\displaystyle }\left[J\right]=\left[\begin{array}{ccc}{\sum }_{n=\mathrm{1}}^{\mathrm{8}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{x}}{x}_{\mathrm{n}}& {\sum }_{n=\mathrm{1}}^{\mathrm{8}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{x}}{y}_{\mathrm{n}}& {\sum }_{n=\mathrm{1}}^{\mathrm{8}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{x}}{z}_n\\ {\sum }_{n=\mathrm{1}}^{\mathrm{8}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{y}}{x}_{\mathrm{n}}& {\sum }_{n=\mathrm{1}}^{\mathrm{8}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{y}}{y}_{\mathrm{n}}& {\sum }_{n=\mathrm{1}}^{\mathrm{8}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{y}}{z}_n\\ {\sum }_{n=\mathrm{1}}^{\mathrm{8}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{z}}{x}_{\mathrm{n}}& {\sum }_{n=\mathrm{1}}^{\mathrm{8}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{z}}{y}_{\mathrm{n}}& {\sum }_{n=\mathrm{1}}^{\mathrm{8}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{z}}{z}_n\end{array}\right].\end{array}$$

We then invert Eq. (A8) to obtain the basis function derivatives in terms of $(x,y,z)$:

$$\begin{array}{}\text{(A11)}& {\displaystyle }\left\{\begin{array}{c}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial x}\\ \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial y}\\ \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial z}\end{array}\right\}=\left[J^{-\mathrm{1}}\right]\left\{\begin{array}{c}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{x}}\\ \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{y}}\\ \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{z}}\end{array}\right\}.\end{array}$$

The LHS of Eq. (A11) contains the spatial derivatives needed to evaluate Eqs. (A3)–(A6).

Equations (A3)–(A6) also contain the viscosity η, which is computed at each quadrature point. In the BP approximation, η is given by Eq. (2); it is a function of the flow factor A and the effective strain rate defined by Eq. (3). We approximate A by its value at the element center. The (squared) effective strain rate ${\dot{\mathit{\epsilon }}}_{\mathrm{e}}^{\mathrm{2}}$ is evaluated at each quadrature point by summing over strain-rate components. The x components are given by

$$\begin{array}{}\text{(A12)}& {\displaystyle \frac{\partial u}{\partial x}}=\sum _{n=\mathrm{1}}^{\mathrm{8}}{\displaystyle \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial x}}{u}_{\mathrm{n}},{\displaystyle \frac{\partial v}{\partial x}}=\sum _{n=\mathrm{1}}^{\mathrm{8}}{\displaystyle \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial x}}{v}_{\mathrm{n}},\end{array}$$

and similarly for the y and z components. The nodal velocities in Eq. (A12) are values from the previous iteration.

We now have the information needed to compute the integrands (A3)–(A6) at quadrature points. To integrate over a hexahedron, we take a weighted sum of the values at each of eight quadrature points. These points are located at reference coordinates $(\widehat{x},\widehat{y},\widehat{z})=(\pm \mathrm{1}/\sqrt{\mathrm{3}},\pm \mathrm{1}/\sqrt{\mathrm{3}},\pm \mathrm{1}/\sqrt{\mathrm{3}})$. To evaluate an integral of the form

$$\begin{array}{}\text{(A13)}& {\displaystyle }\underset{\mathrm{\Omega }}{\int }\mathit{\eta }\left({\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial z}}{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial z}}\right)\mathrm{d}\mathrm{\Omega }\end{array}$$

over element volume Ω, we compute the sum over quadrature points

$$\begin{array}{}\text{(A14)}& {\displaystyle }\sum _{p=\mathrm{1}}^{\mathrm{8}}{w}_p{\mathit{\eta }}_p{\left({\displaystyle \frac{\partial {\mathit{\phi }}_i}{\partial z}}{\displaystyle \frac{\partial {\mathit{\phi }}_j}{\partial z}}\right)}_p\left|J_p\right|,\end{array}$$

where $\left|J\right|$ is the determinant of the Jacobian (Eq. A10). For this choice of quadrature points, each point has a weight w_p=1.

The terms of the element matrices ${\mathbf{K}}_{uu},{\mathbf{K}}_{uv},{\mathbf{K}}_{vu}$, and K_vv are then inserted into the corresponding global matrices ${\mathbf{A}}_{uu},{\mathbf{A}}_{uv},{\mathbf{A}}_{vu}$, and A_vv. This is mostly a matter of bookkeeping. For example, the first row of K_uu corresponds to a particular node of element $(k,i,j)$ (specifically, the node with indices $(k-\mathrm{1},i-\mathrm{1},j-\mathrm{1})$, given the convention for numbering nodes within elements). This row corresponds to a row of the global matrix A_uu, and each of the eight terms in the row of K_uu is associated with a column of A_uu. Glissade determines the correct column and adds the K_uu term to the corresponding term in A_uu. This process proceeds until the code has looped over all the active elements and filled the global matrices.

If written in full, each global matrix would have as many rows and columns as there are active nodes. These matrices, however, are sparse, with a maximum of 27 nonzero terms per row (corresponding to a node and its 26 nearest neighbors in a hexahedral lattice). Glissade therefore assembles and stores arrays of dimension $(\mathrm{27},nz,nx-\mathrm{1},ny-\mathrm{1})$, where nx, ny, and nz are the grid dimensions. The 27 terms of the first array dimension are arranged according to a geometric convention. For example, suppose we are filling columns for the matrix row corresponding to node $(k,i,j)$. Then, by convention, index 1 refers to the node with coordinates $(k-\mathrm{1},i-\mathrm{1},j-\mathrm{1})$, index 14 refers to the node itself (i.e., the diagonal term of the row), and index 27 refers to the node at $(k+\mathrm{1},i+\mathrm{1},j+\mathrm{1})$ (and similarly for the other indices). After assembly, these arrays can be converted to the form required by a particular linear solver.

The remaining assembly consists of evaluating the other terms in Eq. (9) (i.e., the basal and lateral boundary conditions and the gravitational forcing) and implementing Dirichlet boundary conditions, if applicable. We consider these in turn.

A2 Basal boundary conditions

The basal boundary terms in Eq. (9) are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }x:\underset{{\mathrm{\Gamma }}_{\mathrm{B}}}{\int }\mathit{\beta }u{\mathit{\phi }}_{\mathrm{1}}\mathrm{d}\mathrm{\Gamma },\\ \text{(A15)}& {\displaystyle }& {\displaystyle }y:\underset{{\mathrm{\Gamma }}_{\mathrm{B}}}{\int }\mathit{\beta }v{\mathit{\phi }}_{\mathrm{2}}\mathrm{d}\mathrm{\Gamma }.\end{array}$$

The basal face of each cell is a rectangle. To integrate over a rectangle, we sum over four quadrature points lying at $(\widehat{x},\widehat{y})=(\pm \mathrm{1}/\sqrt{\mathrm{3}},\pm \mathrm{1}/\sqrt{\mathrm{3}})$ in a reference square with center (0,0) and vertices $(\pm \mathrm{1},\pm \mathrm{1})$. This reference square is the 2-D analog of the reference cube discussed above. We define four bilinear basis functions on the square (cf. Eq. A7):

$$\begin{array}{}\text{(A16)}& {\displaystyle }\begin{array}{c}{\mathit{\phi }}_{\mathrm{1}}=(\mathrm{1}-\widehat{x})(\mathrm{1}-\widehat{y})/\mathrm{4},\\ {\mathit{\phi }}_{\mathrm{2}}=(\mathrm{1}+\widehat{x})(\mathrm{1}-\widehat{y})/\mathrm{4},\\ {\mathit{\phi }}_{\mathrm{3}}=(\mathrm{1}+\widehat{x})(\mathrm{1}+\widehat{y})/\mathrm{4},\\ {\mathit{\phi }}_{\mathrm{4}}=(\mathrm{1}-\widehat{x})(\mathrm{1}+\widehat{y})/\mathrm{4}.\end{array}\end{array}$$

Given these basis functions and their spatial derivatives, we can compute the Jacobian for the transformation between the reference square and the rectangular cell face using the 2-D versions of Eqs. (A10) and (A11):

$$\begin{array}{}\text{(A17)}& {\displaystyle }\left[J\right]=\left[\begin{array}{cc}{\sum }_{n=\mathrm{1}}^{\mathrm{4}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{x}}{x}_{\mathrm{n}}& {\sum }_{n=\mathrm{1}}^{\mathrm{4}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{x}}{y}_{\mathrm{n}}\\ {\sum }_{n=\mathrm{1}}^{\mathrm{4}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{y}}{x}_{\mathrm{n}}& {\sum }_{n=\mathrm{1}}^{\mathrm{4}}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{y}}{y}_{\mathrm{n}}\end{array}\right],\end{array}$$

$$\begin{array}{}\text{(A18)}& {\displaystyle }\left\{\begin{array}{c}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial x}\\ \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial y}\end{array}\right\}=\left[J^{-\mathrm{1}}\right]\left\{\begin{array}{c}\frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{x}}\\ \frac{\partial {\mathit{\phi }}_{\mathrm{n}}}{\partial \widehat{y}}\end{array}\right\}.\end{array}$$

The integrand at a quadrature point has the form βφ_iφ_j, where the second φ term arises from the finite-element expansion of u at a quadrature point:

$$\begin{array}{}\text{(A19)}& {\displaystyle }u=\sum _{n=\mathrm{1}}^{\mathrm{4}}{u}_{\mathrm{n}}{\mathit{\phi }}_{\mathrm{n}}.\end{array}$$

We determine β at quadrature points from the values at cell vertices:

$$\begin{array}{}\text{(A20)}& {\displaystyle }\mathit{\beta }=\sum _{n=\mathrm{1}}^{\mathrm{4}}{\mathit{\beta }}_{\mathrm{n}}{\mathit{\phi }}_{\mathrm{n}}.\end{array}$$

The integral over a cell is then computed as a sum over quadrature points:

$$\begin{array}{}\text{(A21)}& {\displaystyle }\sum _{\mathrm{p}=\mathrm{1}}^{\mathrm{4}}{w}_{\mathrm{p}}{\mathit{\beta }}_{\mathrm{p}}\left({\mathit{\phi }}_i{\mathit{\phi }}_j{)}_{\mathrm{p}}\right|J_{\mathrm{p}}|,\end{array}$$

where w_p=1 for each point. This procedure yields a 4×4 matrix describing the connections between each vertex and its neighbors in the cell. Since the x term in Eq. (A15) contains u but not v, and the y term contains v but not u, we form 2-D matrices K_uu and K_vv, but not K_uv and K_vu. Each term of K_uu is then inserted into the global matrix A_uu, and similarly for K_vv into A_vv.

This assembly method tends to smooth the β field. If it is desired to resolve sharp discontinuities in β, as in the stream test problem of Sect. 4.3, Glissade supports a local assembly method in which the basal friction at a particular vertex depends on the value of β at that vertex alone.

A3 Lateral boundary conditions

The lateral boundary terms in Eq. (9) are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }x:\underset{{\mathrm{\Gamma }}_{\mathrm{L}}}{\int }p{n}_{\mathrm{1}}{\mathit{\phi }}_{\mathrm{1}}\text{ d}\mathrm{\Gamma },\\ \text{(A22)}& {\displaystyle }& {\displaystyle }y:\underset{{\mathrm{\Gamma }}_{\mathrm{L}}}{\int }p{n}_{\mathrm{2}}{\mathit{\phi }}_{\mathrm{2}}\text{ d}\mathrm{\Gamma }.\end{array}$$

Since these terms are independent of u and v, they contribute to the load vectors b_u and b_v on the right-hand side of Eq. (16). They are integrated over the lateral faces of cells (either grounded or floating) that border the ocean.

The lateral faces bordering the ocean are quadrilaterals that can be mapped to a reference square. The integral over each face is found by summing over four quadrature points. Basis functions are given by Eq. (A16), and the Jacobian of the reference square is found using Eq. (A17). The ice thickness H at each quadrature point is evaluated using

$$\begin{array}{}\text{(A23)}& {\displaystyle }H=\sum _{n=\mathrm{1}}^{\mathrm{4}}{H}_{\mathrm{n}}{\mathit{\phi }}_{\mathrm{n}},\end{array}$$

where the H_n are nodal values interpolated from cell centers.

The integrands have the form pφ, where p is the vertically averaged net pressure normal to the ice edge, given by Eq. (13). The integral of the pressure terms over a lateral face is computed as a sum over quadrature points:

$$\begin{array}{}\text{(A24)}& {\displaystyle }\sum _{p=\mathrm{1}}^{\mathrm{4}}\pm w_{\mathrm{p}}{p}_{\mathrm{p}}\left({\mathit{\phi }}_{\mathrm{i}}{)}_{\mathrm{p}}\right|J_{\mathrm{p}}|,\end{array}$$

where the sign depends on the orientation of the face. The resulting pressure terms are inserted into the load vector (either b_u or b_v, depending on the orientation) in the rows associated with each of the four nodes of the face.

A4 Gravitational driving stress

The gravitational forcing terms in Eq. (9) are

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }x:\underset{\mathrm{\Omega }}{\int }\mathit{\rho }g{\displaystyle \frac{\partial s}{\partial x}}{\mathit{\phi }}_{\mathrm{1}}\text{ d}\mathrm{\Omega },\\ \text{(A25)}& {\displaystyle }& {\displaystyle }y:\underset{\mathrm{\Omega }}{\int }\mathit{\rho }g{\displaystyle \frac{\partial s}{\partial y}}{\mathit{\phi }}_{\mathrm{2}}\text{ d}\mathrm{\Omega }.\end{array}$$

To compute these terms, we evaluate $\partial s/\partial x$ and $\partial s/\partial y$ at each active vertex, typically using Eq. (14) or its y analog.

The integrals in Eq. (A25) are over 3-D elements. Each hexahedral element is mapped to a reference cube as described above. Given $\partial s/\partial x$ at the vertices of a cell, the surface slope terms at quadrature points are

$$\begin{array}{}\text{(A26)}& {\displaystyle \frac{\partial s}{\partial x}}=\sum _{n=\mathrm{1}}^{\mathrm{8}}{\mathit{\phi }}_{\mathrm{n}}{\left({\displaystyle \frac{\partial s}{\partial x}}\right)}_n,{\displaystyle \frac{\partial s}{\partial y}}=\sum _{n=\mathrm{1}}^{\mathrm{8}}{\mathit{\phi }}_{\mathrm{n}}{\left({\displaystyle \frac{\partial s}{\partial y}}\right)}_n,\end{array}$$

where the basis functions φ are given by Eq. (A7) and the spatial derivatives are derived from Eqs. (A10) and (A11). The integral of $\mathit{\rho }g\frac{\partial s}{\partial x}\mathit{\phi }$ over an element is evaluated as a sum over quadrature points:

$$\begin{array}{}\text{(A27)}& {\displaystyle }\sum _{\mathrm{p}=\mathrm{1}}^{\mathrm{8}}{w}_{\mathrm{p}}\mathit{\rho }g{\left({\displaystyle \frac{\partial s}{\partial x}}\right)}_{\mathrm{p}}\left({\mathit{\phi }}_i{)}_{\mathrm{p}}\right|J_{\mathrm{p}}|,\end{array}$$

and similarly for the $\partial s/\partial y$ term. Glissade inserts these terms into the load vectors b_u and b_v.

As described above for the assembly of β terms, this method tends to smooth the surface elevation gradient field, averaging neighbor values of $\partial s/\partial x$ and $\partial s/\partial y$ into the driving stress at each vertex. For problems with sharp variations in surface gradients, we have found the solver to be more robust when the driving stress at each vertex depends on $\partial s/\partial x$ and $\partial s/\partial y$ at that vertex alone. Thus, Glissade also supports a local assembly method for the driving stress.

A5 Dirichlet boundary conditions

Once the matrix has been assembled, it may need to be adjusted for Dirichlet boundary conditions (i.e., prescribed velocity values at certain nodes). A common Dirichlet condition is to set $u=v=\mathrm{0}$ at the bed to enforce a no-slip boundary condition. (A no-slip condition can also be enforced by setting the basal traction coefficient β to a large value, but formally this is not a Dirichlet condition.) Also, it may be desirable to set u and v to observed values at certain locations, as in the Ross Ice Shelf test case (Sect. 4.4).

Suppose that at node $(k,i,j)$ we have u=u_c and v=v_c, where u_c and v_c are prescribed values. Let nr be the row of A_uu associated with this node, and let nc range over the columns with nonzero entries in this row. To enforce the Dirichlet condition, we set ${\mathbf{A}}_{uu}(nr,nc)={\mathbf{A}}_{vv}(nr,nc)=\mathrm{0}$ for all values of nc except nc=nr (the diagonal term); we set ${\mathbf{A}}_{uu}(nr,nr)={\mathbf{A}}_{vv}(nr,nr)=\mathrm{1}$. In addition, we set ${\mathbf{A}}_{uv}(nr,nc)={\mathbf{A}}_{vu}(nr,nc)=\mathrm{0}$ for all nc, since these two matrices do not contain any terms on the diagonal of the full global matrix (i.e., A in Eq. 15). On the right-hand side, we set b_u(nr)=u_c and b_v(nr)=v_c. These operations convert the matrix rows associated with node $(k,i,j)$ to the equations $\mathrm{1}\cdot u=u_{\mathrm{c}},\mathrm{1}\cdot v=v_{\mathrm{c}}$, which clearly have the desired solutions u_c and v_c.

A further step is needed to maintain matrix symmetry, as required for the PCG solver. Consider the term A_uu(nr,nc), where nc is a column associated with a neighboring node. We have already set ${\mathbf{A}}_{uu}(nr,nc)=\mathrm{0}$, so we need to set ${\mathbf{A}}_{uu}(nc,nr)=\mathrm{0}$ to maintain symmetry. The Dirichlet condition is u(nr)=u_c. Thus, we can replace b_u(nc) with ${\mathit{b}}_u\left(nc\right)-{\mathbf{A}}_{uu}(nc,nr)u_{\mathrm{c}}$ and set ${\mathbf{A}}_{uu}(nc,nr)=\mathrm{0}$ without altering the problem. We do this for all the terms in the columns associated with node $(k,i,j)$ (i.e., all the terms multiplied by u_c or v_c in the matrix–vector product). Thus, both the rows and the columns associated with node $(k,i,j)$ are filled with zeros, except for the diagonal term, and the full global matrix remains symmetric.