The Potsdam Parallel Ice Sheet Model (pism-pik) – Part 1: Model Description

We present the Potsdam Parallel Ice Sheet Model (PISM-PIK), developed at the Potsdam Institute for Climate Impact Research to be used for simulations of large-scale ice sheet-shelf systems. It is derived from the Parallel Ice Sheet Model (Bueler and Brown, 2009). Velocities are calculated by superposition of two shallow stress balance approximations within the entire ice covered region: the shallow ice approximation (SIA) is dominant in grounded regions and accounts for shear deformation parallel to the geoid. The plug-flow type shallow shelf approximation (SSA) dominates the velocity field in ice shelf regions and serves as a basal sliding velocity in grounded regions. Ice streams can be identified diagnostically as regions with a significant contribution of membrane stresses to the local momentum balance. All lateral boundaries in PISM-PIK are free to evolve, including the grounding line and ice fronts. Ice shelf margins in particular are modeled using Neumann boundary conditions for the SSA equations, reflecting a hydrostatic stress imbalance along the vertical calving face. The ice front position is modeled using a subgrid-scale representation of calving front motion (Albrecht et al., 2011) and a physically-motivated calving law based on horizontal spreading rates. The model is tested in experiments from the Marine Ice Sheet Model Intercomparison Project (MISMIP). A dynamic equilibrium simulation of Antarctica under present-day conditions is presented in Martin et al. (2011).


Introduction
In order to understand the evolution of ice sheets, especially with respect to their contribution to sea-level rise, there is a need for numerical models which are able to capture the dynamics of sheet-shelf systems as a whole.To this end, there are various types of approaches: Flowline models are both computationally efficient and at the same time very useful for understanding basic processes (e.g., Dupont and Alley, 2005).However, one of the key issues regarding global sea-level rise is the buttressing effect of ice shelves on the sheet which influences the ice flux across the grounding line.In flowline models, this can only be investigated by prescribing a parameterized buttressing strength.Three-dimensional models, by contrast, compute the effect of buttressing for any embayment.
Models based on the Shallow Ice Approximation (SIA, Hutter, 1983) and the Shallow Shelf Approximation (SSA, Morland, 1987;Weis et al., 1999), each recalled in Sect.2, decrease these computational costs.They make use of a scaling analysis of the full Stokes problem, using the fact that ice thickness is small compared to relevant horizontal scales, to enable simulations of whole ice sheets.
SIA-only models (e.g., Greve, 1997;Payne, 1999;Marshall et al., 2000) provide a good approximation for the Published by Copernicus Publications on behalf of the European Geosciences Union.
flow of ice that is frozen to the bedrock.In order to simulate regions of higher velocity such as ice streams, sliding of grounded ice is incorporated in various SIA models by adding a sliding velocity at the base of the ice column which depends directly on the gravitational driving stress.
The SSA has been used to simulate the flow in ice shelves which have no basal friction and thus a different stress regime (e.g., MacAyeal et al., 1996).
The two shallow approximations have been combined using different approaches: some hybrid models (e.g., Huybrechts, 1990;Ritz et al., 2001) solve the SIA and the SSA in distinct map-plane regions, which results in an abrupt change in flow type across the grounding line.Pollard and Deconto (2009) combine both shallow approximations via a heuristic shear stress correction and an additional grounding line velocity correction based on the onedimensional approach by Schoof (2007b).
Here, we present the Potsdam Parallel Ice Sheet Model (PISM-PIK), developed at the Potsdam Institute for Climate Impact Research (PIK), a new hybrid model for marine ice sheets.It combines the two approximations by adding the SIA and the SSA velocity on the whole sheet and shelf system in order to incorporate the different flow regimes in sheet, streams and shelves in a universal manner.
PISM-PIK is based on the Parallel Ice Sheet Model (PISM; Bueler and Brown, 2009), which is a threedimensional thermodynamically-coupled shallow model using a finite-difference discretization.These models are innovative in using the SSA as a sliding law for grounded ice, thereby avoiding discontinuities at the onset of sliding, and provide a framework for consistently modeling sheet-shelf systems.Stream-like flow is modeled by the SSA stress balance and a plastic till model for basal mechanics (Schoof, 2006a).In aiming at modeling the Antarctic ice sheet and shelves (Martin et al., 2011), PIK modifications to PISM have been made particularly with respect to the shelf dynamics.Ice shelves are especially important for an assessment of past and future sea level contributions of Antarctica because of their buttressing effect on the dynamics upstream of the grounding line and its position (De Angelis and Skvarca, 2003;Bamber et al., 2007;Rott et al., 2007;Glasser and Scambos, 2008;Rignot et al., 2008;Goldberg et al., 2009).
In Sect. 2 we describe PISM-PIK with special focus on the modifications compared to PISM.Section 3 shows its performance in experiments from the Marine Ice Sheet Model Intercomparison Project (MISMIP).A summary is given in Sect. 4. The performance of PISM-PIK under present-day boundary conditions for Antarctica is discussed in Martin et al. (2011).

Model description
Our description of the model proceeds from the underlying continuum model to the numerical schemes.The major char-acteristics of the continuum model are adopted from PISM (Sect.2.1), but modifications include changes to the way sliding occurs (Sect.2.2) and major changes in the treatment of marginal boundaries of the ice sheet and shelves, in particular calving fronts .
Boundary data provided to PISM-PIK are bed elevation b, surface mass balance M, geothermal flux G and ocean temperature T o .In simulations for Antarctica such as the dynamic equilibrium simulation in Martin et al. (2011), the surface temperature distribution is a parameterized function of latitude and elevation, and the model is initialized with observed ice thickness data from Lythe et al. (2001) and Le Brocq et al. (2010).In each time step, the model is solved for velocity v, basal melt rate S and calving rate C. Each time step updates the prognostic variables, which are ice temperature T and thickness H .Because the bed elevation b is fixed in time, each update to the ice thickness implies an update to surface elevation h.(For the notation throughout this paper see Table A1.) Each of the model equations is numerically discretized, using a finite difference approach on a rectangular grid.Exactly as described in Bueler and Brown (2009), the new model PISM-PIK runs in parallel using the Portable, Extensible Toolkit for Scientific computation (PETSc) for solving every discretized model equation, including the SSA.Time stepping is explicit and adaptive.Technically speaking, PISM-PIK is implemented as derived classes of the C++ code of the open-source Parallel Ice Sheet Model (PISM), version stable0.2(The PISM authors, 2011).

Field equations and shallow approximations
For completeness we give a short review of the field equations used both in PISM and PISM-PIK.For more details see Bueler and Brown (2009) and Bueler et al. (2007).
The equation of mass continuity describing the evolution of the ice thickness can be derived from incompressibility ∂ µ v µ = 0 and the kinematic equations at the surface and base of the ice.(Throughout this paper, the Einstein summation convention is used.)Here Q is the vertical integral of the horizontal velocity.The mass continuity equation is solved numerically for H .The vertical velocity v z within the ice is also given by incompressibility.
The following shallow equation of conservation of energy includes advection and vertical conduction of heat and a strain dissipation heating term (Greve and Blatter, 2009, Eq. 5.105) The Cryosphere, 5, 715-726, 2011 The strain heating ( ) is a function of the second invariant of the strain rate tensor, which is -modifying the approach by Bueler and Brown (2009) -approximated by an unweighted sum of the strain rate tensors of SIA and SSA.Details about the combination of SIA and SSA in PISM-PIK compared to PISM are given in Sect.2.2.Surface temperature T s and the geothermal flux as well as pressure melting temperature in the case of ice shelves provide boundary conditions for this equation.
Equation ( 3) is solved for the temperature on a zcoordinate system, where the z-coordinate measures the vertical distance above the bedrock for grounded ice and above the ice shelf base for floating ice, so the zero level for z is always the base of the ice.This grid choice, and thermomechanical verification results from PISM using it, are documented in (Bueler et al., 2007).
The flow law is given by Here, ˙ is the strain rate tensor and τ the deviatoric stress tensor with τ ij =σ ij +pδ ij where the pressure p is the isotropic part of the full stress tensor σ .Throughout this paper, the tensor indices i and j denote the two horizontal components of vertically-integrated terms.Following Paterson and Budd (1982), ice softness A(T * ) is given by A(T * )= (3.61×10 −13 )e −6.0×10 4 /( RT * ) , T * ≤263.15K, (1.73×10 +3 )e −13.9×10 4 /( RT * ) , T * >263.15K. (5) The flow law has an inverted form in which the effective viscosity ν depends on the effective strain rate and the temperature.The Stokes stress balance for a viscous fluid (Greve and Blatter, 2009, Sect. 5.1) is approximated using two limiting shallow approximations.
The Shallow Ice Approximation (SIA, Morland and Johnson, 1980;Hutter, 1983) is valid for regions where bottom friction is high enough for the vertical shear stresses to dominate over the horizontal shear stresses and longitudinal stresses.The corresponding velocities are given by The second shallow approximation is the Shallow Shelf Approximation (SSA, Morland, 1987;Weis et al., 1999) Here, τ b = (τ b x ,τ b y ) denotes the basal shear stress.
The vertically-averaged effective viscosity is given by where B is the vertically-averaged ice hardness, ˙ is the strain rate tensor and E SSA is an enhancement factor.
The SSA stress balance is non-local, i.e., when solving for the velocity at a certain grid point, the solution will depend on the whole spatially-distributed stress field.
In an ice shelf where there is zero traction at the base of the ice, the driving stress is exclusively balanced by membrane stresses (those stresses held by viscous deformation), and the basal shear stress is set to zero (τ b =0).The resulting plug flow is described by the SSA, see the upper panel in Fig. 1 (ice shelf ).
Following MacAyeal (1989), the Shallow Shelf Approximation is also used for modeling the fast flow regime in ice streams (see the upper panel in Fig. 1; ice stream).Ice streams have a flow-regime similar to the one in ice shelves (a plug flow) but experience basal resistance.Thus, Eqs. ( 8) and ( 9) are used with nonzero basal friction τ b , which is calculated based on a model for plastic till (Schoof, 2006a) This basal model assumes that till supports applied stresses without deformation until these equal a yield stress τ c .Division by zero is avoided by the addition of a small constant in the denominator (Bueler and Brown, 2009).Yield stress is given by the Mohr-Coulomb model for saturated till, with till cohesion c 0 set to zero (Paterson, 1994) The porewater pressure is parameterized by where λ ranges from 0 to 1 and can be interpreted as water content in the till.A specific parameterization of λ depending on bed elevation and sea level is given for Antarctica in Martin et al. (2011).
www.the-cryosphere.net/5/715/2011/The Cryosphere, 5, 715-726, 2011  Ice sheet with negligible basal sliding.The velocity is dominated by the Shallow Ice Approximation.Friction at the bed leads to shearing in the ice column and a velocity profile as depicted in blue.Ice stream.Both SIA and SSA velocities are relevant for the ice flow.The SSA serves as a sliding law for the shear flow.Ice shelf.The plug flow regime is dominated by the SSA leading to a uniform vertical velocity profile.Lower panel: an example cut through the sheet-stream-shelf transition of Lambert Glacier and Amery Ice Shelf in the equilibrium simulation described in Martin et al. (2011).The model output velocity v=v SIA +v SSA is shown in black.The onset of an ice stream is diagnosed to be where the SSA velocity (red) exceeds the SIA velocity (blue).

Velocity combination, sliding and grounding line motion
Following the basic idea of superposition of SIA and SSA from Bueler and Brown (2009), both SIA and SSA velocities are computed on the whole model domain in PISM-PIK, enabling a smooth transition both from non-sliding ice which is frozen to the bedrock to faster-flowing ice in regions where significant sliding occurs, and from rapidly-sliding ice across the grounding line to floating ice (indicated in Fig. 1).Bueler and Brown (2009, Eqs. 21 and 22) use a weighting function in order to combine the SIA and SSA velocities in the PISM base version where The weighting function f (|v SSA |) ensures a continuous solution of the velocity from the interior of the ice sheet across the grounding line to the ice shelves.It is approximately 1 for small |v SSA | and approximately 0 for large |v SSA | so that the SIA velocity fully governs the overall velocity where the SSA velocity is small (in the interior of the ice sheet).However, the choice of the weighting function is an additional degree of freedom which is not constrained by observations and is therefore not used in PISM-PIK.Instead, the contributions resulting from the two approximations are simply added, still ensuring a smooth transition of the velocity across the grounding line.Thus in PISM-PIK the basal velocities for grounded ice are the SSA velocities v b =v SSA and On ice shelves, experiments with simplified setups have shown that the SIA contribution is negligible due to the low The Cryosphere, 5, 715-726, 2011 www.the-cryosphere.net/5/715/2011/surface gradient so that the dynamics there are dominated by the SSA.In the inner part of the sheet where bottom friction is high enough for vertical shear to dominate over horizontal shear, the SSA contribution is negligible so that the ice velocity there is dominated by the SIA contribution.This SIA/SSA hybrid scheme provides an alternative to SIA-only thermomechanically-coupled sliding models, which allow sliding using evolving basal resistance fields and are generally subject to the failures described qualitatively and quantitatively in Appendix B of Bueler and Brown (2009).In these cases, where sliding is incorporated by adding a basal velocity v b to the SIA velocity v SIA in Eq. ( 7), unbounded vertical velocities arise from (physicallypossible) abrupt changes in basal strength (Fowler, 2001;Bueler and Brown, 2009).The ISMIP-HEINO experiments (Calov et al., 2010), in particular, demonstrate the difficulty such models generally experience with dynamicallyevolving basal resistance fields.Through the superposition of the two velocity fields, each of which are differentiable, unbounded vertical velocities cannot arise with the hybrid scheme used in PISM-PIK.
The basic, mathematical task of model verification is carried out in Bueler and Brown (2009) and Bueler et al. (2007) and references therein.But especially the question whether the dual approximation used here is indeed a valid approximation of the Stokes equations in the transition zone needs to be deferred to other publications.A comparison of the superposition approach to Blatter-Pattyn type models or other "hybrid solvers" like in Pollard and Deconto (2009), Schoof and Hindmarsh (2010), and Goldberg (2011) and performance in ISMIP-HOM would be most interesting as well but is beyond the scope of this paper.
The SIA/SSA hybrid scheme enables the model to simulate qualitatively different flow regimes -fast as well as slow processes.Martin et al. (2011) show that the velocities in an equilibrium simulation of Antarctica span several orders of magnitude, from almost zero close to the ice divides over velocities of a few meters per year in large areas of the interior of the ice sheet to the highest velocities of sliding grounded ice which are associated with stream-like features and reach a few kilometers per year.
We diagnose an ice stream as a region where the SSA velocities are larger than the contribution from the SIA velocities to the vertically averaged overall velocity v, i.e., This definition is purely diagnostic and no different physics are prescribed.It serves as a simple way of identifying the parts of the sheet with significant sliding.The lower panel in Fig. 1 from the dynamic equilibrium simulation discussed in Martin et al. (2011) illustrates how the modeled velocity undergoes a qualitative change from non-sliding over sliding to floating ice.Here, the velocity contributions from SIA (in blue) and SSA (in red) are shown for an example cut through the modeled Lambert Glacier and Amery Ice Shelf.An ice stream is identified as the grounded region where the SSA velocity supersedes the SIA velocity and is understood as the transition region between the sheet without significant sliding where the velocity is almost entirely given by the SIA velocity and the floating shelf where the SSA velocity clearly dominates.
The second merit of the hybrid scheme is the stress transmission across the grounding line.In PISM-PIK, the grounding line is not subject to any boundary conditions.Its position is determined in each time step via a mask which distinguishes grounded ice from floating ice using the flotation criterion The grounding line motion is thus influenced indirectly by the velocities through the ice thickness evolution.Since the SSA velocities are computed non-locally and simultaneously for the shelf and for the sheet, a continuous solution over the grounding line without singularities is ensured and buttressing effects are accounted for.In Sect.3, the reversibility of the grounding line migration on a sloping bed with PISM-PIK is demonstrated in experiments from the Marine Ice Sheet Model Intercomparison Project (MISMIP).The variability of the grounding line under climate forcing is additionally demonstrated in an experiment shown in the supplement of Martin et al. (2011) where a sinusoidal surface temperature forcing is applied to an equilibrium state of Antarctica, resulting in forward and backward motion of the grounding line.

Discretization scheme for mass transport
PISM uses a mass continuity scheme for SIA fluxes which has perfect numerical mass conservation (Bueler et al., 2005).The low-order mass continuity scheme used for SSA fluxes, as described in Bueler and Brown (2009), however, can lead to local mass conservation errors.Keeping track of all mass fluxes and comparing their sum to the change in total ice volume has shown that especially at ice margins problems arise.
We employ a modified scheme in PISM-PIK that is locally mass conserving and that applies to both SIA and SSA velocities in the same way.As in PISM, in general the SSA velocities are computed on the regular grid whereas the SIA velocities are computed on a staggered grid, i.e., on a grid which is shifted by half a grid length compared to the regular grid.For the discretization of the mass transport, the SSA velocities are transferred onto that same staggered grid by averaging over the SSA velocities from the adjacent grid cells.The sum of the SIA and SSA velocities on the staggered grid is the total velocity which is used in a mass-conserving upwind finite difference scheme for the mass continuity Eq. (1): www.the-cryosphere.net/5/715/2011/The Cryosphere, 5, 715-726, 2011 An analogous equation holds for the y-component.This scheme ensures mass conservation since every contribution of mass inflow in one box is necessarily a contribution to mass outflow in an adjacent box.
At the ice front the scheme needs to be slightly modified, since we do not compute velocities in ocean-boxes without ice.To avoid the zero velocity contribution of the ocean box adjacent to the last shelf box, when computing the staggered velocity at the calving front, we use the unstaggered SSA velocity from this last shelf box and the front stress boundary condition.
The properties of this alternative scheme for the pure-SIA mass continuity problem have been tested by comparison to the similarity solution described by Halfar (1983), and the deviations of this solution when using the PISM-PIK scheme are of the same order of magnitude as the ones when using the PISM base version (Bueler et al., 2005).For the SSA, Albrecht et al. (2011) have shown that the analytical solution for the flow line case from Van der Veen (1983) for an ice shelf in equilibrium is better approximated with this alternative mass transport scheme than with the scheme described in Bueler and Brown (2009).It is not necessary to modify the adaptive time-stepping from PISM base version to ensure the numerical stability of this scheme.

Calving front stress boundary condition
At the calving front, the imbalance between verticallyintegrated stresses in the ice and the hydrostatic pressure exerted by the ocean contributes to the ice flow upstream of the calving front.A calving front stress boundary condition is necessary for solving the non-local SSA momentum balance Eqs. ( 8) and ( 9).PISM-PIK therefore incorporates a physical stress boundary condition (Weertman, 1957;Morland and Zainuddin, 1987) describing the stress imbalance at the calving front.By contrast, in PISM there is no such condition, and instead a notional ice shelf extends across the ice-free ocean to the edge of the computational domain.The choice of the ice thickness and viscosity used for such an extension and the effect on the rest of the modeled ice are inadequately understood (MacAyeal et al., 1996;Bueler and Brown, 2009).The stress boundary condition implemented in PISM-PIK reads where z sl denotes the sea-level, H c is the ice thickness at the calving front and n=(n x ,n y ) the horizontal, seawardpointing vector normal to the calving front.
The SSA stress balance Eqs. ( 8) and ( 9) can be expressed in terms of a vertically-integrated stress tensor (Schoof, 2006b) This gives the SSA stress balance a very compact form which is equivalent to Eqs. ( 8) and ( 9) combined: The deviatoric part of the left hand side of the stress boundary condition can then be expressed in terms of the vertically integrated deviatoric stress tensor Neglecting air pressure and assuming cryostatic pressure in the ice we get The right-hand side of these equations, representing the static part of the stress balance τ stat , is computed according to the type of ice front concerned.In PISM-PIK, we distinguish between three types of ice-ocean interfaces, namely shelf calving fronts where marine ice fronts (ice at the coast resting on bedrock below sea level) where and cliffs (ice at the coast resting on bedrock above sea level) where These stress boundary conditions are applied directly at the calving front by replacing certain terms in the discretization of the SSA equations with the respective terms from the boundary condition, making the use of an artificial shelf extension obsolete.The implementation holds for any shape of calving front and is described in detail in Appendix A.

Continuous ice shelf advance and retreat through subgrid parameterization
In order to capture subgrid-scale advance and retreat of the calving front as depicted in Fig. 3 the mechanism introduced by Albrecht et al. ( 2011) is implemented in PISM-PIK.
The subgrid parameterization is a precondition for the application of a continuous calving law like the one described in Sect.2.6 because it models the observed almost vertical cliff-like shape of the calving fronts and preserves realistic driving stresses and hence spreading rates near the ice edge.
Without this mechanism the ice flux into a newly occupied grid cell would be spread out over the entire horizontal domain of that cell, possibly resulting in ice shelf grid cells of only a few meters ice thickness or less.The front of these cells would propagate one grid cell ahead in each time step, leading to an unphysical extension of the ice shelf onto the ocean.
In particular, for an advancing ice front, the discretized version of the mass continuity equation yields a volume increment V for each cell at the ice front that needs to be added to the adjacent ocean grid cell.Without the subgrid parameterization, the volume increment would cover the whole area of this adjacent ocean grid cell and the newly formed shelf would be unrealistically thin.In order to avoid this purely numerical effect, in PISM-PIK a field R is introduced which gives the ratio of ice-covered horizontal area in each gridcell depending on a reference ice thickness H r , equal to the mean ice thickness of the adjacent full shelf grid cells of the previous time step: where a is the area of the grid-cell and R takes values between 0 and 1, R = 1 corresponding to full ice-shelf cells and R = 0 to ice-free ocean cells.Thus, for a cell with 0 < R < 1, the area Ra is covered with ice of the same thickness as the mean of the adjacent full shelf grid cells as illustrated in Fig. 3. Concerning the mass budget for the partially-filled cells, the outflow is determined by the calving rate while the inflow is computed from the ice thickness at the front and the velocities on the staggered grid described in Sect.2.3.
A retreat of the calving front can be modeled in the same manner: first, the calving rate is determined through the calving law described in Sect.2.6.The amount of mass lost through transport by this calving rate is removed by changing the ice volume in the partially-filled cells at the ice front, and, as needed, by transforming adjacent filled cells into partiallyfilled cells.

Calving law
Based on the observation by Doake et al. (1998) and Alley et al. (2008) that the calving rate is linearly-related to the product of near-front thickness, half-width and strain rate, a local first-order law for large-scale ice shelf calving was introduced by Levermann et al. (2011) as a boundary condition to the mass continuity equation.This so-called Eigen Calving law has been implemented in PISM-PIK.
The applied calving rate is based on the eigenvalues ˙ ± of the horizontal strain rate tensor (see Eq. 4).Along most areas of the calving front the corresponding eigen-directions will coincide with the directions parallel to and transverse to the flow (Fig. 4).In regions of divergent flow, where spreading occurs in both principle directions (˙ ± >0), we define the rate of large-scale calving as with K>0 being a proportionality constant.Else, the calving rate is zero.Typically, the maximum spreading rate in an ice shelf can be found along the ice-flow direction downstream of, but close to, the grounding line.If the calving rate were to depend only on ˙ + , it would therefore increase towards the grounding line so that no stable ice-shelf front would develop.Stability arises through the dependence on the product of ˙ + and ˙ − since calving is prohibited where ˙ + ˙ − is negative.
Spreading rate fields and therefore areas of potential calving according to Eq. ( 30) are strongly influenced by the geometry of an shelf ice.In areas where it comes to convergence of ice flow the calving rate is zero whereas it strongly increases where ice flow expands, e.g., at the mouth of an embayment, hindering the shelf to grow outside the embayment.Martin et al. (2011) show that this new dynamically motivated calving law enables PISM-PIK to reproduce realistic calving front positions for many of the ice shelves attached to www.the-cryosphere.net/5/715/2011/The Cryosphere, 5, 715-726, 2011   the Antarctic ice sheet, across a broad range of ice thickness values.
Due to the calving law, an ice bridge connecting the main part of an ice shelf and a wider tongue can be calved off and the thereby separated tongue forms an iceberg of the size of a few grid cells.For these floating pieces with zero total basal resistance the SSA does not have a unique velocity solution (Schoof, 2006b).Therefore they are identified and eliminated, and their volume is reported as lost through calving.
Since no subgrid interpolation is implemented for marine fronts and ice cliffs, mass loss occurs as follows in these cases: in each timestep, a small amount of ice is transported from a grid cell at the marine front or ice cliff into the adjacent ocean grid cell.The hereby emerging shelf cells contain typically very thin ice and therefore calve off immediately.

Experiments from the Marine Ice Sheet Model Intercomparison Project
The numerical performance of PISM-PIK was tested in the context of the Marine Ice Sheet Model Intercomparison Project (MISMIP, Schoof et al., 2009).The results of these simulations allow a comparison with the semi-analytical solution by Schoof (2007a), and give insight into the quality of the numerical treatment of grounding line motion.Note that grounding line migration results directly from the flotation criterion and is thus fully determined by the dynamics described in Sect. 2. We consider only the MISMIP flow-line experiments on a constant slope bed (experiments 1 and 2) in this paper.As PISM-PIK is a three-dimensional model, a flow-line was simulated using periodic boundary conditions in the crossflow direction.Calving occurs at a fixed position in MISMIP so the calving law in Sect.2.6 was not used.The model was run at resolutions of 12 km, 6 km and 3 km.Upon changes in flow parameter A (Eq. 5), the position of steady-state grounding line should closely follow the semi-analytical solution, depicted as a solid black line in Fig. 5a.In order to diagnose the position of the grounding line at subgrid scale, a refinement based on Pattyn et al. (2006) was introduced for the analysis of the MISMIP experiments.
During experiment 1 the softness parameter A is lowered stepwise.During experiment 2 these changes are reversed.Figure 5a shows the modeled grounding line position.It can be seen that (i) the movement of the grounding line is modeled qualitatively correctly at all resolutions, (ii) an offset from the semi-analytical solution is observed in experiment 1, and (iii) the offset is smaller in experiment 2. The offset in experiment 1 decreases with increasing resolution, while in experiment 2 it is independent of the resolution.The dependence on the resolution becomes evident in Fig. 6 which shows the steady-state grounding-line position of experiment 1 for various resolutions.The simulated grounding line positions converge to the semi-analytical solution upon gridrefinement.For a flow-line setup governed by the SSA as considered here, the existence of a shelf should not influence the grounding-line movement in the continuum model (Schoof, 2007a).In the implementation of the problem in PISM-PIK, however, we observe some (numerical) influence of the calving-front treatment on the steady-state results at low resolutions.This is due to the fact that the SSA stress balance equations are solved simultaneously for the entire computational domain encompassing sheet and shelf.
The implementation of the calving front stress boundary condition (see Sect. 2.4) in PISM-PIK distinctly improves the performance compared to PISM which uses a notional shelf extension in ice-free ocean.This shelf extension method generates a less-accurate backstress on the whole ice body influencing the grounding line position.The subgrid interpolation of the calving front as described in Sect.2.5 ensures a controlled shelf growth, while without this treatment a very thin shelf develops immediately.Figure 6 shows that grounding line position performance with a subgrid treatment of the calving front is generally better than without, especially at low resolutions.

Conclusions
The large-scale marine ice sheet model PISM-PIK presented in this paper is based on the Parallel Ice Sheet Model (PISM), with major modifications for ice shelf dynamics.A boundary condition for the SSA equations capturing the stress imbalance at the ice front was incorporated and a novel subgrid  scheme for the advance and retreat of the calving front was introduced.This was combined with a physical calving law which is universal in the sense that it governs the calving process for different shelf geometries.
Grounded ice flow is captured as a direct superposition of velocities from two shallow approximations, SIA and SSA, a modification of the Bueler and Brown (2009) scheme.Both shallow approximations are solved simultaneously in the whole ice area.The new hybrid scheme allows for the treatment of all flow regimes in a consistent manner and permits a diagnostic definition of ice streams.PISM-PIK thus includes a treatment of the transition from vertical-shearing dominated flow, in areas where the ice is frozen to the bed, to plug flow in ice shelves.All this governs the ice flux across the grounding line whose position is determined directly from the flotation criterion, so no boundary conditions are imposed to model its motion.Since the SSA velocities are computed non-locally and simultaneously for shelf and sheet, stress transmission across the grounding line is guaranteed.
The performance of PISM-PIK was tested in context of the Marine Ice Sheet Model Intercomparison Project.In these experiments, the model reproduced the semi-analytically predicted grounding line position dependence on ice softness for the flow-line case, with convergence towards the semianalytical solution with grid refinement.
As a first application, Martin et al. (2011) present a dynamic equilibrium simulation for Antarctica, demonstrating the ability of PISM-PIK to simulate whole sheet-shelf systems.

Appendix A On the discretization of the calving front stress boundary condition
In the following we detail the implementation of the calving front boundary condition into the SSA-scheme.We only describe the scheme for SSA Eq. ( 8).The second equation based on Eq. ( 9) is similar.Discretizing the outer derivatives of Eq. (8) yields The next step is to replace the four expressions of the form 2 ∂{ij } with the respective terms from the boundary condition (24).Using the boundary-indicators a ± ,b ± (see Fig. 7) leads to where a ± ,b ± =0 if one or both of the adjacent cells are icefree, and a ± ,b ± =1 otherwise, and where the value of τ stat  • is determined by one of Eqs. ( 26), ( 27) or ( 28).The calving front is situated at the border between ocean cells and the last completely-filled shelf cells; the stress boundary condition is not applied to partially-filled cells that arise due to the subgrid parameterization described in Sect.2.5.Further discretization yields: For each specific shape of calving front, the respective boundary indicators are set to zero, so that the derivatives of v y in Eq. ( 8) (and of v x in Eq. 9) are partially neglected.

Fig. 1 :
Fig. 1: Superposition of SIA and SSA velocities.Upper panel: Schematic diagram of ice profile and different flow regimes.Throughout the grounded ice region, the ice velocity is determined by linear superposition of the SIA (blue arrows) and the SSA velocities (red arrows).Ice sheet with negligible basal sliding.The velocity is dominated by the Shallow Ice Approximation.Friction at the bed leads to shearing in the ice column and a velocity profile as depicted in blue.Ice stream.Both SIA and SSA velocities are relevant for the ice flow.The SSA serves as a sliding law for the shear flow.Ice shelf.The plug flow regime is dominated by the SSA leading to a uniform vertical velocity profile.Lower panel: An example cut through the sheet-stream-shelf transition of the Lambert Glacier and Amery Ice Shelf in the equilibrium simulation described inMartin et al. (2011).The model output velocity v=v SIA +v SSA is shown in black.The onset of an ice stream is diagnosed to be where the SSA velocity (red) exceeds the SIA velocity (blue).

Fig. 1 .
Fig.1.Superposition of SIA and SSA velocities.Upper panel: schematic diagram of ice profile and different flow regimes.Throughout the grounded ice region, the ice velocity is determined by linear superposition of the SIA (blue arrows) and the SSA velocities (red arrows).Ice sheet with negligible basal sliding.The velocity is dominated by the Shallow Ice Approximation.Friction at the bed leads to shearing in the ice column and a velocity profile as depicted in blue.Ice stream.Both SIA and SSA velocities are relevant for the ice flow.The SSA serves as a sliding law for the shear flow.Ice shelf.The plug flow regime is dominated by the SSA leading to a uniform vertical velocity profile.Lower panel: an example cut through the sheet-stream-shelf transition of Lambert Glacier and Amery Ice Shelf in the equilibrium simulation described inMartin et al. (2011).The model output velocity v=v SIA +v SSA is shown in black.The onset of an ice stream is diagnosed to be where the SSA velocity (red) exceeds the SIA velocity (blue).

Fig. 2 .
Fig. 2. Schematic diagram of forces at the calving front.Hydrostatic pressure (gray arrows) and cryostatic pressure (black arrows) result in a net force in direction of the ocean (blue arrow).The net force is significantly reduced compared to a situation where the hydrostatic pressure is neglected.

Fig. 3 :
Fig. 3: (From Albrecht et al., 2011) Schematic diagram illustrating subgrid-scale advance of the calving front.As a new grid cell i+1 is partially filled with ice its ice thickness is kept constant at the same value as grid cell i. 23

Fig. 3 .Fig. 4 .
Fig. 3.(From Albrecht et al., 2011)  Schematic diagram illustrating subgrid-scale advance of the calving front.As a new grid cell i+1 is partially filled with ice its ice thickness is kept constant at the same value as grid cell i.

Fig. 5 .
Fig. 5. MISMIP experiments on a constant slope bed.(a) Dependence of grounding-line position on rate factor A. Solid black line: semi-analytical prediction.Colored lines: simulation results with PISM-PIK, at x=12, 6, 3 km resolution.(b) Example ice sheet profile.

Fig. 7 :Fig. 7 .
Fig. 7: Stencil for the SSA boundary condition in the case of an ice-filled (white) cell with ice-filled and ice-free (shaded) neighbors.For this specific shape of the calving front the boundary indicators a − ,b + and a −N ,a +N ,b +W have the value 1, whereas a + ,b − and a −S , a +S ,b −W ,b −E ,b +E have value 0.
; v b =v SSA =(v x ,v y ) m s −1 v SIA SIA velocity of ice m s −1 v SSA SSA velocity of ice; v SSA =(v x ,v y ) tensor; σ ij =τ ij −pδ ij Pa=N m

Table A1 .
Table of symbols.