Present stability of the Larsen C ice shelf, Antarctic Peninsula

Abstract

We modelled the flow of the Larsen C and northernmost Larsen D ice shelves, Antarctic Peninsula, using a model of continuum mechanics of ice flow, and applied a fracture criterion to the simulated velocities to investigate the ice shelf’s present-day stability. Constraints come from satellite data and geophysical measurements from the 2008/09 austral summer. Ice-shelf thickness was derived from BEDMAP and ICESat data, and the density–depth relationship was inferred from our in situ seismic reflection data. We obtained excellent agreements between modelled and measured ice-flow velocities, and inferred and observed distributions of rifts and crevasses. Residual discrepancies between regions of predicted fracture and observed crevasses are concentrated in zones where we assume a significant amount of marine ice and therefore altered mechanical properties in the ice column. This emphasizes the importance of these zones and shows that more data are needed to understand their influence on ice-shelf stability. Modelled flow velocities and the corresponding stress distribution indicate that the Larsen C ice shelf is stable at the moment. However, weakening of the elongated marine ice zones could lead to acceleration of the ice shelf due to decoupling from the slower parts in the northern inlets and south of Kenyon Peninsula, leading to a velocity distribution similar to that in the Larsen B ice shelf prior to its disintegration.

Introduction

The recent disintegration of ice shelves on the Antarctic Peninsula is closely linked to climatically controlled southward migration of the −9°C isotherm (Mercer, 1978; Vaughan and others, 2003). The retreat pattern of ice shelves is governed by the overall ice-shelf mass balance including processes such as calving (Skvarca, 1994), disintegration (e.g. Rott and others, 1996), surface and basal melting (Scambos and others, 2000) and also the structure (Glasser and Scambos, 2008; Glasser and others, 2009) and composition of the ice (see Cook and Vaughan, 2009). The Larsen C ice shelf has been thinning over the last two decades (Shepherd and others, 2003), but otherwise does not currently exhibit obvious climatically related signs of retreat (Glasser and others, 2009). However, the −9°C isotherm is presently encroaching upon the ice shelf, which might affect its stability.

Here we present a model of the flow, fracture and stability of peninsular ice shelves. We combine physically based ice-shelf flow modelling (Sandhäger and others, 2005) and linear elastic fracture mechanics (Rist and others, 1996, 1999, 2002), as constrained by satellite and in situ geophysical data, to place quantitative constraints on the present stability of the Larsen C and northernmost Larsen D ice shelves. Although Vieli and others (2006, 2007) and Khazendar and others (2007, 2009) achieved excellent results for ice-shelf rheology with modelling constrained by inversion of satellite-derived velocities, we chose the Sandhäger approach to be more independent from measured velocity data, especially when modelling possible future scenarios. We pay specific attention to the recent hypothesis that mechanical heterogeneity acts to stabilize the Larsen C ice shelf (e.g. Glasser and others, 2009; Holland and others, 2009), but which as yet lacks a quantitative basis. We use as a principal measure of ice-shelf stability the spatial distribution of stress intensity, identifying ice-shelf regions where crevassing or rifting is relatively less or more likely to occur. Mechanically heterogeneous regions, such as the flowbands that originate down-flow of promontories and are subject to basal accretion, would act to reduce local stress intensities and decrease rates of rift propagation, or indeed halt it altogether.

Our modelled flow and fracture distributions are validated against satellite and in situ geophysical field data. These provide a sound basis for future model sensitivity studies of the evolution of the Larsen C ice shelf in a warming climate. Finally, by drawing an analogy to the dynamic evolution of the Larsen B ice shelf prior to its collapse, we comment on the current stability of the Larsen C ice shelf.

Study Area

With an area of approximately 51 000 km² (Cook and Vaughan, 2009), the Larsen C ice shelf is by far the largest on the Antarctic Peninsula (Fig. 1). Located on the eastern side of Graham Land, it is confined by the Jason Peninsula in the north (∼66.5° S) and extends to the Kenyon Peninsula and the Gipps Ice Rise in the south (∼68.5° S). The main ice-shelf flow units originate from Cabinet Inlet, Mill Inlet and Whirlwind Inlet, as well as from the Mobil Oil Inlet, which is characterized by the highest inflow velocities (Fig. 1; Cook and Vaughan, 2009; Glasser and others, 2009). We chose the 2002 configurations of the Larsen C and northernmost Larsen D ice shelves (e.g. Fig. 1) for our modelling study, since the most appropriate satellite datasets were available for that period. The data provide high-quality constraints on our boundary conditions and serve to validate our model outputs.

Fig. 1. MODIS Mosaic of Antarctica (MOA) image of the Larsen C ice shelf (http://nsidc.org/data/nsidc-0280.html) illustrating the Larsen C ice shelf model domain. The two dots indicate the position of the GPS sensors within the 2008/09 field season. They also represent the position of the two crossings of seismic profiles measured to infer density–depth profiles.

Several characteristics of the northeastern and southeastern sectors of the Larsen C ice shelf are particularly relevant to our study. In the northeast, Larsen C is largely sustained by surface accumulation, and numerous smaller inlets and glaciers feed the regions around the Jason Peninsula (Fig. 1; Skvarca, 1994). The ice tends to be heavily crevassed in the transition zones between such feeders and the floating ice shelf (Glasser and others, 2009). The easternmost pinning point of Larsen C is Bawden Ice Rise, generating crevasses that propagate southwards into the ice-shelf body and obliquely to flow, thus delineating nascent icebergs (Fig. 1). In the southeast, an imaginary line, drawn between the eastern tip of the Kenyon Peninsula and Gipps Ice Rise, conceptually separates the Larsen C ice shelf to the north from the Larsen D ice shelf to the south. The entire region is heavily fractured, as manifested in an extensive series of parallel rifts protruding ∼50 km northwards into Larsen C ice shelf and a melange-filled embayment to the south that constitutes the northernmost portion of the Larsen D ice shelf (Cook and Vaughan, 2009). The northernmost portion of the Larsen D ice shelf, and particularly the region between Gipps Ice Rise and Hearst Island, is fed by three major glaciers (Cronus, Casey and Lurabee) and though the ice is heavily crevassed and heterogeneous in its composition, remote-sensing velocities indicate that its flow behaviour might be approximated with a continuum-mechanical approach, similar to the Brunt Ice Shelf (Humbert and others, 2009). Although geographically assigned to Larsen D, we include this part in our model domain in order to investigate the stress distribution of the region between the Kenyon Peninsula and Gipps Ice Rise, which otherwise would have been the domain boundary.

Methods

Model description

The applied continuum-mechanical ice-shelf flow model (Sandhäger and others, 2000, 2005; Grosfeld and Sandhäger, 2004) is based on a finite-difference numerical implementation which simulates progressive ice-shelf evolution as controlled by ice dynamics and variable environmental boundary conditions. It was applied successfully to the Larsen B ice shelf (Sandhäger and others 2005), where simulated velocities were validated by in situ stake measurements. It is beyond the scope of this paper to review in detail the mathematical basis of this ice-shelf model. It is sufficient to emphasize that gravitational driving forces and associated stresses are implemented, while friction at the ice-shelf–ocean boundary and vertical shear strain due to bending forces are neglected (following the approach of MacAyeal and others, 1986). The ice-shelf body is assumed to be in hydrostatic equilibrium, and the horizontal flow velocities are depth- invariant. This is commonly known as the ‘ice-shelf approximation’. Treating the ice front as a perfectly heat-insulating interface, the model equations are solved for ice-flow velocity by considering the balance of forces at an ice edge of idealized rectangular shape (Weertman, 1957). The computation of an approximate solution of the model equations for gridcells of 1.25 km x 1.25 km is performed using the numerical procedures described by Grosfeld and Sandhäger (2004). The ice-shelf model supports the implementation of mechanical decoupling along zones subject to enhanced shear stresses, which was necessary to simulate the flow of Larsen B ice shelf prior to its collapse (Sandhäger and others, 2005). However, decoupling is not observed at Larsen C, so this capability is not required in our model runs.

The model’s primary function is to generate physically appropriate velocity fields under present and future ice-shelf geometries, to allow us to investigate fracture initiation from the stress distribution. We take this approach, rather than using direct observations of surface velocity, for two reasons. Firstly, this allows consistency between analyzing present (for which we have observations) and future (for which we need to predict flow rates) ice-shelf scenarios. Secondly, the observations of surface velocity available to us are suboptimal in terms of continuity and high spatial frequency noise and would be inappropriate for stress calculations without unreasonable spatial filtering. At present, it is not possible to achieve this from velocity observations as these are too noisy (high-frequency spatial variability), but satellite-based observations do allow the model to be validated for current conditions.

Boundary conditions

Ice-flow velocities at the landward boundary of the model domain, including those of the feeding glaciers, were obtained from satellite data. More specifically, the gaps in the RAMP (RADARSAT-1 Antarctic Mapping Project; Jezek, 2002) dataset, particularly near the ice shelf’s grounding line, were filled with new velocities derived from feature tracking between Landsat image pairs (Fig. 2; e.g. Strozzi and others, 2002). Encouragingly, our new velocities agreed well with those from RAMP where there was overlap between the two datasets. Since Landsat image pairs from different years (2000/01, 2006/07) show no significant differences, ice-flow velocities at the grounding line are assumed to have been invariant during this period.

Fig. 2. (a) Flow velocities as published in RAMP (Jezek and others, 2002). (b) Flow velocities derived from Landsat feature tracking. (c) Gridded velocity mosaic (moving average interpolation applied) from the combination of the two datasets.

Ice-shelf geometry, density and flow-rate factor

The grounding zone of the Larsen C ice shelf was determined from a European Remote-sensing Satellite (ERS) tandem mission interferogram; narrow ‘fringes’ readily indicated the position of the zone of tidal flexure (Fig. 3). We defined the seaward boundary of these fringes in our interferometry data as the boundary of the model domain considering only floating parts of the ice shelf unaffected by tidal flexure (Fig. 3).

Fig. 3. Example section of Larsen C showing a comparison of MOA grounding line (http://nsidc.org/data/nsidc-0280.html) and boundary of model domain. Interferogram: ERS-1/2 tandem mission (16/17 November 1995).

The mean temperature of the ice shelf is an essential forcing parameter and its impact on ice rheology is parameterized by the temperature-dependent flow-rate factor B. In previous model studies focusing on the Larsen B ice shelf, located further north, values of B = 420 kPa a^1/3 (Scambos and others, 2000) and values in a wide range from 222 to 601 kPa a^1/3 (Khazendar and others, 2007) have been used. We chose a constant flow-rate factor B of 435 kPa a^1/3 within this range, corresponding to an ice temperature of −12°C (Paterson, 1994), to reflect the more southerly location of the Larsen C ice shelf. Controlled-source seismic measurements (Fig. 1) in the 2008/09 austral summer were inverted for density–depth profiles using the methodology described by King and Jarvis (2007) (Fig. 4). Knowledge of these profiles is essential not only for solving the model equations, but also for improving ice-thickness estimates. We calculated ice thicknesses from Ice, Cloud and land Elevation Satellite (ICESat) elevation data (http://nsidc.org/data/nsidc-0304.html) using the newly derived mean density distribution in the ice column, and were thus able to update the ice-thickness distribution reported by BEDMAP (Lythe and others, 2001). Discrepancies between the two datasets were significant only close to the calving front (Fig. 5). This is unsurprising since ice-shelf basal melting is expected to be enhanced in this region (e.g. Jenkins, 1991).

Fig. 4. Ice density profile derived from inversion of controlled-source seismic data.

Fig. 5. (a) Ice-thickness distribution inferred from integrated BEDMAP (Lythe and others, 2001) radar data and ICESat elevation data (http://nsidc.org/data/nsidc-0304.html). (b) Difference from the thickness distribution published in BEDMAP.

Modelled Ice-Flow Velocities

The modelled flow-velocity field of the Larsen C ice shelf, simulated by adopting the specifications and input parameters listed above, reflects low velocities in the constrained inlets and accelerating flow towards the calving front (Fig. 6a). Maximum ice-flow velocities (∼750 m a⁻¹) dominate near the centre of this front, and velocity gradients are highest close to the grounding zone, particularly at the tips of the Kenyon and Churchill Peninsulas (cf. Figs 1 and 6). Comparison with satellite-derived surface velocities confirms that the velocity pattern of the ice shelf is well reproduced by the model results (Fig. 6). Encouragingly, it is not only the modelled velocity magnitudes that match the observations well, but also the modelled flowline trajectories (Fig. 7). Discrepancies between modelled and observed flowlines are generally minor and are focused on heavily rifted areas (Fig. 7). Further validation of the modelled flow-velocity field comes from two in situ Leica 1200 GPS sensors. These were deployed in the southeastern sector of the Larsen C ice shelf (Fig. 1) for a period of ∼7 weeks (December 2008–February 2009). The inferred mean magnitude of flow velocities during this period was ∼505 m a⁻¹ and thus lies within 5% of the modelled velocities (cf. Figs 1 and 6). The measured velocities are slightly higher than the modelled ones.

Fig. 6. (a) Modelled flow-velocity distribution of the Larsen C ice shelf. (b) Combined velocity from RAMP/feature-tracking grid. (c) Residual between the two datasets (modelled results – observed results). CP: Churchill Peninsula; KP: Kenyon Peninsula.

Fig. 7. Comparison of modelled (yellow) and observed (red) flow trajectories.

Nevertheless, there are some confined regions of the ice shelf where the difference between model and observations is significant (>150 m a⁻¹). The model generally underestimates flow velocities, except for Revelle Inlet south of the Kenyon Peninsula where modelled velocities exceed the observations. Such discrepancies are concentrated at the ice-shelf front and in its southern-sector areas such as down-flow of the eastern tip of the Kenyon Peninsula where major rift systems dominate (Fig. 6c). The inferred discrepancies arise because the model is not currently parameterized for ice-shelf mechanical heterogeneities such as those generated by large rifts and their infill. Rift widening may generate an additional component of ice flow that is not captured by the model. This is probably also the case for the large rifts (50 km) parallel to the central calving front (directly south of Bawden Ice Rise; Fig. 1). These rifts delineate a nascent tabular iceberg that has not yet detached due probably to the presence of softer marine ice (Holland and others, 2009).

In summary, good agreement between predicted and observed flow velocities confirms that our flow model is performing well. The modelled velocity distribution can therefore be considered a reliable basis for calculating stresses and analysing possible crevasse opening within the ice shelf.

Stress Intensity Factor and Crevasse Opening

To identify regions of potential crevasse opening, we applied the two-dimensional fracture criterion of Erdogan and Sih (1963) for the initiation of sharp cracks, first applied to ice by Shyam Sunder and Wu (1990). It was applied successfully to the Filchner–Ronne (Rist and others, 1999) and Larsen B (Bailey and Sammonds, 2007) ice shelves and focuses on the stress intensity factor, K, describing stress intensities in the vicinity of crack tips (Irwin, 1958). K is proportional to the applied stress and the initial half-length, c, of an already existing nuclear flaw within the ice (Rist and others, 1999). As the fracture criterion takes into account pure tensile (mode I) stresses as well as shear stresses in the plane of the crack (mode II), the corresponding stress intensity factors can be written as and , respectively. σ _n represents the stress normal to the crack plane and τ the shear stress, which is again dependent on the friction coefficient, μ.

The fracture criterion reads (for a detailed deduction see Rist and others, 1999):

(1)

with a crack growth initiation angle ϑ ₀, for which the left-hand side of Equation (1) is maximized, the mode I and II stress intensity factors K _I and K _II, and the critical stress intensity K _IC. A crack propagates if the critical threshold of K, the fracture toughness or critical stress intensity, K _IC, is exceeded. K _IC can be approximated as being linearly dependent on ice density (K _IC = 0.257 ρ – 80.7) and typically ranges from 50 kPa m^−1/2 for firn to 150 kPa m^−1/2 for consolidated ice (Rist and others, 1999). Here near-surface density was inferred as ∼512 kg m⁻³ (Fig. 4), yielding K _IC ≈ 51 kPa m^−1/2, which is entirely consistent with the values reported by Rist and others (1999). We therefore adopt a fracture toughness, K _IC,of 50 kPa m^−1/2 in our investigation of the fracture mechanics of the Larsen C ice shelf. For each gridcell, using modelled ice-flow velocities, we derived the principal stresses at the ice surface. Stress intensities were then calculated, following Rist and others (1999) and adopting a friction coefficient, μ, of 0.1 and several initial flaw sizes, c. Stress intensities generally exceed the fracture toughness, K _IC, of ice close to the landward boundary of the model domain (areas shaded in red in Fig. 8, showing the results for c = 0.08 m). Owing to the higher fracture toughness of 150 kPa m^−1/2 for consolidated meteoric ice, the initiation of basal crevasses would need a significantly larger initial flaw size (c = 0.7 m) for the same distribution of areas with crevasse propagation shown in Figure 8.

Fig. 8. Application of the fracture criterion. Shaded areas indicate regions where modelled stress intensities exceed the fracture toughness of the ice and thus support crevasse opening (image: Moderate Resolution Imaging Spectroradiometer (MODIS) 2002). White boxes outline the images shown in Figure 9. CP: Churchill Peninsula; KP: Kenyon Peninsula.

Additional regions for possible crevasse opening (cf. Figs 1 and 8) include: (1) those downstream of Francis Island and the Joerg and Kenyon Peninsulas; and (2) the transition zones between confined inlets and comparatively unconfined shelf ice in the northern sector of Larsen C (labelled northern transition zone in Fig. 8). In all of these regions the ice is subject to enhanced tensile horizontal stresses generated as ice leaves confined inlets and enters the main body of the ice shelf where flow accelerates considerably (Figs 1 and 8). This situation is compounded in the north of Larsen C, where both velocity magnitudes and flow directions undergo marked changes, of which the latter causes a bending moment which may generate additional stresses that promote further crevasse opening in the inlet–shelf transition zones (Fig. 8). The modelled orientation of the maximum tensile stress is perpendicular to the observed crevasses in this ice-shelf region (Fig. 9c). This consistency confirms the quality of the model results as cracks may initially grow in directions that are oblique to the highest tensile stress due to additional shear stress components, but will eventually tend to orientate orthogonally to it (Van der Veen, 1998; Rist and others, 1999).

Fig. 9. Detail of crevasse orientation and critical stress intensity for selected areas. (a) Crevasse orientation (white) in comparison with observed crevasses (black) downstream of Francis Island. (b) Regions where the critical stress intensity factor of 50 kPa m^−1/2 is exceeded. (c, d) Same parameters for the northern transition zone.

In the south of the Larsen C ice shelf, multiple crevasse plumes are advected into the ice shelf on either side of Francis Island and of the Joerg Peninsula (Figs 1 and 9b). Our model results are consistent with further opening of such crevasses downstream of these two promontories, or possibly the growth of other nucleated faults in the advected ice. Indeed, prominent new crevasses open in a flowband in between and down-flow of these promontories (labelled new crevasses in Fig. 8; cf. Fig. 9b). Crevasse orientation is again in agreement with the modelled direction of the maximum tensile stresses (Fig. 9a and b). The Kenyon Peninsula and Table Nunatak, located just off its eastern tip, are causing extensive rifting in the southern sector of Larsen C (Figs 1 and 8). The actual initiation point of the series of rifts, prominent between this peninsula and Gipps Ice Rise (Fig. 1), is located within the grounding zone immediately adjacent to the model domain. However, further crevasse opening is modelled downstream of the Kenyon Peninsula and Table Nunatak (Fig. 8), which is consistent with readily observed rift widening (Fig. 1).

Discussion

Predicted crevasse opening and observed presence of large crevasses are generally in good agreement for the regions downstream of Francis Island and the northern inlets, when an initial flaw size, c, of 0.08 m and a friction coefficient, μ, of 0.1 are assumed. As calculated stress intensities and directions are only weakly dependent on μ (Rist and others, 1999), we are principally concerned with potential variability in c. We found that stress intensity magnitudes scaled with c, while the spatial pattern of modelled stress intensities remained unaffected when c was varied in our model runs. The predicted pattern of crevasse opening (Fig. 8) is therefore unique, and thus reliable as a diagnostic indicator of ice-shelf fracture-mechanical processes. The physical origin, however, remains ambiguous, since the same spatial pattern of crevasse opening as that in Figure 8 could be obtained by systematic covariation of fracture toughness, K _IC, and c (Rist and others, 1999).

There are two prominent ice-shelf regions where significant crevasse opening is predicted but not observed (cf. Figs 1, 8 and 9): (1) downstream (where predicted stress intensity has an ice-shelf wide maximum) and east of the Cole Peninsula in the north; and (2) downstream of the Joerg Peninsula in the south. All observed fractures in these regions are advected into the model domain from the interior of the ice sheet and, as such, are not compliant with our fracture criterion. We believe that these discrepancies between model and observation are evidence of the accretion of marine ice in regions (1) and (2), sustaining softer ice downstream and thus inhibiting crevasse opening, as suggested by Glasser and others (2009) and Holland and others (2009). Since our model is not currently parameterized to allow for such ice-mechanical heterogeneities, model outputs would not be expected to be supported by observations.

Flow modelling and the fracture criterion thus place quantitative constraints on the stabilizing effects of marine ice on the Larsen C ice shelf, which diminishes rates of rift or crevasse propagation. A reduction in marine ice production would therefore result in a weaker ice shelf (Glasser and Scambos, 2008; Holland and others, 2009; Khazendar and others, 2009 (Brunt Ice Shelf)). In the case of Larsen B, regions down-flow of notable promontories in the north and south were inferred to be anomalously mechanically weak, thus allowing unusually high flow-velocity gradients (Sandhäger and others, 2005; Vieli and others, 2006, 2007; Khazendar and others, 2007). Ice-mechanical heterogeneities would as such have supported mechanical decoupling near the northern and southern margins of Larsen B, as well as acceleration of its main body and ultimately, therefore, ice-shelf disintegration. These heterogeneities are reflected in the high spatial variability of the flow parameter, B, published by Khazendar and others (2007). As Larsen C is not currently characterized by Larsen B-style marginal velocity gradients that are anomalously high, it is possible to achieve very good agreement between observed and modelled ice-flow velocity with a constant flow parameter for the entire ice shelf. We may speculate, however, that Larsen C’s velocity distribution may become subject to decoupling of the middle part of the ice shelf, similar to a Larsen B-style dynamic regime, if rates of marine ice production decrease in the lee of the Cole or Joerg Peninsulas. Decoupling of Larsen C’s main body from the northern inlets, as well as from the heavily crevassed region between the Kenyon Peninsula and Gipps Ice Rise, could promote ice acceleration and modification of the ice shelf’s stress regime and, as such, lead to dynamic conditions akin to those that led to disintegration of the Larsen B ice shelf.

Conclusion

The Larsen C ice shelf is inferred to be stable in its current dynamic regime. Ice-mechanical heterogeneities in ice-stream suture zones, sustained by marine ice production down-flow of promontories, have significant stabilizing effects on the ice shelf. Reduction in rates of marine ice production could therefore lead to weakening of suture zones and possibly development of Larsen B-style dynamic conditions prior to its disintegration. This emphasizes the importance of further research into the mechanics of suture zones and their dependence on marine ice provenance, together with thorough quantification of their modification of the ice-shelf stress regime and thus its stability. We plan to extend our continuum-mechanical flow model and fracture criterion to allow for ice-shelf mechanical heterogeneities, supporting predictions of the future evolution of the Larsen C ice shelf with increased confidence.