Model description paper 24 Jun 2020
Model description paper  24 Jun 2020
Description and validation of the icesheet model Yelmo (version 1.0)
 ^{1}Complutense University of Madrid, Madrid, Spain
 ^{2}Geosciences Institute CSICUCM, Madrid, Spain
 ^{3}Potsdam Institute for Climate Impact Research, Potsdam, Germany
 ^{4}Institute for Marine and Atmospheric research Utrecht, Utrecht University, Utrecht, the Netherlands
 ^{5}Laboratoire de Glaciologie, Université Libre de Bruxelles, Brussels, Belgium
 ^{6}Institute of Low Temperature Science, Hokkaido University, Sapporo, Japan
 ^{7}Arctic Research Center, Hokkaido University, Sapporo, Japan
 ^{8}Institute for Geosciences and Environmental Research, Grenoble, France
 ^{1}Complutense University of Madrid, Madrid, Spain
 ^{2}Geosciences Institute CSICUCM, Madrid, Spain
 ^{3}Potsdam Institute for Climate Impact Research, Potsdam, Germany
 ^{4}Institute for Marine and Atmospheric research Utrecht, Utrecht University, Utrecht, the Netherlands
 ^{5}Laboratoire de Glaciologie, Université Libre de Bruxelles, Brussels, Belgium
 ^{6}Institute of Low Temperature Science, Hokkaido University, Sapporo, Japan
 ^{7}Arctic Research Center, Hokkaido University, Sapporo, Japan
 ^{8}Institute for Geosciences and Environmental Research, Grenoble, France
Correspondence: Alexander Robinson (robinson@ucm.es)
Hide author detailsCorrespondence: Alexander Robinson (robinson@ucm.es)
We describe the physics and features of the icesheet model Yelmo, an opensource 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 shallowice and shallowshelf/shelfystream approximations, which makes Yelmo an apt choice for studying a wide variety of problems. Yelmo's main innovations lie in its flexible and userfriendly infrastructure, which promotes portability and facilitates longterm development. In particular, all physics subroutines have been designed to be selfcontained, so that they can be easily ported from Yelmo to other models, or easily replaced by improved or alternate methods in the future. Furthermore, hardcoded 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 icesheet benchmark tests, and we illustrate Yelmo's performance for the Antarctic ice sheet.
The field of continentalscale icesheet 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 longtimescale paleo simulations in part because they are fast but also because they are relatively simple in design, usually relying on lowtech solutions to numerical problems. Most of these models were designed before the era of the highperformance computing cluster, which made it challenging to build models otherwise.
Nowadays, a large number of icesheet 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 higherorder 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 finiteelement/finitevolume 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 icesheet model Yelmo^{1}, which is intended to provide access to complex and robust model physics through an intuitive model design. It is a hybrid icedynamics model that is easy to use and configure. We expect that Yelmo will be useful for longtimescale simulations of the continental ice sheets, coupled climate–icesheet 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 presentday and glacial Antarctic ice sheet (Sect. 7), followed by the conclusions (Sect. 8).
Yelmo has been inspired and largely derived from classical icesheet 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 standalone executable. The strict application of this philosophy drove many design choices and allowed us to develop a robust icesheet 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 “lowtech” 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 communitystandard 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 objectoriented 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 Yelmoobject 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, nonportable 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 selfcontained 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 componentlevel API.
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 icesheet 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 icesheet models, Yelmo uses the Arakawa Cgrid 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 “aanodes”. Velocity components and gradients are calculated on cell edges (“acnodes”), and scalar coefficients, like diffusivity in the SIA approach, are calculated on cell corners (“abnodes”). 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 acnodes), but it also results in greater numerical stability of the model (Macayeal, 1997).
Yelmo requires an evenly spaced Cartesian grid in the horizontal direction, while the vertical component follows a classic sigmacoordinate system (Greve and Blatter, 2009). The vertical axis ζ represents the relative height within the ice sheet, running from ζ=0 at the icesheet base to ζ=1 at the icesheet surface:
where z is the elevation relative to presentday 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 sigmacoordinate 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 acnodes in the vertical and aanodes in the horizontal, while horizontal velocities are calculated on acnodes in the horizontal and aanodes in the vertical. Boundary conditions in a vertical column are applied directly at the ice base and ice surface, which correspond to acnodes (see Fig. 2).
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:
where H is the ice thickness, $\stackrel{\mathrm{\u203e}}{\mathbf{u}}=\left(\stackrel{\mathrm{\u203e}}{u},\stackrel{\mathrm{\u203e}}{v}\right)$ is the depthaveraged 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 depthaveraged 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 icecovered 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 icecovered grid point that has an icefree neighbor, the reference ice thickness of the margin point (H_{ref}) is defined as the minimum thickness of the direct, icecovered 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}}=\mathrm{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):
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 iceshelf growth and thus groundingline advance in transient, glacial simulations but has little impact on the steadystate 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 aanode:
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 aanodes, 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 acnodes, the values of f_{g,acx} and f_{g,acy} are determined by linearly interpolating H_{g} from the two bounding aanodes. When both bounding aanodes are positive ${f}_{\mathrm{g},\mathrm{ac}}=\mathrm{1}$, and when both are negative ${f}_{\mathrm{g},\mathrm{ac}}=\mathrm{0}$. When one aanode is positive (${H}_{{\mathrm{g}}_{+}}$) and one aanode is negative (${H}_{{\mathrm{g}}_{}}$), the grounded fraction on the acnode is determined from linear interpolation:
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
This approach ensures that the surface elevation solution is consistent with the Archimedes flotation criterion on aanodes.
The remaining tasks of the topography component are to diagnose other useful topographic characteristics, such as surface and ice thickness gradients (on acnodes) 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
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}}$),
and the strain rate tensor itself, following index notation, is
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:
Here R is the ideal gas constant; A_{0} and Q_{a} are the temperaturedependent rate factor coefficient and activation energy, respectively (see Greve and Blatter, 2009). E_{f} is a socalled 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:
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,
is solved with a secondorder, 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):
Basal freezeon 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 firstorder 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 agedependent enhancement factor (Greve, 1997b).
3.3 Dynamics
The Yelmo dynamics component is currently representative of a “hybrid” class of icesheet model, treating different modes of ice deformation via a combination of the simplifying shallowice and shallowshelf 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}):
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{\u203e}}{u}={\stackrel{\mathrm{\u203e}}{u}}_{\mathrm{i}}+{u}_{\mathrm{b}}$. To calculate u_{i} and v_{i}, Yelmo uses zeroorder SIA equations:
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{\u25bf}{z}_{\mathrm{s}}$ is the driving stress. In the horizontal plane, the term in brackets is calculated on the abnodes for stability and improved mass conservation (Huybrechts et al., 1996), and then it is staggered onto the acnodes 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 (aanodes). 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):
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 userdefined formulation (power law, regularized Coulomb, etc.), depending on the context (see basal friction description below for details). The depthaveraged (2D) effective viscosity, which is only used for solving the SSA dynamics, is defined as
$\stackrel{\mathrm{\u203e}}{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:
In Yelmo, ${\dot{\mathit{\epsilon}}}_{\mathrm{d}}$ is only used for calculating ${\stackrel{\mathrm{\u203e}}{\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
The depth of the seawater up to the flotation depth, H_{o}, is defined as ${H}_{\mathrm{o}}=\mathrm{min}\left(\mathrm{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 5dev (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^{2} relative error norm (Gagliardini et al., 2013):
where (u_{1},v_{1}) and (u_{0},v_{0}) 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,
The vertical velocity is naturally defined on the acnodes 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
where β represents the basal friction coefficient (with units of Pa yr m^{−1}), 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
Thus in all formulations implemented in Yelmo, the term f(u_{b}) has units of years per meter (yr m^{−1}), 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), powerlaw (e.g. Pattyn, 2017), pseudoplastic powerlaw (e.g. Aschwanden et al., 2013) or regularizedCoulomb (Joughin et al., 2019) formulation. The linear and powerlaw formulations are contained within the pseudoplastic powerlaw formulation, so only the latter and the regularizedCoulomb formulation are needed to represent all four cases.
The pseudoplastic powerlaw formulation (Schoof, 2010; Aschwanden et al., 2013) is
and thus $\mathit{\beta}={c}_{\mathrm{b}}{{u}_{\mathrm{0}}}^{q}{\left{\mathbf{u}}_{\mathrm{b}}\right}^{q\mathrm{1}}$, with the pseudoplastic exponent $q\in \left(\mathrm{0},\mathrm{1}\right)$ and threshold speed u_{0}. This expression results in purely plastic friction for q=0, linear friction for q=1 and powerlaw friction for $\mathrm{0}<q<\mathrm{1}$. With q=1 and u_{0}=1, for example, β=c_{b}, and friction scales linearly with velocity. To obtain the powerlaw formulation used in the original MISMIP experiments (Pattyn et al., 2012), the following parameter values can be prescribed: $q=\mathrm{1}/\mathrm{3}$, u_{0}=1 m yr^{−1} 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
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_{0} 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 Coulombplastic friction, when friction saturates (typically for weak till). When u_{0}=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 icesheet modeling today. However, given the large uncertainty in boundary conditions provided to an icesheet 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 presentday 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_{0} 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 userdefined 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 aanodes. c_{b} is naturally defined on the grid center, while when β=f(u_{b}), the velocity components that are defined on acnodes must be staggered to the grid center. Once β has been calculated using one of the above formulations, it must be staggered to the acnodes for use in the elliptical solver. For purely floating points (i.e., f_{g}=0 at both bounding aanodes) β_{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 acnode.
3.4 Thermodynamics
Thermodynamics in Yelmo is treated in the classical way by solving the following energy conservation equation:
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
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_{0}=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):
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 freezeon of the total available liquid water at the ice base, then the point is treated as a nontemperate ice point.
Yelmo simulates the evolution of the basal water layer thickness in the till following Bueler and van Pelt (2015):
where C_{d} is the prescribed till drainage rate, usually set to C_{d}=0.001 m yr^{−1}. W_{til} is limited to the range $\mathrm{0}\le {W}_{\mathrm{til}}\le {W}_{\mathrm{til},\mathrm{max}}$, where maximum is usually set to ${W}_{\mathrm{til},\mathrm{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, secondorder 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 secondorder central difference scheme. A given column of grid points consists of temperatures defined on the gridcenters (aanodes) and boundary values defined directly at the surface and base of the ice sheet.
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:

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

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

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 secondorder Adams–Bashforth (AB) method (Cheng et al., 2017),
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{\u203e}}{\mathbf{u}}\right)$ is shorthand for $\frac{\partial H}{\partial t}$ as a function of the ice thickness and depthaveraged horizontal velocity field, noting that $\stackrel{\mathrm{\u203e}}{\mathbf{u}}$ is also a function of the ice thickness, material and potentially thermodynamic state of the ice sheet. For this algorithm, β_{1}, β_{2} 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 semiimplicit Adams–Moulton (SAM) method (Cheng et al., 2017),
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:
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 proportionalintegral (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 socalled PI4.2 controller (Söderlind, 2002):
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 secondorder 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 η^{n+1} and η^{n} 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 userprescribed minimum value but smaller than the Courant–Friedrichs–Lewy (CFL) 2D advective limit:
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}=1.9 yr, then the second timestep would likely be Δt_{2}=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_{2} 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.
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 icesheet 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 standalone icesheet 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 objectoriented 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 icesheet 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.
Yelmo has been tested against several icesheet 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 iceshelf 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 dynamicsonly 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 icesheet model dynamics using the SIA solver alone. This test consists of simulating a radially symmetric icesheet 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_{0}=21.2132 km and H_{0}=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 firstorder (p=1.01) numerical convergence with resolution towards the analytical result.
The EISMINT1 moving margin experiment also tests the icesheet model dynamics using the SIA solver alone, with an imposed constant rate factor and diagnosed thermodynamics (i.e., thermodynamics do not impact the icesheet configuration). Radial steadystate surface mass balance and background surface temperature fields are imposed as boundary conditions. Starting from icefree conditions, the ice sheet simulated by Yelmo grows to dynamic and thermodynamic equilibrium within 25 and 100 kyr, respectively. The steadystate summit elevation of Yelmo is 3006.6 m compared to the reported range of 2997.5±7.4 m for socalled “TypeI” discretization models like Yelmo (where diffusivity is staggered to the abnodes). 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.
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{\u03f5}={\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.
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 parallelsided 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, energyconserving results. Yelmo's basal melt rate is essentially identical to the analytical solution for this problem, and its transient behavior is robust.
In contrast, Fig. 6 shows the results of Yelmo's temperature solver for Test B, which simulates a parallelsided 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 pressuremelting 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 lowerresolution case (Δz=10 m), which allows us to conclude that the performance of the temperature solver is robust.
The EISMINT2 benchmark experiments A and F are useful for testing the thermodynamically coupled icesheet 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^{−2} 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 socalled “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).
We also test the capability of the SSA solver and groundingline 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 groundingline advance, given the bedrock is defined as a linear downwardsloping bed. The rate factor is prescribed according to steps that first decrease, allowing groundingline advance, then increase back to the original value. According to theory (Weertman, 1974; Schoof, 2007), only one steadystate 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 icesheet models at coarse resolutions (1 km and greater) are unable to capture proper groundingline 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, downwardsloping 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 powerlaw formulation of Eq. (24) is used with the parameter values $q=\mathrm{1}/\mathrm{3}$, u_{0}=1 m yr^{−1} 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^{−3} yr^{−1}, and the simulation is run for 25 kyr to equilibrate. Next, the rate factor is stepped evenly in logspace every 10 kyr until reaching $A=\mathrm{1}\times {\mathrm{10}}^{\mathrm{19}}$ Pa^{−3} yr^{−1}, 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 groundingline advance is consistent for all resolutions. However, none of the lowestresolution 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 acnode (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 differentresolution simulations are further reduced; however, the system also tends to advance much less given all other conditions are the same.
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, summitlike 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 Eemianage 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 firstorder estimate of the age–depth profile in the ice sheet (Rybak and Huybrechts, 2003).
As further validation of the model's performance, we ran steadystate simulations of the presentday 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^{−2} 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_{0}=100 m yr^{−1} 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}}=\mathrm{min}[\mathrm{1.0},\mathrm{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_{1}=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 steadystate modeled ice distribution.
For the simulation of the presentday state, surface mass balance (SMB) and surface temperature boundary fields were prescribed from a RACMO2.3 simulation driven by ERAInterim data and averaged over 1981–2010 (van Wessem et al., 2018). The iceshelf basal mass balance was set to a spatially constant value of −0.2 m a^{−1} where floating ice exists today and to −2.0 m a^{−1} elsewhere. Figures 10 and 11 show a comparison of the simulations with the observed topography (RTOPO2.1) and the presentday 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^{−1} and 1.9 log(m yr^{−1}), respectively, which fall in the range of other models in the initMIPAntarctica 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 groundingline 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).
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 presentday SMB was maintained, except that points with a low or negative SMB were prescribed with a minimum value of 0.1 m a^{−1}. The ice shelf basal mass balance was set to a spatially constant value of 0 m a^{−1}, 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 presentday simulation. We do not expect this configuration to be realistic, given that isostasy plays no role and a presentdaylike SMB has been imposed. However, this test demonstrates that Yelmo is capable of resolving continentalscale changes in the ice sheet configuration in a plausible way.
We have described the features and physics of the hybrid icesheet–shelf model Yelmo. Yelmo includes the physics to simulate continentalscale ice sheets and floating ice shelves using “shallow” approximations of the ice dynamics. The fully coupled thermomechanical icesheet 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 longtimescale simulations and/or ensembles. It is particularly suited for easy coupling with other models. For example, the simulation of multiple icesheet 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 enthalpybased 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 “depthintegratedviscosity approximation” solver (following, e.g., Goldberg, 2011; Pollard and DeConto, 2012; Lipscomb et al., 2019) in the near future.
Yelmo is maintained as a git repository hosted at https://github.com/palmaice/yelmo (last access: 2 May 2020; palmaice, 2020) under the license GPL3.0. Model documentation can be found at https://palmaice.github.io/yelmodocs/ (last access: 2 May 2020; Yelmodocs, 2020). The exact version of the model, along with the necessary input data, used to produce the results used in this paper is archived on Zenodo (https://doi.org/10.5281/zenodo.3782650, Robinson et al., 2020) and has been tagged in the repository as v1.02.
AR, JAS and MM conceived the model design and features. AR wrote the model code with contributions from the remaining authors. All authors contributed to the model testing and writing the paper.
Heiko Goelzer is a member of the editorial board of the journal.
We would like to thank Mahé Perrette, Christophe Dumas, Gunter Leguy and Bill Lipscomb for valuable discussions about model design that improved Yelmo, Akira Nishida for help with Lis, and Ilaria Tabone and Javier Blasco for extensive model testing at intermediate development points. We are also grateful to the reviewers for helpful comments.
This research has been supported by the Spanish Ministry of Science and Innovation project RIMA (grant no. CGL201785975R). Alexander Robinson was funded by the Ramón y Cajal Programme of the Spanish Ministry for Science, Innovation and Universities (grant no. RYC201620587). Heiko Goelzer has received funding from the program of the Netherlands Earth System Science Centre (NESSC), financially supported by the Dutch Ministry of Education, Culture and Science (OCW) (grant no. 024.002.001). Ralf Greve was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (grant nos. JP16H02224, JP17H06104 and JP17H06323), by the Japanese Ministry of Education, Culture, Sports, Science and Technology (MEXT) through the Arctic Challenge for Sustainability (ArCS) project, and through the Arctic Challenge for Sustainability (ArCS) project (program grant no. JPMXD1300000000).
This paper was edited by Steven Phipps and reviewed by Fuyuki Saito and one anonymous referee.
Albrecht, T., Martin, M., Haseloff, M., Winkelmann, R., and Levermann, A.: Parameterization for subgridscale motion of iceshelf calving fronts, The Cryosphere, 5, 35–44, https://doi.org/10.5194/tc5352011, 2011. a
Albrecht, T., Winkelmann, R., and Levermann, A.: Glacialcycle simulations of the Antarctic Ice Sheet with the Parallel Ice Sheet Model (PISM) – Part 2: Parameter ensemble analysis, The Cryosphere, 14, 633–656, https://doi.org/10.5194/tc146332020, 2020. a
Arakawa, A. and Lamb, V. R.: Computational Design of the Basic Dynamical Processes of the UCLA General Circulation Model, in: General Circulation Models of the Atmosphere, edited by: Chang, J., Vol. 17 of Methods in Computational Physics: Advances in Research and Applications, 173–265, Elsevier, https://doi.org/10.1016/B9780124608177.500094, 1977. a
Aschwanden, A., Aðalgeirsdéttir, G., and Khroulev, C.: Hindcasting to measure ice sheet model sensitivity to initial states, The Cryosphere, 7, 1083–1093, https://doi.org/10.5194/tc710832013, 2013. a, b, c
Brondex, J., GilletChaulet, F., and Gagliardini, O.: Sensitivity of centennial mass loss projections of the Amundsen basin to the friction law, The Cryosphere, 13, 177–195, https://doi.org/10.5194/tc131772019, 2019. a, b
Bueler, E. and Brown, J.: Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model, J. Geophys. Res., 114, F03008, https://doi.org/10.1029/2008JF001179, 2009. a, b
Bueler, E. and van Pelt, W.: Massconserving subglacial hydrology in the Parallel Ice Sheet Model version 0.6, Geosci. Model Dev., 8, 1613–1635, https://doi.org/10.5194/gmd816132015, 2015. a, b
Bueler, E., Lingle, C. S., KallenBrown, J. A., Covey, D. N., and Bowman, L. N.: Exact solutions and verification of numerical model for isothermal ice sheets, J. Glaciol., 51, 291–306, https://doi.org/10.3189/172756505781829449, 2005. a, b, c
Bueler, E., Brown, J., and Lingle, C.: Exact solutions to the thermocoupled shallow ice approximation: effective tools for verification, J. Glaciol., 53, 499–516, https://doi.org/10.3189/002214307783258396, 2007. a, b
Cheng, G., Lötstedt, P., and von Sydow, L.: Accurate and stable time stepping in ice sheet modeling, J. Comput. Phys., 329, 29–47, https://doi.org/10.1016/j.jcp.2016.10.060, 2017. a, b, c, d, e, f, g, h
Cornford, S. L., Martin, D. F., Graves, D. T., Ranken, D. F., Le Brocq, A. M., Gladstone, R. M., Payne, A. J., Ng, E. G., and Lipscomb, W. H.: Adaptive mesh, finite volume modeling of marine ice sheets, J. Comput. Phys., 232, 529–549, https://doi.org/10.1016/j.jcp.2012.08.037, 2013. a, b
Cuffey, K. M. and Paterson, W. S. B.: The physics of glaciers, Academic Press, Burlington, MA, USA, 2010. a, b
Gagliardini, O., Zwinger, T., GilletChaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a newgeneration ice sheet model, Geosci. Model Dev., 6, 1299–1318, https://doi.org/10.5194/gmd612992013, 2013. a, b, c
Gladstone, R. M., Warner, R. C., GaltonFenzi, B. K., Gagliardini, O., Zwinger, T., and Greve, R.: Marine ice sheet model performance depends on basal sliding physics and subshelf melting, The Cryosphere, 11, 319–329, https://doi.org/10.5194/tc113192017, 2017. a, b
Glen, J. W.: The creep of polycrystalline ice, P. Roy. Soc. A Math. Phy., 228, 519–538, 1955. a, b
Goelzer, H., Nowicki, S., Edwards, T., Beckley, M., AbeOuchi, A., Aschwanden, A., Calov, R., Gagliardini, O., GilletChaulet, F., Golledge, N. R., Gregory, J., Greve, R., Humbert, A., Huybrechts, P., Kennedy, J. H., Larour, E., Lipscomb, W. H., Le clec'h, S., Lee, V., Morlighem, M., Pattyn, F., Payne, A. J., Rodehacke, C., Rückamp, M., Saito, F., Schlegel, N., Seroussi, H., Shepherd, A., Sun, S., van de Wal, R., and Ziemen, F. A.: Design and results of the ice sheet model initialisation experiments initMIPGreenland: an ISMIP6 intercomparison, The Cryosphere, 12, 1433–1460, https://doi.org/10.5194/tc1214332018, 2018. a
Goldberg, D. N.: A variationally derived and depthintegrated approximation to a and higherorder glaciological flow model, J. Glaciol., 57, 157–169, https://doi.org/10.3189/002214311795306763, 2011. a, b
Greve, R.: A continuum–mechanical formulation for shallow polythermal ice sheets, Philos. T. R. Soc. SA, 355, 921–974, https://doi.org/10.1098/rsta.1997.0050, 1997a. a, b
Greve, R.: Application of a Polythermal ThreeDimensional Ice Sheet Model to the Greenland Ice Sheet: Response to SteadyState and Transient Climate Scenarios, J. Climate, 10, 901–918, https://doi.org/10.1175/15200442(1997)010<0901:AOAPTD>2.0.CO;2, 1997b. a
Greve, R.: Geothermal heat flux distribution for the Greenland ice sheet, derived by combining a global representation and information from deep ice cores, Polar Data J., 3, 22–36, https://doi.org/10.20575/00000006, 2019. a, b
Greve, R. and Blatter, H.: Dynamics of Ice Sheets and Glaciers, SpringerVerlag, Berlin, 2009. a, b, c, d, e, f, g, h, i, j, k
Greve, R., Wang, Y., and Mügge, B.: Comparison of numerical schemes for the solution of the advective age equation in ice sheets, Ann. Glaciol., 35, 487–494, https://doi.org/10.3189/172756402781817112, 2002. a
Halfar, P.: On the Dynamics and of the Ice and Sheets 2, J. Geophys. Res., 88, 6043–6051, https://doi.org/10.1029/JC088iC10p06043, 1983. a, b
Hoffman, M. J., Perego, M., Price, S. F., Lipscomb, W. H., Zhang, T., Jacobsen, D., Tezaur, I., Salinger, A. G., Tuminaro, R., and Bertagna, L.: MPASAlbany Land Ice (MALI): a variableresolution ice sheet model for Earth system modeling using Voronoi grids, Geosci. Model Dev., 11, 3747–3780, https://doi.org/10.5194/gmd1137472018, 2018. a, b, c, d, e
Huybrechts, P., Instituut, G., and Oerlemans, J.: Evolution of the East Antarctic Ice Sheet: A Numerical Study of ThermoMechanical Response Patterns With Changing Climate, Ann. Glaciol., 11, 52–59, https://doi.org/10.3189/S0260305500006327, 1988. a
Huybrechts, P., Payne, T., AbeOuchi, A., Calov, R., Fabre, A., Fastook, J. L., Greve, R., Hindmarsh, R. C., Hoydal, O., Jóhannesson, T., MacAyeal, D. R., Marsiat, I., Ritz, C., Verbitsky, M. Y., Waddington, E. D., and Warner, R.: The EISMINT benchmarks for testing icesheet models, Ann. Glaciol., 23, 1–12, https://doi.org/10.3189/S0260305500013197, 1996. a, b
Jenkins, A.: A OneDimensional and Model of Ice and ShelfOcean Interaction, J. Geophys. Res., 96, 20671–20677, https://doi.org/10.1029/91JC01842, 1991. a
Joughin, I., Smith, B. E., and Schoof, C. G.: Regularized Coulomb Friction Laws for Ice Sheet Sliding: Application to Pine Island Glacier, Antarctica, Geophys. Res. Lett., 46, 4764–4771, https://doi.org/10.1029/2019gl082526, 2019. a, b, c, d
Kleiner, T., Rückamp, M., Bondzio, J. H., and Humbert, A.: Enthalpy benchmark experiments for numerical ice sheet models, The Cryosphere, 9, 217–228, https://doi.org/10.5194/tc92172015, 2015. a, b, c, d
Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res., 117, F01022, https://doi.org/10.1029/2011JF002140, 2012. a, b
Leguy, G. R., AsayDavis, X. S., and Lipscomb, W. H.: Parameterization of basal friction near grounding lines in a onedimensional ice sheet model, The Cryosphere, 8, 1239–1259, https://doi.org/10.5194/tc812392014, 2014. a
Lipscomb, W. H., Price, S. F., Hoffman, M. J., Leguy, G. R., Bennett, A. R., Bradley, S. L., Evans, K. J., Fyke, J. G., Kennedy, J. H., Perego, M., Ranken, D. M., Sacks, W. J., Salinger, A. G., Vargo, L. J., and Worley, P. H.: Description and evaluation of the Community Ice Sheet Model (CISM) v2.1, Geosci. Model Dev., 12, 387–424, https://doi.org/10.5194/gmd123872019, 2019. a, b, c, d
Ma, Y., Gagliardini, O., Ritz, C., GilletChaulet, F., Durand, G., and Montagnat, M.: Enhancement factors for grounded ice and iceshelf both inferred from an anisotropic ice flow model, J. Glaciol., 56, 805–812, https://doi.org/10.3189/002214310794457209, 2010. a
Macayeal, D.: EISMINT: Lessons in IceSheet Modeling, available at: http://geosci.uchicago.edu/pdfs/macayeal/lessons.pdf (last access: 2 May 2020), 1997. a, b
Martin, M. A., Winkelmann, R., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISMPIK) – Part 2: Dynamic equilibrium simulation of the Antarctic ice sheet, The Cryosphere, 5, 727–740, https://doi.org/10.5194/tc57272011, 2011. a, b
Nishida, A.: Experience in Developing an Open Source Scalable Software Infrastructure in Japan, in: Computational Science and Its Applications – ICCSA 2010, edited by: Taniar, D., Gervasi, O., Murgante, B., Pardede, E., and Apduhan, B. O., 448–462, Springer Berlin Heidelberg, Berlin, Heidelberg, 2010. a
palmaice: Yelmo, GitHub, available at: https://github.com/palmaice/yelmo, last access: 2 May 2020. a
Pattyn, F.: Sealevel response to melting of Antarctic ice shelves on multicentennial timescales with the fast Elementary Thermomechanical Ice Sheet model (f.ETISh v1.0), The Cryosphere, 11, 1851–1878, https://doi.org/10.5194/tc1118512017, 2017. a, b, c
Pattyn, F., Schoof, C., Perichon, L., Hindmarsh, R. C. A., Bueler, E., de Fleurian, B., Durand, G., Gagliardini, O., Gladstone, R., Goldberg, D., Gudmundsson, G. H., Huybrechts, P., Lee, V., Nick, F. M., Payne, A. J., Pollard, D., Rybak, O., Saito, F., and Vieli, A.: Results of the Marine Ice Sheet Model Intercomparison Project, MISMIP, The Cryosphere, 6, 573–588, https://doi.org/10.5194/tc65732012, 2012. a, b, c
Payne, A. J., Huybrechts, P., AbeOuchi, A., Calov, R., Fastook, J. L., Greve, R., Marshall, S. J., Marsiat, I., Ritz, C., Tarasov, L., and Thomassen, M. P.: Results from the EISMINT model intercomparison: The effects of thermomechanical coupling, J. Glaciol., 46, 227–238, https://doi.org/10.3189/172756500781832891, 2000. a, b
Peyaud, V., Ritz, C., and Krinner, G.: Modelling the Early Weichselian Eurasian Ice Sheets: role of ice shelves and influence of icedammed lakes, Clim. Past, 3, 375–386, https://doi.org/10.5194/cp33752007, 2007. a, b
Pollard, D. and DeConto, R. M.: Description of a hybrid ice sheetshelf model, and application to Antarctica, Geosci. Model Dev., 5, 1273–1295, https://doi.org/10.5194/gmd512732012, 2012. a, b, c, d
Quiquet, A., Dumas, C., Ritz, C., Peyaud, V., and Roche, D. M.: The GRISLI ice sheet model (version 2.0): calibration and validation for multimillennial changes of the Antarctic ice sheet, Geosci. Model Dev., 11, 5003–5025, https://doi.org/10.5194/gmd1150032018, 2018. a, b, c, d
Rignot, E., Mouginot, J., and Scheuchl, B.: Ice Flow of the Antarctic Ice Sheet, Science, 333, 1427–1430, https://doi.org/10.1126/science.1208336, 2011. a
Ritz, C., Fabre, A., and Letréguilly, A.: Sensitivity of a Greenland ice sheet model to ice flow and ablation parameters: consequences for the evolution through the last climatic cycle, Clim. Dynam., 13, 11–23, https://doi.org/10.1007/s003820050149, 1997. a, b
Robinson, A., AlvarezSolas, J., Montoya, M., Goelzer, H., Greve, R., and Ritz, C.: Yelmo source code snapshot, Zenodo, https://doi.org/10.5281/zenodo.3782650, 2020. a
Rückamp, M., Greve, R., and Humbert, A.: Comparative simulations of the evolution of the Greenland ice sheet under simplified Paris Agreement scenarios with the models SICOPOLIS and ISSM, Polar Sci., 21, 14–25, https://doi.org/10.1016/j.polar.2018.12.003, 2019. a
Rybak, O. and Huybrechts, P.: A comparison of Eulerian and Lagrangian methods for dating in numerical icesheet models, Ann. Glaciol., 37, 150–158, https://doi.org/10.3189/172756403781815393, 2003. a, b, c, d
Schaffer, J., Timmermann, R., Arndt, J. E., Kristensen, S. S., Mayer, C., Morlighem, M., and Steinhage, D.: A global, highresolution data set of ice sheet topography, cavity geometry, and ocean bathymetry, Earth Syst. Sci. Data, 8, 543–557, https://doi.org/10.5194/essd85432016, 2016. a
Schoof, C.: The effect of cavitation on glacier sliding, P. Roy. Soc. A Math. Phy., 461, 609–627, https://doi.org/10.1098/rspa.2004.1350, 2005. a
Schoof, C.: Ice sheet grounding line dynamics: Steady states, stability, and hysteresis, J. Geophys. Res., 112, F03S28, https://doi.org/10.1029/2006JF000664, 2007. a
Schoof, C.: Icesheet acceleration driven by melt supply variability, Nature, 468, 803–806, https://doi.org/10.1038/nature09618, 2010. a
Seroussi, H., Morlighem, M., Larour, E., Rignot, E., and Khazendar, A.: Hydrostatic grounding line parameterization in ice sheet models, The Cryosphere, 8, 2075–2087, https://doi.org/10.5194/tc820752014, 2014. a
Seroussi, H., Nowicki, S., Simon, E., AbeOuchi, A., Albrecht, T., Brondex, J., Cornford, S., Dumas, C., GilletChaulet, F., Goelzer, H., Golledge, N. R., Gregory, J. M., Greve, R., Hoffman, M. J., Humbert, A., Huybrechts, P., Kleiner, T., Larour, E., Leguy, G., Lipscomb, W. H., Lowry, D., Mengel, M., Morlighem, M., Pattyn, F., Payne, A. J., Pollard, D., Price, S. F., Quiquet, A., Reerink, T. J., Reese, R., Rodehacke, C. B., Schlegel, N.J., Shepherd, A., Sun, S., Sutter, J., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., and Zhang, T.: initMIPAntarctica: an ice sheet model initialization experiment of ISMIP6, The Cryosphere, 13, 1441–1471, https://doi.org/10.5194/tc1314412019, 2019. a
Söderlind, G.: Automatic Control and Adaptive TimeStepping, Numer. Algorithms, 31, 281–310, https://doi.org/10.1023/A:1021160023092, 2002. a
Söderlind, G. and Wang, L.: Adaptive timestepping and computational stability, J. Comput. Appl. Math., 185, 225–243, https://doi.org/10.1016/j.cam.2005.03.008, 2006. a, b
Stearns, L. A. and van der Veen, C. J.: Friction at the bed does not control fast glacier flow, Science, 361, 273–277, https://doi.org/10.1126/science.aat2217, 2018. a
van Wessem, J. M., van de Berg, W. J., Noël, B. P. Y., van Meijgaard, E., Amory, C., Birnbaum, G., Jakobs, C. L., Krüger, K., Lenaerts, J. T. M., Lhermitte, S., Ligtenberg, S. R. M., Medley, B., Reijmer, C. H., van Tricht, K., Trusel, L. D., van Ulft, L. H., Wouters, B., Wuite, J., and van den Broeke, M. R.: Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 2: Antarctica (1979–2016), The Cryosphere, 12, 1479–1498, https://doi.org/10.5194/tc1214792018, 2018. a
Weertman, J.: Stability of the junction of an ice sheet and an ice shelf, J. Glaciol., 13, 3–11, https://doi.org/10.3189/S0022143000023327 , 1974. a
Winkelmann, R., Martin, M. A., Haseloff, M., Albrecht, T., Bueler, E., Khroulev, C., and Levermann, A.: The Potsdam Parallel Ice Sheet Model (PISMPIK) – Part 1: Model description, The Cryosphere, 5, 715–726, https://doi.org/10.5194/tc57152011, 2011. a, b, c
Yelmodocs: Yelmo, General model structure – classes and usage, Yelmodocs, avaialable at: https://palmaice.github.io/yelmodocs/, last access: 2 May 2020. a
The name Yelmo refers to a semidomed, rocky mountain in the Guadarrama Mountains outside of Madrid, Spain.