Marine ice regulates the future stability of a large Antarctic ice shelf

The collapses of the Larsen A and B ice shelves on the Antarctic Peninsula in 1995 and 2002 confirm the impact of southward-propagating climate warming in this region. Recent mass and dynamic changes of Larsen B’s southern neighbour Larsen C, the fourth largest ice shelf in Antarctica, may herald a similar instability. Here, using a validated ice-shelf model run in diagnostic mode, constrained by satellite and in situ geophysical data, we identify the nature of this potential instability. We demonstrate that the present-day spatial distribution and orientation of the principal stresses within Larsen C ice shelf are akin to those within pre-collapse Larsen B. When Larsen B’s stabilizing frontal portion was lost in 1995, the unstable remaining shelf accelerated, crumbled and ultimately collapsed. We hypothesize that Larsen C ice shelf may suffer a similar fate if it were not stabilized by warm and mechanically soft marine ice, entrained within narrow suture zones.

T he stability of an Antarctic ice shelf depends on the balance between the constructive stresses that ensure its integrity and the destructive stresses that compromise it. The spatial distribution and temporal evolution of the first and second principal stresses on Antarctic ice shelves depend on the geometry of the embayment in which the shelf is situated, the spatial and temporal attributes of the ice recharge it receives from feeding glaciers, environmental factors such as the evolution of atmospheric surface and oceanographic basal conditions 1 , and the physical properties of the types of ice of which it is composed 2,3 . On many Antarctic ice shelves, including Larsen C (Figs 1 and 2), tensile stresses dominate over compressional stresses. The normal stress with the highest magnitude, the first principal stress, therefore commonly controls the magnitude and direction of local strain rate and ice-shelf extension, and by implication in uniform strain-rate fields therefore the direction of extensional flow.
Fractures, such as horizontally extending rifts and vertically and horizontally extending surface and basal crevasses can be inherited from feeding glaciers. Fractures may also originate near the grounding line as the ice begins to float, or open up on the shelf in response to bending or shearing forces or thermal processes 4,5 . Such fractures commonly strike orthogonal to ice flow, as typified, for example, on Larsen C ice shelf where several cross-flow zones of major and minor fracturing are present (Fig. 1). Within such zones, rifts and crevasses often form regular sequences that share the same formative conditions 3 , such as, for example, the major rift sequences 'R1' to 'R5' (Fig. 1). Only on rare occasions do fractures intersect the ice-shelf's flow lines at oblique angles, such as in the downflow sectors of sequences R1 and R5 (Fig. 1). If the tensile first principal stress is aligned perpendicular (parallel) to a fracture's strike then that fracture's opening rate is maximized (minimized). As fractures tend to strike perpendicular to flow, we propose that the shelf-wide distribution of the angles between the flow direction and the first principal stress (hereafter 'stress-flow angles') can serve as a first-order criterion on which to judge an ice-shelf's stability. Fundamental to a quantitative prediction of Larsen C's stability in a warming climate is therefore the assessment of the distribution and climatic sensitivity of the stress-flow angles across the ice shelf, and it is this factor that we address in this manuscript.
We adopt our previous approach 4 and model ice-shelf velocities, flow lines and stresses with a continuum-mechanical ice-flow model with spatially uniform ice properties. In subsequent sensitivity studies, we ascertain the impact of spatially variable ice properties and changing model boundary conditions on modelled stress fields and stress-flow angles. More specifically, we follow a three-stage line of argument to demonstrate that marine ice stabilizes Larsen C ice shelf. First, the distribution of stress-flow angles across Larsen C ice shelf is analysed to demonstrate that its frontal portion provides critical restraint for the otherwise unstable central portion. Second, we infer that the rapid disintegration of Larsen B ice shelf between 1995 and 2002 might serve as a plausible blueprint for Larsen C's future demise, because Larsen B was characterized by an equivalent distribution of stress-flow angles pre-collapse. Following a brief introduction to marine ice-bearing suture zones and their anomalous physical properties, we argue, third, that such zones presently prevent the loss of Larsen C's frontal portion and by implication, therefore, stabilize the whole ice shelf.

Results
Stress-flow angles on Larsen C ice shelf. Our first stage of argument considers that both the first and the second principal stresses have highest absolute magnitudes near Larsen C's grounding line and decrease non-uniformly towards the calving front (Fig. 2). Stress-flow angles tend towards 90°near Larsen C's grounding line (Fig. 2a) where feeding glaciers accelerate and spread out laterally as the ice begins to float and basal drag is removed. In the central portion of the ice shelf, the first principal stress aligns with ice-shelf flow (Fig. 2c), and the stress-flow angles thus approach zero (Fig. 2a). Downstream of the embayment, stress-flow angles once more tend towards 90°as the ice shelf spreads out laterally (Fig. 2a,b). Here the tensile first principal stress is oriented along the ice front and parallel to rift zones R1 to R3 (Figs 1 and 2a), favouring a stable ice shelf. Consistently, Larsen C's calving style is presently dominated by infrequent detachment of large tabular icebergs, as demonstrated by the most recent large calving event in 2008 (Fig. 1). With the removal of the stabilizing frontal portion, the largest tensile stress would be oriented perpendicular to the calving front and to the surface and basal crevasses 5-7 that are numerous in this area, thus destabilizing them. Larsen C's frontal portion therefore provides essential restraint for the shelf's central portion that may otherwise be unstable. The transition from compressive to tensile second principal stresses in Larsen C's central to southern frontal portion defines a 'compressive arch' 8 (Fig. 2d)  In this inset, the dark blue stripe overlain on a high-resolution Landsat image marks the area of marine-ice accretion. The newly formed suture zone becomes laterally compacted by the neighbouring converging Trail-Inlet (TI) and Solberg-Inlet (SI) flow units. The lighter purple stripe enclosed by the dark blue stripes traces the meteoric ice contribution from a small glacier on the Joerg Peninsula, as revealed by our ground-penetrating radar data (Fig. 4).
extending from the Kenyon Peninsula in the south towards the Bawden ice rise in the north (Figs 1 and 2a). It was proposed previously that ice-shelf retreat beyond a critical arch may result in rapid disintegration 8 . Because the transition from low to high stress-flow angles is located closer to the calving front than the compressive arch (Fig. 2a), a hypothetical retreating Larsen C ice shelf may become unstable well before the compressive arch is breached.
Stress-flow angles on Larsen B ice shelf. Our second stage of argument considers that the distribution of stress-flow angles on Larsen B ice shelf in 1986 before its collapse (1995-2002, Fig. 3) was similar to those of present-day Larsen C (Fig. 2a), and the stress field in Larsen B's frontal portion was likewise marked by a compressive arch 8 . Larsen B's evolution towards rapid disintegration in 2002 (ref. 1) may therefore offer insights into Larsen C's future stability. Comparable patterns of first principal stresses perpendicular to the flow direction at the grounding line and calving margin, and parallel to flow in the centre, were present before 1995 (ref. 8), and Larsen B calved large tabular icebergs just as Larsen C currently does. In 1995, Larsen B's entire frontal portion calved away (top inset in Fig. 3), resulting in nearzero stress-flow angles at the new calving front and fractureorthogonal tensile stresses that then encouraged the propagation of existing rifts and crevasses. Larsen B's load-bearing capacity subsequently decreased 9 , causing the shelf to crumble by frequent calving of small icebergs and its calving front to continue to retreat ( Fig. 3), concluding in its eventual collapse in 2002 (refs 9,10). Owing to the similarity of shelf-wide stress-angle distributions, Larsen C might similarly disintegrate if its frontal portion, and thus the restraint it provides for the shelf's central portion, is lost. Ongoing preservation of Larsen C's frontal portion is therefore necessary for its stability.
Sensitivity of model outputs to boundary conditions. Larsen C ice shelf has experienced recent mass and dynamic changes that are particularly pronounced in its northern part. These include ice-shelf acceleration 11 , surface lowering due to melt-driven firn compaction 11-14 and thinning of feeding glaciers 15 . We have therefore conducted a series of perturbation experiments with our continuum-mechanical ice-flow model to ascertain the sensitivity of Larsen C's velocity and stress fields to hypothetical changes in the ice-shelf's calving front geometry, the inflow velocities of the feeding glaciers and the thickness of ice shelf ( Supplementary  Fig. 1). A comparison of the first and the second principal stress fields before and after a major calving event that occurred between 2002 and 2008 (see superimposed calving front in Fig. 1) demonstrates only weak sensitivity of Larsen C's principal stress  ARTICLE fields to simulated changes in the geometry of its calving front ( Supplementary Fig. 1). A 20% acceleration of the feeding glaciers would increase the mean velocity of the ice shelf by 12% from 362 m a À 1 to 403 m a À 1 . However, because the spatial velocity gradients across Larsen C ice shelf, and thus the strain rates, remain largely unaffected, the first and second principal stress fields have comparable magnitudes and directions with and without acceleration ( Supplementary Fig. 1). If Larsen C experienced firn compaction or basal melting that spatially averaged 20 m, its mean velocity would decrease byB7% from 362 m a À 1 to 338 m a À 1 . However, once again the first and second principal stress fields with and without thinning have comparable magnitudes (Supplementary Fig. 1).
Despite the recent mass and dynamic changes that Larsen C ice shelf has been experiencing 11-16 , our model sensitivity tests thus reveal that the integrity of Larsen C's stabilizing frontal portion is unlikely to be compromised by mass and dynamic changes in the foreseeable future. Following a brief introduction to ice-shelf suture zones and their anomalous mechanical properties, our third stage of argument considers instead that marine ice-bearing suture zones 2,3 preserve this portion because they prevent rifts from propagating laterally across and coalescing within it; and by implication therefore stabilize the whole of Larsen C ice shelf.
Marine ice in the Larsen C ice shelf. Larsen C, like most other Antarctic Peninsula ice shelves, is principally composed of flowparallel units of meteoric ice that are sustained by feeding glaciers and snow accumulation 1 (Fig. 1), and narrow interstitial suture zones. Suture zones are partially composed of marine ice and commonly appear as smooth flow-parallel bands in satellite imagery 2 (Fig. 1). Prominent suture zones on Larsen C include those originating leeward of the Joerg Peninsula ('J' and red stripe in Fig. 1), Tonkin and Francis Islands (respectively 'TO' and 'FI' in Fig. 1) in the south, and Churchill Peninsula ('C' in Fig. 1) in the north. These zones serve to isolate the prominent areas of fracturing ('R1' to 'R5' in Fig. 1) in the ice-shelf's frontal portion. The presence of marine ice within these suture zones was revealed by airborne radio-echo sounding, substantiated by mathematical modelling of sub-shelf freeze-on (ref. 2) and, within Joerg Peninsula suture zone, delineated at high spatial resolution by our ground-penetrating radar surveys (GPR) undertaken in the 2008/ 09 and 2009/10 austral summers (Fig. 4). Our GPR profiles delineate the base of the meteoric Trail-Inlet and Solberg-Inlet flow units, but cannot detect the suture zone's base 2 (Fig. 4a,b). We therefore used seismic reflection data acquired at P1 (ref. 4) (2008-09; Fig. 4a) and calculations of ice-shelf draft at P2 (refs 4,12) (2009-10; Fig. 4b) to delineate the base of the marine ice within the Joerg Peninsula suture zone. Marine-ice bodies have a temperature similar to the sub-shelf ocean waters from which they are formed ( À 1.5°C to À 2°C), and are therefore anomalously soft 17 . In contrast, meteoric ice-shelf units are much colder because they are derived from feeding glaciers and snow accumulation subject to annually averaged surface temperatures of À 15°C and below. Because the thermal diffusivity of ice is very small, the marked contrast in meteoric versus marine ice temperatures is expected to persist along the entire length of an ice shelf 18 . Warmer marine ice deforms more readily under the same long-term stress loading, imposed upon the ice shelf by gravity-driven flow, than colder meteoric ice 12,19 . A larger proportion of that long-term loading is therefore available to drive elastic fracture in meteoric than in marine ice,   so that warmer marine ice-bearing suture zones are less prone to elastic fracture than surrounding, colder meteoric ice units 20,21 .
Stress-flow angles from observed data. Interferometric Synthetic Aperture Radar (InSAR)-derived flow velocities of Larsen C ice shelf 22 allow, initially, calculation of the spatial distributions of ice-shelf strain rates. Subsequently, the principal stresses and stress-flow angles at the ice-shelf surface are calculated by assuming that the temperature of all shelf ice is equal to Larsen C's annually averaged surface temperature of À 15°C ( Supplementary Fig. 2). Stresses will be higher in colder meteoric ice than in warmer marine ice at a given strain rate, so that marine ice is less likely to fracture 20,21,23 . The assumption of a spatially invariant temperature of À 15°C is consequently violated in major rift zones filled with warmer ice mélange, where strain rates are therefore anomalously high and stress magnitudes overestimated by our calculations (Supplementary Fig. 2). The assumption is also violated in the suture zones containing warm bodies of marine ice, although here shear stress bridging between neighbouring meteoric flow units minimizes the build-up of anomalously high strain rates, so that stress magnitudes are not normally overestimated. Despite these potential limitations, a comparison of Fig. 2a and Supplementary Fig. 2 reveals good overall agreement between the spatial patterns of principal stresses and stress-flow angles reconstructed from modelled and observed data. However, the patterns inferred from observed data ( Supplementary Fig. 2) are relatively noisy and lack prognostic capabilities 4 . We therefore prefer to adopt our previous approach 4 that focuses on analyses of ice-shelf flow velocities modelled with a continuum-mechanical model constrained by spatially uniform ice properties (Fig. 2).

Discussion
Our model sensitivity experiments with simulated softening of marine ice-bearing suture zones and major rifts on Larsen C ice shelf (Fig. 5) reveal that concurrent changes in strain rates within and adjacent to the suture zones act to rotate the trajectories of the first and second principal stresses, respectively, into floworthogonal and flow-parallel directions. Inferred stress rotations have two important consequences. First, higher stress-flow angles tend to be focused within marine ice-bearing suture zones, highlighting their rift-stabilizing potential (Fig. 5). Second, the simulated temperature enhancements within suture and rift zones combine to push the compressive arch landwards, especially in the ice-shelf's northeastern sector, which reduces the likelihood of ice-shelf retreat beyond a critical compressive arch. Together, these two consequences demonstrate that marine ice-bearing suture zones have a critical role in preventing the loss of Larsen C's frontal portion. Suture zones on Larsen C ice shelf tend to become thicker and narrower as they are advected downstream, owing to lateral compaction and the balance between continued accretion and basal melting 2 , as demonstrated by our measurements from the Joerg Peninsula suture zone (Fig. 4). Climatically controlled future ocean warming may accelerate upstream marine ice accretion 2,24,25 , but simultaneously speed-up downstream melting. Melt-induced size reduction of marine ice bodies would decrease their relative influence within a given vertical ice column, reducing that column's ability to resist fracture 26 . Oceanographic modelling is consistent with the presence of two prominent marine ice-bearing suture zones within Larsen B ice shelf before its collapse, originating downflow of Foyn Peninsula and Cape Disappointment 2 . In the absence of direct observations, we speculate accordingly that these zones were either too thin to resist rift propagation or had been weakened over time due to prolonged basal melting. Larsen C's northeastern calving front is placed into specific focus because the stress-flow angles in the northeast already approach zero (Fig. 2a). Here the ice shelf is pinned on the Bawden ice rise, stabilizing the large Churchill Peninsula suture zone and rift zones R4 and R5 (Figs 1 and 2a). Retreat from this ice rise due to suture-and rift-zone weakening would cause Larsen C's northeastern sector to accelerate 27 , and the areal extent of near-zero stress-flow angles to expand towards the south (compare Fig. 2 and Supplementary Fig. 3). Because these consequences would contribute to the de-stabilization of the ice-shelf's frontal portion, future work should identify the role of the ocean in marine ice accretion and in promoting basal melting 13,28,29 .
The presence of variably sized bodies of marine ice has already been inferred for several other Antarctic ice shelves, including the three largest, Filchner-Ronne 30 , Ross 31 and Amery 18,25 . Suture zones within these shelves are also readily identifiable from satellite imagery. Our findings thus highlight the need to map marine ice-bearing suture zones within, determine the oceanographic conditions beneath, and model the evolving stress  regimes of the many Antarctic ice shelves currently experiencing enhanced erosion by ocean warming. A fundamental question that remains as yet unanswered is why Larsen B's entire frontal portion calved away in 1995, presenting the initial trigger for its demise. We recommend therefore that future work should also focus on monitoring the evolution of Antarctic ice shelves' calving fronts in relation to indicative changes in ice flow velocities and crevasse growths. Less feasible at present because of a lack of readily deployed observational techniques, but particularly diagnostic, would be the monitoring of associated changes in the geometry of marine ice inclusions in the shelves' suture zones. Ultimately marine ice-bearing suture zones contribute to regulating Antarctica's future contribution to sealevel rise and must therefore be parameterized in ice-sheet models.

Methods
Model description. Our continuum-mechanical ice-shelf model 32 is implemented numerically using a finite-difference scheme with 1.25 Â 1.25 km grid cells. The model accounts for gravitational driving forces and associated mechanical stresses, and is based on the shallow-shelf approximation. As such it assumes hydrostatic equilibrium and depth-invariant horizontal flow velocities, and neglects friction at the ice-ocean interface and vertical shear-strain changes because of bending forces.
As input, the model requires specification of three-dimensional ice-shelf geometry and the temperature-dependent flow-rate factor B, as well as the inflow velocities at the grounding line boundary of the model domain. The model equations are solved numerically using a force-balance approach 32 , generating as outputs ice-shelf flow velocities and the two-dimensional stress field that is subsequently projected into the principal directions. We previously used the same model implementation to assert the present stability of the Larsen C ice shelf based on the shelf's geometry in 2002 (ref. 4) before a large tabular calving event (Fig. 1). All model runs reported here are previously unpublished, and focus on Larsen C's up-to-date post-calving geometry in 2008 (Figs 1 and 2). The 1.25 Â 1.25 km size of our grid cells is much larger than the width of individual fractures, but smaller than the areal extent of most major fracture zones.
Model implementation of marine ice. Our perturbation experiments implemented the most prominent marine ice-bearing suture zones and rifts on Larsen C ice shelf for the purpose of sensitivity testing (Fig. 5). Adopting previous work 4 , we assumed that meteoric ice units have a mean temperature of À 12°C, increasing from À 15°C at the ice-shelf surface to À 2°C at the base. Our implementation assumes a relatively elevated mean suture-zone ice temperature of À 7°C. Here, marine ice with a temperature of À 2°C represents the bottom half of total suturezone thickness (Fig. 4a), and follows the meteoric ice temperature profile in the suture zone's top half. These values of meteoric and marine ice temperatures are vertically averaged to yield the flow-rate factor B for any given vertical ice column in our 1.25 Â 1.25 km model grid. Our implementation further assumes that the rifts in Larsen C's south-eastern sector (Fig. 5) are filled with an ice mélange 24 that has a vertically constant temperature of À 2°C between the ice-shelf surface and its base.
Acquisition and processing of GPR data. The GPR data presented in Fig. 4  MHz antennas, operated in common offset mode using a snow-scooter-towed sledge assembly. Choice of the 50 MHz antennas provided the best compromise between depth penetration and vertical resolution, following a suite of testing using several antennas of different frequencies in 2008-09. The GPR data were acquired with a 0.8-ns sampling interval, where each trace was acquired eight times and then stacked into a single recorded trace.
A compromise between towing speed as well as required window length and number of stacks for each trace implied that one trace was recorded Bevery 3-4 m along the survey lines. In 2008-09, a handheld GPS was linked to the GPR system, locating each individual recorded radar trace with a planimetric precision of approximately±5 m. In 2009-10, precise (±0.1 m) planimetric and height location of the antennas was achieved with a differential Leica System 1200 GPS. The raw radar data were processed using the commercial Reflex-W package. Standard techniques were applied, including de-wow, automatic gain control and band-pass filter, as well as adjustment for noticeable surface topography in 2009-10 immediately downstream of the Joerg Peninsula (Figs 1 and 4). The data are presented as un-migrated profiles because the meteoric-marine ice interfaces are more readily distinguishable in this format. Travel times were converted to depth assuming a depth-averaged radar velocity of 0.175 m ns À 1 based on commonmidpoint surveys conducted in both field seasons 4,5 . The use of the older and lower-power PE-100 GPR system in 2009-10 resulted in much noisier data ( Fig. 4b) than those collected in 2008-09 (Fig. 4a), although marine ice bodies are still readily distinguishable (Fig. 4b).
Determination of suture-zone thicknesses. Our seismic reflection data were acquired in walk-away mode in 2008-09 (ref. 4) using a Geometrics-Geode-based system and explosive shots deployed in shallow drill holes, yielding an estimate of the suture-zone thickness at P1 (Figs 1 and 4a). No seismic data were available in 2009-10. Instead the suture-zone thickness at P2 (Figs 1 and 4b) was calculated using (a) the precise differential GPS data available for P2, under the assumption of hydrostatic equilibrium and a mean density profile derived from the 2008/09 seismic data 4 ; and (b) a recently published method for the estimation of firn-air content on Larsen C 12 . Results from both methods agree to within 0.5%.