Creep deformation and buttressing capacity of damaged ice shelves: theory and application to Larsen C ice shelf

. Around the perimeter of Antarctica, much of the ice sheet discharges to the ocean through ﬂoating ice shelves. The buttressing provided by ice shelves is critical for modulating the ﬂux of ice into the ocean, and the presently observed thinning of ice shelves is believed to be reducing their buttressing capacity and contributing to the acceleration and thinning of the grounded ice sheet. However, relatively little attention has been paid to the role that fractures play in the ability of ice shelves to sustain and transmit buttressing stresses. Here, we present a new framework for quantifying the role that fractures play in the creep deformation and buttressing capacity of ice shelves. We apply principles of continuum damage mechanics to derive a new analytical relation for the creep of an ice shelf that accounts for the softening inﬂuence of fractures on longitudinal deformation using a state damage variable. We use this new analytical relation, combined with a temperature calculation for the ice, to partition an inverse method solution for ice shelf rigidity into independent solutions for softening damage and stabilizing backstress. Using this new approach, ﬁeld and remote sensing data can be utilized to monitor the structural integrity of ice shelves, their ability to buttress the ﬂow of ice at the grounding line, and thus their indirect contribution to ice sheet mass balance and global sea level. We apply this technique to the Larsen C ice shelf using remote sensing and Operation IceBridge data, ﬁnding damage in areas with known crevasses and rifts. Backstress is highest near the grounding line and upstream of ice rises, in agreement with patterns observed on other ice shelves. The ice in contact with the Bawden ice rise is weakened by fractures, and additional damage or thinning in this area could diminish the backstress transmitted upstream. We model the consequences for the ice shelf if it loses contact with this small ice rise, ﬁnding that ﬂow speeds would increase by 25 % or more over an area the size of the former Larsen B ice shelf. Such a perturbation could potentially destabilize the northern part of Larsen C along pre-existing lines of weakness, highlighting the importance of the feedback between buttressing and fracturing in an ice shelf.


Introduction
The majority of the Antarctic ice sheet drains to the ocean through floating ice shelves (Barkov, 1985), most of which are contained in embayments or run aground against ice rises, ice rumples or islands. These pinning points buttress the flow of neighboring grounded ice (Thomas, 1979) and influence the position of the grounding line, where the ice detaches from the bed and becomes afloat in the ocean (Thomas, 1979;Schoof, 2007;Rignot et al., 2008;Gagliardini et al., 2010;Pritchard et al., 2012;Gudmundsson, 2013). Thus ice shelves play a major role in modulating the mass balance and contribution to sea level rise of the Antarctic ice sheet. This influence was brought into sharp focus following the collapse of the Larsen B ice shelf in 2002, after which the tributary glaciers that fed the shelf accelerated three-to eightfold Scambos et al., 2004) with sustained dynamic thinning and retreat ongoing (Rott et al., 2011).
Ice shelves are thinning in many sectors of Antarctica, primarily a result of ocean-driven basal melting (Rignot and Jacobs, 2002;Pritchard et al., 2012;Rignot et al., 2013). At the same time, the outlet glaciers and ice streams that flow into these ice shelves are thinning (Rignot et al., 2008;Pritchard et al., 2009), a correlation attributed to the diminished buttressing provided by the ice shelves (e.g. Rignot et al., 2008;Pritchard et al., 2009;Gagliardini et al., 2010). The implied reasoning is that as ice shelves thin, the lateral surface area over which shear stress can be transmitted is reduced, which decreases the total buttressing force provided by the shelf (e.g. . Reduced ice shelf buttressing can cause retreat of the grounding line and increased ice sheet mass loss Schoof, 2007).
However, the thinning of an ice shelf also makes it more susceptible to fracture (Shepherd et al., 2003). Any reduction in resistive backstress, from thinning or weakening of the ice, increases the net longitudinal tensile stress in the ice (Jezek, 1984;Rist et al., 2002), which may allow new fractures to form and existing fractures to penetrate deeper into or across the ice. During ice shelf retreat, stabilizing features such as ice rises and ice rumples can transition into being nucleation points for fractures and actually contributing to the destabilization of the shelf (Doake and Vaughan, 1991). Thus the decreased buttressing provided by thinning ice shelves could be due to -or compounded by -mechanical weakening associated with fractures.
Given their sensitivity to atmospheric and oceanic changes and their role in modulating the discharge of grounded ice into the ocean, and thus to global sea level, there is a need for improved models and observations that link ice shelf buttressing to mechanical weakening. Buttressing has been identified as one of the least understood processes in marine ice sheet dynamics (e.g. Schoof, 2007), which has motivated a number of recent numerical studies of the relationship between ice shelf buttressing and grounding line stability for a marine ice sheet (Goldberg et al., 2009;Gagliardini et al., 2010;Favier et al., 2012;Gudmundsson, 2013). However, in spite of the widespread recognition of the importance of fractures for the flow and stability of ice shelves (e.g Jezek et al., 1985;Doake and Vaughan, 1991;Vaughan, 1993;van Der Veen, 1998a;Scambos et al., 2000;Kenneally and Hughes, 2004;Larour et al., 2004a;Glasser et al., 2009;Khazendar et al., 2009;Jansen et al., 2010;Albrecht and Levermann, 2012;Luckman et al., 2012;McGrath et al., 2012), little attention has been given to the impact of fractureinduced weakening on ice shelf buttressing.
Recently, first steps have been taken to incorporate the physical effects of fractures into large-scale ice sheet models (Albrecht and Levermann, 2012;Borstad et al., 2012). Here, we build on these recent efforts by introducing a theoretical framework for quantifying the buttressing capacity and mechanical integrity of an ice shelf that may be weakened by fractures. We apply the principles of continuum damage mechanics to the standard momentum balance equations for the creep deformation of an ice shelf to account for the rheological weakening caused by fractures using a scalar damage variable. In this context, damage is a state variable that accounts for the influence of fractures on longitudinal deformation visible at the surface of the ice shelf. In addition to characterizing the likelihood of calving and the overall stability of an ice shelf (e.g. Borstad et al., 2012), damage feeds back with the ability of the ice to support and transmit backstress. We calculate independent solutions for damage and backstress from the results of an inverse method solution for the spatial rheology of an ice shelf combined with an estimate of the bulk temperature of the ice. Using this framework, the magnitude of damage and backstress for an ice shelf can be quantified and monitored through time using remote sensing data.
We begin with a review of the analytical theory for ice shelf creep developed by Weertman (1957) and generalized by Thomas (1973a). A linear mapping from continuum damage mechanics theory is then applied to transform this depthaveraged creep equation into one that depends on damage. We apply this new framework using data for the Larsen C ice shelf covering the period 2006-2009, concluding with a discussion on possible extensions and additional applications of this new approach.
2 Background Weertman (1957) derived an analytical expression for the creep of an idealized ice shelf from first principles combined with the empirical relation governing the power law creep of glacier ice. For an ice shelf free to deform in only one direction, taken here as the x direction, with y and z the lateral and vertical coordinates, respectively, and with the origin at sea level, the analytical solution for ice shelf creep iṡ whereε xx is the longitudinal strain rate (generally assumed to be independent of depth), ρ i is the density of ice, g is gravitational acceleration, H is the ice thickness, B is the ice rigidity, ρ w is the density of seawater, n is the flow law exponent, and overbars indicate depth-integrated quantities (Weertman, 1957). For convenience, we write this aṡ where = ρ i (1 − ρ i /ρ w ). The buoyancy-driven longitudinal stress in an ice shelf is 1/2 gH 2 and arises due to the density difference between ice and seawater. Vertical variations in ice density and ice rigidity can be accounted for by writing the rigidity term as where b and s are the vertical coordinates representing the base and surface of the ice shelf, respectively (Thomas, 1973a). Hereafter we omit the overbars and assume that all relevant quantities are depth-integrated unless otherwise specified. Thomas (1973a) generalized the theory of Weertman (1957) to account for all horizontal components of the strain rate tensor as well as the presence of backstress due to the ice shelf running aground against ice rises, ice rumples or lateral margins. This generalized relation for the creep of an ice shelf can be written aṡ where σ b > 0 is the scalar backstress and α =ε yy /ε xx and β =ε xy /ε xx represent the contributions of lateral and shear strain rate, respectively (Thomas, 1973a). The flow is considered converging for α < 0 and diverging for α > 0, with parallel flow for α = 0 . For α = −2 the definition of the effective stress gives 1/2 gH = σ b , therefore the equation is insoluble (the right-hand side is 0/0). Note thatε xx is sensitive to errors in α for strongly converging flow where α < −1 . In Eq. (5) the strain rate is negative in areas where σ b > 1/2 gH , which occurs just upstream of ice rises (Thomas, 1973b), or where the flow is strongly compressive longitudinally or laterally and thus α < −2. However, where the longitudinal strain rate is tensile, which covers the majority of a typical ice shelf, the following simplified form of Eq. (5) can be used, where θ represents the contribution of lateral and shear strain rate terms (Thomas, 1973a) Equation (6) is the form most commonly cited in the literature (e.g. Thomas and MacAyeal, 1982;Jezek, 1984;Jezek et al., 1985;Rist et al., 2002;, though it must be kept in mind that it is only a specialized case of Eq. (5) and care must be taken to ensure that the appropriate relation (Eq. 5) is used in regions where the flow has a compressive component, especially if Eq. (5) is solved for a term other than the strain rate.
Equations (5) and (6) are derived from the relationship for the effective stress and are therefore invariant under coordinate transformation (Jezek et al., 1985). The x direction is commonly taken as the flow direction (e.g. Thomas, 1973a, b;Thomas and MacAyeal, 1982), though it is often convenient to choose the x direction coincident with the largest positive (tensile) principal strain rate, for which β = 0 (e.g. Thomas, 1973b;MacAyeal and Holdsworth, 1986). For an ice shelf free to creep in only one direction, such as a shelf embayed between parallel frictionless walls,ε yy = 0 anḋ ε xy = 0 and therefore α = β = 0, θ = 1/8, and the unidirectional solution of Weertman (1957) is obtained (Eq. 2). For an unbounded ice shelf free to creep in any horizontal direction, α = 1, β = 0 and the unconfined solution of Weertman (1957) is recovered.
The backstress term σ b can be partitioned into individual terms representing the influence of lateral or marginal shear stress between an ice shelf and its embaying walls, backstress caused by locally grounded areas of ice, or resistance caused by the compressive confluence of neighboring tributary glaciers in an ice shelf (Thomas, 1973b. Under idealized scenarios it is possible to formulate explicit expressions for these individual resistance terms, but in general the backstress needs to be solved for in terms of the other known variables in Eqs. (5) or (6) (e.g. Thomas and MacAyeal, 1982).
Several studies have incorporated principles of fracture mechanics into the analytical framework of Thomas (1973a). Jezek (1984) independently calculated backstress on the Ross Ice Shelf in areas where the penetration height of bottom crevasses was known from radar sounding experiments, assuming that the difference between measured crevasse height and that predicted by fracture mechanics was due to the stabilizing presence of backstress. This development introduced a crack length explicitly into the formulation of Thomas (1973a), though it only applies in a field of closely spaced crevasses. Rist et al. (2002) and Kenneally and Hughes (2004) similarly invoked backstress to rectify the discrepancy between measured and predicted basal crevasse penetration heights, this time using Linear Elastic Fracture Mechanics (LEFM). In this approach, the crack length enters the framework implicitly through the implied backstress. These studies formulate the full stress in an ice shelf to calculate crack initiation and arrest, whereas only the deviatoric stress is represented in the creep equations, a distinction that could be important in interpreting patterns of fractures or damage in an ice shelf and in distinguishing between the presence of fractures and the influence of fractures on buoyancy-driven creep.
These techniques for using basal crevasse penetration heights to constrain the backstress only apply to areas where basal crevasses are present. Furthermore, they rely on fracture mechanics assumptions such as linear elasticity, homogeneity, and fully brittle failure that may not hold, or may only hold over a short time period when the crack first forms www.the-cryosphere.net/7/1931/2013/ The Cryosphere, 7, 1931-1947, 2013 and the deformation is predominantly elastic (e.g. MacAyeal and . Nevertheless, the motivation for accounting for the influence of fractures in the theoretical framework of Thomas (1973a) is sound, given our current understanding of the importance of fractures on the evolution, stability, and collapse of ice shelves. Within this context, a theory such as continuum damage mechanics can fill a longstanding need by providing a framework that accounts for the influence of fractures of any type, anywhere in an ice shelf.

Damage theory
Continuum damage mechanics is a theory that accounts for the effects of fractures on material behavior while maintaining a continuum representation of the material (e.g. Lemaitre, 1996). The spatial scale of interest in mechanical modeling is typically the macroscale response of a material or structure, which may be orders of magnitude larger than the fractures that exist at the microscale. Even though fractures are by nature local phenomena, and their spatial extent may be small relative to the macroscale response of interest, their effect on deformation or strain can be measured at the macroscale (Murakami and Ohno, 1981). Ignoring these features simply because they cannot be resolved by a model is not prudent. In damage mechanics, the effects of fractures at the microscale are accounted for through the definition of effective material properties or state variables without the need to resolve these features individually. Ice shelves typically have spatial extents on the order of 10-100 km, and most numerical ice shelf models discretize the governing equations over a spatial scale on the order of 0.1-1 km. Surface crevasses can penetrate up to 40-50 m deep in cold ice, less in temperate ice, and more in the presence of surface meltwater (van der Veen, 1998b), and basal crevasses can propagate much further, even as high as sea level (e.g. Jezek, 1984). These scales are smaller than can be resolved with a typical ice shelf model, moreover it would not be feasible to represent explicitly in three dimensions the mechanics of the hundreds to thousands of crevasses that may exist over an entire ice shelf. Most ice shelf models use a two-dimensional depth-integrated formulation of the governing equations (Shallow Shelf Approximation or SSA, MacAyeal, 1989) which precludes explicit modeling of features that occupy only a fraction of the ice thickness.
Damage mechanics is not limited to representing microscale fractures exclusively, however. Damage models are also appropriate for representing the propagation of macroscopic fractures. This is especially the case for heterogeneous materials, where the macroscopic crack tip is inherently surrounded by a zone of distributed microcracking as a result of material or structural heterogeneity. Modeling such fractures by the propagation of a damage field that smears the heterogeneity over an appropriate length scale can ensure the cor-rect dissipation of mechanical energy where a singular crack model may not (Bažant and Jirásek, 2002).
For an ice shelf, therefore, the propagation of throughthickness rifts can also be represented using damage mechanics, though it is also possible to explicitly represent such macroscopic fractures as discontinuous features in a numerical model (Larour et al., 2004a, b). However, rifts may not behave as classical brittle cracks with well-defined stress singularities near their tips, as must be assumed in order to apply a linear fracture theory such as Linear Elastic Fracture Mechanics (LEFM). Instead, there appears to be a length scale of about 1 km that defines the distributed zone of microcracking and strain localization surrounding the crack tip (Bassis et al., 2007). This large distributed zone of cracking, known as the "fracture process zone" in the parlance of quasi-brittle fracture mechanics, engenders nonlinearity into the fracture scaling, necessitating either a quasi-brittle (nonlinear) representation of the fracture (e.g., for snow, Borstad and Mc-Clung, 2013) or the use of a nonlocal damage model (e.g. Borstad and McClung, 2011).
Modeling fractures using either a nonlinear fracture theory or a nonlocal damage model requires certain conditions be met to ensure that the results are not sensitive to the chosen discretization scheme. This is typically represented as a restriction on the maximum element size. A common rule of thumb is that the element size should be no larger than about one third of the process zone size to satisfy mesh objectivity (Bažant, 2005). This ensures that damage does not localize down to the size of a single element and that the direction of propagation is not influenced by the orientation of the mesh or grid. A properly formulated nonlocal model can also represent the initiation of a fracture from a smooth boundary without the need for a prescribed flaw or notch (Borstad and McClung, 2011;Bažant, 2005).
A final advantage of damage mechanics over brittle fracture mechanics is that it works within the existing constitutive framework of the material or structure. A state damage variable is introduced via a linear mapping applied to the relevant governing equations, in this case the relations for the creep of a floating ice shelf. We begin by discussing the nature and consequences of the chosen mapping scheme and then derive a new relation for the creep of ice shelves that depends on damage.

Linear mapping between physical and effective spaces
Damage is a state variable introduced to achieve a desired linear mapping between the actual physical state of the material, which may be fractured, damaged, or otherwise heterogeneous, and an effective state that is compatible with a homogeneous, continuum representation of the applicable governing equations. The nature of this mapping, and therefore the physical definition of the damage variable, depends on the chosen equivalence scheme between these two reference states. A common choice of equivalence scheme is that of strain (or strain rate) equivalence, which states that the actual material state and the effective state are mapped such that they achieve the same strain (or strain rate) under the same stress or loading. This equivalence scheme is commonly used for damage models applied to polycrystalline ice (e.g. Pralong and Funk, 2005; Duddu and Waisman, 2012) as well as cohesive snow (Borstad and McClung, 2011). Pralong et al. (2006) demonstrated that this equivalence scheme is thermodynamically consistent for damage as well as healing (damage reversal) of ice. We adopt this equivalence scheme here, which leads to the following definition of an effective stress, whereσ is the effective (damage-dependent) stress, σ is the Cauchy stress, and D is the scalar (isotropic) damage variable (e.g. Borstad et al., 2012;Duddu and Waisman, 2012;Pralong and Funk, 2005;Lemaitre, 1996). Damage ranges between zero, for fully intact ice, to one for ice which has fully failed or lost all load bearing capacity. When applying the strain equivalence principle to an arbitrary volume element, damage has a physical interpretation as the loss of load bearing surface area within the element. For scalar (isotropic) damage, damage is interpreted as representing the weakest cross section of the element since this section governs the ultimate load bearing capacity of the element as a whole. For fully viscous (long timescale) deformation, this type of simple area reduction due to fractures likely dominates the material response over the shape or orientation of the cracks (Murakami and Ohno, 1981), which favors a scalar representation of damage. However, damage can also be represented as a tensor to account for varying damage on different orthogonal planes, which can account for any anisotropy induced by the orientation of fractures in a material. Wu and Mahrenholtz (1993) applied such an anisotropic damage model to creep rupture data for polycrystalline ice and found that the model did not fit the experimental data any better than an isotropic model. However, insufficient experimental data exist for properly calibrating an anisotropic model and comparing it to a scalar model (Duddu and Waisman, 2012). Nonetheless, at the spatial scale of an ice shelf, and over long timescales associated with predominantly creep deformation, it remains to be demonstrated whether an anisotropic model is necessary to account for the orientation of crevasses and rifts. For simplicity we adopt the scalar isotropic damage mapping of Eq. (8), but note that a tensorial approach could be derived by replacing the scalar D with a tensor and carrying on with the proceeding derivation. Since we will be applying the linear mapping defined by Eq. (8) to a depth-integrated equation (Eq. 5), the simple geometric interpretation of damage described above is unfortunately lost. Instead, damage becomes a phenomenological parameter of the resulting equation.

Derivation of damage-dependent creep relation
Following a choice of equivalence principles and associated linear mappings between actual and effective spaces, the damage model is derived by substituting the effective stress σ in place of the Cauchy stress σ anywhere that the Cauchy stress arises in the applicable governing equations, in this case the creep relations for an ice shelf derived by Weertman (1957) and Thomas (1973a). Here we apply the linear mapping to the deviatoric stress, since this is the stress that governs the creep of an ice shelf and that appears in the solutions of Weertman (1957) and Thomas (1973a). However, formally the mapping applies to the full Cauchy stress tensor, which implies that damage will also map onto an effective pressure term, as demonstrated in a glaciological context by Pralong and Funk (2005). Since we do not calculate the pressure in the present analysis, and for brevity, we omit the full derivation here. Instead, we begin with an expression for the longitudinal strain rate of an ice shelf derived by Thomas (1973a), where σ xx is the (depth-integrated) longitudinal deviatoric stress. Substitution of the effective stress given by Eq. (8) in place of σ xx leads to Note thatε =ε by definition when applying the strain equivalence principle, therefore we omit the tilde here. The deviatoric longitudinal stress is given by (Thomas, 1973a) where F is the magnitude of backforce in opposition to the flow at a given point and the backstress is defined as σ b = F /H . Substituting Eq. (11) into Eq. (10) and simplifying leads tȯ This relation represents the creep of a damaged ice shelf, in a simplified form analogous to Eq. (6). The general form is analogous to Eq. (5), with B replaced by (1−D)B. This new relation accounts for areas which are fractured and which therefore serve to destabilize the shelf as well as for any buttressing stresses which act to stabilize the shelf. As with Eq. (5) caution near the grounding line where the assumption of hydrostatic equilibrium may not be satisfied and where vertical shearing may also be important. Furthermore, damage in this context does not take into account the possible influence of fractures on the response of an ice shelf to flexural stresses or modes of deformation other than longitudinal. This distinction may be important in the vicinity of the ice front where bending stresses are present or if widespread melt ponding induces flexural stresses in an ice shelf (e.g. MacAyeal and . According to Eq. (12), the strain rate of a damaged ice shelf is proportional to (1 − D) −n , indicating the sensitivity of the flow to the initiation or evolution of damage. For a moderate level of damage of D = 0.2, the strain rate approximately doubles compared to the solution for undamaged ice (assuming n = 3 and no change in backstress). For D = 0.5, a level approaching the threshold damage for calving identified by Borstad et al. (2012), the strain rate increases by a factor of 8. Thus damage can account for large increases in strain rate that would not be compatible with lower order influences such as warming of the ice or development of preferential crystal fabric, at least over short distances.
When solving Eq. (12) for D, care must be taken with the sign ofε xx since the strain rate can be negative (Thomas, 1973a). Eq. (13) solved for D is thus For negative (compressive) strain rates, the numerator of Eq. (13) is also negative since the condition σ b > 1/2 gH will hold, ensuring that D ≤ 1. For a consistent set of observations defining the right-hand side of Eq. (13), the second term will be less than one, ensuring D ≥ 0. Recall that when α = −2, the equality 1/2 gH = σ b holds by definition, making the equation insoluble.
Fractures are more likely to originate in areas of high tension, for which |α| < 2, or in areas of strong shear where lateral and longitudinal strain rates are of similar magnitude and thus α > −2. Providing that backstress is less than the driving stress, a condition which also holds over the majority of an ice shelf, we can write the expression for damage in the simplified form

Backstress and damage calculation
Given velocity and thickness data as well as a calculation of ice temperature or ice rigidity, there are still generally two unknowns in Eqs. (13) and (14): damage and backstress. As has been confirmed in previous studies that implicitly enforced D = 0 everywhere, backstress is present over much of an ice shelf (e.g. Thomas and MacAyeal, 1982). A number of analytical expressions for the backstress term have been proposed for idealized cases of lateral resistance or grounding at an ice rise or ice rumple (Thomas, 1973a, b), which may allow damage to be calculated independently from backstress for simple geometries. An alternative method of determining backstress, that we introduce here for the first time, is to use an inverse control method to provide an independent calculation of the ice shelf rheology as a starting point. The analytical model (Eq. 12) is then used, in combination with a calculation of the bulk temperature of the ice shelf, to partition the inverse method solution into separate and independent solutions for damage and backstress. Thus two constraints (the inverse solution for rigidity and the temperatureparameterized rigidity) are used to solve for the two unknowns, damage and backstress. We outline this approach by first noting that, since D is dimensionless, the term (1 − D)B in Eq. (12) can be aggregated into a bulk ice rigidity term B i (with the subscript i standing here for "inverse"). This aggregate term is obtained from an inverse method that seeks a spatially variable rigidity field (B i ) that minimizes the misfit between modeled and observed surface velocities (e.g. Larour et al., 2005). Anomalously low values of B i obtained from inverse methods, that is, values of B i that are lower than would be expected for pure ice at that appropriate bulk temperature of the shelf, are typically associated with fractured or damaged ice (Larour et al., 2005;Vieli et al., 2007;Khazendar et al., 2007Khazendar et al., , 2011. Conversely, anomalously high values of B i are typically found in areas where backstress is present. To partition the inverse method solution into damage and backstress, we first mask the inverse solution such that This masked inverse solution B i − contains information about both the temperature of the ice and any areas of softened or weakened ice, but not backstress, which comes from the analytical creep model using B i − as an input. The term B(T ) comes from a standard temperature parameterization for ice rigidity (e.g. Cuffey and Paterson, 2010) and determines the level of rigidity appropriate for an unconfined ice shelf of the given temperature. To calculate backstress, B i − is substituted into Eq. (12) in place of (1 − D)B and the equation is solved for σ b : To calculate damage, Eq. (16) is substituted into Eq. (13), which can then be simplified to temperature, but not on each other. The principal source of uncertainty in damage is the uncertainty in the appropriate value of ice rigidity B(T ) given the (typically unknown) temperature of the ice and possibly factors such as ice fabric, impurities, cavities, water content, etc. The uncertainty in backstress is further compounded by errors in observed or modeled strain rates.

Application to Larsen C ice shelf
We apply this new theory to assess the damage and buttressing capacity of the Larsen C ice shelf. Draining an area roughly five times larger than the Larsen B ice shelf before it disintegrated in 2002, Larsen C is the largest remaining ice shelf on the Antarctic Peninsula (Khazendar et al., 2011).
Though modeling studies suggest that the ice shelf appears to be stable at present (Jansen et al., 2010;Khazendar et al., 2011), remote sensing observations have indicated progressive thinning and surface elevation lowering of the ice shelf over the last three decades (Shepherd et al., 2003;Fricker and Padman, 2012). The northern half of Larsen C has accelerated since 2000 (Khazendar et al., 2011), leading to speculation about possible destabilization of the ice shelf in the coming years. The influence of fractures, crevasses and rifts on Larsen C has been inferred in a number of previous studies. Observations of fractures in the shelf have been linked to locally enhanced strain rates (Rack et al., 2000;Jansen et al., 2010). Khazendar et al. (2011) noted that rifts and other fractures were linked to large spatial variations in the inferred rheology of the shelf. Several recent studies of the ice shelf have observed basal crevasses along longitudinal transects of the shelf using ground-penetrating radar (McGrath et al., 2012;Luckman et al., 2012), indicating that basal crevasses are likely widespread on the ice shelf.

Surface velocity
We use horizontal surface velocity observations for the period 2006-2009 for the ice shelf (Rignot et al., 2011b), generated at a sample spacing of 150 m for Larsen C (Fig. 1) vs. 900 m in the remainder of Antarctica. For Larsen C, the data come primarily from the Japanese Aerospace Exploration Agency's Advanced Land Observing Satellite (ALOS) PALSAR, processed using speckle tracking (Mouginot et al., 2012). The median error in velocity magnitude is about 5 m yr −1 , which translates into a strain rate error on the order of 3 × 10 −4 yr −1 (Rignot et al., 2011b). The median error in velocity direction is 1.7 • , with higher error near the grounding line.
It is instructive to look at the ratios α =ε yy /ε xx and β =ε xy /ε xx to assess the flow of the ice shelf. Figure 2a shows areas of lateral convergence (α < 0) and divergence (α > 0) of the flow. Between the Bawden and Gipps ice rises downstream of Churchill and Hollick-Kenyon Peninsulas, the flow of the ice shelf is predominantly divergent. Other areas of divergent flow are evident closer to the grounding line, often after glaciers have passed through points of constriction and begin expanding laterally into the shelf. After a short distance, however, α becomes negative as neighboring tributary glaciers meet downstream of major promontories and provide lateral resistance to each other's viscous expansion into the shelf, thus causing compressive lateral strain rates (ε yy < 0). The high magnitude of α upstream of the Bawden and Gipps ice rises indicates the importance of these features to the force balance of the shelf. A sharp transition from positive to negative α occurs just upstream of the Bawden ice rise, a consequence of the longitudinal strain rate becoming negative near the ice rise due to the strong resistance encountered.
The magnitude of shear strain rates are substantially higher than longitudinal strain rates in the vicinity of the grounding line where glaciers flow past grounded lateral margins, especially along Churchill and Hollick-Kenyon Peninsulas (Fig. 2b). Gudmundsson (2013) noted similar results using a numerical model, namely that the largest (in magnitude) deviatoric stresses were along the shear margins of an ice sheet/ice shelf system. For Larsen C, the relative influence of shear diminishes with distance into the shelf and with distance from the ice rises.

Ice thickness
The thickness of the ice shelf is taken from the Bedmap2 data set (Fretwell et al., 2013), which for ice shelves is based on the data of Griggs and Bamber (2011), except for a zone of exclusion within 5 km of the grounding line that was filled, where available, with airborne thickness measurements. The seaward extent of the thickness data corresponds to the boundaries of the ice shelf in the 2003/2004 MODIS Mosaic of Antarctica (Scambos et al., 2007). Between 2003Between /2004Between and 2009, the ice front of Larsen C in some places extended beyond these earlier limits, but the ice front used in the following calculations was drawn to stay within the MOA-derived limits.

Temperature calculation
A steady-state temperature solution for the ice shelf is calculated as a function of surface and basal ice temperature and the mass flux at the base of the ice shelf. Basal melting rates (Fig. 3a) are computed using the Massachusetts Institute of Technology general circulation model (MITgcm) with a three equation thermodynamic representation of the freezing/melting process in the cavity below the ice shelf (as in Schodlok et al., 2012). The model domain is derived from that of the Estimating the Circulation and Climate of the Ocean, Phase II (ECCO2) project (Menemenlis et al., 2008), but with higher-resolution horizontal grid spacing of ∼ 1 km and 60 vertical levels, as in Schodlok et al. (2012). The bathymetry below the ice shelf is derived in part from NASA Operation IceBridge data (Cochran and Bell, 2012).
The bulk temperature of the ice shelf is determined by first calculating an analytical steady-state temperature profile through the thickness of the ice, taking into account vertical advection and diffusion of heat into the base of the ice as a function of surface and basal ice temperature and calculated basal melting rates (Holland and Jenkins, 1999). This solution requires the assumption that surface accumulation balances basal melting such that the vertical velocity of the shelf is constant (Holland and Jenkins, 1999). This assumption is reasonable for Larsen C, as the median rates of basal melting (Fig. 3a) and surface accumulation  are nearly identical. The surface temperature of the shelf was taken from the Regional Atmospheric Climate Model (RACMO) mean air temperature results from the period 1980-2004 (van den Broeke and van Lipzig, 2004). Although horizontal advection is not explicitly accounted for in the analytical solution, the surface temperature is reduced by 3 • C everywhere to approximately account for the horizontal advection of colder continental ice into the shelf, a tuning that was optimized for the neighboring Larsen B ice shelf (Sandhäger et al., 2005) and which we adopt here for simplicity. The basal ice temperature is assumed to be −2 • C everywhere, approximately the pressure melting temperature of the ice. The resulting analytical temperature profile is then depth-integrated and used to determine the ice rigidity B(T ) following a standard temperature parameterization (Cuffey and Paterson, 2010). We note that due to the nonlinearity in the temperature-rigidity relationship, the rigidity calculated in this manner (B(T )) is not equivalent to the vertically integrated rigidity B(T (z)). However, the difference between the two calculations is small compared to the uncertainty associated with the temperature of the ice, which is the largest overall source of uncertainty in our calculations. It is more straightforward to represent this uncertainty in sensitivity analyses using B(T ± T ), as will be seen below.

Inversion for ice rigidity
We calculate bulk ice rigidity B i using an inverse control method (MacAyeal, 1993;Rommelaere and MacAyeal, 1997) MacAyeal, 1989) that minimizes a cost function measuring the misfit between modeled and observed surface velocity. A partial differential equation constrained optimization algorithm is used that calculates the gradient of the cost function with respect to B and then updates B using a steepest-descent approach. The initial value of rigidity B • provided to the algorithm is calculated as (1−D)B(T ), where B(T ) is the ice rigidity parameterized from the temperature calculation and D is calculated using Eq. (13) with σ b = 0, truncated where necessary to ensure that damage remains within physically acceptable limits (D ∈ [0, 1]). We find that this procedure for specifying the initial state for the control method leads to better fits to the velocity data than providing a constant initial value across the ice shelf, the technique typically applied in previous studies (e.g. Larour et al., 2005;Khazendar et al., 2007Khazendar et al., , 2011.
In an ensemble of model runs starting from different initial states, the average misfit between modeled and observed velocity was minimized when specifying the initial rigidity as described above compared to using any of a wide variety of uniform initial states.

Basal melting and ice shelf temperature
The computed basal melting rates (Fig. 3a) are highest (up to 8 m yr −1 ) near the grounding line, similar to the rates modeled by Holland et al. (2009) near the grounding line. They are also similar to the rates inferred near the grounding line from remote sensing data and modeled surface mass balance by Rignot et al. (2013). The computed mean and median melting rates are 1.3 m yr −1 and 0.7 m yr −1 , respectively, higher than the mean of 0.4 m yr −1 calculated by Rignot et al. (2013). Computed melting is generally higher in the southern part of the shelf, in contrast to the pattern of freezing inferred by Rignot et al. (2013) in the south. The model computes small rates of freezing in the north of the shelf as well as in suture zones downstream of major promontories, similar to the pattern modeled by Holland et al. (2009) and consistent with inferred softer ice in these areas by Khazendar et al. (2011).
We note that, even though Cochran and Bell (2012) measured a maximum gravity anomaly associated with Bawden ice rise, this ice rise is not represented in the bathymetry data set used by the ocean circulation model. For this reason, the melt rates and thermal state of the ice shelf may not be well represented in the vicinity of this ice rise. Rignot et al. (2013) calculated concentrated basal melting of several meters per year near the Bawden ice rise, and though we compute basal melting near Bawden it is not as high nor as concentrated.
The depth-integrated temperature for Larsen C is shown in Fig. 3b. The pattern of temperature follows that of the basal melt rates, as expected. The ice in the southern portion of the shelf is colder, reflecting higher melting rates. The ice is also expected to be colder further south as a result of advection of colder grounded ice into the shelf. We assume an uncertainty of 3 • C in our calculated bulk temperature, which is the standard deviation of the temperature calculation over the entire ice shelf. Given the lack of direct temperature measurements of the ice shelf, this assigned uncertainty is somewhat arbitrary but is chosen more to test the sensitivity of the damage calculation to variations in temperature than to quantify the actual skill of the temperature calculation for representing the true thermal state of the ice shelf.

Damage and backstress
The damage calculated from Eq. (17) is shown in Fig. 4. Figure 4a shows the best estimate for damage using the assumed temperature of the ice (Fig. 3b) to determine the ice rigidity (B(T )). Figure 4b and c show damage assuming an ice shelf with a depth-integrated temperature uniformly warmer and colder by 3 • C, respectively, all else the same. Figure 5a shows the backstress calculated using Eq. (16), with longitudinal strain rates and the parameters α and β taken from modeled velocities resulting from the inverse calculations rather than the observed velocities, which effectively filters out spikes or artefacts in the observational data and produces a smoother backstress field than would be produced if strain rates were calculated directly from observed velocity data. Figure 5b shows a normalized measure of backstress using a parameter analogous to the buttressing parameter f introduced by Dupont and Alley (2005), defined here as Where f = 0, the ice shelf is considered "unbuttressed" and the viscous stress available for thinning of the ice shelf is entirely unopposed. The shelf is considered "fully buttressed" at f = 1, at which point the longitudinal strain rate is zero.
To consider the case of f > 1, where the ice shelf is considered "overbuttressed" and the longitudinal strain rate is compressive (negative), it is necessary to rewrite Eq. (18) by substituting the more general expression for σ b (Eq. 16).

Consequences of diminished backstress
The consequences of losing contact with the Bawden and Gipps ice rises was investigated in model experiments. The inversion for ice rigidity indicated anomalously stiff ice in the vicinity of both ice rises, consistent with the inversion results of Khazendar et al. (2011). The backstress provided by these ice rises was removed by reducing the inverted ice rigidity down to a level appropriate for the assumed temperature of the ice, similar to the masking described by Eq. (15)  within a 60 km radius of the ice rise. The same was done within a 70 km radius of the Gipps ice rise. These radii were chosen based on visual inspection of the inferred rigidity pattern. Elsewhere the original (unmasked) inverted rigidity B i was used. The diagnostic equations for ice shelf creep were then solved using the common SSA approximation (MacAyeal, 1989) to determine the instantaneous speedup of ice flow if contact with these ice rises was abruptly lost. Figure 6a and c shows the impact on the flow of Larsen C if contact with the Bawden ice rise is lost. The flow velocity near the ice rise would increase up to 200 m yr −1 , a 50 % increase, with an area nearly the size of the Larsen B ice shelf experiencing a speedup in excess of 100 m yr −1 , a 25 % increase. A similar order of magnitude change, over a similar but more disjointed area, would occur if the ice shelf lost contact with the Bawden and Gipps ice rises simultaneously ( Fig. 6b and d). Only velocity changes 20 m yr −1 greater than observed are plotted in Fig. 6 to clearly isolate the speedup signal above the "noise" associated with the misfit between modeled and observed velocity in the inversion solution for the ice rigidity. Nearly all of the nodes in the model (95 %) had a misfit of 20 m yr −1 or less, so we only focus on velocity change greater than this in Fig. 6.   The Cryosphere, 7, 1931-1947, 2013 www

Thermal state of the ice shelf
The spatial pattern of temperature calculated for the shelf appears to capture the tributary structure of the ice shelf to the extent that the bathymetry beneath the ice shelf (which influences the ocean circulation and thus melting rate pattern) represents an extension of the bedrock structure upstream of the grounding line that directs the flow into the shelf. We note that the presence of narrow bands of marine ice identified by Holland et al. (2009) are not captured in the steadystate temperature calculations here, at least not beyond areas where the ocean circulation model indicates direct freezing. It would be reasonable to expect warmer depth-averaged temperature along longitudinal bands that contain marine ice, therefore we may be overestimating damage in areas where marine ice is present. While it would be possible to manually prescribe warmer temperature along these inferred longitudinal bands, our primary purpose here is to calculate backstress and damage given an independent temperature estimate for the shelf and demonstrate sensitivity of the results to variations in temperature. The large uncertainty we assign to our calculation for ice temperature reflects both the absence of actual temperature measurements against which to validate our results as well as the fact that the calculations neglect horizontal advection. The thermal state of the ice advected horizontally across the grounding line will be important in determining the temperature of the ice shelf. We may be overestimating the tem-perature in some areas and underestimating it in others, depending on the rate of basal melting or freezing that in turn determines the thickness of the thermal boundary layer at the base of the ice column (e.g. Sergienko et al., 2013).
For a 3 • C uncertainty in the bulk temperature of the ice, the ice rigidity B(T ) varies in the range of 8-29 %, with a mean of 15 %. This uncertainty likely dominates over that in the inverted ice rigidity B i , therefore from Eq. (17) the mean uncertainty in damage is 15 %. The uncertainty is higher for the warmer areas of the shelf owing to the nonlinearity in the temperature-rigidity relationship as the melting temperature is approached (Cuffey and Paterson, 2010). Thus much of the northern half of the shelf as well as areas downstream of the major promontories of the shelf where marine ice accumulates have higher uncertainty, whereas colder meteoric flow units have less. The fact that the uncertainty in damage is directly proportional to, and dominated by, the uncertainty in ice temperature suggests that field studies of ice shelves should include measurements of the ice temperature whenever possible. Better knowledge of the thermal state of ice shelves will be necessary to distinguish between the influence of temperature vs. fracture in observations of ice flow and for making predictions of the mechanical integrity of ice shelves in the future.

Comparison of damage with Operation IceBridge altimetry
In Fig. 4, many of the rifts visible in the MOA image appear damaged. For those rifts that are damaged, the spatial agreement between damage and rift location is good. However, not all rifts are damaged, and not all damaged rifts are equally damaged. Numerous rifts have high levels of damage, approaching the effect of traction free cracks with vanishing viscosity (D → 1) in a continuum sense. Other rifts have intermediate damage, and some appear entirely undamaged and are deforming coherently and indistinguishably from the surrounding ice shelf. The explanation for this behavior likely has to do with the stabilizing influence of mélange filling in between the rift flanks (Khazendar and Jenkins, 2003). Figure 7a shows a close up of an area of rifts that originate downstream of Hollick-Kenyon Peninsula and proceed toward the Gipps ice rise. The black line indicates a NASA Operation IceBridge flight in 2009. ATM surface elevation data from this flight are plotted in Fig. 7b and compared against the damage calculation along the same track. The ATM data have a 40 m spacing along track and an 80 m sample width, and are plotted as elevation above the GL04c geoid for consistency with the Bedmap2 thickness data (Fretwell et al., 2013)   , which shows damage (red, left) and IceBridge ATM surface elevation above mean sea level (black, right) using the GL04c geoid (Fretwell et al., 2013).
The rifts in Fig. 7 show varying levels of damage. Damage decreases and the surface elevation between the rift flanks increases for rifts labeled 1 to 3 in Fig. 7. Rift 1 is the youngest of the transect and has relatively little mélange to stabilize it. With age and further advection with the shelf, the surface elevation (a proxy for mélange thickness) increases for Rifts 2 and 3 and subsequently the damage decreases. Rift 4 appears only moderately damaged despite the lowest surface elevation recorded along the transect. This can probably be explained by the kinked shape of this fracture and some measure of stress shielding given its proximity to Rift 5, the most damaged rift and the rift that is currently interacting with the Gipps ice rise. Rift 6 is the oldest of the transect and has the thickest mélange of all, contributing to its behavior as if completely healed.

Damage patterns
Aside from rifts, damage is consistently highest near the grounding line, presumably the result of the large stresses -and thus fracture -induced as the ice shelf adjusts toward flotation. Areas of damaged ice extend from the grounding line for long distances into the ice shelf in several places. The patterns of damage downstream of the major promontories of the shelf are similar to the patterns of softer ice inferred by Khazendar et al. (2011). The magnitude of damage in these areas is lower than the damage associated with rifts, an expected result since the crevasses in these areas, primarily basal crevasses, occupy only a fraction of the full ice thickness (McGrath et al., 2012;Luckman et al., 2012).
Damage downstream of promontories decreases with distance into the shelf, in many places reducing to zero before reaching the ice front. This indicates that the influence of fractures on longitudinal deformation is diminishing with distance from the grounding line. This is perhaps a counterintuitive result, as longitudinal traces of fractures are visible from the grounding line all the way to the ice front in many parts of the ice shelf (Glasser et al., 2009). However, this pattern is consistent with inverse method results for ice rigidity for Larsen C (Khazendar et al., 2011) as well as Larsen B (Khazendar et al., 2007;Vieli et al., 2007) that show softer ice in the lee of promontories transitioning into stiffer ice closer to the ice front. A direct inversion for damage on Larsen B also showed a similar pattern (Borstad et al., 2012).
There are several possible explanations for this behavior. The first possibility is that the calculated temperature field for the ice shelf is incorrect. Since the damage pattern is directly linked to the assumed temperature of the ice, any errors in temperature will introduce errors in the magnitude and pattern of damage. The steady-state temperature calculation used here ignores horizontal advection, and no in situ temperature data exist to check the validity of our calculations. The true temperature field of the ice shelf may have a more banded structure longitudinally as a result of horizontal advection, with a bulk temperature that would depend strongly on the presence and thickness of a marine ice layer. To produce damage more continuously from the grounding line to the ice front along trains of known crevasses, the temperature of the ice would need to grow progressively colder -more so than we already calculate -moving away from the grounding line.
If the calculated temperature of the ice shelf is reasonable, it is still possible for the influence of fractures to diminish as they advect with the shelf. In some areas, marine ice is likely accumulating at the tips of basal crevasses and progressively healing the fractures as they advect with the shelf. However, this vanishing influence of fractures does not necessarily imply that the fractures are fully healing and disappearing. When a fracture first forms, it contributes to damage in two ways; first, by reducing the available load bearing surface area, and second, by causing a stress concentration at the tip of the fracture (Murakami and Ohno, 1981). Both mechanisms would contribute to increasing the mean stress over the remaining load bearing cross section. When a basal crevasse forms near the grounding line, for example, the fracture tip is initially sharp and therefore both contributions to damage are present. With subsequent advection into the shelf, the crack tip is blunted by creep deformation and the stress concentration ahead of the crack tip relaxes until eventually only the The Cryosphere, 7, 1931-1947, 2013 www.the-cryosphere.net/7/1931/2013/ area reduction contributes to damage. Thus the inferred damage would decrease even if the absolute crack length does not change.
A further possibility for explaining undamaged areas of the ice shelf where fractures are known or assumed to be present is that the definition of the effective stress given by Eq. (8) is inappropriate or incomplete. According to Eq. (8), at failure (D = 1) the stress (and thus strain rate) tend to infinity, whereas in creep rupture experiments they tend to finite values (Murakami and Ohno, 1981). One method to account for this observation is to define two different effective stress measures, one for use in the evolution of damage and one for use in the creep relation. In this approach, an effective stress of the form of Eq. (8) is used for damage evolution, whereas an effective stress with D replaced by cD where 0 < c < 1 is used for the creep relation (e.g. Eq. 10), ensuring finite stress and strain rate at failure (Murakami and Ohno, 1981). This would explicitly distinguish between a reduced level of damage that represents the influence of fractures vs. the full damage that represents the presence of fractures. The damage solution in Fig. 4 would represent the former, whereas thus the full damage representing the latter may link more continuously to the ice front in a manner consistent with observations of longitudinal trains of fractures (e.g. Glasser et al., 2009).

Backstress patterns
The magnitude and pattern of backstress for Larsen C (Fig. 5), with higher backstress near the grounding line, at the confluence of neighboring tributary glaciers, and upstream of ice rises, is consistent with findings for other ice shelves (Thomas, 1973b;Thomas and MacAyeal, 1982;Kenneally and Hughes, 2004). The Bawden and Gipps ice rises contribute locally high backstress to the shelf. Just upstream of these ice rises, the ice shelf is overbuttressed (f > 1) and the flow is compressive longitudinally. The high backstress associated with these ice rises is also consistent with the high values of ice rigidity inferred in these locations by Khazendar et al. (2011). Approaching the McDonald ice rumples on the Brunt ice shelf, Thomas (1973b) found that the backstress decreased approximately inversely with distance from the ice rumples (σ b ∝ 1/r 1.07 ). Along a flowline intersecting the Bawden ice rise, we find a nearly equivalent relationship (σ b ∝ 1/r 1.05 ).
The backstress associated with the Gipps ice rise covers a larger area than that of the Bawden ice rise (Fig. 5). The ice in this part of the shelf is more heterogeneous and fractured, with a mélange of icebergs and sea ice filling in the open areas to the east of Hollick-Kenyon Peninsula (e.g. Fig. 7b). In addition to encountering the longitudinal resistance of the ice rise, the flow is also laterally compressive as it approaches the ice rise. For these reasons, it is not possible to discern a similar relationship of backstress vs. distance for the Gipps ice rise.
Another noteworthy aspect of Fig. 5 is the backstress associated with the confluence of neighboring tributary glaciers in the ice shelf. Elevated backstress is evident in many areas where streamlines converge downstream of promontories, a result of lateral compression as neighboring glaciers flow together. Thomas (1973b) described this as a "bottleneck" effect, finding a similar pattern of backstress on the Amery ice shelf. In Cabinet Inlet, the high backstress is likely the result of both shear stress along Churchill Peninsula and with a sort of bottleneck effect as the ice approaches the end of the peninsula and turns the corner to begin flowing toward the Bawden ice rise.
Isolated patches of grounded ice near the grounding line determined from differential InSAR data (Rignot et al., 2011a) are highlighted in red in Fig. 5. Just beyond the end of both Churchill and Hollick-Kenyon Peninsulas, small grounded spots coincide with the origin of rifts. Near the grounding line in a few places, these grounded spots appear to contribute to the inferred backstress. In many places these spots coincide closely with the locations of basal crevasse formation, as evidenced by the longitudinal features that originate at these locations. These small areas of grounded ice may be very important to the flow and the stability of the ice shelf if they are responsible for the kinds of longitudinal features that coincide with, and are typically thought to be the cause of, the termination of rift propagation (e.g. Glasser et al., 2009).

Weakening contact with Bawden ice rise?
Given its small size and relief, the Bawden ice rise appears to be the most likely location for the ice shelf to become ungrounded and destabilized in the near future. If the ice shelf were to lose contact with this ice rise, the flow of the entire northern half of the ice shelf would be affected ( Fig. 6a and c). The ice adjacent to the Bawden ice rise is already quite damaged (Fig. 4a). Over the period 2003-2008 the thickness change of the northern part of the ice shelf was small , though the flow velocity increased by 15 % or more between 2000 and 2008 (Khazendar et al., 2011). Progressive weakening of the ice in the vicinity of the ice rise, and associated reduction in backstress transmitted upstream, may be a contributing factor in this acceleration. This could indicate that the ice rise is transitioning from acting to stabilize the ice shelf to acting as an "indentation wedge" that generates or exacerbates fractures, possibly contributing to the eventual destabilization of the ice front as witnessed for other ice shelves (e.g. Doake and Vaughan, 1991).
Further weakening of the ice adjacent to the ice rise would diminish the backstress transmitted upstream, with significant consequences for the local stability of the shelf. An abrupt velocity increase of the type modeled in Fig. 6 could destabilize the rifts that are already present in the vicinity of Jason Peninsula, possibly leading them to propagate southward into the ice shelf. These rifts currently terminate along a band of marine ice identified by Holland et al. (2009). If the rifts were to remain stabilized by this marine ice after the stress perturbation associated with loss of contact with the Bawden ice rise, then the northern part of the ice shelf may remain intact. However, if the stress was sufficient for the rifts to propagate through this marine ice band, these rifts could link up with the wide train of basal crevasses that originate at the tip of Churchill Peninsula and extend all the way to the ice front (McGrath et al., 2012). This could potentially cause a large portion of the northern ice shelf to calve off, approaching the size of the 2002 Larsen B collapse. Key questions in this scenario are how fractures are stabilized along longitudinal features that contain marine ice (e.g. Glasser et al., 2009) and the limits to this stabilizing influence, questions which remain open. While the consequences of ungrounding of the small Bawden ice rise are speculative at present, the ice rise is a significant pinning point for the Larsen C ice shelf and deserves further monitoring and study.
Note that the grounding line is not affected in the perturbation simulations presented here. This is likely the result of the high levels of backstress near the grounding line associated with lateral convergence and shearing along promontories (Fig. 5). A backstress perturbation originating at the ice front, at least for the configuration of Larsen C, does not appear to matter much to the grounding line. However, this may not be a generalizable result, as some ice shelves likely have pinning points that connect more directly to the grounding line.

Damage evolution
In addition to providing a straightforward framework for monitoring the mechanical integrity and buttressing capacity of ice shelves using remote sensing data, the types of damage calculations presented here will serve as important benchmarks for testing and optimizing dynamic damage evolution models. Dynamic damage models have been developed and calibrated at smaller spatial and temporal scales, and their introduction into large-scale ice sheet models will eventually allow for modeling calving, rifting, retreat and disintegration of ice shelves. The damage mechanics formulation used in the present paper is derived from only the momentum balance equations for diagnostic calculations of ice velocity. Simulating the advection and evolution of damage in time requires a separate dynamic damage function coupled with the equations for heat and mass balance for an ice shelf. A number of established functional frameworks for modeling damage evolution have been applied to polycrystalline ice (e.g. Pralong and Funk, 2005;Duddu and Waisman, 2012) as well as analogous high homologous temperature materials that obey power law creep (e.g. Lemaitre, 1996). Though the validity of such formulations at the scale of modeling an entire ice shelf remains to be confirmed, the benchmarks provided by the methods in this study will aid with the devel-opment and testing of such an evolution model within the framework of a large-scale ice sheet model.

Conclusions
The analytical theory for the creep deformation of a floating ice shelf was extended using continuum damage mechanics to account for the softening influence of fractures on longitudinal deformation. A scalar damage variable was introduced to represent the influence of fractures on an ice shelf while maintaining a continuum representation of the material physics, an advantage for implementation in largescale models. This framework is fully generalizable to any ice shelf, with resistive backstress accounting for areas where the strain rate is lower than the limit for an unconfined ice shelf and damage accounting for areas of fracture-induced softening. A control method inversion for ice rigidity combined with the damage theory allows independent damage and backstress calculations for an ice shelf.
Using field and remote sensing data, this new framework enables monitoring the structural integrity of ice shelves, their ability to buttress the flow of ice at the grounding line, and thus their indirect contribution to ice sheet mass balance and global sea level. For the Larsen C ice shelf, we find damage in areas with visible crevasses and rifts. The level of damage of rifts is modulated by the thickness of the mélange that fills in between the rift flanks, with the thickest mélange reducing damage to zero in older rifts. The Bawden and Gipps ice rises appear to be critical pinning points that stabilize the ice front in its current configuration. The ice in contact with the Bawden ice rise is mechanically weak, and further weakening would lead to a significant increase in flow velocity and possible destabilization of the northern part of the ice shelf.