Parameterization for Subgrid-scale Motion of Ice-shelf Calving Fronts

A parameterization for the motion of ice-shelf fronts on a Cartesian grid in finite-difference land-ice models is presented. The scheme prevents artificial thinning of the ice shelf at its edge, which occurs due to the finite resolution of the model. The intuitive numerical implementation diminishes numerical dispersion at the ice front and enables the application of physical boundary conditions to improve the calculation of stress and velocity fields throughout the ice-sheet-shelf system. Numerical properties of this subgrid modification are assessed in the Potsdam Parallel Ice Sheet Model (PISM-PIK) for different geometries in one and two horizontal dimensions and are verified against an analytical solution in a flow-line setup.


Introduction
Ice shelf fronts are predominantly observed to have an almost vertical cliff-like shape with a typical ice thickness of a few hundred meters (idealized sketch in Fig. 1a).Bending of this ice wall imposes strong tensile and shear stresses close to the terminus and promotes crevassing (Reeh, 1968;Scambos et al., 2009).Calving icebergs are cut off from the shelf along intersecting crevasses (Kenneally and Hughes, 2006) and are swept away onto the open ocean where they melt.
As a precondition for the computational treatment of calving processes and for imposing the correct boundary conditions and thereby properly computing the stress field within the shelf, we focus here on the subgrid-scale motion of the ice front on a fixed rectangular grid.One challenge arises through the finite resolution when the ice front advances seaward.Without special treatment the ice flux into a newly Correspondence to: A. Levermann (anders.levermann@pik-potsdam.de) occupied grid cell is spread out over the entire horizontal domain of the grid cell (Fig. 1b).Thus, the finite-difference ice-transport scheme (here first-order upwind) can produce grid cells of only a few meters ice thickness (or even less).In numerical models these cells are considered as floating iceshelf grid cells whose front propagates one grid cell ahead at each time step.This is generally faster than the motion of the actual moving ice-shelf margin and has no proper physical basis.The model hence produces a situation in which the ice front has no sharp vertical profile as it should have but an unphysical extension in the direction of the open ocean.In such a situation also the corresponding ice-thickness gradient which drives the ice flow is unrealistic.The dispersion effect depends mainly on the time step and spatial discretization length.It should be distinguished from the numerical diffusion of an upwind mass-transport scheme, which is often applied in finite difference models and takes the form of an additional diffusion term due to the asymmetry of the scheme.
The dispersion effect of the ice-thickness discontinuity at the front of the ice shelf is also purely numerical and does not agree with observations nor is it consistent with the underlying physical equation.Here we present a parameterization of ice-front motion on a Cartesian grid that avoids this undesirable phenomenon.There are elaborate concepts such as Immersed Boundary Methods (e.g., Mittal and Iaccarino, 2005) or Sharp Interface Methods (e.g., Marella et al., 2005), which were developed for moving boundaries in turbulent flow simulations.Our subgrid method for the slow motion of an ice front enables the application of a proper boundary condition for the stress field within the ice shelf denoted here as "calving front boundary condition" or CFBC (Weertman, 1957).Furthermore, most of the recent iceberg calving theories (e.g., Warren, 1992;Kenneally and Hughes, 2002;Van der Veen, 2002;Benn, 2007) require a steep ice wall, which is guaranteed by the proposed concept.Calving rates can be applied adequately.The paper is organized in three main parts.In Sect. 2 the features of PISM-PIK that are directly relevant for the calving front are briefly summarized.Section 3 introduces the proposed subgrid-parameterization of ice-front motion in the flow-line case and its generalization to flow in two horizontal dimensions.In Sect. 4 the parameterization is tested in simulations with PISM-PIK for a flow-line setup as well as for the Larsen and Ross Ice Shelves.We conclude in Sect. 5.

Model description
The parameterization for subgrid-scale ice-front motion, introduced in Sect.3, is applied in the Potsdam Parallel Ice Sheet Model, PISM-PIK, which is based on the thermomechanically coupled open-source Parallel Ice Sheet Model (PISM stable 0.2 by Bueler et al., 2008).Within the model, the stress balance for a floating ice shelf with negligible basal friction is computed according to the Shallow Shelf Approximation (SSA, Morland, 1987;MacAyeal, 1989;Weis et al., 1999) on a fixed rectangular grid.Solving the stressbalance equations in SSA with appropriate boundary conditions yields vertically integrated velocities, which are used for horizontal ice-transport.A full description of the model is provided by Winkelmann et al. (2010) and its performance in a setup of the Antarctic ice sheet under present-day boundary conditions is discussed by Martin et al. (2010).Here we summarize some aspects relevant for the parameterization.
The mass-transport scheme is particularly important for the ice-front motion.It approximates the ice-flux equation.In order to illustrate the general idea we restrict ourselves to the one-dimensional (flow-line) case with V , H , v x and a being ice volume, thickness, velocity and area of a grid cell (variables summarized in Table 1).
For simplicity we ignore surface and bottom mass balance.
In the vicinity of a discontinuity like a propagating ice-shelf front the appropriate numerical discretization is an upwind transport scheme.PISM base code (Bueler and Brown, 2009) uses a combination of an upwind and a centered scheme in the SSA region which does not conserve the total numerical ice mass.In PISM-PIK we introduce a first-order upwind scheme on a staggered grid which is based on the finite volume method (as generally discussed in Morton and Mayers, 2005).At the ice-front boundary the scheme has to be adjusted since there are no ice velocities on the open ocean.
In accordance with the applied conservative upwind scheme we get for the flux through the boundary (with terminal ice thickness H c and terminal velocity v c ) into the adjacent grid cell on the seaward side of the ice-shelf front The Courant-Friedrichs-Lewy criterion (CFL, Courant et al., 1928Courant et al., , 1967) ) guarantees numerical stability, i.e., the volume increment advected to the ocean grid cell is always smaller than the ice volume of the last shelf cell, This will play a role in the discussion of the treatment of residual ice volume in Sect.3. The solution of the SSA equations as a linear second order elliptical boundary-value problem requires boundary conditions.The boundary condition for the calving front (CFBC) has been implemented in most finite-difference ice-sheet and ice-shelf models in a simplified way.The ice shelves are extended artificially beyond the actual ice front with an ice thickness extrapolated from the existing ice shelf or are simply reduced to 1 m thickness (MacAyeal et al., 1996;Ritz et al., 2001).At the rectilinear boundaries of the computational domain the artificial calving front is always perpendicular to one of the two coordinate axes (x and y with indices i and j) and the dynamic Neumann boundary condition is easily imposed.Generally, for an ice front facing to the direction The Cryosphere, 5, [35][36][37][38][39][40][41][42][43][44]2011 www.the-cryosphere.net/5/35/2011/ At these positions the hydrostatic pressure term of the boundary condition (right-hand side of the equations) substitutes the velocity gradients used in the SSA equations (with ν as vertically averaged effective viscosity).In PISM-PIK we apply the dynamic boundary condition for each shelf grid cell facing the ocean to at least one side (for details see Winkelmann et al., 2010), Our subgrid parameterization guarantees a steep calving front and hence yields the correct stress balance.
In order to test the general idea of calving front advance and to find steady-state front positions using our parameterization we apply a simple calving condition that has been used in a number of previous model studies (Ritz et al., 2001;Peyaud et al., 2007;Paterson, 1994).It is based on the fact that observed ice thicknesses at calving fronts in Antarctica vary mostly between 150 and 250 m.We thus eliminate ice in any grid cell that ( 1) is located at the calving front and (2) has ice thickness less than a critical threshold H cr .The results are qualitatively independent of the specific choice of H cr for which we use a value of 250 m throughout the paper.

Parameterization
In this section we describe the subgrid parameterization of ice-front motion for both the flow-line case (one horizontal dimension) and generalize to ice flux in two horizontal dimensions.For the discretized ice-shelf model with clear-cut terminal cliff at grid cell [i] (as illustrated in Fig. 1a) the discretized flux equation Eq. ( 2) yields an ice-volume increment to be added to the adjacent ocean grid cell [i+1] with horizontal area a.The corresponding volume increment V i+1 is, without our scheme, associated with a thin ice layer of thickness H i+1 = V i+1 /a covering the whole surface of the grid cell.In our implementation a volume increment is associated with a slab ice block adjacent to the cliff (Fig. 2a).
To that end, we define a field that has the value R=1 on shelf grid cells and R=0 on ice-free ocean grid cells.Scalar values 0<R<1 in grid cells at the interface between shelf and icefree ocean are associated with the ratio of ice covered horizontal area in a grid cell with a defined reference ice thickness H r to the total grid-cell area a, which is calculated as where V is the current ice volume within the partially-filled grid cell [i+1].Thus, the slab ice block of ice thickness H r covers an area a R of the grid cell.In the flow-line case we choose this reference value to be equal to the ice thickness H r ≡H c of the adjacent shelf cell of the previous time step.In a dynamic simulation a positive ice flux through the boundary yields a new volume increment at each time step that is added to the boundary grid cell [i+1].Consequently, the area in the grid cell covered with ice and hence the ratio R increases.An increasing R is interpreted as an advancing front within a grid cell, i.e., on subgrid scale.When the ice volume in grid cell [i+1] exceeds the threshold V lim =aH r , this ice-shelf grid cell is considered to be full with R=1.In the next time step the adjacent grid cell [i+2] can start growing in ice volume (Fig. 2b).During this forward propagation of the front boundary the terminal ice thickness H c might change.Hence, the ratio R changes since H r is updated each time step even if there is no ice flux through the boundary during that particular time step.Retreat of the ice-shelf front in response to a continuous calving rate (e.g., Benn, 2007, Eq. 1) can be treated in a similar fashion as the ice advance.In the following experiments we use a simple calving condition as prescribed in Sect. 2. A generalization to other calving laws is straight forward.In such a situation, a partially filled grid cell is drained with a negative flux Q − = −c H r , where c is the magnitude of the calving rate.If the grid cell is empty with V i+1 ≤0, further ice loss acts onto the adjacent ice-shelf grid cell with R i =1+V i+1 /(a H c )≤1. Analogous to the CFL-limited numerical propagation speed of the front also the retreat is thereby restricted to at most one grid cell per time step for the simple calving rule as well as for more sophisticated calvingrate parameterizations.Negative ice volumes are set to zero after this procedure.
The generalization of Eq. ( 1) to two-dimensional horizontal ice volume flux is simply with a generalized CFL-criterion as in the PISM base code where ε is a small factor to avoid division by zero.As an example, grid cell [i,j ] in Fig. 3 borders two ice-shelf cells with velocity components directed to this grid cell on the ice-free ocean.Here, the reference ice thickness H r is the average of the ice thicknesses of those two neighboring shelf cells (or better a flux-weighted average).The volume flux through the two boundaries together with the volume of the previous time step adds up to the new volume V i,j .Herewith the ratio R i,j of ice-covered area in this grid cell is evaluated.When a subgrid ice front advances and a grid cell at the boundary is considered to become full (R=1) according to the reference thickness H r it is possible that some residual ice volume remains unaccounted for A convenient way to treat this remaining volume is to simply omit it (variant 1).Obviously this has the disadvantage of an artificial ice loss but the advantage that it does not interfere with the ice-shelf dynamics upstream of the moving ice front because H c is properly represented.Hence, the imposed CFBC, which is evaluated at the ice-shelf front and depends sensitively on the boundary ice thickness H c (Eq. 4), enables the accurate computation of velocities according to the SSA throughout the ice shelf.In order to conserve numerical ice mass and to still keep the numerical treatment as simple as possible, the residual ice mass can alternatively be equally redistributed to the neighboring grid cells on the ice-free ocean (variant 2).For these adjacent cells H r must be determined.Using adaptive time steps according to the CFL-stability criterion (Eq. 3) guarantees that the size of the volume increments V advected to the next ocean grid cell is limited.In the model of an unconfined ice shelf (e.g., flow-line case) a special problem occurs because the largest velocities are typically found at the evolving front.There, the advection of the ice-thickness discontinuity as in Eq. ( 2) has maximum propagation speed of one grid length x per time step t.Hence, we have max(V res )=a H c , which is redistributed to the adjacent grid cell [i+2].If we choose H r =H c this cell at the interface between ice shelf and ice-free ocean is completely filled within one time step (R i+2 =1), and the ice shelf evolves with an ice wall at the front that does not decrease in ice thickness, which has a strong impact on the dynamics throughout the ice shelf.In order to avoid this problem for variant 2 we reduce the reference thickness to H r,red and make a linear guess according to the analytical solution (Eq.11), which is described in the next section.The expected slope at the front depends mainly on the power law of the marginal ice thickness H c and typical constant ice parameter C and boundary values The Cryosphere, 5, 35-44, 2011 www.the-cryosphere.net/5/35/2011/This solves the problem of the unrealistic thick ice wall at the front for adaptive time steps.Note that even an inaccurate guess for the reference ice thickness jeopardizes neither mass conservation nor the basic idea of the parameterization.

Application in numerical simulations
As mentioned before, the membrane-stress balance in SSA is a non-local boundary-value problem and its over-all solution is controlled by the boundary conditions.Thus the introduction of a numerical method that alters the boundaries requires verification against an analytical solution.It is a robust feature of unconfined ice shelves that the thinning rate becomes smaller with increasing distance from the grounding line.In the model the position of the calving front of a certain ice thickness (e.g., 250 m) is sensitive to this ice thickness itself and hence small variations in ice transport have a noticeable effect.This makes the flow-line setup with applied calving rule to be a strong sensitivity test of the proposed parameterization.
For a first assessment we apply a simple ice-shelf setup, where the flow is one-dimensional in the sense that all quantities perpendicular to the flow line are constant and only the unidirectional spreading of ice is considered.We apply periodic boundary conditions at the lateral boundaries, which is associated with an infinitely broad unconfined ice shelf.At the upstream boundary, ice thickness and velocity are prescribed to 600 m and 300 m/yr.The bathymetry can be cho-sen to any arbitrary value deep enough to fulfill the floatation condition.Neither accumulation nor melting are taken into account here.There is no thermocoupling since a constant ice hardness of B 0 =1.9×10 8 Pa s 1/3 according to MacAyeal et al. (1996) is used.
The model solution is compared with the following analytical solution of the flow-line case.We choose the x-axis as the direction of the main ice flow with constant ice inflow Q 0 =v x,0 H 0 and with vanishing transversal components.In the flow-line case the stress equilibrium equations in SSA simplify considerably.Since ice is treated as a non-linearly viscous, isotropic fluid with a constitutive relation of Arrhenius-Glen-Nye form (Paterson, 1994) the equations can be integrated and rearranged with constant hardness B 0 =A −1/n 0 and flow-law exponent n=3.The solution for the spreading rate was first found by Weertman (1957) to be Inserting this into the ice-thickness Eq. ( 1), we obtain after integration the ice thickness and velocity profiles for the steady state (Van der Veen, 1999), There is no ice-shelf front considered in the analytical steady state, since lim x→∞ H (x)=0.But when appropriate boundary conditions are applied we can assume that the modeled transient profile is congruent to the analytical profile up to the advancing front at position x c .This theoretical position of the free boundary at time t can be derived from integration of Q 0 t= x c 0 H (x )dx , which yields The implementation of the calving front boundary condition (CFBC) in PISM-PIK was an essential step on the path towards reliable velocity distributions within the ice shelf.In the PISM standard code an ice-shelf extension scheme is used, where the product of effective viscosity and ice thickness is held constant.This is problematic since the calculation of the SSA yields velocities on the ice-free ocean.Due to the non-locality of the SSA equations this will generally influence velocities on the ice shelf.Its effect on the iceshelf propagation is shown in comparison to a model result with applied CFBC in the flow-line case (Fig. 4).For this experiment the subgrid parameterization (variant 1) and the simple calving rule are applied.When using the standard extension scheme with vanishing velocities at the boundary of the computational domain, the experiment shows that the velocity profile underestimates the analytical solution by more than 100 m/yr in the outer regions of the ice shelf (>100 km).Consequently, with constant ice flux across the profile, the ice thickness profile can be expected to be too thick, here about 60 m at the terminus.When applying the CFBC, however, the profile of the velocity is steeper and on average only 9 m/yr less than the the exact solution profile.Thus, calving front ice thickness of 250 m is reached at a position 10 km close to the position expected from the analytical solution (Fig. 4a, dashed).Iceshelf velocities are calculated independently of velocities on the ice-free ocean, which enhances the performance.
The CFBC gives best results when the calving front has a rectangular shape (in side view, as in Fig. 2).This shape is guaranteed with the examined subgrid treatment at the calving front.Without the subgrid parameterization though, we observe a disperion of the steep ice front with grid cells of very thin ice at the front (Fig. 1b).Hence, the CFBC is applied at a false position for a false terminal ice thickness H c .Accordingly, the velocity calculation gives false results throughout the whole ice shelf.In the following experiment we assess the influence of resolution on the steady state ice shelf for the flow-line setup with applied CFBC.We describe a three-step refinement path and divide the computational box of 505 km length into 51, 101 and 201 grid cells.When ice calves off at ice thickness 250 m the resulting profiles (Fig. 5a) show that ice fronts stabilize at distances between 160 km from the upstream boundary for coarse resolution ( x=10 km) and 145 km on a fine grid ( x=2.5 km).The latter matches the expected calving front position of 145 km calculated from the analytical solution (Eq.11).The resulting velocity distributions increase monotonically in downstream direction with largest values at the terminus.For the coarse-resolution case this value reaches 730 m/yr, while for the fine-resolution case the terminal velocity equals the analytical value of 720 m/yr (Fig. 5b).The deviation of the front position is probably due to the transport scheme which becomes more inaccurate for coarse resolution (truncation error in Taylor approximation), especially in the steep region close to the upstream boundary.Since ice thicknesses are computed on the regular grid half a grid cell (several kilometers) upstream of the defined staggered velocities, slightly stronger mass fluxes can be assumed than for finer resolution.Furthermore, the used calving rule is quite rough since ice of whole grid cells is cut off.Nevertheless, for all tested resolutions a good match of the profiles can be seen upstream of the calving front.Thus, at the expected calving front position (vertical dashed), calculated velocities underestimate the analytical value by less than 1%.Fig. 6.Transient flow-line ice thickness and velocity profiles in lateral view after 300 yr of time evolution, calculated using three different variants of residual mass treatment.The model result with adaptive time stepping is shown in blue, i.e., t≈6-10 yr.The result with fixed short time step t=1 yr is plotted in green.Magenta dashed is the analytical profile after evolution time t=300 yr.A resolution of 101×5 km is used.CFBC and subgrid parameterization are not applied for the control case variant 0.
For the transient part of the simulations in the flow-line case, when the calving front propagates downstream initiated at the boundary, different effects of the two variants of treatment of residual ice volumes are revealed.For comparison a simulation is shown where none of the two described variants of subgrid calving front treatment is used (denoted as "variant 0"), so basically the PISM extension scheme with vanishing velocities at the boundary of the computational domain.In this case, the propagating calving front suffers from strong numerical dispersion (ice thickness declines by 250 m over a distance of 80 km) especially for small time steps about 10 times shorter than adaptive time steps (Fig. 6a, green, t=1 yr).Thus, grid cells of very thin ice occur in the terminus region and the related velocities decrease in flow-line direction (Fig. 6b, green) influenced by velocity calculation in the region of the ice-free ocean and the by the shape of the dispersed front.For adaptive time steps though (Fig. 6a, blue) the front is much steeper, but a rather small numerical effect at the front is observed, so called "wiggles", which will be discussed later.Upstream of the front, in both transient cases the computed ice thickness profiles overestimate the analytical solution, which can be expressed in terms of the coefficient of determination, which is for adaptive time steps r 2 = 0.81 and for short time steps slightly better (value 1 means perfect match).Consequently, with application of calving at a certain ice thickness, a steady-state front position far beyond the analytical one can be anticipated.This is analogous to the first experiment result (Fig. 4) where the shelf-extension scheme was applied.
If we use the subgrid parameterization and cut off the occurring residual ice volumes (variant 1) we get accurately shaped profiles according to the analytical solution both in the transient phase (Fig. 6c, d) and in steady state (Fig. 5) with a coefficient of determination of r 2 >0.97 close to 1.The big advantage of this procedure is the rectangular shape of the calving front without any disturbing wiggles throughout the whole transient phase.This leads to a proper application of the CFBC and accurate velocity profiles (Fig. 6d).Variant 1 produces a certain mass loss, which can be easily reported and discussed (it is not caused by the actual transport scheme).The mass loss is negligible for small time steps (green) but it is quite large in a shelf propagating with maximal time steps according to the CFL-criterion (blue), which can be seen in the deviation of the front position in relation to the analytical front.
The mass-conserving variant 2 yields a very accurate profile (green) for short time steps ( t=1 yr) with r 2 = 0.98.Also the front position is very close to the transient analytical solution (magenta dashed), with a deviation of not more than one grid cell length.Generally, adaptive time stepping is used in order to decrease computational cost.In this special case CFL stability condition yields t= x/v c ≈ 6-10 yr and round-off errors cannot be efficiently damped in the front region.Thus, we can observe during time evolution small wiggles at the calving front (Fig. 6e, blue profile).In this case www.the-cryosphere.net/5/35/2011/The Cryosphere, 5, 35-44, 2011 the CFBC evaluated at the last shelf grid cell leads to underestimated velocity values along the ice shelf (Fig. 6f, blue), which gives a smaller r 2 = 0.84.Hence, the ice shelf is on average 37 m too thick and the front lags behind the analytically calculated front position (magenta dashed).
The subgrid-parameterization of ice-front motion with both variants of residual mass treatment is designed to be applied in two-dimensional and realistic setups as for Larsen A and B Ice Shelf as shown in Fig. 7 and for the Ross Ice Shelf in Fig. 8. Along the smooth ice front the R-field has values of 0≤R<1, while the ice shelf with values R=1 is shaded in light gray and grounded areas in dark gray.The figures show a steady state snapshot with applied continuous physical calving rate, but details are discussed elsewhere (Levermann et al., 2011).
In a two-dimensional realistic setup of a confined ice shelf or even in a setup of the Antarctic ice sheet with several ice shelves attached the adaptive time steps are determined according to the maximal velocity magnitude of the whole computational domain (Eq.7).The generalized CFLcriterion is used to limit the amount of residual ice mass, which is redistributed equally to the neighbor grid cells in an unphysical way with regard to the physical flow across the boundary.We could apply a more rigorous criterion, but simulations in realistic setups confirm that the error is small.The maximal flux through the boundary for a certain adap-  tive time step occurs for a single pair of cells, which is located probably at the ice front with distance from confinements, whereas along the rest of the ice-shelf front velocities are lower than the maximal value.Numerical damping is relatively efficient in most front regions as well as in the inner parts of the ice shelf.Hence, transient phenomena like wiggles at the front are rarely observed.

Discussion and conclusions
In this paper we presented a numerical method that enables the subgrid motion of ice-shelf fronts in a finite-difference model.This prevents the steep margin from being numerically dispersed and allows for a proper application of the Neumann boundary condition for the approximated stressbalance calculations.Flow-line simulations with the Potsdam Parallel Ice Sheet Model (PISM-PIK) for different resolution have been compared with the exact analytical solutions.The modification of the transport scheme at the icefront boundary, which implicates a redistribution of residual ice volumes at this moving front has been assessed and discussed.The proposed procedure opens the way to an appropriate determination and application of calving rates in realistic models of combined ice-sheet/ice-shelf dynamics.
In the simple flow-line case a very good match of the calculated velocity profile and the analytical solution confirm a correct application of the CFBC.Comparison of a simulation with a realistic setup of the Ross Ice Shelf against observed velocity data from the Ross Ice Shelf Geophysical and Glaciological Survey (RIGGS, Thomas et al., 1984;Bentley, 1984) verifies the proper implementation also in two horizontal dimensions even when the boundary line is not straight (not discussed in this paper).An additional benefit in using the implemented CFBC is that velocity calculation is independent of information outside of the ice-shelf boundaries.Thus, velocities on the ice-free ocean can be ignored, which simplifies the calculation of SSA-velocities and reduces computational cost.Very well approximated steadystate ice thickness and velocity profiles with a coefficient of determination of r 2 >0.99 are observed for all of the three tested resolution.The steady-state front position converges to the analytical solution for resolution refinement.In simulations of the Antarctic ice sheet (e.g., Martin et al., 2010), coarse resolution of about 20 km grid length or more are used and marginally overestimated mass fluxes can be expected here in the ice-stream and shelf region.
The subgrid parameterization of ice-shelf-front motion implicates the handling of residual ice-volume increments that arise from the restriction R≤1 of the ice coverage ratio.We show in transient simulations that variant 1 yields very accurate flow-line profiles for both ice thickness and velocity, although the modeled front lags behind the analytical front due to the cut-off of residual ice volumes.We use this variant for the application of calving rates that depend sensitively on the velocity field in the vicinity of the front (e.g., Levermann et al., 2011).Thereby, the residual ice volumes are reported as additional mass loss, although they are not physically motivated.Note that in the flow line case with adaptive time stepping (Fig. 6c, blue) variant 1 can produce large mass losses.In a more realistic case, however, these mass losses are far smaller, comparable to the flow-line case of shorter time steps (as in Fig. 6c, green) since the residual ice volumes are generally much smaller and the cut-off occurs less often.The mass-conservative variant 2, however, yield accurate results and does not increase computational cost distinguishably since the CFL-criterion limits the size of the volume increments and the numerical redistribution of residual ice volumes to adjacent grid cells on the ocean is generally executed only once each time step.

Fig. 1 .
Fig. 1.(a) Schematic of a discretized ice shelf is shown in lateral view with decreasing ice thickness in positive i-direction.The local calving front is located at the interface between the last shelf grid cell [i] and the adjacent open ocean cell [i+1].(b)In every time step a volume increment is calculated for each grid cell according to the scheme approximating Eq. (1).Thus, in every time step, the marginal cliff moves one grid cell further into the open ocean and may thin out relatively fast.

Fig. 2 .
Fig. 2. Lateral view of discretized ice shelf with subgrid-scale parameterization.(a) The volume increment in grid cell [i+1] is associated with a slab of ice of the same ice thickness as the adjacent full shelf grid cell [i].(b) When the grid cell [i+1] at the calving front is full of ice according to the associated reference thickness, the next following cell [i+2] can gain in ice volume.

Fig. 3 .
Fig.3.A bird's eye view of the ice shelf calving front approximated by a rectangular mesh grid.Gray shaded area denotes ice-free ocean, the ice shelf area is white with exemplary velocity vectors defined on the regular grid.Velocity components directed to the open ocean are shown in green.

Fig. 4 .
Fig.4.Profiles of ice thickness (top) and velocity (bottom) in flowline case in lateral view after evolution time 1000 yr.Fixed Dirichlet boundary on the left and calving front on the right hand side of the computational domain.Blue: the result with applied CFBC already in steady state.Green: with shelf-extension scheme and still advancing.Black dashed is the expected profile from analytical solution in steady state when cut off at 250 m ice thickness.In both cases a resolution of 101×5 km, adaptive time stepping and variant 1 of residual mass treatment are used.

Fig. 5 .
Fig. 5. Steady-state flow-line ice thickness and velocity profiles in lateral view calculated with different resolution but constant size of computational domain.Calving front position (145 km) and terminal velocity (720 m/yr) as expected from the analytical solution with constant ice hardness for 250 m terminal ice thickness are shown as black dashed line.Adaptive time steps and variant 1 of residual ice treatment are used.

Fig. 7 .Fig. 7 .
Fig. 7. Realistic steady-state model simulation of Larsen A and B Ice shelf (light gray) with grounded parts (dark gray) and the ice-free ocean (white).Values of R at the propagating ice shelf front are colored.

Fig. 8 .
Fig. 8. Realistic steady-state model simulation of Ross Ice shelf (light gray) with grounded parts (dark gray) and the ice-free ocean (white).Values of R at the propagating ice shelf front are colored.

Fig. 8 .
Fig. 8. Snapshot of a realistic steady state model simulation of Ross Ice shelf (light gray) with grounded parts (dark gray) and the icefree ocean (white).Values of R at the propagating ice shelf front are colored.

Table 1 .
Table of variables and abbreviations.