Introduction
The evolution of ice shelves is controlled by many factors (Reference Jacobs, Hellmer, Doake, Jenkins and FrolichJacobs and others, 1992), one of which is the episodic calving of large tabular icebergs. Although much progress has been made in iceberg monitoring, the mechanisms behind calving are still poorly documented and understood (Reference ReehReeh, 1968; Reference HughesHughes, 1998, p. 196). Some of the factors that are not well known include the timing of calving events, the origin of the large rifts that are precursors to a calving event, their mode and rate of propagation, the factors controlling their propagation and how this chain of mechanisms could be affected by climate change.
The objective of this work is to study the evolution of rifts prior to a large calving event. Our area of interest is Hemmen Ice Rise (HIR) on the Ronne Ice Shelf, Antarctica. HIR is surrounded by an extensive field of rifts that are at the origin of large calving events (Reference Hartl, Thiel, Wu, Doake and SieversHartl and others, 1994; Reference Rignot and MacAyealRignot and MacAyeal, 1998). On 13 October 1998, a tabular iceberg 145 × 50 km2 in size broke off from the ice shelf. The iceberg detached along a pre-existing line of weakness, namely a rift emanating from HIR. This event provides a unique opportunity to study the mechanics of rift propagation because extensive satellite imagery was collected before (and after) the event.
The data consist of synthetic aperture radar (SAR) images spanning a 6 year time period prior to the major calving event. The radar images are processed interferometrically to obtain both the horizontal and vertical velocities of the ice shelf around and near the rifts. Time series of radar amplitude imagery are also used to measure the propagation of the rifts and study their time evolution.
The discussion presents a series of results. First, we study the influence of the tidal motion of the Ronne Ice Shelf on the propagation of rifts. Second, we examine the contribution to their propagation from the horizontal flow of the ice shelf. Third, we compare the results with a viscous model of ice and with linear elastic fracture mechanics (LEFM). Finally, we discuss how these results will help to develop a better modelling of fracture processes and iceberg calving on an ice shelf.
Study Area
Hemmen Ice Rise, on the Ronne Ice Shelf, lies on the eastern flank of Berkner Island (Figs 1 and 2). This small ice rise (20 km long by 3 km wide) is at the origin of a vast field of rifts, some up to 40 km long and up to 5 km wide. These rifts originate along the margins of HIR. Numerous crevasses form in the shear margins of the ice rise. Surface crevasses are clearly visible in the radar imagery, but there are probably also numerous bottom crevasses (Reference Van der VeenVan der Veen, 1998b). Eventually, some of these crevasses develop into rifts, i.e. deep incisions through the entire ice-shelf thickness, which create major lines of weakness in the ice shelf. The rifts then propagate westward, parallel to the ice-shelf front, at rates of about 1000 m a-1, which is comparable to the ice-shelf velocity. On 13 October 1998, the rift marked 3 in Figure 1a became unstable and rapidly propagated through the ice shelf, giving birth to tabular iceberg A38. This large calving event was soon followed by a series of other calving events in the same region (http://www.ccrs.nrcan.gc.ca/ccrs/rd/apps/marine/ice/calv98_e.html).
ERS-1 geocoded SAR images of the study area around HIR near Berkner Island, on the Ronne Ice Shelf, (a) ERS-1 amplitude radar image obtained in February 1992. (b) Corresponding interferogram. The fringe rate has been chosen so that each fringe corresponds to ΔV = 10 m a-1, where V is calculated according to Equation (1). The direction of north, the satellite track direction and the range direction are indicated with arrows. Small arrows near the rifts show the direction of the velocity vector of the ice-shelf flow. The propagation rates of rifts are given with double arrows. Profile A-B in the right panel is plotted in Figure 4a. ©European Space Agency 1992.
Sequence of RADARSAT-1 amplitude SAR images during the calving of iceberg A38 on 13 October 1998.
The complex pattern of ice motion around HIR prior to calving is discussed in detail by Reference Rignot and MacAyealRignot and MacAyeal (1998). The ice shelf flows around HIR and along Berkner Island towards the east. The presence of large rifts north of HIR releases the stresses in the ice shelf, modifying the ice flow. The ice-shelf trajectory deviates to the northeast between HIR and the ice front. Yet the presence of melange between the rifts (Reference MacAyeal, Rignot and HulbeMacAyeal and others, 1998) tends to maintain cohesion to the rifts. In Figure 1a, the ice melange between rifts 1 and 2 is connected/attached to the ice rise, and therefore constrains its motion, whereas ice between rifts 3, 4 and 5 is free to flow northwards. This flow to the north contributes to the opening of the rifts and to building of stresses at the rupture tip of rifts, ultimately resulting in a propagation of these rifts into the ice shelf. Similarly, this flow may help to transform crevasses along the shear margins of HIR into new rifts. The appearance of these rifts follows a remarkable periodicity of about 10 km, or 10 years, between rifts.
Observations
The data consist of interferograms derived from European Remote-sensing Satellite-1 (ERS-1) and RADARSAT-1 radar images between 1992 and 1997. These interferograms measure the deformation of ice accumulated over, respectively, 9 and 24 days. Combining these interferograms in different ways, it is possible to separate the continuous process of creep deformation from the more cyclic motion due to ocean tides. It is also possible to evaluate changes in the geometry of the rifts and their rates of propagation, through the analysis of radar amplitude images. Table 1 summarizes the images used and the combinations employed.
SAR images used to study the period preceding calving of iceberg A38
Geometric observations
The propagation of the active rifts is measured from SAR amplitude images acquired in 1992, 1996, 1997 and 1998. The radar amplitude images are geocoded at a sample spacing of 50 m on a polar stereographic grid. We measure the propagation rates by feature tracking along the rifts. The results are given inTable 2.
Propagation of rifts, 1992-98
Rifts 4 and 5 are inactive: no significant propagation or shortening of the rifts is observed in the data. This result is confirmed by the fact that imprints of rifts 4 and 5 are still visible on iceberg A38 (Fig. 2b). The inactivity of these rifts may be due to the thick melange filling them (Reference Rignot and MacAyealRignot and MacAyeal, 1998). Towards the ice front, the melange is older and thicker. This layer transmits stresses from one side of the rift to the other, thus acting as a coherent material (Reference MacAyeal, Rignot and HulbeMacAyeal and others, 1998). Furthermore, these are the rifts most distant from HIR, and hence least influenced by lateral shear from HIR. These conditions make it harder for the rifts to propagate.
Rifts 1-3 are actively propagating from 1992 to 1998, but the time series of measurements is not dense enough to determine whether the propagation is continuous or proceeds through growth/arrest increments. The latter process is more likely to occur, as suggested by Reference Parsons, Snellen, Muggeridge, Cocks and PonterParsons and others (1989) and Reference DeFranco and DempseyDeFranco and Dempsey (1994), although their conclusions concerned samples of sea ice in a stable crack- growth configuration, which may not apply to glacier ice. Nevertheless, Table 2 shows that rifts 2 and 3 have a similar pattern of propagation, with most propagation occurring between 1992 and 1996. Rift 1, unlike rifts 2 and 3, did not propagate between 1992 and 1996. The propagation took place mainly between February 1996 and October 1997. This is consistent with the position of rift 1 along the HIR margins: it is located (Figs 1a and 2a) at a point where the bay sides diverge. Reference SandersonSanderson (1979) has shown that the state of strain rate reaches a maximum in the ice shelf at this point, which facilitates rapid propagation of rifts.
Dynamics of the ice shelf
Interferograms are used to characterize the ice-shelf dynamics. Different combinations of interferograms are used to obtain both the horizontal creep flow of the ice shelf and the vertical deviation due to tidal perturbations.
Small perturbations to the ice flow
Interferograms are built from pairs of images separated by 24 days for RADARSAT-1 (1997) and 9 days for ERS-1 (1992). Figure 1b shows the 1992 interferogram and the corresponding geocoded amplitude image. Radar phase differences are converted into a velocity, V, using (Reference RignotRignot, 1996):
where ϕij is the flattened interferometric phase between images i and j (with the influence of topography and interferometric baseline removed), taken at time t i and tj , Zi and Z j are the vertical positions of the ice shelf (positive upwards) at time t i and tj , λ is the wavelength of the radar (λ = 5.6 cm), ψ is the angle between the local vertical of the ice shelf and the radar illumination direction and Vx and Vz are the horizontal and vertical velocity of the ice shelf, with the x axis parallel to the ground and perpendicular to the satellite track. V is displayed in Figure 1b.
We introduce a first-order development of Vx and Vz in Equation (1), to precisely evaluate small perturbations to V, which are of interest here.
where 
is the creep deformation of the ice shelf, and (δV
X
, δ
V
y
, δV
z
) a small perturbation to it. Equation (4) is taken from Reference Joughin, Kwok and FahnestockJoughin and others (1998) and takes into account the influence of the slope tan(α) on the velocity vector, assuming that ice flows parallel to the ice surface.
Horizontal displacement of the ice shelf
We will first study the evolution of the horizontal creep flow 
which is the dominant term in V (Equation (5)) Figure 1b shows V around HIR in 1992. Figure 3a and b show a close-up of V around the rupture tip of rift 3 in 1992 and 1997, respectively. Figure 4a shows a profile of V across the five main rifts (profile A-B from Fig. 1b) in 1992, and Figure 4b shows a profile along rift 3 of the velocity differential ΔV between opposite flanks of the rift in 1992 and 1997.
Interferograms for rift 3 in 1992 (a) and 1997 (b).The fringe rate for V is the same in 1992 and 1997 (2 m a-1 per cycle).
(a) Profile A-B (Fig 1b) of Vacross rifts 1-5 in 1992. Each shift in velocity, ΔV, is marked with the corresponding rift number. (b) Evolution of the velocity differential ΔV along rift 3 in 1992 and 1997.
Between 1992 and 1997, the direction of the radar differs by 9.9° modulo 180° (because ERS-1 and RADARSAT-1 were viewing in opposite directions); therefore the 1992 and 1997 interferograms measure nearly the same projection of the velocity vector, so that their comparison reveals changes in flow velocity. Note that the direction of the x axis is almost perpendicular to rift 3, which facilitates the evaluation of opening rates.
The main feature that can be observed in Figures 1b and 3 is the offset in fringe rate across rifts. There is a difference in projected velocity between the two flanks of the rift. The corresponding velocity difference is a measure of the opening rate of the rift, which is a good indicator of the propagation of the rupture tips. This velocity differential is observed on rifts 1-3 (Fig. 4a), which are actively propagating (Table 2), but not on rifts 4 and 5.
Figure 4b shows the velocity differential ΔV along rift 3 in 1992 and 1997. The velocity differential increases linearly from the rupture tip, until it reaches a threshold where it becomes independent of the distance, d tip, to the rupture tip. This shows two distinct zones: one of pure rotation (20 km long in 1992, 30 km long in 1997) and a subsequent zone where the flanks of the rift are moving away from each other at a velocity of 55 m a-1 in 1992 and 62 m a-1 in 1997. This implies that the active size of the rift to be taken into account in any fracture theory is smaller than the actual size of the rift because the part where ΔV is independent of d tip is not influencing the opening rate of the rift and is not participating in the state of stress around the rupture tip.
Tidal oscillations
As shown by Reference RignotRignot (1996), the tidal signal can be isolated by applying a double-difference technique with three or more radar images. This technique is applied here with three images acquired in February 1992. Calling the three images, separated by 3 days, 1, 2 and 3, the double differencing of the phases leads to:
This double difference only depends on the tidal displacement and the incidence angle if we assume that the creep motion of the ice shelf is steady and continuous, which is a good approximation here. Figure 5 shows the corresponding interferogram. Each fringe represents an incremental change in elevation of 3.4 cm in 3 days. On a 10 km scale, the average tidal elevation differential ranges from 35 cm. For instance, the ice constrained between rifts 2 and 3 is 3 cm higher than the ice constrained between rifts 3 and 4.
ERS-1 double-difference interferogram in 1992. Each fringe represents a tidal elevation increment of 3.4 cm.
Such a tidal oscillation could trigger a mode III fracture (Reference AndersonAnderson, 1995). Across the rifts, the tidal offset is at most half a fringe, or 1.7 cm, but on either side of them the phase differential varies very slowly. The offset between each side of the rift diminishes toward the tip and no singularity is observed at the tip. Locally the contribution of the tidal signal to V in Equation (5) is negligible. On a scale of 10 km and more, the tidal elevation differential is about 85 m a-1, which represents one-tenth of 
and therefore cannot be ignored in the interpretation of simple difference interferograms. No differential interferometry was available in 1997, but the study of a similar double difference in 1996 (Reference Rignot, Padman, MacAyeal and SchmeltzRignot and others, 2000) reached the same conclusions. For the rest of the study, the influence of the tidal signal along the rifts is therefore ignored.
Vertical deformation of the ice shelf
We are now discussing the possibility of vertical motion of the ice shelf other than tidal oscillation. As can be seen in Equation (5), 
contribute to V to the first order. This combination makes it difficult to isolate and evaluate the vertical motion δVz
.
The contribution of the slope to V can be calculated. From our digital elevation model (Reference Bamber and BindschadlerBamber and Bindschadler, 1997), we have tan(α) = 2 × 10-4 in absolute value. With in the order of 1000 m a-1, the resulting contribution 
in absolute value.
ERS-1 is a right-looking SAR and RADARSAT-1 was left-looking during the first Antarctic Mapping Mission. This allows us to separate -δVx
from δV
z
cot(ψ) and 
in Equation (5). As Table 1 shows, the ERS-1 interferogram (3069-2940) uses the first image (orbit 3069) as reference; the RADARSAT-1 interferogram (985210195) uses the second (orbit 9852). The time-spans of these two interferograms are therefore opposite in sign. If remains the same in 1992 and 1997, and ignoring the tidal oscillation contribution along the rifts, we have in 1992:
and in 1997:
where the x axis corresponds to the across-track direction in 1992. Note that ψ remains constant to the first order between 1992 and 1997 (ψ = 23.4° in 1992 and ψ = 28° in 1997).
is unchanged between Equations (7) and (8) because of the opposite-sign time difference combined with the opposite-sign direction of view. In contrast, δV
z
cot(ψ) and change sign between 1992 and 1997 because the z axis remains unchanged but the time-spans have opposite sign in the interferograms.
Figure 3 shows distinctive patterns along the flanks and at the tip of the rifts. These patterns are better detected in the derivative of the signal taken in the direction of view. In Figure 6, the background horizontal creep is homogeneous and the residual deformation pattern around the rifts is more visible.
Derivative of the velocity (strain rate) taken in the direction of view of the radar in 1992 (a) and 1997 (b). Aclose-up of rifts 5 and 3 in 1997 is shown in Figure 7.
Note in Figure 3 that the patterns of deformation near the rupture tip are opposite in sign between 1992 and 1997. In 1992, the fringes curve towards the tip, whereas in 1997 the fringes curve outwards. Equations (7) and (8) show that this type of reversal is characteristic of either a slope-induced horizontal motion or a vertical deformation, but could not be due to the background creep of the ice . Here, it can only be a vertical motion. The motion observed could not be due to the snow accumulation on the surface, because the patterns are localized along the rifts and at the rupture tips.
Figure 7a and c (left parts) show V
for rifts 5 and 3 in 1997. The horizontal contribution has been removed by subtracting a linear fit on each side of the rift. In each case, bands of upwards positive velocity are present on each side of the rifts. Profiles taken along each rift (Fig. 7b and d) show that these bands are 1 km wide, with peak values of 0.8 m a-1. At a distance of 500 m from the rifts, decreases strongly to become negative.
Close-up of rifts 5 (a, b) and 3 (c, d) in 1997. (a) Vertical velocity maps observed in 1997 (left) and deduced from the linear viscous bending model (right). Profiles 1 and 2 are plotted in (b) together with the corresponding model profiles. The same layout is plotted for rift 3 in (c) and (d). Rift 3 has been artificially rotated with relation to Figure 6.
Towards the tip of the rift, a negative singularity occurs only for active rifts 1-3. The tip patterns are similar in this case. In the case of inactive rifts 4 and 5, the bands disappear towards the tip, leaving no observable singularity. Similar bands are observed at the ice front. This suggests that this deformation could be associated with a bending of the ice shelf, similar to the one observed on ice fronts prior to the calving of icebergs. We cannot rule out a local variation of the slope along the rifts and around the rupture tips, but this local variation would have to be opposite on each side of a rift, which is unlikely. We therefore choose to model these perturbations as vertical motion induced by a bending of the ice shelf.
Viscous Bending of the Ice Shelf
To explain the vertical deformation of the ice shelf shown in Figure 7, a model developed by Reference ReehReeh (1968) is generalized to the two-dimensional case. The deformation of the plate is controlled by the imbalance between hydrostatic pressure in the ice and hydrostatic pressure in the water on the faces of the rift, which generates a vertical motion of the ice plate. Because of the time-scales involved (years), the analysis is made using a linear viscous regime. For a complete explanation of the concepts, we refer the reader to Reference ReehReeh (1968). The model and the finite-element implementation are presented in detail in the Appendix.
Figure 8a shows the evolution of Vz along a profile perpendicular to the rift, according to the model. Perpendicular to the rift, VZ reaches a maximum V zmax about 500 m from the rift. The maximum decreases in time, and moves towards the rift. V zmax is responsible for the narrow bands observed along the rifts. Figure 8b shows the evolution of this maximum in time for different viscosities μ. V zmax decreases asymptotically with time, and the rate of decrease is inversely proportional to μ.
One-dimensional model for the linear viscous bending of an ice shelf (a) Modelled vertical velocity profile perpendicular to a rift. Viscosity μ is 1014 Pa s. (b) Evolution of modelled V zmax vs time for different viscosities. The observed vertical velocity peaks are also plotted. The different profiles are used to evaluate the best-fit viscosity of the ice shelf.
The generalization of the model to a plate allows for the interpretation of the observed rupture tip patterns. Both a static (no propagation) and dynamic (active rift) modelling of the rupture have been carried out. In the dynamic case, the mesh is artificially propagated at the rupture tip at each time increment. The geometric configuration of the plate used is shown in Figure 9.
Geometric configuration used in the finite-element model of linear viscous bending and in the LEFM double-cantilevered beam propagation model. L = 30 km, a = 10 km, e = 6 km and h = 500 m. For the finite-element model, elevation and slope are constrained to 0 far from the rift (red boundaries) and a bending moment is applied on the flanks of the rift (green boundaries). For clarity, the mesh displayed has only 6000 elements. The actual mesh used in the computations has 60 000 elements.
Static case
The static case is used to adjust the parameters of the model to achieve the best fit with the observations. Several calculations have been conducted with different widths, e. A good fit is obtained for e > 3000 m. To choose the viscosity, the velocity maps were used to compute strain rates and evaluate the real viscosity using:
where σ′ is the stress deviator, εe the effective strain rate, 
the strain rate and B the flow constant. μ is computed from (Reference MacAyeal, Rignot and HulbeMacAyeal and others, 1998):
We use 
The viscosities obtained from the observations are in the range 1013-1015 Pa s, with the largest values recorded far from the rifts, and the smallest near the rifts. If the estimated viscosities are taken as parameters in our model, the results remain the same as those obtained with a constant viscosity μ = 1014 Pa s.
Another way of evaluating the best-fit viscosity is to calculate the evolution of V zmax in time for different viscosities (Fig. 8b).We compare this result with the observed values of V zmax for different rifts. The age of each rift is estimated from the distance to HIR and the mean horizontal ice-shelf velocity. The best fit is obtained for μ = ¼ 1014 Pa s. Note that rift 4 deforms more rapidly. We will return to this later.
From now on, we carry out the calculations with μ = 1014 Pa s. The model results, together with observations for rift 5, are shown in Figure 7a. Figure 7b shows a comparison between the model and observations on two profiles selected perpendicular to the rift. The following parameters are used: L = 30 km, 2e = 12 km, ρwater = 1023 kgm-3, ρ ice = 900kg m-3 and h = 500m. L is the length of the plate, 2e the width, ρwater the water density, ρice the ice density and h the ice thickness (Fig. 9). The ice density takes into account the effect of low-density firn layers in the surface of the ice shelf. The computation spans a time period of 40 years.
Figure 7a and b show a good fit between observation and model. There is no tip deformation pattern observed, which confirms that rift 5 is not propagating. This result is consistent with the observed propagation rates in Table 2. Near the edges of the rift, the model predicts a strong decrease in V z . This cannot be verified in the observations because the signal is not reliable in this area (because of multiple reflections in the rift). Nevertheless, the observed V z decreases near the edges, which is consistent with the model.
Dynamic case
In the dynamic model, the crack is artificially propagated. At each time-step, one cell of the mesh in front of the rupture tip of the rift is split in two. A static analysis of the linear viscous bending is done, iterating on the step before. Therefore, our model does not take into account any dynamic effects during the propagation of the rift. Consequently, our model is compatible with: (1) a propagation continuous in time; and (2) a stick-slip type of propagation with jump- arrest increments of the rift (Reference Parsons, Snellen, Muggeridge, Cocks and PonterParsons and others, 1989), as long as each of the jump-arrest increments is inferior in length to the rupture tip element characteristic size.
In our model, the value of μ and the size of the rupture tip cell determine the average in time of the propagation rate of the rift. The results are shown in Figure 7c and d, together with the observations for rift 3 in 1997. The propagation rate used in the model is 1000 m a-1.
Rift 3 exhibits the same VZ positive bands along the rift, but the deformation at the tip is different. Compared with the observations in Figure 7a and c, a wide deformation pattern develops at the tip. Positive bands of vertical deformation surround the rupture tip, and a strong negative singularity develops at the tip. The model explains these conditions, which suggests that this is an active rift.
Some differences between the model and the observations appear in profile 2 of Figure 7d. The magnitude of VZ is larger in the model near the tip of the rift. This difference could be due to the assumptions made about μ: μ is independent of temperature and stress fields, which could vary substantially towards the rupture tip. Although it is in principle possible to deduce the propagation rate of rifts from the magnitude of tip deformation, the process is complex and we have not pursued this aspect.
Elastic Fracture Singularity
The vertical deformation pattern along rifts has been interpreted as a vertical viscous bending of the ice shelf. We now evaluate the concentration of stress predicted by the LEFM. This mode I rupture yields a horizontal deformation as the rift propagates.
LEFM predicts the displacement differential across a rift in mode I (Reference HellanHellan, 1984) in the polar (r, θ) coordinate system to be:
where Δu is the differential displacement (Fig. 9), K I is the stress intensity factor and G is the shear modulus of ice (G = E/2(1 + v) with v the Poisson’s ratio).
The rift is propagating at a rate v = ∂r/∂t; therefore:
E is the Young’s modulus. We use K I = 0.2 MPa m-1/2 (Reference Mulmule and DempseyMulmule and Dempsey, 2000), v = 1000 ma-1 and E = 9 × 109 Pa (Reference Petrenko and WhitworthPetrenko and Whitworth, 1999, p. 39) to obtain Δv = 5.610-10 r -1/2.
In order to observe Δv with interferometric synthetic aperture radar (InSAR), Δv needs to exceed 0.1m a-1, which means r ≤ 4 cm. This is much less than the resolution cell (7 m) of the radar, hence Δv is not observable in the interferograms. Even it it were detectable, the signal would be masked by the viscous deformation of ice (Δv viscous ≫ Δv elastic).
Nevertheless, as shown by Reference Mulmule and DempseyMulmule and Dempsey (2000), the behaviour of large cracks in tabular ice plates of size ≥200 m should be well explained by LEFM. Although ice around HIR is mainly in creep regime, rifts propagate at rates superior or comparable to the background ice flow (Table 2). In this case, LEFM is expected to yield a good prediction of the propagation rate of rifts. This theory was applied successfully to study crevasse penetration (Reference Van der VeenVan der Veen, 1998a, Reference Van der Veenb).
From 1994 to 1998, the propagation of rifts 1-3 appears to be stable. The configuration of rifts 2 and 3 is similar to a double-cantilevered beam. This configuration is known to be stable under displacement-controlled propagation (Reference Dempsey, Bentley and SodhiDempsey and others, 1986; Reference Parsons, Snellen, Muggeridge, Cocks and PonterParsons and others, 1989). We adopt here a similar modelling for the rifts. Figure 9 shows the geometry used.
The displacement differential Δu between the two flanks is (Efunda Engineering Fundamentals, http://www/efunda.com/formulae/solid_mechanics/beams/casestudy_bc_cantilever.cfm):
where Δu is the distance between the two flanks of a rift, at distance x from the tip, I is the moment of inertia of the beam (I = he 3/12), e the distance between two rifts and a the length of the rift. According to LEFM, the action of external forces on a rift can be reduced to the action of a stress σ on the flanks of the rift (Reference AndersonAnderson, 1995, p. 66). This stress is due to a distribution of water and lithostatic pressure on the flanks of a rift, and a distribution of longitudinal stress due to the ice flow. The assumption is made here that σ is constant, because of the relatively slow variation in the background creep of the ice shelf (Fig. 6). The strain energy due to this repartition of stress on the flanks of the rift is:
The driving force ϕ is derived from (Reference AndersonAnderson, 1995, p. 43):
which yields:
σ can be eliminated using Equation (12) at x = a:
where α = 4EI/h .
If we derive Equation (15) in time and replace a by ϕa 4/Δu(a)2, we obtain:
where v = da/dt
If the propagation is assumed to be continuous in time, the driving stress ϕ is equal to the resistance R(a) of the ice shelf at each time increment (Reference Dempsey, Bentley and SodhiDempsey and others, 1986):
R(a) is unknown in this study. A common assumption is that R is constant and independent of a, as in the case of a brittle fracture, in which case:
This formulation assumes that the propagation is a continuous process. Extensive laboratory experiments on ice have shown, however, that under displacement control a stick- slip propagation occurs. To take this type of propagation into account, we define, following Reference Parsons, Snellen, Muggeridge, Cocks and PonterParsons and others (1989), an initiation fracture toughness Kc
and a crack arrest intensity factor K
a. In mode I fracture, toughness is linked to resistance (Reference HellanHellan, 1984, p. 58) by 
and 
where we define R
c and R
a as the initiation fracture resistance and crack arrest resistance of the ice shelf, respectively. R
a is a material characteristic, whereas R
c depends on the loading rate. Defining α = R
c/R
a = (K
c/K
a)2, and considering as before a brittle fracture (dR
c/da = 0, dR
a/da = 0), we obtain:
Tests carried out by Reference Parsons, Snellen, Muggeridge, Cocks and PonterParsons and others (1989) on sea ice yield α = 1.54, which gives c = 1.78. Introducing a stick- slip mechanism of propagation yields a faster propagation rate. The case α = 1 gives c = 2, which is the continuous case where K c = K a.
Equations (19) and (20) can be used to evaluate the propagation rate from observations. The quantity 
is evaluated at the beginning of the rift using the horizontal velocity maps obtained from the interferograms. The results are shown in Figure 10. The observed positions of active rifts 1-3 between 1992 and 1998 are compared with positions determined using propagation rates evaluated from Equation (19). The calculated propagation rates overestimate the observations in 1992 but are broadly correct in 1997. The misfit in 1992 could be removed by using α ≤ 1. This would mean that the crack initiation toughness was inferior to the crack arrest toughness, which is not realistic. Another explanation could be that the configuration of the rifts in 1992 differs from that in 1997. In 1992, rifts 1-3 were less far open than in 1997: the double-cantilevered beam configuration may not be suitable for the modelling of rift propagation in 1992.
Evolution of the propagation of rifts compared with a LEFMmodel, based on a double-cantilevered beam configuration. (a) corresponds to rift1, (b) to rift 2 and (c) to rift 3.Two different propagation rates have been evaluated in1992 and1997 using Equation (16), except for rift 1, which did not have the required double-cantilever beam configuration in 1992.We plotted the observed propagation distance of the rifts as observed in the SAR amplitude images, and overlaid our predicted curves.
Discussion and Conclusions
This study shows that very subtle patterns of deformation are revealed by InSAR along rupture tips of rifts in an ice shelf. In the case of the Ronne Ice Shelf, a fair degree of precision was reached in the observation of tidal oscillations. It was possible to conclude that elastic fracture in mode III is unlikely. Yet the fatigue associated with long-term tidal oscillations (22 000 cycles in 30 years) could play an important role in the propagation of rifts, which we did not investigate further.
The observed vertical deformation was explained using a linear viscous plate bending model. This viscous model has simplifying assumptions that we discuss here in light of the model/data comparison:
-
1. The ice shelf is assumed to behave like a plate. At the tip and near the edges of rifts this hypothesis is valid: the vertical deformations computed were <10% of the ice- shelf thickness;
-
2. The dependence of viscosity on temperature variations and stress distribution throughout the thickness was not modelled. Thus, the model will not predict any variation of the bending through the ice-shelf thickness. These variations could be large because of the strong gradient of temperature throughout the ice shelf;
-
3. The thickness of the ice shelf was assumed to be constant, which is a good approximation except possibly along the ice front and along the rifts, because of the presence of bottom melting. High bottom melt rates along rifts could modify the thickness distribution and therefore alter the viscous bending process. According to the linear viscous bending model, an 80% linear decrease in thickness starting 1 km from the rift would decrease V Zmax by 55% at any given time. Measuring the ice-shelf thickness along the rifts is critical to improve the fit of the model;
-
4. The influence of calving was ignored: no limit was put on the bending angle reached by the edges of rifts. This situation is not realistic. As discussed by Reference ReehReeh (1968), the surface longitudinal stress produced by bending of an ice shelf leads to fracture of the ice and release of icebergs. The anomaly spotted in Figure 8b for rift 4 could be due to a calving event. The presence of debris ice in the rifts, as discussed by Reference Rignot and MacAyealRignot and MacAyeal (1998), confirms this hypothesis. Taking calving into account would increase the predicted vertical velocity, unless the viscosity used in the model is reduced;
-
5. The boundary conditions along the flanks of the rifts were determined assuming that the sections remained vertical during deformation. This assumption may not hold at the tip of a rift in the dynamic case. The presence of a vertical deformation singularity implies high bending angles, which lower the bending moment on the flanks of rifts. Our deformations may therefore have been overestimated.
Overall, the viscous model gives good qualitative results, but the simplifying assumptions discussed above preclude inversion of propagation rates from tip deformation patterns. The differences between model and observations in Figure 7d (profile 2) show that linking V Zmax at the tip of the rift to the propagation rate is difficult. Yet the state of propagation of the rift (active or inactive) is well indicated by the deformation pattern recorded at the tip of the rift.
The linear viscous bending model implies the use of continuum mechanics in a quasi-static approach to rift propagation. Here, the dynamic propagation of rifts was modelled by an artificial splitting of elements at the rupture tip of the rift. This procedure allows our linear viscous bending model to be valid in case of continuous propagation of the rift, or in case of stick-slip propagation, as long as each slip is smaller than the rupture tip element
Similarly, in our LEFM model of rupture tip propagation, a stick-slip propagation approach is used that converges toward a continuous propagation mode if the correction factor c tends towards 2. From Equations (18) and (20), propagation rates depend on the evolution of the resistance of the ice shelf. The instability of rifts will result from an evolution of the resistance R(a). Typically, resistance to crack growth increases as crack length increases (Reference Dempsey, Bentley and SodhiDempsey and others 1986), which means that the double-cantilevered beam configuration is stable. The onset of final rupture therefore stems from either a modification of the loading configuration or atypical evolution of the resistance. A study of the geometric configuration prior to the rupture indicates that the double-cantilevered beam configuration is valid before the calving event. Further studies should be conducted to assess the evolution of the resistance of the ice shelf R(a) in order to correctly evaluate its importance in triggering the final rupture.
Acknowledgements
This work was performed at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration, Cryospheric Sciences Program. We thank the European Space Agency, the VECTRA project and Alaska SAR Facility for distributing the radar data employed in this study. We thank the reviewers and especially N. Reeh for their valuable input into our work and their excellent recommendations, and the editor, H. Rott, for his help in the reviewing process. Finally, we thank M. Schmeltz for her assistance with image processing and finite-element modelling.
Appendix: Vertical Linearviscous Bending Model
Only the modifications used to generalize the Reeh one-dimensional viscous beam model to a two-dimensional linear viscous plate model are presented here. The geometry used for the creation of the model is shown in Figure 9. The dimensions of the plate model are chosen in order to correspond to the configuration of rift 3.
Reference ReehReeh (1968) used the Kirchhoff hypothesis (Reference HuangHuang, 1988) for thin beams. We used the Mindlin-Reissner hypothesis for thick plates to develop our model. This was done to avoid continuity problems between finite elements in our formulation.
We start from the well-known equilibrium equations of a floating plate:
where Mx , My , Mxy are the moments over a section of the ice shelf, Z is the elevation of the ice shelf, Qx , Qy are the vertical shear forces, ρw is the density of the water and g the acceleration of gravity.
A linear viscous model is chosen because the characteristic times of the phenomena observed are of the order of a year. Elevation Z, moments Mx, My, Mxy and vertical shear forces Qx, Qy are linked according to Reference NadaiNadai (1963):
where μ and ρi are the viscosity and the density of ice, h is the thickness of the ice shelf and θx , θy are respectively the x and y rotation of a vertical section of the ice shelf.
To solve the problem using a finite-element formulation, we use the following boundary conditions far from the rift:
A rift is modelled by two flanks separated by a very small distance. On the flanks, the repartition of water and litho- static pressure yields a bending moment (Reference ReehReeh, 1968):
with 
the normal and tangent unitary vectors to the rift faces, 
and d
i = ρi/ρw the ice density.
On each flank, the shear stress due to the water friction is nil:
with 
Equations (A1-A8), together with the boundary conditions, are solved using finite-element methods. The discretization in time was done using a Euler implicit scheme that is always stable. The mesh was refined (60 000 elements) to take into account the singularity of the stress field around the tip of the rift (Fig. 9). The validity of our model was checked against Reeh’s computations. The vertical deformation on a profile chosen perpendicular to the rift was identical to the deformation at a calving ice front (Fig. 8).
For the dynamic case, we artificially split a mesh element at the tip of the rift to simulate its propagation. Different propagation rates were considered and in each case the results were stationary after 20-30 iterations.
