Improvements in one-dimensional grounding-line parameterizations in an ice-sheet model with lateral variations

The use of a boundary-layer parameterization of buttressing and ice flux across grounding lines in a two-dimensional ice-sheet model is improved by allowing general orientations of the grounding line. This and another modification to the model's grounding-line parameterization are assessed in two settings: a narrow fjord-like domain (MISMIP+), and in future simulations of West Antarctic ice retreat under RCP8.5-based climates. The new modifications are found to have significant effects on the 10 fjord results, which are now within the envelopes of other models in the MISMIP+ intercomparison. In contrast, the modifications have little effect on West Antarctic retreat, presumably because dynamics in the wider major Antarctic basins are adequately represented by the model's previous simpler one-dimensional formulation. As future grounding lines retreat across very deep bedrock topography in the West Antarctic simulations, buttressing is weak and deviatoric stress measures exceed the ice yield stress, implying that structural failure at these grounding lines would occur. We suggest that these grounding-line 15 quantities should be examined in similar projections by other ice models, to better assess the potential for future structural failure.


Introduction
Analytic boundary-layer treatments of buttressing and ice velocities across grounding lines (e.g., Schoof, 2007) are usually 1-D, i.e., formulated with one horizontal dimension along the flowline and no lateral variations. In ice-sheet models with two 20 horizontal dimensions, such formulations can be used to prescribe the approximate flow across grounding lines. In our previous work (Pollard and DeConto, 2012;DeConto and Pollard, 2016), this was done simply by applying the 1-D expressions at individual one-grid-cell-wide segments separating pairs of grounded and floating cells, so that the orientation of each single-cell "grounding-line" segment is parallel to either the x or the y axis. Although this is consistent with the one-dimensional character of the formulation in Schoof (2007), it neglects the actual orientation of the wider-scale grounding line, and results in non-25 isotropic buttressing amounts for x vs. y directions.
Here we implement a more rigorous, isotropic treatment of grounding-line buttressing and ice flow, by applying the 1-D expressions to normal flow across a more realistic grounding-line orientation that is not constrained to one or the other grid axes.
In principle this is more rigorous than the previous single-cell treatment, and is expected to improve model results. We also implement a modification to improve behavior with strong buttressing, and a minor improvement in crevasse depths, as 30 described below.
Two types of experiments are used to assess the above modifications. First, simulations are performed for a fjord-like glacier confined to a relatively narrow channel, as in the MISMIP+ intercomparison (Cornford et al. 2020). Because of the confining https://doi.org/10.5194/gmd-2020-131 Preprint. Discussion started: 25 May 2020 c Author(s) 2020. CC BY 4.0 License. lateral boundaries and a central bedrock depression, grounding lines in these simulations have large two-dimensional curvatures, and provide a good test for the changes implemented here. Second, much larger-scale simulations of future ice retreat in West

35
Antarctica are performed, forced by warming climates corresponding to the extreme RCP8.5 greenhouse gas emissions scenario.
In section 2, the modifications to the buttressing, grounding-line flux and crevasse-depth parameterizations are described in detail. Section 3 presents results for the narrow fjord-like MISMIP+ experiments, and section 4 presents results of the West Antarctic future simulations, comparing results for different combinations of the new model modifications. In section 5, deviatoric stresses at grounding lines in West Antarctic simulations (without hydrofracturing or cliff-failure physics) are 40 examined, as they retreat across very deep bedrock topography in central West Antarctica in future centuries, to assess the potential for structural failure that could lead to very rapid disintegration of the remaining ice. In an Appendix, two additional more speculative modifications to the model's grounding-line flux parameterization are described.

45
As described in Pollard and DeConto (2012), the primary grid in the finite-difference ice-sheet model is the h-grid, with ice thicknesses (h) defined at the center of each cell. Ice in each h-grid cell is either floating in the ocean or is grounded, depending on the ice thickness, bedrock elevation, and sea level. At the grid-cell level, the boundary between floating and grounded-ice regions consists of piecewise-linear segments at the edges between pairs of h-grid cells, with each edge parallel to the x or y axis (Fig. 1a). The model uses an Arakawa-C grid, in which horizontal u and v velocities are staggered half a grid cell in the x and y 50 directions respectively, so each segment separating floating or grounded h-cells has a u or v velocity defined at its mid-point, as indicated in Fig. 1a. (The model performs a sub-grid interpolation that refines the grounding-line position between each pair of h-grid cells and does not coincide with the cell edge between them, but that does not affect the material presented here). The model ice dynamics uses a hybrid combination of vertically integrated shallow ice and shallow shelf approximations (SIA, SSA), with the seaward ice flux at grounding lines imposed as a boundary condition according to an analytical expression 60 relating ice flux to ice thickness (Schoof, 2007): https://doi.org/10.5194/gmd-2020-131 Preprint. Discussion started: 25 May 2020 c Author(s) 2020. CC BY 4.0 License.
where qg is the ice flux and Ug is the ice velocity across the grounding line, and h is ice thickness at the grounding line. ρi and ρw are the densities of ice and ocean water, respectively, and g is gravitational acceleration. A is the rheological coefficient and n is 65 the exponent for ice deformation. C is the coefficient and m is the exponent for basal sliding (Schoof, 2007), written as Cs and ms in Pollard and DeConto (2012). The term θ in (1a) represents buttressing by ice shelves, i.e., the amount of back stress caused by pinning points or lateral forces on the ice shelf further downstream. The buttressing factor θ is defined as the ratio of vertically averaged horizontal deviatoric stress normal to the grounding line, relative to its value if the ice shelf was freely floating with no lateral constraints and no back stress. (The latter free-floating value is always extensional, balancing the difference between the 70 column-mean hydrostatic ice pressure at the grounding line with the smaller mean horizontal component of ocean-water pressure on the ice shelf. Pinning points or lateral forces on the ice shelf reduce this value towards zero, i.e., less extensional and more compressive, so θ = 1 for unbuttressed grounding lines and diminishes towards 0 as buttressing increases.) The analysis for grounding-line flux and buttressing in Schoof (2007) is limited to one-dimensional flowline geometry. In our standard model (Pollard and DeConto, 2012), Eq. (1) is applied across individual one-grid-cell-wide segments separating pairs of 75 grounded and floating grid cells, so that the orientation of each single-cell "grounding-line" segment is parallel to either the x or the y axis, as sketched in Fig. 1a. In the standard model, the buttressing factors θu and θv in the x and y directions respectively are: where η is the non-linear strain-dependent ice viscosity, and the numerators in (2a,b) are 2x the deviatoric stress (times ice 80 thickness h) in the x or y directions.
Although this previous treatment of θ is consistent with the one-dimensional character of the formulation in Schoof (2007), it does not capture the wider-scale orientation of the real grounding line, which does not actually run along the "staircase" singlecell segments as in Fig. 1a. Furthermore, it results in non-isotropic θ values for the u and v staggered-grid velocities.
An alternative treatment of θ is described below, and sketched in Fig. 1b. The new method allows for general grounding-line 85 orientations running at an angle to the grid axes, and applies the ice flux given by (1) in a direction normal to this grounding line.
First, an estimate of the grounding-line orientation is needed, that represents a spatial smoothing of the boundaries of nearby cells. A simple algorithm is used, as follows.
(i) Consider all grid cells within a given radius Rc of the location in question (xc, yc), and take the average of the x and y coordinates of cells with ocean or floating ice (not grounded ice), (xo, yo). If this radius extends beyond the domain boundaries, virtual points are used with their grounded or floating property equal to that extended normally from the domain boundary.
(ii) Then the normal to the grounding line (in the direction towards the ocean) is (xoxc, yo -yc). The length of this vector is normalized to 1 meter, and is called (nx,ny) below. This orientation is used in the calculation of N, the net deviatoric stress normal to the grounding line. The equations below follow 100 Gudmundsson (2013, his Eqs. 2, 6 and 12).

110
Then the buttressing factor θ is given by The denominator is the net normal deviatoric stress that would result for a freely floating and completely unbuttressed ice shelf (or a vertical ice face with no ice shelf at all). θ is used in Eq.
(1) to obtain the imposed ice flux qg and velocity Ug normal to the grounding line.

115
Finally, Ug (in the direction (nx,ny)) is resolved into its x-and y-axis components: These velocity components are imposed in the final SSA solution at each time step, at staggered u or v-grid points as appropriate located at the mid points between pairs of grounded and floating h-grid cells. (It is easy to show that this decomposition of Ug 120 onto the u and v-grids results in the physically correct net flux across the actual grounding line, averaged over many u and v-grid points).
As well as entering in the Schoof grounding-line flux Eq. (1a), the buttressing factor θ also affects the effects of crevasses and hydrofracturing in grounding-zone cliff-failure . These physics are not enabled for all MISMIP+ runs and most of the Antarctic runs below.

Grid-cell weighting of imposed grounding-line velocities
If the buttressing factor θ falls to zero or below in (6), this corresponds to net zero or compressive horizontal deviatoric stress normal to the grounding line, and zero or compressive (negative) strain in the direction of flow. However, its use in the Schoof formation for grounding-line ice velocity (1) predicts zero velocity for θ = 0, and becomes invalid for θ < 0, as noted by Reese et al. (2018). In our application, θ in Eq. (1a) is restricted to be greater or equal to 0, but can still unrealistically cause the Schoof

130
velocity Ug to become very small or zero if θ falls to 0.
To avoid this problem, here we modify the condition (Box 9 in Pollard and DeConto, 2012) that selects the grid-cell location at which the Schoof grounding-line velocity components ug or vg from (7) are imposed. First, we describe this condition in the standard model. The following is written with all variables and indices for u velocities on the u-grid, but it applies equally to v velocities on the v-grid.

135
The condition affects whether ug is applied at the u-grid point separating a pair of floating vs. grounded h-grid cells, or at the next downstream u-grid point. The choice depends on the difference ugugrid, where ugrid is the "grid-solution" velocity from a preliminary solution of the SSA equations without Schoof constrains, as mentioned above. The grounding line is assumed to be retreating if ug >> ugrid, and be advancing if ug << ugrid. If the former, ug is imposed at the "actual" grounding-line location between the pair of grounded vs. floating h-grid points, and if the latter, ug is imposed at the next downstream (floating) location.

140
This rule is ad-hoc, but yields reasonable grounding-line migration in idealized tests and real-world simulations. It works because in probable advancing situations with relatively small ug, it is applied at the downstream edge of the first floating h-grid cell, generally allowing it to thicken towards grounding; and in probable retreating situations with relatively large ug, it is applied at the downstream edge of the last grounding h-grid cell, generally allowing it to thin towards flotation.
To implement this condition, a weighting factor w (0 to 1) is used: where the normalizing velocity Unorm = (10 5 m 2 yr -1 ) / hg, and hg is the ice thickness at the grounding line. In early model versions described in Pollard and DeConto (2012, Box 9), this was imposed as a binary choice with w either 0 or 1, i.e., as if Unorm in (8) https://doi.org/10.5194/gmd-2020-131 Preprint. Discussion started: 25 May 2020 c Author(s) 2020. CC BY 4.0 License.
was very small. In more recent model versions the gradual weighting in (8) has been used, which was found to yield smoother behavior in some tests.

150
The SSA equations for (u,v) over the whole domain are solved with the following impositions at the two grounding-zone points mentioned above, denoted here by u-grid indexes iau between a grounded vs. floating h-grid pair, and ibu at the next u-grid point downstream (all on the x-axis, with the y-axis index j not shown). In the standard model used to date:

155
In the SSA sparse matrix and right-hand side of the equation for u at grid point (iau, j), the first quantity of the pair in Eq. (9a) is used to weight u towards ug with weight w, and the second quantity indicates that all other quantities (which would yield the non-Schoof "grid solution") are multiplied by 1-w. The same is done for u at point (ibu, j), using Eq. (9b) with w and 1-w switched.
Here, we introduce a slight modification to the weighting, aimed at avoiding the problem of ug becoming zero as θ falls to zero in Eq. (1):

160
At iau: At ibu: The only difference from the standard model is the use in (10b) of ugrid(iau) at point ibu instead of ug. This term is weighted by 1w, so is significant in probable advancing situations with small ug (w close to 0); the use of ugrid(iau) at ibu still allows the first floating h-grid cell to thicken sufficiently towards grounding. In the old model, as ug falls to zero and w falls to zero, the imposed 165 velocity at ibu in (9b) also falls to zero. In the new model, the imposed velocity in (10b) does not tend to zero, which is more reasonable. Instead, the divergence ∂u/∂x between iau and ibu tends to zero, as it should for strong buttressing causing θ to fall towards zero.

Calculation of crevasse depths
For all runs in this paper, an improvement is made in the parameterization of crevasse depths, used both in "normal" calving and 170 also in the cliff-failure physics . Crevasse depths are set to the Nye-depth (at which total horizontal stress is zero for surface crevasses, or is equal to water pressure for basal crevasses; Nye, 1957;Jezek, 1984;Nick et al., 2010).
Previously, the divergence (∂u/∂x + ∂v/∂y) was used along with ice viscosity as a simple estimate of the horizontal deviatoric stress . Here, this is replaced by the maximum principal deviatoric stress (Turcotte and Schubert, 1982), calculated from the strain rates and viscosity. This is a small improvement "in principle". It has no effect in the idealized fjord

175
MISMIP+ experiments for which calving is disabled, and has little effect on the West Antarctic simulations (not shown).
As a first test of the modifications above, we use the MISMIP+ experiments (Cornford et al., 2020). These simulate glacier flow in a rectangular fjord-like channel, and involve significant two-dimensional curvatures of grounding lines. The channel is 80 km wide, with bedrock generally sloping downstream and an ice shelf flowing into the ocean. There is a bedrock depression at mid-  (Cornford et al, 2020). Also shown is the "control" grounding line after spinup at year 0 (thick black line). Axes scales are in km. The first 80 km of the channel is not shown. Note the nearly 3x stretching of the channel width relative to the length. With our standard model, i.e., without any of the modifications above, the grounding-line variations were significantly faster and 195 larger than other higher-order, higher-resolution models in the intercomparison (Cornford et al, 2020). Here, we perform various sensitivity experiments with one or more of the modifications above: A. the previous model with single-cell "staircase" grounding lines (Eq. 2).
C. the 2-D orientation, and the improved grid-cell weighting in section 2.2 (Eq. 10). envelopes of other higher-order, higher-resolution models in the MISMIP+ intercomparison for Ice1 (shown as background 205 shading in Fig. 3). This suggests that the modifications above are real physical improvements to our model.

Results: West Antarctic simulations
To test the modifications in real-world scenarios at larger scales than the idealized fjord experiments above, we simulate retreat of the West Antarctic ice sheet due to future climate warming. The climate forcing follows that in DeConto and Pollard (2016) for the extreme RCP8.5 greenhouse gas emissions scenario, with atmospheric temperatures and precipitation from regional climate model simulations, and oceanic temperatures from a transient future simulation with the NCAR CCSM4 global climate 245 model (Shields et al., 2017). The ice sheet is initialized to modern observed (Fretwell et al., 2013), and run from 1950 CE for 500 years. A nested domain is used spanning West Antarctica with a polar stereographic grid of 10 km resolution, and with lateral boundary conditions supplied by an earlier continental-scale run.
The mechanisms of hydrofracturing and cliff-failure DeConto and Pollard, 2016) are disabled in the main simulations below, so the future collapse of West Antarctica is relatively slow and driven mainly by sub-ice-shelf oceanic melt 250 and ductile processes as in other models (Feldmann and Levermann, 2015;Golledge et al., 2015;Arthern and Williams, 2017).
This provides a better test of the modifications above, without the overall retreat being dominated by more drastic retreat mechanisms.

Results: potential for structural failure at West Antarctic grounding lines
As grounding lines retreat across central West Antarctica in the RCP8.5-driven simulations above, they encounter very deep 290 bathymetry with depths of ~1 to 2.5 km below sea level, especially in the Bentley Subglacial Trench (Fig. 10). Simple vertically integrated force balance calculations (Bassis and Walker, 2012;Pollard et al., 2015) and vertically resolved modeling (Bassis and Jacobs, 2013;Ma et al., 2017;Schlemm and Levermann, 2019;Benn et al., 2019;Parizek et al, 2019;cf. Clerc et al., 2019) suggest that ice columns at such deep grounding lines, if unbuttressed or only weakly buttressed by ice shelves or mélange, will be structurally unstable, with deviatoric stresses exceeding the material yield stress of the ice. Once initiated, structural "cliff" 295 failure would be expected to propagate extremely rapidly into ice upstream of the grounding line, only stopping when shallower bathymetry is reached, or if buttressing increases somehow.
where τx´x´ is the depth-averaged normal deviatoric stress at the grounding line (in direction x´ to distinguish it from the model's x axis). Note that this applies equally to grounding lines with ice shelves, and to ice cliffs at grounding lines without an ice shelf (for which θ =1). The crevasse depths ds (surface) and db (basal) are Nye-depths as described in section 2.3 above, and depend on 310 principal deviatoric stress. Their sum ds + ds = θ h / 2, where h is the ice thickness .
Assuming x´ is close to the horizontal principal stress direction, the quantity 2τx´x´ is a good approximation for the difference in the two principal stresses in the x´and z plane, which is reported in laboratory experiments as a measure of ice yield strength, typically around ~1 MPa (Bassis and Walker, 2012). Several other considerations may modify this value and the concept of a uniform ice yield strength itself (Parizek et al., 2019;Clerc et al., 2019), including deformation unique to cliffs such as slumping 315 and torques, ice cohesion and modes of failure depending on depth, and importantly, the amount of pre-existing fractures, buried crevasses, bubbly ice and/or cm-scale grain sizes, as opposed to relatively pristine ice with small (~mm-scale) grain sizes. Ice with extensive pre-existing damage is prevalent in most ice cores and presumably throughout Antarctica, and has yield strengths around ~1 MPa, much weaker than pristine ice; Parizek et al. (2019) and Clerc et al. (2019) agree that maximum heights of subaerial ice cliffs (above sea level, with ~9 times that below sea level) are approximately 100 to 200 m for pre-damaged ice, and ~500 m for pristine ice. In our diagnosis below, we assume that central West Antarctic ice is typically pre-damaged (or if it is not already, it will likely become so as the rapidly retreating grounding line approaches from the north), and so assume an ice yield strength of around ~1 MPa.
In Fig. 11, θ and the relevant deviatoric stress measure (2τx´x´ from Eq. 11) are plotted at grounding lines for the simulation without hydrofracturing or cliff failure (type 2), and for the model version C with both new modifications. For modern, 2τx´x´ is 325 far below 1 MPa at all grounding lines, as it should be as no significant structural failure is observed today. At +400 years into the run (~2350 CE), when the retreating central West Antarctic grounding lines are beginning to encounter deep (>1 km) bathymetry, the surviving ice shelves shown in Fig. 8 still provide some buttressing, and most buttressing factors are well below 1; (even though these ice shelves are too short and thin to reach distant pinning points, the lateral curvature in their flow produces back stress). Most 2τx´x´ values are somewhat below 1 MPa, indicating that extensive structural failure is unlikely. However, by +500 years (~2450 CE), central West Antarctic grounding lines experience even deeper bathymetry, and many θ values are at or close to 1 (weakly buttressed or essentially unbuttressed). Many 2τx´x´ magnitudes are at or exceed 1 MPa, indicating that structural failure of these grounding-line columns would occur.

340
The modifications described and tested above in (i) 2-D orientation of the grounding line in calculating the buttressing factor and imposed normal ice flow, and (ii) grid-cell weighting of imposed grounding-line velocities, are physically reasonable. The first modification more rigorously represents the true geometry of the grounding line.
In the idealized fjord-like MISMIP+ experiments, which involve strong 2-D curvature of grounding lines in a narrow channel, the modifications have significant effects on the model's grounding-line variations, bringing them in line with those of other 345 higher-order higher-resolution models in the MISMIP+ intercomparison (Cornford et al., 2020). This suggest that the modifications are real improvements to the model physics.
In contrast, the modifications have relatively little effect in large-scale simulations of future rapid West Antarctic ice retreat. This is presumably because of the larger lateral scales of major West Antarctic basins, so that grounding-line retreat in these basins is more one-dimensional in character, and better represented by the simpler "staircase" grounding-line treatment of the standard 350 model.
The more rigorous treatments of grounding-line orientation and buttressing factors allow us to better diagnose the force balance at grounding lines in the West Antarctic simulations, to see if structural failure could occur in a future with unmitigated greenhouse-gas warming. We find that when grounding lines reach very deep central West Antarctic regions (~1 to 2.5 km below sea level) after about 500 years, ice-shelf buttressing is weak and the deviatoric stress measures widely exceed the ice yield 355 stress, implying that structural failure would occur at these grounding lines. In that case, a runaway disintegration could be initiated, with structural failure propagating very rapidly into the remaining grounded ice (Schlemm and Levermann, 2019), which in the absence of renewed buttressing would continue until shallower bathymetry is reached to the south.
Several other ice sheet-shelf models have performed similar projections of future West Antarctic retreat (e.g., Feldmann and Levermann, 2015;Golledge et al., 2015;Arthern and Williams, 2017), some with higher order and/or higher resolution than 360 ours. We suggest it would be beneficial to examine these grounding-line quantities in other model simulations, to more robustly assess the danger of structural failure in future centuries under RCP8.5-like climate warming.
Apart from ice shelves, another potential source of buttressing is from mélange. Huge amounts of floating ice debris (mélange) would be generated in front of the retreating ice fronts in the above scenarios. In major Greenland fjords today such as Jakobshavn and Helheim, mélange is considered to provide significant back stress on the glacier calving front, at least in winter 365 (e.g., Burton et al., 2018). However, in one study using a heuristic continuum model of mélange  Other processes that could reduce the deep bathymetry encountered by future grounding lines are bedrock rebound under the reduced ice load, and less gravitational attraction of the ocean by the receding ice (Gomez et al., 2015). The West Antarctic 370 simulations here include the first process, using a relatively simple ELRA (Elastic Lithosphere Relaxing Asthenosphere) bed model (Pollard and DeConto, 2012), and the rebound of the modern bathymetry (Fig. 10) under the central grounding lines after 400 to 500 years (Fig. 11) is minor. However, recent geophysical data indicate very low mantle viscosities below parts of West Antarctica (Heeszel et al., 2016;Barletta et al., 2018), which could produce faster rebound and shallower bathymetry by the time grounding lines retreat into central regions. Work to develop Earth-sea level models with laterally varying properties and ice-375 ocean gravitational interaction, and couple them with ice-sheet models, is ongoing (Gomez et al., 2018;Powell et al., 2020).

Appendix
A few recent analytical studies have investigated aspects of boundary-layer treatments of grounding-line zones, going beyond Schoof (2007), such as Haseloff and Sergienko (2018) and Sergienko and Wingham (2019). Here we briefly test two modifications to our grounding-line flux implementation that are more heuristic and speculative than those in the main paper.

380
They are roughly motivated by the recent work although they cannot represent it directly and address different aspects.

A1. Strain softening
This modification addresses the presumed underestimate of strain softening in the grounding zone in a purely 1-D flowline treatment such as Schoof (2007

390
(e.g., Thoma et al., 2014), where LM is either LM N or LM N , and A and n are the rheological coefficient and exponent respectively appearing in the Schoof formula Eq. 1a. In our implementation η is computed using the velocity solution of the previous iteration (Pollard and DeConto, 2012), at the last grounded cell adjacent to the grounding line.
This is not rigorous because the Schoof analysis incorporates the 1-D dependence (A1) in its derivation, and not (A2). However the modification to A in (A4) is at least in the right direction (increasing the ice flux across the grounding line), and may be useful as a crude approximation.

400
Sergienko and Wingham (2019) found that in ice streams with high basal sliding coefficients, the boundary-layer expansion of Schoof (2007) is not valid, and can overestimate the flux of ice across the grounding zone. Following on from that paper, the ratio of the newly calculated flux to the Schoof-calculated flux, in idealized tests for small basal slopes, ranges from ~0.6 to 1 but can be much smaller for steeper slopes (O. Sergienko, personal communication, 2020). This analysis cannot be represented by modifications in our model. However, we can crudely estimate the possible effect of such changes for small basal slopes at least,

405
by simply reducing all imposed grounding-line velocities in (1) by a constant factor, i.e., multiplying Ug given by (1b) by a factor 0.6.

A3. Effects on results
The effects of applying each of the modifications described above are shown here. Figs   Code and Output Availability. The ice sheet model code is available on request from the corresponding author 435 (pollard@essc.psu.edu). That and selected model output will be available at Penn State's Data Commons, http://www.datacommons.psu.edu.
Author contributions. DP and RD conceived the project and design. DP performed coding and simulations and wrote the manuscript with input from RD.
Competing interests. The authors declare that they have no conflict of interest.