IMBAY – a multi-approximation 3 D ice-dynamics model for comprehensive applications : model description and examples

Glaciers and ice caps exhibit currently the largest cryospheric contributions to sea level rise. Modelling the dynamics and mass balance of the major ice sheets is therefore an important issue to investigate the current state and the future response of the cryosphere in response to changing environmental conditions, namely global warming. This requires a powerful, easy-to-use, versatile multi-approximation ice dynamics model. Based on the well-known and established ice sheet model of Pattyn (2003) we develop the modular multi-approximation thermomechanic ice model R IMBAY , in which we improve the original version in several aspects like a shallow ice–shallow shelf coupler and a full 3D-groundingline migration scheme based on Schoof’s (2007) heuristic analytical approach. We summarise the full Stokes equations and several approximations implemented within this model and we describe the different numerical discretisations. The results are cross-validated against previous publications dealing with ice modelling, and some additional artificial set-ups demonstrate the robustness of the different solvers and their internal coupling. RIMBAY is designed for an easy adaption to new scientific issues. Hence, we demonstrate in very different set-ups the applicability and functionality of R IMBAY in Earth system science in general and ice modelling in particular.


Introduction
According to the Fourth Assessment Report (AR4) of the Intergovernmental Panel on Climate Change (IPCC) (IPCC, 2007) it is unequivocal, that Earth's climate is warming since about 1850.This trend has been observed e.g. in rising air and ocean temperatures, in increased snow and ice melting, and in a rising sea level.According to more recent publications (e.g.Church et al., 2011;Rahmstorf et al., 2012) the trends estimated even for the worst scenarios of the AR4 are already reached or surpassed.Therefore, the imminent climate change will have profound impact on society (e.g.Hanson et al., 2011).
However, none of the complex numerical Earth system models (ESMs) in the IPCC report, used to compute the future climate trends, include the possible climate feedbacks of the large ice sheets in Greenland and Antarctica, resulting in large uncertainties for the global mean sea-level predictions.
These ice sheets play a crucial role in the Earth's hydrological cycle as they store about 75 % of the Earth's fresh water.In general, ice sheets accumulate mass from snow precipitation, which is compacted and finally transformed into ice.It follows the gravitational force and flows downhill from summits towards the ice sheet margin.However, this simplified view gets much more complex as different flow regimes exist within ice sheets (Fig. 1): the ice sheet's homogeneity is disturbed by nunataks and fast flowing ice streams; at the base, subglacial lakes and a hydrological network alternates the basal boundary conditions of the ice sheet; and at the edges ice shelves interact with the ocean by massive melting and Published by Copernicus Publications on behalf of the European Geosciences Union.
Fig. 1.Sketch illustrating several aspects/components to be considered in ice sheet modelling (adapted after Sandhäger, 2000).oped since 2009.Although the underlying Higher Order 70 Model (HOM) and Full Stokes (FS)-physics remained basically unchanged, a Shallow Shelf Approximation (SSA) solver has been added to calculate the horizontally averaged velocities of ice streams and ice shelves.Additionally, the numerical solver implementation, the discretisation, the cou-75 pling between different solvers and the user interface have been improved in many aspects since it diverted from the original model.Keeping in mind that ice models have to deal with many different geophysical settings and boundary conditions, it is challenging to design a computer code which 80 is able to fulfil these needs for a large variety of users and applications.RIMBAY has been designed to be easy applicable to new scenarios, easy to extend and with clear interfaces to couple it with existing codes.
This paper is structured as follows: First, we clarify in 85 section 2 the sometimes imprecise usage of the term model, before we present in section 3 the mathematical equations and several approximations founding the mathematical background of RIMBAY.Thereafter, we describe the numerical finite-difference implementation of these equations and how 90 they can be solved with existing numerical solvers for linear differential equations in section 4. Some more details about the code-implementation are given in section 5, before we present some idealised example-applications of RIMBAY, with a main focus on cross-validation with previously pub-95 lished ice-model results and an example of internal codecoupling in section 6.Finally, we demonstrate in section 7 the wide spectrum of applications RIMBAY is already used for by several users.

Multi-approximation ice sheet/shelf model RIMBAY 100
The term model is used in several ways in Earth system science, which can be sometimes confusing.Therfore, we first define what we understand as model, or to be more precise between which types of model we distinguish: -Equations form the mathematical model describing the 105 fundamental relationship between the relevant values of interest (e.g., velocity, temperature, and viscosity).In our context, these equations are mostly coupled differential equations which can not be solved analytically.-These equations are solved with a computer, which re-110 quires a discretisation of the equations.This can be done in several distinct ways, depending on the demand of accuracy, stability, convergence properties, and resources (memory usage and computational coast).We refer to this as the numerical model.

115
-This numerical model has to be translated into a computer language (mostly a high-level programming language like Matlab, Fortran, C, or C++).It is common sense to refer to this computer program as model, too.
We use the expression code or the implementation to 120 specify the lines forming this (sometimes compiled binary) program.-Finally, the code is applied to answer a specific scientific question (e.g. the contribution to sea level rise) of a specific domain (e.g., whole Antarctica or a subregion 125 like the area of the Pine Island Glacier), or to study processes (e.g. the impact of basal water on ice dynamics) and the sensitivity to parameters or boundary conditions (e.g.geothermal heat flux, bedrock topography or ice thickness distribution).These applications of a com-130 puter program are often called model, too.We refer to these applications as experiments or scenarios.
In general, we use the term RIMBAY for the implementation of the discretised equation, and therefore the compiled binary code, which includes not only the mathematical model, but also a sophisticated command-line interpreter and inputoutput interfaces for an easy usage.RIMBAY is distributed with a suit of example-and reference-scenarios and several additional programs (mainly based on the bash-script language) providing several options to visualize the computed 140 results with the Generic Mapping Tools (GMT, Wessel and Smith, 1998;Wessel et al., 2013).In the following sections, we elaborate on these different model types and how they are used in RIMBAY.

Mathematical model 145
The mathematical field equations are based upon the conservation of mass, momentum, and energy with the (constant) density ρ, the velocity vector v = (v x , v y , v z ) = (u, v, w), the gravitational acceleration g = Fig. 1.Sketch illustrating several aspects/components to be considered in ice sheet modelling (adapted after Sandhäger, 2000).
iceberg calving.Therefore, a numerical model has to deal with many different aspects of an ice sheet (and ice shelf) to represent it's complex dynamic behaviour adequately and to improve future projections or hindcasts for palaeoclimatology.
During the last years great efforts have been undertaken to improve existing ice models and to incorporate them into coupled climate models (e.g.Rutt et al., 2009;Gillet-Chaulet et al., 2012;Levermann et al., 2012).Here, we present the Revised Ice Model Based on frAnk pattYn, the multi-approximation ice sheet/ice shelf model RIMBAY.This model is originally based on the higher-order numerical ice-flow model of Pattyn (2003), which has been tested and applied to many scenarios (e.g.Pattyn, 2002;Pattyn et al., 2004;Pattyn, 2008Pattyn, , 2010)).RIMBAY itself has been developed since 2009.Although the underlying higher order model (HOM) and full Stokes (FS)-physics remained basically unchanged, a shallow shelf approximation (SSA) solver has been added to calculate the horizontally averaged velocities of ice streams and ice shelves.Additionally, the numerical solver implementation, the discretisation, the coupling between different solvers and the user interface have been improved in many aspects since it diverted from the original model.Keeping in mind that ice models have to deal with many different geophysical settings and boundary conditions, it is challenging to design a computer code which is able to fulfil these needs for a large variety of users and applications.RIMBAY has been designed to be easy applicable to new scenarios, easy to extend and with clear interfaces to couple it with existing codes.
This paper is structured as follows: first, we clarify in Sect. 2 the sometimes imprecise usage of the term model, before we present in Sect. 3 the mathematical equations and several approximations founding the mathematical background of RIMBAY.Thereafter, we describe the numerical finite-difference implementation of these equations and how they can be solved with existing numerical solvers for linear differential equations in Sect. 4. Some more details about the code implementation are given in Sect.5, before we present some idealised example applications of RIMBAY, with a main focus on cross-validation with previously published ice-model results and an example of internal code coupling in Sect.6.Finally, we demonstrate in Sect.7 the wide spectrum of applications RIMBAY is already used for by several users.

Multi-approximation ice sheet/shelf model RIMBAY
The term model is used in several ways in Earth system science, which can be sometimes confusing.Therefore, we first define what we understand as model and, to be more precise, between which types of models we distinguish.
-Equations form the mathematical model describing the fundamental relationship between the relevant values of interest (e.g.velocity, temperature, and viscosity).
In our context, these equations are mostly coupled differential equations which can not be solved analytically.
-These equations are solved with a computer, which requires a discretisation of the equations.This can be done in several distinct ways, depending on the demand of accuracy, stability, convergence properties, and resources (memory usage and computational coast).We refer to this as the numerical model.
-This numerical model has to be translated into a computer language (mostly a high-level programming language like Matlab, Fortran, C, or C++).It is common sense to refer to this computer program as model, too.We use the expression code or the implementation to specify the lines forming this (sometimes compiled binary) program.
-Finally, the code is applied to answer a specific scientific question (e.g. the contribution to sea level rise) of a specific domain (e.g.whole Antarctica or a subregion like the area of the Pine Island Glacier), or to study processes (e.g. the impact of basal water on ice dynamics) and the sensitivity to parameters or boundary conditions (e.g.geothermal heat flux, bedrock topography or ice thickness distribution).These applications of a computer program are often called model, too.We refer to these applications as experiments or scenarios.
In general, we use the term RIMBAY for the implementation of the discretised equation, and therefore the compiled binary code, which includes not only the mathematical model, but also a sophisticated command-line interpreter and inputoutput interfaces for an easy usage.RIMBAY is distributed with a suit of example -and reference -scenarios and several Geosci.Model Dev., 7, 1-21, 2014 www.geosci-model-dev.net/7/1/2014/additional programs (mainly based on the bash-script language) providing several options to visualise the computed results with the generic mapping tools (GMT) (GMT, Wessel and Smith, 1998;Wessel et al., 2013).In the following sections, we elaborate on these different model types and how they are used in RIMBAY.

Mathematical model
The mathematical field equations are based upon the conservation of mass, momentum, and energy: with the (constant) density ρ, the velocity vector v = (v x , v y , v z ) = (u, v, w), the gravitational acceleration g = (0, 0, −g), the two dimensional stress tensor τ , the (potential) temperature θ, the heat capacity c, the thermal conductivity κ, and the internal frictional heating Q i .In the following we consider Cartesian coordinates, with the vertical coordinate z upwards and neglect acceleration.In case of an incompressible fluid with a constant density the continuity equation (conservation of mass) follows as The stress tensor τ is split into a deviatoric part τ and an isotropic pressure, which is defined as the negative trace of the stress tensor: where 1 symbolises the identity matrix.

Equation of motion
Because velocities in ice sheet/shelf modelling are rather small, acceleration can be ignored and the momentum equation can be written as According to Paterson (1994), the constitutive equation for polycrystalline ice links the deviatoric stresses to the strain rates, applying the effective viscosity η, which can be described by the Glen-type flow law (e.g.Cuffey and Paterson, 2010): with n = 3, the pressure-corrected ice temperature θ = θ + αp, with a constant α = 9.8 × 10 −4 K Pa −1 (Greve and Blatter, 2009), and the effective strain rate (valid for incompressibility as ˙ xx + ˙ yy + ˙ zz = 0 follows from Eq. 4) The temperature dependent rate factor A( θ ) is parameterised according to the Arrhenius relationship after Hooke (1981) or Paterson and Budd (1982).Combining Eq. ( 6) and Eq. ( 7) we get the so-called full Stokes (FS) equations for ice modelling: surface S to the height z (Van der Veen and Whillans, 1989;Pattyn, 2008): Here, the first term in Eq. ( 12) describes the hydrostatic part and R zz the resistive part, sometimes also referred to as vertical resistive longitudinal stress.
Depending on the scientific issue, several approximations of Eq. ( 10) might be reasonable, which are described in the following subsection.

Higher-order approximation
The HOM approximation of Pattyn (2003) applies the hydrostatic approximation, by neglecting the resistive stress R zz in Eqs. ( 10)-( 12) for the vertical velocity and the vertical normal stress.These are only relevant (but still almost two orders of magnitude below the other normal stress and shear stress components, Pattyn, 2000) where the ice flow regime changes, as in the vicinity of ice margins or ice divides.Additionally, ignoring the horizontal derivatives of the vertical velocity in Eq. ( 10 for the horizontal velocities.The vertical velocity at depth z can be derived by integrating the continuity equation Eq. ( 4) from the base B vertically:

Shallow shelf or shelfy stream approximation
A second common approximation is the shallow shelf approximation or shelfy stream approximation (SSA).This assumes that the horizontal velocity is depth-independent ( ∂u ∂z = ∂v ∂z = 0), which is the case for ice shelf regions and fast flowing ice streams decoupled from the ground.Integrating Eq. ( 14) through the ice from the base B to the surface S, and defining U and V as the vertically integrated velocities leads to (e.g.Morland, 1987;MacAyeal, 1989;Pattyn, 2010 where the basal shear stress τ bi retards the otherwise unhampered flow on bedrock till.It can be expressed in terms of the basal friction parameter β 2 and the horizontal velocity: A thorough derivation of Eq.16 can be found in Greve and Blatter (2009).Both, the shelfy stream approximation and the shallow shelf approximation are expressed by Eq. ( 16).The only difference is, that for an ice shelf or above a subglacial lake β 2 is zero, while it might reach several thousand Pa a m −1 for a slippery bedrock, which especially applies to basal lubricated areas.As a rule of thumb above dry bedrock a value of β 2 = 25 000 Pa a m −1 would correspond to a typical frictional stress of about 100 kPa (Paterson, 1994) if a velocity of about 4 m a −1 is assumed (Thoma et al., 2012).Finally, because of the lacking vertical shear stresses Eq. ( 9) reduces to

Shallow ice approximation
The most rigid approximation is the shallow ice approximation (SIA, which is a reasonable simplification for large ice bodies, when the horizontal length scale is much larger than the ice thickness (e.g.Hutter, 1983).Assuming that the horizontal derivation of the vertical velocity is much smaller than the vertical derivation of the horizontal velocity ( ∂w ∂x ∂u ∂z ) and applying the hydrostatic approximation (which reduces the vertical momentum balance to the hydrostatic term) we derive Basically, this approximation decouples the horizontal velocities, allowing local solutions for the velocity field, instead of a much more complex and time-consuming implicit solver.
The numerical resources of this SIA are so low (compared to any other approximations) that it is still widely used (and useful) for many applications.

Boundary conditions
Several boundary conditions have to be formulated to solve the different approximations of the equation of motion.
1. We apply a stress-free surface boundary condition: with the normal vector n s orthogonal to the surface.
2. For the horizontal velocities at the ice base, we apply either -a no-slip condition for the tangential velocities -a Weertman-type sliding law (e.g.Paterson, 1994;Cuffey and Paterson, 2010), linking the sliding velocity with the basal shear stress: The basal drag is defined as the sum of all basal resistive forces (Van der Veen and Whillans, 1989;Pattyn, 2003): with τ ij = τ ij (B).In case of the SIA these equations simplify to -or a stress free base when a substantial amount of water is present, like in the case of subglacial lakes and ice shelves; this implies β 2 = 0.
3. For the vertical velocity at the base, we apply a kinematic boundary condition: with the basal melt rate ṁB .
4. At lateral boundaries of the model domain, we apply either -zero ice thickness (H = 0), -Dirichlet boundary conditions with fixed velocities.The no-slip condition (u = 0), which would imply frozen ice at nunataks, is a special case of this.
-A Neumann free-slip boundary condition: at ice-nunatak edges, with the unit vector n ⊥ orthogonal to the edge, or -a (dynamic) Neumann boundary condition for an ice shelf-ocean interface (e.g.Greve and Blatter, 2009;Joughin et al., 2009;Pattyn, 2010) with the outward-pointing unit vector (n x , n y ), which is perpendicular to the (vertical) ice shelf front.
-or periodic boundary conditions.
These equations are converted in terrain following σ coordinates by applying with the ice thickness H and the surface height S.This coordinate transformation leads to additional metric terms in the equations, which are described in detail in Pattyn (2003) or Greve and Blatter (2009).The advantage is, that the vertical coordinate ranges from σ = 0 at the surface to σ = 1 at the ice base, independent of the local ice thickness and the bedrock elevation.

Temperature calculation
Assuming a constant thermal conductivity κ, the temperature evolution (Eq. 3) can be divided into an advective, a diffusive and a source term: www.geosci-model-dev.net/7/1/2014/Geosci.Model Dev., 7, 1-21, 2014 Neglecting horizontal diffusion and assuming that the internal heat source results mainly from the ice deformation (Paterson, 1994) we obtain with the effective deviatoric stress τ (defined similar to the effective strain rate in Eq. 9) and The boundary conditions applied to solve this thermodynamic equation are -the mean air temperature at the surface of the ice body, -a Dirichlet boundary condition according to the pressure melting point of ice (e.g.Paterson, 1994), θ = −8.7 × 10 −4 K m −1 H , at the ice base when the ice is floating (like above subglacial lakes and for ice shelves), and -a Neumann boundary condition at the base b for grounded ice: with the basal stress τ b = τ 2 bx + τ 2 by and the geothermal heat flux G.

Ice sheet evolution
Integration of Eq. ( 1) from the base B to the surface S leads to an equation for the ice evolution.Defining the ice thickness H = S − B, accounting for melting or accumulation at the surface and/or base and assuming a constant ice density ρ we get with the mass balance (in m a −1 ) is defined as Basal freezing can be implemented by negative basal melting.
4 Numerical model

Linear and non-linear solvers
The coupled pair of equations for the horizontal velocity field for the FS, HOM, and SSA equations (Eqs. 10,14,16) depend on the strain-rate dependent viscosity (Eq.8), resulting in a non-linear problem.In the full Stokes case, the horizontal velocities depend also on the vertical velocity.However, these problems can be solved iteratively as indicated in Fig. 2.
In the full Stokes case, the vertical velocity w can be estimated from the continuity Eq. ( 4), imposing kinematic boundary condition at the lower ice surface (including melt rates).
According to Pattyn (2003), it is sufficient to solve the system of linear equations for u and v successively, instead of solving both equations at once.In general, we iteratively solve where l is the iteration, A nm contains the coefficients of the left-hand side of the relevant equation to solve, while b n is the forcing term on the right-hand side of the equation.The placeholder x m symbolises the horizontal velocities u ij or v ij (Eqs.10, 14, or 16), the potential temperature θ ij (Eq.28), or the ice thickness H ij (Eq.31), respectively.The indices n and m symbolise the consecutively numbered grid nodes from Two methods to solve the linear system of Eq. ( 33) are available within RIMBAY: first, the fast and efficient biconjugate gradient method with a Jacobian preconditioner (linbcg) from the numerical recipes (NR) (Press et al., 2007); second, the Library of Iterative Solvers (LIS) from Nishida (2010).
The Library of Iterative Solvers (LIS) provides a bunch of preconditioners and solvers, including the recommendable generalised minimal residual (gmres) method, which can also be applied to solve non-symmetric matrices.For both methods the effective compressed row storage (CRS) sparse matrix method is used as a default to store the elements of the matrix A nm .However, for the LIS, the modified sparse row (MSR) format is implemented, too.Comparisons with respect to the calculated velocities have shown -the differences for the two storage formats (CRS vs. MSR) are negligible, -the differences between the linbcg solver from Press et al. ( 2007) and the very same preconditioner/solver combination from the Library of Iterative Solvers are negligible, but the solver of Press et al. ( 2007) needs less computational resources.
-if specific preconditioner-solver combinations converge, the difference between different combinations are negligible.
Summarised, if a solution of the linear system can be computed with a reasonable accuracy, the results can be trusted.6 M. Thoma et al.: Description of the ice flow model RIMBAY Fig. 2. Sequence of iteratively solved variables within RIMBAY.In the SSA-case the product of Hη is calculated, instead of the viscosity η, only.The grayish highlighted variables are calculated only in the FS-case.(Within the main loop, the vertical velocity needs only to be calculated in the non FS-cases.)get with the mass balance (in m/a) is defined as Basal freezing can be implemented by negative basal melt-
4 Numerical model

Linear and non-linear solvers
The coupled pair of equations for the horizontal velocity field for the FS, HOM, and SSA equations (Eq. 10, Eq. 14, and Eq.16) depend on the strain-rate dependent viscosity (Eq.8), resulting in a non-linear problem.In the full Stokes case, the horizontal velocities depend also on the vertical velocity.However, these problems can be solved iteratively as indicated by Fig. 2. In the full Stokes case, the vertical velocity w can be estimated from the continuity equation 4, imposing kinematic boundary condition at the lower ice surface (including melt rates).
According to Pattyn (2003), it is sufficient to solve the system of linear equations for u and v successively, instead of solving both equations at once.In general, we iteratively solve where l is the iteration, A nm contains the coefficients of the left-hand side of the relevant equation to solve, while b n is the forcing term on the right-hand side of the equation.The placeholder x m symbolizes the horizontal velocities u ij or v ij (Eq. 10, Eq. 14, or Eq.16), the potential temperature θ ij (Eq.28), or the ice thickness H ij (Eq.31), respectively.The indices n and m symbolize the consecutively numbered grid Two methods to solve the linear system of Eq. 33 are available within RIMBAY: First, the fast and efficient biconjugate gradient method with a Jacobian preconditioner (linbcg) from the Numerical Recipes (NR) (Press et al., 2007), second 230 the Library of Iterative Solvers (LIS) from Nishida (2010).The Library of Iterative Solvers (LIS) provides a bunch of preconditioners and solvers, including the recommendable Generalized Minimal RESidual (gmres) method, which can also be applied to solve non-symmetric matrices.For both 235 methods the effective compressed row storage (CRS) sparse matrix method is used as a default to store the elements of the matrix A nm .However, for the LIS, the modified sparse row (MSR) format is implemented, too.Comparisons with respect to the calculated velocities have shown Summarised, if a solution of the linear system can be computed with a reasonable accuracy, the results can be trusted.
In general, we suggest to start with the faster linbcg algorithm (from Press et al., 2007) and to switch to the gmres solver with a Jacobian or ILU (Incomplete LU decomposi-255 tion) preconditioner if it should fail.When solving the linearized equations for the horizontal velocity, the viscosity η might vary over a few orders of magnitude, which requires a sophisticated convergence scheme.Hence, a simple Picard iteration might fail.There-260 fore, Pattyn (2003) extended this scheme by the unstable manifold correction (UMC), introduced by Hindmarsh and Payne (1996), which results in a proper convergence of the solution.In RIMBAY the UMC is applied in the SSA, HOM, and FS solvers.In general, we suggest to start with the faster linbcg algorithm (from Press et al., 2007) and to switch to the gmres solver with a Jacobian or ILU (Incomplete LU decomposition) preconditioner if it should fail.
When solving the linearised equations for the horizontal velocity, the viscosity η might vary over a few orders of magnitude, which requires a sophisticated convergence scheme.Hence, a simple Picard iteration might fail.Therefore, Pattyn (2003) extended this scheme by the unstable manifold correction (UMC), introduced by Hindmarsh and Payne (1996), which results in a proper convergence of the solution.In RIMBAY the UMC is applied in the SSA, HOM, and FS solvers.

Discretisation
When equations are discretised, it is important to realise where exactly the individual variables are located.This is quite simply defined for the unstaggered Arakawa A-grid (e.g.Arakawa and Lamb, 1977;Purser and Leslie, 1988) where all variables are located in the very same grid position.However, sometimes a different approach has numerical advantages.Besides the traditional (unstaggered) A-grid, the staggered Arakawa C-grid is optionally available in RIM-BAY for the SIA and SSA solvers.On the Arakawa C-grid, the horizontal velocities are defined in between the thickness (and viscosity) nodes as illustrated in Fig. 3.

Ice Sheet evolution
As an example, we formulate the implemented discretisation of the ice sheet evolution (Eq.31) explicitly for the two different grids.Additionally, the detailed discretisation on the C-grid of the SSA equation of motion (Eq.16) is given in Appendix B. quite simply defined for the unstaggered Arakawa A-Grid (e.g., Arakawa and Lamb, 1977;Purser and Leslie, 1988) 270 where all variables are located in the very same grid position.However, sometimes a different approach has numerical advantages.Besides the traditional (unstaggered) A-Grid, the staggered Arakawa C-Grid is optionally available in RIM-BAY for the SIA-and SSA-solvers.On the Arakawa C-Grid,

275
the horizontal velocities are defined in between the thickness (and viscosity) nodes as illustrated in Fig. 3.

Ice Sheet evolution
As an example, we formulate the implemented discretisation of the ice sheet evolution (Eq.31) explicitly for the two dif-280 ferent grids.Additionally, the detailed discretisation on the C-grid of the SSA equation of motion (Eq.16) is given in appendix B.

C-Grid
For the C-Grid, where the velocities are defined inbetween thickness nodes, the equation of the ice sheet evolution (Eq.31) can be written as an implicit first order finite difference equation as Rearrang values, l in the fol linear so The co the inter mulated grid cell plicitly i grid cell.
If the domain, these by edge of boundari 0) edge ∆t 2 Vi,j ∆y

C-grid
For the C-grid, where the velocities are defined in-between thickness nodes, the equation of the ice sheet evolution (Eq.31) can be written as an implicit first order finite difference equation as Rearranging Eq. ( 34) with respect to the five discrete H t+1 values, located and numbered as indicated in Fig. 4, results in the following coefficients for the sparse matrix A nm of the linear solver: Rearranging Eq. 34 with respect to the five discrete H t+1values, located and numbered as indicated by Fig. 4, results in the following coefficients for the sparse matrix A nm of the linear solver: These coefficients represent the non-zero elements of each single row n for each ij-element of the matrix A nm , with C n3 indicating the central node at (i,j) and b n indicates the forcing term on the right-hand-side.The coefficients derived in the last subsection are valid for the interior of the ice.Boundary conditions have to be formulated at the edges of the ice sheet.Open boundaries for grid cells adjacent to ocean or ice-free land are simply implicitly implemented by assuming H = 0 at the respective grid cell.
If the ice adjoins a nunatak or a lateral end of the model domain, closed boundary conditions are applied.We define These coefficients represent the non-zero elements of each single row n for each i, j element of the matrix A nm , with C n3 indicating the central node at (i,j ) and b n the forcing term on the right-hand side.
The coefficients derived in the last subsection are valid for the interior of the ice.Boundary conditions have to be formulated at the edges of the ice sheet.Open boundaries for grid cells adjacent to ocean or ice-free land are simply implicitly implemented by assuming H = 0 at the respective grid cell.
If the ice adjoins a nunatak or a lateral end of the model domain, closed boundary conditions are applied.We define these by setting the velocity (and thus the flux) of ice over the edge of the specific grid cell, to zero.For example, closed boundaries at the eastern (U i,j = 0) and southern (V i−1,j = 0) edges would result in C n1 = C n4 = 0 and in Eq. ( 35).

A-grid
For the the A-grid a pure advective scheme to solve Eq. ( 31) would be numerically problematic, because of the velocitypressure gradient coupling.Hence, we decompose the equation into a weighted advective and diffusive part by applying the identity (∇H +∇B)(∇S) −1 = 1, derived from a simple gradient formulation of S = H + B (surface elevation S equals ice thickness H plus ice bottom B): with f ad = 1 for pure diffusion and f ad = 0 in case of pure advection.With the definition of the non-linear (because it depends on the solution H ) diffusion vector The finite difference formulation of Eq. ( 37) as well as the coefficients C nm for the sparse linear matrix for the interior and the boundary conditions are given in the Appendix A.
In general, it would be appropriate to apply the diffusive equation (with f ad = 1), because a Lax method has to be used to numerically stabilise the advective part of Eq. ( 37) (see Appendix A).Unfortunately, this adds numerical dissipation (numerical diffusion) and results in a time-step dependence of the solution.However, if the ice body contains ice shelves and/or ice divides with flat areas, the reciprocal value of (∇S) −1 becomes very large and counteracts the stabilising effect of the otherwise stable diffusive implementation.Despite this problem of exchanging stability towards convergence (with respect to decreasing time steps) this approach has been discussed in some applications (e.g.Pattyn et al., 2006;Docquier et al., 2011).
As an alternative to overcome the restrictions involved with the numerical representation with respect to (∇S) −1 in Eq. ( 36), we implemented a mass conserving (time step independent) upwind scheme, based on Eq. ( 34).Averaging the horizontal velocities from their central (A-grid) location towards the grid-cell edges according to U c i,j = 1 2 (U i,j + U i,j +1 ) and V c i,j = 1 2 (V i,j + V i+1,j ) leads to and the following coefficients for the sparse matrix A nm : Geosci.Model Dev., 7, 1-21, 2014 www.geosci-model-dev.net/7/1/2014/ b n =H t i,j + ṁi,j t. (39)

Ice sheet-ice shelf coupling and grounding line flux
The mechanically correct way of coupling a ice sheet system with an ice shelf system would be a FS approach.According to Pattyn et al. (2013) a horizontal resolution of less than 0.5 km is necessary to capture the grounding line (GRL) migration accurately.This however, is computationally costly and inefficient, especially for major parts of the ice sheet and ice shelf, which are at a large distance from the GRL where reduced physics is sufficient (see Eqs. 16 and 18, respectively).Either a finite element discretisation (as in the Elmer/ice model or the Ice Sheet System Model (ISSM), e.g.Zwinger et al., 2007;Larour et al., 2012) or a finite volume approach are necessary to implement FS physics in a reasonable way.Another approach to increase the grid resolution in specific regions of interest are adaptive grids (e.g.Gladstone et al., 2010;Cornford et al., 2013), although they have (to our knowledge) not been applied to HOM or FS physics, yet.For coarse resolution finite difference models (with grid sizes beyond one kilometre), Pollard and DeConto (2009) and Pollard and DeConto (2012) suggested a heuristic approach, based on the semi-analytical grounding line flux solution, derived by Schoof (2007): with the longitudinal stress τ xx just downstream of the grounding line, and the unbuttressed stress τ f = 0.5 ρgh g (1 − ρ/ρ Ocean ).The grounding line flux Q S x is estimated from the ice thickness h g at the interpolated sub-grid grounding line position.Before the ice evolution Eq. ( 30) is solved, the Schoof flux (estimated on a sub-grid scale) constrains the flux across the grounding line by correcting the previous estimated velocity, located on a discrete grid node according to With the ice thickness H at the last grounded node (the model's grounding line) and H float the ice thickness at the first floating node downstream.The distinction depends on the relation between the analytical Schoof-flux Q S i and the modelled flux through the last gridded node i then more ice is transported into the ice shelf and the grounding line retreats or stays constant, the velocity at the grounding line is corrected according to Eq. ( 41).If Q S i < Q M i then less ice is transported into the ice shelf and more ice is kept in the ice sheet, the grounding line advances and the velocity of the first floating node is corrected according to Eq. ( 41).A detailed description of this method is given in Docquier et al. (2011), Pollard and DeConto (2012), and Pattyn et al. (2012).To avoid unrealistic velocity steps, we additionally apply a conservative 2D-Gaussian filter to the grounding line nodes to smooth the resulting velocity field.

General Information
The RIMBAY code is mainly2 written in C++ and has about 30 000 (mostly) well documented lines.For historical reasons the code is not completely object oriented yet, the majority is organised into classes, and the number of global variables (which should be avoided as much as possible in any code) is close to zero.A reasonable degree of code separation into several C++ classes, allows an easy maintenance of the code.Well-defined interfaces (public methods of the C++ classes) enable an easy extension of the code for upcoming developments in ice modelling and/or further reaching applications (see Sect. 7).
The GNU build system3 (also known as the autotools) is a suite of programming tools designed to assist in making source-code packages portable to many Unix-like systems.It generates system-and environment-dependent Makefiles automatically and attends dependencies between different source (and header) files.Thanks to the GNU build system, RIMBAY has been compiled and tested successfully on several different Unix-platforms without any code adjustments.To distribute, develop, and maintain RIMBAY we use the distributed revision control system monotone4 , which keeps track of any changes within the code and provides a sophisticated automatic merging of development branches.
One of the main programming paradigms for RIMBAY is that the very same (compiled) code has to run every single (previous successfully tested) scenario without any code editing and/or recompiling.To achieve this, RIMBAY is started with command-line arguments and loads the specific scenario from parameter files and (if requested) optionally from a netcdf file, too.The well established netcdf output format of RIMBAY ensures that the computed results can subsequently be post-processed with the desired software packages, if the supplied GMT-bash scripts (Wessel and Smith, 1998;Wessel et al., 2013) (included in the RIMBAY-monotone database) should not be sufficient.
The RIMBAY code comes with a test suite containing nearly 50 different scenarios.These small and fast-running scenarios are designed to ensure that future model developments do not interfere with previous results.

Solver coupling
The coupling of SIA and SSA at the grounding line is realised by applying depth averaged velocities from the SIA solver as a Dirichlet boundary condition for the SSA solver.This transition can be located either at the last grounded node (the numerical GRL) or several grid nodes inside the ice sheet.In the latter case a transition zone (or grounding zone) is defined by a region where the solutions of the SIA and the SSA solvers are interpolated.
If the HOM/FS should not be applied to the whole model domain (which might be reasonable to save computational time), one or more region(s) of interest can be defined.In that case, the resource-consuming HOM/FS solver is limited to these regions only, while the faster SIA and SSA solvers are applied elsewhere and provide the Dirichlet boundary conditions for the HOM/FS solver (see example in Sect.6.4).

Validation
The implementation of the different mathematical models (SIA, SSA, HOM, and FS) to calculate the horizontal velocity field are validated separately in this subsection.Additionally we show that the solver for the ice sheet evolution and the solver coupling produces reasonable results.The temperature evolution and thermomechanical coupling is not reconsidered here.Although the solvers have been revised, their results are identical to those published by Pattyn (2003) and Thoma et al. (2012).

SIA solver
The A-grid implementation of the SIA within RIMBAY is mainly identical to those of Pattyn (2003) and has already been validated successfully against the moving-margin Eismint benchmark described in Huybrechts et al. (1996) within Pattyn (2003).Here, we compare the estimated ice thicknesses, derived with the A-grid (Type-II according to Huybrechts et al., 1996) and the C-grid RIMBAY implementation for the fixed-and moving-margin benchmark experiments, with results published by Huybrechts et al. (1996) and Bueler et al. (2005).Figure 5 shows that the A-grid implementation produces results very close to the reference, while the C-grid implementation results in a 0.38 % larger ice thickness.Considering the very different discretisations (compare

405
The A-Grid implementation of the SIA within RIMBAY is mainly identical to those of Pattyn (2003) and has already been validated successfully against the moving-margin Eismint benchmark described in Huybrechts et al. (1996) within Pattyn (2003).Here, we compare the estimated ice thick-410 nesses, derived with the A-Grid (Type-II according to Huybrechts et al., 1996) and the C-Grid RIMBAY implementaion for the fixed-and moving margin benchmark experiments, with results published by Huybrechts et al. (1996) and Bueler et al. (2005).Fig. 5 shows that the A-Grid implementa-415 tion produces results very close to the reference, while the C-Grid implementation results in a 0.38% larger ice thickness.Considering the very different discretisations (compare Eq.34 and Eq.A2) of the ice evolution equation, this is acceptable.420

SSA-Solver
The A-and the C-Grid implementations of the SSA are compared with a diagnostic tabular iceberg experiment of Jansen et al. (2005).In this experiment, the horizontal velocity field of a rectangular iceberg with a constant thickness of 250 m 425 and an isothermal temperature of −20 • C is calculated.The viscosity is calculated according to Eq.8 with n = 3 and a temperature dependent rate factor given by the Arrhenius relationship after Paterson and Budd (1982)  Eqs.34 and A2) of the ice evolution equation, this is acceptable.

SSA solver
The A-and the C-grid implementations of the SSA are compared with a diagnostic tabular iceberg experiment of Jansen et al. (2005).In this experiment, the horizontal velocity field of a rectangular iceberg with a constant thickness of 250 m and an isothermal temperature of −2 • C is calculated.The viscosity is calculated according to Eq. ( 8) with n = 3 and a temperature dependent rate factor given by the Arrhenius relationship after Paterson and Budd (1982).Our horizontal velocities are calculated on a 1 km grid and are in close agreement with those presented by Jansen et al. (2005) (Fig. 6a).Additionally, we rotate the iceberg, to demonstrate the independence of the model results from the iceberg's orientation within the rectangular grid (Fig. 6b-e).This test is essential for the modelling of evolving ice sheet fronts, which are rarely aligned with the grid orientation in real geometries.
A much more complex proof-of-concept is shown in Fig. 7.This artificially constructed geometry with a grid resolution of 2 km features -a non-constant ice thickness; -two discontinuous areas, which are solved simultaneously by the numerical solver; -a quite complex shaped ice-water front with corners, tongues, and an inlet at x ≈ 200 km.
-The brown areas in Fig. 7  modelling, possibly coupled with an atmosphere and ocean model in an Earth-System model approach.For reversibility tests and transient experiments the 5 km and 10 km resolutions were considered, only.In order to overcome the problem of capturing grounding line migration in coarse res-490 olutions, we apply the heuristic rule described in section 4.4.Despite the rather coarse resolution (compared with most of the other 15 participating numerical ice models), RIMBAY accomplished the velocity-field comparison of the diagnostic experiment and in particular the reversibility test for ground-495 ing line migration.The latter is a prerequisite for modelling the Antarctic Ice Sheet, which is sourrounded by ice shelves along more than half of its coast line.However, as a consequence of the imposed heuristic grounding line condition, which is not valid for a compressive flow, RIMBAY is (as all 500 other A-HySSA models participated in the MISMIP3d intercomparison are) incapable of reproducing a grounding line retreat at free-slip walls.Despite the partly discontinuous grounding line retreat, resulting from the rather coarse resolution, we conclude that RIMBAY is capable of simulating 505 large ice bodies with attached ice shelves for different climate conditions.-Additionally, a small nunatak (with an area of 10 km× 5 km = 50 km 2 ) located within the left iceberg with free-slip boundary conditions is added.
The modelled velocity pattern is consistent with the expectations, which are -higher velocities at higher ice fronts, -zero velocities at no-slip boundaries, and -a reduced, orthogonally orientated velocity field at free-slip boundaries.
The difference between the A-grid and the C-grid (not shown) are negligible.Therefore, we conclude that the SSAsolver implementation produces reasonable and robust results, even for complex geometries.

SIA-SSA solver coupling and GRL migration
The results of the Marine Ice Sheet Model Intercomparison Project (MISMIP) (Pattyn et al., 2012) are a good benchmark test for the capability of a coupled ice sheet/shelf model to predict grounding line migrations.In this 2D flow-line experiment the position of the GRL is compared with the boundary layer theory of Schoof (2007).We applied RIM-BAY with different horizontal resolutions and a transition zone of 100 km, imposing the heuristic condition according to the method described in Sect.4.4.Figure 8 indicates that the semi-analytical steady state grounding-line positions, according to the boundary layer theory of Schoof (2007), are in general reproduced well with RIMBAY.However, some delayed movements happen because of numerical issues in this idealised set-up.
Recently, RIMBAY participated in the extended 3D variant of the MISMIP, which investigates the grounding line response to external forcings (Pattyn et al., 2013).We performed different scenarios with comparable coarse resolutions between 2 and 20 km, because our main focus was on the applicability of these approximations with respect to www.geosci-model-dev.net/7/1/2014/Geosci.Model Dev., 7, 1-21, 2014 large-scale modelling, possibly coupled with an atmosphere and ocean model in an Earth system model approach.For reversibility tests and transient experiments only the 5 km and 10 km resolutions were considered.In order to overcome the problem of capturing grounding line migration in coarse resolutions, we apply the heuristic rule described in Sect.4.4.Despite the rather coarse resolution (compared with most of the other 15 participating numerical ice models), RIMBAY accomplished the velocity-field comparison of the diagnostic experiment and in particular the reversibility test for grounding line migration.The latter is a prerequisite for modelling the Antarctic Ice Sheet, which is surrounded by ice shelves along more than half of its coast line.However, as a consequence of the imposed heuristic grounding line condition, which is not valid for a compressive flow, RIMBAY is (as all other A-HySSA models that participated in the MISMIP3d intercomparison) incapable of reproducing a grounding line retreat at free-slip walls.Despite the partly discontinuous grounding line retreat, resulting from the rather coarse resolution, we conclude that RIMBAY is capable of simulating large ice bodies with attached ice shelves for different climate conditions.

HOM and FS solvers
The numerical core for the HOM solver is very similar to the original implementation of Pattyn (2003), validated in the the Ice Sheet Model Intercomparison Project for higher order models (ISMIP-HOM) experiments (Pattyn et al., 2008).The FS implementation is basically an extension of this code and has originally been published in Pattyn (2008) for a linear rheology (with n = 1 in Eq. 8) and successfully been expanded for non-linear theologies (with n = 3) by Thoma et al. (2010Thoma et al. ( , 2012)).The results of these specific codefragments are already published, hence we do not present any additional validation of the FS solution here.However, we present two coupled SIA-FS-SSA experiments to demonstrate the flexibility of RIMBAY, with respect of a nested HOM/FS domain within a SIA-SSA domain.
The first experiment, is simply an extension of the original MISMIP experiment discussed in Sect.6.3 and Fig. 8.In addition to the coupled SIA-SSA solver, we modelled an area of 250 km in the vicinity of the GRL with the HOM and FS solver, respectively.We confined this test to the highest viscosity applied in Pattyn et al. (2012) (1/A = 2.1544 × 10 23 s Pa 3 ), where the modelled grounding line position is at about 1058.7 km (Fig. 8).Applying the HOM solver in the vicinity of the grounding line, after the SIA-SSA model reached a steady state, results in a slight retreat of about 8 km (which is below the grid size of 10 km) to 1051.3 km.This result is in close agreement with the theoretical position (1051.9km) according to the boundary layer theory of Schoof (2007).Switching from the HOM solver to the FS solver, however, does not change the GRL position anymore

HOM-and FS-Solver
The numerical core for the HOM-solver is very similar to the original implementation of Pattyn (2003), validated in the the 510 Ice Sheet Model Intercomparison Project for Higher-Order Models (ISMIP-HOM) experiments (Pattyn et al., 2008).The FS implementation is basically an extension of this code and has originally been published in Pattyn (2008) for a linear rheology (with n = 1 in Eq. 8) and successfully been ex-515 panded for nonlinear theologies (with n = 3) by Thoma et al. (2010Thoma et al. ( , 2012)).The results of these specific code-fragments are already published, hence we do not present any additional validation of the FS solution here.However, we present two coupled SIA-FS-SSA experiments to demonstrate the flexi-  significantly, neither does an extension of the HOM/FS domain from 250 km up to 900 km.
In the second experiment, the bedrock is downward sloping with a central trough: the horizontal resolution is 5 km and the accumulation is set to ṁac = 0.5 m a −1 (Fig. 9a).A Weertman-type sliding law (Eq.20) is applied as basal boundary condition, modified with an additional basal sliding reduction in the model's domain centre according to Geosci.Model Dev., 7, 1-21, 2014 www.geosci-model-dev.net/7/1/2014/ with C = 10 7 Pa m −1/3 s 1/3 , m = 1 3 , x j = 300 km, x i = 100 km, y j = 100 km, and y i = 10 km (this reduction is similar to those applied in Pattyn et al., 2013).The ice is surrounded by nunataks in the south, west, and north and an ocean in the east.We apply dS/dx = 0 at the western ice divide, free-slip boundary conditions along the southern and northern nunataks and the dynamic boundary conditions according to Eq. ( 25) at the ice-ocean boundary.First, the model is run with the coupled SIA-SSA solver and a transition zone of 50 km (imposing the heuristic condition outlined in Sect.4.4 at the grounding line) until a steady state is reached.Within the transition zone, the SIA and the SSA solutions for the velocity field are interpolated.The final steady state of this control experiment is shown in Fig. 9a, indicating the ice's geometry as well as the vertically averaged horizontal velocity.
Thereafter, (first) the HOM solver and (later) the FS solver are applied to the region, indicated in Fig. 9.As a result the grounding line advantages from 356 km to 398 km (HOM) and 408 km (FS solver), respectively, in this synthetic experiment.Pattyn et al. (2012Pattyn et al. ( , 2013) ) and Drouet et al. (2013) already discussed the limitations of the SIA-SSA approximations with respect to grounding line migration and pointed out that a high spatial resolution would be necessary to map the whole dynamic behaviour of transient states in ice sheet models.
However, in ice sheet/ice shelf models on a continental or even global scale and on long timescales (millennia), a high spatial resolution (below 10 km) and a FS solver, which consumes significant more computational resources than the SIA-SSA approximations, might be too ambitious at present.With respect to the large uncertainties of other atmospheric and ocean modelling issues, like boundary conditions and parameterisations, the drawback of the SIA-SSA approximation might be tolerable for most large-scale applications.

Conclusions
We have shown, that RIMBAY is capable of reproducing results of previously published experiments and benchmark tests (upper part of Table 1).In addition, RIMBAY has already been successfully applied in many very different scenarios in recent years.These applications range from highresolution FS modelling of ice flow across subglacial lakes (Thoma et al., 2010(Thoma et al., , 2012) ) in studies concerning the interaction between ice sheet, ice shelf and the ocean with a coupled SIA-SSA solver (Determann et al., 2012(Determann et al., , 2014;;Pattyn et al., 2013), to coupling RIMBAY with the Community Earth System Models (COSMOS) (Barbi et al., 2013) and the Viscoelastic Lithosphere and Mantle model (VILMA) (Konrad et al., 2013(Konrad et al., , 2014)), which calculates the isostatic adjustment of a spherical Earth to (ice-)surface loads.
Two additional modules are implemented within RIM-BAY, broadening its versatility: first, the water layer concept, developed by Goeller et al. (2013a), providing a sophisticated concept for the evolution of a large-scale subglacial hydrological network, which interacts with the ice sheet by modifying the basal boundary conditions.Second, a sub-grid scale Lagrangian-tracer module, allowing to track tracer propagation through the ice, which assists with the interpretation of the origin and age of ice cores (Sutter et al., 2013).All mentioned applications and modules are summarised in Table 1.
Several numerical ice flow model codes have been developed over the past years.In particular, the more recent FE approaches to solve the FS equations, as implemented in the ISSM (Larour et al., 2012;Seroussi et al., 2012) and ELMER/ICE (e.g.Zwinger et al., 2007;Gillet-Chaulet et al., 2012), are very promising.Their ability to adjust the spatial grid resolution with respect to the area under investigation is very useful.Their main drawbacks are currently the computational resources needed, which prohibit longterm projections on millennial timescales, and that (at least to the author's knowledge) there is no general concept of moving grids, which could adjust along a migrating area of investigation.To our knowledge only ISSM has the potential of coupling models of varying orders of complexity (SIA/SSA/HOM/FS) (Seroussi et al., 2012).
Well-known FD thermomechanically coupled ice sheet models are the Community Ice Sheet Models (CISM, based on GLIMMER; Payne, 1999;Rutt et al., 2009;Bougamont et al., 2011;Lemieux et al., 2011), the Parallel Ice Sheet Model (PISM; Martin et al., 2011;Winkelmann et al., 2011), the PennState3D ice-sheet/shelf model (PenState3D; Pollard and DeConto, 2012) and the SImulation COde for POLythermal Ice Sheets (SICOPOLIS; Greve, 1995;Sato and Greve, 2012).A recent overview about the most up-to-date ice sheet models is given by Bindschadler et al. (2013).All these models have proven their flexibility in several applications.However, none of these coupled SIA-SSA models has the option to simulate selected domains (e.g. the vicinity of grounding lines or ice streams) within a larger (e.g.continental scale) area with a (potentially migrating) FS approach.
With RIMBAY we provide a versatile open-source ice dynamics model to the scientific community.The code is not parallelised yet (apart from the LIS), and can be run even on common single-processor Linux or Unix systems.RIM-BAY can be redistributed and/or modified under the terms of the GNU General Public License as published by the Free Software Foundation, either 3rd (or any later) version of the license5 .RIMBAY fills a gap between several demands of the ice sheet modelling community, because it combines (a) the simplicity of a finite difference model, which can be run on a single processor, with (b) the option to model selected regions with a HOM or FS model, and (c) the potential to apply  hese coefficients represent the non-zero elements of a sinle row n for each ij -element of the matrix A nm , while b n epresents the right-hand side of Eq.A2.
Boundary conditions have to be formulated at the edges of he ice body.In case of open boundaries, the unknown value i+1 is virtually extrapolated from the interior.
n example for the ice thickness H is illustrated in Fig- re A1.Substitution of the ice thickness H, the velocity v i , he diffusion D i , and the ice bottom B on the eastern edge acording to Eq. A5 as well as on the southern edge according o ψ j−1 = 2ψ j − ψ j−1 in Eq.A2 results in the highlighted odifications If the ice adjoins a nunatak, closed b are applied.In this case the ice thicknes velocity u i+ 1 /2 , and the diffusion D i+ 1 /2 Alternatively we can express the unknow and B i−1 = ern edge results in the highlighted modifi to Eq. A2: These coefficients represent the non-zero elements of a single row n for each ij element of the matrix A nm , while b n represents the right-hand side of Eq. (A2).
Boundary conditions have to be formulated at the edges of the ice body.In case of open boundaries, the unknown value ψ i+1 is virtually extrapolated from the interior.

Defining
The lateral boundary condition along the ice shelf front is defined as free slip, and hence ∂U ∂y = 0 ∂V ∂x = 0, (B10) or N: U i+1,j = U i,j ⇒ γ u 1 = 0 in C u 7 and C u 5 , S: U i−1,j = U i,j ⇒ γ u 2 = 0 in C u 3 and C u 5 , E: V i,j +1 = V i,j ⇒ γ v 1 = 0 in C v 6 and C v 5 , W: V i,j −1 = V i,j ⇒ γ v 2 = 0 in C v 4 and C v 5 .In case of no slip boundary conditions, the flow on the boundary is zero, and hence u i+1,j + u i,j 2 = 0 v i+1,j + v i,j 2 = 0 (B11) or N: u i+1,j = −u i,j ⇒ γ u 1 = 0 in C u 7 and 2 in C u 5 , S: u i−1,j = −u i,j ⇒ γ u 2 = 0 in C u 3 and 2 in C u 5 , E: v i,j +1 = −v i,j ⇒ γ v 1 = 0 in C v 6 and 2 in C v 5 , W: v i,j −1 = −v i,j ⇒ γ v 2 = 0 in C v 4 and 2 in C v 5 .2012) * : about 30 percent of the whole domain (the lake and its vicinity) is calculated with a FS solver, note that there is no ice evolution in this experiment, only the ice temperature evolves over time.
with the basal tangential stress component τ b = τ • n b and the normal vector n b orthogonal to the ice base, the basal friction coefficient C, and the basal friction exponent m.

240-
the differences for the two storage formats (CRS vs. MSR) are negligible, the differences between the linbcg solver from Press et al. (2007) and the very same preconditioner/solver combination from the Library of Iterative Solvers are 245 negligible, but the solver of Press et al. (2007) needs less computational resources.if specific preconditioner-solver combination converge, the difference between different combinations are negligible.250

Fig. 2 .
Fig. 2.Sequence of iteratively solved variables within RIMBAY.In the SSA case the product of H η is calculated, instead of the viscosity η, only.The grayish highlighted variables are calculated only in the FS case.(Within the main loop, the vertical velocity needs only to be calculated in the non-FS cases.)

Fig. 3 .
Fig. 3. Location of nodes on a C-Grid.The location of H-, η-, and θ-nodes are indicated by dots while the location of the horizontal velocities are indicated by arrows.The stars indicate certain intergrid nodes used in the numerical implementation.The red color indicates corresponding nodes of the central i,j-node and the colorcoded increments (∆...) at the edges refer to the corresponding grid node distances.

Fig
Fig. 4. R first order

Fig. 3 .
Fig. 3. Location of nodes on a C-grid.The location of H , η, and θ nodes are indicated by dots while the location of the horizontal velocities are indicated by arrows.The stars indicate certain inter-grid nodes used in the numerical implementation.The red colour indicates corresponding nodes of the central i,j node and the colourcoded increments ( . ..) at the edges refer to the corresponding grid node distances.

Fig. 4 .
Fig. 4. Relative positions and numbering of nodes for the implicit first order finite difference formulation of the ice evolution (Eq.31).

Fig. 6 .Fig. 7 .
Fig. 6.Modelled horizontal velocity for a synthetic iceberg in different orientations.The enlarged inlet shows exemplarily the orientation of the normal vectors at the ice shelf front in green (compare Eq.25).

Fig. 6 .Fig. 6 .Fig. 7 .
Fig. 6.Modelled horizontal velocity for a synthetic iceberg in different orientations.The enlarged inlet shows exemplarily the orientation of the normal vectors at the ice shelf front in green (compare Eq. 25).

Fig. 7 .
Fig. 7. Modelled horizontal velocity for two synthetic floating ice structures of complex geometries.Nunataks are indicated in brown.At the southern (lower) edge no-slip boundary conditions are applied, at the northern edge and at the ice rise in the left ice body free-slip boundaries are valid.

Fig. 9 .
Fig. 9. Geometry for the experiment described in section 6.4.a) Bedrock topography and ice geomtery.The horizontal ice velocity is plotted on top of the ice sheet surface; the magenta and red lines indicate the interpolated (sub-grid scale) GRL-positions for the coupled SIA/SSA, the HOM-(dotted) and the FS-(solid) solution, respectively; the black rectangle indicates the region, where the FS-solver is applied.Additionally, the basal friction parameter β 2 (according to Eq. 20) is shown.b) Profile along y = 100 km.The dashed black lines indicate the area where the HOM and FS solutions are calculated, respectively; the red lines indicate the shape of the corresponding ice geometry for the HOM-solution (dotted) and FS-solution (solid).

Fig. 9 .
Fig. 9. Geometry for the experiment described in Sect.6.4.(a) Bedrock topography and ice geomtery.The horizontal ice velocity is plotted on top of the ice sheet surface; the magenta and red lines indicate the interpolated (sub-grid scale) GRL-positions for the coupled SIA-SSA, the HOM (dotted) and the FS (solid) solutions, respectively; the black rectangle indicates the region, where the FS solver is applied.Additionally, the basal friction parameter β 2 (according to Eq. 20) is shown.(b) Profile along y = 100 km.The dashed black lines indicate the area where the HOM and FS solutions are calculated, respectively; the red lines indicate the shape of the corresponding ice geometry for the HOM solution (dotted) and FS solution (solid).

,
Fig. B1.Relative positions and numbering of nodes for the SSA.
Relative positions and numbering of nodes for the implicit first order finite difference formulation of the ice evolution (Eq.31).

Table 1 .
Pattyn (2008)(2010)ications of the ice sheet/shelf model RIMBAY.The FS experiments appear as validation as well, as there is no explicit benchmark available.The two FS validations differ, asThoma et al. (2010)extended the originally linear flow law used byPattyn (2008).

Table 2 .
Acronyms.thenumericalmodel (with small effort) to new synthetic or realistic scenarios.Typical computational times needed by RIMBAY for some representative applications are given in Appendix C.Based on the specific needs, RIMBAY can be applied easily to new scientific issues and is easy extensible because of well-structured interfaces.It provides a broad spectrum of applicability and functionality and could therefore contribute to solve the pressing questions of global climate change.M. Thoma et al.: Description of the ice fl Geosci.Model Dev., 7, 1-21, 2014www.geosci-model-dev.net/7/1/2014/ Relative positions and numbering of nodes for the SSA.

Table C1 .
Typical CPU wall-clock times consumed to simulate typical domains with the specified numerical complexities.