Description and validation of the ice-sheet model Yelmo (version 1.0)

We describe the physics and features of the ice-sheet model Yelmo, an open-source project intended for collaborative development. Yelmo is a thermomechanical model, solving for the coupled velocity and temperature solutions of an ice sheet simultaneously. The ice dynamics are currently treated via a “hybrid” approach combining the shallow-ice and shallow-shelf/shelfy-stream approximations, which makes Yelmo an apt choice for studying a wide variety of problems. Yelmo's main innovations lie in its flexible and user-friendly infrastructure, which promotes portability and facilitates long-term development. In particular, all physics subroutines have been designed to be self-contained, so that they can be easily ported from Yelmo to other models, or easily replaced by improved or alternate methods in the future. Furthermore, hard-coded model choices are eschewed, replaced instead with convenient parameter options that allow the model to be adapted easily to different contexts. We show results for different ice-sheet benchmark tests, and we illustrate Yelmo's performance for the Antarctic ice sheet.

1 Introduction

The field of continental-scale ice-sheet modeling started with a handful of pioneering models (e.g., Huybrechts et al., 1988; Ritz et al., 1997; Greve, 1997a). These models were computationally efficient for the resources available at the time. Typical grid resolutions were on the order of 20–40 km, and generally the shallow ice approximation (SIA) was used to solve the ice dynamics. These classic models have been most useful for long-timescale paleo simulations in part because they are fast but also because they are relatively simple in design, usually relying on low-tech solutions to numerical problems. Most of these models were designed before the era of the high-performance computing cluster, which made it challenging to build models otherwise.

Nowadays, a large number of ice-sheet models exist, supported by a growing and active community of developers. Models today represent a broad spectrum of approaches that incorporate different levels of physical complexity and computational ingenuity. These models include hybrid approaches that heuristically combine the SIA with the shallow shelf approximation (SSA) (e.g., Bueler and Brown, 2009; Winkelmann et al., 2011; Pollard and DeConto, 2012; Pattyn, 2017; Quiquet et al., 2018) and higher-order approximations (e.g., Goldberg, 2011; Cornford et al., 2013; Hoffman et al., 2018; Lipscomb et al., 2019), including full Stokes solutions (e.g., Larour et al., 2012; Gagliardini et al., 2013). Newer models often feature finite-element/finite-volume methods (e.g.,  Larour et al., 2012; Gagliardini et al., 2013; Hoffman et al., 2018) or adaptive mesh refinement (Cornford et al., 2013), which allows simulation of complex terrain and very high resolution where it is needed (e.g., at the grounding line in Antarctica). While more complex models are driving advances in our understanding of the physics and relevant processes of ice sheets over a range of timescales, simpler and thus faster methods are still required to understand the evolution of the ice sheets on multimillennial timescales.

Here we introduce the ice-sheet model Yelmo¹, which is intended to provide access to complex and robust model physics through an intuitive model design. It is a hybrid ice-dynamics model that is easy to use and configure. We expect that Yelmo will be useful for long-timescale simulations of the continental ice sheets, coupled climate–ice-sheet modeling, ensemble simulations and uncertainty studies, as well as for teaching. Below, we first describe the model design (Sect. 2), followed by the physics (Sect. 3), timestepping approach (Sect. 4) and application programming interface (API; Sect. 5). Then, we present results for several benchmark experiments to validate the model performance (Sect. 6) and simulations of the present-day and glacial Antarctic ice sheet (Sect. 7), followed by the conclusions (Sect. 8).

2 Model design

Yelmo has been inspired and largely derived from classical ice-sheet models that have been used successfully for many years – with the most in common with GRISLI (Ritz et al., 1997; Quiquet et al., 2018) and SICOPOLIS (Greve, 1997a, 2019). However, in contrast to many models, Yelmo was designed from scratch to run as a modular library that can be called by other programs rather than as a stand-alone executable. The strict application of this philosophy drove many design choices and allowed us to develop a robust ice-sheet model library with a clear API that would be difficult to develop in an ad hoc way later. Thus, developing this framework was a primary reason to build a new model, rather than continuing the development of other active projects such as GRISLI and SICOPOLIS.

Yelmo is written in Fortran 2003, which provides continuity from previous code bases and supports the fact that clarity and readability of the code are important features. Like SICOPOLIS and other models, we have opted for “low-tech” solutions whenever possible, meaning that internally coded routines are preferred, and, thus, the external dependencies of the model are kept to a minimum. This ensures that the algorithms used remain accessible and easily changeable. Nonetheless, Yelmo has two key dependencies: the NetCDF library for convenient community-standard input/output capability and the Library of Iterative Solvers for linear systems (https://www.ssisc.org/lis/, last access: 2 May 2020) (Lis, Nishida, 2010), which is used for solving the elliptical SSA equations. The latter can be compiled with OpenMP parallelization active, which can speed up this computationally intensive step.

Yelmo has been designed to be user friendly (i.e., straightforward to understand), accessible, portable and adaptable. These features were facilitated by the design choice to separate what we call the “model accounting” from the model physics itself and by following an object-oriented approach. There are no global variables in Yelmo (except for a few global constants related to the general physics of the planet being simulated), which means that variables and parameters are saved together in containers (called derived types in Fortran) specific to each of the main Yelmo components, such as dynamics or thermodynamics, as described in the sections below. These containers make up the individual components of the overall Yelmo object (itself a container) that contains all of the variables, parameters and information needed to simulate a given domain of the ice sheet. Multiple instances of the Yelmo object can therefore be defined in a program (e.g., one Yelmo-object instance for Greenland and one for Antarctica), and each one will operate fully in isolation from the others. This is the model accounting, which is of a specific design built into Yelmo and is thus the only part not easily portable to other models.

The model physics, meanwhile, consists of subroutines that are fully portable and, whenever possible, only rely on native data types (e.g., scalar, vector and array). In other words, the specific, non-portable design structure of the Yelmo object does not contaminate the physics subroutines, since the necessary variables and parameters are always passed as arguments. This approach requires that all input and output to subroutines must be defined as arguments. Each argument must further always be given an intent characteristic (in Fortran, the intent of an argument can be one of IN, OUT or INOUT), which ensures that only the variables destined for output from the routine can be modified inside it. This approach not only aides debugging and provides programmatic safety but also provides a clear blueprint to users of what each subroutine does. Most importantly, the subroutines are thus fully self-contained and can be used in other programs and contexts, as long as the correct arguments can be provided.

Concerning the model accounting, the Yelmo object contains all parameters and variables needed to run a given domain. For clarity and convenience, it has been divided into four components: topography, dynamics, material properties and thermodynamics (Fig. 1). Each component has an associated set of functions to load parameters, allocate and initialize the variables, update the variables (i.e., the actual physics calculation step), and finally terminate the instance of the component at the end of the program. This pattern is followed for all four components and represents the component-level API.

Figure 1 Overview of the Yelmo model structure highlighting state variables in the four components: topography, dynamics, material and thermodynamics, as well as the boundary conditions required to run the model. The thick black border for boundary variables indicates that these fields are never modified internally by Yelmo, while the components with a thin black border or dashed line are allowed to be modified depending on the context. When, for example, the topography is updated (dashed line), no other components are allowed to be modified (solid lines).

Each component contains variables and parameters necessary for the calculation of its specific physics; however each component also relies on the variables defined in other components since the ice sheet is a highly coupled and nonlinear system. The benefit of the somewhat artificial division of components made here is that the use of INTENT statements ensures that variables of a given component can only be modified in the corresponding module. For example, when the update subroutine of the topography module is called, only the object containing topography variables is defined as INTENT=INOUT, while the objects containing dynamics, material and thermodynamics variables are all defined as INTENT=IN. Analogous to the design of the physics subroutines, the use of intent statements here not only makes the model blueprint clear but also enforces consistency with the overall design of the Yelmo structure. The hope is that this will not only make the model more user friendly but that it will also naturally lead to more disciplined model development in the future.

In addition to the four components that contain prognostic and diagnostic model variables, the Yelmo object includes a boundary component, which defines all fields that Yelmo requires as input from external sources (Fig. 1). These fields can be obtained from other coupled models or simply by loading data; however Yelmo does not make any assumptions about their source. The boundary component is defined as INTENT=IN in all modules, so that Yelmo does not have the right to modify them internally. This conceptual isolation of the ice-sheet model serves to ensure that coupling with other models is as straightforward as possible, because it is clear by design which variables should be provided to Yelmo as boundary conditions. This is a key feature of Yelmo in comparison with many other models.

Yelmo also makes use of a working precision variable, which allows for the model to be compiled with any real precision. For most applications, single precision (32 bits) is sufficient. Double precision (64 bits) gives equivalent results for the tests we have made. Nonetheless, this choice is left open to the user.

In terms of model physics, each component of Yelmo was built to work independently, in the sense that a given component is agnostic to the methods used to calculate variables from other components. For example, the temperature and velocity fields are used by the material component without any knowledge about the physics and numerical approximations behind them. This means that sometimes simplifying assumptions cannot be used, even though they may be valid in some cases (such as assuming that the strain rate is only due to SIA terms where the ice sheet is frozen to the bed). However, the benefit is that typically the most general solutions possible have been implemented for each component. Thus, when the physics of one component is changed or upgraded, it is likely that the other modules will not require any modification.

Grid information is also stored in the main Yelmo object, and a single grid is defined for use with all components. Like many ice-sheet models, Yelmo uses the Arakawa C-grid staggering approach (Arakawa and Lamb, 1977) extended to 3D, as shown in Fig. 2. Scalar variables, such as temperature, are defined at the cell centers, which in Yelmo are designated as “aa-nodes”. Velocity components and gradients are calculated on cell edges (“ac-nodes”), and scalar coefficients, like diffusivity in the SIA approach, are calculated on cell corners (“ab-nodes”). The specific numerical discretization of the finite difference equations largely follows the approach of Macayeal (1997). The advantage of this approach is that it benefits from the natural staggering that occurs when calculating gradients (e.g., the surface slope is naturally defined on the ac-nodes), but it also results in greater numerical stability of the model (Macayeal, 1997).

Figure 2 Yelmo staggered grid definition and nomenclature. The (a) horizontal grid assumes constant resolution in the x and y directions, while in the (b) vertical variable resolution is allowed. With any given cell defined as a 3D box, scalar variables are calculated on cell centers (aa-nodes), velocities are calculated on cell faces (ac-nodes, edges in 2D), and scalar coefficients are calculated on cell edges (ab-nodes, corners in 2D). Figure design adapted from Hoffman et al. (2018).

Yelmo requires an evenly spaced Cartesian grid in the horizontal direction, while the vertical component follows a classic sigma-coordinate system (Greve and Blatter, 2009). The vertical axis ζ represents the relative height within the ice sheet, running from ζ=0 at the ice-sheet base to ζ=1 at the ice-sheet surface:

$$\begin{array}{}\text{(1)}& \mathit{\zeta }\left(z\right)=h\left(z\right)/H,\end{array}$$

where z is the elevation relative to present-day sea level, h(z) is the ice thickness up to the elevation z within the ice sheet and H is the total ice thickness. Yelmo can be defined with any specified number of vertical grid points, which can be unevenly spaced. Typically, we have set n_z=20, and the ζ axis is defined with higher resolution near the base and surface of the ice sheet, which is important for resolving thermodynamics and ages accurately. Use of the sigma-coordinate system simplifies the numerics of an evolving domain in the vertical direction and inherently results in higher resolution for grid points with less ice thickness (Greve and Blatter, 2009). Vertical velocities are calculated on ac-nodes in the vertical and aa-nodes in the horizontal, while horizontal velocities are calculated on ac-nodes in the horizontal and aa-nodes in the vertical. Boundary conditions in a vertical column are applied directly at the ice base and ice surface, which correspond to ac-nodes (see Fig. 2).

3 Model physics

Yelmo solves for two prognostic variables using coupled equations of mass and energy conservation: the ice thickness (2D field) and ice temperature (3D field). Velocity (3D vector field) is diagnosed from approximations of ice flow assuming a nonlinear flow law. These equations are described in the subsections below, along with additional considerations related to each component. For more details on the derivation of the equations, thorough explanations can be found in various references (Greve and Blatter, 2009; Cuffey and Paterson, 2010) and thus are not repeated here.

3.1 Topography

The evolution of the ice thickness in the model is determined from mass conservation:

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial H}{\partial t}}=-\mathbf{▿}\cdot H\stackrel{\mathrm{‾}}{\mathbf{u}}+\dot{a}+{\dot{b}}_{\mathrm{g}}+{\dot{b}}_{\mathrm{f}}-\dot{c},\end{array}$$

where H is the ice thickness, $\stackrel{\mathrm{‾}}{\mathbf{u}}=\left(\stackrel{\mathrm{‾}}{u},\stackrel{\mathrm{‾}}{v}\right)$ is the depth-averaged horizontal velocity, $\dot{a}$ is the surface mass balance, ${\dot{b}}_{\mathrm{g}}$ and ${\dot{b}}_{\mathrm{f}}$ are the basal mass balance for grounded and floating ice, respectively, and $\dot{c}$ is the calving rate at floating ice margins. In Yelmo, in order to obtain more accurate mass balance accounting, the advection of ice and source contributions are treated separately as follows. First, a forward Euler explicit method (or optionally an upwind implicit method) is used to solve for the ice thickness without accounting for $\dot{a}$, ${\dot{b}}_{\mathrm{g}}$, ${\dot{b}}_{\mathrm{f}}$ or $\dot{c}$. The depth-averaged horizontal velocity is obtained from the dynamics component from the previous iteration (see Timestepping below). Next the mass balance terms $\dot{a}$, ${\dot{b}}_{\mathrm{g}}$ and ${\dot{b}}_{\mathrm{f}}$ are applied. It should be noted that the basal mass balance of floating ice is a boundary variable for Yelmo (i.e., it is obtained externally and passed to Yelmo), while the basal mass balance of the grounded ice is calculated internally as part of the thermodynamics solver (see Thermodynamics section below).

Yelmo also includes special treatment of grid points at the floating margin of the ice sheet, by making a distinction between ice-covered grid points that are totally and partially filled following Albrecht et al. (2011) and Lipscomb et al. (2019). This is done in a relatively simple, yet effective way to avoid artificially thin ice thickness at the ice margin. For each floating ice-covered grid point that has an ice-free neighbor, the reference ice thickness of the margin point (H_ref) is defined as the minimum thickness of the direct, ice-covered neighbors. This represents the minimum ice thickness for which the cell can be considered completely ice covered. The fraction of ice cover is then defined as $f_{\mathrm{ice}}=\min \left(H/H_{\mathrm{ref}},\mathrm{1}\right)$. Whenever f_ice<1, the grid cell is considered dynamically inactive, which ensures zero ice flux through the downstream edge of a partially filled margin cell. In this way, the ice cell can be filled with ice from upstream, and when the threshold of f_ice=1 is reached, the ice shelf can advance.

In the final mass conservation step, calving $\dot{c}$ is treated at the floating ice margins. Currently, a simple threshold method has been implemented, as well as a threshold+flux method (Peyaud et al., 2007). In both methods, the calving rate applied to the ice sheet is defined following Lipscomb et al. (2019):

$$\begin{array}{}\text{(3)}& \dot{c}={\displaystyle \frac{H_{\mathrm{ref}}-H}{{\mathit{\tau }}_{\mathrm{c}}}},\end{array}$$

where τ_c is the characteristic calving time, usually set to 1–10 years, and H_ref is the margin ice thickness as defined above. Setting τ_c to higher values facilitates ice-shelf growth and thus grounding-line advance in transient, glacial simulations but has little impact on the steady-state distribution of ice shelves for the present day. This calving rate is applied only when the ice thickness of ice margin points is below a threshold value (simple threshold method) or when the ice thickness is below a threshold value and the upstream flux is not sufficient to return the ice thickness to above the threshold (threshold + flux method). For paleo simulations the latter is our preferred method, as it allows for more robust ice shelf advance (Peyaud et al., 2007).

Once the ice thickness has been updated, Yelmo diagnoses whether the ice should be grounded or floating. To facilitate this step, the height above flotation as measured in ice thickness, i.e., how close a grid point is to the Archimedes flotation criterion, is calculated on each aa-node:

$$\begin{array}{}\text{(4)}& H_{\mathrm{g}}=H-{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{sw}}}{\mathit{\rho }}}\max \left(z_{\mathrm{sl}}-z_{\mathrm{b}},\mathrm{0}\right),\end{array}$$

where ρ is the ice density and ρ_sw the seawater density, and z_sl and z_b are the boundary fields of sea level and bedrock elevation, respectively. H_g can thus be positive, zero or negative. When H_g is positive, the ice thickness exceeds the flotation criterion and is considered grounded, while when H_g is zero or negative, the ice is considered floating.

Yelmo also calculates the grounded fraction of each grid point, f_g. On aa-nodes, f_g is only assigned binary values to maintain consistency with the overall grid definition: zero when H_g≤0 or one when H_g>0. However, on ac-nodes, the values of f_g,acx and f_g,acy are determined by linearly interpolating H_g from the two bounding aa-nodes. When both bounding aa-nodes are positive $f_{\mathrm{g},\mathrm{ac}}=\mathrm{1}$, and when both are negative $f_{\mathrm{g},\mathrm{ac}}=\mathrm{0}$. When one aa-node is positive ($H_{{\mathrm{g}}_{+}}$) and one aa-node is negative ($H_{{\mathrm{g}}_{-}}$), the grounded fraction on the ac-node is determined from linear interpolation:

$$\begin{array}{}\text{(5)}& f_{\mathrm{g},\mathrm{ac}}={\displaystyle \frac{H_{{\mathrm{g}}_{+}}}{(H_{{\mathrm{g}}_{+}}-H_{{\mathrm{g}}_{-}})}}.\end{array}$$

Alternatively, it is possible to calculate f_g via subgrid bilinear interpolation of H_g to intermediate points to determine the grounded area fraction. However, this operation is more computationally intensive, and we find that in practice, the simple linear interpolation method is sufficient.

The surface elevation (z_s) is calculated following Pattyn (2017) as

$$\begin{array}{}\text{(6)}& z_{\mathrm{s}}=\max \left[z_{\mathrm{b}}+H,z_{\mathrm{sl}}+(\mathrm{1}-{\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{sw}}}})H\right].\end{array}$$

This approach ensures that the surface elevation solution is consistent with the Archimedes flotation criterion on aa-nodes.

The remaining tasks of the topography component are to diagnose other useful topographic characteristics, such as surface and ice thickness gradients (on ac-nodes) and topographic masks.

3.2 Material

The material component of Yelmo handles the calculation of the rate factor, the strain rate tensor and effective strain rate, the effective viscosity, and, optionally, the age of the ice. Essentially, the material variables make the link between thermodynamics and dynamics, since the rate factor depends on temperature, and the strain rate depends on velocity. No distinction is made between the type of approximation used to solve the dynamics here – rather all equations follow from the more general hydrostatic approximation (Greve and Blatter, 2009).

The effective viscosity, used to determine strain heating in the thermodynamics component, is calculated as

$$\begin{array}{}\text{(7)}& \mathit{\eta }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\dot{\mathit{\epsilon }}}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}-n}{\mathrm{2}n}}\left(A^{-\mathrm{1}/n}\right),\end{array}$$

where $\dot{\mathit{\epsilon }}$ is the effective strain rate, n is the Glen's Flow law exponent (Glen, 1955; Greve and Blatter, 2009), typically set to n=3, and A is the rate factor. The effective strain rate is given by the second invariant of the strain rate tensor (${\dot{\mathit{\epsilon }}}_{\mathrm{ij}}$),

$$\begin{array}{}\text{(8)}& \dot{\mathit{\epsilon }}={\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\dot{\mathit{\epsilon }}}_{\mathrm{ij}}{\dot{\mathit{\epsilon }}}_{\mathrm{ij}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}$$

and the strain rate tensor itself, following index notation, is

$$\begin{array}{}\text{(9)}& {\dot{\mathit{\epsilon }}}_{\mathrm{ij}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial u_i}{\partial x_j}}+{\displaystyle \frac{\partial u_j}{\partial x_i}}\right),i,j=\mathrm{1},\mathrm{2},\mathrm{3}.\end{array}$$

The rate factor, $A\left(x,y,z\right)$, can be prescribed to be a constant value or calculated as a function of ice temperature following an Arrhenius equation:

$$\begin{array}{}\text{(10)}& A\left(T^{\prime }\right)=E_{\mathrm{f}}{A}_{\mathrm{0}}{e}^{-Q_{\mathrm{a}}/R{T}^{\prime }}.\end{array}$$

Here R is the ideal gas constant; A₀ and Q_a are the temperature-dependent rate factor coefficient and activation energy, respectively (see Greve and Blatter, 2009). E_f is a so-called enhancement factor, which is used to approximate the effect of anisotropic flow. In Yelmo, it is possible to specify different values of the enhancement factor for different flow regimes (shear, stream and shelf). The shelf value is prescribed anywhere ice is floating, while the inland value of E_f is a weighted average between the shear and stream value with the weighting given by a diagnosis of the vertical shearing fraction at any given point:

$$\begin{array}{}\text{(11)}& f_{\mathrm{shr}}={\displaystyle \frac{\left({\dot{\mathit{\epsilon }}}_{\mathrm{xz}}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{\mathrm{yz}}^{\mathrm{2}}\right)}{{\dot{\mathit{\epsilon }}}^{\mathrm{2}}}}.\end{array}$$

Typical values of the enhancement factor for the shearing, streaming and shelf regime are E_shr=3.0, E_strm=1.0 and E_shlf=0.7, respectively (Ma et al., 2010).

In addition, it is possible to track the deposition time (i.e., age) or other conservative tracers of the ice using an Eulerian tracer advection model. The general 3D advection equation of a conservative variable X,

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial X}{\partial t}}=-u{\displaystyle \frac{\partial X}{\partial x}}-v{\displaystyle \frac{\partial X}{\partial y}}-w{\displaystyle \frac{\partial X}{\partial z}},\end{array}$$

is solved with a second-order, upwind explicit method. The ice surface boundary condition must be imposed. When tracing the ice deposition time, the ice surface boundary condition is X(t)=t. At the ice base, an initial deposition time is prescribed to be several thousand years before the start of the simulation; however this plays little role in the resulting vertical profile of deposition times. When ice is melting at the base $\left(\dot{b}<\mathrm{0}\right)$, the following flux boundary condition is defined (Rybak and Huybrechts, 2003):

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial X}{\partial t}}=-u_{\mathrm{b}}{\displaystyle \frac{\partial X}{\partial x}}-v_{\mathrm{b}}{\displaystyle \frac{\partial X}{\partial y}}-\dot{b}{\displaystyle \frac{\partial X}{\partial z}}.\end{array}$$

Basal freeze-on is assumed to be negligible. It is well known that Eulerian solvers lose accuracy towards the base of the ice sheet, and therefore this method can only be considered to give a first-order estimate of a conservative tracer (Greve et al., 2002; Rybak and Huybrechts, 2003). It can nonetheless be useful for diagnosing the age of ice, in order to know the timescale of different dynamic properties or to, e.g., impose an age-dependent enhancement factor (Greve, 1997b).

3.3 Dynamics

The Yelmo dynamics component is currently representative of a “hybrid” class of ice-sheet model, treating different modes of ice deformation via a combination of the simplifying shallow-ice and shallow-shelf approximations (SIA and SSA, respectively). In the following, the description of the dynamics equations follows closely the notation and definitions of Greve and Blatter (2009) and Pollard and DeConto (2012).

Yelmo treats the horizontal velocity $u\left(x,y,z\right)$ and $v\left(x,y,z\right)$ as the sum of transport via internal shear (u_i, v_i) and basal sliding (u_b, v_b):

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{\displaystyle }u& {\displaystyle }=u_{\mathrm{i}}+u_{\mathrm{b}}\\ {\displaystyle }v& {\displaystyle }=v_{\mathrm{i}}+v_{\mathrm{b}}.\end{array}\end{array}$$

Here, and analogously for v, u_b(x,y) is vertically constant, and $u_{\mathrm{i}}\left(x,y,z_{\mathrm{b}}\right)=\mathrm{0}$, where the subindex “b” here represents the basal boundary of the ice sheet. It also holds that in the vertical average (denoted by a bar), $\stackrel{\mathrm{‾}}{u}={\stackrel{\mathrm{‾}}{u}}_{\mathrm{i}}+u_{\mathrm{b}}$. To calculate u_i and v_i, Yelmo uses zero-order SIA equations:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{u}_i\left(z\right)& =-\left[\mathrm{2}|\mathit{▿}{\mathit{\tau }}_{\mathrm{d}}{|}^{(n-\mathrm{1})}\underset{z_{\mathrm{b}}}{\overset{z}{\int }}A{\left(z_{\mathrm{s}}-z^{\prime }\right)}^n\mathrm{d}{z}^{\prime }\right]{\mathit{\tau }}_{\mathrm{d},\mathrm{x}}\\ v_i\left(z\right)& =-\left[\mathrm{2}|\mathit{▿}{\mathit{\tau }}_{\mathrm{d}}{|}^{(n-\mathrm{1})}\underset{z_{\mathrm{b}}}{\overset{z}{\int }}A{\left(z_{\mathrm{s}}-z^{\prime }\right)}^n\mathrm{d}{z}^{\prime }\right]{\mathit{\tau }}_{\mathrm{d},\mathrm{y}},\end{array}\end{array}$$

where u_i(z) and v_i(z) are the horizontal components of the SIA velocity as a function of depth at a given location, A is the material rate factor of the ice, which is obtained from the material component (Eq. 10), n is the Glen's Flow law exponent (Glen, 1955; Greve and Blatter, 2009), and ${\mathit{\tau }}_{\mathrm{d}}=\left({\mathit{\tau }}_{\mathrm{d},\mathrm{x}},{\mathit{\tau }}_{\mathrm{d},\mathrm{y}}\right)=\mathit{\rho }gH\mathit{▿}{z}_{\mathrm{s}}$ is the driving stress. In the horizontal plane, the term in brackets is calculated on the ab-nodes for stability and improved mass conservation (Huybrechts et al., 1996), and then it is staggered onto the ac-nodes where it is multiplied with the driving stress. In the vertical plane, the horizontal velocities are calculated at the vertical center of each grid point (aa-nodes). Following Bueler and Brown (2009), we use the SSA solution to calculate the transport implied by sliding at the base (i.e., in regions of ice streams and floating ice shelves):

$$\begin{array}{}\text{(16)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial x}}& \left[H{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{d}}\left(\mathrm{4}{\displaystyle \frac{\partial u_{\mathrm{b}}}{\partial x}}+\mathrm{2}{\displaystyle \frac{\partial v_{\mathrm{b}}}{\partial y}}\right)\right]+{\displaystyle \frac{\partial }{\partial y}}\left[H{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{d}}\left({\displaystyle \frac{\partial u_{\mathrm{b}}}{\partial y}}+{\displaystyle \frac{\partial v_{\mathrm{b}}}{\partial x}}\right)\right]\\ & ={\mathit{\tau }}_{\mathrm{d},\mathrm{x}}-{\mathit{\tau }}_{\mathrm{b},\mathrm{x}}\\ {\displaystyle \frac{\partial }{\partial y}}& \left[H{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{d}}\left(\mathrm{4}{\displaystyle \frac{\partial v_{\mathrm{b}}}{\partial y}}+\mathrm{2}{\displaystyle \frac{\partial u_{\mathrm{b}}}{\partial x}}\right)\right]+{\displaystyle \frac{\partial }{\partial x}}\left[H{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{d}}\left({\displaystyle \frac{\partial u_{\mathrm{b}}}{\partial y}}+{\displaystyle \frac{\partial v_{\mathrm{b}}}{\partial x}}\right)\right]\\ & ={\mathit{\tau }}_{\mathrm{d},\mathrm{y}}-{\mathit{\tau }}_{\mathrm{b},\mathrm{y}},\end{array}\end{array}$$

where $\left({\mathit{\tau }}_{\mathrm{b},\mathrm{x}},{\mathit{\tau }}_{\mathrm{b},\mathrm{y}}\right)=-\mathit{\beta }\left(u_{\mathrm{b}},v_{\mathrm{b}}\right)$ (or in vector notation ${\mathit{\tau }}_{\mathrm{b}}=-\mathit{\beta }{\mathbf{u}}_{\mathrm{b}}$) is the basal stress due to friction. The basal friction coefficient β is set to zero for floating ice shelves and can otherwise be set to a constant value or follow another user-defined formulation (power law, regularized Coulomb, etc.), depending on the context (see basal friction description below for details). The depth-averaged (2D) effective viscosity, which is only used for solving the SSA dynamics, is defined as

$$\begin{array}{}\text{(17)}& {\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{d}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\stackrel{\mathrm{‾}}{B}{\left({\dot{\mathit{\epsilon }}}_{\mathrm{d}}^{\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{\mathrm{0}}^{\mathrm{2}}\right)}^{\frac{\mathrm{1}-n}{\mathrm{2}n}}\end{array}$$

$\stackrel{\mathrm{‾}}{B}=\frac{\mathrm{1}}{H}{\int }_{z_{\mathrm{b}}}^{z_{\mathrm{s}}}{A}^{-\mathrm{1}/n}\mathrm{d}z$ is the vertically averaged ice hardness, ${\dot{\mathit{\epsilon }}}_{\mathrm{d}}$ is the 2D effective strain rate and ${\dot{\mathit{\epsilon }}}_{\mathrm{0}}^{\mathrm{2}}$ is a small regularization factor for avoiding a potential singularity when velocity gradients are zero. The 2D effective strain rate is calculated as a reduced form of the second invariant of the strain rate tensor (Eq. 9) that does not include vertical shear terms:

$$\begin{array}{}\text{(18)}& {\dot{\mathit{\epsilon }}}_{\mathrm{d}}^{\mathrm{2}}={\left({\displaystyle \frac{\partial u_{\mathrm{b}}}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial v_{\mathrm{b}}}{\partial y}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\partial u_{\mathrm{b}}}{\partial x}}{\displaystyle \frac{\partial v_{\mathrm{b}}}{\partial y}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial u_{\mathrm{b}}}{\partial y}}+{\displaystyle \frac{\partial v_{\mathrm{b}}}{\partial x}}\right)}^{\mathrm{2}}.\end{array}$$

In Yelmo, ${\dot{\mathit{\epsilon }}}_{\mathrm{d}}$ is only used for calculating ${\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{d}}$, while the 3D effective strain rate is calculated from the full strain rate tensor in the material component (see Sect. 3.2). Calculating the full tensor during the iterative SSA solution procedure would be much more computationally expensive, while the 2D effective strain rate is already sufficient for the vertically integrated SSA equations (Pollard and DeConto, 2012).

The stress boundary condition imposed at the floating ice front, following Winkelmann et al. (2011) and Greve and Blatter (2009), is

$$\begin{array}{}\text{(19)}& \begin{array}{rl}H{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{d}}& \left(\mathrm{4}{\displaystyle \frac{\partial u}{\partial x}}+\mathrm{2}{\displaystyle \frac{\partial v}{\partial y}}\right)n_x+H{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{d}}\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)n_y\\ & =\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }g{H}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{sw}}g{H}_{\mathrm{o}}^{\mathrm{2}}\right)n_x\\ H{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{d}}& \left(\mathrm{4}{\displaystyle \frac{\partial v}{\partial y}}+\mathrm{2}{\displaystyle \frac{\partial u}{\partial x}}\right)n_y+H{\stackrel{\mathrm{‾}}{\mathit{\eta }}}_{\mathrm{d}}\left({\displaystyle \frac{\partial v}{\partial x}}+{\displaystyle \frac{\partial u}{\partial y}}\right)n_x\\ & =\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\rho }g{H}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\rho }}_{\mathrm{sw}}g{H}_{\mathrm{o}}^{\mathrm{2}}\right)n_y.\end{array}\end{array}$$

The depth of the seawater up to the flotation depth, H_o, is defined as $H_{\mathrm{o}}=\min \left(\max \left(z_{\mathrm{sl}}-z_{\mathrm{b}},\mathrm{0}\right),\frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{sw}}}H\right)$. This is the depth of the ocean directly adjacent to the ice sheet, which acts to reduce the outward pressure at the floating ice margin. In contrast to Winkelmann et al. (2011), this boundary condition is not currently used in Yelmo for grounded ice, where Eq. (16) applies.

The SSA equations are nonlinear, elliptical, partial differential equations with nonlocal solutions. Yelmo uses Lis for the numerical solution using the biconjugate gradient method. The subroutine to discretize the equations and to call Lis was ported from the latest SICOPOLIS version 5-dev (Greve, 2019; Rückamp et al., 2019) and subsequently modified for model design choices in Yelmo. We use a Picard iteration method to account for the nonlinear dependence of the effective viscosity (η_d) and, potentially, the basal friction coefficient (β) on velocity. Convergence of the SSA solution is tested using the L² relative error norm (Gagliardini et al., 2013):

$$\begin{array}{}\text{(20)}& {\mathit{\delta }}_{u,v}={\displaystyle \frac{\mathrm{2}\sqrt{\sum {\left(u_{\mathrm{1}}-u_{\mathrm{0}}\right)}^{\mathrm{2}}+\sum {\left(v_{\mathrm{1}}-v_{\mathrm{0}}\right)}^{\mathrm{2}}}}{\sqrt{\sum {\left(u_{\mathrm{1}}+u_{\mathrm{0}}\right)}^{\mathrm{2}}+\sum {\left(v_{\mathrm{1}}+v_{\mathrm{0}}\right)}^{\mathrm{2}}}}},\end{array}$$

where (u₁,v₁) and (u₀,v₀) are the velocity solutions for the current and previous iterations, respectively, and the sum is made over all grid points with nonzero velocity being considered by the SSA solver. By default, we consider a convergence limit of ${\mathit{\delta }}_{u,v}={\mathrm{10}}^{-\mathrm{2}}$, which is typically achieved within 1–10 iterations, depending on the context. This limit can be specified by the user.

The result of solving the above equations is the hybrid, 3D horizontal velocity field (u, v). The vertical velocity w can then be diagnosed by applying a kinematic boundary condition at the base and integrating the continuity equation for incompressible flow (Greve and Blatter, 2009), from z_b to z,

$$\begin{array}{}\text{(21)}& w\left(z\right)=\dot{b}-\left(u_{\mathrm{b}}{\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial x}}+v_{\mathrm{b}}{\displaystyle \frac{\partial z_{\mathrm{b}}}{\partial y}}\right)-\underset{z_{\mathrm{b}}}{\overset{z}{\int }}\left({\displaystyle \frac{\partial u}{\partial x}}+{\displaystyle \frac{\partial v}{\partial y}}\right)\mathrm{d}{z}^{\prime }.\end{array}$$

The vertical velocity is naturally defined on the ac-nodes in the vertical plane, analogous to the horizontal velocity in the horizontal plane. The above dynamics update results in a 3D hybrid velocity field (u, v, w).

Basal frictional stress, as it appears in the SSA elliptical equations, is defined as

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

where β represents the basal friction coefficient (with units of Pa yr m⁻¹), which can be defined in several ways. β is prescribed to be zero for floating ice and otherwise can be set to a constant or a spatially varying field, and, depending on the formulation used, it can depend on velocity itself. For this reason, we also define c_b as the bed friction coefficient, which we consider to only provide information about conditions at the physical bed (the nature of basal sediments, basal hydrology, effective pressure, etc.), independent of velocity. In the model, therefore, β is defined as

$$\begin{array}{}\text{(23)}& \mathit{\beta }=c_{\mathrm{b}}f\left({\mathbf{u}}_{\mathrm{b}}\right).\end{array}$$

Thus in all formulations implemented in Yelmo, the term f(u_b) has units of years per meter (yr m⁻¹), and the coefficient c_b has units of pascals (Pa), which helps to facilitate its physical interpretation.

Most commonly, β is defined using a linear (e.g. Quiquet et al., 2018), power-law (e.g. Pattyn, 2017), pseudo-plastic power-law (e.g. Aschwanden et al., 2013) or regularized-Coulomb (Joughin et al., 2019) formulation. The linear and power-law formulations are contained within the pseudo-plastic power-law formulation, so only the latter and the regularized-Coulomb formulation are needed to represent all four cases.

The pseudo-plastic power-law formulation (Schoof, 2010; Aschwanden et al., 2013) is

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

and thus $\mathit{\beta }=c_{\mathrm{b}}{u_{\mathrm{0}}}^{-q}{\left|{\mathbf{u}}_{\mathrm{b}}\right|}^{q-\mathrm{1}}$, with the pseudo-plastic exponent $q\in \left(\mathrm{0},\mathrm{1}\right)$ and threshold speed u₀. This expression results in purely plastic friction for q=0, linear friction for q=1 and power-law friction for $\mathrm{0}<q<\mathrm{1}$. With q=1 and u₀=1, for example, β=c_b, and friction scales linearly with velocity. To obtain the power-law formulation used in the original MISMIP experiments (Pattyn et al., 2012), the following parameter values can be prescribed: $q=\mathrm{1}/\mathrm{3}$, u₀=1 m yr⁻¹ and $c_{\mathrm{b}}=\mathrm{3.165176}\times {\mathrm{10}}^{\mathrm{4}}$ Pa.

Alternatively, the regularized Coulomb law (Schoof, 2005; Brondex et al., 2019; Joughin et al., 2019) is defined as

$$\begin{array}{}\text{(25)}& {\mathit{\tau }}_{\mathrm{b}}=-c_{\mathrm{b}}{\left({\displaystyle \frac{\left|{\mathbf{u}}_{\mathrm{b}}\right|}{\left|{\mathbf{u}}_{\mathrm{b}}\right|+u_{\mathrm{0}}}}\right)}^q{\displaystyle \frac{{\mathbf{u}}_{\mathrm{b}}}{\left|{\mathbf{u}}_{\mathrm{b}}\right|}},\end{array}$$

and thus $\mathit{\beta }=c_{\mathrm{b}}{\left(\left|{\mathbf{u}}_{\mathrm{b}}\right|+u_{\mathrm{0}}\right)}^{-q}{\left|{\mathbf{u}}_{\mathrm{b}}\right|}^{q-\mathrm{1}}$. Again q is the nonlinear exponent, and u₀ is an empirical threshold speed that dictates the transition from Coulomb friction when cavitation effects dominate at the base (typically for a hard bed) to Coulomb-plastic friction, when friction saturates (typically for weak till). When u₀=0 or q=0, purely plastic friction is recovered.

The merits and physical basis of the different possible friction formulations and nonlinear exponents are still under active debate (Aschwanden et al., 2013; Stearns and van der Veen, 2018; Brondex et al., 2019; Joughin et al., 2019), and all of the above formulations are used in ice-sheet modeling today. However, given the large uncertainty in boundary conditions provided to an ice-sheet model, which include bedrock topography, sediment composition and distribution, basal hydrology and its temporal evolution, etc., it is clear that the use of any formulation will rely on empirical tuning. Also, as noted above, different choices for the friction exponents or threshold values can reduce a given formulation to another. Although modeling studies have shown that all four cases above can produce realistic velocity fields of the present-day ice sheets (e.g. Goelzer et al., 2018; Joughin et al., 2019), it remains to be seen how the choice of friction formulation may impact transient changes in the ice sheet.

For these reasons, we have chosen to implement the friction formulations in the most general way possible in the code, with essentially two free parameters: q as a nonlinear exponent and u₀ as a threshold speed. Meanwhile, c_b is a 2D field that can be set to a constant value or a spatially and/or temporally varying field based on, e.g., whether the ice is frozen to the bed or temperate, on till strength (Bueler and van Pelt, 2015), effective pressure, or other user-defined criteria. As mentioned above, a Picard iteration method is used to solve for basal friction, η_d and the SSA velocity solution until convergence of the velocity solution is reached.

β and c_b are initially defined on aa-nodes. c_b is naturally defined on the grid center, while when β=f(u_b), the velocity components that are defined on ac-nodes must be staggered to the grid center. Once β has been calculated using one of the above formulations, it must be staggered to the ac-nodes for use in the elliptical solver. For purely floating points (i.e., f_g=0 at both bounding aa-nodes) β_ac=0, and for purely grounded points, β_ac is the average of the two neighbors. At the grounding line, Yelmo allows several options to handle staggering. These include simple averaging, taking the upstream value of β, taking the downstream value of β or taking the weighted average based on the grounded fraction of the ac-node.

3.4 Thermodynamics

Thermodynamics in Yelmo is treated in the classical way by solving the following energy conservation equation:

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{\partial T}{\partial t}}={\displaystyle \frac{k}{\mathit{\rho }c}}{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial z^{\mathrm{2}}}}-u{\displaystyle \frac{\partial T}{\partial x}}-v{\displaystyle \frac{\partial T}{\partial y}}-w{\displaystyle \frac{\partial T}{\partial z}}+{\displaystyle \frac{\mathrm{\Phi }}{\mathit{\rho }c}},\end{array}$$

where k and c are the ice thermal conductivity and specific heat capacity, respectively. The evolution of the ice temperature T is driven by vertical diffusion, horizontal and vertical advection, and internal strain heating due to ice shearing, Φ, where

$$\begin{array}{}\text{(27)}& \mathrm{\Phi }=\mathrm{4}\mathit{\eta }{\dot{\mathit{\epsilon }}}^{\mathrm{2}}.\end{array}$$

Horizontal diffusion is assumed to be negligible (Greve and Blatter, 2009). At the air–ice interface (i.e., the ice surface), the ice temperature is prescribed via the input boundary temperature field T_s, limited to a maximum value of T₀=273.15 K. At the base of floating ice, the ice temperature is prescribed to be the expected freezing temperature of seawater as a function of depth (Jenkins, 1991), except near the grounding line, where the temperature is prescribed to be the pressure melting point of ice. At the base of grounded ice, when the ice temperature is below the pressure melting point, the vertical gradient of temperature is prescribed as $\partial T/\partial z=-Q_{\mathrm{geo}}/k$, where the geothermal heat flux (Q_geo) is provided as a boundary field to Yelmo. If the temperature at the ice base reaches the pressure melting point, then the temperature is set to the pressure melting point, and the basal mass balance is diagnosed as (Cuffey and Paterson, 2010):

$$\begin{array}{}\text{(28)}& {\dot{b}}_{\mathrm{g}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }L}}\left(Q_{\mathrm{b}}+k{\displaystyle \frac{\partial T}{\partial z}}{|}_{\mathrm{b}}+Q_{\mathrm{geo}}\right),\end{array}$$

where ${\dot{b}}_{\mathrm{g}}$ is the basal mass balance of grounded ice (negative for melting), L is the latent heat of fusion for ice, Q_b is the basal heat production to due sliding and $\frac{\partial T}{\partial z}{|}_{\mathrm{b}}$ is the ice temperature gradient at the base. Yelmo calculates ${\dot{b}}_{\mathrm{g}}$, which is a model output, in contrast to ${\dot{b}}_{\mathrm{f}}$ (basal mass balance of floating ice), which is prescribed in Yelmo as a boundary condition. Once the ice base is temperate (i.e., at the pressure melting point), it will remain so as long as $W_{\mathrm{til}}-\left(\frac{{\mathit{\rho }}_{\mathrm{w}}}{\mathit{\rho }}{\dot{b}}_{\mathrm{g}}\right)\mathrm{\Delta }t>\mathrm{0}$, where ${\dot{b}}_{\mathrm{g}}$ is used from the previous timestep, and W_til is the water layer thickness in the till beneath the ice sheet. In other words, if it is expected that an energy deficit will result in freeze-on of the total available liquid water at the ice base, then the point is treated as a non-temperate ice point.

Yelmo simulates the evolution of the basal water layer thickness in the till following Bueler and van Pelt (2015):

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\partial W_{\mathrm{til}}}{\partial t}}=-{\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_w}}{\dot{b}}_{\mathrm{g}}-C_{\mathrm{d}},\end{array}$$

where C_d is the prescribed till drainage rate, usually set to C_d=0.001 m yr⁻¹. W_til is limited to the range $\mathrm{0}\le W_{\mathrm{til}}\le W_{\mathrm{til},\max }$, where maximum is usually set to $W_{\mathrm{til},\max }=\mathrm{2}$ m. This approach allows for W_til to maintain consistency with the thermodynamic state of the ice sheet at all times. It does not include horizontal transport, as this could potentially be treated by an external basal hydrology model. It is also possible to disable calculation of W_til inside of Yelmo and instead consider it as a boundary variable. However, given the adaptive timestepping approach used by Yelmo, we have found that updating W_til internally at each timestep helps to avoid artificial oscillations that may develop otherwise when the thermodynamics and basal friction are coupled.

Equation (26) is solved with an implicit method in the vertical direction, while the horizontal advection is solved separately applying an explicit, second-order upwind forward Euler method. This separation allows the energy conservation in the vertical to be solved as a 1D column model. The discretization of vertical diffusion follows the form presented by Hoffman et al. (2018), while the discretization of vertical advection follows a second-order central difference scheme. A given column of grid points consists of temperatures defined on the grid-centers (aa-nodes) and boundary values defined directly at the surface and base of the ice sheet.

4 Timestepping

Yelmo makes use of a predictor–corrector (PC) method combined with adaptive timestepping to balance speed and stability, following the method developed by Cheng et al. (2017). This approach requires calculating the ice thickness twice per timestep, while all other variables can be calculated only once per timestep. Applying a PC method significantly improves the accuracy of the solution compared to a simple forward Euler timestepping method. Furthermore, it facilitates the calculation of a stability metric at each timestep that can be used to evaluate model performance and forms the basis of a robust adaptive timestepping approach (Cheng et al., 2017). A given timestep therefore consists of three parts:

1. Predictor step. The topography component (namely the ice thickness) is predicted using the dynamics, material and thermodynamics solutions from previous timesteps.

2. Update step. Using the predicted topography solution, the dynamics, material and thermodynamics components are then also updated.

3. Corrector step. Using the updated dynamics, material and thermodynamics component solutions, the topography component is finally calculated again, starting from the ice thickness solution of the previous timestep.

In Yelmo, the predictor step is calculated via the second-order Adams–Bashforth (AB) method (Cheng et al., 2017),

$$\begin{array}{}\text{(30)}& H_{n+\mathrm{1}}^{\star }=H_n+\mathrm{\Delta }{t}_n\left[{\mathit{\beta }}_{\mathrm{1}}f\left(H_n,{\stackrel{\mathrm{‾}}{\mathbf{u}}}_n\right)+{\mathit{\beta }}_{\mathrm{2}}f\left(H_{n-\mathrm{1}},{\stackrel{\mathrm{‾}}{\mathbf{u}}}_{n-\mathrm{1}}\right)\right],\end{array}$$

where H^⋆ is the predicted ice thickness, Δt is the timestep, and ${\mathit{\beta }}_{\mathrm{1}}=\mathrm{1}+\frac{{\mathit{\zeta }}_t}{\mathrm{2}}$, ${\mathit{\beta }}_{\mathrm{2}}=-\frac{{\mathit{\zeta }}_t}{\mathrm{2}}$ and ${\mathit{\zeta }}_t=\frac{\mathrm{\Delta }{t}_n}{\mathrm{\Delta }{t}_{n-\mathrm{1}}}$. The labels n, n−1 and n+1 indicate the current, previous and next timestep, respectively. Here, $f\left(H,\stackrel{\mathrm{‾}}{\mathbf{u}}\right)$ is shorthand for $\frac{\partial H}{\partial t}$ as a function of the ice thickness and depth-averaged horizontal velocity field, noting that $\stackrel{\mathrm{‾}}{\mathbf{u}}$ is also a function of the ice thickness, material and potentially thermodynamic state of the ice sheet. For this algorithm, β₁, β₂ and ζ_t are timestep dependent, but the subscript n has been dropped for clarity. Once $H_{n+\mathrm{1}}^{\star }$ has been calculated, the other components are updated, and finally the corrector step is then calculated via the semi-implicit Adams–Moulton (SAM) method (Cheng et al., 2017),

$$\begin{array}{}\text{(31)}& H_{n+\mathrm{1}}=H_n+{\displaystyle \frac{\mathrm{\Delta }{t}_n}{\mathrm{2}}}\left[f\left(H_{n+\mathrm{1}}^{\star },{\stackrel{\mathrm{‾}}{\mathbf{u}}}_{n+\mathrm{1}}\right)+f\left(H_n,{\stackrel{\mathrm{‾}}{\mathbf{u}}}_n\right)\right],\end{array}$$

where H_n+1 is the corrected ice thickness for the next timestep.

For the AB–SAM timestepping method, Cheng et al. (2017) derived the following expression for the leading term of the local truncation error:

$$\begin{array}{}\text{(32)}& {\mathit{\tau }}_{n+\mathrm{1}}={\displaystyle \frac{{\mathit{\zeta }}_t\left(H_{n+\mathrm{1}}-H_{n+\mathrm{1}}^{\star }\right)}{\left(\mathrm{3}{\mathit{\zeta }}_t+\mathrm{3}\right)\mathrm{\Delta }{t}_n}}.\end{array}$$

The local truncation error is valuable for diagnosing the performance of the model and can be used as an indicator of numerical stability. For a small enough timestep, $H_{n+\mathrm{1}}^{\star }$ and H_n+1 should be indistinguishable, and ${\mathit{\tau }}^{n+\mathrm{1}}\sim \mathrm{0}$. However, as the timestep increases, the local truncation error will also increase.

An adaptive timestepping approach based on a proportional-integral (PI) controller method is therefore used to maximize the timestep while maintaining the truncation error below a specified threshold (Cheng et al., 2017; Söderlind and Wang, 2006). Defining the maximum truncation error over all grounded grid points as η=max|τ|, the next timestep is calculated using the so-called PI4.2 controller (Söderlind, 2002):

$$\begin{array}{}\text{(33)}& \mathrm{\Delta }{t}_{n+\mathrm{1}}={\left({\displaystyle \frac{\mathit{ϵ}}{{\mathit{\eta }}^{n+\mathrm{1}}}}\right)}^{\left(k_I+k_p\right)}{\left({\displaystyle \frac{\mathit{ϵ}}{{\mathit{\eta }}^n}}\right)}^{-k_p}\mathrm{\Delta }{t}_n,\end{array}$$

where ϵ is the target tolerance, and $k_I=\mathrm{2}/\mathrm{10}$ and $k_p=\mathrm{1}/\mathrm{10}$ are reasonable control parameters for the second-order AB–SAM timestepping method used here (Söderlind and Wang, 2006). This algorithm ensures that the timestep increases when η<ϵ and decreases when η>ϵ. The use of both ηⁿ⁺¹ and ηⁿ helps to avoid rapid fluctuations in the timestep, which improves model stability and results in a predictable timestep size as a function of the target tolerance.

For practical purposes, the timestep is further treated as follows. The timestep must be larger than a user-prescribed minimum value but smaller than the Courant–Friedrichs–Lewy (CFL) 2D advective limit:

$$\begin{array}{}\text{(34)}& \mathrm{\Delta }{t}_{\mathrm{cfl}}=C_{\mathrm{cfl}}\max {\left|{\displaystyle \frac{\stackrel{\mathrm{‾}}{u}}{\mathrm{\Delta }x}}+{\displaystyle \frac{\stackrel{\mathrm{‾}}{v}}{\mathrm{\Delta }y}}\right|}^{-\mathrm{1}},\end{array}$$

where the maximum is taken over all grid points and C_cfl=1.0. Furthermore, the adaptive timestep is adjusted to ensure that the model stays synchronized with the frequency that Yelmo is being called externally. We found that the latter requirement often results in highly uneven timestepping; e.g., if Yelmo is called with a timestep of Δt_tot=2.0 yr and the first adaptive timestep is determined to be Δt₁=1.9 yr, then the second timestep would likely be Δt₂=0.1 yr. To avoid this possibility and increase stability, the condition is imposed that if any given adaptive timestep is predicted in the range of $\mathrm{0.5}\mathrm{\Delta }{t}_{\mathrm{tot}}<\mathrm{\Delta }t<\mathrm{\Delta }{t}_{\mathrm{tot}}$, then Δt=0.5Δt_tot. In this example, this condition would ensure that $\mathrm{\Delta }{t}_{\mathrm{1}}=\mathrm{\Delta }{t}_{\mathrm{2}}=\mathrm{1.0}$ yr, unless Δt₂ needed to be smaller for stability. Finally, if the maximum local truncation error η is larger than a specified threshold for any given integration, then the integration is discarded and repeated with a progressively smaller timestep until η diminishes and stability is restored, or the timestep reaches the minimum allowed value.

5 Model interface

The Yelmo model interface is designed to be clear and simple but also flexible. In its essence, there are three main model functions: yelmo_init to initialize the model variables, yelmo_update to perform the ice-sheet model calculations for a given timestep and yelmo_end to terminate the Yelmo object (free it from memory).

The first subroutine, yelmo_init, is used to load parameters, initialize variables in memory (i.e., allocate arrays) and, optionally, to initialize the topographic state variables (ice thickness, masks, etc.). No other variables are initialized here in the sense of being populated with data values, which is left to the user. An additional, optional helper function can be used, yelmo_init_state, which populates the remaining model variables in the material, thermodynamics and dynamics components. This initialization step is separated from that of topography because in practice, sometimes boundary variables (e.g., surface temperature) need the surface elevation as input in order to be determined. In contrast, the remaining variables, namely dynamics and thermodynamics, often rely on boundary variables to be initialized. Thus, a typical initialization sequence for a stand-alone ice-sheet model simulation could first call yelmo_init, then load or calculate boundary variables and then call yelmo_init_state to finalize the initialization of all Yelmo variables. After this sequence, the Yelmo state should be consistent with running the model for one timestep with the prescribed boundary conditions and a fixed topography. If the model will be initialized from a restart file, then these data are loaded in each case based on parameter choices.

The next subroutine, yelmo_update, is used to advance the model state to a new timestep. Any modifications to boundary variables are left to the user externally, and Yelmo expects that the boundary conditions are valid for this timestep. The subroutine does not take any arguments to modify the model behavior – rather, all model configuration choices are specified in the parameters of the Yelmo object itself. These are initially loaded from a parameter file in the call to yelmo_init; however, it is possible to modify any parameter values during simulations, allowing for changing the model configuration transiently depending on the experimental setup. An additional optional subroutine, yelmo_update_equil, is available to facilitate equilibration. This routine effectively calls yelmo_update for a specified time window with unchanging boundary conditions and allows for the temporary modification of some key model parameters (such as the maximum allowed adaptive timestep and the maximum allowed SSA velocity).

The last subroutine, yelmo_end, simply removes the Yelmo object from memory (i.e., all domain variables are deallocated). After calling yelmo_end, it is possible to reinitialize the Yelmo object via yelmo_init, for example, in order to test a different grid resolution or other configuration.

There are several input/output routines defined for Yelmo. yelmo_write_init can be used to initialize a NetCDF model output file with the axes of model dimensions defined from the Yelmo object and writing of static fields like domain masks. The writing of model output for individual timesteps is left to the user to maintain flexibility, as most programs require specific fields to be written (examples can be found in the test programs included with the code – see further below). In addition, yelmo_restart_write will create a NetCDF file and write all Yelmo fields as a snapshot, which can be used to restart the model (loading of a restart file can be activated with parameter choices).

As mentioned above, given the object-oriented approach, it is possible to run multiple Yelmo domains in one program. Each domain must be initialized separately via yelmo_init and the variable fields populated with initial values; then separate calls to yelmo_update are needed during timestepping, and finally each object should be terminated at the end of the program via yelmo_end. With this structure, minimum modification of another model, like a global climate model is needed, to incorporate online ice-sheet evolution or to simulate an ensemble of ice sheets in one program. Furthermore, all fields are directly accessible within the main program to facilitate coupling. For example, the 2D array of surface elevation of the topography component of the Yelmo Antarctica domain could be referenced as yelmo_ant%tpo%now%z_srf. While it is clear that the nesting of several containers (derived types) results in a rather long variable reference, it is unambiguous and straightforward to use.

6 Model validation and benchmarks

Yelmo has been tested against several ice-sheet model validation tests and benchmarks in wide use today. These include the Halfar dome experiment (Halfar, 1983; Bueler et al., 2005), the EISMINT1 (Huybrechts et al., 1996) and EISMINT2 (Payne et al., 2000) model intercomparison experiments that test uncoupled and coupled dynamics–thermodynamics, respectively, and MISMIP (Pattyn et al., 2012) for ice-shelf dynamics, among others. By design, many of these experiments allow isolation of specific model features for testing. When the model passes more complex benchmark tests, the simpler experiments are somewhat redundant (if the model passes a coupled thermodynamics–dynamics benchmark, the model should necessarily also be able to pass a dynamics-only benchmark). However, it should be noted that in the process of model development, all tests prove to be extremely useful. The results of all tests will not be reported here, but several are highlighted below to demonstrate that Yelmo performs well.

The Halfar dome experiment, a specific case of the more general Test B of Bueler et al. (2005), tests the ice-sheet model dynamics using the SIA solver alone. This test consists of simulating a radially symmetric ice-sheet with zero mass balance and resting on a flat bed, deforming under gravitational stress. The analytical solution is known at every time, allowing a direct comparison of the simulation to the desired result. The simulation parameters consist of the margin radius and dome elevation, in this case set to the values suggested by Halfar (1983): R₀=21.2132 km and H₀=707.1 m – see Bueler et al. (2005) for further details. Figure 3 shows the root mean square error (RMSE) of the simulation with the analytical result after 200 years for a range of model resolutions. Yelmo demonstrates first-order (p=1.01) numerical convergence with resolution towards the analytical result.

Figure 3 Root mean square error (RMSE) of the Halfar dome experiment after 200 yr simulated by Yelmo compared to the analytical solution versus model resolution. The value of p=1.01 indicates the order of convergence as the resolution increases.

The EISMINT1 moving margin experiment also tests the ice-sheet model dynamics using the SIA solver alone, with an imposed constant rate factor and diagnosed thermodynamics (i.e., thermodynamics do not impact the ice-sheet configuration). Radial steady-state surface mass balance and background surface temperature fields are imposed as boundary conditions. Starting from ice-free conditions, the ice sheet simulated by Yelmo grows to dynamic and thermodynamic equilibrium within 25 and 100 kyr, respectively. The steady-state summit elevation of Yelmo is 3006.6 m compared to the reported range of 2997.5±7.4 m for so-called “Type-I” discretization models like Yelmo (where diffusivity is staggered to the ab-nodes). The basal temperature relative to the pressure melting point (i.e., homologous temperature) at the summit simulated by Yelmo is −13.37 ^∘C, which lies within the EISMINT1 range of $-\mathrm{13.40}\pm \mathrm{0.56}$ ^∘C. These and other relevant statistics are given in Table 1.

Table 1 Yelmo performance in the EISMINT1 moving margin experiment (“Moving”), as well as in the EISMINT2 experiments A and F. Where available, metrics with the ensemble mean and standard deviation from the original benchmark experiments are also provided for comparison.

We also use the EISMINT1 moving margin experiment to demonstrate the capability of the adaptive timestepping approach in Yelmo. By setting the tolerance parameter ϵ, Yelmo automatically adjusts the timestep to maintain the truncation error in ice thickness η around this value. Figure 4a shows the time series of the adaptive timestep used by Yelmo for a 25 kyr simulation for different resolutions. The timestep exhibits oscillations around a mean value, which is typical for such a PI approach (Cheng et al., 2017). When the timestep grows larger, the truncation error increases. This leads to a reduction in the timestep, and the error decreases. Figure 4b shows the mean timestep used by Yelmo over the last 10 kyr of the simulation versus model resolution. Given a tolerance of $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{2}}$, Yelmo's mean timestep is Δt=6.96 yr, Δt=1.59 yr, Δt=0.24 yr and Δt=0.06 yr for resolutions of 50, 25, 10 and 5 km, respectively. As expected, the timestep must be reduced for higher resolutions. These results are in line with those of Cheng et al. (2017) for the same experiment (Δt=12.4 yr for 60 km resolution). It should be noted that the truncation error increases nonlinearly as a function of the timestep, so setting a higher tolerance does not translate directly into a larger timestep.

Figure 4 (a) Adaptive timestepping for the EISMINT1 moving margin experiment. Time series of the timestep used by Yelmo for grid resolutions of 50 km, 25 km and 10 km and a tolerance of $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{2}}$ and (b) the mean adaptive timestep in the time range of 15–25 kyr versus model resolution. Separate lines in (b) show results for different values of the tolerance parameter ϵ.

Next we validate the thermodynamics component first by performing the benchmark experiments Test A and Test B of Kleiner et al. (2015). In contrast to an enthalpy solver, Yelmo uses a temperature solver that assumes all water produced in the ice column drains directly to the bed, and so temperate ice in the vertical column has no water content. In cases where water content of up to 3 % could be present in the basal layers of the column, Yelmo's solver would be inaccurate. Nonetheless, we expect that the temperature solver should be sufficiently accurate to simulate ice sheets on long timescales and large spatial domains. Figure 5 shows the performance of Yelmo's temperature solver for Test A of Kleiner et al. (2015), which simulates a column of ice in a parallel-sided slab with no horizontal advection and no internal strain heating that undergoes warming and subsequent cooling at the surface. In this case, no water content should develop in the vertical column, so a temperature and enthalpy solver should give identical, energy-conserving results. Yelmo's basal melt rate is essentially identical to the analytical solution for this problem, and its transient behavior is robust.

Figure 5 Time series of the (a) basal temperature, (b) basal melt rate and (c) basal water layer corresponding to the thermodynamic benchmark experiment Exp. A (Kleiner et al., 2015). The analytical solution (thick, light-red line) for the basal melt rate is compared to Yelmo results (black lines). Not that where the Yelmo results are not visible, they overlap with the analytical solution.

In contrast, Fig. 6 shows the results of Yelmo's temperature solver for Test B, which simulates a parallel-sided slab on a sloping bed with a prescribed horizontal velocity and strain heating profile in steady state. In this case, water is generated in the basal layers of the ice column (see Fig. 6c); however Yelmo cannot reproduce this solution. Nonetheless, because temperature is limited to the pressure-melting point, the simulated ice temperature profile is in full agreement with the analytical profile. This is true both for a very high resolution case (Δz=0.5 m) and a lower-resolution case (Δz=10 m), which allows us to conclude that the performance of the temperature solver is robust.

Figure 6 Steady-state vertical profiles of (a) enthalpy, (b) temperature and (c) water content for the thermodynamic benchmark experiment Exp. B (Kleiner et al., 2015). The analytical solution (thick, light-red lines) is compared to Yelmo results for a vertical resolution of Δz=0.5 m (n_z=400, black lines) and Δz=10 m (n_z=20, light green lines). The vertical grey line in (b) shows the pressure-melting point as prescribed in this experiment. Note that where the analytical solution is not visible, it overlaps with the Yelmo results.

The EISMINT2 benchmark experiments A and F are useful for testing the thermodynamically coupled ice-sheet model with SIA dynamics like in EISMINT1. The experiments are identical to the EISMINT1 moving margin experiment, except the resolution is doubled (25 km), and the surface temperature is prescribed to be independent of ice thickness. Experiment A prescribes a summit temperature of 238.15 K, while experiment F is 15.00 K colder, which promotes an increase in the region of ice frozen to the bedrock. The statistics for these experiments are listed in Table 1 as well. Figure 7 shows the basal homologous temperature distribution for experiments A and F. Yelmo produces temperature patterns in both experiments, which are consistent with both the benchmark results (Payne et al., 2000) and other more recent models (e.g., Bueler et al., 2007; Hoffman et al., 2018). Axial symmetry, assessed by comparing the basal temperature field with a mirror of itself along the x or y axis, is maintained to a precision of 10⁻² K. This symmetry is not critical to realistic applications, but a lack of at least axial symmetry in this test is often indicative of numerical artifacts. In experiment F, Yelmo produces the so-called “cold spokes”, which have been shown to be related to internal strain heating in regions of steep gradients in ice thickness and are largely numerical in nature (Bueler et al., 2007).

Figure 7 Steady-state, basal homologous temperature (^∘C) distribution after 100 kyr obtained by Yelmo in EISMINT2 (a) experiment A and (b) experiment F. Areas that have reached the pressure-melting point have been shaded grey. The contour lines represent ice thickness at 500 m intervals up to 3500 m.

We also test the capability of the SSA solver and grounding-line treatment by running the MISMIP protocol experiments (Pattyn et al., 2012). Particularly, MISMIP EXP 1 (advance) and EXP 2 (retreat) are useful for testing the reversibility of grounding-line advance, given the bedrock is defined as a linear downward-sloping bed. The rate factor is prescribed according to steps that first decrease, allowing grounding-line advance, then increase back to the original value. According to theory (Weertman, 1974; Schoof, 2007), only one steady-state grounding line position should exist for each step – i.e., the ice sheet should advance and retreat symmetrically without showing hysteresis. It is now well known, however, that ice-sheet models at coarse resolutions (1 km and greater) are unable to capture proper grounding-line migration, even when subgrid parameterizations to mimic higher resolution are applied (Seroussi et al., 2014; Gladstone et al., 2017).

In the MISMIP experiment performed here, the linear, downward-sloping bedrock is defined in the x direction as z_b=720.0–778.5(x∕750.0) with x in kilometers and z_b in meters. The bedrock elevation does not change in the y direction, which extends to ±50 km to allow the simulation of a symmetric ice stream flowing in the positive x direction. The power-law formulation of Eq. (24) is used with the parameter values $q=\mathrm{1}/\mathrm{3}$, u₀=1 m yr⁻¹ and $c_{\mathrm{b}}=\mathrm{3.165176}\times {\mathrm{10}}^{\mathrm{4}}$ Pa. The rate factor is initially prescribed to be $A=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{16}}$ Pa⁻³ yr⁻¹, and the simulation is run for 25 kyr to equilibrate. Next, the rate factor is stepped evenly in log-space every 10 kyr until reaching $A=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{19}}$ Pa⁻³ yr⁻¹, and then the rate factor is increased in the same way until returning to the original value.

Figure 8 shows results for this MISMIP experiment with Yelmo at different grid resolutions, ranging from 20 km down to 2.5 km, and with different treatments of basal friction near the grounding line. When the default model setup is used, with no special treatment at the grounding line, the grounding-line advance is consistent for all resolutions. However, none of the lowest-resolution simulations show grounding line retreat as the rate factor increases again. At a resolution of 5 km, some minor grounding line retreat can be seen, and for 2.5 km, the model is more successful at retreating though it remains 400 km from the target. In contrast, when the basal friction β is scaled at the grounding line by the grounded fraction of the ac-node (f_g,ac), the hysteresis is greatly reduced. The 5 km simulation retreats to within 200 km of the original position, and the 2.5 km simulation retreats to within 100 km of the original position, thus showing convergence to the correct solution with resolution. With this setup, even the 10 and 20 km simulations retreat significantly. In a third case, the basal friction is also linearly scaled to zero as the ice sheet approaches flotation (Leguy et al., 2014; Gladstone et al., 2017). In this case, the hysteresis and differences between different-resolution simulations are further reduced; however, the system also tends to advance much less given all other conditions are the same.

Figure 8 Yelmo performance in the MISMIP bedrock advance and retreat simulations on a linear sloping bed. Panel (a) shows the imposed rate factor A, with 10 kyr steps of decreasing and then increasing values. Panels (b), (c) and (d) show the grounding line position evolution for each of three model configurations, respectively: “Default” is the standard model setup, with no special treatment of friction at or near the grounding line; “Subgrid” uses the grounded fraction at the grounding line to scale the basal friction; and “Scaling” applies both the grounded fraction and imposes a linear reduction in basal friction as the ice sheet approaches flotation. Separate simulations were run for resolutions ranging from 20 km down to 2.5 km.

Yelmo’s Eulerian conservative tracer model is validated with a simulation of ice age in an idealized configuration against the analytical solution presented by Rybak and Huybrechts (2003). In this case, summit-like conditions are imposed, in that horizontal advection is neglected, and the vertical velocity is assumed to decrease linearly with depth. Figure 9 shows the solution with Yelmo as compared to the analytical result. For a nominal vertical resolution of n_z=30 points and single or double precision, the age tracer gives errors in the range of 0.2 %–0.5 % over most of the column of the ice sheet, with higher errors at the base. Increasing the vertical resolution to n_z=50 points decreases the error by an order of magnitude, while using n_z=30 with higher resolution at the base of the ice sheet allows a similar reduction in error with significant computational savings. For Eemian-age ice in such simplified conditions, the latter case gives an uncertainty of less than 1 kyr. It is expected that the error would increase for more realistic 3D domains; however the Eulerian age solver can be used for a first-order estimate of the age–depth profile in the ice sheet (Rybak and Huybrechts, 2003).

Figure 9 Analytical age–depth profile for idealized summit compared to Yelmo Eulerian age tracing. (a) The ice age relative to present day is shown for the analytical solution (thick, light-red line) and Yelmo with a resolution of n_z=30 with a linear vertical axis and compiled at double precision (black line). (b) The associated relative error is given for this case (black line), as well as for a higher resolution of n_z=50 and for n_z=30 compiled at single precision (grey lines), and finally for a resolution of n_z=30 compiled at double precision but with exponentially increasing resolution at the base instead of a linear axis (green line).

7 Antarctica

As further validation of the model's performance, we ran steady-state simulations of the present-day and glacial Antarctic ice sheet. These simulations, run at 32 km resolution, have been deliberately simplified to include the minimum complexity necessary to simulate the ice sheet without additional external components. There was no active isostasy model, and geothermal heat flux was set to 50 mW m⁻² everywhere. Basal friction followed a linear law, where $\mathit{\beta }=\left(c_f{\mathit{\lambda }}_{\mathrm{b}}/u_{\mathrm{0}}\right)\left(\mathit{\rho }gH\right)$ with u₀=100 m yr⁻¹ used as a scaling constant. We prescribed c_f=0.15 (unitless) for most of the domain, except for ad hoc adjustments in specific regions to improve the match with observations. This was additionally scaled by an exponential function of bedrock elevation: ${\mathit{\lambda }}_{\mathrm{b}}=\min [\mathrm{1.0},\exp ((z_{\mathrm{b}}-z_{\mathrm{1}})/(z_{\mathrm{1}}-z_{\mathrm{0}})\left)\right]$, analogous to the approach of (Martin et al., 2011). We set z₁=250 m everywhere and $z_{\mathrm{0}}=-\mathrm{2000}$ m for West Antarctic Ice Sheet (WAIS) regions feeding the Ronne ice shelf and $z_{\mathrm{1}}=-\mathrm{200}$ m elsewhere. Friction was scaled by the grounded fraction at the grounding line, but no additional scaling is applied. The enhancement factor parameters set for these simulations were E_shr=2.5, E_strm=0.7 and E_shlf=0.5. The bedrock topography and initial ice thickness were prescribed from the RTOPO2.1 dataset (Schaffer et al., 2016), after which the model ran for 50 kyr, reaching a steady-state modeled ice distribution.

For the simulation of the present-day state, surface mass balance (SMB) and surface temperature boundary fields were prescribed from a RACMO2.3 simulation driven by ERA-Interim data and averaged over 1981–2010 (van Wessem et al., 2018). The ice-shelf basal mass balance was set to a spatially constant value of −0.2 m a⁻¹ where floating ice exists today and to −2.0 m a⁻¹ elsewhere. Figures 10 and 11 show a comparison of the simulations with the observed topography (RTOPO2.1) and the present-day observed velocity (Rignot et al., 2011). With this relatively simple model setup, it is possible to obtain reasonable agreement with observations. The root mean square errors (RMSEs) in ice thickness, velocity and log(velocity) are 320 m, 270 m yr⁻¹ and 1.9 log(m yr⁻¹), respectively, which fall in the range of other models in the initMIP-Antarctica intercomparison project (Seroussi et al., 2019). The simulated ice sheet is thinner than the observed ice sheet over large parts of East Antarctica, with a broad positive bias near the South Pole (Fig. 11). The margins of the ice sheet are the most difficult to match, in particular, the grounding-line positions of the large ice shelves, leading to larger biases in these regions. This pattern is quite consistent with other studies (e.g., Martin et al., 2011; Quiquet et al., 2018; Albrecht et al., 2020). Overall, the dome configuration, slow deformational speeds and even most ice streams as they penetrate inland are well represented by the model (Figs. 10 and 11).

Figure 10 (a) Antarctica present-day ice sheet configuration and surface velocities from observations, compared to (b) a steady-state simulation with Yelmo. (c) In contrast, Antarctica glacial configuration and surface velocities simulated by Yelmo. Simulations were performed at 32 km resolution. The colors show surface velocity in meters per year, and the dark grey contours show surface elevation in 500 m intervals (thick lines correspond to 1000, 2000 and 3000 m above sea level). The black line shows the grounding-line position.

Figure 11 (a) Simulated present-day ice-thickness minus observations and (b) simulated versus observed ice surface velocity. The colors in (a) show the ice thickness difference in meters, and the dark grey contours show surface elevation in 500 m intervals (thick lines correspond to 1000, 2000 and 3000 m above sea level). The black line shows the grounding-line position. In (b), the dark red line indicates a perfect correlation between model and observations.

We use the same setup with modified boundary conditions to simulate a configuration resembling that of a deep glacial period like the Last Glacial Maximum. The surface temperature was set to 10 ^∘C colder, and the present-day SMB was maintained, except that points with a low or negative SMB were prescribed with a minimum value of 0.1 m a⁻¹. The ice shelf basal mass balance was set to a spatially constant value of 0 m a⁻¹, and sea level was lowered by 120 m. In this case, the grounded ice sheet advances until the continental shelf break and thickens inland (Fig. 10). A similar structure of ice streams can be seen, due to the topographic dependence of β, but their speed is greatly reduced compared to those of the present-day simulation. We do not expect this configuration to be realistic, given that isostasy plays no role and a present-day-like SMB has been imposed. However, this test demonstrates that Yelmo is capable of resolving continental-scale changes in the ice sheet configuration in a plausible way.

8 Conclusions and future work

We have described the features and physics of the hybrid ice-sheet–shelf model Yelmo. Yelmo includes the physics to simulate continental-scale ice sheets and floating ice shelves using “shallow” approximations of the ice dynamics. The fully coupled thermomechanical ice-sheet model has been validated against several benchmark tests and has been shown to simulate the dynamic configuration of the Antarctic ice sheet well.

Yelmo is expected to be useful for long-timescale simulations and/or ensembles. It is particularly suited for easy coupling with other models. For example, the simulation of multiple ice-sheet domains with independent parameter configurations coupled to a global climate model can be achieved in a simple and straightforward way. Also, given that the subroutines representing the physics of the model have been isolated from the “model accounting”, it is possible to test individual model components in different contexts easily. This should facilitate future model development and comparison of different methods.

The model framework has been designed to facilitate the incorporation of new and different physics. Thus, this initial release of Yelmo lays the foundation for several future developments. These may include more advanced calving and basal friction schemes, as well as improved treatment of the grounding line. We also plan to transition to an enthalpy-based thermodynamics solver; however this will require an adaptive vertical axis to be able to map the height of transition between temperate and cold ice accurately. We also plan to implement a variationally derived “depth-integrated-viscosity approximation” solver (following, e.g., Goldberg, 2011; Pollard and DeConto, 2012; Lipscomb et al., 2019) in the near future.