New Mass-Conserving Bedrock Topography for Pine Island Glacier Impacts Simulated Decadal Rates of Mass Loss

High-resolution ice ﬂow modeling requires bedrock elevation and ice thickness data, consistent with one another and with modeled physics. Previous studies have shown that gridded ice thickness products that rely on standard interpolation techniques (such as Bedmap2) can be inconsistent with the conservation of mass, given observed velocity, surface elevation change, and surface mass balance, for example, near the grounding line of Pine Island Glacier, West Antarctica. Using the BISICLES ice ﬂow model, we compare results of simulations using both Bedmap2 bedrock and thickness data, and a new interpolation method that respects mass conservation. We ﬁnd that simulations using the new geometry result in higher sea level contribution than Bedmap2 and reveal decadal-scale trends in the ice stream dynamics. We test the impact of several sliding laws and ﬁnd that it is at least as important to accurately represent the bedrock and initial ice thickness as the choice of sliding law.


Introduction
Numerical models of marine ice streams-that is, ice streams flowing over bedrock lying well below sea level such as the major glaciers of West Antarctica-are accurate if they can properly represent the dynamics of the region around the grounding line. The grounding line itself separates ice thick enough to contact bedrock from the floating ice shelves and takes a leading-order role in the evolution of the entire ice stream. The ice flux across the grounding line is strongly dependent on the thickness of the ice there (Schoof, 2007), and this flux, balanced with the integrated flux upstream, determines the rate of sea level change from the ice sheet. As a result, having an accurate representation of the ice sheet geometry and in particular the bedrock elevation, at a sufficiently high resolution, is essential.
However, there is often a disparity in resolution between the ice sheet geometry products used as input, and the models themselves. By "ice sheet geometry product" we mean a compilation of bedrock elevation and ice thickness data that offers complete coverage of some region, typically with data posted on a uniform grid. While models may operate on subkilometer resolutions, required to accurately model grounding line migration Pattyn et al., 2013), geometry products are derived from radar flight paths, which tend to be at least 5-10 km apart, and so must employ some kind of interpolation. In addition, previous studies have shown that products based on conventional interpolation techniques, such as kriging (e.g., Bedmap2, covering Antarctica) can be inconsistent with the conservation of mass, given observations of velocity, surface mass balance (SMB), and surface elevation change (dh∕dt) . For Pine Island Glacier, a fast-flowing ice stream in the Amundsen Sea Embayment, West Antarctica, Bedmap2, along with velocity observations, produces a thickening tendency of order 100 m yr −1 in the region of the grounding line, which is not observed in the pattern of dh∕dt or SMB (Rignot et al., 2014). As one of the major sources of sea level contributions from the Antarctic ice sheet, this catchment is the focus of our study. In this paper, we model the Pine Island Glacier catchment using two different bed geometries: Bedmap2 and a geometry created using the principle of mass conservation, that is, a mass-conserving initial thickness and bedrock along the lines of Millan et al. (2017), Morlighem et al. (2011), andRignot et al. (2014). We aim to demonstrate the importance of having a consistent representation of the bed, particularly in the vicinity of the grounding line, to modeling short-term ice stream dynamics.

Mass-Conserving Bed
We solve a two-dimensional advection-diffusion problem approximating the principle of mass conservation to obtain a field of ice thickness covering Pine Island Glacier, given dense fields of observed velocity, dh∕dt, and modeled SMB (RACMO2), together with sparse ice thickness data measured along airborne radar flight lines. Further details of the algorithm and data used can be found in the supporting information. In fast-flowing regions, the ice thickness, h mc (mc for mass conserved), is subtracted from the upper ice surface, s, provided by Bedmap2, to give the lower ice surface, b mc . For the grounded ice, b mc is equivalent to the bedrock topography, r mc . In slow-flowing regions, Bedmap2 geometry (h bm , r bm ) is used (with a velocity-weighted transition between the two geometries). Radar thickness measurements from the problematic region near the grounding line-where the data as supplied indicate a submarine ridge but may do so erroneously (Rignot et al., 2014)-are removed, allowing mass conservation to fill in the gap.
To obtain a sub-ice shelf bathymetry, we use r bm , although we make some minor alterations. The ice shelf cavity has been partially surveyed using an Autonomous Underwater Vehicle (AUV) (Jenkins et al., 2010) and information about the bathymetry from this has been incorporated in Bedmap2 (Fretwell et al., 2013). However, the AUV was unable to get close to the grounding line in the northern part of the ice shelf, and in Bedmap2 this region is particularly shallow as it is interpolated between a ridge underneath the ice shelf and the shallow ridge at the grounding line. Therefore, we recompute the bathymetry by interpolating between our new deeper grounding line and the extent of the AUV bathymetry surveys.

Simulations
In order to test the geometry, we run 50 year simulations of Pine Island Glacier for both the Bedmap2 geometry (h bm , r bm ) and our new mass-conserved geometry (h mc , r mc ). We use the BISICLES ice sheet model, a vertically integrated stress model based on Schoof and Hindmarsh (2010) and described by Cornford et al. (2013), which includes simplified vertical shear in the effective viscosity but neglects it in the mass flux. It relies upon adaptive mesh refinement to allow for high-resolution modeling (up to 250 m in the simulations reported here) near the grounding line and shear margins, while avoiding the computational cost of a uniformly fine grid.
First, for each geometry, we initialize BISICLES to find unknown two-dimensional fields of the sliding law coefficient (C) and viscosity stiffening factor ( ) by minimizing an objective function, so as to reduce the misfit between modeled speeds and velocity observations . We then estimate a background melt rate function, M 0 (x, y, t), constructed to hold the ice shelf close to steady state given the present-day thickness and velocity, but evolving so that enhanced melt rates migrate with the grounding line. This method is described in detail by Cornford et al. (2015). For the purpose of this paper, we will briefly describe the representation of the sliding law.
At the base of the ice sheet b, basal traction b (or more formally defined, basal shear stress) is the tangential (shear) component of the stress tensor at the bed | z=b , and can be found using where u b is the tangential velocity to the bed plane and N is the effective pressure, approximated here as follows a sliding law where the ice is grounded on bedrock (with C(x, y) and m being parameters of the chosen law) and is zero where the ice is floating. We test the impact of a number of sliding laws by setting m = 1 for a linear viscous law; m = 1∕3 and m = 1∕9 for nonlinear laws. As well as these traditional power law relationships between velocity and stress, we also test a Coulomb-limited sliding law after Tsai et al. (2015) where = 0.5. This law builds Coulomb friction behavior into the sliding dynamics near the grounding line, where it is thought that a glacial till layer influences the sliding (Tsai et al., 2015). For regions where ice thickness above flotation is small, the Coulomb sliding applies (equation (2), whereas farther upstream, where ice is thicker, the power law applies (equation (1). Here we use m = 1∕3 with the corresponding C(x, y).
The two geometries and four sliding laws give a total of eight configurations. From the initial state, we run each of them forward in time to 50 years. The data sets available for use in the initialization are from varying epochs, and as such the starting point of the forward simulations is imprecise. For example, the Bedmap2 geometry and grounding line position is relevant to 1996, whereas the velocity observations cover 2007-2009 . No additional melt rate forcing beyond that given by M 0 (x, y, t) is applied over the 50 years.

Results
The main difference between r mc and r bm is the removal of a topographic rise in the vicinity of the 1996 grounding line. This feature has a distinctive effect on all of the r bm simulations: an extensive region of negative flux divergence just upstream from the grounding line ( Figure 1a). The reason for this is straightforward: ice flowing through this region experiences little acceleration, so that ∇ ⋅ (uh) ≈ u ⋅ ∇h. Put another way, given that ice is not being stretched as it passes through this region, thicker ice upstream is simply transported downstream to replace thinner ice. By removing suspect ice thickness data here and applying mass conservation, the ice is thicker, and therefore, given the same surface elevation, the bed is deeper. This removes the highly negative signal in the flux divergence in the grounded ice that is at odds with the observed surface elevation change (Figure 1b). A large region of negative flux divergence does remain in the r mc simulations, but it is within the ice shelf, and so can be realistically countered by oceanic melt rates.
Note the difference in the grounding line position between r bm and r mc (Figure 1). This is due to the differences between the point of flotation determined by the upper surface elevation and computed thickness. The sea level contribution from the Pine Island Glacier catchment is greater in the new geometry simulations (h mc , r mc ) than the Bedmap2 simulations (h bm , r bm ), independent of which sliding law is used. Table 1 lists the mean rate of mass above flotation lost over the course of the 50 year simulations. Simulations carried out with (h bm , r bm ) place their mean rates in the range 0.04-0.07 (mm yr −1 sea level equivalent), compared to 0.08-0.10 (mm yr −1 sea level equivalent) with (h mc , r mc ). Table 1 also shows that while mass loss decreases as m → 0 in the Bedmap2 simulations, that is, as the sliding law tends toward a yield stress rule, mass loss increases as m → 0 with the mass-conserving bed.
As well as a greater mean rate of mass loss, the mass-conserving bed produces a qualitatively different time varying signal. Figure 2b shows that the rate of sea level contribution tends to decrease over the course of the r bm simulations. This is in contrast to the power sliding law simulations using r mc , which experience a sharp increase in the rate at ∼30 years. This acceleration lasts for approximately 5 years before the rate slowly subsides again, although it does not entirely decay to previous rates within the simulated period.
The r mc Coulomb sliding law simulation has an early acceleration in the rate of sea level rise, peaking at ∼10 years, after which the rate generally decreases for the rest of the simulation. The resulting total sea level contribution is at least 20% greater for the Coulomb law simulation, using r mc , compared to the other sliding laws. For the r bm simulations, while the Coulomb sliding law does contribute more to sea level than other sliding laws, the difference between them is less pronounced.
Grounding line migration occurs in all simulations (Figure 3), with the particular geometry determining the extent and pattern of retreat. There is more grounding line retreat in the r mc simulations than the r bm simulations. In the r bm simulations, the central portion of the fast-flowing trunk quickly ungrounds from the protruding topographic rise (Figure 3c). However, the grounding line then remains in a stable position for the rest of the simulation. In contrast, in the r mc simulations, the central portion of the grounding line progressively retreats throughout the 50 years to a distance approximately 8 km upstream of the initial position (Figure 3b).
The different sliding laws tend to have little influence on the resulting pattern of retreat. That said, in the Coulomb law simulations, the portions of the grounding line peripheral to the main channel tend to experience more retreat than the power sliding laws.

Discussion
The motivation for creating a new geometry of Pine Island Glacier using mass conservation was to remove the apparently erroneous "grounding line wedge" present in Bedmap2. This has been achieved (Figures 3b  and 3c), and as a result, we have also removed the thickening tendency near the grounding line (Figure 1). A likely explanation for the erroneous thickness measurements, used by Bedmap2, is that strong returns from water-filled basal crevasses were picked out in the radargrams, rather than the "true" base of the ice (Rignot et al., 2014). These crevasses occur when the ice is near flotation, in the transition zone between two stress regimes across the grounding line (Jacobel et al., 2014). The thicker ice found here in our new geometry, compared to Bedmap2, is consistent with the assumption that the interpreted basal radar return is a minimum estimate of ice thickness, but in reality the thickness is likely greater. The ice surface in the grounding zone is heavily crevassed, which impedes access for ground truth measurements. The direct results of the mass conservation scheme, namely, the thickness and bedrock maps, are broadly similar to other studies that have used mass conservation to find ice thickness. For example, Rignot et al. (2014), and more recently Millan et al. (2017), used the optimization method of Morlighem et al. (2011), which produces a bed topography 150 to 350 m deeper than Bedmap2 near the 1996 grounding line. Nias et al. (2016) used a relaxation scheme, where the ice thickness was altered iteratively during the model initialization process, in an attempt to reduce the flux divergence anomaly. This also produced a bed topography up to 215 m deeper than Bedmap2 in the region of the grounding line wedge. However, elsewhere in the trunk of Pine Island Glacier, it was often shallower than Bedmap2, by 100 to 200 m. Given the assumption that the true topography is at least as deep as that indicated by radar measurements, due to englacial reflections, this seems less likely than the new topography generated here. While the results of these studies exhibit some similar features to our geometry compared to Bedmap2, we have gone on to assess the impact of the geometry error on modeled ice stream dynamics.
The modification to the geometry affects the dynamical response of Pine Island Glacier on decadal timescales. The difference in model results between the two bed topographies used is due to the presence, in the case of r bm , or absence (r mc ), of the topographic rise near the initial grounding line. The grounding line in the r bm simulations does not retreat off this "grounding line wedge" over the course of the 50 year simulations (Figure 3c), which is at odds with observations (Rignot et al., 2014). No such obstacle to retreat exists in the r mc simulations, and so the grounding line progressively retreats throughout the simulations (Figure 3b), which is more consistent with the observed motion (Rignot et al., 2014).
The retreat over r mc is not uniform in time. Between 20 and 30 years, the grounding line jumps upstream by approximately 4 km (Figure 3b). This is reflected in the surface elevation change. For example, Figure 4 demonstrates that for the m = 1∕3 simulation, this period of more rapid retreat coincides with a short period of enhanced thinning, concentrated near the grounding line and propagating upstream. A similar pattern, albeit with slightly different onset times, occurs for the other power sliding laws.
Enhanced thinning coincides with the step increase in the rate of sea level contribution (Figure 2b). The maximum thinning rates persists for just one year, before subsiding to previous rates. However, the rate of sea level contribution is not so quick to recover, particularly for the m = 1∕3 and m = 1∕9 power sliding laws, where the high rates of sea level rise are sustained for the remaining duration of the simulations (Figure 2b).
Similar short-term variability in thinning rates have been observed in satellite data from Pine Island Glacier over the last couple of decades. Konrad et al. (2017) collated data from numerous satellite altimeter missions to give a continuous record of thinning rates in Pine Island Glacier, and its neighbors. From ∼2002 to 2010 there was an acceleration in the thinning rates closest to the grounding line, to a maximum of ∼6 m yr −1 , in 2007. They found that the enhanced thinning rates rapidly propagated upstream, at a rate of 10-12 km yr −1 along the central flow lines of Pine Island Glacier. In 2010 the observed thinning rates near the grounding line dropped suddenly back to previous rates of ∼1 m yr −1 .
These observations, combined with our model results, indicate a mechanism for internal variability within Pine Island Glacier. Previous work shows that grounding line retreat can be sensitive to relatively low-amplitude topographic features (Durand et al., 2011;Nias et al., 2016;Sun et al., 2014). Here we find that the retreat at 30 year coincides with movement of the grounding line across a small circular depression in the bed, only approximately 20-30 m in amplitude and 4 km across. The localized ungrounding of this depression results in strong dynamic thinning that propagates upstream. While the thinning rates return to their previous magnitudes, once the grounding line stabilizes again on the locally prograde slope on the upstream side of the depression, the rate of sea level contributions does not. Sun et al. (2014) found that the height of the grounding line wedge in Bedmap2 only needs to change on the order of 10 m to alter the timing of the onset of retreat off it by decades. They suggest that once the retreat off this ridge does occur, grounding line retreat is then less sensitive to noise in the bed topography. However, our r mc results show that on decadal timescales, small features can impact on grounding line retreat, thinning rates and sea level contributions. Perhaps averaged over longer timescales, the influence of such features diminishes.
Whether the shallow depression that causes the signal in our r mc simulations is actually present in reality, we cannot confirm without further measurements, although the application of mass conservation in our method implies that it is. Despite this uncertainty, our results demonstrate that the subtlest of features can have a significant impact on short-term variability in thinning rates, which can have a longer term impact on rates of sea level rise (Figure 2b). This demonstrates the importance of having an accurate representation of the geometry when modeling ice streams.
Our simulations show that accurately representing the geometry of the ice sheet is at least as important as the choice of sliding law. For the two geometries, the relationships between power law and mass loss are the inverse of one another: as the sliding law becomes more plastic (m → 0), mass loss from the r mc simulations increases, whereas it decreases for r bm . Joughin et al. (2010) tested a range of different sliding law behaviors, with m → 0, in the Pine Island Glacier catchment. They found that the observed speedup at the grounding line 1996-2009 was best reproduced by a mixed model for sliding, with m = 1∕3 in hard-bedded areas and plastic behavior (m → 0) where till is likely, and a model with plastic conditions across the ice stream. The linear viscous law (m = 1) was unable to reproduce the instantaneous response to an ungrounding. The results of the studies discussed above suggest that our simulations using m = 1∕9 are more likely to represent true sliding behavior than the m = 1 or 1∕3 simulations. This is supported by a series of seismic surveys at multiple sites on the ice stream and its tributaries, which suggest that a weak till layer is ubiquitous across Pine Island Glacier (Brisbourne et al., 2017). This has implications for mass loss, with the two m = 1∕9 simulations resulting in the greatest disparity between the two geometries, with the mean rate of sea level contributions for r mc being more than double that of r bm (Table 1).
The effect of sliding law exponent on stability is a profoundly local one, depending on bedrock geometry. For example, neighboring Thwaites Glacier has been shown to stabilize as m → 0 (Parizek et al., 2013), similarly to our r bm simulations. More plastic beds tend to quickly transfer changes at the grounding line up glacier. In our r mc simulations, as m → 0, mass loss across the grounding line is enhanced, as thinning propagates more quickly up glacier (thinning rates of 2 m yr −1 spreads up the central trunk at a rate of 6.7 km yr −1 for m = 1, compared to 9.1 km yr −1 for m = 1∕9), whereas the stabilizing effect of the grounding line wedge in r bm is more effective as m → 0.
The Coulomb sliding law, in which basal shear stress is proportional to the effective pressure near the grounding line, removes the discontinuity in basal traction between the ice stream and ice shelf. This results in a greater sensitivity to perturbations compared to traditional power laws (Brondex et al., 2017;Tsai et al., 2015), which is seen here, with the resulting mass loss, for either bed, being greater than that of the power law simulations. This effect is more pronounced for r mc , owing to the more retrograde slope at the grounding line. Even the Coulomb sliding law for r bm does not produce as much mass loss as any of the r mc simulations, again demonstrating the importance of having an accurate representation of ice sheet geometry.

Conclusion
We have modeled decadal-scale variability in the dynamics of Pine Island Glacier using maps of ice thickness and bedrock elevation created using the principle of mass conservation. Compared to simulations using Bedmap2 geometry, which are unable to produce observed grounding line retreat from the 1996 position, these simulations are an improvement. Subtle variation in the geometry near the grounding line can trigger a response in the ice sheet that is felt hundreds of kilometers upstream. This is demonstrated by the thinning signal caused by retreat over a low-amplitude feature in the bed topography.
These results highlight the importance of using accurate and highly resolved bed topography in marine ice sheet modeling. The impact of the bed topography is at least as great as the influence of the particular form of the sliding law used. As m → 0, the difference in simulated mass loss between the two geometries increases.