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

. Ice ﬂow forced by gravity is governed by the full Stokes (FS) equations, which are computationally expensive to solve due to the nonlinearity introduced by the rheol-ogy. 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 10 3 years or longer. The shallow ice approximation (SIA) and shallow shelf approximation (SSA) are commonly used but are ac-curate 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

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 10 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(IPCC, 2014. 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 (Ritz et al., 2015). 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 (Alley et al., 2005;Phillips et al., 2010;Stokes et al., 2015). However, the computation time and memory required for an FS model to be applied to ice sheets restricts simulations to sub-millennial timescales Gladstone et al., 2012a;Nowicki et al., 2013;Seddik et al., 2012Seddik et al., , 2017Joughin et al., 2014). Therefore, approximations of the FS equations are employed for simulations over long timescales, such as the shallow ice approximation (SIA; Hutter, 1983), the shallow shelf approximation (SSA; Morland, 1987;MacAyeal, 1989), Blatter-Pattyn (Pattyn, 2003), and hybrid models (Hindmarsh, 2004;Bernales et al., 2017). 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 Published by Copernicus Publications on behalf of the European Geosciences Union.  Pattyn et al., 2012). 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 (Feldmann et al., 2014). It should be noted that the experiments in MISMIP3d consisted of laterally extruded idealized 2-D geometries with quite small sideward disturbances, and MISMIP+ (Asay-Davis et al., 2016) 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 (Gladstone et al., 2018).
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, Seroussi et al. (2012) coupled FS and SSA, in the framework of the Ice Sheet System Model (ISSM; Larour et al., 2012). 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 (Ahlkrona et al., 2016) 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 continentalscale simulations of marine ice sheets. However, the coupling is targeted at conducting paleo-simulations, for which much larger ice shelves have been present (Jakobsson et al., 2016;Nilsson et al., 2017). 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. 2, together with the boundary conditions. Memory and performance estimates of an FS-SSA coupling, independent of the specific coupling implemented, are provided in Sect. 2.3. Section 3 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. 4. The simulation of a 3000-year long cycle of grounding line advance and retreat (described in Sect. 4.2.2) 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 divergencefree: where 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 (Hutter, 1982). Therefore, the conservation of linear momentum under the action of gravity g can be described by where ∇ is the gradient operator, p pressure, η viscosity, ρ ice density, and g gravity. Letting σ denote the stress tensor, pressure p is the mean normal stress (p = −1/3 i σ ii ) and D(u) is the strain rate tensor, related by where I is the identity tensor. Together, Eqs. (1) and (2) are called the full Stokes equations. Observations by Glen (1952) suggest that the viscosity depends on temperature T and the effective strain rate D(u):  (Paterson, 1994). 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. (4), 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. 2.1.

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 (Morland, 1987;MacAyeal, 1989). 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. 2.2). Then, the conservation of linear momentum, Eq.
(2), simplifies to where the subscript h represents the components in the x-y plane, η the vertically integrated viscosity, H the thickness of the ice shelf, and z s the upper ice surface; see where w is eliminated using incompressibility; Eq. (1). The SSA equations are still nonlinear through η, 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.

Boundary conditions and time evolution
The coupling is applied to a marine ice sheet, with bedrock lying (partly) below sea level (see Fig. 1), 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 (u, N ), that relates the basal velocity u b and effective pressure N to the basal shear stress as where t i are the vectors spanning the tangential plane, n is the normal to the bed, and a b describes basal refreezing or melt. A sliding law suggested by Budd et al. (1979) is assumed, which depends on u b and the height above buoyancy z * such that Here, the sliding parameter β is constant in time and space. In line with Gladstone et al. (2017), instead of modeling N, a hydrostatic balance is assumed to approximate z * , implying a subglacial hydrology system entirely in contact with the ocean: where z b is the lower ice surface and ρ w the water density and the sea level is at z = 0. Equation (11) implies that z * equals zero when the flotation criterion (Archimedes' principle) is satisfied, i.e., where

Ice-ocean interface
As soon as the seawater pressure p w at the ice base z 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 Durand et al. (2009). At the ice-ocean interface, the tangential friction is neglected (f (u, N ) ≡ 0 in Eq. 8) and and σ · n = 0 above sea level (z > 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, σ · n = 0) and ice base at z s and z b behave as free surfaces according to where a s/b is the accumulation (a s/b > 0) or ablation (a s/b < 0) in meter ice equivalent per year, at the surface or base, respectively. By vertical integration of the incompressibility condition, Eq. (1), w can be eliminated using Leibniz integration rule and substituting the free surface equations (Eq. 14), which yields the thickness advection equation: where u and 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 u 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., Helanow and Ahlkrona, 2018) and there is no need for specialized iterative solution techniques to solve the so-called saddle-point problem that the FS equations pose (see Benzi et al., 2005). Suppose that the computational domain is discretized with N z nodes regularly placed in the z direction and N h nodes in a horizontal footprint mesh and is decomposed into two parts ( SSA and FS ; see Fig. 1). The fraction of nodes in SSA is denoted as θ with 0 < θ < 1. The number of nodes in FS is then approximately (1 − θ )N h N z , and in SSA , it is θ N 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 2N h (θ + 2(1 − θ )N z ). 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 h (θ + 3(1 − θ )N z ). The memory requirement for a physical domain in d dimensions reduces to when part of the domain is modeled by the SSA equations. The memory requirements for mesh-related quantities reduce to q mesh = 1 − θ + θ/N z in both 2-D and 3-D. The quotients q var and q mesh are close to 1 − θ if N z 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 4C FS N h N z in the whole domain, for FS 4C FS (1 − θ )N h N z in FS , and for SSA 2C SSA θ N h in SSA . The coefficients C FS and C 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 ass = 1 − θ + C SSA θ/2C FS N z . If C FS ≈ C SSA , then the reduction is approximately as in Eq. (16). The same conclusion holds in 2-D. Therefore, the reduction in that part is estimated to be similar to the reduction in Eq. (16).

Method for coupling FS and SSA
All equations are solved in Elmer/Ice  using the FEM. First the velocity 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 u 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 (Baiocchi et al., 1993), the recommended stabilization method in Gagliardini and Zwinger (2008). First, the coupling for a given geometry is presented, followed by the coupled surface evolution, both summarized in Algorithm 1.
The FS domain FS contains the grounded ice and a part of the shelf around the grounding line; see Fig. 1. The SSA domain SSA is restricted to a part of the ice shelf and starts at the coupling interface x c at the first basal mesh nodes located at least a distance d GL from the grounding line x GL , such that

Boundary conditions at the coupling interface
Horizontal gradients of the velocity are not neglected in the SSA equations (unlike in the SIA; Hutter, 1983). 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., Ahlkrona et al., 2016). 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 x c ; we assume that the cryostatic pressure acts on FS at x c : where n is normal to the coupling interface x c . The FS velocity at x c provides a Dirichlet inflow boundary condition to the SSA equations. Then, the Neumann boundary condition in Eq. (18) has to be adjusted based on the ice flow as calculated for SSA . This is done using the contact force denoted by f SSA , as explained below. The SSA equations are linearized and, by means of FEM, discretized. This leads to a matrix representation Au = b, where u is the vector of unknown variables (here, horizontal SSA velocities). In FEM terminology, the vector b that describes the forces driving or resisting ice flow is usually called the body force and A the system matrix . In Elmer/Ice, Dirichlet conditions for a node i are prescribed by setting the ith row of A to zero, except for the diagonal entry which is set to be unity, and b i is set to have the desired value (Råback et al., 2018). For an exact solution of Au = b, the residual f = Au − b is zero. If we instead use the system matrix A 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, f SSA = A SSA u SSA − b SSA is the vertically integrated contact force and needs to be scaled by the ice thickness H . In Elmer/Ice, f SSA is mesh dependent and needs to be scaled by the horizontal mesh resolution ω as well. For 2-D configurations, ω = 1. Using f SSA instead of explicitly calculating the stress is advantageous since it is extremely cheap to find the contact force if A SSA is stored.
To summarize the boundary conditions at x c , for FS, an external pressure is applied: where f SSA := 0 in the first iteration (for its derivation, see Appendix A). For SSA, a Dirichlet inflow boundary condition, provides the coupling to the FS solution. Here we take the u FS at z b , but any z can be chosen since x c should be located such that u FS (x c , z) hardly varies with z. Every iteration, f SSA , and u FS (x c , z b ) are updated until convergence up to a tolerance ε c .

Surface evolution
The surface evolution is calculated differently in the two domains FS and SSA . Equation (14) is applied to FS for the evolution of z s and z b , avoiding assuming hydrostatic equilibrium beyond the grounding line, since the flotation criterion is not necessarily fulfilled close to the grounding line (Durand et al., 2009). The thickness advection equation, Eq. (15), is used for SSA , which is advantageous since the ice flux q = H u SSA is directly available (because u 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 s and z b for SSA is then calculated from the flotation criterion, Eq. (12). At x c , H SSA = H FS is applied as a boundary condition to the thickness equation. First the surface evolution is solved for FS ; then 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 FS , SSA , or at the coupling interface x c , based on the user-defined d GL . Technically, the domain decomposition is based on the use of passive elements implemented in the overarching Elmer code (Råback et al., 2018), which allow for deactivating and reactivating elements. An element in FS is passive for the SSA solver, which means that it is not included in the global matrix assembly of A SSA and vice versa. Two kinds of iterations are involved, since computing either u FS,k or u SSA,k for the kth coupled iteration also requires Picard iteration by the nonlinearity in the viscosity. As the experiments will show, calculating u 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 u SSA,k , instead of continuing until the convergence tolerance ε P is reached, since it is inefficient to solve very accurately for u FS,k if the boundary condition at x c is not yet accurate.
Despite interrupting the Picard iteration, the final solution includes a converged FS solution since the coupled tolerance ε c is reached. Picard iteration for u SSA,k is always continued until convergence since the computation time is negligible compared to FS.
An element may switch from SSA to FS , for example during grounding line advance. Then, the coupled iteration either starts with the initial condition for u FS if the element is in FS for the first time, or the latest u FS (t) 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.

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 Greve and Blatter (2009). It consists of a 200 km long ice shelf (see Fig. 2), with a horizontal inflow velocity u(0, z) = 100 m yr −1 and a calving front at x = 200 km where the hydrostatic pressure as exerted by the sea water is applied. The shelf thickness linearly decreases from 400 m at x = 0 to 200 m at x = 200 km; z b and z s follow from the flotation criterion (Eq. 12). 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 A1.
Three models are applied to this setup, FS-only, SSA-only, and the coupled model, for which the horizontal velocities are denoted as u FS , u SSA , and u c , respectively. The relative node-wise velocity differences between u SSA and u 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 coupling location is fixed at x c = 100 km, as no grounding line is present to relate x c to. In the first coupled 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. 5).

Three-dimensional ice shelf ramp
The 2-D ice shelf ramp is extruded along the y axis (see Fig. 3). On both lateral boundaries at y = 0 and 20 km, u · n = 0. All other boundary conditions remain identical to the 2-D case, and the coupling interface is located halfway x c = (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 u · n = 0 at both lateral sides yields a negligible flow in the y direction (v FS < 10 −8 m yr −1 ). Despite differences in the models, the relative difference in u is below 1.5 %. Running the experiment with The maximum relative difference between u FS and u 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 %.

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.

Diagnostic MIS experiment
The domain starts with an ice divide at x = 0, where u = 0, and terminates at a calving front at x = L = 1800 km. An equidistant grid with grid spacing x = 3.6 km is used. Other values of constants and mesh characteristics are specified in Table A2. Gagliardini et al. (2016) showed that resolving grounding line dynamics with an FS model requires very high mesh resolution around the grounding line. However, Gladstone et al. (2017) showed that the friction law assumed in this study (see Sect. 2.2.1) reduces the mesh sensitivity of the FS model compared to the Weertman friction law assumed in Gagliardini et al. (2016), allowing the coarse mesh used here. The bedrock (m) is negative below sea level and is given by Basal melt is neglected, and the surface accumulation a s (m yr −1 ) is a function of the distance from the ice divide: www.geosci-model-dev.net/11/4563/2018/ Geosci. Model Dev., 11, 4563-4576, 2018 This experimental setup is almost equivalent to Gladstone et al. (2017), 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 (Gladstone et al., 2012b). 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 Gladstone et al. (2017) is performed, starting from a uniform slab of ice (H = 300 m), applying the accumulation in Eq. (22) 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 10 −5 ) is reached. This removal of buttressing leads to grounding line retreat from 871.2 to 730.8 km (Fig. 4).
Again, FS-only, SSA-only, and the coupled model are applied to this setup. Where u FS ≥ 5 m yr −1 , the relative difference between u FS and u SSA is below 1.8 %. The velocity u c is given in Fig. 4, with d GL = 30 km such that 58 % of the nodes in the horizontal footprint mesh are located inside SSA (θ = 0.58). The coupled model converges after 27 FS iterations on the restricted domain FS , compared to 24 Picard iterations in the FS model. The relative difference between u FS and u c is below 0.5 % (Fig. 4); this small difference shows that d GL = 30 km is sufficient. For this configuration, 4 % of the FS nodes are located between x GL and x c with d GL = 30 km; hence, decreasing d GL does not affect the proportion of nodes in FS significantly. Therefore, d GL is kept equal to 30 km for the prognostic experiment.

Prognostic MIS experiment
The prognostic experiment aims to verify model reversibility as in Schoof (2007). Starting from the steady-state geometry, the ice temperature T is lowered over a period of 500 years from −10 to −30 • C and back according to The resulting change in A (see Eq. 4) induces a grounding line advance and retreat and changes SSA by Eq. (17). Afterwards, T = −10 • 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 c and u FS after 3000 years is 10 m yr −1 (shown in Fig. 5), corresponding to a relative difference of 1.6 %. The time evolution of x GL , u b (x GL ), H (x GL ), and the grounded volume V g is shown in Figs. 6 and 7. In general, u 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 GL = 1036.8 km in the FS model and x GL = 1044 km in the coupled model. The FS model returns back to the original x GL = 730.8 km, but the coupled model yields x GL = 734.4 km, an offset of one grid point. The maximum difference in thickness is 1 %. After 3000 years, V g still decreases but the relative difference is below 10 −5 between two time steps.
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 A in Table 1), and assembly time per FS iteration is similar for the coupled and FS model. Only 5 % of the computation time is Geosci. Model Dev., 11,[4563][4564][4565][4566][4567][4568][4569][4570][4571][4572][4573][4574][4575][4576]2018 www.geosci-model-dev.net/11/4563/2018/   used to solve the linearized FS system (t s in Table 1). For all coupled simulations, assembling and solving the SSA matrix (t SSA ) takes 4 %-6 %. All the time that is left will be called overhead, t 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 o and t A , yielding a 10 % decrease in computation time with respect to the FS model. This speedup comes from a lower number of FS  (Table 1) and a slight decrease in the time used to solve the linearized FS system (13 % lower than the time that the FS model takes).

Discussion
The presented coupling is dynamic, since the coupling interface x c changes with grounding line changes, but the distance d GL that defines x 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 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 (Jouvet, 2016). 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 SSA = 0.03t FS ) and from the performance estimates in Sect. 4.2. This is due to an inefficient matrix assembly. The assembly of the system matrix A FS restricted to FS currently takes at least as much time as the assembly for the full domain , even though the domain FS is much smaller than ; in Eq. (13), θ = 0.5 for the ice shelf ramp and θ = 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 (Råback et al., 2018), 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 (Ahlkrona et al., 2016), 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 SSA . One is free to choose between the two different FS solvers in Elmer/Ice (see Gagliardini et al., 2013) 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 FS scales with the size of 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 (Gladstone et al., 2017).
In experiments involving areas where SIA is applicable, this new FS-SSA model can be combined with the ISCAL method in Ahlkrona et al. (2016) that couples SIA and FS in Elmer/Ice. This mixed model is motivated by paleosimulations, but reducing computational work by the combination of multiple approximation levels is also convenient for parameter studies, ensemble simulations, and inverse problems.
Code availability. 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 https://doi.org/10.5281/zenodo.1202407 (van Dongen et al., 2018). 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).
Use the definition of σ and the divergence theorem to rewrite Eq. (A1): The operation A : B denotes the sum i,j A ij B ij . The test function v vanishes on the inflow boundary i , has a vanishing normal component on the bedrock boundary b , and lives in the Sobolev space The space V 0 has this form because the boundary conditions on i and b are of Dirichlet type. Furthermore, there is a lateral boundary for FS ∈ R 3 , where the normal component also vanishes (v| · n = 0), and we assume a vanishing Cauchy-stress vector for unset boundary conditions to velocity components, such that the integral over vanishes. Then, the boundary integral in Eq. (A2) consists of a sum of the remaining boundary terms: The open SSA domain SSA ∈ R 2 , coupled to FS ∈ R 3 , has the boundary SSA = SSAint ∪ CF ∪ where SSAint is adjacent to FS and partly coinciding with FSint (but of one dimension less) and CF where f g = ρgH ∇ h z s and ∇ h is the horizontal gradient operator. The boundary condition on SSAint is the Dirichlet condition (Eq. 20), and the force due to the water pressure at the calving front CF is f CF , as in Eq. (13) Apply the divergence theorem to Eq. (A8) to obtain A mesh is constructed to cover FS and SSA with nodes at x i . In the finite element solution of Eq. (A9), the linear test function v i ∈ W 0 is non-zero at x i and zero in all other nodes. The integral over SSAint vanishes when v ∈ W 0 . The finite element solution u h of Eqs. (A6) and (A9) satisfies It follows from Eq. (A9) that with a test function v i ∈ W that is non-zero on SSAint and the solution u h from Eq. (A10) x i ∈ SSA ∪ CF ∪ SSAint .
The first integral in Eq. (A12) corresponds to (A SSA u SSA ) i in Sect. 3.1 and b SSAi to the second and third integrals. By www.geosci-model-dev.net/11/4563/2018/ Geosci. Model Dev., 11, 4563-4576, 2018 Eq. (A10), the right-hand side of Eq. (A12) vanishes for all x i in SSA and on CF , but for a node on the internal boundary, x i ∈ SSAint , the force f SSA from the ice due to the state u h in SSA is obtained. The internal pressure in the ice in SSA is assumed to be cryostatic as in Eq. (18). The total force on FSint consists of one component due to the state u h at SSAint and one due to the cryostatic pressure there. Let * SSA denote the mesh on SSA , which is extruded in the z direction. The common boundary between FS and * SSA is FSint , and let f * SSA be the stress force there, independent of z. Since and the corresponding strong form of the boundary condition at FSint is σ · n = H −1 f SSA − ρg(z s − z)n; cf. Eq. (19). Thus, by computing the residual as in Eq. (19), the two finite element solutions in FS and SSA are coupled together at the common boundary FSint and SSAint .