Introduction

Grounding line migration is crucial when examining the dynamics of marine ice-sheet–shelf systems (Reference Vieli and PayneVieli and Payne, 2005; Reference Pattyn, Huyghe, De Brabander and De SmedtPattyn and others, 2006). It is especially important considering that large areas of the West Antarctic ice sheet might be subject to potentially self-accelerated ice loss via the marine ice-sheet instability (Reference Joughin and AlleyJoughin and Alley, 2011). Ice of ~ 3.3 m sea-level equivalent is grounded below sea level on bedrock that is downsloping landward (Reference Bamber, Riva, Vermeersen and LeBrocqBamber and others, 2009), and is thereby potentially subject to the instability. While in one-dimensional (1-D) ice flow the grounding line is found not to be stable (Reference SchoofSchoof, 2007), the limitation of the 1-D case disregards stabilizing effects of the ice shelves (Reference Dupont and AlleyDupont and Alley, 2005; Reference GudmundssonGudmundsson, 2013). Recent observations of thinning, acceleration and grounding line retreat of the ice in parts of West Antarctica (Reference RignotRignot and others, 2008; Reference Pritchard, Ligtenberg, Fricker, Vaughan, Van den Broeke and PadmanPritchard and others, 2012) underline the demand on ice-sheet models to feature a realistic representation of grounding line motion.

In previous model intercomparison exercises, the shallow approximation models often failed to reproduce grounding line migration (Reference PattynPattyn and others, 2012), especially for relatively low resolutions of 10–20 km at which continental-scale simulations are generally performed (Reference Huybrechts and WoldeHuybrechts and De Wolde, 1999; Reference Pollard and DeContoPollard and DeConto, 2009; Reference Greve, Saito and Abe-OuchiGreve and others, 2011; Reference MartinMartin and others, 2011). Semi-analytical solutions for simplified geometries derived from boundary-layer theory based on the shallow-shelf approximation were used as benchmark (Reference SchoofSchoof, 2007). Other studies addressed the transient grounding line behavior of different flowline ice-sheet models (Reference DrouetDrouet and others, 2013) and the performance of different grounding line parameterizations (Reference Gladstone, Payne and CornfordGladstone and others, 2010). Numerical models generally have to tackle the delicate situation found at the grounding line: there is an abrupt change in the surface gradient and basal roughness, which strongly influences grounding line dynamics.

Here we analyze the resolution-dependent performance of the shallow-ice approximation (SIA)/shallow-shelf approximation (SSA) hybrid Parallel Ice Sheet Model (PISM) in reproducing grounding line motion in a model set-up and experiment sequence according to the three-dimensional (3-D) Marine Ice Sheet Model Intercomparison Project (MISMIP3d) (Reference PattynPattyn and others, 2013). This project aims to examine the capability of ice-sheet models to reproduce advance and subsequent retreat of the grounding line, i.e. its reversibility, in a 3-D ice-sheet–shelf system which is temporarily perturbed by a local basal lubrication. Our study is an update of PISM’s performance in the MISMIP3d experiments using a modified computation scheme for the driving stress at the grounding line. The migration of the grounding line is compared for two model versions of PISM (i.e. with and without applying a subgrid interpolation of the grounding line) to the performance of the full-Stokes, high-resolution, finite-element model Elmer/Ice (Reference Favier, Gagliardini, Durand and ZwingerFavier and others, 2012; Reference GagliardiniGagliardini and others, 2013). Comparison is done for PISM experiments on multiple spatial resolutions, ranging from ∆x = 1 to 16.67 km on a regular rectangular grid.

Model

The Parallel Ice Sheet Model (PISM) is an open-source, thermomechanically coupled, 3-D, finite-difference model (Reference Bueler and BrownBueler and Brown, 2009; www.pism-docs.org). It uses a superposition on the SIA and the SSA of the stress balance to calculate velocities in grounded ice (Reference WinkelmannWinkelmann and others, 2011). Since SSA velocities are used as basal velocities for grounded parts of the ice, a smooth transition of the velocity field between the grounded and floating regimes is ensured (Reference MartinMartin and others, 2011). SSA velocities and ice thickness are co-located on a fixed rectangular grid, while effective viscosity is defined on a staggered grid.

The model used in this study is based on PISM version stable 0.5. Since its release the code has been further developed including the implementation on a linear interpolation scheme for the grounding line position, referred to as ‘LI’ in Reference Gladstone, Payne and CornfordGladstone and others, (2010), where grounding line position is indicated by a fraction of gridcell length. Extension of the flowline ‘LI’ scheme to two horizontal dimensions is done by first applying it to the x and y directions separately. The fractional positions span a rectangle (visualized as overlap in Fig. 1) which is associated with the effective grounded-to-floating area ratio for each cell. Accordingly, basal resistance is proportionally increased/reduced in gridcells that are partially grounded/floating. This parameterization attenuates the discontinuity in basal resistance usually found when going from the last grounded into the first floating cell. Interpolated basal melt is neglected here. The grounding line interpolation has already been used in PISM’s contribution to MISMIP3d, but later the computation of the driving stress at the grounding line was modified: instead of a centered difference scheme across the grounding line, two one-sided differences upstream and downstream of the grounding line are now used to calculate the driving stress in the last grounded and first floating cell, respectively. This modification contributes significantly to PISM’s improved capability to represent grounding line reversibility. It should be noted that driving stress is not interpolated on a subgrid scale (which might be the case in other models). Here the ice thickness, evaluated at the center of each gridcell, determines via the flotation criterion whether the ice inside the gridcell is floating or grounded. Considering grounded cells that are neighbors of floating cells and vice versa, the above-mentioned driving-stress scheme is applied. All changes to the code are documented via the open-access revision management software github (https://github.com/pism), so reproducibility of experiments with certain model versions (e.g. code version of MISMIP3d contribution vs present study) is ensured.

Fig. 1. The two-dimensional expansion of the linear interpolation of grounding-line position on a regular grid (gray), where all data are co-located in cell centers indicated by red stars, is illustrated in three steps: (a) Linear interpolation in x –direction between cell centers of grounded and floating points reveals expected grounding line position (blue arrow). This distance is used to define a rectangular fraction of the cell area, which tends to be grounded (blueish). (b) Analogously in the y –direction, this procedure yields gridcells, which are partly covered by a blueish rectangle. (c) Overlap of the two panels is shown in dark blue rectangles, corresponding to the effective fraction of a gridcell area, and is expected to be grounded in the used formulation. A real grounding line could be located along the dashed line.

The experiments and model set-ups in this study follow the outline of the reversibility experiment of MISMIP3d (Reference PattynPattyn and others, 2013; also parameters therein) which consists of three sub-experiments (referred to as ‘Stnd’, ‘P75S’ and ‘P75R’). The model domain stretches from –800 to 800 km in the x –direction (flow direction) and from –50 to 50 km in the y –direction, with a bed elevation b = –100 – |x |/1000 (x and b in m, positive above sea level) which does not vary in the y –direction and is constant over time. Horizontal boundary conditions are a free-slip wall at y = ± 50 km and a calving front position which may not exceed y = ± 700 km. Due to the symmetry of the set-up, i.e. two symmetry axes along x = 0 (ice divide) and y = 0, respectively, only the positive ranges of the x and y axes are considered in the following.

Experiments

In the first experiment (Stnd) an ice-sheet–shelf-system is grown from an initial state of uniform ice thickness of 500 m and run into equilibrium (30 000 model years, rate of relative volume change smaller than 10^–5 a^–1, position of grounding line no longer varies; Fig. 2). For the subsequent experiment (P75S) the basal resistance is decreased locally by introducing a Gaussian-shaped perturbation at the center of the ice sheet’s grounding line (y = 0, axis of symmetry). The perturbed system is then run for 100 model years in which the grounding line is expected to advance and curve as a response to the basal lubrication. For the third experiment (P75R) the basal resistance is reset to its original constant field and the model is run for another 500 model years. According to theoretical calculations by Reference SchoofSchoof, (2007), the grounding line should exhibit reversibility, which in this setup means a retreat to its original position and again taking the shape of a straight line. The sequence of the three experiments described above is carried out on a rectangular mesh for seven different spatial resolutions: ∆x = ∆y = {1, 2, 4, 5, 10, 12.5, 16.67} km. In this study, two PISM versions are run, one without (model A) and one with (model B) applied subgrid grounding line interpolation. In both versions, the above-mentioned modified scheme for driving-stress computation at the grounding line is used. To determine changes with respect to PISM’s performance in the earlier MISMIP3d intercomparison, model A may be compared to TAL3/7 and model B to TAL1/5 in Reference PattynPattyn and others, (2013).

Fig. 2. Evolution of grounding line position (a, b) and ice volume (c, d) during 30 000 years of model spin-up for the two different model versions: (a, c) model A and (b, d) model B.

We also compare our results with results from Elmer/Ice, a full-Stokes finite-element model. It served as a reference for the MISMIP3d intercomparison, being run on very high resolution (minimal ∆x = 0.05 km) and clearly exhibiting the grounding line reversibility proposed by theory.

Resolution-Dependent Performance of Pism

Without subgrid interpolation (model A)

During spin-up (Stnd) the grounding line in model A evolves smoother and reaches a position farther downstream as resolution is increased (Fig. 2a). The resulting steady-state grounding line positions x _Stnd,_P range from 483 to 634 km.The ice volume increases towards an equilibrium value which is larger the finer the grid (highest and lowest resolution differ by a factor of ~ 1.4; Fig. 2c). The response of the grounding line to the perturbation in run P75S, i.e. its advance, which is most pronounced at y = 0, is stronger the lower the resolution (Fig. 3). Only for ∆x = 1 and 2 km (the two highest resolutions) does the grounding line retreat in experiment P75R to or even beyond its original position, hence exhibiting reversibility. For all other resolutions the resulting grounding line is located significantly downstream of its original position, being farther downstream with decreased resolution.

Fig. 3. Anomaly of the grounding line position for the PISM version without subgrid grounding line interpolation (model A): P75S – Stnd in red, P75R – Stnd in blue. Anomalies for Elmer/Ice are shown in gray (minimal ∆x = 0.05). The plots are analogous to Figure 5 of Reference PattynPattyn and others, (2013) (showing absolute values). Lower-right panel shows the difference in steady-state grounding line positions (Stnd) between PISM and Elmer/Ice for all tested resolutions. Vertical gray dashed line denotes the semi-analytical grounding line position according to Reference SchoofSchoof, (2007), calculated for the MISMIP3d set-up and parameters.

Regarding the transient response of the grounding line to the perturbation (Fig. 4), the advance/retreat for resolutions ∆x ≤ 5 km mainly occurs in the first 20 model years of the run, being smoothest for the highest resolution. For lower resolutions the grounding line migration is more step-like. It is generally known that spatial resolution is a key element in modeling grounding line migration (Reference Vieli and PayneVieli and Payne, 2005; Reference Docquier, Perichon and PattynDocquier and others, 2011; Reference Gladstone, Payne and CornfordGladstone and others, 2012; Reference PattynPattyn and others, 2013).

Fig. 4. Time-dependent position of the grounding line at the symmetry axis y = 0 (upper curves, light colors) and at the free-slip wall y = 50 km (lower curves, dark colors) during (P75S, red) and after perturbation (P75R, blue) for the PISM version without subgrid grounding line (model A). Time-span is 100 model years for both experiments. The plots are analogous to Figure 6 of Reference PattynPattyn and others, (2013).

Improvement through subgrid interpolation (model B)

For the evaluation of the subgrid grounding line positions in the x –direction the grounded-to-floating area ratio within cells along the grounding line is used. This ratio is provided by the model-computed two-dimensional subgrid interpolation of the grounding line (Fig. 1).

Regarding model B the evolution of the interpolated grounding line position during spin-up (Stnd) differs only slightly across resolutions (Fig. 2b). The resulting steady-state positions x _Stnd,_P range from 595 to 607 km (giving a width of variation an order of magnitude smaller than the model A runs; lower right panels of Figs 3 and 5). Consequently, grounding lines are close to the semi-analytical position for the SSA (Reference SchoofSchoof, 2007) of x _sa = 606.8 km, calculated for the MISMIP3d set-up and parameters. This is supported by the fact that in the simulations the SSA dominates the hybrid scheme (ratio SIA to SSA velocities <0.15 for the major part of the area covered by grounded ice). The increase in ice volume towards an equilibrium value is almost independent of resolution (Fig. 2d). After perturbation (P75S) the subgrid grounding lines show an advance of 10–26 km at the center line, and a slight retreat in the vicinity of the side margin for most resolutions. The grounding line positions after 500 model years of experiment P75R are very close to or even match their initial position and shape.

Fig. 5. Anomaly of the grounding line position as in Figure 1 for the PISM version with applied subgrid grounding line interpolation(modelB). Note different scale on x –axis.

The interpolated grounding line migrates near-continuously in time but also exhibits some variability (Fig. 6). Decreasing the resolution, oscillations between two or more grounding line positions (caused by cells in which the ice is flipping between being grounded and floating) occur more often. The magnitude of these jumps increases with grid size, but they are below gridcell length and therefore significantly smaller than the advance/retreat itself. Ignoring these numerical artifacts, the grounding line migration is similar even in magnitude to the full-Stokes solution. The fast response of the grounding line to perturbation is similar to that of other hybrid SIA/SSA and pure SSA models participating in MISMIP3d (maximum grounding line position reached after ~ 30 model years)

Fig. 6. Time-dependent position of the grounding line for the PISM version with applied subgrid grounding line interpolation (model B). Color coding same as in Figure 4, but here blue/red curves show the 11 year running mean of the yearly data underlain in gray. Note different scale on y –axis

Quantifying model response times

The analysis of the transient model response in MISMIP3d showed that the participating models differ especially in terms of response time (Reference PattynPattyn and others, 2013). To quantify response times and thus deliver a well-defined criterion to compare between different models, response functions of the form

(1)

are applied here. Equation (1) describes a system in which the response R (change in grounding line position) to a perturbation (local basal lubrication) will asymptotically approach an upper value, i.e. the response magnitude A (eventual grounding line position after 100 model years). The response time τ then gives the time after which ~ 63% of the full response is reached, allowing a quantification of the system’s sensitivity to the perturbation with respect to time. A and τ are obtained using least squares for fitting response function R to the curve of transient grounding line advance at y = 0 (yearly data) in run P75S (Fig. 7).

Fig. 7. Response function (dashed, black) according to Eqn (1) fitted to the curve of the time-dependent interpolated position of the grounding line at the symmetry axis y = 0 (run P75S, 11 year running mean in red, yearly data underlain in gray, same curves as in Fig. 6). Response time and magnitude A are given for each resolution.

The resulting response times confirm that PISM’s response is generally faster than Elmer/Ice’s (exception ∆x = 12.5 km). Regarding the three highest resolutions, for which PISM’s transient response is most consistent, the response time (τ = 6.2–7.5 years) is about half of Elmer/Ice’s response time (τ = 12.8 years).

Discussion and Conclusions

We find that, using an updated model version of PISM, the reversibility of the grounding line is captured throughout all resolutions covered here. Relative positions of advanced and retreated grounding lines are comparable to Elmer/Ice’s performance, qualitatively as well as quantitatively (Fig. 5). Especially for low resolutions, this is a major improvement compared to PISM’s performance in MISMIP3d (no reversibility at ∆x = 16.67 km). As far as we know, this is the first time that a shallow model has captured grounding line reversibility at low resolutions without applying a flux correction at the grounding line. We attribute this improvement to the combined application of the modified driving-stress calculation and the subgrid grounding line interpolation. The combination of both features further leads to a steady-state (Stnd) grounding line position being (1) much closer to Reference SchoofSchoof’s, (2007) semi-analytic position in SSA and (2) much less resolution-dependent than in MISMIP3d. Both schemes were easy to implement in PISM and are computationally cheap. Hence we think that other models might also benefit from the combined use of both features.

Without using the subgrid scheme, the relative advance of the grounding line is still qualitatively comparable to Elmer/Ice’s performance; however, reversibility only exists at a grid size of ∆x ≤ 2 km. In addition, the steady-state grounding line positions are highly dependent on resolution, which also holds for the response times of advance and retreat, revealing a more and more step-like behavior with decreased resolution. This is much more consistent among tested resolutions when applying the interpolation scheme. However, the small incremental steps of grounding line migration allow for more variability. These numerical artifacts, which are also seen in other models using a subgrid interpolation (Reference PattynPattyn and others, 2013), become stronger for lower resolutions. Considering the experimental set-up, this is not surprising as the model grid size of 16.67 km is of the same order of magnitude as the expected maximum grounding line position change of ~ 18 km. Hence, the maximum relative advance and the transient response in this specific experiment need to be interpreted with caution for the very low resolutions. PISM’s transient response is comparable to but less pronounced in terms of magnitude than in the other SIA/SSA hybrid models participating in MISMIP3d, which, in contrast to PISM, adjust the ice flux at the grounding line (Reference SchoofSchoof, 2007). Response time and response magnitude of PISM are closest to the performance of the pure SSA models. This is plausible as the SIA plays only a minor role throughout the experiments in our model (i.e. SSA velocities are dominant). Using response functions to quantify and compare model response times (Fig. 7) reveals that for the three highest resolutions PISM responds about twice as fast as Elmer/Ice (~6.5 years and 13 years, respectively). However, comparability is limited as both models start from different initial states (in terms of grounding line position and hence also ice geometry). Since we find that the initial grounding line position strongly affects model response, the above numbers have to be interpreted with caution.

We conclude that PISM’s performance in the MISMIP3d reversibility experiments improves on earlier results (Reference PattynPattyn and others, 2013), especially for resolutions of ∆x = 10 km and lower when using a modified driving-stress computation at the grounding line in addition to a subgrid grounding line interpolation. The abrupt change in driving stress and basal friction across the grounding line finds a much better numerical representation, and grounding line migration is captured more realistically. This supports PISM’s application to the Antarctic ice sheet in regional as well as continental-scale set-ups, also at resolutions of 10 km and lower.