Introduction

The recent disintegration of ice shelves on the Antarctic Peninsula is closely linked to climatically controlled southward migration of the −9°C isotherm (Reference MercerMercer, 1978; Reference VaughanVaughan and others, 2003). The retreat pattern of ice shelves is governed by the overall ice-shelf mass balance including processes such as calving (Reference SkvarcaSkvarca, 1994), disintegration (e.g. Reference Rott, Skvarca and NaglerRott and others, 1996), surface and basal melting (Reference Scambos, Hulbe, Fahnestock and BohlanderScambos and others, 2000) and also the structure (Reference Glasser and ScambosGlasser and Scambos, 2008; Reference GlasserGlasser and others, 2009) and composition of the ice (see Reference Cook and VaughanCook and Vaughan, 2009). The Larsen C ice shelf has been thinning over the last two decades (Reference Shepherd, Wingham, Payne and SkvarcaShepherd and others, 2003), but otherwise does not currently exhibit obvious climatically related signs of retreat (Reference GlasserGlasser 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 (Reference Sandhäger, Rack and JansenSandhäger and others, 2005) and linear elastic fracture mechanics (Reference Rist, Sammonds, Murrell, Meredith, Oerter and DoakeRist and others, 1996, Reference Rist1999, Reference Rist, Sammonds, Oerter and Doake2002), 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 Reference Vieli, Payne, Du and ShepherdVieli and others (2006, Reference Vieli, Payne, Shepherd and Du2007) and Reference Khazendar, Rignot and LarourKhazendar and others (2007, Reference Khazendar, Rignot and Larour2009) 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. Reference GlasserGlasser and others, 2009; Reference Holland, Corr, Vaughan, Jenkins and SkvarcaHolland 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² (Reference Cook and VaughanCook 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; Reference Cook and VaughanCook and Vaughan, 2009; Reference GlasserGlasser 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; Reference SkvarcaSkvarca, 1994). The ice tends to be heavily crevassed in the transition zones between such feeders and the floating ice shelf (Reference GlasserGlasser 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 (Reference Cook and VaughanCook 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 (Reference Humbert, Kleiner, Mohrholz, Oelke, Greve and LangeHumbert 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

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; Reference JezekJezek, 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. Reference Strozzi, Luckman, Murray, Wegmuller and WernerStrozzi 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 (Reference JezekJezek 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 (Reference Scambos, Hulbe, Fahnestock and BohlanderScambos and others, 2000) and values in a wide range from 222 to 601 kPa a^1/3 (Reference Khazendar, Rignot and LarourKhazendar 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 (Reference PatersonPaterson, 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 Reference King and JarvisKing 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 (Reference Lythe and VaughanLythe 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. Reference JenkinsJenkins, 1991).

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

Fig. 5. (a) Ice-thickness distribution inferred from integrated BEDMAP (Reference Lythe and VaughanLythe 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 (Reference Holland, Corr, Vaughan, Jenkins and SkvarcaHolland 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 Reference Erdogan and SihErdogan and Sih (1963) for the initiation of sharp cracks, first applied to ice by Reference Shyam Sunder and WuShyam Sunder and Wu (1990). It was applied successfully to the Filchner–Ronne (Reference RistRist and others, 1999) and Larsen B (Reference Bailey and SammondsBailey and Sammonds, 2007) ice shelves and focuses on the stress intensity factor, K, describing stress intensities in the vicinity of crack tips (Reference Irwin and FlüggeIrwin, 1958). K is proportional to the applied stress and the initial half-length, c, of an already existing nuclear flaw within the ice (Reference RistRist 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 Reference RistRist 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 (Reference RistRist 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 Reference RistRist 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 Reference RistRist 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 (Reference Van der VeenVan der Veen, 1998; Reference RistRist 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 μ (Reference RistRist 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 (Reference RistRist 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 Reference GlasserGlasser and others (2009) and Reference Holland, Corr, Vaughan, Jenkins and SkvarcaHolland 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 (Reference Glasser and ScambosGlasser and Scambos, 2008; Reference Holland, Corr, Vaughan, Jenkins and SkvarcaHolland and others, 2009; Reference Khazendar, Rignot and LarourKhazendar 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 (Reference Sandhäger, Rack and JansenSandhäger and others, 2005; Reference Vieli, Payne, Du and ShepherdVieli and others, 2006, Reference Vieli, Payne, Shepherd and Du2007; Reference Khazendar, Rignot and LarourKhazendar 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 Reference Khazendar, Rignot and LarourKhazendar 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.