Description and validation of the ice-sheet model Nix v1.0

We present a physical description of the ice-sheet model Nix v1.0, an open-source project intended for collaborative development. Nix is a two-dimensional (flowline combined with a vertical dimension) thermomechanical model written in C and C++ that simultaneously solves for the momentum balance equations, mass conservation and temperature evolution. Nix's velocity solver includes a hierarchy of Stokes approximations: Blatter–Pattyn, depth-integrated higher order and shallow shelf. The grounding-line position is explicitly solved by a moving coordinate system that avoids further interpolations. The model can be easily forced with any external boundary conditions. Nix has been verified for standard test problems, showing versatility from regular machines (lightweight memory allocation) to high-performance computing (multi-threading capabilities). Resolutions below 0.1 km are attainable even with minimal computational resources: Nix's serial run finalizes within hours on a single CPU. Here we show results for a number of benchmark experiments from the Marine Ice Sheet Intercomparison Project (MISMIP) and assess grounding-line migration with an overdeepened bed geometry. Lastly, we further exploit the thermomechanical coupling by designing a suite of experiments where the forcing is a physical variable, unlike previously idealized forcing scenarios where ice temperatures are implicitly fixed via an ice rate factor. That is, we use atmospheric temperature and oceanic temperature anomalies to assess model hysteresis behaviour with active thermodynamics. Our results show that hysteresis in an overdeepened bed geometry is similar for atmospheric and oceanic forcings. Notably, the classical hysteresis loop is widened for both forcing scenarios (i.e. atmospheric and oceanic) if the ice sheet is thermomechanically active as a result of the internal feedback among ice temperature, stress balance and viscosity. These results show that a temperature-dependent ice viscosity provides inertia and stability to the ice sheet, regardless of the particular external forcing applied. In summary, Nix combines rapid computational capabilities with a Blatter–Pattyn stress balance fully coupled to a thermomechanical solver, not only validating against established benchmarks but also offering a powerful tool for advancing our insight into ice dynamics and grounding-line stability.

1 Introduction

Marine ice sheets, such as the present-day West Antarctic Ice Sheet (WAIS), are of particular interest for the glaciological community and have been fundamental objects of study in the last few decades (Payne et al., 2000; Pattyn et al., 2008; Bamber et al., 2009; Feldmann and Levermann, 2015; Shepherd et al., 2018; Martin et al., 2019; Rignot et al., 2019; Robel et al., 2019; Pattyn and Morlighem, 2020; Garbe et al., 2020; DeConto et al., 2021; Joughin et al., 2021; Hill et al., 2023). Since their bedrock lies mostly below sea level, they are prone to rapid changes (Bentley, 1998), leading a number of authors to question their stability (e.g. Bamber et al., 2009; Mouginot et al., 2014; Paolo et al., 2015; Feldmann and Levermann, 2015; Shepherd et al., 2018; Rignot et al., 2019; Robel et al., 2019; Pattyn and Morlighem, 2020; Garbe et al., 2020; Joughin et al., 2021; Hill et al., 2023). A complete WAIS collapse would imply 3–5 m of sea-level rise (Bentley, 1998; Fretwell et al., 2013), leaving the future of the WAIS a key uncertainty for sea-level projections.

An accurate numerical description of the grounding line is thus fundamental for the reliability of such projections. A number of attempts have been made in the past to simulate grounding-line migration within marine ice-sheet models. Weertman (1974) and Thomas and Bentley (1978) proposed that no stable steady states of the grounding line could be found on inland-sloping or retrograde beds. Hindmarsh (1993) later introduced the possibility of “neutral equilibrium” under the premise that the equilibrium position is continuous, and hence there exists an infinite number of equilibrium configurations. More recently, Vieli and Payne (2005) assessed the influence of numerical details and discretization on the dynamics of the grounding line, concluding that a reliable method of treating grounding-line migration within numerical ice-sheet models was unknown. Later studies confirmed the possibility of numerical artefacts (Pattyn et al., 2006; Hindmarsh, 2006; Schoof, 2006a, b, 2011), in agreement with the early works of Weertman (1974) and Thomas and Bentley (1978). Even so, Vieli and Payne (2005) and Pattyn et al. (2006) hypothesized the possibility of neutral equilibrium, first introduced by Hindmarsh (1993). The analytical approach of Schoof (2007a) based on asymptotic expansions eventually concluded that these results were numerical artefacts appearing for certain parameter regimes. Only in the absence of basal sliding has the possibility of non-unique steady states been raised (Nowicki and Wingham, 2008).

Amidst the lack of a reliable model of grounding-line migration, the first Marine Ice Sheet Intercomparison Project (MISMIP; Pattyn et al., 2012) shed light on the agreement of modelling efforts to describe the grounding-line motion and assessed the appropriateness of numerical schemes. The authors proposed a set of benchmark experiments on an idealized two-dimensional bed geometry, concluding that moving-grid models are the most reliable choice from a numerical perspective as the grounding line is part of the solution and no interpolations are required.

MacAyeal and Barcilon (1988) notably showed that a two-dimensional free-floating shelf has no effect on the dynamics of the grounded ice upstream of it (later underlined by Schoof, 2007a). As a result, a boundary condition can be directly imposed at the grounding line that is solely dependent on the ice thickness therein, irrespective of the particular shape or the dynamics of the shelf. A correct description of a two-dimensional marine ice sheet thus relies on an appropriate formulation of the grounded ice dynamics, especially near the terminus position where ice streaming is generally found.

For a comparison with the semi-analytical solutions of Schoof (2007a), ice streaming (i.e. fast-flowing ice due to basal sliding) becomes a necessary condition given that the boundary layer theory assumes rapid sliding near the grounding line. Ice streams are in fact a distinct feature of ice sheets with no counterpart in other geophysical thin-film flows. These regions of rapidly flowing ice exhibit velocities up to 3 orders of magnitude faster than those of the usual glacial ice, yet they only account for a small fraction of the total ice-sheet area (e.g. less than 5 % of the Antarctic Ice Sheet; Bamber et al., 2000). Even so, it is important to represent them correctly to evaluate ice outflow discharge, ice-sheet sensitivity and overall stability.

The rapid flow of ice streams fails to be explained by vertical shearing of ice. In other words, friction at the bed is typically smaller than the driving stress predicted by a lubrication approximation (Whillans and van der Veen, 1997; Joughin et al., 2004). Rather, high ice-stream velocities are caused by the deformation of meltwater-saturated, weak subglacial till (Alley et al., 1986; Blankenship et al., 1986; Engelhardt et al., 1990), thus consistent with geophysical studies showing that basal sliding is fundamentally a sort of Coulomb slip connected with the mechanical failure of plastic till (e.g. Tulaczyk, 1999).

Schoof (2006b) later extended the work to depth-integrated viscous flows used in three-dimensional ice-sheet models. That is, a variational formulation of the two-dimensional shallow-shelf approximation (SSA) equations is given without assumptions on the extension of the sliding domain. In fact, as noted by the author, sliding regions must be determined as part of the solution and are consequently not known a priori. Notably, a solvability condition was also derived (as in Schoof, 2006b) to guarantee the existence of physical solutions. Strictly speaking, if the till is too weak such that the total momentum of applied forces is greater than the maximum momentum of frictional force about a given point, then no solutions are expected to exist.

A variational formulation entails strong consequences from both a physical and a mathematical point of view. Particularly, it eludes explicit manipulation of the unknown sliding domain extension, additionally provides a numerical method for solving the ice flow problem, and ensures the well posedness of the SSA non-linear elliptic equations since they can be derived from a convex and functional bounded below (Schoof, 2006b). However, the time-evolving system of the SSA stress balance coupled to the advection equation is not yet known to be mathematically well posed (Bueler and Brown, 2009).

More recently, Goldberg (2011) derived a higher-order stress approximation using variational methods with similar accuracy to the Blatter–Pattyn momentum equations (Blatter, 1995; Pattyn, 2003), though differences are particularly notable for resolutions below 20 km. The velocity solver was first adapted for multimillennial three-dimensional ice-sheet models such as CISM (Lipscomb et al., 2019), where this depth-integrated velocity approximation was referred to as DIVA. Nevertheless, the DIVA solver had previously been used in continental-scale models by Arthern et al. (2015) and Arthern and Williams (2017). The numerical stability of this solver was systematically studied by Robinson et al. (2022), who showed that the DIVA solver outperformed the remaining solvers in terms of both model performance and the representation of the ice-flow physics itself.

The appropriate stress balance treatment is merely one of the challenges associated with ice streaming and grounding-line stability. Understanding the mechanisms governing its temporal variability also remains a major obstacle, particularly with the aim of developing models of ice-sheet dynamics (Robel et al., 2013). Given the broad range of ice-flow speeds observed in real ice sheets (e.g. Shepherd and Wingham, 2007; Truffer and Fahnestock, 2007; Vaughan and Arthern, 2007), numerical simulations of these rapidly flowing bands are a well-known difficulty, partially due to the fact that fast grounded ice flow is a combination of sliding over a hard and soft bed and shear deformation of the basal. Moreover, high-quality spatially distributed ice observations of near-base conditions are rare, and constraining models in fast-flowing regions becomes challenging (Bueler and Brown, 2009). Various modelling approaches have been considered to correctly describe the large complexity of ice-stream dynamics. Tulaczyk et al. (2000) found that subglacial hydrology yields multiple modes of ice-stream flow in a highly reduced model. Parameterizations of observed small-scale phenomena (e.g. drainage networks) were later considered by coupling a flow-band model and a simple hydrological model (Bougamont et al., 2003; Bougamont and Tulaczyk, 2003). Another flow-band model was employed by van der Wel et al. (2013), additionally introducing a dynamic drainage model.

Two-dimensional models have helped tremendously to understand ice-sheet dynamics from both a theoretical (e.g. Weertman, 1974; Hindmarsh, 1993, 1996; Chugunov and Wilchinsky, 1996; Schoof, 2005, 2006b, 2007a, b, 2011) and a modelling perspective. Numerous authors have contributed to the latter, thus demonstrating the practical use of a two-dimensional setup. Hindmarsh and le Meur (2001) assessed the dynamical processes involved in the retreat of marine ice sheets, with a particular interest in the WAIS during the Last Glacial Maximum. Haseloff and Sergienko (2018) later considered the effect of buttressing on grounding-line dynamics, thus corroborating the findings of existing numerical studies that the stability of confined marine ice sheets is influenced by ice-shelf properties. Other two-dimensional ice-sheet models additionally employ real bedrock geometry sections. This is the case of Pattyn et al. (2006), who studied the role of transition zones in marine ice-sheet dynamics, and Jamieson et al. (2012), where ice-stream stability was investigated on a reverse bed slope. The realism of the two-dimensional setup can also account for glacial isostatic adjustment. To illustrate this, Payne (1995) studied limit cycles in the basal thermal regime of ice sheets considering a constant diffusivity of the asthenosphere. More recently, Bassis et al. (2017) investigated how Heinrich events are triggered by ocean forcing and modulated by isostatic adjustment, though the viscosity dependency on temperature was not considered. Other examples of simplified physics that neglect thermochemical coupling and instead focus on attribution exercises of anthropogenic-induced ice-sheet retreat (e.g. Christian et al., 2022) are consequently biased by unrealistic constant temperatures in both space and time. Lastly, even though ice shelves are not explicitly resolved in two-dimensional models, the potential role of buttressing can also be considered via a parameterization (e.g. Dupont and Alley, 2005; Schoof, 2007a; Jamieson et al., 2012; Robel et al., 2014, 2019).

Despite the extensive research on the topic, important questions regarding the particular effect of thermodynamics remain unanswered. Specifically, it is unclear whether marine ice sheets have discrete steady-surface profiles if ice temperatures can freely evolve in time and what the potential implications would be for the hysteresis behaviour in overdeepened bed geometries. In ice-streaming regions, ice flow occurs mostly along one main direction, thus becoming the preferred axis across which lateral variations are negligible. It is a common approach to reduce the number of horizontal dimensions to the main flow direction so as to minimize computing time while allowing for realistic applications. Nonetheless, the thermal state of the ice and the potential oceanic forcing are still fundamental to understanding the future evolution of ice sheets that have not been considered in low-dimensional models.

In line with the above, we herein introduce the two-dimensional ice-sheet model Nix. Unlike previous two-dimensional models, the default setup consists of a Blatter–Pattyn stress balance fully coupled to a thermodynamical solver that accounts for both vertical and along-flow horizontal heat advection, as well as vertical diffusion. Note that other configurations are also possible given the independent structure of Nix functionalities. In other words, the user can select a particular stress balance description from a hierarchy of Stokes approximations (i.e. Blatter–Pattyn, depth integrated viscosity approximation, SSA or SIA) and optionally solve the associated heat problem. The paper is structured as follows. We begin by describing the technical model design (Sect. 2). We then describe the physical approximations (Sect. 3) and numerics (Sect. 4) of the model. Several benchmarks and idealized experiments are presented in Sect. 5. A thorough discussion is given in Sect. 6, and the work is summarized in Sect. 7.

2 Model design

Nix is open-source software, available under the Creative Commons Attribution 4.0 International License. The model has been derived from scratch with a clear application programming interface (API). It is written in C and C++ for efficiency and extremely fast computing (see Sect. 7) and is readily available to run in any high-performance computing cluster. There are two key dependencies: NetCDF (Rew and Davis, 1990; Brown et al., 1993) and Eigen (Jacob and Guennebaud, 2010) libraries. The former handles tasks for convenient community-standard input/output capability, whereas the latter serves to define vectors, matrices and further necessary computations (Fig. 1).

Figure 1 Overview of Nix's modular structure. Each colour represents a C++ class: dynamics, material, topography, thermodynamics and boundary conditions. The Python wrapper is a user-friendly option, and the code can be compiled without any additional dependencies in any standard high-performance computing cluster.

The current design offers a user-friendly Python wrapper module that handles directory management and compilation, though it can be compiled and run independently. The exact version used to produce the results of this work is archived in a persistent Zenodo repository (Moreno-Parada et al., 2023), while the latest version can be accessed on GitHub at https://github.com/d-morenop/nix (last access: 27 March 2025).

Nix users can optionally select parallel computing (supported by the Eigen library) simply by enabling OpenMP (Graham et al., 2006) on the employed compiler, potentially convenient for high resolutions in the Blatter–Pattyn approximation, where large sparse matrices must be inverted. Nix is extremely well optimized, thus achieving peak performance with only eight logical threads. Both efficient memory allocation and optimized computation combined allow for high-resolution simulations (Δx<0.1 km) that can even be executed on a regular PC. Section 7 thoroughly elaborates on the parallelization details and optimal resource choices. Moreover, it is also possible to use Eigen's matrices, vectors and arrays for fixed size within CUDA kernels (Nickolls et al., 2008).

As a result, Nix presents itself as a remarkably versatile model combining usage simplicity, low computational costs and high-order physics at extremely high resolution (Δx<0.1 km). Moreover, the practical value of the Nix ice-sheet model lies not only in its high-resolution performance, but also in the gap it fills within the model hierarchy spectrum: a thermodynamically coupled and computationally inexpensive two-dimensional model solving for the higher-order Blatter–Pattyn stress balance.

3 Model physics

In this section, the fundamental equations of the model are described. Generally speaking, we consider an ice slab of two spatial dimensions (i.e. horizontal and vertical) by coupling a particular choice of stress balance, the advection equation and the associated heat problem.

Our system evolves thermodynamically in time through three main processes concerning heat propagation: vertical diffusion, horizontal and vertical advection, and internal deformation of the ice. Viscosity is thus dependent on both the strain rate and the temperature. With respect to dynamics, basal friction can be parameterized by three distinct formulations (linear, power law and regularized Coulomb). Additionally, basal friction captures the thermal state of the base through a two-valued friction coefficient encapsulating frozen and thawed bedrock.

3.1 The Blatter–Pattyn approximation

Ice sheets and glaciers are generally described as an incompressible fluid with a low-Reynolds-number flow. Conservation of momentum is ensured through the Stokes equations, a quasi-static stress description where inertial and advective terms are neglected due to the slow movement of the ice.

The typical ice-sheet geometry allows us to further simplify the Stokes flow equations by defining an aspect ratio ε. Given the characteristic length scales for the horizontal and vertical dimensions, ε≪1 (e.g. Greve and Blatter, 2009). Simply by keeping terms of order 𝒪(ε) in the Stokes equations, the Blatter–Pattyn model (Blatter, 1995; Pattyn, 2003) arises with a hydrostatic approximation error of 𝒪(ε²) (Dukowicz et al., 2010; Schoof and Hindmarsh, 2010). This first-order approximation forms an elliptic coercive problem, which is significantly easier to solve than the intricate saddle-point problem of the full-Stokes system.

For the purpose of this work, we consider two spatial dimensions: horizontal x and vertical z. This reduces the computational time considerably and allows for extremely high spatial resolutions (Δx≤0.1 km) while explicitly accounting for the vertical gradients in ice viscosity and velocity. The two-dimensional version of the Blatter–Pattyn model can be written as

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial }{\partial x}}\left(\mathrm{4}\mathit{\eta }{\displaystyle \frac{\partial u}{\partial x}}\right)+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\eta }{\displaystyle \frac{\partial u}{\partial z}}\right)=\mathit{\rho }g{\displaystyle \frac{\partial h}{\partial x}},\end{array}$$

where ρ is the ice density, g is the gravitational acceleration, η(x,z) is the effective viscosity, h(x) is the surface elevation and u(x,z) is the ice velocity.

The problem is subjected to a set of boundary conditions. Nix considers potential friction at the base of the ice and a free surface on the upper boundary (Veen and Whillans, 1989). In terms of velocity gradients, the free surface condition can be expressed as (Pattyn, 2003)

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial u}{\partial z}}=\mathrm{4}{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{\partial h}{\partial x}},\end{array}$$

and the basal drag is defined as the sum of all resistive forces:

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial u}{\partial z}}=\mathrm{4}{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{\partial b}{\partial x}}+{\displaystyle \frac{\mathit{\tau }\left(u\right)}{\mathrm{2}\mathit{\eta }}},\end{array}$$

for the base (z=b) in the presence of potential drag τ(u) (see Sect. 3.4 for a thorough description on basal friction). A stress-free base can be obtained simply by setting τ=0 in Eq. (3).

We further assume an ice divide at one end of the domain (x=0), where $\partial u/\partial x=\mathrm{0}$, and hydrostatic equilibrium at the shelf–ocean boundary (x=L), where the water pressure balances the longitudinal stress gradient. The full problem thus takes the following succinct form (hereinafter, subscripts denote partial differentiation):

$$\begin{array}{}\text{(4)}& \left\{\begin{array}{ll}{\left(\mathrm{4}\mathit{\eta }{u}_x\right)}_x+{\left(\mathit{\eta }{u}_z\right)}_z=\mathit{\rho }g{h}_x,& x\in \mathcal{I},z\in \mathcal{L},\\ u_z=\mathrm{4}{u}_x{h}_x,& x\in \mathcal{I},z\in \partial {\mathcal{L}}^{+},\\ \mathrm{2}\mathit{\eta }{u}_z=\mathrm{8}\mathit{\eta }{u}_x{b}_x+\mathit{\tau },& x\in \mathcal{I},z\in \partial {\mathcal{L}}^{-},\\ u_x=\mathrm{0},& x=\mathrm{0},z\in \mathcal{L},\\ \mathrm{8}\mathit{\eta }{u}_x=\mathit{\rho }g{H}^{\mathrm{2}}-{\mathit{\rho }}_{\mathrm{w}}g{z}^{\mathrm{2}},& x=L,z\in \mathcal{L},\end{array}\right.\end{array}$$

where ρ_w is the water density; H is the ice thickness evaluated at the grounding line x=L; and the $\partial {\mathcal{L}}^{\pm }$ symbols denote the upper and lower vertical boundaries, respectively.

3.2 The depth-integrated viscosity approximation

We now briefly describe the mathematical problem underlying the depth-integrated viscosity approximation (DIVA) stress balance in Cartesian coordinates (Goldberg, 2011; Lipscomb et al., 2019).

As for the Blatter–Pattyn model, we consider only one horizontal dimension, leaving the nature of DIVA/SSA equations unaltered. This allows us to regard our model as a longitudinal section of a three-dimensional ice stream:

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial }{\partial x}}\left(\mathrm{4}\stackrel{\mathrm{‾}}{\mathit{\eta }}H{\displaystyle \frac{\partial u}{\partial x}}\right)+\mathit{\tau }\left(u\right)=\mathit{\rho }gH{\displaystyle \frac{\partial h}{\partial x}}.\end{array}$$

Since the stress balance is also a second-order partial differential equation on the velocity, we again need two boundary conditions. Analogously to the Blatter–Pattyn approximation, we assume an ice divide at one end. At the other end of the domain, the problem is subjected to a dynamic boundary condition that accounts for the balance between cryostatic and hydrostatic pressures. Thus, we can express the DIVA/SSA boundary problem in the following compact form:

$$\begin{array}{}\text{(6)}& \left\{\begin{array}{ll}{\left(\mathrm{4}\stackrel{\mathrm{‾}}{\mathit{\eta }}H{u}_x\right)}_x+\mathit{\tau }\left(u\right)=\mathit{\rho }gH{h}_x,& x\in \mathcal{I},\\ u=\mathrm{0},& x=\mathrm{0},\\ \mathrm{8}\stackrel{\mathrm{‾}}{\mathit{\eta }}{u}_x=\mathit{\rho }g{H}^{\mathrm{2}}-{\mathit{\rho }}_{\mathrm{w}}g{D}^{\mathrm{2}},& x=L,\end{array}\right.\end{array}$$

where D is the distance from the sea surface to the bottom of the ice.

Equation (5) is an elliptic non-linear differential equation. In the purely SSA form (neither velocity nor viscosity dependency on z, i.e. $u=\stackrel{\mathrm{‾}}{u}$ and $\mathit{\eta }=\stackrel{\mathrm{‾}}{\mathit{\eta }}$), it constitutes the simplest form of longitudinal stress balance derivable from the Stokes model (Bueler and Brown, 2009). We can then solve for the velocity u(x) by integrating Eq. (6) if the functions H(x), $\stackrel{\mathrm{‾}}{\mathit{\eta }}\left(x\right)$, h(x) and τ(u) are known. We have implemented an implicit algorithm so as to numerically integrate the one-dimensional DIVA/SSA equation (see integrating scheme description in Appendix A2).

3.3 The advection coupling

Given that all models presented herein provide a quasi-static description of the ice flow, the stress balance does not determine the temporal evolution of the system but rather represents an equilibrium state for a particular ice thickness H(x) and viscosity η(x,z) configuration. The temporal evolution is generally considered by coupling the stress balance to the advection equation:

$$\begin{array}{}\text{(7)}& H_t+{\left(\stackrel{\mathrm{‾}}{u}H\right)}_x=S\left(x\right),\end{array}$$

where S(x) is the surface mass balance. Given that Eq. (7) is a first-order equation, we only need one boundary condition $H(x=\mathrm{0},t)$ and the consequent initial condition $H(x,t=\mathrm{0})$.

We now couple Eqs. (6) and (7) to study the evolution of the ice thickness H(x,t) governed by the advection equation, where the velocity field u(x,z) satisfies the stress balance imposed by a particular choice of the Stokes approximation. That is, our problem takes the following mathematical form:

$$\begin{array}{}\text{(8)}& \left\{\begin{array}{ll}{H}_t+{\left(\stackrel{\mathrm{‾}}{u}H\right)}_x=S\left(x\right),& x\in \mathcal{I},t>\mathrm{0},\\ {\left(\mathrm{4}\mathit{\eta }{u}_x\right)}_x+{\left(\mathit{\eta }{u}_z\right)}_z=\mathit{\rho }g{h}_x,& x\in \mathcal{I},z\in \mathcal{L},\\ H=H_{\mathrm{0}},& x\in \mathcal{I},t=\mathrm{0},\\ h_x=\mathrm{0},& x=\mathrm{0},t>\mathrm{0},\\ u=\mathrm{0},& x=\mathrm{0},t>\mathrm{0},\\ \mathrm{8}\mathit{\eta }{u}_x=\mathit{\rho }g{H}^{\mathrm{2}}-{\mathit{\rho }}_{\mathrm{w}}g{z}^{\mathrm{2}},& x=L,z\in \mathcal{L},t>\mathrm{0}.\end{array}\right.\end{array}$$

From a purely physical perspective, Eq. (8) describes a fluid membrane of variable thickness driven by its own weight that evolves in time due to advection.

3.4 Basal friction

Basal shear stress can generally be expressed as a function of the sliding velocity u_b and the effective pressure N; i.e. ${\mathit{\tau }}_{\mathrm{b}}=f(u_{\mathrm{b}},N$). The physical properties of the material over which the ice may potentially slide can correspond either to a hard bedrock flow (e.g. Weertman, 1957) or to a Coulomb plastic rheology (e.g. Tulaczyk et al., 1998). Moreover, the influence of the sliding velocity on τ_b is often represented by a power friction law, although a regularization term u₀ – accounting for local properties of the bed – has been shown to outperform such a power law in both pressurized ice experiments (Zoet and Iverson, 2020) and observations (Minchew et al., 2018; Stearns and van der Veen, 2018; Joughin et al., 2019)

As a result, Nix can calculate the basal shear stress (i.e. basal drag) via two independent formulations: a pseudo-plastic power law (Schoof, 2010; Aschwanden et al., 2013) and the regularized-Coulomb law (Schoof, 2005; Gagliardini et al., 2007; Joughin et al., 2019). The former reads

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

where $u_{\mathrm{0}}=\mathrm{100}\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$ and c_b is a spatially variable friction coefficient defined below. We focus on two particular cases of the pseudo-plastic law based upon the choice of the exponent q, namely the linear law (q=1; e.g. Quiquet et al., 2018) and the purely plastic law (q=0).

On the other hand, the regularized-Coulomb formula is given by

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

behaving as a power law for small sliding velocities (i.e. u_b<u₀) while always yielding a bounded friction value for arbitrarily high velocities (i.e. u_b≫u₀). Following Zoet and Iverson (2020), we set $q=\mathrm{1}/\mathrm{5}$ and $u_{\mathrm{0}}=\mathrm{100}\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}$ by default to ensure a reasonable transition to the steady-state shear stress supported by the till bed. The same study empirically established that q remains unaffected by variations in the detailed bed surface geometry.

The basal drag coefficient β is usually defined as

$$\begin{array}{}\text{(11)}& \mathit{\beta }=c_{\mathrm{b}}\left(x\right)N,\end{array}$$

where N=ρgH is the overburden pressure exerted by the ice column and c_b(x) is a coefficient that reflects the bedrock characteristics.

Nevertheless, for simplicity and consistency with prior benchmark experiments such as MISMIP (Pattyn et al., 2012), the model also allows us to represent basal friction as

$$\begin{array}{}\text{(12)}& {\mathit{\tau }}_{\mathrm{b}}=C{u}^q.\end{array}$$

With the chosen value of $C=\mathrm{7.624}\times {\mathrm{10}}^{\mathrm{6}}\mathrm{Pa}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}{\mathrm{s}}^{\mathrm{1}/\mathrm{3}}$ and $q=\mathrm{1}/\mathrm{3}$, a sliding velocity of about 35 m yr⁻¹ yields a basal shear stress of 80 kPa.

3.5 Thermodynamics

The ice temperature in the flow line depends on the two spatial dimensions x and z (horizontal and vertical, respectively) along with time (i.e. $\mathit{\theta }=\mathit{\theta }(x,z,t)$). Heat transfer is further considered to occur due to vertical diffusion, both horizontal and vertical advection, and internal heat deformation. Energy conservation is ensured in a classical approach by a balance equation that neglects horizontal diffusion (Greve and Blatter, 2009):

$$\begin{array}{}\text{(13)}& \left\{\begin{array}{ll}\mathit{\rho }c{\mathit{\theta }}_t=k{\mathit{\theta }}_{zz}-\mathit{\rho }c\left(u{\mathit{\theta }}_x+w{\mathit{\theta }}_z\right)+\mathrm{\Phi },& \\ & x\in \mathcal{I},z\in \mathcal{L},t>\mathrm{0},\\ \mathit{\theta }={\mathit{\theta }}_{\mathrm{0}},& x\in \mathcal{I},z\in \mathcal{L},t=\mathrm{0},\\ {\mathit{\theta }}_z=-G/k,& x\in \mathcal{I},z=\partial {\mathcal{L}}^{-},t>\mathrm{0},\\ \mathit{\theta }={\mathit{\theta }}_L,& x\in \mathcal{I},z=\partial {\mathcal{L}}^{+},t>\mathrm{0},\end{array}\right.\end{array}$$

where k is the ice conductivity, c is the specific heat capacity, $\mathrm{\Phi }=\mathrm{4}\mathit{\eta }{\dot{\mathit{\epsilon }}}^{\mathrm{2}}$ denotes the internal strain heating, G is the geothermal heat flow, θ₀ is the initial temperature profile and θ_L is the surface ice temperature. The $\partial {\mathcal{L}}^{\pm }$ symbols denote the upper and lower vertical boundaries, respectively.

The energy balance is discretized using an upwind scheme with a forward Euler step and centred differences for the spatial derivatives (see Appendix A4 for a detailed description).

3.6 Viscosity

We consider Glen's flow law (Glen, 1995; Nye, 1957) to relate the shear stress, the ice temperature and the pressure of isotropic polycrystalline ice. Formation of anisotropic fabric is considered via a flow enhancement factor.

As shown in Sect. 3.1, the Blatter–Pattyn stress balance equations define the effective viscosity as

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

where B is the ice hardness, n=3 is the exponent in Glen's flow law, ${\dot{\mathit{\epsilon }}}^{\mathrm{2}}$ is the effective strain rate and ${\dot{\mathit{\epsilon }}}_{\mathrm{0}}^{\mathrm{2}}$ is a regularization factor to elude potential singularities when velocity gradients are zero. Notably, for a two-dimensional model with explicit thermodynamics, the viscosity expression further simplifies the expression of B and ${\dot{\mathit{\epsilon }}}^{\mathrm{2}}$:

$$\begin{array}{}\text{(15)}& {\displaystyle }B=& {\displaystyle }A(\mathit{\theta }{)}^{-\mathrm{1}/n},\text{(16)}& {\displaystyle }{\dot{\mathit{\epsilon }}}^{\mathrm{2}}=& {\displaystyle }{\left({\displaystyle \frac{\partial u}{\partial x}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial u}{\partial z}}\right)}^{\mathrm{2}},\end{array}$$

where n=3 is the Glen flow exponent. A(θ) is the rate factor and follows Arrhenius' law:

$$\begin{array}{}\text{(17)}& A\left(\mathit{\theta }\right)=A_{\mathrm{0}}E{e}^{-Q/R\mathit{\theta }},\end{array}$$

where A₀ and Q are the temperature-dependent rate factor coefficient and activation energy, respectively (Greve and Blatter, 2009). E_f is the so-called enhancement factor, commonly used to approximate the effect of anisotropic flow. It is possible to specify different values of the enhancement factor for different flow regimes (shear or stream). Typical values of the enhancement factor for the shearing and streaming regimes are E_shr=3.0 and E_strm=0.7, respectively (Ma et al., 2010). Here we use a default value of E=1.0 for both.

For the vertically integrated stress balance models (i.e. DIVA and SSA), Eqs. (14) and (15) are modified slightly by computing the vertically averaged quantities $\stackrel{\mathrm{‾}}{\mathit{\eta }}$ and $\stackrel{\mathrm{‾}}{B}$ following the generic formula $\stackrel{\mathrm{‾}}{f}=\frac{\mathrm{1}}{H}{\int }_b^{h}fdz$.

3.7 Grounding line

Nix aims to simulate the flow of a sliding ice sheet. Since the longitudinal stress at the grounding line x=L is simply a function of the ice thickness therein, i.e. H(x=L) for a two-dimensional ice sheet (Schoof, 2007a), the behaviour of grounded ice and the location of the grounding line itself are completely independent of the floating part.

Neither the potential distinct shapes of the ice shelf (e.g. due to sub-shelf melting) nor the calving affect the dynamics of grounded ice. Thus, the flotation condition and the stress condition (Eq. 8) can be considered boundary conditions at the grounding line. These two conditions are in fact sufficient to study the ice thickness evolution and the grounding-line migration.

Following Hindmarsh (1996), an explicit expression for the grounding-line migration rate $\dot{L}$ can be readily obtained from a total differentiation of the flotation condition:

$$\begin{array}{}\text{(18)}& \dot{L}\equiv {\displaystyle \frac{\mathrm{d}L}{\mathrm{d}t}}={\displaystyle \frac{\mathit{\varrho }{D}_t+{\left(\stackrel{\mathrm{‾}}{u}H\right)}_x-S}{H_x-\mathit{\varrho }{D}_x}},\end{array}$$

where D is the water depth at the grounding line and $\mathit{\varrho }={\mathit{\rho }}_{\mathrm{w}}/\mathit{\rho }$ is the water-to-ice density ratio, respectively.

More recent studies suggest that the maximum terminus thickness is bounded by the yield strength of ice τ_c (Bassis and Walker, 2012; Bassis and Jacobs, 2013). Hence, a maximum ice thickness at the terminus occurs when the stress exceeds the depth-integrated strength of ice:

$$\begin{array}{}\text{(19)}& H^{\text{max}}={\displaystyle \frac{{\mathit{\tau }}_c}{\mathit{\rho }g}}+\sqrt{{\left({\displaystyle \frac{{\mathit{\tau }}_c}{\mathit{\rho }g}}\right)}^{\mathrm{2}}+\mathit{\varrho }{D}^{\mathrm{2}}},\end{array}$$

thus constraining the terminus thickness such that $H(x=L,t)\le H^{\text{max}}$.

This approach eludes semi-empirical parameterizations of the calving (as in Meier and Post, 1987; Van der Veen, 1996, 2002; Nick et al., 2009, 2010) and further provides a lower bound on the rate of grounding-line advance (Bassis et al., 2017). Combining the continuity equation and the material derivative of H^max (Eq. 19), an expression for the rate of advance/retreat of the terminus can be readily obtained:

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\mathrm{d}L}{\mathrm{d}t}}\ge {\displaystyle \frac{H_t}{H_x^{\text{max}}-H_x}},\end{array}$$

at x=L. The negative sign indicates retreat.

The inequality in Eq. (20) is analogous to the grounding-line migration derived for a marine ice sheet by Schoof (2007a, b). Particularly, if H^max is given by the flotation condition, Eq. (20) exactly reproduces the grounding-line position derived by Schoof (2007a) (Bassis et al., 2017).

3.8 Sub-shelf melting parameterization

Oceanic melting beneath ice shelves is the main driver of the current mass loss of the Antarctic Ice Sheet. For this reason, Nix considers various melting parameterizations, such as simple scaling with far-field thermal driving (e.g. Favier et al., 2019).

We adhere to local yet physically based parameterizations based on ocean circulation models (Grosfeld et al., 1997). That is, the linear dependency can be expressed as

$$\begin{array}{}\text{(21)}& M={\mathit{\gamma }}_T\mathit{\varrho }{\displaystyle \frac{c_{\text{po}}}{L_{\mathrm{i}}}}\left(T-T_{\mathrm{0}}\right),\end{array}$$

where γ_T is the heat exchange velocity, T₀ is a reference temperature, c_po is the specific heat capacity of the ocean mixed layer and L_i is the latent heat of fusion of ice.

This linear formulation, with a constant exchange velocity γ_T, assumes a circulation in the ice-shelf cavity that is independent of the ocean temperature. This assumption is not supported by modelling (Holland et al., 2008; Donat-Magnin et al., 2017) or observational studies (Jenkins et al., 2018), which suggest a larger circulation in response to a warmer ocean, subsequently increasing melt rates. One manner to account for this positive feedback is by considering a quadratic dependency (Holland et al., 2008):

$$\begin{array}{}\text{(22)}& M={\mathit{\gamma }}_T{\left(\mathit{\varrho }{\displaystyle \frac{c_{\text{po}}}{L_{\mathrm{i}}}}\right)}^{\mathrm{2}}{\left(T-T_{\mathrm{0}}\right)}^{\mathrm{2}}.\end{array}$$

These two parameterizations have been employed in numerous studies (e.g. review in Asay-Davis et al., 2017; Favier et al., 2019). This melt rate is included as an additional frontal ablation term in the ice flux computation (Eq. A16), and it amounts to an additional outflow of ice beyond the grounding-line velocity provided by the stress balance. By default, Nix uses a quadratic parameterization.

4 Model numerics

4.1 Moving-grid transformation

Nix uses a nonuniform moving spatial grid that explicitly solves the grounding-line position (Fig. 2). By default, the grid point distribution yields higher resolution near the grounding line following a polynomial or an exponential law (details in Appendix A). Evenly spaced grids are also possible by setting the polynomial order to 1.

As already noted by Pattyn et al. (2012), moving-grid models are presumably the best choice in two-dimensional models from a numerical perspective, as the grounding-line position L(t) is part of the solution and no interpolations are required. Given that neither the terminus position L(t) (i.e. the grounding line) nor the ice thickness H(x,t) is fixed in time, we adopt a moving grid to trace their positions:

$$\begin{array}{}\text{(23)}& \mathit{\sigma }={\displaystyle \frac{x}{L\left(t\right)}},\mathit{\zeta }={\displaystyle \frac{z-b\left(x\right)}{H(x,t)}},\mathit{\tau }=t,\end{array}$$

thus mapping the time-dependent intervals $\mathrm{0}\le x\le L\left(t\right)$ and $\mathrm{0}\le z\le H(x,t)$ into fixed ones $\mathrm{0}\le \mathit{\sigma }\le \mathrm{1}$ and $\mathrm{0}\le \mathit{\zeta }\le \mathrm{1}$. The variable τ is merely introduced to distinguish partial derivatives defined by keeping both σ and ζ constant (as opposed to keeping x and z constant).

Figure 2 Nix's staggered grid definition follows an Arakawa C scheme (Arakawa and Lamb, 1977). The number of grid points in the horizontal r and vertical p dimensions are fixed in time (see Appendix A). As we employ a moving grid, the position of the last horizontal point $(r-\mathrm{1}/\mathrm{2})$ explicitly tracks the grounding line L(t). The grid spacing in the vertical Δζ_j and horizontal Δσ_i axes can be spatially dependent as Nix allows for nonuniform grids. Ghost points required to satisfy the boundary condition at the ice divide are noted in grey on the edge.

As a result, the corresponding derivatives contain additional terms (upon application of the Leibniz rule):

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\partial }{\partial z}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{H}}{\displaystyle \frac{\partial }{\partial \mathit{\zeta }}},\text{(25)}& {\displaystyle \frac{\partial }{\partial x}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{L}}{\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}-{\displaystyle \frac{\mathrm{1}}{H}}\left[\left(\mathrm{1}-\mathit{\zeta }\right){\displaystyle \frac{\partial b}{\partial \mathit{\sigma }}}+\mathit{\zeta }{\displaystyle \frac{\partial H}{\partial \mathit{\sigma }}}\right]{\displaystyle \frac{\partial }{\partial \mathit{\zeta }}},\text{(26)}& {\displaystyle \frac{\partial }{\partial t}}& {\displaystyle }={\displaystyle \frac{\partial }{\partial \mathit{\tau }}}-{\displaystyle \frac{\mathit{\sigma }}{L}}{\displaystyle \frac{\partial L}{\partial \mathit{\tau }}}{\displaystyle \frac{\partial }{\partial \mathit{\sigma }}}-{\displaystyle \frac{\mathit{\zeta }}{H}}{\displaystyle \frac{\partial H}{\partial \mathit{\tau }}}{\displaystyle \frac{\partial }{\partial \mathit{\zeta }}}.\end{array}$$

For simplicity and analogously for the Blatter–Pattyn approximation, the advection equation coupled with the SSA/DIVA stress balance can be written in terms of the new variables. Thus, Eq. (8) reads

$$\begin{array}{}\text{(27)}& \left\{\begin{array}{ll}L{H}_{\mathit{\tau }}-\mathit{\sigma }\dot{L}{H}_{\mathit{\sigma }}+{\left(uH\right)}_{\mathit{\sigma }}=LS(\mathit{\sigma },\mathit{\tau }),& \mathit{\sigma }\in \stackrel{\mathrm{̃}}{\mathcal{I}},\mathit{\tau }>\mathrm{0},\\ {\left(\mathrm{4}\stackrel{\mathrm{‾}}{\mathit{\eta }}H{u}_{\mathit{\sigma }}\right)}_{\mathit{\sigma }}+\mathit{\tau }{L}^{\mathrm{2}}=L^{\mathrm{2}}\mathit{\rho }g{h}_{\mathit{\sigma }}& \mathit{\sigma }\in \stackrel{\mathrm{̃}}{\mathcal{I}},\\ H=H_{\mathrm{0}},& \mathit{\sigma }\in \stackrel{\mathrm{̃}}{\mathcal{I}},\mathit{\tau }=\mathrm{0},\\ H_{\mathit{\sigma }}=\mathrm{0},& \mathit{\sigma }=\mathrm{0},\mathit{\tau }>\mathrm{0},\\ u=\mathrm{0},& \mathit{\sigma }=\mathrm{0},\\ \mathrm{4}\mathit{\eta }{u}_{\mathit{\sigma }}=\left(\mathit{\rho }g{H}^{\mathrm{2}}-{\mathit{\rho }}_{\mathrm{w}}g{D}^{\mathrm{2}}\right)L/\mathrm{2},& \mathit{\sigma }=\mathrm{1},\end{array}\right.\end{array}$$

where $\stackrel{\mathrm{̃}}{\mathcal{I}}\in [\mathrm{0},\mathrm{1}]$ is the transformed interval and the subscripts denote partial differentiation.

Likewise, the third evolution equation that determines the behaviour of our system (i.e. the energy balance, Eq. 13) can be readily obtained in terms of our new variables:

$$\begin{array}{}\text{(28)}& \left\{\begin{array}{ll}\mathit{\rho }c\left[L{\mathit{\theta }}_{\mathit{\tau }}-\mathit{\sigma }\dot{L}{\mathit{\theta }}_{\mathit{\sigma }}-\mathit{\zeta }L{H}_{\mathit{\tau }}{\mathit{\theta }}_{\mathit{\zeta }}/H\right]=& \\ kL{\mathit{\theta }}_{\mathit{\zeta }\mathit{\zeta }}/H^{\mathrm{2}}-\mathit{\rho }cu\left[{\mathit{\theta }}_{\mathit{\sigma }}-(b_{\mathit{\sigma }}+\mathit{\zeta }{H}_{\mathit{\sigma }}){\mathit{\theta }}_{\mathit{\zeta }}/H\right]+L\mathrm{\Phi },& \\ & \mathit{\sigma }\in \stackrel{\mathrm{̃}}{\mathcal{I}},\mathit{\zeta }\in \stackrel{\mathrm{̃}}{\mathcal{L}},\mathit{\tau }>\mathrm{0},\\ \mathit{\theta }={\mathit{\theta }}_{\mathrm{0}},& \mathit{\sigma }\in \stackrel{\mathrm{̃}}{\mathcal{I}},\mathit{\zeta }\in \stackrel{\mathrm{̃}}{\mathcal{L}},\mathit{\tau }=\mathrm{0},\\ {\mathit{\theta }}_z=-G/k,& \mathit{\sigma }\in \stackrel{\mathrm{̃}}{\mathcal{I}},\mathit{\zeta }=\partial {\stackrel{\mathrm{̃}}{\mathcal{L}}}^{-},\mathit{\tau }>\mathrm{0},\\ \mathit{\theta }={\mathit{\theta }}_L,& \mathit{\sigma }\in \stackrel{\mathrm{̃}}{\mathcal{I}},\mathit{\zeta }=\partial {\stackrel{\mathrm{̃}}{\mathcal{L}}}^{+},\mathit{\tau }>\mathrm{0},\end{array}\right.\end{array}$$

where the transformed intervals are again denoted by $\stackrel{\mathrm{̃}}{\mathcal{I}}$ and $\stackrel{\mathrm{̃}}{\mathcal{L}}$, respectively.

4.2 Spatial integration

4.2.1 Implicit scheme and Picard iteration

The lateral boundary condition is in fact non-trivial to implement using an explicit scheme (e.g. a shooting-like method) since it depends on the first spatial derivative of the velocity at the terminus position σ=1, which might lead to convergence issues. Nix thus includes an alternative velocity solver based on an implicit discretization scheme of all stress balance models described in Sect. 3 (numerical details in Appendix A).

To account for the potential non-linearity in the velocity as a consequence of the viscosity and basal friction τ(u), the implicit solver uses a initial guess τ₀ and η₀ and then enters a Picard iteration (see Theorem 2.2 in Teschl, 2012). A solution is hence obtained when the convergence criterion

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\left|\right|u^n-u^{n-\mathrm{1}}\left|\right|}{\left|\right|u^n\left|\right|}}<{\mathit{\varphi }}_{\text{tol}}\end{array}$$

is satisfied. The tolerance ϕ_tol can be set by the user, but the default value is 10⁻⁶.

For the Blatter–Pattyn approximation, a sparse matrix must be solved in each Picard iteration. To do so, we apply a biconjugate gradient stabilized method (commonly known as BiCGSTAB, van der Vorst, 1992) with an incomplete preconditioner (ILUT). In contrast, the DIVA/SSA approximation solely requires solving a tridiagonal matrix at each Picard iteration step, where the ice viscosity is updated. A tridiagonal solver algorithm is implemented as a subroutine within Nix to avoid additional external dependencies (see Appendix A).

4.3 Time integration

Once the velocity field u(x,z) is obtained for a given set of boundary conditions and a particular ice thickness initial distribution H(σ,τ₀), we can compute the time evolution of the latter as a consequence of the advection imposed by $\stackrel{\mathrm{‾}}{u}\left(x\right)$ and the surface mass balance S(x,t) (Eq. 7). Thus, this coupled system, formed by the momentum conservation and the continuity equation (Eq. 27), is fully integrated in two steps: first, a spatial integration to obtain the velocity (where the ice viscosity is known) and, second, a forward time integration to determine the new ice thickness. Lastly, the energy balance equation is integrated to compute the new temperature field.

Specifically, for a given initial ice thickness distribution H(x,t₀), the stress balance equation is spatially integrated, thus yielding the velocity u(x,z). Then, the solution u(x,z) (at t₀) and H(x,t₀) allow us to integrate the continuity equation forward in time, consequently obtaining $H(x,t_{\mathrm{0}}+\mathrm{\Delta }t)$. Additionally, this new ice thickness distribution yields $\mathit{\theta }(x,z,t_{\mathrm{0}}+\mathrm{\Delta }t)$, thus constituting a self-consistent iterative method.

5 Methods and experimental setup

Prior to any comprehensive description of the results, we must test whether Nix is capable of reproducing the benchmark tests of the Marine Ice Sheet Model Intercomparison Project (MISMIP; Pattyn et al., 2012). To this end, we perform all three MISMIP experiments: relaxation to steady state on a downward-sloping bed (Exp. 1), reversal of parameters (Exp. 2) and hysteresis on an overdeepening bed (Exp. 3). The aim of Exp. 1 was to show that there should be a single stable equilibrium profile on a downward-sloping bed. A backward parameter relaxation in Exp. 2 was intended to demonstrate that grounding-line positions should be identical during advance and retreat, as steady states are unique. Exp. 3 was designed to assess whether ice-sheet models exhibit hysteresis behaviour and has become a benchmark for testing the capability of numerical models to simulate grounding-line migration.

First, we adopt the exact same problem definition so as to perform a one-to-one comparison. Next, we run an ensemble of simulations to address the question of whether the hysteresis with respect to model parameter variations found in MISMIP Exp. 3 is still present even if the thermal state of the ice can evolve in time (as opposed to the idealized constant ice factor A set in Pattyn et al., 2012). Fixing A uniquely determines a constant ice temperature, since A(T) is a bijective function of temperature. We therefore impose an atmospheric forcing (i.e. the ice surface boundary condition) that spans a wide range of realistic temperatures. As the geothermal heat flux provides a positive energy contribution, we expect a different thermal equilibrium profile for each imposed surface temperature. This yields a different viscosity field for each scenario, consequently leading to a different equilibrium velocity. As noted by Sergienko et al. (2013), the temperature profile is mostly determined by horizontal advection in streaming regions, thus bringing forward a strongly non-linear feedback worthy of attention.

Lastly, we force the system via ocean temperature anomalies with respect to a reference value T₀ so that $\mathrm{\Delta }{T}_{\text{oce}}=T_{\text{oce}}-T_{\mathrm{0}}$ while keeping the air temperature constant throughout the simulation. We then convert these temperature anomalies into sub-shelf melting at the grounding line (e.g. Favier et al., 2019) by computing any of the parameterizations described in Sect. 3.8. Even though the air temperature is held constant (i.e. the boundary condition of our heat problem), the thermal state of the ice may evolve as both the thickness and the extent are perturbed by the changing sub-shelf melting at the grounding line. Our particular ocean forcing consists of steps of 0.5 °C evenly spaced in time by 30 kyr to ensure equilibration, from ΔT_oce=0 to a maximum applied anomaly of 7 °C. Then, we reverse the forcing to recover the unperturbed state (i.e. zero anomaly).

It is worth noting that the basal friction remains identical to that in the MISMIP experiments for both the atmospheric and the oceanic forcings. This means that no additional dependency of friction on temperature or hydrology is considered.

6 Results

All simulations shown herein ran at a horizontal resolution of Δx=2 km and 35 vertical layers. The particular mesh employed is unevenly spaced, with an increasing density of points toward the base and near the grounding line, following an exponential relation (see Appendix A). The different experimental setups are summarized in Table 1.

Table 1 Nix suite of experiments. The first row replicates the MISMIP benchmark tests, whereas MISMIP+therm explores the hysteresis behaviour of a thermomechanically active ice sheet in two different forcing scenarios: atmospheric and oceanic.

6.1 MISMIP benchmark experiments

As a performance test for the Nix ice-sheet model, simulations (Fig. 3) fairly reproduce the results shown by models that employ a stretched coordinate system like ours.

Figure 3 (a, c) Ice-sheet extent. (b, d) Grounding-line position as a function of the MISMIP forcing A for three independent Stokes approximations: SSA, DIVA and Blatter–Pattyn. The grey line represents the analytical solution at equilibrium from Schoof (2007a), with the solid line indicating a stable branch and the dashed line indicating an unstable branch. The markers represent Nix's results after the equilibration time given in Pattyn et al. (2012). Bed geometries correspond to Experiments 1 and 2 (both advance and retreat) in the first row and Experiment 3 in the second row.

Figure 3b shows both the advancing and the retreating phase in Pattyn et al. (2012) (Experiments 1 and 2, respectively), with equilibrium grounding-line positions that coincide and points that thus overlap. Additionally, we find no stable equilibrium states for the inland-sloping bed of Exp. 3 (Fig. 3c), in agreement with the theoretical considerations by Weertman (1974) and Schoof (2007a). That is, as we gradually decrease A (i.e. increasing ice viscosity), the grounding-line position advances across the downward-sloping bed until the upward-sloping region is reached (Fig. 3c). The ice flux at the grounding line then continues to increase as A decreases. As illustrated in Fig. 3d, when the ice flux is large enough that a stable solution exists beyond the unstable region (on the right-hand side of the bedrock peak in Fig. 3c), the grounding line traverses the upward-sloping sector, reaching a new stable solution.

Additionally, we repeat the three MISMIP experiments using the more sophisticated velocity solvers available in Nix: DIVA and Blatter–Pattyn. A direct inspection of Fig. 3b and d reveals that the solutions are nearly identical to those of the simpler SSA version, for both the downward-sloping and the overdeepening beds. The hysteresis is particularly well captured in all three Stokes approximations. Thus, we use the SSA solver in the remainder of the current work to minimize computational costs, unless otherwise stated.

6.2 MISMIP and thermodynamics

To exploit the fact that Nix is fully coupled with a thermodynamic solver, we further investigate the equilibrium states (Schoof, 2007a) when the system is forced via two different forcings: air temperatures T_air and ocean temperature anomalies ΔT_oce. Both describe more realistic conditions with slight variations. The former implies the same underlying perturbation mechanism (as for the idealized rate factor A forcing): temperature changes within the ice modify its viscosity so that the grounding line migrates to reach a new equilibrium position. Nevertheless, when forcing the system with ocean temperature anomalies ΔT_oce while keeping the air temperature constant, we perturb the system via an additional outflow term at the grounding line. By separately studying each mechanism, we can determine whether a marine-terminating ice sheet might undergo hysteresis under different forcings.

6.2.1 Air temperature T_air forcing

With the aim of building a thermomechanically active version of the MISMIP experiments, the natural choice is to convert the idealized ice rate factor forcing in MISMIP into temperatures (via Arrhenius' law, Eq. 17) and then use them explicitly as a forcing of the new experimental setup. For a more sophisticated forcing, a vertical dependency of the temperature is further considered via a lapse rate. The default setup in Nix accounts for adiabatic conditions (${\mathrm{\Gamma }}_{\mathrm{a}}=-\mathrm{9.8}\mathrm{°}\mathrm{C}{\mathrm{km}}^{-\mathrm{1}}$), though a moist lapse rate is also available (${\mathrm{\Gamma }}_{\mathrm{m}}=-\mathrm{5.0}\mathrm{°}\mathrm{C}{\mathrm{km}}^{-\mathrm{1}}$) (e.g. Stone and Carlson, 1979).

The particular atmospheric forcing is imposed at the sea level T_air(t), as shown in Fig. 4a. Starting from warm conditions T_air=0 °C, it reaches a minimum value of −30 °C in gradual steps that last 40 kyr each to ensure thermal quasi-equilibrium. Nonetheless, lower temperatures are present near the ice divide as the surface extends far above the sea level, wherein the lapse rate correction becomes relevant. This experiment reflects the insulating effect of the ice sheet as the forcing eventually reaches the initial atmospheric temperature but the grounding line does not retreat (Fig. 4b and f). It is not possible to make a one-to-one comparison (Schoof, 2007a), given the different physical descriptions of the system. Nonetheless, it is illustrative to represent the grounding-line position as a function of the ice temperature therein, evaluated at two different depths, and the vertical mean (Fig. 4b). The near-base temperature closely matches that of the theoretical prediction by the boundary layer. The lower layers appear to be shifted to the right, as the effect of the warmer surface becomes relevant.

Figure 4 Overdeepened bed experiment forced with atmospheric temperatures. (a) Forcing time series: atmospheric temperatures T_air(t). Each black symbol represents three snapshots of particular interest: the initial state (triangle, warmest conditions), coldest forcing conditions (star) and final state (square, same exact atmospheric conditions as the beginning). (b) Grounding-line position as a function of the ice temperature evaluated at two different depths (near the base and surface) and vertical mean. Note that the initial and end states differ strongly in ice extent even though the atmospheric forcing is identical. The grey line represents analytical results in the absence of thermodynamics (i.e. imposed rate factor A), following Schoof (2007a).

In terms of the hysteresis behaviour, the jump over the retrograde region of the bed geometry occurs for a near-base temperature of −30 °C, as predicted by the semi-analytical counterpart (see solid grey line, Fig. 4). Even so, when the forcing returns to the initial value, the grounding line does not retreat back to its original position and remains advanced (black square in Fig. 4b). In contrast, for shallower ice layers, the warming branch extends far from the analytical results, thus showing larger bistability against atmospheric temperature changes. As for the near-base layer, when the initial forcing state is eventually retrieved, the ice sheet extends beyond the bedrock peak and does not retreat. This is also well captured in Fig. 5 by comparing panels b and f, knowing that both are equilibrium states with identical forcing. These results differ strongly when using a downward-sloping bed geometry, as shown in Fig. 5a, c and e, where thermodynamics is also active but hysteresis is not present, as the bedrock does not present retrograde regions.

Figure 5 Ice-sheet extent and temperature distribution for the prograde (a, c, e) and the overdeepened (b, d, f) bed geometries. Each row represents a snapshot given by each symbol in Fig. 4 (triangle, star and square, respectively): the initial warm state (a, b), coldest atmospheric conditions (c, d) and final atmospheric configuration (e, f). Note that the overdeepening bed geometry yields a final ice-sheet profile extended far beyond the initial state even though the boundary conditions are identical, thus exhibiting hysteresis. Colours indicate the ice temperature at the given time.

Lastly, it is worth noting how close the temperature is to the melting point, particularly on the right-hand side of the base in all panels due to the combined contribution of geothermal heat flux and frictional heat dissipation (Fig. 5a, b, e and f), favoured by a warmer atmospheric temperature. This is unlike Fig. 5c and d, where the lower surface temperature perturbs the entire temperature profile, thus only partially cooling the ice-sheet base. It must be stressed that near the grounding line, there is a reduction in the basal temperature as a result of considerably thinner ice, thus providing less thermal insulation from the colder surface.

6.2.2 Ocean temperature anomalies ΔT_oce forcing

In this configuration, we apply ocean temperature anomalies with respect to a reference value T₀ so that $\mathrm{\Delta }{T}_{\text{oce}}=T-T_{\mathrm{0}}$. We then convert these temperature anomalies into sub-shelf melting at the grounding line (e.g. Favier et al., 2019) by using any of the parameterizations shown in Sect. 3.8.

We first perform two identical hysteresis experiments forced by ΔT_oce that solely differ on the thermodynamic treatment of the ice: an idealized fixed ice rate factor A (Fig. 6, blue curve) and a more realistic active thermodynamic scenario with a constant boundary condition T_air (Fig. 6, red curve). Results show that active thermodynamics considerably widens the width of the hysteresis loop. This behaviour resembles that obtained for the atmospherically forced simulations (Fig. 4), where we find a larger extent of the cooling branch compared to that of the semi-analytical solutions (grey line, Fig. 4).

Figure 6 (a) External forcing time series: ocean temperature anomalies ΔT_oce(t). The time duration of each step equals 30 kyr. (b) Hysteresis experiments for the overdeepened bed geometry forced via slowly varying ocean temperature anomalies ΔT_oce(t). Blue: constant ice rate factor $A={\mathrm{10}}^{-\mathrm{26}}{\mathrm{Pa}}^{\mathrm{3}}\mathrm{s}$. Red: active thermodynamics A=f(T) with fixed boundary condition $T_{\text{air}}=-\mathrm{40}\mathrm{°}\mathrm{C}$. Each forcing step is run for 30 kyr to ensure quasi-equilibrium (solid dots). A quadratic sub-shelf parameterization is employed in both scenarios. Heat exchange velocity parameter $\mathit{\gamma }={\mathrm{10}}^{-\mathrm{3}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$.

In addition to the experiments carried out to assess the importance of thermodynamics for the hysteresis behaviour of a marine-terminating ice sheet, we have performed a sensitivity study to quantify the differences caused by parameter uncertainty (e.g. Favier et al., 2019), particularly with respect to the heat exchange velocity γ.

Figure 7 illustrates the high sensitivity that stems from the heat exchange velocity parameter γ. It is worth noting that the retreat is much more sensitive to the particular γ choice than the later advance as the anomalies approach zero. That is, all intermediate values advance at $\mathrm{\Delta }{T}_{\text{oce}}=+\mathrm{1.5}\mathrm{°}\mathrm{C}$. In contrast, the retreat occurs for a wider range of temperature anomalies from +4.5 to +6.5 °C for $\mathit{\gamma }=\mathrm{2.1}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$, respectively.

Figure 7 Sensitivity tests. As in Fig. 6a but for (a) the constant ice rate factor $A={\mathrm{10}}^{-\mathrm{26}}{\mathrm{Pa}}^{\mathrm{3}}\mathrm{s}$. (b) Active thermodynamics with fixed T_air. Values of the heat exchange velocity parameter γ are given in ${\mathrm{10}}^{-\mathrm{5}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ and fall within the spanned range in Favier et al. (2019).

For a quadratic sub-shelf parameterization (Eq. 22), the retreat takes place at $\mathrm{\Delta }{T}_{\text{oce}}=+\mathrm{2.5}\mathrm{°}\mathrm{C}$ for the highest heat exchange velocity calibrated in Favier et al. (2019), i.e. $\mathit{\gamma }=\mathrm{100}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. Even the lowest parameter value presented in the same work also shows a retreat if the ocean temperature anomalies reach $\mathrm{\Delta }{T}_{\text{oce}}=+\mathrm{6.5}\mathrm{°}\mathrm{C}$. Multiple different values of γ advance back at nearly the same particular forcing value ΔT_oce 7, whereas the retreat happens at significantly different values. To illustrate this, we can take the hysteresis loops corresponding to $\mathit{\gamma }=\mathrm{20}\times {\mathrm{10}}^{-\mathrm{5}}$, $\mathrm{25}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{35}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. They retreat at $\mathrm{\Delta }{T}_{\text{oce}}=+\mathrm{5.5}$, +5.0 and +4.5 °C, respectively, whereas the advance takes place at precisely the same anomaly value $\mathrm{\Delta }{T}_{\text{oce}}=+\mathrm{2.5}\mathrm{°}\mathrm{C}$.

7 Model scalability and performance

This section illustrates the overall performance achieved by the Nix model in terms of computation speed for the MISMIP benchmark experiments (Pattyn et al., 2012). Moreover, both strong scalability and weak scalability are tested for the same experimental setup under the most sophisticated solver available in the model: the Blatter–Pattyn stress balance.

Nix's computational speed is assessed for each of the three velocity solvers available: SSA, DIVA and Blatter–Pattyn. An additional convergence study is detailed in Appendix B. Our domain is now restricted to a linear prograde slope for simplicity (i.e. Experiment 1 in Pattyn et al., 2012). Horizontal resolution is gradually increased, reaching 0.05 km.

Figure 8 shows the dependency of computational speed on the total number of horizontal grid points. The SSA and DIVA solvers exhibit a similar performance for all resolutions, whereas the higher-order Blatter–Pattyn solver is nearly 2 orders of magnitude slower. This is simply a result of the corresponding sparse-matrix solution in the latter, which is significantly slower than the tridiagonal solver applied for both SSA and DIVA.

Figure 8 Nix's computational speed for the three solvers available. The double x axis represents the number of horizontal grid points n and the corresponding resolution at the grounding line Δx. Note that Nix allows for an unevenly spaced stretched grid that explicitly tracks the grounding-line position.

It is challenging to give a one-to-one comparison with other existing models, given the differences in dimensionality. Generally, other ice-sheet models solving for higher-order momentum balance coupled with active thermodynamics are comprehensive three-dimensional models. To give an estimation, the MALI ice-sheet model (Hoffman et al., 2018) control simulations average 5.26 simulated years per wall-clock hour. In contrast, MISMIP experiments run with Nix reach ∼10⁵ simulated years per wall-clock hour on average. Thus, there is a 5 order of magnitude difference in terms of computational time.

Nix's parallel implementation is particularly of interest for the Blatter–Pattyn stress balance. The velocity solution relies on the biconjugate gradient stabilized (BiCGSTAB) algorithm, where multi-threading can be exploited if the model is compiled with OpenMP enabled. Nix employs a row-major sparse-matrix format to achieve the best performance. By default, the tolerance and maximum number of iterations over the sparse linear system are set to $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{6}}$ and N=10³, respectively. Note that ϵ should not be mistaken for the non-linear Picard convergence criterion ϕ_tol related to viscosity convergence. Tolerance values above 10⁻⁴ give rise to numerical instabilities, which hamper viscosity convergence.

As a synthetic test for scalability, we restrict to a linear prograde slope (i.e. Experiment 1 in Pattyn et al., 2012) and assess both the strong and the weak scalability of the model. The former fixes the problem size (total number of grid points) while increasing the number of threads (CPUs) available to the sparse solver. For the latter, there is a constant CPU load for different grid sizes. To illustrate the strong scalability, an evenly spaced grid with 10⁵ horizontal points and 10³ vertical layers is considered. The necessary RAM to run such a job is ∼54 GB. The results of scalability are characterized by the speed-up $t_{\mathrm{1}}/t_n$, where t_n represents the time required to complete the job using n threads (e.g. Gagliardini et al., 2013). Equivalently, the efficiency is defined as $t_{\mathrm{1}}/\left(n{t}_n\right)$.

The results of the strong scalability are shown in Fig. 9. Parallelization through Eigen library multi-threading allows for an overall acceleration factor of ∼2. This result is further tested under permutations of varying optimization levels during compilation with the OpenMP flag (i.e. O1, O2 and O3) and the maximum number of iterations over the linear problem N=10¹, 10² and 10³. Minor differences are found among experiments with different optimization levels, in agreement with other more complex models such as ISSM (Fischler et al., 2022). Irrespective of both the particular number of maximum iterations and the optimization choice, the efficiency rapidly drops as the number of threads increases.

Figure 9 Acceleration and efficiency for strong-scalability experiments. The maximum number of iterations N in the sparse linear problem is given as a coloured legend. Line styles denote the three levels of optimization provided by OpenMP during compilation (O1, O2 and O3, in increasing order).

Analogously, Nix's efficiency under weak scalability is depicted in Fig. 10. Beyond four threads, the performance reaches values below 25 %. There are slight differences among optimization levels during parallel compilation. Higher levels, such as O2 and O3, appear to outperform O1 as the total number of linear iterations increases, simply as a result of a larger total number of computations.

Figure 10 Efficiency for weak-scalability experiments. The maximum number of iterations N in the sparse linear problem is given as a coloured legend, ranging from 10¹ to 10⁶. Line styles denote the three levels of optimization provided by OpenMP during compilation (O1, O2 and O3, in increasing order).

In view of the scaling results (Figs. 9 and 10), two conclusions can be drawn. First, the serial execution of Nix's problem is extremely well optimized for the problem size of interest. This further supports the applicability of Nix as a lightweight model that can run at extremely high resolution with minimum computational resources. Second, users can reduce the simulation wall time by a factor of ∼2 by setting eight threads (CPU) and a conservative optimization level provided by OpenMP (i.e. O1 flag).

8 Discussion

The results of the MISMIP benchmark experiments are successful given the good agreement between our numerical solution and the semi-analytical work of Schoof (2007a) (Fig. 3). From a modelling perspective, our grounding-line position is slightly shifted upstream, like in other moving-grid models shown in Pattyn et al. (2012). Nevertheless, a test of sensitivity to spatial resolution shows an asymptotic convergence toward the semi-analytical solution (Fig. B1, Appendix B), thus providing robustness to our results. A further comparison among the Nix velocity solvers shows excellent agreement on the equilibrium solutions for both bed geometries studied herein: downward sloping and overdeepening.

For active thermodynamics and air temperature forcing, the corresponding temperature range spanned by the MISMIP ice rate factor A(T) does not yield a full advance/retreat of the ice sheet. This can be understood by the insulating effect of the ice sheet. The MISMIP idealized forcing with varying rate factors simultaneously modifies the ice viscosity over the entire domain, whereas the real temperatures given by an active thermomechanical solver will adjust to the new surface temperature (Fig. 4). This further means that the impact on the viscosity is weaker as there are other heat sources, such as the basal friction dissipation and the geothermal heat flow.

The stability of the system is accordingly perturbed, as shown in Fig. 6. In particular, the bistability of the system is increased, in the sense that a larger range of perturbation values has two stable solutions. To illustrate this, for a rate factor $A={\mathrm{10}}^{-\mathrm{26}}{\mathrm{Pa}}^{\mathrm{3}}\mathrm{s}$, the oceanic anomaly perturbation ranges from $\mathrm{\Delta }{T}_{\text{oce}}=+\mathrm{3}$ to +7 °C, whereas for the thermomechanically active scenario (i.e. temperature-dependent viscosity), the bistability solution is found for a wider range from $\mathrm{\Delta }{T}_{\text{oce}}=+\mathrm{1.0}$ to +8.0 °C. This can be understood as the result of a thermal adjustment that occurs when the temperature of the ice can evolve in time. As the forcing changes over time (i.e. the ocean temperature anomalies), the ice flux at the grounding line is modified and the inland ice thickness is perturbed accordingly. The new ice thickness distribution implies a slightly different temperature solution, and the viscosity is consequently modified. Knowing that the viscosity field determines the velocity solution via the stress balance, we therefore find a clear feedback that allows the ice sheet to adjust if thermodynamics are active.

It must be noted that this stability study employs a time duration in each forcing step of 30 kyr to elude transient responses (Fig. 6b). This allows the ice geometry to reach a steady state as these experiments are intended to be quasi-equilibrium simulations (similar to MISMIP). If the time between steps was reduced, the hysteresis loop would then be a transitory response given that the ice temperature may not adjust to the new geometry. The added value of the thermomechanical coupling relies not only on the transitory response (i.e. thermal inertia), but also on the perturbed stability of the quasi-equilibrium hysteresis loop. In other words, the thermomechanical coupling already determines the stable regions in a quasi-steady description and not only through the effects of thermal inertia. Our goal here is to first show this more fundamental behaviour. Further work is needed to assess the relevance of potential transient responses.

It is worth noting the fundamental role of vertical advection if the system is to be forced with air temperatures. The magnitude of vertical advection and its vertical dependency determine the temperature distribution within the ice (Moreno-Parada et al., 2024). This mechanism strongly determines the equilibrium profiles as it modifies the overall outflow of ice via the viscosity dependency on temperatures (i.e. Arrhenius law in Eq. 17).

We must also discuss the calibrated values of γ_T (Favier et al., 2019). The impact on the hysteresis behaviour is interesting as it does not imply a symmetric effect on the retreat/advance of the ice sheet (Fig. 7). Strictly speaking, the temperature anomaly necessary for the ice sheet to retreat is far more sensitive to the particular γ value than the anomaly necessary to advance. One potential explanation for this interesting behaviour is the very nature of the melting/calving parameterization at the grounding line. Since M grows with ΔT_oce (Eqs. 21 and 22), the difference among M values for a fixed γ increases with the temperature anomaly. Hence, knowing that the advance occurs when the ice flux value reaches a certain value at the grounding line, the temperature anomaly range covered by different γ that yields such melting/calving values is smaller (as ΔT_oce→0). Moreover, if thermodynamics is considered (Fig. 7b), the required calving at the grounding line to retreat is generally larger than that for a fixed rate factor (Fig. 7b). In fact, for sufficiently low values of γ, the thermomechanically active ice sheet never retreats. The thermal behaviour of the ice thus provides additional inertia in the sense that the ice sheet is less prone to changing its current state, thus endowing the system with higher stability. These results are not exclusive of the oceanic forcing, as also shown in Fig. 4, where the ice sheet does not collapse to a retreated position when the perturbation vanishes.

Besides Nix's extremely fast computations, marginal parallel speed-up is noticeable for the Blatter–Pattyn solver (Figs. 9 and 10). Given the scalability analysis in Sect. 7, this result is consistent for all permutations of the total maximum number of linear iterations and the three different levels of OpenMP optimization during compilation. There are a number of potential causes underlying this behaviour. First, Nix's most complicated calculations are carried out by the velocity solver, and it relies on the BiCGSTAB algorithm implemented in the Eigen library. Such a method allows for multi-threading if OpenMP is enabled during compilation, and it has certain limitations concerning sparse-matrix–vector multiplication. This operation is fundamental knowing that the Eigen library performs a static partitioning of the sparse matrix. That is, rows are divided among the available threads leveraged by OpenMP, and each thread then processes a contiguous block of rows to minimize the synchronization overhead. Since the number of non-zero entries per row remains constant, the partitioning approach seemed justified. Nevertheless, Eigen does not allow for more complex partitioning strategies, such as graph-based approaches (Catalyurek and Aykanat, 1999), corner partitioning (Wolf et al., 2008) or block-cyclic partitioning (e.g. Aboelaze et al., 1991; Petitet and Dongarra, 1999), to minimize the total communication volume while keeping the computation balanced across compute nodes.

A second potential cause of marginal parallel performance concerns overhead. Since the BiCGSTAB method computes two matrix–vector operations per iteration, there is an increased communication overhead compared to that of simpler solvers. For small- to medium-sized problems, the overhead of creating/managing threads may negate some potential speed-up. This is of particular relevance for the problem at hand: a two-dimensional ice-sheet model. The Blatter–Pattyn stress balance yields a two-dimensional velocity field that is considerably less computationally expensive than the three-dimensional counterpart. To test the overhead hypothesis, we repeated the scalability performance test (not shown here) with an algorithm that only computes one matrix–vector multiplication per iteration: the conjugate gradient method (e.g. Shewchuk, 1994). Results show that both methods have similar efficiency when parallelization is enabled, suggesting that overhead is not the main cause of marginal speed-up. This further supports the hypothesis that sparse-matrix–vector multiplication is the limiting component when enabling parallelization. We propose that Eigen simplicity of a row-wise partitioning of the sparse matrix is not enough to yield optimal parallelization. Future work is needed to explore more advanced approaches of the partitioning strategy.

Lastly, the results presented herein show that active thermodynamics perturbs the hysteresis loop and the overall stability of an ice sheet. The particular grid over which the equations are discretized does not alter this behaviour. Consequently, our results are not numerical artefacts of the selected mesh. Therefore, we do not expect our results to change for a different grid discretization. More precisely, we expect the same physical behaviour as long as the ice viscosity varies with temperature changes, irrespective of the chosen grid. The exact grounding-line position may differ for a different grid, yet the physical mechanism underlying this mechanism remains unperturbed. It is thus expected that other models with distinct meshes will exhibit a similar mechanism to that observed in Nix simulations.

9 Conclusions

The thermomechanically coupled two-dimensional model Nix has been presented and thoroughly described. There are a number of novelties compared to other two-dimensional models: a stress balance given by the Blatter–Pattyn approximation, a fully coupled thermodynamics solver, explicit calculation of the grounding line by a stretched coordinate system, adaptive time stepping and potential melting/calving at the grounding line. The Nix ice-sheet model is designed to study a complex system with minimum resources: within a few hours, a 100 m resolution higher-order stress balance coupled with full thermodynamics is feasible even on a regular laptop.

First, Nix's performance was tested by reproducing Experiments 1, 2 and 3 from MISMIP benchmarks (Pattyn et al., 2012). Results were further compared to semi-analytical solutions (Schoof, 2007a), yielding excellent agreement for all Stokes approximations available in Nix: SSA, DIVA and Blatter–Pattyn. In general, our grounding-line position slightly underestimates the boundary layer results in Schoof (2007a) and all moving-grid models participating in MISMIP. The well-known hysteresis behaviour in Experiment 3 is also captured.

The complexity of the system was further increased by solving the associated heat problem. This allows us to investigate to which extent the hysteresis behaviour under parameter variations is perturbed as a result of a temperature-dependent viscosity. In so doing, we designed two different suites of experiments regarding the forcing variables of the system: air temperatures T_air and ocean temperature anomalies ΔT_oce.

When forcing with air temperatures, the hysteresis loop width is widened and the system exhibits larger bistability as the dynamics are perturbed via the ice viscosity dependency on temperature. This is a direct result of the insulating effect of the ice, thus preventing the colder atmospheric temperatures from travelling downward from the uppermost layers. In an idealized overdeepened bed geometry, it is necessary to reach an air temperature of −40 °C at sea level (provided an adiabatic lapse rate dependency with height is present) for the grounding line to advance beyond the bedrock local maximum.

If the system is instead forced by ocean temperature anomalies (i.e. melting/calving at the grounding line), we find that the hysteresis behaviour also persists. Notably, the ocean temperature anomaly at which the ice sheet retreats depends on the particular heat exchange parameter. Our results show that the quadratic parameterization retreats at $\mathrm{\Delta }{T}_{\text{oce}}=+\mathrm{8.0}\mathrm{°}\mathrm{C}$ of the temperature anomaly. The system advances back to its unperturbed state at $\mathrm{\Delta }{T}_{\text{oce}}=+\mathrm{1.0}\mathrm{°}\mathrm{C}$ for the quadratic parameterization.

These results show that a temperature-dependent ice viscosity provides inertia and stability to the ice sheet, regardless of the particular external forcing applied (i.e. oceanic or atmospheric). Therefore, the thermal state of the ice must be accounted for, even for low-dimensional or conceptual models, particularly if the potential sudden retreat of an ice sheet is to be assessed.

Appendix A: Discretization schemes and nonuniform grids

We herein elaborate on a thorough description of the discretization schemes of our flowline model where we follow the ordinary notation $q({\mathit{\sigma }}_i,{\mathit{\zeta }}_j,{\mathit{\tau }}_n)\equiv q_{i,j}^n$.

The position in the spatio-temporal coordinates is then given by ${\mathit{\sigma }}_i={\sum }_{k=\mathrm{0}}^i\mathrm{\Delta }{\mathit{\sigma }}_k$, ${\mathit{\zeta }}_i={\sum }_{k=\mathrm{0}}^i\mathrm{\Delta }{\mathit{\zeta }}_k$ and τ_n=nΔτ. Note that Nix allows for a nonuniform spatial grid where the spacing between two consecutive points follows a desired distribution (polynomial or exponential). This yields high resolutions near the grounding line while minimizing the total number of grid points. The horizontal index $i\in {\mathcal{W}}_r=\mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },r\mathit{\}}$, where r is the number of points in which the horizontal axis is divided. Likewise, the vertical index follows $j\in {\mathcal{W}}_p=\mathit{\{}\mathrm{0},\mathrm{1},\mathrm{2},\mathrm{\dots },p\mathit{\}}$, where p is the number of vertical layers.

Nix's finite-difference discretization considers unevenly spaced grids, commonly used in the glaciological community, where higher resolutions are desired near the base, while minimizing the required number of points to reduce computational costs. A new coordinate system is thus built considering two types of nonuniform grid spacing: polynomial and exponential. For any normalized variable ξ (such as σ and ζ, Eq. 23), a new nonuniform grid $\stackrel{\mathrm{̃}}{\mathit{\xi }}$ can be obtained by a power law transformation:

$$\begin{array}{}\text{(A1)}& \stackrel{\mathrm{̃}}{\mathit{\xi }}={\mathit{\xi }}^n,\end{array}$$

where n is the spacing order, or an exponential map:

$$\begin{array}{}\text{(A2)}& \stackrel{\mathrm{̃}}{\mathit{\xi }}={\displaystyle \frac{e^{s\mathit{\xi }}-\mathrm{1}}{e^s-\mathrm{1}}},\end{array}$$

where s is the spacing factor for the exponential grid and $\stackrel{\mathrm{̃}}{\mathit{\xi }}$ is substituted by the variables σ and ζ. In this study, we have employed n=2 and s=2.

A1 Blatter–Pattyn stress balance discretization

The discretization is straightforward for an Arakawa C grid. The position of the grounding line L(t) is located on the velocity grid (following Vieli and Payne, 2005). Thus, if the horizontal axis is divided into r points, the ice thickness grid ranges between $i=\mathrm{0},\mathrm{1},\mathrm{\dots },r-\mathrm{1}$, whereas the velocity grid (staggered) indexes read $i=\mathrm{1}/\mathrm{2},\mathrm{3}/\mathrm{2},\mathrm{\dots },r-\mathrm{1}/\mathrm{2}$. The fractional index implies that the point $(i+\mathrm{1}/\mathrm{2},j)$ lies between (i,j) and $(i+\mathrm{1},j)$ and analogously for the vertical index j.

The Blatter–Pattyn stress balance can be written as

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{2}}{L^{\mathrm{2}}\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}}}\left[\begin{array}{c}{\mathit{\eta }}_{i+\mathrm{1},j}{\displaystyle \frac{u_{i+\mathrm{3}/\mathrm{2},j}-u_{i+\mathrm{1}/\mathrm{2},j}}{\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{3}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}}}\\ -{\mathit{\eta }}_{i,j}{\displaystyle \frac{u_{i+\mathrm{1}/\mathrm{2},j}-u_{i-\mathrm{1}/\mathrm{2},j}}{\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i-\mathrm{1}/\mathrm{2}}}}\end{array}\right]+\\ & {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\left(H_i\right)}^{\mathrm{2}}\mathrm{\Delta }{\mathit{\zeta }}_{j+\mathrm{1}/\mathrm{2}}}}\left[\begin{array}{c}{\mathit{\eta }}_{i,j+\mathrm{1}}{\displaystyle \frac{u_{i,j+\mathrm{3}/\mathrm{2}}-u_{i,j+\mathrm{1}/\mathrm{2}}}{\mathrm{\Delta }{\mathit{\zeta }}_{j+\mathrm{3}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\zeta }}_{j+\mathrm{1}/\mathrm{2}}}}\\ -{\mathit{\eta }}_{i,j}{\displaystyle \frac{u_{i,j+\mathrm{1}/\mathrm{2}}-u_{i,j-\mathrm{1}/\mathrm{2}}}{\mathrm{\Delta }{\mathit{\zeta }}_{j+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\zeta }}_{j-\mathrm{1}/\mathrm{2}}}}\end{array}\right]=\\ & \mathit{\rho }g{\displaystyle \frac{h_{i+\mathrm{1}}-h_i}{L\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}}}.\end{array}\end{array}$$

We thus have a linear system of r×p unknowns that can be solved by applying standard linear algebraic solvers. For each time step, we build a matrix of coefficients with dimension (rp)×(rp):

$$\begin{array}{}\text{(A2)}& \underset{\left(rp\right)\times \left(rp\right)}{\mathrm{A}}\cdot \underset{\left(rp\right)\times \left(\mathrm{1}\right)}{\mathrm{u}}=\underset{\left(rp\right)\times \left(\mathrm{1}\right)}{\mathrm{F}}.\end{array}$$

Since our discretization stencil includes six points, $(i+\mathrm{3}/\mathrm{2},j)$, $(i+\mathrm{1}/\mathrm{2},j)$, $(i-\mathrm{1}/\mathrm{2},j+\mathrm{1})$, $(i,j+\mathrm{3}/\mathrm{2})$, $(i,j+\mathrm{1}/\mathrm{2})$ and $(i,j-\mathrm{1}/\mathrm{2})$, we obtain a sparse matrix that allows for optimized inversion. For r=500 and p=25, only 0.048 % of the coefficient matrix is a nonzero entry:

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}& \left(\begin{array}{cccccccccc}& & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & & \\ {\mathit{\alpha }}_{i-\mathrm{2}M,j}& \mathrm{\dots }& {\mathit{\alpha }}_{i-M,j}& \mathrm{\dots }& {\mathit{\alpha }}_{i,j-\mathrm{1}/\mathrm{2}}& {\mathit{\alpha }}_{i,j+\mathrm{1}/\mathrm{2}}& \mathrm{\dots }& {\mathit{\alpha }}_{i+M,j}& \mathrm{\dots }& {\mathit{\alpha }}_{i+\mathrm{2}M,j}\\ & & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & & \\ & & & & & & & & & \end{array}\right)\\ & \cdot \left(\begin{array}{c}{u}_{i-\mathrm{3}/\mathrm{2},j}\\ \mathrm{⋮}\\ u_{i-\mathrm{1}/\mathrm{2},j}\\ \mathrm{⋮}\\ u_{i,j-\mathrm{1}/\mathrm{2}}\\ u_{i,j+\mathrm{1}/\mathrm{2}}\\ \mathrm{⋮}\\ u_{i+\mathrm{1}/\mathrm{2},j}\\ \mathrm{⋮}\\ u_{i+\mathrm{3}/\mathrm{2},j}\end{array}\right)=\left(\begin{array}{c}{F}_{i-\mathrm{3}/\mathrm{2},j}\\ \mathrm{⋮}\\ F_{i-\mathrm{1}/\mathrm{2},j}\\ \mathrm{⋮}\\ F_{i,j-\mathrm{1}/\mathrm{2}}\\ F_{i,j+\mathrm{1}/\mathrm{2}}\\ \mathrm{⋮}\\ F_{i+\mathrm{1}/\mathrm{2},j}\\ \mathrm{⋮}\\ F_{i+\mathrm{3}/\mathrm{2},j}\end{array}\right).\end{array}\end{array}$$

The boundary conditions are then needed to evaluate the velocity at the edge of the domain. For the free surface we apply the following three-point derivative scheme:

$$\begin{array}{}\text{(A4)}& u_{i,r-\mathrm{1}/\mathrm{2}}={\displaystyle \frac{\mathrm{4}{u}_{i,r-\mathrm{3}/\mathrm{2}}-u_{i,r-\mathrm{5}/\mathrm{2}}+\mathrm{2}\mathit{\chi }\mathrm{\Delta }{\mathit{\zeta }}_{r-\mathrm{3}/\mathrm{2}}}{\mathrm{3}}},\end{array}$$

where ξ is evaluated right below the uppermost layer (i.e. $j=r-\mathrm{3}/\mathrm{2}$):

$$\begin{array}{}\text{(A5)}& \mathit{\chi }=\mathrm{4}{\displaystyle \frac{u_{i+\mathrm{1}/\mathrm{2},r-\mathrm{3}/\mathrm{2}}-u_{i-\mathrm{1}/\mathrm{2},r-\mathrm{3}/\mathrm{2}}}{\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i-\mathrm{1}/\mathrm{2}}}}{\displaystyle \frac{h_{i+\mathrm{1}}-h_i}{\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}}}.\end{array}$$

The basal boundary condition is thus obtained analogously with a non-zero drag. In this case, we discretize with the following simple two-point scheme:

$$\begin{array}{}\text{(A6)}& u_{i,\mathrm{1}/\mathrm{2}}=u_{i,\mathrm{3}/\mathrm{2}}-\mathit{\mu }\mathrm{\Delta }{\mathit{\zeta }}_{i,\mathrm{1}/\mathrm{2}},\end{array}$$

where μ is computed right above the basal boundary (i.e. $j=\mathrm{3}/\mathrm{2}$):

$$\begin{array}{}\text{(A7)}& \mathit{\mu }=\mathrm{4}{\displaystyle \frac{u_{i+\mathrm{1}/\mathrm{2},\mathrm{3}/\mathrm{2}}-u_{i-\mathrm{1}/\mathrm{2},\mathrm{3}/\mathrm{2}}}{\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i-\mathrm{1}/\mathrm{2}}}}{\displaystyle \frac{b_{i+\mathrm{1}}-b_i}{\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}}}.\end{array}$$

With these expressions, the velocity field is then complete. Lastly, at the ice divide and the grounding line, the velocity is obtained from symmetry and hydrostatic equilibrium arguments, respectively (elaborated on in Appendix A2).

A2 DIVA/SSA stress balance discretization

The discretization is straightforward for a staggered grid. The position of the grounding line L(t) is located on the staggered grid (following Vieli and Payne, 2005). Thus, if the domain is divided into n points, the ice thickness grid ranges between $i=\mathrm{0},\mathrm{1},\mathrm{\dots },r-\mathrm{1}$, whereas the velocity grid (staggered) indexes read $i=\mathrm{1}/\mathrm{2},\mathrm{3}/\mathrm{2},\mathrm{\dots },r-\mathrm{1}/\mathrm{2}$.

The SSA stress balance can then be written as¹

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{2}}{L^{\mathrm{2}}\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}}}\\ & \times \left[{\mathit{\eta }}_{i+\mathrm{1}}{H}_{i+\mathrm{1}}{\displaystyle \frac{u_{i+\mathrm{3}/\mathrm{2},j}-u_{i+\mathrm{1}/\mathrm{2},j}}{\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{3}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}}}-{\mathit{\eta }}_i{H}_i{\displaystyle \frac{u_{i+\mathrm{1}/\mathrm{2},j}-u_{i-\mathrm{1}/\mathrm{2},j}}{\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i-\mathrm{1}/\mathrm{2}}}}\right]\\ & -{\mathit{\beta }}_{i+\mathrm{1}/\mathrm{2}}^{\mathrm{2}}{u}_{i+\mathrm{1}/\mathrm{2}}\\ & =\mathit{\rho }g{\displaystyle \frac{H_{i+\mathrm{1}}+H_i}{\mathrm{2}}}{\displaystyle \frac{h_{i+\mathrm{1}}-h_i}{L\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}}},\end{array}\end{array}$$

where the friction coefficient β² reads

$$\begin{array}{}\text{(A9)}& {\mathit{\beta }}_{i+\mathrm{1}/\mathrm{2}}=C{u}_{i+\mathrm{1}/\mathrm{2}}^{m-\mathrm{1}}\end{array}$$

so that τ_b=β²u.

Assuming our domain is divided into n points, the corresponding tridiagonal matrix is built at every time step as (where we have dropped the super-index to lighten the notation)

$$\begin{array}{}\text{(A10)}& \left(\begin{array}{ccccc}{B}_{\mathrm{0}}& C_{\mathrm{0}}& & & \\ A_{\mathrm{1}}& B_{\mathrm{1}}& C_{\mathrm{1}}& & \\ & A_{\mathrm{2}}& B_{\mathrm{2}}& \mathrm{\ddots }& \\ & & \mathrm{\ddots }& \mathrm{\ddots }& C_{r-\mathrm{2}}\\ & & & A_{r-\mathrm{1}}& B_{r-\mathrm{1}}\end{array}\right)\left(\begin{array}{c}{u}_{\mathrm{1}/\mathrm{2}}\\ u_{\mathrm{3}/\mathrm{2}}\\ \mathrm{⋮}\\ u_{r-\mathrm{3}/\mathrm{2}}\\ u_{r-\mathrm{1}/\mathrm{2}}\end{array}\right)=\left(\begin{array}{c}{F}_{\mathrm{0}}\\ F_{\mathrm{1}}\\ \mathrm{⋮}\\ F_{r-\mathrm{2}}\\ F_{r-\mathrm{1}}\end{array}\right),\end{array}$$

where $A_{\mathrm{0}}=C_{n-\mathrm{1}}=\mathrm{0}$ by definition.

The non-zero entries of the matrix and the inhomogeneous term read

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}{\displaystyle }{A}_i& {\displaystyle }={\mathit{\gamma }}_i{\mathit{\eta }}_i{H}_i,\\ {\displaystyle }{B}_i& {\displaystyle }=-{\mathit{\gamma }}_i\left({\mathit{\eta }}_{i+\mathrm{1}}{H}_{i+\mathrm{1}}+{\mathit{\eta }}_i{H}_i\right)-{\mathit{\beta }}_i^{\mathrm{2}},\\ {\displaystyle }{C}_i& {\displaystyle }={\mathit{\gamma }}_i{\mathit{\eta }}_{i+\mathrm{1}}{H}_{i+\mathrm{1}},\\ {\displaystyle }{F}_i& {\displaystyle }=\mathit{\rho }g{\displaystyle \frac{H_{i+\mathrm{1}}+H_i}{\mathrm{2}}}{\displaystyle \frac{h_{i+\mathrm{1}}-h_i}{L\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}}},\end{array}\end{array}$$

where ${\mathit{\gamma }}_i=\mathrm{2}/(L\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}{)}^{\mathrm{2}}$.

For the edge of the matrix (i.e. $i=n-\mathrm{1}$), we use the following values:

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}{\displaystyle }{A}_{r-\mathrm{1}}& {\displaystyle }={\mathit{\eta }}_{r-\mathrm{1}}{H}_{r-\mathrm{1}},\\ {\displaystyle }{B}_{r-\mathrm{1}}& {\displaystyle }=-{\mathit{\gamma }}_{r-\mathrm{1}}\left({\mathit{\eta }}_{r-\mathrm{1}}{H}_{r-\mathrm{1}}\right)-{\mathit{\beta }}_{n-\mathrm{1}}^{\mathrm{2}},\\ {\displaystyle }{C}_{r-\mathrm{1}}& {\displaystyle }=\mathrm{0},\\ {\displaystyle }{F}_{r-\mathrm{1}}& {\displaystyle }=\mathit{\rho }g{H}_{r-\mathrm{1}}{\displaystyle \frac{h_{r-\mathrm{1}}-h_{r-\mathrm{2}}}{L\mathrm{\Delta }{\mathit{\sigma }}_{r-\mathrm{3}/\mathrm{2}}}}.\end{array}\end{array}$$

For the boundary values, we set²

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}{\displaystyle }{u}_{\mathrm{1}/\mathrm{2}}& {\displaystyle }=-u_{\mathrm{3}/\mathrm{2}},\\ {\displaystyle }{u}_{r-\mathrm{1}/\mathrm{2}}& {\displaystyle }=u_{r-\mathrm{3}/\mathrm{2}}+{\displaystyle \frac{L\mathrm{\Delta }{\mathit{\sigma }}_{r-\mathrm{3}/\mathrm{2}}}{\mathrm{8}{\mathit{\eta }}_{r-\mathrm{1}}}}\left(\mathit{\rho }g{H}_{r-\mathrm{1}}^{\mathrm{2}}-{\mathit{\rho }}_{\mathrm{w}}g{D}^{\mathrm{2}}\right),\end{array}\end{array}$$

where D is the bedrock depth below sea level, the first equality follows from symmetry arguments at the ice divide (i=1) and the second equality implies that the ice momentum is equated by the hydrostatic pressure of the water.

A3 Advection discretization

For the advection equation we also chose an implicit scheme for numerical stability:

$$\begin{array}{}\text{(A14)}& \begin{array}{rl}{\displaystyle \frac{H_i^{n+\mathrm{1}}-H_i^n}{\mathrm{\Delta }{\mathit{\tau }}^n}}& ={\mathit{\sigma }}_i{\displaystyle \frac{{\dot{L}}^n}{L^n}}{\displaystyle \frac{H_{i+\mathrm{1}}^{n+\mathrm{1}}-H_{i-\mathrm{1}}^{n+\mathrm{1}}}{\left(\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i-\mathrm{1}/\mathrm{2}}\right)}}\\ & -{\displaystyle \frac{\mathrm{2}\left(q_{i+\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}-q_{i-\mathrm{1}/\mathrm{2}}^{n+\mathrm{1}}\right)}{L^n\left(\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i-\mathrm{1}/\mathrm{2}}\right)}}+S_i^n,\end{array}\end{array}$$

where the ice flux is defined as

$$\begin{array}{}\text{(A15)}& q_{i+\mathrm{1}/\mathrm{2}}=u_{i+\mathrm{1}/\mathrm{2}}{\displaystyle \frac{H_{i+\mathrm{1}}+H_i}{\mathrm{2}}}.\end{array}$$

However, at the grounding line we use

$$\begin{array}{}\text{(A16)}& q_{r-\mathrm{1}/\mathrm{2}}=u_{r-\mathrm{1}/\mathrm{2}}{H}_{r-\mathrm{1}}+M,\end{array}$$

where M represents a flux anomaly at the terminus driven by variable ocean conditions. In real glaciers, these anomalies could be driven by variable calving, submarine melt, or a combination of these factors.

The advection equation can be rewritten as

$$\begin{array}{}\text{(A17)}& A_i{H}_{i-\mathrm{1}}^{n+\mathrm{1}}+B_i{H}_i^{n+\mathrm{1}}+C_i{H}_{i+\mathrm{1}}^{n+\mathrm{1}}=F_i\end{array}$$

so that the corresponding matrix is also tridiagonal:

$$\begin{array}{}\text{(A18)}& \begin{array}{rl}{\displaystyle }{A}_i& {\displaystyle }={\mathit{\gamma }}_i^n\left({\mathit{\sigma }}_i\dot{L}-u_{i-\mathrm{1}}\right),\\ {\displaystyle }{B}_i& {\displaystyle }=\mathrm{1}+{\mathit{\gamma }}_i^n\left(u_i-u_{i-\mathrm{1}}\right),\\ {\displaystyle }{C}_i& {\displaystyle }={\mathit{\gamma }}_i^n\left(-{\mathit{\sigma }}_i\dot{L}+u_i\right),\\ {\displaystyle }{F}_i& {\displaystyle }=H_i^n+\mathrm{\Delta }{\mathit{\tau }}^n{S}_i^n,\end{array}\end{array}$$

where ${\mathit{\gamma }}_i^n=\mathrm{\Delta }{\mathit{\tau }}^n/\mathrm{2}(\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i-\mathrm{1}/\mathrm{2}})L^n$.

As the ice divide is located at i=1 (note that we start counting at i=0), the boundary condition then reads

$$\begin{array}{}\text{(A19)}& H_{\mathrm{0}}^n=H_{\mathrm{2}}^n,\end{array}$$

since σ=1 is a symmetry axis.

A4 Grounding-line scheme

The terminus position L (i.e. the grounding line) is not fixed in time. Direct discretization of Eq. (18) in terms of σ coordinates leads to

$$\begin{array}{}\text{(A20)}& \begin{array}{rl}{\dot{L}}^n& \equiv {\displaystyle \frac{\mathrm{d}L}{\mathrm{d}\mathit{\tau }}}\\ & ={\displaystyle \frac{\left\{\begin{array}{c}-L^n\mathrm{\Delta }{\mathit{\sigma }}_{r-\mathrm{1}/\mathrm{2}}{S}_{r-\mathrm{1}}^n\\ +\mathrm{2}\left(q_{r-\mathrm{1}/\mathrm{2}}^n-q_{r-\mathrm{3}/\mathrm{2}}^n\right)/\left(\mathrm{\Delta }{\mathit{\sigma }}_{r-\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{r-\mathrm{3}/\mathrm{2}}\right)\end{array}\right\}}{H_{r-\mathrm{1}}^n-H_{r-\mathrm{2}}^n+\mathit{\varrho }\left(b_{r-\mathrm{1}}^n-b_{r-\mathrm{2}}^n\right)}}.\end{array}\end{array}$$

A5 Thermodynamics discretization scheme

Unlike previous discretizations, the temperature field $\mathit{\theta }(\mathit{\sigma },\mathit{\zeta },\mathit{\tau })$ has an additional dependency on the vertical coordinate ζ that brings a higher degree of complexity (Eq. 24).

The energy balance (Eq. 28) is discretized using an upwind scheme with a forward Euler step and a centred difference for the spatial derivatives. The lengthy expression reads

$$\begin{array}{}\text{(A21)}& \begin{array}{rl}& \mathit{\rho }c\left[{\displaystyle \frac{{\mathit{\theta }}_{i,j}^{n+\mathrm{1}}-{\mathit{\theta }}_{i,j}^n}{\mathrm{\Delta }{\mathit{\tau }}^n}}-{\mathit{\sigma }}_i{\displaystyle \frac{{\dot{L}}^n}{L^n}}{\displaystyle \frac{{\mathit{\theta }}_{i+\mathrm{1},j}^n-{\mathit{\theta }}_{i-\mathrm{1},j}^n}{\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i-\mathrm{1}/\mathrm{2}}}}+\right.\\ & -\left.{\displaystyle \frac{{\mathit{\zeta }}_{i,j}}{H_i^n}}{\displaystyle \frac{H_i^{n+\mathrm{1}}-H_i^{n-\mathrm{1}}}{\mathrm{2}\mathrm{\Delta }{\mathit{\tau }}^n}}{\displaystyle \frac{{\mathit{\theta }}_{i,j+\mathrm{1}}^n-{\mathit{\theta }}_{i,j-\mathrm{1}}^n}{\mathrm{\Delta }{\mathit{\zeta }}_{i,j+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\zeta }}_{i,j-\mathrm{1}/\mathrm{2}}}}\right]\\ & ={\displaystyle \frac{k}{{\left(H_i^n\right)}^{\mathrm{2}}}}{\displaystyle \frac{{\mathit{\theta }}_{i,j+\mathrm{1}}^n-\mathrm{2}{\mathit{\theta }}_{i,j}^n+{\mathit{\theta }}_{i,j-\mathrm{1}}^n}{{\left(\mathrm{\Delta }{\mathit{\zeta }}_{i,j+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\zeta }}_{i,j-\mathrm{1}/\mathrm{2}}\right)}^{\mathrm{2}}}}+\\ & -\mathit{\rho }c{\displaystyle \frac{u_i^n}{L^n}}\left[{\displaystyle \frac{{\mathit{\theta }}_{i+\mathrm{1},j}^n-{\mathit{\theta }}_{i-\mathrm{1},j}^n}{\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i-\mathrm{1}/\mathrm{2}}}}\right.\\ & -\left({\displaystyle \frac{b_{i+\mathrm{1}}-b_{i-\mathrm{1}}}{\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i-\mathrm{1}/\mathrm{2}}}}\right.\\ & \left.+{\mathit{\zeta }}_{i,j}{\displaystyle \frac{H_{i+\mathrm{1}}^n-H_{i-\mathrm{1}}^n}{\mathrm{\Delta }{\mathit{\sigma }}_{i+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\sigma }}_{i-\mathrm{1}/\mathrm{2}}}}\right)\\ & \left.\times {\displaystyle \frac{{\mathit{\theta }}_{i,j+\mathrm{1}}^n-{\mathit{\theta }}_{i,j-\mathrm{1}}^n}{H_i^n\left(\mathrm{\Delta }{\mathit{\zeta }}_{i,j+\mathrm{1}/\mathrm{2}}+\mathrm{\Delta }{\mathit{\zeta }}_{i,j-\mathrm{1}/\mathrm{2}}\right)}}\right]+{\mathrm{\Phi }}_i^n.\end{array}\end{array}$$

A6 Adaptive time stepping

We take an adaptive time-stepping approach to enhance the computational performance of the flowline model. Unlike the proportional–integral (PI) methods, we exploit the fact that Picard's iteration already computes a metric to determine convergence. Thus, without additional calculations, we let the new time step evolve within a range set by the user [Δt_min,Δt_max] with a quadratic dependency on the error:

$$\begin{array}{}\text{(A22)}& \begin{array}{rl}\mathrm{\Delta }\stackrel{\mathrm{̃}}{t}& =\left(\mathrm{1}-{\left({\displaystyle \frac{\text{min}\left[\mathit{\epsilon }\left(t\right),{\mathit{\varphi }}_{\text{pic}}\right]}{{\mathit{\varphi }}_{\text{pic}}}}\right)}^{\mathrm{2}}\right)\left(\mathrm{\Delta }{t}_{\text{max}}-\mathrm{\Delta }{t}_{\text{min}}\right)\\ & +\mathrm{\Delta }{t}_{\text{min}},\end{array}\end{array}$$

where ϕ_pic is the tolerance in Picard's iteration and ε(t) is the error in the current iteration defined as ${\mathit{\epsilon }}^i=\left|\right|u^i-u^{i-\mathrm{1}}\left|\right|$ (Smedt et al., 2010). If the solution has not converged in the given time step (i.e. ε>ϕ_pic), Eq. (A22) ensures that the time step is set to the minimum value.

Then, we apply some relaxation to provide stability and avoid spurious oscillations:

$$\begin{array}{}\text{(A23)}& \mathrm{\Delta }t=\mathit{\alpha }\mathrm{\Delta }t+(\mathrm{1}-\mathit{\alpha })\mathrm{\Delta }\stackrel{\mathrm{̃}}{t},\end{array}$$

where α=0.7. Finally, we ensure that the time step is no greater than the Courant–Friedrichs–Lewy (CFL) condition:

$$\begin{array}{}\text{(A24)}& \mathrm{\Delta }t=\text{min}\left[\mathrm{\Delta }t,\mathrm{\Delta }{t}_{\text{CFL}}\right].\end{array}$$

Appendix B: Convergence

Nix's fast computation capabilities allow us to reach extremely high resolutions of Δx≤0.1 km. We describe herein the behaviour of the main variables of glaciological interest upon reducing the grid size. Given that Nix tries to simulate marine-terminating ice sheets, all variables are evaluated at the grounding line. For simplicity, the convergence experiments are performed over a linear prograde slope as described in Experiment 1 of MISMIP (Pattyn et al., 2012).

Figure B1 shows that ice velocities and thickness largely change as the resolution is increased. The former gradually decreases its value, converging to the analytical result of ∼ 700 m yr⁻¹ for Δx≤0.5 km. The ice thickness at the grounding line exhibits an opposite dependency with resolution, also converging to the analytical counterpart for Δx≤0.5 km. The total outflow of ice at the grounding line q(L) is not a monotonic function of the grid size. Since the grounding-line position is dictated by the ice flux value, the Nix model first advances with increasing resolution and ultimately retreats. It is worth noting that, despite the fact that ice velocity, thickness and fluxes are significantly close to the analytical solution, Nix overestimates the grounding-line position by roughly 6 % for the highest-resolution experiment Δx=0.06 km.

Figure B1 Convergence study with the Nix model. From top to bottom: (a) ice velocity, (b) terminus position, (c) ice thickness and (d) ice flux. All variables are evaluated at the grounding line. The dashed grey lines illustrate the values corresponding to the analytical counterpart (Schoof, 2007a). The double x axis denotes the total number of horizontal grid points n and the corresponding spatial resolution Δx given the stretched coordinate transformation.

In summary, Nix's velocity and ice thickness solutions converge to the analytical results elaborated on by Schoof (2007a). As is expected for any numerical model, a certain bias is found for high resolutions (≤0.13 km). Nonetheless, such a description would require solving a higher-order stress balance instead, for which the absence of an analytical solution complicates the comparison. This exercise shows that the grounding-line position for coarser grids might be close to that of the highest resolution for the wrong reason: an overestimation of the ice velocities combined with an underestimation of the ice thickness.