Dynamically coupling full Stokes and shallow shelf approximation for marine ice sheet flow using Elmer/Ice (v8.3)

Abstract

Ice flow forced by gravity is governed by the full Stokes (FS) equations, which are computationally expensive to solve due to the nonlinearity introduced by the rheology. Therefore, approximations to the FS equations are commonly used, especially when modeling a marine ice sheet (ice sheet, ice shelf, and/or ice stream) for ${\mathrm{10}}^{\mathrm{3}}$ years or longer. The shallow ice approximation (SIA) and shallow shelf approximation (SSA) are commonly used but are accurate only for certain parts of an ice sheet. Here, we report a novel way of iteratively coupling FS and SSA that has been implemented in Elmer/Ice and applied to conceptual marine ice sheets. The FS–SSA coupling appears to be very accurate; the relative error in velocity compared to FS is below 0.5 % for diagnostic runs and below 5 % for prognostic runs. Results for grounding line dynamics obtained with the FS–SSA coupling are similar to those obtained from an FS model in an experiment with a periodical temperature forcing over 3000 years that induces grounding line advance and retreat. The rapid convergence of the FS–SSA coupling shows a large potential for reducing computation time, such that modeling a marine ice sheet for thousands of years should become feasible in the near future. Despite inefficient matrix assembly in the current implementation, computation time is reduced by 32 %, when the coupling is applied to a 3-D ice shelf.

Introduction

Dynamical changes in both the Greenland and Antarctic ice sheets are, with medium confidence, projected to contribute 0.03 to 0.20 m of sea level rise by 2081–2100 . The main reason for the uncertainty in these estimates is a limited understanding of ice dynamics. Thus, there is a great need for improvement of ice dynamical models . The gravity-driven flow of ice is described by the full Stokes (FS) equations, amended by a nonlinear rheology described by Glen's flow law. Model validation is required over centennial to millennial timescales to capture the long response time of an ice sheet to external forcing . However, the computation time and memory required for an FS model to be applied to ice sheets restricts simulations to sub-millennial timescales . Therefore, approximations of the FS equations are employed for simulations over long timescales, such as the shallow ice approximation SIA;, the shallow shelf approximation SSA;, Blatter–Pattyn , and hybrid models .

Any ice sheet model accounting for ice shelves needs to resolve grounding line dynamics (GLD). Despite many recent efforts, modeling GLD still poses a challenge in numerical models, as illustrated by the wide range of results obtained in the Marine Ice Sheet Model Intercomparison Project MISMIP;. In MISMIP3d, GLD differ between FS models and SSA models, with discrepancies attributed to so-called higher-order terms which are neglected in SSA models but included in FS models . Based on these model intercomparisons, it is advised to use models that include vertical shearing to compute reliable projections of ice sheet contribution to sea level rise . On the other hand, it is not entirely clear how much of the difference in GLD is due to the different numerical treatment of the grounding line problem in shallow models. An updated version of the hybrid SIA/SSA Parallel Ice Sheet Model (PISM) that uses a modified driving stress calculation and subgrid grounding line interpolation showed GLD comparable to an FS model . It should be noted that the experiments in MISMIP3d consisted of laterally extruded idealized 2-D geometries with quite small sideward disturbances, and MISMIP$+$ may give more insight on realistic situations. Additionally, there is a recent publication that sheds new light on a possible problem with the setup of MISMIP experiments .

Solving the FS equations over large spatiotemporal domains is still infeasible. However, solvers combining approximations (e.g., SIA or SSA) with the FS equations allow the simulation of ice dynamics over long time spans without introducing artifacts caused by the application of approximations in parts of the domain where they are not valid. For instance, coupled FS and SSA, in the framework of the Ice Sheet System Model ISSM;. They apply the tiling method which includes a blending zone of FS and SSA. Their result looks promising with respect to both accuracy and efficiency but is limited to diagnostic experiments. The Ice Sheet Coupled Approximation Levels (ISCAL) method couples SIA and FS by a nonoverlapping domain decomposition that dynamically changes with time. ISCAL is implemented in Elmer/Ice , an open-source finite element software for ice sheet modeling. Here, we present a novel coupling between FS and SSA, also by the implementation of a nonoverlapping domain decomposition in Elmer/Ice. The domain decomposition changes dynamically with grounding line advance and retreat. GLD are modeled with FS and coupled to SSA on the ice shelf via boundary conditions. The equations discretized by the finite element method (FEM) are solved iteratively, alternating between the FS and the SSA domain, until convergence is reached.

The extent of present-day ice shelves is limited to approximately 10 % of the area of Antarctica . Therefore, one may question the reduction in computational work by applying SSA to model ice shelves in continental-scale simulations of marine ice sheets. However, the coupling is targeted at conducting paleo-simulations, for which much larger ice shelves have been present . In that case, a large part of the interior of a marine ice sheet is modeled with SIA, SSA is applied to the ice shelves, and the FS domain is restricted to ice streams and areas around the grounding line.

An overview of the FS and SSA equations governing ice sheet and shelf dynamics in three dimensions (3-D) is presented in Sect. , together with the boundary conditions. Memory and performance estimates of an FS–SSA coupling, independent of the specific coupling implemented, are provided in Sect. . Section  describes the coupled FS–SSA model, hereafter “coupled model”. The coupling is applied to a conceptual ice shelf ramp and marine ice sheet in Sect. . The simulation of a 3000-year long cycle of grounding line advance and retreat (described in Sect. ) shows the robustness of the coupling.

Governing equations of ice flow

Ice is considered to be an incompressible fluid, such that mass conservation implies that the velocity is divergence-free: $$\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0},$$ where $\mathit{u}=(u,v,w{)}^T$ describes the velocity field of the ice with respect to a Cartesian coordinate system $(x,y,z{)}^T$, where $z$ is the vertical direction. For ice flow, the acceleration term can be neglected in the Navier–Stokes equations . Therefore, the conservation of linear momentum under the action of gravity $\mathit{g}$ can be described by $$-\mathrm{\nabla }p+\mathrm{\nabla }\cdot \left(\mathit{\eta }(\mathrm{\nabla }\mathit{u}+(\mathrm{\nabla }\mathit{u}{)}^T)\right)+\mathit{\rho }\mathit{g}=\mathbf{0},$$ where $\mathrm{\nabla }$ is the gradient operator, $p$ pressure, $\mathit{\eta }$ viscosity, $\mathit{\rho }$ ice density, and $\mathit{g}$ gravity. Letting $\mathit{\sigma }$ denote the stress tensor, pressure $p$ is the mean normal stress ($p=-\mathrm{1}/\mathrm{3}{\mathrm{\Sigma }}_{\mathrm{i}}{\mathit{\sigma }}_{ii}$) and $\mathbf{D}\left(\mathit{u}\right)$ is the strain rate tensor, related by $$\mathit{\sigma }=\mathrm{2}\mathit{\eta }\mathbf{D}\left(\mathit{u}\right)-p\mathbf{I}=\mathit{\eta }(\mathrm{\nabla }\mathit{u}+(\mathrm{\nabla }\mathit{u}{)}^T)-p\mathbf{I},$$ where $\mathbf{I}$ is the identity tensor. Together, Eqs. () and () are called the full Stokes equations. Observations by suggest that the viscosity depends on temperature $T$ and the effective strain rate $D\left(\mathit{u}\right)$: $$\begin{array}{cc}& {\displaystyle }& {\displaystyle }\mathit{\eta }(\mathit{u},T)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathcal{A}\left(T{)}^{-\frac{\mathrm{1}}{n}}D\right(\mathit{u}{)}^{\frac{\mathrm{1}-n}{n}},& {\displaystyle }& {\displaystyle }D\left(\mathit{u}\right)=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\left({\displaystyle \frac{\partial u}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial v}{\partial y}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial w}{\partial z}}\right)}^{\mathrm{2}}\right)}{\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left({\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial u}{\partial z}}+{\displaystyle \frac{\partial w}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial v}{\partial z}}+{\displaystyle \frac{\partial w}{\partial y}}\right)}^{\mathrm{2}}\right)},\end{array}$$ where Glen's exponent $n=\mathrm{3}$. The fluidity parameter $\mathcal{A}$ increases exponentially with temperature as described by the Arrhenius relation . This represents a thermodynamically coupled system of equations. However, in the current study, we focus on the mechanical effects and a uniform temperature is assumed. Due to the velocity dependence of the viscosity in Eq. (), the FS equations form a nonlinear system with four coupled unknowns, which is time consuming to solve. Therefore, many approximations to the FS equations have been derived in order to model ice sheet dynamics on long timescales; see Sect. .

Shallow shelf approximation

Floating ice does not experience basal drag, hence all resistance comes from longitudinal stresses or lateral drag at the margins. For ice shelves, the SSA has been derived by dimensional analysis based on a small aspect ratio and surface slope . This dimensional analysis shows that vertical variation in $u$ and $v$ is negligible, such that $w$ and $p$ can be eliminated by integrating the remaining stresses over the vertical and applying the boundary conditions at the glacier surface and base (described in Sect. ). Then, the conservation of linear momentum, Eq. (), simplifies to $${\mathrm{\nabla }}_{\mathrm{h}}\cdot \left(\mathrm{2}\stackrel{\mathrm{‾}}{\mathit{\eta }}\left({\mathbf{D}}_{\mathrm{h}}\left(\mathbf{u}\right)+\text{tr}\left({\mathbf{D}}_{\mathrm{h}}\right(\mathbf{u}\left)\right)\mathbf{I}\right)\right)=\mathit{\rho }gH{\mathrm{\nabla }}_{\mathrm{h}}{z}_{\mathrm{s}},$$ where the subscript h represents the components in the $x$–$y$ plane, $\stackrel{\mathrm{‾}}{\mathit{\eta }}$ the vertically integrated viscosity, $H$ the thickness of the ice shelf, and $z_{\mathrm{s}}$ the upper ice surface; see Fig. . The effective strain rate in Eq. () simplifies to $$\begin{array}{cc}& {\displaystyle }& {\displaystyle }{D}_{\mathrm{h}}\left(\mathit{u}\right)={\displaystyle }& {\displaystyle }\sqrt{{\left({\displaystyle \frac{\partial u}{\partial x}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial v}{\partial y}}\right)}^{\mathrm{2}}+{\displaystyle \frac{\partial u}{\partial x}}{\displaystyle \frac{\partial v}{\partial y}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}{\left({\displaystyle \frac{\partial u}{\partial y}}+{\displaystyle \frac{\partial v}{\partial x}}\right)}^{\mathrm{2}}},\end{array}$$ where $w$ is eliminated using incompressibility; Eq. (). The SSA equations are still nonlinear through $\stackrel{\mathrm{‾}}{\mathit{\eta }}$, but since $w$ and $p$ are eliminated and vertical variation in $u$ and $v$ is neglected, the 3-D problem with four unknowns is reduced to a 2-D problem with two unknowns. Therefore, the SSA model is less computationally demanding than FS. The horizontal velocities are often of main interest, for example when results are validated by comparison to observed horizontal surface velocity. If desirable, the vertical velocity can be computed from the incompressibility condition.

Overview of the notations and domain decomposition for a conceptual marine ice sheet. The vertical scale is exaggerated. The sea level at $z=\mathrm{0}$ is dashed blue and the interface between the FS and SSA domains is solid red. The bed elevation is denoted by $b$, the coupling interface by $x_{\mathrm{c}}$, and the grounding line by $x_{\mathrm{GL}}$. The distance between $x_{\mathrm{c}}$ and $x_{\mathrm{GL}}$, defined in Eq. (), is denoted by $d_{\mathrm{GL}}$.

[Figure omitted. See PDF]

Boundary conditions and time evolution

The coupling is applied to a marine ice sheet, with bedrock lying (partly) below sea level (see Fig. ), and involves boundaries in contact with the bedrock, ocean and atmosphere. The only time dependency is in the evolution of the free surfaces.

Bedrock

Where the ice is grounded (in contact with the bedrock), the interaction of ice with the bedrock is commonly represented by a sliding law $f(\mathit{u},N)$, that relates the basal velocity ${\mathit{u}}_{\mathrm{b}}$ and effective pressure $N$ to the basal shear stress as $$\begin{array}{}& {\displaystyle }& {\displaystyle }({\mathit{t}}_{\mathrm{i}}\cdot \mathit{\sigma }\cdot \mathit{n}{)}_{\mathrm{b}}=f(\mathit{u},N)\mathit{u}\cdot {\mathit{t}}_{\mathrm{i}},i=\mathrm{1},\mathrm{2},& {\displaystyle }& {\displaystyle }(\mathit{u}\cdot \mathit{n}{)}_{\mathrm{b}}+a_{\mathrm{b}}=\mathrm{0},\end{array}$$ where ${\mathit{t}}_{\mathrm{i}}$ are the vectors spanning the tangential plane, $\mathit{n}$ is the normal to the bed, and $a_{\mathrm{b}}$ describes basal refreezing or melt. A sliding law suggested by is assumed, which depends on ${\mathit{u}}_{\mathrm{b}}$ and the height above buoyancy $z_{*}$ such that $${\displaystyle }f(\mathit{u},N)=-\mathit{\beta }\left|{\mathit{u}}_{\mathrm{b}}{|}^{\frac{\mathrm{1}}{n}-\mathrm{1}}{z}_{*}\right(N).$$ Here, the sliding parameter $\mathit{\beta }$ is constant in time and space. In line with , instead of modeling $N$, a hydrostatic balance is assumed to approximate $z_{*}$, implying a subglacial hydrology system entirely in contact with the ocean: $$z_{*}\left(H\right)=\left\{\begin{array}{ll}H& \text{if }{z}_{\mathrm{b}}\ge \mathrm{0},\\ H+z_{\mathrm{b}}{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{\mathit{\rho }}}& \text{if }{z}_{\mathrm{b}}<\mathrm{0},\end{array}\right.$$ where $z_{\mathrm{b}}$ is the lower ice surface and ${\mathit{\rho }}_{\mathrm{w}}$ the water density and the sea level is at $z=\mathrm{0}$. Equation () implies that $z_{*}$ equals zero when the flotation criterion (Archimedes' principle) is satisfied, i.e., where $$z_{\mathrm{s}}=\left(\mathrm{1}-{\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{w}}}}\right)H,z_{\mathrm{b}}=-{\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{w}}}}H.$$

Ice–ocean interface

As soon as the seawater pressure $p_{\mathrm{w}}$ at the ice base $z_{\mathrm{b}}$ is larger than the normal stress exerted by the ice at the bed, the ice is assumed to float. For a detailed description of the implementation of the contact problem at the grounding line in Elmer/Ice, see . At the ice–ocean interface, the tangential friction is neglected ($f(\mathit{u},N)\equiv \mathrm{0}$ in Eq. ) and $$\mathit{\sigma }\cdot \mathit{n}=-p_{\mathrm{w}}\mathit{n}\text{ where }{p}_{\mathrm{w}}\left(z\right)=-{\mathit{\rho }}_{\mathrm{w}}gz\text{ if }z\le \mathrm{0},$$ and $\mathit{\sigma }\cdot \mathit{n}=\mathrm{0}$ above sea level ($z>\mathrm{0}$). Calving at the seaward front of the ice shelf is not explicitly modeled, but the length of the modeling domain is fixed and ice flow from the shelf out of the domain is interpreted as a calving rate.

Surface evolution

Ice surface (assumed stress-free, $\mathit{\sigma }\cdot \mathit{n}=\mathrm{0}$) and ice base at $z_{\mathrm{s}}$ and $z_{\mathrm{b}}$ behave as free surfaces according to $${\displaystyle \frac{\partial z_{\mathrm{s}/\mathrm{b}}}{\partial t}}+u_{\mathrm{s}/\mathrm{b}}{\displaystyle \frac{\partial z_{\mathrm{s}/\mathrm{b}}}{\partial x}}+v_{\mathrm{s}/\mathrm{b}}{\displaystyle \frac{\partial z_{\mathrm{s}/\mathrm{b}}}{\partial y}}=w_{\mathrm{s}/\mathrm{b}}+a_{\mathrm{s}/\mathrm{b}},$$ where $a_{\mathrm{s}/\mathrm{b}}$ is the accumulation ($a_{\mathrm{s}/\mathrm{b}}>\mathrm{0}$) or ablation ($a_{\mathrm{s}/\mathrm{b}}<\mathrm{0}$) in meter ice equivalent per year, at the surface or base, respectively. By vertical integration of the incompressibility condition, Eq. (), $w$ can be eliminated using Leibniz integration rule and substituting the free surface equations (Eq. ), which yields the thickness advection equation: $${\displaystyle \frac{\partial H}{\partial t}}+{\displaystyle \frac{\partial H\stackrel{\mathrm{‾}}{u}}{\partial x}}+{\displaystyle \frac{\partial H\stackrel{\mathrm{‾}}{v}}{\partial y}}=a_{\mathrm{s}}-a_{\mathrm{b}},$$ where $\stackrel{\mathrm{‾}}{u}$ and $\stackrel{\mathrm{‾}}{v}$ are the vertically integrated horizontal velocities.

Memory and performance estimates of an FS–SSA coupling

The reduction in the memory required for an FS–SSA coupling by domain decomposition, compared to an FS model, can be estimated. This estimate is independent of the specific implementation of the coupling between the domains and concerns only the most ideal implementation in which no redundant information is stored. The main advantage of the SSA model is that ${\mathit{u}}_{\mathrm{SSA}}$ is independent of $z$, such that the SSA equations can be solved on a part of the domain with a mesh of 1 dimension fewer. Besides that, there are fewer unknowns since $p$ and $w$ are eliminated. An additional advantage of eliminating $p$ is that the resulting system is mathematically easier to discretize and solve. In particular, difficulties related to a stable choice for the basis functions for the pressures and velocities are avoided see, e.g., and there is no need for specialized iterative solution techniques to solve the so-called saddle-point problem that the FS equations pose see.

Suppose that the computational domain $\mathrm{\Omega }$ is discretized with $N_z$ nodes regularly placed in the $z$ direction and $N_{\mathrm{h}}$ nodes in a horizontal footprint mesh and is decomposed into two parts (${\mathrm{\Omega }}_{\mathrm{SSA}}$ and ${\mathrm{\Omega }}_{\mathrm{FS}}$; see Fig. ). The fraction of nodes in ${\mathrm{\Omega }}_{\mathrm{SSA}}$ is denoted as $\mathit{\theta }$ with $\mathrm{0}<\mathit{\theta }<\mathrm{1}$. The number of nodes in ${\mathrm{\Omega }}_{\mathrm{FS}}$ is then approximately $(\mathrm{1}-\mathit{\theta })N_{\mathrm{h}}{N}_z$, and in ${\mathrm{\Omega }}_{\mathrm{SSA}}$, it is $\mathit{\theta }{N}_{\mathrm{h}}$, neglecting shared nodes on the boundary. For a 3-D physical domain, SSA has two unknowns ($u$ and $v$) and FS has four unknowns ($u$, $v$, $w$, and $p$). Hence, the memory needed to store the solution with a coupled model is proportional to $\mathrm{2}{N}_{\mathrm{h}}(\mathit{\theta }+\mathrm{2}(\mathrm{1}-\mathit{\theta }\left)N_z\right)$. For a 2-D simulation in the $x-z$ plane, where FS has three unknowns and SSA only one, the memory is proportional to $N_{\mathrm{h}}(\mathit{\theta }+\mathrm{3}(\mathrm{1}-\mathit{\theta }\left)N_z\right)$. The memory requirement for a physical domain in $d$ dimensions reduces to $$\begin{array}{cc}{\displaystyle }& {\displaystyle }{q}_{\mathrm{var}}={\displaystyle \frac{\text{coupled model memory}}{\text{FS model memory}}}=\mathrm{1}-\mathit{\theta }+{\displaystyle \frac{\mathit{\theta }}{(\mathrm{5}-d)N_z}},\\ & {\displaystyle }& {\displaystyle }d=\mathrm{2},\mathrm{3},\end{array}$$ when part of the domain is modeled by the SSA equations. The memory requirements for mesh-related quantities reduce to $q_{\mathrm{mesh}}=\mathrm{1}-\mathit{\theta }+\mathit{\theta }/N_z$ in both 2-D and 3-D. The quotients $q_{\mathrm{var}}$ and $q_{\mathrm{mesh}}$ are close to $\mathrm{1}-\mathit{\theta }$ if $N_z\mathit{\gtrsim }\mathrm{10}$.

The computational work is more difficult to estimate a priori since it depends on the implementation of the coupling. The dominant costs are for the assembly of the finite element matrices, the solution of the nonlinear equations, and an overhead for administration in the solver. The work to assemble the matrices grows linearly with the number of unknown variables. Suppose that this work for FS in 3-D is $\mathrm{4}{C}_{\mathrm{FS}}{N}_{\mathrm{h}}{N}_z$ in the whole domain, for FS $\mathrm{4}{C}_{\mathrm{FS}}(\mathrm{1}-\mathit{\theta })N_{\mathrm{h}}{N}_z$ in ${\mathrm{\Omega }}_{\mathrm{FS}}$, and for SSA $\mathrm{2}{C}_{\mathrm{SSA}}\mathit{\theta }{N}_{\mathrm{h}}$ in ${\mathrm{\Omega }}_{\mathrm{SSA}}$. The coefficients $C_{\mathrm{FS}}$ and $C_{\mathrm{SSA}}$ depend on the basis functions for FS and SSA and the complexity of the equations. The reduction in assembly time for the matrix is $q_{\mathrm{ass}}=\mathrm{1}-\mathit{\theta }+C_{\mathrm{SSA}}\mathit{\theta }/\mathrm{2}{C}_{\mathrm{FS}}{N}_z$. If $C_{\mathrm{FS}}\approx C_{\mathrm{SSA}}$, then the reduction is approximately as in Eq. (). The same conclusion holds in 2-D. Therefore, the reduction in that part is estimated to be similar to the reduction in Eq. ().

[Figure omitted. See PDF]

Method for coupling FS and SSA

All equations are solved in Elmer/Ice using the FEM. First the velocity $\mathit{u}$ (using FS or SSA) is solved for a fixed geometry at time $t$. The mesh always has the same dimension as the physical modeling domain, but ${\mathit{u}}_{\mathrm{SSA}}$ is only solved on the basal mesh layer, after which the solution is re-projected over the vertical axis. Then, the geometry is adjusted by solving the free surface and thickness advection equations using backward Euler time integration. The nonlinear FS and SSA equations are solved using a Picard iteration. The discretized FS equations are stabilized by the residual free bubbles method , the recommended stabilization method in . First, the coupling for a given geometry is presented, followed by the coupled surface evolution, both summarized in Algorithm 1.

The FS domain ${\mathrm{\Omega }}_{\mathrm{FS}}$ contains the grounded ice and a part of the shelf around the grounding line; see Fig. . The SSA domain ${\mathrm{\Omega }}_{\mathrm{SSA}}$ is restricted to a part of the ice shelf and starts at the coupling interface ${\mathit{x}}_{\mathrm{c}}$ at the first basal mesh nodes located at least a distance $d_{\mathrm{GL}}$ from the grounding line ${\mathit{x}}_{\mathrm{GL}}$, such that $$\begin{array}{cc}{\displaystyle }\left|\right|\mathit{x}-{\mathit{x}}_{\mathrm{GL}}\left|\right|& {\displaystyle }:=\sqrt{(x-x_{\mathrm{GL}}{)}^{\mathrm{2}}+(y-y_{\mathrm{GL}}{)}^{\mathrm{2}}+(z-z_{\mathrm{GL}}{)}^{\mathrm{2}}}\\ & {\displaystyle }& {\displaystyle }\ge d_{\mathrm{GL}}\text{ for all }\mathit{x}\text{ in }{\mathrm{\Omega }}_{\mathrm{SSA}}.\end{array}$$

Boundary conditions at the coupling interface

Horizontal gradients of the velocity are not neglected in the SSA equations unlike in the SIA;. Thus, not only FS and SSA velocities but also their gradients have to match, in order to allow a coupling of the two. Therefore, one cannot solve one system of equations independently, for use as an input to the other system, as done for a one-way coupling e.g.,. Instead, the coupling of FS and SSA is solved iteratively, updating the interaction between FS and SSA velocities in each iteration to obtain mutually consistent results. SSA-governed ice shelf flow is greatly influenced by the inflow velocity from the FS domain. Therefore, we start the first iteration of the coupled model by solving the FS equations. A boundary condition is necessary at ${\mathit{x}}_{\mathrm{c}}$; we assume that the cryostatic pressure acts on ${\mathrm{\Omega }}_{\mathrm{FS}}$ at ${\mathit{x}}_{\mathrm{c}}$: $${\mathit{\sigma }}_{\mathrm{FS}}\cdot \mathit{n}({\mathit{x}}_{\mathrm{c}},z)=\mathit{\rho }g(z_{\mathrm{s}}-z)\mathit{n},$$ where $\mathit{n}$ is normal to the coupling interface ${\mathit{x}}_{\mathrm{c}}$. The FS velocity at $x_{\mathrm{c}}$ provides a Dirichlet inflow boundary condition to the SSA equations. Then, the Neumann boundary condition in Eq. () has to be adjusted based on the ice flow as calculated for ${\mathrm{\Omega }}_{\mathrm{SSA}}$. This is done using the contact force denoted by ${\mathit{f}}_{\mathrm{SSA}}$, as explained below.

The SSA equations are linearized and, by means of FEM, discretized. This leads to a matrix representation $\mathbf{A}\mathit{u}=\mathit{b}$, where $\mathit{u}$ is the vector of unknown variables (here, horizontal SSA velocities). In FEM terminology, the vector $\mathit{b}$ that describes the forces driving or resisting ice flow is usually called the body force and $\mathbf{A}$ the system matrix . In Elmer/Ice, Dirichlet conditions for a node $i$ are prescribed by setting the $i$th row of $\mathbf{A}$ to zero, except for the diagonal entry which is set to be unity, and ${\mathit{b}}_{\mathrm{i}}$ is set to have the desired value . For an exact solution of $\mathbf{A}\mathit{u}=\mathit{b}$, the residual $\mathit{f}=\mathbf{A}\mathit{u}-\mathit{b}$ is zero. If we instead use the system matrix ${\mathbf{A}}_{\mathrm{SSA}}$ obtained without the Dirichlet conditions being set, the resulting residual is equal to the contact force that would have been necessary to produce the velocity described by the Dirichlet boundary condition. Since the SSA equations are vertically integrated, ${\mathit{f}}_{\mathrm{SSA}}={\mathbf{A}}_{\mathrm{SSA}}{\mathit{u}}_{\mathrm{SSA}}-{\mathit{b}}_{\mathrm{SSA}}$ is the vertically integrated contact force and needs to be scaled by the ice thickness $H$. In Elmer/Ice, ${\mathit{f}}_{\mathrm{SSA}}$ is mesh dependent and needs to be scaled by the horizontal mesh resolution $\mathit{\omega }$ as well. For 2-D configurations, $\mathit{\omega }=\mathrm{1}$. Using ${\mathit{f}}_{\mathrm{SSA}}$ instead of explicitly calculating the stress is advantageous since it is extremely cheap to find the contact force if ${\mathbf{A}}_{\mathrm{SSA}}$ is stored.

To summarize the boundary conditions at ${\mathit{x}}_{\mathrm{c}}$, for FS, an external pressure is applied: $${\mathit{\sigma }}_{\mathrm{FS}}\cdot \mathit{n}({\mathit{x}}_{\mathrm{c}},z)=\mathit{\rho }g(z-z_{\mathrm{s}})\mathit{n}+{\displaystyle \frac{{\mathit{f}}_{\mathrm{SSA}}\left({\mathit{x}}_{\mathrm{c}}\right)}{\mathit{\omega }H}},$$ where ${\mathit{f}}_{\mathrm{SSA}}:=\mathrm{0}$ in the first iteration (for its derivation, see Appendix ). For SSA, a Dirichlet inflow boundary condition, $${\mathit{u}}_{\mathrm{SSA}}\left({\mathit{x}}_{\mathrm{c}}\right)={\mathit{u}}_{\mathrm{FS}}({\mathit{x}}_{\mathrm{c}},z_{\mathrm{b}}),$$ provides the coupling to the FS solution. Here we take the ${\mathit{u}}_{\mathrm{FS}}$ at $z_{\mathrm{b}}$, but any $z$ can be chosen since ${\mathit{x}}_{\mathrm{c}}$ should be located such that ${\mathit{u}}_{\mathrm{FS}}({\mathit{x}}_{\mathrm{c}},z)$ hardly varies with $z$. Every iteration, ${\mathit{f}}_{\mathrm{SSA}}$, and ${\mathit{u}}_{\mathrm{FS}}({\mathit{x}}_{\mathrm{c}},z_{\mathrm{b}})$ are updated until convergence up to a tolerance ${\mathit{\epsilon }}_{\mathrm{c}}$.

Surface evolution

The surface evolution is calculated differently in the two domains ${\mathrm{\Omega }}_{\mathrm{FS}}$ and ${\mathrm{\Omega }}_{\mathrm{SSA}}$. Equation () is applied to ${\mathrm{\Omega }}_{\mathrm{FS}}$ for the evolution of $z_{\mathrm{s}}$ and $z_{\mathrm{b}}$, avoiding assuming hydrostatic equilibrium beyond the grounding line, since the flotation criterion is not necessarily fulfilled close to the grounding line . The thickness advection equation, Eq. (), is used for ${\mathrm{\Omega }}_{\mathrm{SSA}}$, which is advantageous since the ice flux $\mathit{q}=H{\mathit{u}}_{\mathrm{SSA}}$ is directly available (because ${\mathit{u}}_{\mathrm{SSA}}$ does not vary with $z$) and no vertical velocity is needed. Moreover, only one time-dependent equation is solved instead of one for the lower and one for the upper free surface. The evolution of the surfaces $z_{\mathrm{s}}$ and $z_{\mathrm{b}}$ for ${\mathrm{\Omega }}_{\mathrm{SSA}}$ is then calculated from the flotation criterion, Eq. (). At ${\mathit{x}}_{\mathrm{c}}$, $H_{\mathrm{SSA}}=H_{\mathrm{FS}}$ is applied as a boundary condition to the thickness equation. First the surface evolution is solved for ${\mathrm{\Omega }}_{\mathrm{FS}}$; then ${\mathrm{\Omega }}_{\mathrm{SSA}}$ follows.

The algorithm

The iterative coupling for one time step is given by Algorithm 1. First, the shortest distance $d$ to the grounding line is computed for all nodes in the horizontal footprint mesh at the ice shelf base. Then, a mask is defined that describes whether a node is in ${\mathrm{\Omega }}_{\mathrm{FS}}$, ${\mathrm{\Omega }}_{\mathrm{SSA}}$, or at the coupling interface ${\mathit{x}}_{\mathrm{c}}$, based on the user-defined $d_{\mathrm{GL}}$. Technically, the domain decomposition is based on the use of passive elements implemented in the overarching Elmer code , which allow for deactivating and reactivating elements. An element in ${\mathrm{\Omega }}_{\mathrm{FS}}$ is passive for the SSA solver, which means that it is not included in the global matrix assembly of ${\mathbf{A}}_{\mathrm{SSA}}$ and vice versa.

Two kinds of iterations are involved, since computing either ${\mathit{u}}_{\mathrm{FS},k}$ or ${\mathit{u}}_{\mathrm{SSA},k}$ for the $k$th coupled iteration also requires Picard iteration by the nonlinearity in the viscosity. As the experiments will show, calculating ${\mathit{u}}_{\mathrm{FS},k}$ dominates the computation time in the coupled model. The coupled model is therefore more efficient if the total number of FS Picard iterations (the sum of FS Picard iterations over all coupled iterations) decreases. This is accomplished by limiting the number of FS Picard iterations before continuing to compute ${\mathit{u}}_{\mathrm{SSA},k}$, instead of continuing until the convergence tolerance ${\mathit{\epsilon }}_{\mathrm{P}}$ is reached, since it is inefficient to solve very accurately for ${\mathit{u}}_{\mathrm{FS},k}$ if the boundary condition at ${\mathit{x}}_{\mathrm{c}}$ is not yet accurate. Despite interrupting the Picard iteration, the final solution includes a converged FS solution since the coupled tolerance ${\mathit{\epsilon }}_{\mathrm{c}}$ is reached. Picard iteration for ${\mathit{u}}_{\mathrm{SSA},k}$ is always continued until convergence since the computation time is negligible compared to FS.

An element may switch from ${\mathrm{\Omega }}_{\mathrm{SSA}}$ to ${\mathrm{\Omega }}_{\mathrm{FS}}$, for example during grounding line advance. Then, the coupled iteration either starts with the initial condition for ${\mathit{u}}_{\mathrm{FS}}$ if the element is in ${\mathrm{\Omega }}_{\mathrm{FS}}$ for the first time, or the latest ${\mathit{u}}_{\mathrm{FS}}\left(t\right)$ computed in this element, before switching to SSA.

Numerical experiments

To validate the coupled model, we first verify for a conceptual ice shelf ramp that solutions obtained with the coupled model resemble the FS velocity in 2-D and 3-D. Then the coupled model is applied to a 2-D conceptual marine ice sheet (MIS). Whenever “accuracy of the coupled model” is mentioned, this refers to the accuracy of the coupled model compared to the FS model. Investigating the accuracy of the FS model itself is outside the scope of this study. No convergence study of the FS model with respect to discretization in either time or space is performed. Instead, equivalent settings are used for the FS and coupled model, such that they can be compared, and the FS model is regarded as a reference solution.

Ice shelf ramp

Two-dimensional ice shelf ramp

A simplified test case is chosen for which the analytical solution to the SSA equations exists in 2-D as described in . It consists of a 200 km long ice shelf (see Fig. ), with a horizontal inflow velocity $u(\mathrm{0},z)=\mathrm{100}$ m yr${}^{-\mathrm{1}}$ and a calving front at $x=\mathrm{200}$ km where the hydrostatic pressure as exerted by the sea water is applied. The shelf thickness linearly decreases from 400 m at $x=\mathrm{0}$ to 200 m at $x=\mathrm{200}$ km; $z_{\mathrm{b}}$ and $z_{\mathrm{s}}$ follow from the flotation criterion (Eq. ). By construction, the SSA model is expected to be a good approximation of the FS model. The domain is discretized by a structured mesh and equidistant nodes on the horizontal axis and extruded along the vertical to quadrilaterals. All constants used and all mesh characteristics are specified in Table .

Three models are applied to this setup, FS-only, SSA-only, and the coupled model, for which the horizontal velocities are denoted as $u_{\mathrm{FS}}$, $u_{\mathrm{SSA}}$, and $u_{\mathrm{c}}$, respectively. The relative node-wise velocity differences between $u_{\mathrm{SSA}}$ and $u_{\mathrm{FS}}$ stay below 0.02 % in the entire domain. However, computing time for the SSA solution only takes 3 % of that of the FS solution, which is promising for the potential speedup of the coupled model.

The horizontal velocity $u_{\mathrm{c}}$ (m yr${}^{-\mathrm{1}}$) and node-wise difference $|u_{\mathrm{FS}}-u_{\mathrm{c}}|/u_{\mathrm{FS}}\cdot \mathrm{100}$ (%) in the coupled solution for the 2-D ice shelf ramp. The vertical scale is exaggerated 100 times. The ice thickness ranges from 400 to 200 m.

[Figure omitted. See PDF]

The coupling location is fixed at $x_{\mathrm{c}}=\mathrm{100}$ km, as no grounding line is present to relate $x_{\mathrm{c}}$ to. In the first coupled iteration, $u_{\mathrm{c}}(x_{\mathrm{c}},z_{\mathrm{b}})=\mathrm{100}$ m yr${}^{-\mathrm{1}}$, while in the final solution $u_{\mathrm{FS}}(x_{\mathrm{c}},z_{\mathrm{b}})=\mathrm{4805}$ m yr${}^{-\mathrm{1}}$. The cryostatic pressure applied to ${\mathrm{\Omega }}_{\mathrm{FS}}$ at $x_{\mathrm{c}}$ buttresses the ice flow completely and the force imbalance of the hydrostatic pressure at the calving front does not yet influence the velocity $u_{\mathrm{c}}$ in ${\mathrm{\Omega }}_{\mathrm{FS}}$. In the second iteration, when ${\mathit{f}}_{\mathrm{SSA}}$ is applied, a maximum difference of only 0.3 % between $u_{\mathrm{FS}}$ and $u_{\mathrm{c}}$ in the entire domain remains. The coupling converges after three iterations; the velocity $u_{\mathrm{c}}$ and relative difference compared to FS are shown in Fig. . Convergence of the coupled model requires 31 FS Picard iterations compared to 35 for FS-only. However, assembly time per FS iteration almost doubles in the coupled model compared to the FS model, and assembly time dominates the computational work in this simplified 2-D case. Therefore, the coupled model needs almost twice as much computation time as the FS model. This issue is due to the use of passive elements and is addressed in the “Discussion” section (Sect. ).

Three-dimensional ice shelf ramp

The 2-D ice shelf ramp is extruded along the $y$ axis (see Fig. ). On both lateral boundaries at $y=\mathrm{0}$ and 20 km, $\mathit{u}\cdot \mathit{n}=\mathrm{0}$. All other boundary conditions remain identical to the 2-D case, and the coupling interface is located halfway ${\mathit{x}}_{\mathrm{c}}=(\mathrm{100},y)$ km. First the solutions of the FS and SSA model in Elmer/Ice will be compared before applying the coupled model.

The limited width of the domain (20 km) in combination with the boundary condition $\mathit{u}\cdot \mathit{n}=\mathrm{0}$ at both lateral sides yields a negligible flow in the $y$ direction ($v_{\mathrm{FS}}<{\mathrm{10}}^{-\mathrm{8}}$ m yr${}^{-\mathrm{1}}$). Despite differences in the models, the relative difference in $u$ is below 1.5 %. Running the experiment with the SSA model takes only 0.8 % of the time needed to run it with the FS model.

The maximum relative difference between ${\mathit{u}}_{\mathrm{FS}}$ and ${\mathit{u}}_{\mathrm{c}}$ is 1.4 %, which is of the same order of magnitude as the velocity difference between FS and SSA. The mean assembly time per FS iteration is 6 % higher than in the FS-only model, but the solution time decreases by 55 %. Convergence of the coupled model requires 30 FS iterations compared to 27 for FS-only. The total computation time decreases by 32 %.

Horizontal velocity $u_{\mathrm{c}}$ (m yr${}^{-\mathrm{1}}$) from coupled model for the 3-D ice shelf ramp with ${\mathit{x}}_{\mathrm{c}}=(\mathrm{100},y)$ km.

[Figure omitted. See PDF]

Marine ice sheet

First, a diagnostic MIS experiment is performed in 2-D to compare velocities for the initial geometry. After one time step, velocity differences between the coupled and FS models yield geometric differences. In prognostic experiments, velocity differences can therefore be due to the coupling and to the different geometry for which the velocity is solved. Computation times for the FS and coupled model are presented for the prognostic case only.

The coupled velocity $u_{\mathrm{c}}$ (m yr${}^{-\mathrm{1}}$) and relative difference $|u_{\mathrm{FS}}-u_{\mathrm{c}}|/u_{\mathrm{FS}}\cdot \mathrm{100}$ (%), for the diagnostic MIS experiment. The bedrock is shaded in grey, $x_{\mathrm{GL}}=\mathrm{730.8}$ km, and $x_{\mathrm{c}}=\mathrm{763.2}$ km (the mesh resolution yields $\left|\right|{\mathit{x}}_{\mathrm{c}}-{\mathit{x}}_{\mathrm{GL}}|{|}_{\mathrm{h}}=\mathrm{32.4}$ km). The vertical scale is exaggerated 100 times with an ice thickness ranging from 1435 to 296 m.

[Figure omitted. See PDF]

Diagnostic MIS experiment

The domain starts with an ice divide at $x=\mathrm{0}$, where $u=\mathrm{0}$, and terminates at a calving front at $x=L=\mathrm{1800}$ km. An equidistant grid with grid spacing $\mathrm{\Delta }x=\mathrm{3.6}$ km is used. Other values of constants and mesh characteristics are specified in Table . showed that resolving grounding line dynamics with an FS model requires very high mesh resolution around the grounding line. However, showed that the friction law assumed in this study (see Sect. ) reduces the mesh sensitivity of the FS model compared to the Weertman friction law assumed in , allowing the coarse mesh used here. The bedrock (m) is negative below sea level and is given by $$b\left(x\right)=\mathrm{200}-\mathrm{900}{\displaystyle \frac{x}{L}}.$$ Basal melt is neglected, and the surface accumulation $a_{\mathrm{s}}$ (m yr${}^{-\mathrm{1}}$) is a function of the distance from the ice divide: $$a_{\mathrm{s}}\left(x\right)={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{w}}}{\mathit{\rho }}}{\displaystyle \frac{x}{L}}.$$ This experimental setup is almost equivalent to , except that they applied a buttressing force to the FS equations. It is possible to parameterize buttressing for the SSA equations as well through applying a sliding coefficient . This was not done here as it may introduce a difference between the FS and SSA models that is unrelated to the coupling.

The diagnostic experiments are run on a steady-state geometry computed by the FS model. First, the experiment “SPIN” in is performed, starting from a uniform slab of ice ($H=\mathrm{300}$ m), applying the accumulation in Eq. () for 40 kyr, such that a steady state is reached. The geometry yielded from these SPIN runs (which include buttressing) is used in simulations without buttressing until a new steady state (defined as a relative ice volume change below ${\mathrm{10}}^{-\mathrm{5}}$) is reached. This removal of buttressing leads to grounding line retreat from 871.2 to 730.8 km (Fig. ).

Again, FS-only, SSA-only, and the coupled model are applied to this setup. Where $u_{\mathrm{FS}}\ge \mathrm{5}$ m yr${}^{-\mathrm{1}}$, the relative difference between $u_{\mathrm{FS}}$ and $u_{\mathrm{SSA}}$ is below 1.8 %. The velocity $u_{\mathrm{c}}$ is given in Fig. , with $d_{\mathrm{GL}}=\mathrm{30}$ km such that 58 % of the nodes in the horizontal footprint mesh are located inside ${\mathrm{\Omega }}_{\mathrm{SSA}}$ ($\mathit{\theta }=\mathrm{0.58}$). The coupled model converges after 27 FS iterations on the restricted domain ${\mathrm{\Omega }}_{\mathrm{FS}}$, compared to 24 Picard iterations in the FS model. The relative difference between $u_{\mathrm{FS}}$ and $u_{\mathrm{c}}$ is below 0.5 % (Fig. ); this small difference shows that $d_{\mathrm{GL}}=\mathrm{30}$ km is sufficient. For this configuration, 4 % of the FS nodes are located between $x_{\mathrm{GL}}$ and $x_{\mathrm{c}}$ with $d_{\mathrm{GL}}=\mathrm{30}$ km; hence, decreasing $d_{\mathrm{GL}}$ does not affect the proportion of nodes in ${\mathrm{\Omega }}_{\mathrm{FS}}$ significantly. Therefore, $d_{\mathrm{GL}}$ is kept equal to 30 km for the prognostic experiment.

Computation times for the MIS simulation of 3000 years with FS-only and the coupled model. Model C$i$ denotes the coupled model, with $i$ being the maximum number of nonlinear FS iterations per coupled iteration; C5–C9 are omitted for brevity. The assembly time for ${\mathbf{A}}_{\mathrm{FS}}$ is denoted as $t_{\mathrm{A}}$. All relative computation times are given in percentage of the total time $t_{\mathrm{tot}}$. The number of FS and coupled iterations are averaged over the time steps.

Model	$t_{\mathrm{A}}$ (%)	$t_{\mathrm{s}}$ (%)	No. FS iter.	$t_{\mathrm{o}}$ (%)	$t_{\mathrm{SSA}}$ (%)	No. coupled iter.	$t_{\mathrm{tot}}$ (cpu s)
FS	75	5	5.0	20	–	–	48 641
C10	68	4	4.6	25	4	2.7	49 724
C4	61	3	3.7	31	5	2.9	44 143
C3	59	3	3.6	33	5	3.1	44 040
C2	56	3	3.4	36	5	3.4	44 334
C1	49	3	3.2	43	6	4.2	47 135

Absolute difference $|u_{\mathrm{FS}}-u_{\mathrm{c}}|$ (m yr${}^{-\mathrm{1}}$) after 3000 years. The vertical scale is exaggerated 100 times. The ice thickness ranges from 1445 to 296 m.

[Figure omitted. See PDF]

Prognostic MIS experiment

The prognostic experiment aims to verify model reversibility as in . Starting from the steady-state geometry, the ice temperature $T$ is lowered over a period of 500 years from $-\mathrm{10}$ to $-\mathrm{30}$ ${}^{\circ }$C and back according to $$\begin{array}{cc}{\displaystyle }& {\displaystyle }T\left(t\right)=-\mathrm{10}(\mathrm{2}-\cos (\mathrm{2}\mathit{\pi }t/\mathrm{500}\left)\right){}^{\circ }\text{C}\\ & {\displaystyle }& {\displaystyle }\text{for }\mathrm{0}\le t\le \mathrm{500}\text{years.}\end{array}$$ The resulting change in $\mathcal{A}$ (see Eq. ) induces a grounding line advance and retreat and changes ${\mathrm{\Omega }}_{\mathrm{SSA}}$ by Eq. (). Afterwards, $T=-\mathrm{10}$ ${}^{\circ }$C for 2500 years. Mass balance forcing is kept constant throughout. The length of one time step is 1 year.

The maximum difference between $u_{\mathrm{c}}$ and $u_{\mathrm{FS}}$ after 3000 years is 10 m yr${}^{-\mathrm{1}}$ (shown in Fig. ), corresponding to a relative difference of 1.6 %. The time evolution of $x_{\mathrm{GL}}$, $u_{\mathrm{b}}\left(x_{\mathrm{GL}}\right)$, $H\left(x_{\mathrm{GL}}\right)$, and the grounded volume $V_{\mathrm{g}}$ is shown in Figs.  and . In general, $u_{\mathrm{b}}$ is slightly higher in the coupled model, with a maximum difference of 5.3 % in the entire experiment. The grounding line advances to $x_{\mathrm{GL}}=\mathrm{1036.8}$ km in the FS model and $x_{\mathrm{GL}}=\mathrm{1044}$ km in the coupled model. The FS model returns back to the original $x_{\mathrm{GL}}=\mathrm{730.8}$ km, but the coupled model yields $x_{\mathrm{GL}}=\mathrm{734.4}$ km, an offset of one grid point. The maximum difference in thickness is 1 %. After 3000 years, $V_{\mathrm{g}}$ still decreases but the relative difference is below ${\mathrm{10}}^{-\mathrm{5}}$ between two time steps.

Time evolution of $x_{\mathrm{GL}}$ (red) and $u_{\mathrm{b}}\left(x_{\mathrm{GL}}\right)$ (blue) with solid lines for FS and dashed lines for the coupled model.

[Figure omitted. See PDF]

To investigate the efficiency of the coupled model, the simulation is performed with 10 different settings, where the maximum number of FS iterations per coupled iteration is varied from 1 to 10. The assembly of the FS matrix takes 75 % of the computation time of the FS model (see $t_{\mathrm{A}}$ in Table ), and assembly time per FS iteration is similar for the coupled and FS model. Only 5 % of the computation time is used to solve the linearized FS system ($t_{\mathrm{s}}$ in Table ). For all coupled simulations, assembling and solving the SSA matrix ($t_{\mathrm{SSA}}$) takes 4 %–6 %. All the time that is left will be called overhead, $t_{\mathrm{o}}$, which includes launching solvers, i.e., allocating memory space for vectors and matrices, the surface evolution, and solvers for post-processing. As expected, the total number of FS iterations is the smallest when just performing one FS Picard iteration per coupled iteration. However, the model then changes between solvers more often, meaning that both overheads and the time to solve the SSA model increase. It turns out that a limit of three FS Picard iterations per coupled iteration balances minimizing $t_{\mathrm{o}}$ and $t_{\mathrm{A}}$, yielding a 10 % decrease in computation time with respect to the FS model. This speedup comes from a lower number of FS Picard iterations (Table ) and a slight decrease in the time used to solve the linearized FS system (13 % lower than the time that the FS model takes).

Time evolution of $H_{\mathrm{GL}}=H\left(x_{\mathrm{GL}}\right)$ and grounded volume $V_{\mathrm{g}}$ with solid lines for FS and dashed lines for the coupled model.

[Figure omitted. See PDF]

Discussion

The presented coupling is dynamic, since the coupling interface ${\mathit{x}}_{\mathrm{c}}$ changes with grounding line changes, but the distance $d_{\mathrm{GL}}$ that defines ${\mathit{x}}_{\mathrm{c}}$ has to be chosen such that the FS velocity at the interface is almost independent of $z$. In the experiments described in Sect. 4, this is already the case at the grounding line. We propose that further studies let ${\mathrm{\Omega }}_{\mathrm{SSA}}$ be determined automatically, for example, by a tolerance for the vertical variation in the horizontal velocities, which should be close to zero in order to allow for a smooth coupling to SSA. Another option is to use a posteriori error estimates based on the residual .

The current implementation in Elmer/Ice does not give as much speedup as expected from computation times of the FS- and SSA-only models for the ice shelf ramp ($t_{\mathrm{SSA}}=\mathrm{0.03}{t}_{\mathrm{FS}}$) and from the performance estimates in Sect. 4.2. This is due to an inefficient matrix assembly. The assembly of the system matrix ${\mathbf{A}}_{\mathrm{FS}}$ restricted to ${\mathrm{\Omega }}_{\mathrm{FS}}$ currently takes at least as much time as the assembly for the full domain $\mathrm{\Omega }$, even though the domain ${\mathrm{\Omega }}_{\mathrm{FS}}$ is much smaller than $\mathrm{\Omega }$; in Eq. (13), $\mathit{\theta }=\mathrm{0.5}$ for the ice shelf ramp and $\mathit{\theta }=\mathrm{0.58}$ for the diagnostic MIS experiment. Since the assembly time dominates the total solution time in simple 2-D simulations, this is problematic. The inefficient assembly is caused by the use of passive elements implemented in the overarching Elmer code , which allow the de- and reactivation of elements. A passive element is not included in the global matrix assembly, but every element must be checked to determine if it is passive. The inefficient assembly can be overcome by implementing the coupling on a lower level, hard-coded inside the FS solver. This was done for the coupling of SIA and FS in ISCAL , which showed significant speedup when restricting the FS solver to a smaller domain. However, using passive elements is more flexible, since the coupling is independent of the solver used to compute velocities outside ${\mathrm{\Omega }}_{\mathrm{SSA}}$. One is free to choose between the two different FS solvers in Elmer/Ice see or to apply ISCAL. The latter is irrelevant in the experiments presented here since both the grounded and floating ice experiences low basal drag, and SIA is not capable of representing ice stream and shelf flow. Only a preliminary 3-D experiment is performed here, since the current implementation is not sufficiently efficient to allow extensive testing in 3-D. If the coupling is implemented efficiently such that the time spent on solving the FS equations on the restricted domain ${\mathrm{\Omega }}_{\mathrm{FS}}$ scales with the size of ${\mathrm{\Omega }}_{\mathrm{FS}}$, the computational work is expected to decrease significantly (see Sect. 4.2).

Conclusions

We have presented a novel FS–SSA coupling in Elmer/Ice, showing a large potential for reducing the computation time without losing accuracy. At the coupling interface, the FS velocity is applied as an inflow boundary condition to SSA. Together with the cryostatic pressure, a depth-averaged contact force resulting from the SSA velocity is applied as a boundary condition for FS. The main finding of this study is that the two-way coupling is stable and converges to a velocity that is very similar to the FS model in the tests on conceptual marine ice sheets, and it yields a speedup in 3-D.

In diagnostic runs, the relative difference in velocity obtained from the coupled model and the FS model is below 1.5 % when applying SSA at least 30 km seaward from the grounding line. During a transient simulation, where the coupling interface changes dynamically with the migration of the grounding line, the coupled model is very similar to the FS model, with a maximum difference of 5.3 % in basal velocity at the grounding line. An offset of 3.6 km remains in the reversibility experiment in Sect. 4.3, which is within the range of the expected resolution dependence for FS models .

In experiments involving areas where SIA is applicable, this new FS–SSA model can be combined with the ISCAL method in that couples SIA and FS in Elmer/Ice. This mixed model is motivated by paleo-simulations, but reducing computational work by the combination of multiple approximation levels is also convenient for parameter studies, ensemble simulations, and inverse problems.

The code of Elmer/Ice is available at https://github.com/ElmerCSC/elmerfem/tree/elmerice (last access: 13 November 2018). An example of the coupling is provided at 10.5281/zenodo.1202407 . The version of the Elmer/Ice code that includes the coupling discussed in the paper can be accessed by using the hash qualifier linked to the commit of the coupling code at https://github.com/ElmerCSC/elmerfem/archive/ba117583defafe98bb6fd1793c9c6f341c0c876.zip (last access: 13 November 2018).

Derivation of the interface boundary condition

The boundary condition in Sect.  between the FS and the SSA domains is derived following a standard procedure in FEM using the weak formulation of the equations. Let ${\mathrm{\Omega }}_{\mathrm{FS}}\in {\mathbb{R}}^d$, $d=\mathrm{2},\mathrm{3}$ denote the open FS domain in two or three dimensions with the boundary ${\mathrm{\Gamma }}_{\mathrm{FS}}$. After multiplying Eq. () with a test function $\mathit{v}$ and integrating over the domain ${\mathrm{\Omega }}_{\mathrm{FS}}$, the weak form of Eq. () is $$-\underset{{\mathrm{\Omega }}_{\mathrm{FS}}}{\int }\mathit{v}\cdot \left(\mathrm{\nabla }\cdot \mathit{\sigma }\right)=\underset{{\mathrm{\Omega }}_{\mathrm{FS}}}{\int }\mathit{\rho }\mathit{v}\cdot \mathit{g}.$$ Use the definition of $\mathit{\sigma }$ and the divergence theorem to rewrite Eq. (): $$\underset{{\mathrm{\Omega }}_{\mathrm{FS}}}{\int }\mathit{\eta }\mathbf{D}\left(\mathit{u}\right):\mathbf{D}\left(\mathit{v}\right)-\underset{{\mathrm{\Omega }}_{\mathrm{FS}}}{\int }p\mathrm{\nabla }\cdot \mathit{v}=\underset{{\mathrm{\Omega }}_{\mathrm{FS}}}{\int }\mathit{\rho }\mathit{v}\cdot \mathit{g}+\underset{{\mathrm{\Gamma }}_{\mathrm{FS}}}{\int }\mathit{v}\cdot \mathit{\sigma }\cdot \mathit{n}.$$ The operation $\mathbf{A}:\mathbf{B}$ denotes the sum ${\sum }_{i,j}{A}_{ij}{B}_{ij}$. The test function $\mathit{v}$ vanishes on the inflow boundary ${\mathrm{\Gamma }}_{\mathrm{i}}$, has a vanishing normal component on the bedrock boundary ${\mathrm{\Gamma }}_{\mathrm{b}}$, and lives in the Sobolev space $\left[W^{\mathrm{1},\mathrm{1}/n+\mathrm{1}}\right({\mathrm{\Omega }}_{\mathrm{FS}}){]}^d$ , i.e., $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathit{v}\in {\mathcal{V}}_{\mathrm{0}}=\mathit{\{}\mathit{v}\in [W^{\mathrm{1},\mathrm{1}/n+\mathrm{1}}\left({\mathrm{\Omega }}_{\mathrm{FS}}\right){]}^d|\mathit{v}{|}_{{\mathrm{\Gamma }}_{\mathrm{i}}}=\mathrm{0},\\ & {\displaystyle }& {\displaystyle }\mathit{v}{|}_{{\mathrm{\Gamma }}_{\mathrm{b}}}\cdot \mathit{n}=\mathrm{0}\mathit{\}}.\end{array}$$ The space ${\mathcal{V}}_{\mathrm{0}}$ has this form because the boundary conditions on ${\mathrm{\Gamma }}_{\mathrm{i}}$ and ${\mathrm{\Gamma }}_{\mathrm{b}}$ are of Dirichlet type. Furthermore, there is a lateral boundary ${\mathrm{\Gamma }}_{\mathrm{\ell }}$ for ${\mathrm{\Omega }}_{\mathrm{FS}}\in {\mathbb{R}}^{\mathrm{3}}$, where the normal component also vanishes ($\mathit{v}{|}_{{\mathrm{\Gamma }}_{\mathrm{\ell }}}\cdot \mathit{n}=\mathrm{0}$), and we assume a vanishing Cauchy-stress vector for unset boundary conditions to velocity components, such that the integral over ${\mathrm{\Gamma }}_{\mathrm{\ell }}$ vanishes. Then, the boundary integral in Eq. () consists of a sum of the remaining boundary terms: $$\underset{{\mathrm{\Gamma }}_{\mathrm{FS}}}{\int }\mathit{v}\cdot \mathit{\sigma }\cdot \mathit{n}=\sum _{i=\mathrm{1}}^{d-\mathrm{1}}\underset{{\mathrm{\Gamma }}_{\mathrm{b}}}{\int }f\mathit{u}\cdot {\mathit{t}}_{\mathrm{i}}\mathit{v}\cdot {\mathit{t}}_{\mathrm{i}}-\underset{{\mathrm{\Gamma }}_{\mathrm{w}}}{\int }{p}_{\mathrm{w}}\mathit{n}\cdot \mathit{v}+\underset{{\mathrm{\Gamma }}_{\mathrm{FSint}}}{\int }\mathit{v}\cdot \mathit{\sigma }\cdot \mathit{n},$$ given by the boundary conditions on ${\mathrm{\Gamma }}_{\mathrm{b}}$ in Eqs. () and (), on the ocean boundary ${\mathrm{\Gamma }}_{\mathrm{w}}$ in Eq. (), and the internal boundary ${\mathrm{\Gamma }}_{\mathrm{FSint}}$ between the FS and the SSA domains. The force $\mathit{\sigma }\cdot \mathit{n}$ on ${\mathrm{\Gamma }}_{\mathrm{FSint}}$ is determined by the SSA solution.

The open SSA domain ${\mathrm{\Omega }}_{\mathrm{SSA}}\in {\mathbb{R}}^{\mathrm{2}}$, coupled to ${\mathrm{\Omega }}_{\mathrm{FS}}\in {\mathbb{R}}^{\mathrm{3}}$, has the boundary ${\mathrm{\Gamma }}_{\mathrm{SSA}}={\mathrm{\Gamma }}_{\mathrm{SSAint}}\cup {\mathrm{\Gamma }}_{\mathrm{CF}}\cup {\mathrm{\Gamma }}_{\mathrm{\ell }}$ where ${\mathrm{\Gamma }}_{\mathrm{SSAint}}$ is adjacent to ${\mathrm{\Omega }}_{\mathrm{FS}}$ and partly coinciding with ${\mathrm{\Gamma }}_{\mathrm{FSint}}$ (but of one dimension less) and ${\mathrm{\Gamma }}_{\mathrm{CF}}$ is at the calving front. Let $\mathbf{B}$ have the elements $$\begin{array}{cc}{\displaystyle }& {\displaystyle }{B}_{\mathrm{11}}=\mathrm{4}\stackrel{\mathrm{‾}}{\mathit{\eta }}{\displaystyle \frac{\partial u}{\partial x}}+\mathrm{2}\stackrel{\mathrm{‾}}{\mathit{\eta }}{\displaystyle \frac{\partial v}{\partial y}},B_{\mathrm{12}}=B_{\mathrm{21}}=\stackrel{\mathrm{‾}}{\mathit{\eta }}{\displaystyle \frac{\partial u}{\partial y}}+\stackrel{\mathrm{‾}}{\mathit{\eta }}{\displaystyle \frac{\partial v}{\partial x}},\\ & {\displaystyle }& {\displaystyle }{B}_{\mathrm{22}}=\mathrm{2}\stackrel{\mathrm{‾}}{\mathit{\eta }}{\displaystyle \frac{\partial u}{\partial x}}+\mathrm{4}\stackrel{\mathrm{‾}}{\mathit{\eta }}{\displaystyle \frac{\partial v}{\partial y}},\end{array}$$ when $d=\mathrm{3}$. If $d=\mathrm{2}$, then $B=\mathrm{4}\stackrel{\mathrm{‾}}{\mathit{\eta }}\partial u/\partial x$. Then the SSA equations Eq. () can be written as follows: $${\mathrm{\nabla }}_{\mathrm{h}}\cdot \mathbf{B}={\mathit{f}}_{\mathrm{g}},$$ where ${\mathit{f}}_{\mathrm{g}}=\mathit{\rho }gH{\mathrm{\nabla }}_{\mathrm{h}}{z}_{\mathrm{s}}$ and ${\mathrm{\nabla }}_{\mathrm{h}}$ is the horizontal gradient operator. The boundary condition on ${\mathrm{\Gamma }}_{\mathrm{SSAint}}$ is the Dirichlet condition (Eq. ), and the force due to the water pressure at the calving front ${\mathrm{\Gamma }}_{\mathrm{CF}}$ is ${\mathit{f}}_{\mathrm{CF}}$, as in Eq. () but integrated over $z$. Define the two test spaces $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\mathcal{W}=\mathit{\{}\mathit{v}\in [W^{\mathrm{1},\mathrm{1}/n+\mathrm{1}}\left({\mathrm{\Omega }}_{\mathrm{SSA}}\right){]}^{d-\mathrm{1}}|\mathit{v}{|}_{{\mathrm{\Gamma }}_{\mathrm{\ell }}}\cdot \mathit{n}=\mathrm{0}\mathit{\}},\\ & {\displaystyle }& {\displaystyle }{\mathcal{W}}_{\mathrm{0}}=\mathit{\{}\mathit{v}\in \mathcal{W}|\mathit{v}{|}_{{\mathrm{\Gamma }}_{\mathrm{SSAint}}}=\mathrm{0}\mathit{\}}.\end{array}$$ Multiply Eq. () by $\mathit{v}\in {\mathcal{W}}_{\mathrm{0}}$ and integrate. The weak form of Eq. () is $$\underset{{\mathrm{\Omega }}_{\mathrm{SSA}}}{\int }\mathit{v}\cdot ({\mathrm{\nabla }}_{\mathrm{h}}\cdot \mathbf{B})=\underset{{\mathrm{\Omega }}_{\mathrm{SSA}}}{\int }\mathit{v}\cdot {\mathit{f}}_{\mathrm{g}}.$$ Apply the divergence theorem to Eq. () to obtain $$\begin{array}{cc}{\displaystyle }& {\displaystyle }-\underset{{\mathrm{\Omega }}_{\mathrm{SSA}}}{\int }{\mathrm{\nabla }}_{\mathrm{h}}\mathit{v}:\mathbf{B}+\underset{{\mathrm{\Gamma }}_{\mathrm{SSA}}}{\int }\mathit{v}\cdot \mathbf{B}\cdot \mathit{n}=-\underset{{\mathrm{\Omega }}_{\mathrm{SSA}}}{\int }{\mathrm{\nabla }}_{\mathrm{h}}\mathit{v}:\mathbf{B}+\\ & {\displaystyle }& {\displaystyle }+\underset{{\mathrm{\Gamma }}_{\mathrm{CF}}}{\int }\mathit{v}\cdot {\mathit{f}}_{\mathrm{CF}}+\underset{{\mathrm{\Gamma }}_{\mathrm{SSAint}}}{\int }\mathit{v}\cdot {\mathit{f}}_{\mathrm{SSA}}=\underset{{\mathrm{\Omega }}_{\mathrm{SSA}}}{\int }\mathit{v}\cdot {\mathit{f}}_{\mathrm{g}}.\end{array}$$ A mesh is constructed to cover ${\mathrm{\Omega }}_{\mathrm{FS}}$ and ${\mathrm{\Omega }}_{\mathrm{SSA}}$ with nodes at ${\mathit{x}}_{\mathrm{i}}$. In the finite element solution of Eq. (), the linear test function ${\mathit{v}}_{\mathrm{i}}\in {\mathcal{W}}_{\mathrm{0}}$ is non-zero at ${\mathit{x}}_{\mathrm{i}}$ and zero in all other nodes. The integral over ${\mathrm{\Gamma }}_{\mathrm{SSAint}}$ vanishes when $\mathit{v}\in {\mathcal{W}}_{\mathrm{0}}$. The finite element solution ${\mathit{u}}_{\mathrm{h}}$ of Eqs. () and () satisfies $$\begin{array}{cc}{\displaystyle }& {\displaystyle }-\underset{{\mathrm{\Omega }}_{\mathrm{SSA}}}{\int }{\mathrm{\nabla }}_{\mathrm{h}}{\mathit{v}}_{\mathrm{i}}:\mathbf{B}\left({\mathit{u}}_{\mathrm{h}}\right)+\underset{{\mathrm{\Gamma }}_{\mathrm{CF}}}{\int }{\mathit{v}}_{\mathrm{i}}\cdot {\mathit{f}}_{\mathrm{CF}}-\underset{{\mathrm{\Omega }}_{\mathrm{SSA}}}{\int }{\mathit{v}}_{\mathrm{i}}\cdot {\mathit{f}}_{\mathrm{g}}=\mathrm{0},\\ & {\displaystyle }& {\displaystyle }{\mathit{x}}_{\mathrm{i}}\in {\mathrm{\Omega }}_{\mathrm{SSA}}\cup {\mathrm{\Gamma }}_{\mathrm{CF}}.\end{array}$$ It follows from Eq. () that with a test function ${\mathit{v}}_{\mathrm{i}}\in \mathcal{W}$ that is non-zero on ${\mathrm{\Gamma }}_{\mathrm{SSAint}}$ and the solution ${\mathit{u}}_{\mathrm{h}}$ from Eq. () $$\begin{array}{}& {\displaystyle }& {\displaystyle }\underset{{\mathrm{\Gamma }}_{\mathrm{SSAint}}}{\int }{\mathit{v}}_{\mathrm{i}}\cdot {\mathit{f}}_{\mathrm{SSA}}=& {\displaystyle }& {\displaystyle }=\underset{{\mathrm{\Omega }}_{\mathrm{SSA}}}{\int }{\mathrm{\nabla }}_{\mathrm{h}}{\mathit{v}}_{\mathrm{i}}:\mathbf{B}\left({\mathit{u}}_{\mathrm{h}}\right)-\underset{{\mathrm{\Gamma }}_{\mathrm{CF}}}{\int }{\mathit{v}}_{\mathrm{i}}\cdot {\mathit{f}}_{\mathrm{CF}}+\underset{{\mathrm{\Omega }}_{\mathrm{SSA}}}{\int }{\mathit{v}}_{\mathrm{i}}\cdot {\mathit{f}}_{\mathrm{g}},& {\displaystyle }& {\displaystyle }{\mathit{x}}_{\mathrm{i}}\in {\mathrm{\Omega }}_{\mathrm{SSA}}\cup {\mathrm{\Gamma }}_{\mathrm{CF}}\cup {\mathrm{\Gamma }}_{\mathrm{SSAint}}.\end{array}$$ The first integral in Eq. () corresponds to $({\mathbf{A}}_{\mathrm{SSA}}{\mathit{u}}_{\mathrm{SSA}}{)}_{\mathrm{i}}$ in Sect.  and $b_{\mathrm{SSAi}}$ to the second and third integrals. By Eq. (), the right-hand side of Eq. () vanishes for all ${\mathit{x}}_{\mathrm{i}}$ in ${\mathrm{\Omega }}_{\mathrm{SSA}}$ and on ${\mathrm{\Gamma }}_{\mathrm{CF}}$, but for a node on the internal boundary, ${\mathit{x}}_{\mathrm{i}}\in {\mathrm{\Gamma }}_{\mathrm{SSAint}}$, the force ${\mathit{f}}_{\mathrm{SSA}}$ from the ice due to the state ${\mathit{u}}_{\mathrm{h}}$ in ${\mathrm{\Omega }}_{\mathrm{SSA}}$ is obtained. The internal pressure in the ice in ${\mathrm{\Omega }}_{\mathrm{SSA}}$ is assumed to be cryostatic as in Eq. (). The total force on ${\mathrm{\Gamma }}_{\mathrm{FSint}}$ consists of one component due to the state ${\mathit{u}}_{\mathrm{h}}$ at ${\mathrm{\Gamma }}_{\mathrm{SSAint}}$ and one due to the cryostatic pressure there. Let ${\mathrm{\Omega }}_{\mathrm{SSA}}^{\ast }$ denote the mesh on ${\mathrm{\Omega }}_{\mathrm{SSA}}$, which is extruded in the $z$ direction. The common boundary between ${\mathrm{\Omega }}_{\mathrm{FS}}$ and ${\mathrm{\Omega }}_{\mathrm{SSA}}^{\ast }$ is ${\mathrm{\Gamma }}_{\mathrm{FSint}}$, and let ${\mathit{f}}_{\mathrm{SSA}}^{\ast }$ be the stress force there, independent of $z$. Since ${\int }_{z_{\mathrm{b}}}^{z_{\mathrm{s}}}{\mathit{f}}_{\mathrm{SSA}}^{\ast }={\mathit{f}}_{\mathrm{SSA}}$ at ${\mathrm{\Gamma }}_{\mathrm{FSint}}$, we have ${\mathit{f}}_{\mathrm{SSA}}^{\ast }=H^{-\mathrm{1}}{\mathit{f}}_{\mathrm{SSA}}$. Let ${\mathit{v}}_{\mathrm{i}}$ be a test function on ${\mathrm{\Omega }}_{\mathrm{FS}}\bigcup {\mathrm{\Omega }}_{\mathrm{SSA}}^{\ast }$ which is non-zero on ${\mathrm{\Gamma }}_{\mathrm{FSint}}$ and zero in all other nodes. Then the weak form of the force balance at ${\mathrm{\Gamma }}_{\mathrm{FSint}}$ is $$\begin{array}{cc}{\displaystyle }& {\displaystyle }\underset{\mathrm{FSint}}{\int }{\mathit{v}}_{\mathrm{i}}\cdot \mathit{\sigma }\cdot \mathit{n}=\underset{\mathrm{FSint}}{\int }{\mathit{f}}_{\mathrm{SSA}}^{\ast }\cdot {\mathit{v}}_{\mathrm{i}}-\underset{\mathrm{FSint}}{\int }\mathit{\rho }g(z_{\mathrm{s}}-z)\mathit{n}\cdot {\mathit{v}}_{\mathrm{i}}\\ & {\displaystyle }& {\displaystyle }=\underset{\mathrm{FSint}}{\int }{H}^{-\mathrm{1}}{\mathit{f}}_{\mathrm{SSA}}\cdot {\mathit{v}}_{\mathrm{i}}-\underset{\mathrm{FSint}}{\int }\mathit{\rho }g(z_{\mathrm{s}}-z)\mathit{n}\cdot {\mathit{v}}_{\mathrm{i}},\end{array}$$ and the corresponding strong form of the boundary condition at ${\mathrm{\Gamma }}_{\mathrm{FSint}}$ is $$\mathit{\sigma }\cdot \mathit{n}=H^{-\mathrm{1}}{\mathit{f}}_{\mathrm{SSA}}-\mathit{\rho }g(z_{\mathrm{s}}-z)\mathit{n};$$ cf. Eq. (). Thus, by computing the residual as in Eq. (), the two finite element solutions in ${\mathrm{\Omega }}_{\mathrm{FS}}$ and ${\mathrm{\Omega }}_{\mathrm{SSA}}$ are coupled together at the common boundary ${\mathrm{\Gamma }}_{\mathrm{FSint}}$ and ${\mathrm{\Gamma }}_{\mathrm{SSAint}}$.

Numerical values of the constants used in the ice shelf ramp experiment. Since the shelf is afloat, there is no sliding at the base.

Parameter	Symbol	Value	Unit
Ice density	$\mathit{\rho }$	900	kg m${}^{-\mathrm{3}}$
Water density	${\mathit{\rho }}_{\mathrm{w}}$	1000	kg m${}^{-\mathrm{3}}$
Gravitational acceleration	$g$	9.81	m s${}^{-\mathrm{2}}$
Fluidity parameter	$\mathcal{A}$	${\mathrm{10}}^{-\mathrm{16}}$	Pa${}^{-\mathrm{3}}$ yr${}^{-\mathrm{1}}$
Number elements	$N_z$	10
	$N_x$	120
	$N_y$	10
Picard convergence tolerance	${\mathit{\epsilon }}_{\mathrm{P}}$	${\mathrm{10}}^{-\mathrm{3}}$
Coupled convergence tolerance	${\mathit{\epsilon }}_{\mathrm{c}}$	${\mathrm{10}}^{-\mathrm{4}}$

Numerical values of the constants for the MIS experiment.

Parameter	Symbol	Value	Unit
Ice density	$\mathit{\rho }$	910	kg m${}^{-\mathrm{3}}$
Water density	${\mathit{\rho }}_{\mathrm{w}}$	1000	kg m${}^{-\mathrm{3}}$
Gravitational acceleration	$g$	9.81	m s${}^{-\mathrm{2}}$
Sliding parameter	$\mathit{\beta }$	$\mathrm{7}\times {\mathrm{10}}^{-\mathrm{6}}$	MPa m${}^{-\mathrm{4}/\mathrm{3}}$ yr${}^{\mathrm{1}/\mathrm{3}}$
Temperature	$T$	$-\mathrm{10}$	${}^{\circ }$C
Number elements	$N_z$	11
	$N_x$	500
Picard convergence tolerance	${\mathit{\epsilon }}_{\mathrm{P}}$	${\mathrm{10}}^{-\mathrm{4}}$
Coupled convergence tolerance	${\mathit{\epsilon }}_{\mathrm{c}}$	${\mathrm{10}}^{-\mathrm{4}}$
Time step	d$t$	1	year

NK, ECHvD, RSWvdW, MBvG, PL, and LvS designed the study. ECHvD implemented the coupling and carried out the numerical simulations, with support from TZ and GC. ECHvD drafted the paper with support from NK, and all authors contributed to the final version.

The authors declare that they have no conflict of interest.

Acknowledgements

This work has been supported by FORMAS grant 214-2013-1600 to Nina Kirchner. Thomas Zwinger's contribution was supported by the Academy of Finland (grant number 286587). The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the PDC Center for High Performance Computing at KTH. We are grateful to Mika Malinen, Peter Råback, and Juha Ruokolainen for advice in developing the coupling, to Rupert Gladstone for providing the setup as in , and to Felicity Holmes, Guillaume Jouvet, and Daniel Farinotti for their feedback on a draft of the paper. We wish to acknowledge the constructive comments of two anonymous reviewers, which contributed to improving the paper. Edited by: Didier Roche Reviewed by: two anonymous referees