Introduction
Long, >100 km, transverse rifts near the fronts of ice shelves become the planes along which tabular icebergs calve. Calving is an important ice mass loss mechanism around the Antarctic. The large tabular icebergs A-43 and A-44, which calved from the western front of the Ronne Ice Shelf (RIS), Antarctica, in May 2000, contained an ice mass of ∼2075 × 1012 kg, approximately equivalent to 1 year of accumulation across the entire Antarctic, 2326 × 1012 kg a−1 (Reference Giovinetto and ZwallyGiovinetto and Zwally, 2000). In the RIS (Figs 1 and 2) such rifts are observed to begin as shear-margin fractures along the lateral boundaries of outlet streams feeding the shelf. Fractures may propagate laterally as they are advected through the shelf, responding to changes in the local stress field. The local stress field depends on the glaciological (far-field) stress, the geometry of the fracture and the geometries of nearby fractures. Inhomogeneity of ice material properties may also play a role. It may be important to understand how relatively small fractures evolve into large ice-shelf rifts so that the calving mass loss mechanism may be accounted for in predictive ice-sheet models, either directly or by parameterization.
Feature map digitized from the MOA. Evans Ice Stream tributaries bounded by suture zones observable in the MOA are labeled for reference in the text.
Fractures in the EIS outflow. Locations of experiments described are indicated by number on the right. Upstream fractures are inferred to be early evolutions of large rifts near the shelf front (Reference Fahnestock, Scambos, Bindschadler and KvaranFahnestock and others, 2000).
Existing calving criteria, developed for tidewater calving glaciers, rely on empirical correlations between quantities such as water depth and front location (Reference MeierMeier and others, 1994; Reference Van der VeenVan der Veen, 1996) with an incomplete representation of the physical processes involved in iceberg production. Similar developments for the large, embayment-filling ice shelves associated with marine ice sheets remain elusive. At present, calving criteria tend to focus on longitudinal strain rates, an approach that implicitly emphasizes vertical propagation (Reference ReehReeh, 1968; Reference Fastook and SchmidtFastook and Schmidt, 1982; Reference HughesHughes, 1983; Reference AlleyAlley and others, 2008; Reference Nick, van der Veen and VieliNick and others, 2008). The work reported here suggests a different view is warranted, one that embraces fracture propagation in the horizontal plane, at least in the case of large embayed ice shelves, such as the Ronne and Ross Ice Shelves.
Fractures are an observable effect of stress in a material. They initiate at pre-existing cracks or flaws and propagate when the stress intensity at a crack tip is large enough to overcome the fracture toughness of the material. Once propagation is initiated, it continues until the stress intensity at the propagating tip falls below the value required to overcome the fracture toughness. Considering only fracture geometry and stress intensity, the longer a fracture becomes, the more likely it is to continue propagating in a given remote stress field (Reference LawnLawn, 1993). Yet fractures in glacier ice (and other geologic materials) are observed to have finite lengths. The explanation for this is often that interactions between adjacent fractures reduce stress intensity at the fracture tip and limit propagation. Interaction might be in the region of adjacent tips (Reference Cotterell and RiceCotterell and Rice, 1980) or through the development of a stress shadow cast by a larger fracture (Reference Pollard and AydinPollard and Aydin, 1988). Alternatively, the fracture may be driven by an internal pressure (i.e. a fluid), in which case the fluid pressure drops as the fracture extends. In the RIS, spatial variation in material properties may also be important.
Crack tip arrest is often observed to coincide with structural boundaries, such as suture zones between ice from adjacent ice-stream tributaries (visual inspection of Figs 1 and 2; also Reference LeDouxLeDoux, 2007; Reference GlasserGlasser and others, 2009; Reference Holland, Corr, Vaughan, Jenkins and SkvarcaHolland and others, 2009). Here we address the observed relationship quantitatively. There are two possible roles for suture zones in crack tip arrest. First, enhanced shear within these zones may reduce the mode I stress intensity at fracture tips. Second, past shear experienced by this ice may modify its mechanical properties (fracture toughness, elasticity constants). Near the surface, cold air ponds in crevasses, cooling the ice and making it more brittle, while at depth the ice may be relatively warm and less brittle due to its shear history (Reference Harrison, Echelmeyer and LarsenHarrison and others, 1998). The ice fabric may have a preferred crystal orientation due to its strain history. Once a crack tip has arrested at a suture zone, microfracturing in the process zone (to dissipate energy when propagation is not possible) may further modify ice rheology (Reference BassisBassis and others, 2007). Propagation may reinitiate in response to that modification or by advection into a more favorable stress field. Where fractures are observed to propagate through a structural boundary, they often go on to arrest at another boundary. To investigate fracture evolution and evaluate the role of structural boundaries in that process, we model fracture propagation in the ice shelf using linear elastic fracture-mechanics theory and a boundary element numerical method.
Feature Map of the Ronne Ice Shelf
For this study we digitized a feature map of the Ronne–Filchner Ice Shelf using the MODIS Mosaic of Antarctica (MOA) image map (T. Haran and others, nsidc.org/data/nsidc-0280.html). The mosaic is a composite of individual MODIS (Moderate Resolution Imaging Spectroradiometer) images taken between 20 November 2003 and 29 February 2004, with an effective resolution of 125 m. The range of sun illumination angles in the composite image makes it possible to trace many features, including large fractures, streak lines and disturbed ice associated with relict ice-stream shear margins now advecting through the shelf. Fractures appear either as sharp, sunlit/shadowed faces or as shadows cast within sagging snow bridges (Reference Merry and WhillansMerry and Whillans, 1993). Small fractures, or those well covered by snow, are not visible. The closely spaced fractures in a chaotic shear margin leave a signature in the shadows cast along the margin track (Reference Merry and WhillansMerry and Whillans, 1993). Slope breaks, such as at the coasts of ice rises and at the grounding line, are also visible as shadows. Features in the map presented here were digitized at the dark/light boundary.
In this work we examine fractures advecting downstream from the outlet of Evans Ice Stream (EIS), in the western region of the RIS (Figs 1 and 2). The orientations of streak lines and relict shear margins in this region agree well with modern velocity azimuths, indicating steady flow over the lifetimes of features embedded within the flow. Thus, change in fracture geometry with distance downstream may be equated with fracture evolution.
Suture zones that form where tributaries to the EIS outlet merge persist into the ice shelf and can be traced to the shelf front (we identify these sutures by letter for ease of reference; Fig. 1). Each suture zone is characterized by a lateral change in fracture pattern. The regional character of fracture patterns in the EIS outflow changes as the ice advects from the grounding line to the shelf front. These changes may be due to rotation of fractures with respect to principal stress directions, changes in the stress field itself and, perhaps, also inhomogeneity in ice properties. Fractures originating from suture zone D, in the easternmost part of the EIS outflow, advect downstream to become large transverse rifts near the shelf front.
Four major fracture-pattern changes take place with distance downstream along suture zone D (Fig. 2). (1) Immediately downstream of the grounding line is a visually complicated region characterized by arcuate crevasses, intersecting crevasses and relict shear margins. The shear margin fractures become less prominent with distance from the grounding line. (2) Downstream of Fowler Peninsula, closely spaced (1–2 km) upstream-pointing fractures develop, overprinting older crevasses. These fractures arrest where they approach other upstream-pointing fractures associated with a suture zone in ice from an adjacent EIS tributary. (3) About 160 km downstream of Fowler Peninsula fractures propagate east, normal to the suture zone, arresting at the suture zone between ice from EIS and from Carlson Inlet. (4) Finally, ∼200 km from the shelf front, suture zone D is breached by a small number of fractures that propagate west as extensions of the fractures that arrested their eastward propagation at the Evans/Carlson suture zone. Near the shelf front, the longest fractures exceed 100 km.
Fracture-Mechanics Model
Theory
Linear elastic fracture mechanics provides the physical and mathematical framework for studying fracture propagation in the ice shelf. Near the ice surface on timescales of days to months, glacier ice behaves elastically and undergoes brittle fracture at sufficiently large stresses (cf. Reference Nath and VaughanNath and Vaughan, 2003; Reference Sergienko, MacAyeal and BindschadlerSergienko and others, 2009). Fracture propagation is treated as a mixed-mode opening (mode I) and sliding (mode II) phenomenon. In the present work, test fractures that mimic observed geometries are investigated using a remote (glaciological) stress field derived from observed ice-shelf flow. We compute equilibrium fracture geometry for a given stress state and do not compute the rate of propagation.
We are interested in the horizontal propagation of long fractures in the ice shelf, so we take advantage of the large ratio between horizontal and vertical length scales (tens of kilometers to <1 km) and assume a plane-strain condition. This requires that fracture geometry, material properties and boundary conditions, including boundary tractions, do not vary in the vertical direction, which, in turn, requires us to assume that material properties are constant across a fracture and that tidal displacement does not affect one side of a fracture more than the other. We consider these assumptions reasonable, except in the near-shelf-front environment.
The governing equations for plane strain in an elastic material are presented by many workers (e.g. Reference MuskhelishviliMuskhelishvili, 1963). Here we provide a brief overview. The static stress equilibrium equations are
in which the horizontal coordinates i, j ∈ x, y are used and the stress tensor, σij , is symmetric. The constitutive relationship between stresses and strains, εij , for an isotropic linear elastic material under plane strain (Hooke’s law) is
which δij is the Kronecker delta and the repeated indices in εkk imply summation. The shear modulus, μ, and Poisson’s ratio, v, are independent. In the case of an incompressible material, Poisson’s ratio is 0.5 and Equation (2) is simpler.
The compatibility equation, which ensures continuity following deformation, is
Writing Equation (3) in terms of stresses requires a rewritten stress equilibrium
and the definition of linear elasticity (Equation (2)). The result is
It is convenient to rewrite the compatibility equation using a scalar quantity. The Airy stress function, ψ(x, y), is chosen such that
in which case, the biharmonic,
is produced. The solution to the plane-strain problem lies in finding an appropriate function, ψ(x, y).
A fracture will propagate only when the loading stress at one or both tips is large enough to overcome the toughness of the material. Thus, to determine whether or not a fracture will propagate, and , if so, the propagation direction, we must be able to evaluate the stress field near the tip. This, in turn, depends on the geometry of the crack tip. Assuming a ‘sharp slit’ geometry, we follow the stress intensity approach developed by Reference IrwinIrwin (1957) and expanded by many others (Reference Erdogan and SihErdogan and Sih, 1963; Reference Paris and SihParis and Sih, 1965; Reference Kanninen and PopelarKanninen and Popelar, 1985; Reference Pollard, Segall and AtkinsonPollard and Segall, 1987; Reference LawnLawn, 1993). Close to the fracture tip (but at a distance slightly above 0 to avoid the singularity), stresses at an orientation θ relative to the fracture plane and a distance r from the fracture tip are found to vary as r −1/2 and are amplified by a scalar quantity called the stress intensity factor (Reference Rice and LiebowitzRice, 1968):
in which Km is the stress intensity factor for mode m. Reference Pollard, Segall and AtkinsonPollard and Segall (1987) found the r −1/2 relationship to be suitable for r <10% the fracture half-length, a. Using the approach of Reference Erdogan and SihErdogan and Sih (1963), the stress intensity factor may be written in terms of the loading stress, σm :
in which the constant γ is defined by the fracture geometry. Superposition allows addition of stress intensity factors for a given mode.
The loading stress required in Equation (9) is not straightforward for many geometries, so displacement near the fracture tip, which varies as r 1/2, is used instead. Maximum displacement, D maxm , near an active (propagating) tip is related to the loading stress for mode m:
in which E is Young’s modulus. If the maximum displacement can be evaluated, then Equation (10) may be used together with Equation (9) to compute the stress intensity at the fracture tip:
Here, our interest is mixed-mode propagation. We employ the maximum circumferential stress criterion, which assumes that the fracture propagates at an angle θ from the tip in a direction normal to the maximum circumferential stress,
(Reference Erdogan and SihErdogan and Sih, 1963; Reference Ingraffea and AtkinsonIngraffea, 1987):
in which the subscript on K indicates the mode of loading. For propagation, the quantity
in Equation (12) must be equal to a critical stress intensity factor at which the material fails, known as the fracture toughness, K
IC. Dividing Equation (12) by the fracture toughness yields a propagation criterion
When the right-hand side of Equation (13) is ≥1, propagation occurs. The angle θ 0 is determined using the condition that propagation is along the plane in which shear stress is zero
which yields a useful solution:
Analytical solutions for in-plane stresses and stress intensity factors are available for straightforward geometries and loadings of individual and sets of fractures (e.g. Reference SihSih, 1973; Reference Tada, Paris and IrwinTada and others, 2000). A numerical solution is required for the more complicated scenarios of interest here (cf. Reference Crouch and StarfieldCrouch and Starfield, 1983).
Numerical model
We use the displacement discontinuity boundary element method (BEM) to solve the elastic boundary value problem (see Appendix) and have modified the ‘Frac2D’ codes developed by A.L. Thomas (unpublished information, 1991) for this purpose. The method allows easy computation of stresses in the region of a fracture or set of fractures using their associated displacements and evaluates whether or not a given fracture will propagate and, if so, the propagation direction. Boundary stresses are simulated by evaluating the cumulative (and unknown) influence of the shear and normal displacement discontinuities of all elements on each element. Once boundary stresses are computed, stresses can be evaluated for any location within the model domain, defined using an observation grid.
We validated our implementation of the numerical scheme using a set of experiments with increasingly complicated boundary conditions (Reference LeDouxLeDoux, 2007). A simple test is shown here, in which a single fracture is placed within a small box boundary and stresses are simulated. A uniform left-lateral shear stress, σxx = 0, σyy = 0, σxy = −0.5 MPa, yields an expected stress field near the fracture (Fig. 3). Fracture propagation is computed from the simulated stresses. This approach is widely used in geomechanics (e.g. Reference Sempere and MacdonaldSempere and Macdonald, 1986; Reference Olson and PollardOlson and Pollard, 1989; Reference DahmDahm, 2000).
Mean stresses in a box under left-lateral shear, σxx = 0, σyy = 0, σxy = −0:5 MPa, with one fracture boundary. The contour interval is 10 MPa; a dashed contour indicates negative values.
Model implementation
In each experiment, a model domain is defined that includes fractures of interest and a region boundary that allows the remote stress field to be well resolved. Starter fractures may trace existing fractures or may be simplified for the purpose of testing a specific hypothesis. Simple rectangular geometries were used for the region boundaries. The element length is 1–2 km for the region boundary and 0.5–1 km for fracture boundaries, reflecting their more complicated geometries.
Remote stresses and elasticity constants are required to initialize the model for each experiment. The remote (glaciological) stresses used to specify boundary conditions are computed on a 1 km grid using the inverse form of Glen’s flow law, with a uniform rate factor of 760 kPa a1/3 (an average value for the western part of the RIS; Reference Larour, Rignot, Joughin and AubryLarour and others, 2005). Surface velocities from interferometric synthetic aperture radar (InSAR) (personal communication from I. Reference Larour, Rignot, Joughin and AubryJoughin, 2005) are used to compute infinitesimal strain rates as gradients of the velocity field. The original velocity data are smoothed using a running average in a 20 km square in order to remove a small-wavelength artifact in one velocity component that is enhanced by the strain-rate calculation. Ice thickness is from the BEDMAP dataset (Reference Lythe and VaughanLythe and others, 2000). Deviatoric stress components are computed following the engineering hydrodynamics practice of taking two-thirds the ice thickness, H (Reference Lambe and WhitmanLambe and Whitman, 1969), so the overburden pressure used is 2/3 ρgH (Fig. 4). Regional stresses, M (Table 1), averaged from the region boundary are prescribed as parameters to the model, so they represent a constant remote-stress value at all locations within the model domain.
(a) Observed surface speed, with contour interval 100 m a−1, and (b) computed principal deviatoric stresses in the EIS outflow. The contour interval is 100 m a−1. Lighter black curves in the stress map indicate extension, and heavier gray curves indicate compression. Sutures and fracture locations are from Figure 1. Locations of experiments described here are indicated by number.
Large-scale remote stresses and coefficients for the model experiments. M are regional stresses,
are boundary scalings, v is Poisson’s ratio, E is Young’s modulus and K
IC is the fracture toughness
Boundary stresses are initialized using the stress field derived from observation. The initialization is performed with only the region boundary; fracture boundaries are added to the domain after this procedure. Principal stresses are scaled and elasticity constants are tuned for each region boundary, so gridded glaciological stresses resemble ice-shelf stresses (on a 500 m grid), according to a root-meansquare comparison between model- and velocity-derived values of mean stresses, computed as (σxx + σyy )/2. Scaling of the boundary stresses is always required to obtain a satisfactory fit. Calibrated elasticity constants used in various model domains range between 7500 and 8000 MPa for the modulus of elasticity and between 0.27 and 0.29 for Poisson’s ratio, similar to values reported in the literature (Reference MellorMellor, 1975; Reference HutterHutter, 1983; Reference RistRist and others, 1999, Reference Rist, Sammonds, Oerter and Doake2002; Reference King and JarvisKing and Jarvis, 2007).
Once a region boundary is calibrated, a fracture boundary (or boundaries) is added to the domain and the effect of the fracture(s) on the stress field is evaluated. Principal stresses for the fracture boundary may be further scaled to generate propagation. We evaluate this requirement by shortening our starter ‘test’ fractures relative to observed fractures and tuning the fracture boundary scaling such that the observed propagation is reproduced. The tuning parameters for region,
, and fracture,
, are scalar quantities so principal stress directions, and thus fracture propagation directions, are not affected by this procedure. In any given experiment, fracture tips may be active (allowed to propagate) or held fixed, depending on the topic being addressed. In the simple, single-mode case, fracture propagation is evaluated by comparison of the modeled stress intensity factor, Km
(Equation (11)), at the fracture tip and the fracture toughness, K
IC. Where the stress intensity factor is >K
IC, the tip must propagate. For mixed-mode propagation, fracture growth depends on K
I, K
II and K
IC via the propagation criterion in Equation (13). Estimates of fracture toughness for glacier ice range from ∼0.1 to 0.3 Mpa m1/2 (Reference Rist, Sammonds, Oerter and DoakeRist and others, 2002); a value of 0.3 Mpa m1/2 is used in all but one of the experiments reported here.
When propagation of a fracture tip is indicated, a boundary element is concatenated to the fracture boundary with an appropriate orientation (Equation (15)) and length computed as half the length of the boundary element. Following the addition of a new element, boundary stresses, and thus stress intensity and propagation direction, may change. If the stress field is unchanged, each fracture increment begins at the previous element tip and grows in the same orientation. If the stress field changes, the orientation of the fracture increment adjusts to minimize shear loading. The iterative process is used only to predict fracture geometry; there is no time component in the model. A range of fracture geometries are produced during model iteration, including no propagation, propagation that extends to the edge of the model domain (the region boundary), and episodic arrest and reactivation due to change in the stress field via change in fracture length, by propagation of the opposite tip.
Experiments
The importance of structural boundaries in crack tip arrest is investigated by modeling the propagation behavior of fractures at distinct steps in their downstream evolution toward the calving front. We create a set of test fractures located in transitional zones from one dominant fracture pattern to another, using observed fracture geometries. In each case, the model domain is initialized as described above. Model results are plotted using contoured mean stresses in order to clearly demonstrate how the presence of the fracture modifies the stress field and concentrates stresses at the tips of a mixed-mode fracture.
As we do not have any particular insight into the physical characteristics of suture zones with respect to stress intensity, we do not attempt to resolve these features and instead pose the question: would a fracture with a specified geometry in a given remote stress field propagate were a structural boundary not present? If the model fracture propagates but we do not observe propagation in the ice shelf, then we may conclude that the suture zone plays an important role in observed crack tip arrest. If the fracture does not grow, then the stress conditions must not be favorable for propagation and we may draw no conclusion about the role of the suture zone.
Experiment 1
Closely spaced, 1–2 km apart, subparallel fractures develop along suture zone D downstream of Fowler Peninsula (Figs 1 and 2). The fractures have a distinct geometry, with upstream-pointing left-lateral (relative to ice flow) segments. As they advect downstream, they develop short, 1 km, transverse right-lateral segments. The left-lateral segments of the fractures rotate slightly and their upstream-pointing segments become less visible in the MOA. About 160 km downstream from Fowler Peninsula, the transverse segments of these fractures propagate eastward, arresting at the suture zone between ice from EIS and Carlson Inlet. The propagation direction is similar to the orientation of the original transverse segments. With increasing distance downstream, longer fractures remain visible while shorter fractures become less visible in the MOA, suggesting that the features are inactive and filling with snow.
The model is used to test the hypothesis that tip arrest is due to the presence of the Evans/Carlson suture zone and to investigate interaction between adjacent fractures. Relatively large compressive stresses at the eastern (right-lateral) tips of the test fractures in the remote stress field suggest that conditions in this region may be favorable for propagation. Downstream of the experiment location, fractures with similar geometries are observed to have propagated eastward and arrested at the suture zone.
Propagation of two adjacent test fractures, 1α and 1β, is examined. The geometry of test fracture 1α is identical to an observed fracture. Test fracture 1β, 3.4 km downstream of test fracture 1α, also follows an observed fracture trace but is shortened by 5 km at its eastern tip. Shortening allows us to evaluate the model’s ability to reproduce the observed geometry, as well as test the hypothesis that the suture zone causes tip arrest. Tips of the upstream-pointing fracture segments are set to be inactive. This allows the complete fracture geometry to be used in computing near-field stresses while maintaining the non-propagating tips at their observed locations.
Calibration follows the procedure outlined above (Table 1; Fig. 5). We evaluate scaling of the fracture boundary stresses by determining the value,
, required to initiate propagation for a given value of the fracture toughness, K
IC. Using K
IC = 0.3 Mpa m1/2 and
, the right-lateral end of fracture boundary 1β propagates once to relieve shear stress, producing a characteristic kink geometry, and then arrests (Fig. 6b).
Experiment 1. (a) Mean stresses from observation and (b, c) two different calibrations of elastic properties for the domain using
. The contour interval is 50 kPa. The preferred calibration (b) uses v = 0.29 and E = 8000 MPa. The calibration in (c) uses v = 0.28 and E = 8000 MPa and produces an unobserved gradient from top to bottom of the domain. Fractures are displayed to orient the reader, but are not included in stress calibrations. Model simulated stresses are interpolated to a 500 m grid for comparison with stresses computed from the observed velocity field. The root-mean-square difference in the mean stress is 85 kPa between (a) and (b) and is 185 kPa between (a) and (c).
Experiment 1 mean stresses and propagation geometries for two different values of
. The contour interval is 50 kPa. Dashed lines represent additional observed geometry. Squares indicate propagation increments. (a) Initial mean stresses considering presence of fracture boundaries. (b) With
the long fracture propagates once to relieve shear stress, producing a kink at the fracture tip. (c) The fracture propagates with
. Step 9 in the calculation is shown. (d) MOA image of the fractures investigated here and in experiment 2.
Using
induces the observed propagation of the artificially shortened test fracture 1β for K
IC = 0.3 Mpa m1/2 in this domain (Fig. 6c). The shorter 1α is in the stress shadow of the longer 1β and does not propagate. The transverse tip of test fracture 1β propagates in the same direction as the observed fracture but, unlike the observed fracture, the model fracture propagates through the Evans/ Carlson suture zone. Thus, the suture zone appears to be important to the observed pattern of fracture tip arrest at this structural boundary.
Experiment 2
The upstream-pointing segments of the fractures propagating from suture zone D become less distinct with distance downstream, presumably as the relict features fill with snow. Western tips of the open (transverse) fracture segments coincide with suture zone D. Transverse segments begin as short features, only a few kilometers in length, and grow eastward to a length of 10–12 km, arresting at the Evans/ Carlson suture zone (Figs 2 and 4). While the features originated in this suture zone far upstream, the zone appears to be a barrier to westward propagation of still-active transverse segments of the fractures.
Two questions arise from the observed fracture evolution. First, does suture zone D in fact limit westward fracture propagation? Second, is limiting westward propagation (as we did in experiment 1) important to eastward propagation of these fractures? We investigate these questions using a 5 km test fracture, 2α, that retains the geometry of experiment 1 but without the upstream-pointing segment (Fig. 7). Eliminating the upstream-pointing segment aligns the western side of the fracture with principal compressive stresses, a situation that should be favorable for propagation. The region boundary is modified so the upstream boundary is closer to the single transverse fracture. Calibration follows the usual procedure. The best calibration is obtained with
and elasticity constants v = 0.28 and E = 8000 MPa (2a in Table 1). Boundary stresses are scaled along the fracture boundary to initiate propagation and both tips are active.
Experiment 2 mean stresses and propagation geometries. (a) Initial mean stresses for calibrated elasticity parameters (2a in Table 1). Test fracture 2α is displayed to orient the reader. Observed values are shown in Figure 5a. The contour interval is 50 kPa. (b) Step 6 in the calculation, eastward propagation of the right-lateral tip is limited. (c) Step 18 in the calculation, left-lateral tip propagates westward, an unobserved phenomenon. Propagation of the right-lateral tip is promoted by westward lengthening of the fracture. (d) An experiment with a shorter initial fracture length and the same remote stresses (2b in Table 1) required K IC = 0.1 MPa m−1/2 for propagation
Sustained propagation of the western tip of the fracture boundary is generated with K
IC = 0.3 MPa m1/2 and
, while propagation of the eastern tip is episodic. The eastern tip experiences one episode of eastward propagation when the fracture length reaches ∼8 km (due to westward propagation) and a second episode when the fracture length reaches ∼10 km (again due to westward propagation). Eastward propagation follows the observed propagation direction.
Stresses at the western end of the fracture boundary are large enough to yield propagation through suture zone D, contrary to observation. This implies that the suture zone inhibits propagation, supporting the hypothesis under consideration. Unlike experiment 1, in which the eastern tip of the propagating test fracture extended to the edge of the model region, the test fracture in this experiment does not propagate eastward without significant (unobserved) westward lengthening. From this, we conclude that eastward propagation requires inactive western tips, again supporting the hypothesis that inhomogeneities in the ice are important to fracture propagation in the ice shelf. Additional experiments with longer fracture geometries downstream of this location, not reported here, yield similar results (Reference LeDouxLeDoux, 2007).
The calculation is repeated with a shorter, 2.6 km, test fracture boundary, 2β, and a slightly narrower region boundary (Fig. 7d; 2b in Table 1). Propagation initiates at the western tip of the test fracture with K
IC = 0.1 Mpa m1/2 and =
. As the calculation proceeds, the western tip propagates ∼1 km beyond suture zone D before initiating a single episode of propagation at the eastern tip to release accumulated shear stress along the fracture. The western tip then continues to propagate as a mode I fracture.
Experiment 3
Westward propagation of fractures previously arrested at suture zone D eventually produces >100 km long fractures near the shelf front (Fig. 2). As the fractures investigated in experiments 1 and 2 approach the front, right-lateral tips propagate a short distance eastward to the Evans/Carlson suture zone, where they arrest. Some fractures then propagate westward, through several EIS tributary suture zones.
The region boundary is placed in the area of two ∼70 km long subparallel fractures spaced ∼15 km apart (Figs 4 and 8). The upstream fracture is used to define our test fracture boundary and to evaluate the results of our calculations. The test fracture is ∼95 km from the shelf front, and the downstream end of the region boundary is ∼80 km from the front. Principal stresses in this region are relatively uniform, except near the western shear margin and very close to the front (Fig. 4). Following an initial experiment, we translate the fracture boundary upstream in order to find a best fit between modeled and observed fracture geometries.
Experiment 3a mean stresses and propagation geometries. The contour interval is 50 kPa. (a) Observed and (b) calibrated mean stresses. Test fracture 3α is shown to orient the reader. (c) Initial mean stresses considering presence of the fracture boundary. (d) Mean stresses and fracture geometry at step 36 in the calculation. The right-lateral tip is at the region boundary. Were the domain to be expanded laterally, the left-lateral tip of this very long fracture would continue propagating with similarly poor representation of the observed fracture geometry.
Initialization is challenging for the large region boundary required to embrace long fractures. In particular, care is required with regard to effects of shear in the west and correct representation of the small stress gradients in the eastern two-thirds of the domain (Fig. 8a and b). The best calibration is obtained with region boundary stress scaling
and elasticity constants v = 0:27 and E = 7500 MPa (3 in Table 1). Note that the magnitude of the resulting stress field is on average ∼100 kPa, or ∼7% lower than observed mean stresses.
We first consider propagation of a single fracture originating in ice from the easternmost Evans tributary. The initial test fracture boundary, 3α, follows the easternmost 20 km of an observed 70 km fracture with its eastern tip at the Evans/ Carlson suture. When both tips are active the eastern tip propagates eastward to the edge of the model domain (Fig. 8d). This behavior is not observed and we repeated the experiment with the eastern tip fixed as non-propagating (not shown). With K
IC = 0.3 Mpa m1/2,
, the test fracture propagates westward, as observed, past suture zones between ice from several EIS tributaries, but the modeled geometry deviates significantly from the observed geometry. This suggests that the fracture did not propagate in this stress field, rather in a stress field somewhere upstream.
The effect of advection through a non-uniform glaciological stress field can be treated in an ad hoc way by translating the test fracture boundary upstream and varying its initial length, now referred to as test fracture boundary 3β. This procedure does indeed improve modeled fracture geometry (Fig. 9). Other combinations of upstream translation and modification of the initial boundary length could produce similar improvements in the modeled fracture geometry. A comparison of mean stress magnitudes near both fracture tips (Fig. 9b and c) provides additional quantitative support that suture zones are important. Here the left-lateral tip of a long fracture propagates into an unfavorable stress field near the shelf margin, and the fracture arrests when the right-lateral tip is unable to propagate. As the left-lateral tip approaches a region of lower glaciological stresses, stresses near the right-lateral tip become increasingly concentrated (due to longer fracture length and greater glaciological stresses). To test the importance of fracture length in this region, we shortened test fracture boundary 3β and translated it only 5.7 km upstream from test fracture 3α. The test fracture did not propagate.
Experiment 3b mean stresses and propagation geometries, using test fracture 3β translated 8.9 km upstream from 3α and lengthened. The glaciological and initialized mean stresses are shown in Figure 8. The contour interval is 50 kPa. (a) Initial mean stresses considering the presence of the fracture boundary. (b) Step 12 in the calculation. (c) Step 42 in the calculation. (d) Result of an additional experiment with fracture boundary translated only 5.7 km upstream from 3α and shortened from 3β; no propagation occurred.
Experiment 4
The westward-propagating tips of the long fractures examined in experiment 3 pass relict fractures in ice between the A and B sutures (Fig. 4). These short fractures do not propagate through adjacent shear margins and disappear from view in the MOA downstream of this region (Fig. 2). We return to the original experiment 3 set-up to examine interaction between the long propagating fractures and the shorter, apparently relict, features. Here the model is used to investigate the roles of fracture geometry and tip interaction due to modification of near-field stresses. None of our calculations with the additional fracture boundary reproduce observed geometries.
A long fracture boundary, 4α, that traces the observed fracture in experiment 3a is used and a new fracture boundary, 4β, with a shallow V shape is added to the domain (Fig. 10). The 4α boundary is longer than in experiment 3, a change that favors correct reproduction of the observed geometry (as discussed in experiment 3). The shape of fracture boundary 4β represents the intersection of fractures that propagated eastward from suture zone A and westward from suture zone B. Here we show the outcomes of calculations with two different 4β geometries. The right-lateral tip of 4α is held fixed while all others are free to propagate.
Experiment 3b mean stresses and propagation geometries, using test fracture 3β translated 8.9 km upstream from 3α and lengthened. The glaciological and initialized mean stresses are shown in Figure 8. The contour interval is 50kPa. (a) Initial mean stresses considering the presence of the fracture boundary. (b) Step 12 in the calculation. (c) Step 42 in the calculation. (d) Result of an additional experiment with fracture boundary translated only 5.7 km upstream from 3α and shortened from 3β; no propagation occurred.
As the calculation begins, the left-lateral end of the 4β fracture boundary propagates over one iteration to relieve shear stress. Propagation of 4β begins in earnest as the left-lateral end of the 4α boundary approaches 4β and its right-lateral tip is affected by the associated stress concentration in the region around the propagating 4α tip. As 4β propagates (Fig. 10c), it modifies principal stress directions in the near field and a characteristic shape for interacting, propagating fracture tips develops (Reference Olson and PollardOlson and Pollard, 1989; Reference Cruikshank, Zhao and JohnsonCruikshank and others, 1991).
In a second calculation, a kink is added to the right-lateral end of the 4β boundary. The result is similar to the previous calculation, in that the short fracture does not propagate until the left-lateral tip of the long fracture approaches. The final geometry of 4β is, however, very different as a result of the initial kink at the right-lateral end of the boundary (Fig. 10d). In both cases, the final geometry of the long fracture is similar to that obtained in experiment 3 because the near-field modification due to the short fracture is minor: the right-lateral tip of the short fracture is distant. Propagation of the 4β fractures is not observed in the ice shelf. Again, we conclude that the suture zones bounding these features inhibit crack tip propagation.
Discussion
A small subset of the many transverse fractures present in the ice shelf become planes along which tabular icebergs calve. The experiments presented here address why this should be the case, but the underlying problem of tabular iceberg production remains. Detailed documentation of fracture propagation preceding a major calving event is limited (Reference Lazzara, Jezek, Scambos, MacAyeal and van der VeenLazzara and others, 1999). One reasonably well-observed event is the calving of iceberg B-15 from the eastern front of the Ross Ice Shelf in March 2000. In that event, large fractures, like the near-front fractures examined here, connected with fractures generated by shear at the corner of the ice-shelf front to form the 300 km long and 40 km wide iceberg. Both the pre-existing fractures and new fractures forming at the front corner of the ice shelf were necessary, but neither alone was sufficient to generate the massive calving event. Only near the front do fractures grow to the lengths required to produce large icebergs. We take this as support of our emphasis on horizontal propagation, but the importance of vertical propagation to the fractures studied here can be evaluated by calculating stress intensity at a vertical tip in the near-front environment.
Following Reference Van der VeenVan der Veen (1998) and Reference Scambos, Hulbe, Fahnestock and BohlanderScambos and others (2000), the stress intensity for mode I opening in a dry fracture is the sum of the effect of tensile stress and the effect of the overburden, which tends to close the fracture. The first component
is in which f(λ) is an empirically derived function of the ratio of crevasse depth, d, to ice thickness and τyy is the deviatoric stress normal to the fracture plane. The lithostatic component
is in which ρ i and ρ s represent the ice and surface snow densities, respectively, C is a densification constant taken to be 0.02 m−1 and G is a numerically derived function of γ = b/d and λ. Using the same parameter values as Reference Scambos, Hulbe, Fahnestock and BohlanderScambos and others (2000), along with ice thickness and stresses for the fracture in experiment 3, and a fracture depth of 50 m (an estimate following Reference Fricker, Bassis, Minster and MacAyealFricker and others, 2005), the resulting stress intensity factor at vertical tips at this depth is negative, indicating that vertical propagation is not favorable in this region. While vertical propagation (either from the surface or the base) is clearly required at the moment of calving, it is not what sets in place the fundamental plane-view geometry that must be followed when large icebergs calve.
Conclusion
Visual inspection of high-resolution satellite imagery reveals that surface crevasses in Antarctic ice shelves are abundant, but that they remain relatively short over most of their residence time and appear to arrest at structural boundaries in the ice (Reference LeDouxLeDoux, 2007; Reference GlasserGlasser and others, 2009). One explanation for the limited length is inhibition of propagation due to stress shadows cast by adjacent fractures. We demonstrate this effect in experiment 1. However, some fractures must be long for others to be inhibited, so an explanation is still required for the relatively short lengths of most fractures. We demonstrate other characteristic interaction patterns in experiment 4. A second possibility is that the glaciological stress fields are simply too small to promote propagation. We address this by shortening our test fractures relative to observed fracture lengths and reproducing observed geometries. Once propagation is initiated, it tends to continue because stress intensity increases with fracture length. On occasion, for shorter transverse fractures in regions of larger shear stresses (Figs 7d and 9d), glaciological stress fields relative to the fracture geometry may simply be too small. A third possibility is that inhomogeneity in ice properties plays a fundamental role in governing fracture length. We demonstrate this effect in experiments 1, 2a and 3, in which modeled glaciological stresses would have been sufficient to drive propagation.
Suture zones between tributaries to the Ronne Ice Shelf are clearly associated with the arrest of propagating fracture tips. This notion is supported anecdotally by visual inspection of satellite imagery. It is supported quantitatively by our numerical model experiments. In each experiment, we pose the question: would a propagating fracture arrest in its observed location in the absence of inhomogeneity in ice properties? The answer is consistently no. Because the models use far-field glaciological stresses and homogeneous material properties (i.e. variations due to the presence of a suture zone are not treated) we conclude that suture zones play an important role in fracture mechanics within the ice shelf, in particular by promoting crack tip arrest.
Acknowledgements
The work was supported by NSF OPP awards No. 0125754 and No. 0440670. The authors thank I. Joughin for the velocity data. Comments from two reviewers helped clarify the text.
Appendix
We use a boundary element method (BEM) to solve the elastic boundary value problem. BEMs require discretization of only the domain boundary and fractures, resulting in a smaller set of equations than would arise from other approaches (e.g. a finite-element method). In particular, we use the displacement discontinuity method, a BEM developed to handle slit-like openings or thin fractures (Reference CrouchCrouch, 1976; Reference Crouch and StarfieldCrouch and Starfield, 1983). For the present work, we have modified the Frac2D software developed by A.L. Thomas (unpublished information, 1991), following Reference Crouch and StarfieldCrouch and Starfield (1983). We provide an overview of the method here. A more complete treatment is given by Reference LeDouxLeDoux (2007).
A fracture can be viewed as a set of point disturbances within a stress field. The remote, or glaciological, stress field is modified in the near field by these disturbances. Our goal is to investigate that modification in a finite region within what we assume to be an infinite plane (horizontal ice-shelf dimensions are large compared with fracture dimensions). In the numerical scheme, the set of disturbances are midpoints of straight-line boundary elements along a fracture trace. Analytical expressions are written to describe the effect of each point disturbance on the stress field within the finite region. Linearity of the governing equations allows the complete solution to be formed by summation of individual solutions.
The model boundaries are a contour defining the finite region containing a fracture or set of fractures and the traces of fractures themselves (Fig. 11). Boundary conditions consist of tractions (shear and normal) and displacements on the elements. It is convenient to define a local coordinate system,
,
, for each element, rotated an angle β from a fixed external system x, y, with its origin at the element midpoint. Shear and normal stresses, σ
s and σ
n, hereafter referred to as ‘boundary stresses’, are defined at the midpoints of the elements:
The mean stress is taken to be the remote stress associated with viscous deformation of the ice. Removing the remote stress leaves residual, near-field boundary stresses, b s and b n. Stress state and material properties are assumed to be constant along each boundary element.
Schematic of boundary element setup, modified from Reference Crouch and StarfieldCrouch and Starfield (1983). (b) is the fracture boundary shown in (a). Local reference frames
are used for each boundary element.
Each boundary element along a fracture trace can be thought of as two coincident surfaces. A displacement discontinuity is the movement of one surface relative to the other. The displacement is assumed to be constant along the length of a given boundary element. Displacement discontinuities,
and
, are
in which
and
represent components of the displacement in the local reference frame and the superscripts indicate positive and negative sides of the fracture. Displacement discontinuities are assumed continuous everywhere except across the fracture, and change in displacement is assumed to be constant for any given boundary element. By convention, negative
(normal to the element) indicates that the surfaces are moving away from each other (opening mode) and negative
(shear) indicates that the positive fracture surface is moving to the left relative to the negative fracture surface (Fig. 11).
The solutions for constant displacement along a line element in an infinite plane (Reference CrouchCrouch, 1976; Reference Crouch and StarfieldCrouch and Starfield, 1983) are displacements
and stresses
The function
is
for an element with half-length a.
Boundary stresses are simulated by evaluating the cumulative unknown shear and normal displacement discontinuities due to all boundary elements (fractures and the region boundary) on all other boundary elements. Each boundary element, j, exerts a stress influence on each other boundary element, i. The length, orientation and location of an element j relative to element i are used to compute influence coefficients:
Each coefficient represents the shear or normal stress acting on the midpoint of the ith element due to a constant unit shear or normal displacement at the midpoint of the jth element.
The relationship between boundary stresses, bn and bs, influence coefficients and displacement discontinuities is
Using the principle of superposition, the cumulative influence of the shear and normal displacement discontinuities of all j elements on the ith element is
for a domain of N elements. The result is a set of 2N simultaneous linear equations for the 2N unknowns for each element i. In practice, fracture geometries are not simple and each set of equations requires a series of coordinate transformations.
The model is used to simulate fracture propagation in a stress field that evolves with the growing fracture. Starting from an initial stress field and fracture geometry, we iterate on the system of governing equations until a steady state is found between the model stress field and fractures, that is, stress intensity at fracture tips falls and remains below that required for propagation.
