Introduction
Of the two main processes that remove mass from ice shelves (basal melt and iceberg calving), iceberg calving is the less understood. Basal melt rates have been estimated for the major Antarctic ice shelves (Reference Rignot and JacobsRignot and Jacobs, 2002; Reference Rignot and ThomasRignot and Thomas, 2002; Reference Joughin and PadmanJoughin and Padman, 2003) and directly measured and/or modeled on ice shelves such as the Amery Ice Shelf (AIS) (Reference Hemer, Hunter and ColemanHemer and others, 2006). Iceberg calving, however, is of a sporadic nature, making it more difficult to quantify. Large iceberg-calving events tend to occur as part of a natural cycle where the front advances (by ice flow) beyond its embayment walls and then retreats (by iceberg calving). Recurrence intervals between large calving events vary between ice shelves, but typically are on the order of several decades (Reference BuddBudd, 1966; Reference Scambos, Hulbe and FahnestockScambos and others, 2003). While such events are necessary to maintain the mass balance of the Antarctic ice sheet, the concern is that they may become more frequent in response to atmospheric and oceanic warming. This concern is particularly valid because of the potential of iceberg calving to remove large amounts of ice very rapidly. Therefore, small perturbations in calving frequency can translate into large changes in the mass balance of the ice shelves and thus significantly accelerate trends in deglaciation.
Icebergs detach when one or several fractures near the calving front isolate an iceberg. For ice shelves, large fractures that penetrate the entire ice-shelf thickness (termed ‘rifts’) tend to originate along sharp pointed segments of the ice margin, pinning points or suture zones between ice streams (Reference Lazzara, Jezek, Scambos and MacAyealLazzara and others, 1999; Reference Fricker, Young, Allison and ColemanFricker and others, 2002). There is currently limited knowledge of the processes driving rift propagation, and numerical ice-sheet models either do not incorporate calving at all, or treat iceberg calving using untested ‘calving laws’ (Reference Benn, Charles and MottramBenn and others, 2007). Without proper treatment of the calving process in models, any current quantitative prediction of the response of an ice sheet to climate change is seriously compromised.
This study aims to remedy our limited understanding of the rifting process through a detailed field study of a propagating rift on the AIS. This study extends the previous study by Reference Bassis, Coleman and FrickerBassis and others (2005) in which a network of eight seismometers and six GPS (global positioning system) receivers were deployed around the tip of one of the rifts studied by Reference Fricker, Young, Coleman and BassisFricker and others (2007) on the AIS during the 2002/03 austral summer field season. The pilot study revealed that rift propagation is episodic, with recurrence intervals of 10–24 days. Each burst, deduced from seismic swarms (detected with seismometers) coincident with rapid rift widening (detected with GPS), lasted ∼1 to 4 hours. Between bursts there was a relatively constant background rate of seismicity. The source of the detected seismic events clustered around the rift tip, providing compelling evidence that the seismicity was related to rift propagation. However, with only one field season of measurements and only three propagation events, it was difficult to know whether these observations were representative of typical behavior and how varied the recurrence intervals between propagation events are, a detail that has important implications for the process by which stress accumulates and is released at the rift tip (Reference Shimazaki and NakataShimazaki and Nakata, 1980). In addition, the small number of seismic stations and low sampling frequency made it difficult to accurately locate the source of the seismic events. This made it impossible to determine if seismicity around the rift tip is diffuse, evidence of a large area (or process zone) in which strain energy is dissipated, or highly localized, as is predicted by some ‘thin-strip’ process-zone models. For example, the Dugdale–Barenblatt model predicts a thin zone of cohesion ahead of the rift tip (Reference BarenblattBarenblatt, 1959; Reference DugdaleDugdale, 1960).
In this paper, we extend the results of Reference Bassis, Coleman and FrickerBassis and others (2005) with data from two additional field seasons, during which we deployed a more densely distributed network of seismometers and GPS receivers. The network was specifically designed to improve our ability to identify and locate rift-related icequakes. These measurements were conducted on the AIS over the austral summers 2004/05 and 2005/06. Although one of the goals of this study was to confirm the results from 2002/03, we also wanted to extend the time series of rupture events to provide a better statistical sample of (1) recurrence times between swarms and (2) the number of events in each swarm and magnitude of rift widening during each swarm. We also aimed to map the spatial extent and pattern of fracture events. Although study of the spatio-temporal pattern of seismicity is a common technique in earthquake physics, it has never before been applied to ice-shelf rifting.
Location of study
Our study site is located near the front of the AIS (Fig. 1). Historic shipboard sightings and imagery dating back to the Corona mission show that the last major calving event from the AIS occurred in late 1963 or early 1964 when a 10000 km2 iceberg detached, removing close to one-fifth of the total surface area of the shelf (Reference BuddBudd, 1966; Reference Fricker, Young, Allison and ColemanFricker and others, 2002). Since then, the AIS has steadily advanced and is expected to return to its most advanced position in the mid-2020s, when another major calving event is predicted (Reference Fricker, Young, Allison and ColemanFricker and others, 2002). Satellite imagery has revealed that in the process of ice-shelf re-advancement, several rifts have formed near the front of the AIS.
Upper right: Map showing the location of the Amery Ice Shelf in East Antarctica. Left: MODIS image acquired on 18 January 2006 showing the Amery Ice Shelf and the ‘Loose Tooth’ rift system. The rift is currently propagating into a suture zone formed where two ice streams merged (black curve). Lower right: LANDSAT image acquired on 18 December 2002 showing a close up of the Loose Tooth rift system. The L1–T1–T2 triple junction was first observed in 1995. Since then L1 has widened but has not increased in length.
Our efforts have focused on monitoring a system of these rifts that form the outline of an intermediate-sized iceberg (30 km by 30 km) termed the ‘Loose Tooth’. The associated active rift system (Fig. 1, bottom right panel) consists of two longitudinal-to-flow rifts ∼30 km apart that initiated ∼20 years ago (L1 and L2) and two transverse-to-flow rifts (T1 to the west and T2 to the east) that initiated at the tip of L1, forming a triple junction first observed in 1995. Both rifts T2 and T1 initiated in the transition zone where transverse-to-flow strain rates begin to exceed longitudinal-to-flow strain rates (Reference Young and HylandYoung and Hyland, 2002).
Rift T2 currently propagates at about 4 md−1 (Reference Fricker, Young, Coleman and BassisFricker and others, 2007), a rate that is approximately equivalent to the flow speed of ice in this region (Reference Bassis, Coleman and FrickerBassis and others, 2005). When T2 connects with L2, an iceberg 900 km2 in area will calve. L2 has not propagated at all since monitoring began in 1996 but has widened by several kilometers. Analysis of historical ice-front positions shows a Loose-Tooth-sized indent in the ice front that preceded the major calving event of 1963/64 (Reference Fricker, Young, Allison and ColemanFricker and others, 2002). This suggests that the Loose Tooth rift system may repeatedly form as part of the cycle of advance and retreat of the ice shelf (Reference Fricker, Young, Allison and ColemanFricker and others, 2002). It also suggests that the detachment of the Loose Tooth may be a precursor to another major calving episode. For these reasons our field campaign was concentrated near the tip of rift T2, whose propagation will ultimately determine the timing of the release of the Loose Tooth iceberg.
In the field we observed that the rift tip was not well defined, but instead tapered off into a series of micro- and mesoscale cracks (Fig. 2a–c). During a helicopter reconnaissance flight, we identified a position near the end of this network of fractures, where no cracks were visible on the surface of the ice shelf, as the location of the rift tip. This location was typically several kilometers ahead of where the rift tip was identified in satellite images, due to errors in geolocation of the image and limitations due to the image pixel size (∼250 m). During the helicopter flights we also noticed that the interior of the rift was filled with large blocks of ice that had fallen in from the sides (shown in Fig. 2b and e). Some of these blocks were almost completely covered with snow. Between field seasons we saw substantial variation in the amount and thickness of the debris that filled the rift. We were not able to tell whether this was because of active rearrangement of blocks within the rift or seasonal variation in the snow that covered the rift. We also noticed that the rift propagates through a field of longitudinal-to-flow (and normal-to-rift) crevasses (indicated with arrows on Fig. 2). Early in the field seasons snow bridges obscured many of these. Later on, many of the snow bridges appeared to sag and collapse, revealing the full extent of the crevassing. Another feature that we noticed was that the northern/seaward side of the rift was uplifted ∼1 to 2 m higher than the southern side, with the amount of uplift decreasing towards the tip (Fig. 2d).
(a) Photograph taken from a helicopter looking west towards the triple junction Illustrating the field of transverse crevasses. (b) Photograph taken looking east towards the rift tip showing how the rift tapers off towards the tip. The transverse crevasses are also clearly evident in the photograph as are blocks of ice that appear to have fallen into the rift from the rift walls. (c) Photograph taken near the rift tip looking east showing the width of the rift. (d) Photograph looking north at the rift ∼2 km from the visible rift tip showing the uplift of the northern (ocean) side. (e) Photograph looking east showing the mixture of snow and ice blocks that fills the rift.
Survey design
Our seismic network was designed to detect high-frequency icequakes generated by ice-fracturing events near the tip of rift T2. For this reason, we installed observing stations in a dense network around the rift tip. As the experiment progressed over different field seasons, our observing equipment was modified slightly, but, in general, our observation stations each operated one seismometer and one GPS receiver. During the 2002/03 field season we had six stations, each with a single vertical component L-4C seismometer recording at 10 Hz (0.1 s) and a dual-frequency GPS receiver recording at 0.033 Hz (30 s), plus two additional seismometer-only stations (Reference Bassis, Coleman and FrickerBassis and others, 2005, fig. 1). In the 2004/05 and 2005/06 seasons we increased the number of stations to 12, each with a three-component L-28 seismometer digitized with a Quanterra Q330 data logger and a dual-frequency GPS receiver recording at 0.5 Hz (2 s), which recorded for 52 and 81 days, respectively. Hereafter we refer to the 2002/03 field season as ‘season 1’, the 2004/ 05 field season as ‘season 2’ and the 2005/06 field season as ‘season 3’. To facilitate comparison of timing between field seasons, we also reference days of our survey relative to day of year (DOY) 332.
A network of stations that encloses the rift tip enables more accurate location of seismic events. Therefore, to maximize the chance that the rift tip would be enclosed by our network, we placed instruments both far ahead of and behind where the surface expression of the rift tip was visible. We also deployed a line of stations normal to the rift at the tip (stations a–f in Fig. 3). These stations enable us to determine how rapidly the rift is widening and compare rift widening between station pairs on opposite sides of the rift, and, as a way of detecting differences between the seaward side of the rift and the landward side, we had stations on the same side of the rift with comparable baselines. Other stations were placed at intermediate distances to form triangles along the length of the rift from which we could compute horizontal strain rates. The geometry of the network (Fig. 3) was maintained between seasons 2 and 3, but the center of the network was translated to account for both ice flow and rift propagation so the rift tip remained close to the center of the network at the beginning of each field season.
Diagram showing the relative positions of observing stations (gray triangles) relative to the tip of rift T2, indicated by a gold star for seasons 2 and 3. Observing stations were placed relative to where the tip of rift T2 was observed in the field.
Most of the L-28 seismometers operated continuously during deployment. During season 2 the aluminum poles on which the GPS antenna were placed melted out and tipped over during the first 20 days because they had not been placed deep enough in the snow. This problem was mitigated for season 3 by placing poles deeper into the snow and switching from aluminum to wooden poles. Unfortunately, we still observed some artifacts in the GPS signal, possibly related to motion of the GPS poles or to frost lenses forming on the GPS antennae after day 30 during season 3. Because we cannot unambiguously separate real glaciological signals from those caused by the motion of the poles, we exclude the GPS data of season 2 from this analysis and also GPS data after day 30 of season 3.
The similarity in network geometry facilitates comparisons between seismicity detected during seasons 2 and 3. However, the distances between stations were not exactly the same nor did we reoccupy the exact same geographical sites. There is an additional uncertainty because the distance between the ‘true’ rift tip position and the center of our network may differ, depending on how accurately the tip could be determined in the field. The low sampling frequency of our seismometers (due to storage limitations) of 10 Hz during season 1, combined with the anti-aliasing filter used for analysis, cut our bandwidth to frequencies of 1–5 Hz. Unfortunately this barely overlapped with the response of the L-28 seismometers deployed during seasons 2 and 3, with a natural period of 5 Hz. This, in combination with different network geometry, makes a quantitative comparison with season 1 more challenging.
Seismic Processing and Results
Passive seismology is a powerful tool to study ice-shelf rifting; not only can seismic measurements be used to determine the spatial and temporal distribution of ice-fracturing events, but also the elastic waves excited by these events contain information about the rupture process itself. Exploiting this information requires some processing. The general processing strategy we use is: (1) generate a catalog of all seismic arrival times at each station; (2) associate arrivals at different stations with a common event; and (3) use the differential travel times between stations to locate the source of the seismic event. Once a catalog of events has been created, the locations of the events can be used to exclude seismicity that is not related to the rift, resulting in a population of events related to ice-shelf rifting. The spatio-temporal distribution of this population of events contains information about the ice-shelf rifting process. The details of our implementation of the algorithms to process the seismic data are given in the Appendices. Briefly, the basis of our detection algorithm is the STA/LTA algorithm described by Reference WithersWithers and others (1998) and we use a grid-location algorithm with a constant velocity to locate events. In our seismic analysis, we consider only seismic events detected at four or more stations. This eliminates spurious associations caused by single-station noise. It also enables us to locate all the events used in our analysis and to determine whether they have an origin within our study region or outside it (e.g. small calving events from the ice front).
Temporal distribution of seismicity
We detected over 10 000 individual seismic arrivals during seasons 2 and 3, with a larger number of events detected at the two stations closest to the rift tip (c and d) than at stations further away from the rift tip (e.g. stations a and f; see Fig. 3). We restricted our analysis to those events with arrivals visible at four or more stations – the subset of events that we could locate. This reduced the dataset to 4160 events during season 2 and 5330 events during season 3, corresponding to a mean rate of seismicity of 80 events per day during season 2 and 63 events per day during season 3, a net decrease in the overall seismicity of roughly 10% and marginally above the statistical uncertainty in our detection measurements.
To examine fluctuations in seismicity over shorter time periods, we computed the number of events in 2 hour time bins (upper two panels in Fig. 4). During season 2 there was a relatively constant background seismicity punctuated by three large swarms of seismicity on days 30, 53 and 68. Each swarm lasted 1–3 hours and consisted of 120–175 events. In contrast, there is little clear swarm-like activity during season 3. We detect one clear swarm (on day 16) along with a smaller bump that might be a smaller swarm around day 52. Both of these swarms contain fewer events (100 and 40, respectively) than the swarms observed during season 2. In particular, the second swarm during season 3 is much smaller than any of the three swarms observed during season 2.
Temporal distribution of seismicity. Upper two panels: histograms of the number of icequakes detected at four or more stations (bin size = 2 hours) for seasons 2 and 3. Lower two panels: temporal distribution of the relative magnitudes of seismic events. Seismic swarms are emphasized with shaded bars. There are a few smaller spikes of seismicity, hinting at the presence of smaller swarms that are not as obvious as the main swarms within the seismic record.
We next looked at temporal variations in the relative magnitude of events, determined from the log10 of the peak amplitude of the waveform with the largest amplitude (neglecting attenuation and geometric spreading). The variation in relative magnitudes with time is shown in the lower two panels of Figure 4. Swarms appear to consist of a mixture of events of all sizes. However, large events did not occur exclusively during swarms; there were a few large-magnitude events that occurred between swarms. For example, there are a few large events before and after the first swarm during season 2. Some of these large-magnitude events appear to correspond to smaller spikes in the seismicity. These smaller spikes may be evidence of the presence of smaller seismic swarms that are less well demarcated within the data (denoted in Fig. 4 with arrows). We also found that the main swarms did not consist of a main ‘shock’ of much larger magnitude than all of the other events followed by ‘tails’ of smaller aftershocks or preceded by foreshocks, as is often seen in earthquake studies. There is, however, some evidence of accelerating seismicity leading up to some of the swarms. The swarm on day 30 during season 2, in particular, seems to be accompanied by an increase in the magnitude of events leading up to the swarm and then a decrease in the magnitude of events following the swarm. This suggests that while ice-shelf rifting shares some of the features of earthquake rupture, there are also significant differences between the two processes.
Spatial distribution of seismicity
The large number of events we detected at four or more stations enabled us to produce the first maps of the spatial distribution of rift-related seismicity (Fig. 5). The relative magnitude of events is denoted by the size of each circle (determined, as before, from the log10 of the peak amplitude of the waveform with the largest amplitude). The seismic velocity that minimized the sum of the squares of the residuals at all stations for all events was 2250 ms−1, similar to both the shear wave speed in meteoric ice and the velocity in the upper 50 m of the ice shelf (Reference McMahon and LackieMcMahon and Lackie, 2006). This indicates that either (1) seismic events have shallow hypocenters or (2) we detect little P-wave energy from the waveforms, a result that seems consistent with the results of Reference Stuart, Murray, Brisbourne and StylesStuart and others (2005).
Map showing events (blue circles) that were located with respect to the network geometry for season 2 (left panel) and season 3 (right panel). The size of each circle is proportional to the log10 of the peak amplitude of the seismogram at the station closest to the event. Gray triangles indicate the locations of the seismometers. Dashed lines indicate regions where we may be observing propagation of normal-to-rift crevasses.
We initially expected that seismicity would be concentrated in a small annulus or disk ahead of the rift tip. Instead, the pattern of seismicity for both seasons 2 and 3 indicates that there is a highly clustered line of seismicity along the length of the rift axis. Surprisingly, this cluster of seismicity has a long tail of events extending several kilometers back from the rift tip towards the triple junction. During season 3, there is a large cluster of seismic events located ∼500m ahead of the rift and offset from the main rift-centric cluster of seismicity. This cluster of seismicity is accompanied by an increase in diffuse seismicity ahead of the rift tip. An additional feature present in both seasons is swaths of seismicity oriented approximately in a normal-to-rift direction. We strongly suspect that these swaths are related to the location of some of the normal-to-rift crevasses that we observed in the field.
To test whether this spatial pattern of seismicity was caused by the geometry of our network, we cross-validated our results by relocating all events, each time with a different station omitted, for all 12 stations; this is akin to the so-called ‘jack-knife’ method in statistical estimation. This did not significantly affect the horizontal positions of the events. However, it did alter many of the source depths we derived by up to several hundreds of meters. Unfortunately, this suggests that our network does not have adequate resolution to determine depths of events within the ice shelf. For this reason, we restrict our analysis to horizontal positions. We did not try to identify focal mechanisms for seismic events because it appears that the first detected arrival may be a shear wave.
Distribution of event magnitudes
One of the most ubiquitous features of earthquake seismicity is the power-law size–frequency distribution of event magnitudes. The frequency of earthquakes with magnitudes not less than M can generally be fitted according to the Gutenberg–Richter relationship
where A and b are constants (b is typically ∼1), M is the earthquake magnitude and N is the number of events per day with magnitude not less than M. Equation (1) implies that there are many more small events than large events. Previously, we only considered events detected at four or more stations. This puts a lower limit on the size of events we can detect that makes it difficult to accurately assess the size–frequency distribution. Instead, we calculated the size–frequency distribution for each station for both seasons 2 and 3. Three representative examples are shown in Figure 6. As we are unable to locate events detected at fewer than four stations, we cannot correct for geometric spreading and attenuation and thus may overestimate the true size–frequency distribution. For all stations (and both field seasons) we find power-law behavior over three orders of magnitude. For each station the value of b is ∼1, similar to earthquakes, ranging between 0.8 and 1.2, and reasonably consistent between field seasons. Despite neglecting geometric spreading and attenuation we can say that, at least qualitatively, like earthquakes, there are considerably fewer large events than small events and there does not appear to be a characteristic size of rupture events.
Size–frequency distribution for icequakes determined from three different stations for season 2, where M is the relative magnitudes of events and 
is the number of events per day with magnitude not less than M. All three stations show the same power-law behavior with slopes (b) of ∼1.
GPS Processing and Results
GPS processing
GPS data were processed using the Geodatic Inc. Real-Time Dynamics (RTD) software package with instantaneous network positioning (Reference Bock, Nikolaidis, de Jonge and BevisBock and others, 2000), previously applied to monitoring earthquakes (Reference Nikolaidis and BockNikolaidis and others, 2001) and volcanoes (Reference Mattia, Rossi, Guglielmino and AloisiMattia and others, 2004). Because of the small extent of the network, we used independent L1 and L2 phase GPS observations, effectively neglecting differential atmospheric effects (Reference Bock, Nikolaidis, de Jonge and BevisBock and others, 2000). We used a satellite elevation cut-off of 10°. Since the entire network is in motion, we allowed all station coordinates to freely adjust. An advantage of kinematic processing is that common motion due to ice flow (lateral) and tides (vertical) is removed. This means that the time series of GPS positions are determined relative to one of the sites within the network. Using this method we obtained relative positions between the GPS receivers at every epoch (every 2 s). We then chose to express the positions of the sites relative to site c (Fig. 3), generating a time series of positions relative to site c for each of the other sites. For seasons 2 and 3, the higher sampling rate and shorter baseline distances yielded solutions with less noise than for season 1 (2 s cf. 30 s). After processing, we first removed a linear trend and then discarded all outliers greater than 3.5 times the interquartile range of the detrended baseline lengths. The censored time series of positions was then smoothed using a 30 min Gaussian running mean to remove high-frequency variability in the GPS processing.
Horizontal strain rates
To determine strain rates we tessellated our network of GPS sites into individual triangles using Delaunay triangulation (e.g. Reference Schroeder and ShephardSchroeder and Shephard, 1988). A typical example, shown in Figure 7, shows the displacement between station pairs is dominated by a linear trend, indicating that the background strain rate within the ice is relatively constant. We then computed relative velocities between stations using a least-squares fit to the time series of positions determined from our GPS receivers. The velocity gradients, calculated from the slope of the fit between station pairs, were then used to compute strain rates for each triangle individually, using standard geodetic techniques (see, e.g., Reference MalvernMalvern, 1969). The eigenvalues and eigenvectors of the matrix representation of the strain-rate matrix give the principal strain rates and principal axes (shown in Table 1). Standard errors computed for the strain rates were then computed for each triangle (see, e.g., Reference TaylorTaylor, 1997, ch. 8). The errors shown in Table 1 are quite small, typically several orders of magnitude smaller than the strain rates. The standard error is almost certainly an underestimate of the true error, as the strain rate may vary across each triangle and each GPS-derived position measurement is not independent of every other position measurement because errors due to the ionospheric or atmospheric correction may be correlated. However, it is clear that the errors in the strain rate are relatively small compared to the estimated strain rates.
Upper panel: displacement of site e relative to site c for season 3. Lower panel: residual of the displacement of site e relative to site c after a linear trend has been removed. Transient signals, which may be related to motion of the poles on which the GPS antennae were mounted, are visible after about day 30 and are highlighted with a gray box. This signal is incoherent across the network, suggesting that it is an artifact rather than a real glacio-logical signal.
Principal strain rates calculated for each triangle. E 11 is the argest principal strain rate and E 22 is the smallest. θ is the angle that the largest principal axis makes with the east-west axis. The formal error in the calculation is probably an underestimate of the true error. The largest uncertainty in the direction of the principal axes was 2°
All but two of the strain rates were extensional with typical magnitudes ranging between ∼10–5 and 10–4 d–1, about an order of magnitude larger than strain rates near the front of the AIS away from the rift tip (Reference Young and HylandYoung and Hyland, 2002). Triangles e–k–f and d–e–k show a small compressional strain rate in a direction approximately normal to the rift axis which may be related to widening of the rift. The principal axes of the strain rates are displayed graphically in Figure 8. The largest principal strain rates occur in triangles anchored on both sides of the rift (triangles d–k–c and g–d–c). In these triangles, the principal strain rate is aligned in a direction close to rift-normal. In contrast, those triangles that outline an area that is entirely on one side of the rift have principal axes that are close to rift-parallel, suggesting that normal-to-rift stress is primarily accommodated by rift widening. Overall strain rates determined from triangles that cross the rift are about an order of magnitude larger than those from triangles situated entirely on one side. This indicates that the rift is widening at a rate that is larger than the glaciological spreading rate and may explain the small compression evident in triangles e–k–f and d–e–k. In addition, there is a rotation in the orientation of the principal axes of the strain rates counterclockwise in the stations ahead of the rift tip (see, e.g., triangles j–i–a, g–j–i, g–d–e), a feature consistent with numerical models of the strain-rate field near the tip of rift T2 (personal communication from R. Warner and J. Court, 2007), but at odds with the symmetric strain fields often observed at crack tips (Reference LawnLawn, 1993). The rotation in the principal strain rates near the tip of rift T2 is probably related to the rigid-body rotation of the Loose Tooth as rift T2 lengthens and rift L1 opens. The rigid-body rotation may act as a vorticity source in the velocity field.
Map of the orientation of the principal axes of the strain-rate tensor for season 3 (black lines). The length of each axis is proportional to the magnitude of the strain rate. The triangles outlined in gray show the station triangulation used to calculate each strain-rate tensor. The approximate rift tip (shown as a black square) was located between stations c and d. An approximate outline of the rift (filled gray region) is included to show the orientation of each triangle used to compute a strain rate relative to the rift.
Transient strains
The baseline extensions for each station pair are dominated by linear trends due to the constant background spreading of the ice, the magnitude and orientation of which is approximately given by the strain rates discussed in the previous subsection. However, the lower panel of Figure 7 shows that, in addition to the large-scale flow of ice, there is substantial temporal variation in the strain rate. We found that the variations in the strain rate after day 30 became incoherent across the network and thus may be due to processing artifacts introduced by, for example, the ‘settling’ of the GPS poles after temperatures began to rise or frost forming lenses on the GPS antenna. Because we cannot be sure if this signal is of real glaciological origin, we discarded all data after day 30 from the rest of our analysis.
To examine transient signals in baseline extension prior to day 30, we computed rift-normal and rift-parallel displacement for each GPS epoch (i.e. every 2 s interval). We then removed the large-scale motion by subtracting a linear fit to the long-term trend, resulting in a time series of detrended epoch-by-epoch displacements for each station pair, as was done in the lower panel of Figure 7. Figure 9 shows 10 days surrounding the first swarm in season 3. (Recall that we have rejected GPS data from season 2 because the poles on which the GPS antennae were mounted fell over.) All three baselines that contain stations on both sides of the rift (pairs c–d, c–e, c–g) show a sharp jump of about 1 cm in the width of the rift, coincident with the swarm of seismicity. The coincidence of the widening with the seismic swarm and coherence of the signal across several station pairs indicates that this signal, unlike that after day 30, is real and not a processing artifact. The larger noise in some of the displacements is caused by the longer baselines between station pairs, and the small diurnal oscillations in the signal are probably due to a combination of multipath effects and aliasing vertical tidal displacements into the horizontal direction (GPS poles are not perfectly vertical). However, this signal is small compared to the size of the rift widening.
Residual after removing a linear trend in the displacement between station pairs plotted along with a histogram of the seismicity (2 hour bin size) for 10 days around the first swarm in season 3. The jump in the transverse-to-rift component of the displacement coincides with the onset of the seismic swarm and only occurs in station pairs that span the rift.
The displacement appears to be entirely in the normal-to-rift direction (red), with little if any motion visible in the parallel-to-rift direction (black). In contrast, none of the baselines of station pairs that are on the same side of the rift (c–b, c–j) show any evidence of rapid widening during this time period. The magnitude of the signal is about 1 cm, large compared to the root-mean-square of the signal (0.3 cm) and, over 2 hours, gives a rough velocity ∼50 mm h−1. The displacement appears to have the same magnitude for all station pairs that span the rift and occurs simultaneously for all station pairs, including c–g which extends far ahead of the visible rift tip. Furthermore, the rift-widening signal appears to take place over a longer duration than that of the swarm in seismicity. This suggests that the seismicity is unlikely to be caused by blocks falling into the rift and wedging it open, as, if that were the case, we would expect the widening to precede the seismic swarm. These results are in excellent agreement with those determined during season 1 when we found three seismic swarms, each of which was accompanied by rapid rift widening that occurred over 1–2 hours with ∼1 cm magnitude. However, the increased number of GPS receivers allowed us to identify the signal in more station pairs and the increased sampling rate enabled us to increase the temporal resolution of this signal.
Discussion
Results from all three field seasons are summarized in Figure 10. We found three swarms of seismicity during seasons 1 and 2 and one smaller swarm during season 3. For each seismic swarm for which we have GPS data we find that the swarms are coincident with ∼1 cm of rift widening. The persistence of seismic swarms from season to season, accompanied by rapid rift widening, confirms our initial conclusion from season 1, that rift propagation is episodic (Reference Bassis, Coleman and FrickerBassis and others, 2005). The higher sampling rate of the GPS and greater sensitivity of our seismic network shows that the widening appears to last for several hours after the swarm ends, indicating that it may be partially related to post-propagating viscous relaxation processes within the rift or ice shelf.
Histogram showing the number of events per 3 hour bin for (a) season 1, (b) season 2 and (c) season 3. Each of the major swarms detected is shaded in gray for emphasis.
We also found that both the number of swarms and the number of events within each swarm decreased during season 3, resulting in a lower total rate of seismicity. One explanation for the decreased seismicity is that intrinsic or statistical variations in propagation rate – and hence seismicity – are not resolved by the short duration of our field campaigns. To test this it would be necessary to obtain continuous measurements over a longer period of time, ideally a full year. Another possibility is that the ‘true’ propagation rate of the rift is decreasing and this is reflected in the lower rate of seismicity we observed. To examine this further, we compared our field results with rift-propagation rates derived from Multi-angle Imaging Spectroradiometer (MISR) imagery for rifts T1 and T2 (Reference Fricker, Young, Coleman and BassisFricker and others, 2007). Figure 11 shows the rift lengths of T1 and T2 from Reference Fricker, Young, Coleman and BassisFricker and others (2007) extended through January 2006. Also shown is a quadratic fit to the time series of rift lengths for both rifts. Measurements are only available between October and March, when there is sufficient sunlight, causing a data gap over each winter. There is an overall decrease in the propagation rate for T2, which intriguingly corresponds to a simultaneous increase in rift T1’s propagation rate. This decrease in T2’s propagation rate is consistent with the results from our field studies, and the simultaneous increase in T1’s propagation rate suggests that the two rifts are not propagating independently of each other. In summary, there may be longer-term variations in seismicity but it appears that the decrease in seismicity we observed is related to changes in the long-term rate of rift propagation.
Time series of rift lengths for T1 and T2 derived from MISR. Upper panel shows rift length of T2. Lower panel shows rift length of T1. The error bars shown are more conservative than those used by Reference Fricker, Bassis and MinsterFricker and others (2005). Instead of using 1 pixel, we used the maximum decrease in rift length between successive points within each season (i.e. we assumed a measurement indicating that the rift closes due to measurement error).
The rate of long-term rift propagation is, at least partly, controlled by the background glaciological stress within the ice. This stress is changing as the rift lengthens, the ice front advances and the Loose Tooth rotates within its embayment. However, the timing of the decrease in seismicity also corresponds to the rift propagating into a suture zone between two flow bands, formed where two ice streams merged 500 km upstream (shown in Fig. 1). At the base of the ice shelf, along this same suture zone, lies a thick band of marine ice (Reference Fricker, Young, Allison and ColemanFricker and others, 2002). Since they form at the boundaries of merging ice streams, suture zones may contain many pre-existing fractures and shear bands that result in a decrease in the stress concentrated ahead of the rift tip, leading to a decrease in rift-propagation rate. Such a decrease in rift-propagation rate when the rift tip enters a suture zone is consistent with observations for T1, which in winter 2002 emerged through a suture zone and had a sudden spurt of growth (Reference Fricker, Young, Coleman and BassisFricker and others, 2007). It is possible that the cluster of events seen ahead of the T2 rift tip during season 3 is evidence that the stress is concentrated on the other side of the suture zone and the rift will begin to propagate more rapidly again once the rift jumps from its present position to the location of the cluster of seismicity. Another possible explanation for the slowdown of T2 is the increasing thickness of marine ice encountered in the rift’s path. If the thickness of marine ice is controlling the rate of propagation then we expect the rift would continue to slow down until the thickness of the marine ice decreases again. This leaves us with two competing hypotheses for the decrease in rift-propagation rate, with different predictions for the trend in rift-propagation rates over the next few years. As the band of marine ice does not exactly coincide with the suture zone, continued monitoring of propagation rates will enable us to directly test if either of these two hypotheses is correct or if the larger-scale stress within the ice primarily determines the speed at which rifts propagate.
Why is propagation episodic?
Our discovery that rift propagation is episodic is not very surprising. Many geophysical fracturing processes also display episodic behavior (e.g. earthquakes). One possibility is that propagation events might be ‘triggered’ by short-timescale environmental forcing (e.g. tides, winds, ocean swell). In a companion paper we will show that the timing of the propagation events does not correspond to any obvious environmental variable. There we will present calculations that show most environmental variables exert a stress that is small in comparison with the internal glaciological stress of the ice. Our observations, in combination with order-of-magnitude calculations, indicate that if environmental variables play a role in rift propagation, it is subtle.
Another possibility is that the episodic bursts of propagation are related to the crystal structure of ice. Fracture experiments on ice are known to produce fracture jump–arrest cycles in which cracks propagate forward in discrete jumps (Reference Rist, Sammonds and OerterRist and others, 2002). It is therefore tempting to conclude that the episodic bursts of propagation that we observe are a larger-scale manifestation of the fracture propagation seen in laboratory experiments. This interpretation, however, is problematic because rifts are tens of kilometers long and the effect of inhomogeneities such as grain boundaries and brine pockets that are responsible for the laboratory-scale jump–arrest behavior have an increasingly small effect at large length scales. In linear elastic fracture mechanics (LEFM), fracture occurs when the stress concentrated near the tip of a crack exceeds a material property, called the fracture toughness (denoted K IC). Mathematically, the critical stress (σ crit) at which fracture occurs decreases as the square root of the crack length (L)
If we further assume that material inhomogeneities result in a spatially variable fracture toughness, we can express the fracture toughness as the sum of a constant fracture toughness, 
with a smaller variable term, Δ. Then Equation (2) becomes
Figure 12 shows the critical stress as a function of rift length for a variety of values of Δ for K IC = 150 kPa m1/2 (Reference Rist, Sammonds and OerterRist and others, 2002). Even for variation in the fracture toughness as high as 100%, the stress necessary to propagate a 5 km long rift is less than 5 kPa. In comparison, typical values for the glaciological stress are on the order of hundreds of kilopascals (Reference PatersonPaterson, 1994). This illustrates two things. First, LEFM predicts that small-scale material inhomogeneities should have an increasingly small effect on rift propagation for large rifts. Second, large fractures should be unstable, a familiar result of LEFM (Reference LawnLawn, 1993).
Minimum stress necessary to propagate a rift as a function of rift length as determined from Equation (2) for K IC = 150 kPa m1/2 and Δ ranging from 0 to 150 kPa m1/2.
It is possible that the rift is loaded in such a way that the stress on the rift decreases with increasing rift length, resulting in a stable loading configuration. Reference Larour and RignotLarour and others (2004) suggest the loading should be modeled as a double cantilevered beam. Instead, we suggest a more realistic explanation for episodic propagation is that the rift propagates by the formation of a series of micro- and mesoscale cracks that extend ahead of the rift tip. As more and more of these cracks form and propagate they eventually coalesce, extending the length of the rift (a schematic diagram of this process is shown in Fig. 13). Conceptually, this explains the background seismicity and swarms that we observed with the seismometers, as well as the slow (over the course of several hours) rift widening that we observed with the GPS receivers. It also follows a qualitative pattern similar to that observed in subcritical crack growth in laboratory experiments (Reference Atkinson and MeredithAtkinson and Meredith, 1987) and postulated for lithospheric rift propagation (Reference Floyd, Tolstoy and MutterFloyd and others, 2002) and for crevasse formation (Reference WeissWeiss, 2004). In some experiments in rock, acoustic emissions caused by the initiation of small micro-cracks are observed over a wide region (Reference Lockner, Byerlee, Kuksenko, Ponomarev and SidorinLockner and others, 1992). At the beginning of the process single, isolated micro-cracks form. These micro-cracks grow and multiply, resulting in an increase in the density of cracks per unit volume. Eventually, a critical crack density is reached and the micro-cracks begin to coalesce into a single main fracture. Furthermore, the production of micro-cracks reduces the stress concentrated near the tip of the cracks. This has led some researchers to propose scaling laws where the fracture toughness scales with crack length, resulting in potentially stable propagation even for very large fractures (e.g. Reference Floyd, Tolstoy and MutterFloyd and others, 2002).
Sketch of rift propagation by micro- and mesoscale crack initiation. In the first stage (a), a series of small and medium-sized cracks initiate around the rift tip. When the density of these cracks approaches a critical value, these smaller cracks begin to merge and create new rift surface (b), leading to the seismic swarms that we identified (c).
What controls the recurrence intervals of episodic propagation?
Earthquake mechanics may provide an appropriate analogy to rifting that can be exploited to determine what controls the intervals in episodic propagation. For example, Reference Shimazaki and NakataShimazaki and Nakata (1980) propose four different models through which stress accumulates and is released through earthquake rupturing events. These four models correspond to conditions in which: (1) both the magnitude and time of major rupture events are predictable; (2) the magnitude of each major rupture event is predictable but the time between events is variable; (3) the time between major rupture events is predictable but the magnitude is not; and (4) neither the time nor the magnitude is predictable. Although these models were proposed in the context of shear cracks, taking a liberal definition of slip to include our rift-widening events, we can qualitatively apply these models of earthquake recurrence intervals to intervals between rift-propagation events. In contrast to earthquakes, the recurrence relation between episodic propagation events is short (i.e. several weeks) and this enables us to immediately begin testing various models of stress accumulation. If we take the rift-widening signal as a measure of the size of each event, then each event is approximately the same size, ∼1 cm. Figure 14 shows slip (i.e. rift widening) as a function of time for seasons 1 and 2, where we have assumed that the slip between rift-widening events is small. As we have only three events it is difficult to conclusively rule out any of the models. However, with the limited data we have, the best-fitting model is slip-predictable (model 2). The slopes of the best-fitting line for each of the seasons are statistically identical. This implies that the critical stress at which each propagation event occurs is constant, but that the stress drop varies between events. Unfortunately, this fit is only marginally better than any of the other three models. Nonetheless, our results suggest that we would detect enough propagation events to decide which model is appropriate if we were able to deploy instruments over a full annual cycle. This would help substantially in understanding how stress accumulates and is released at the tip of ice-shelf rifts.
Recurrence times of rift-widening events for seasons 1 and 2. The dashed black line indicates the best-fitting slip-predictable model for season 1, while the dashed blue line indicates the best-fitting slip-predictable model for season 2. Slopes for the two seasons are statistically identical.
Conclusions
In this study we have confirmed our earlier discovery from the 2002/03 field season that rift propagation is episodic (Reference Bassis, Coleman and FrickerBassis and others, 2005). Swarms, containing 100–175 seismic events, accompanied by rapid rift widening of ∼1 cm are seen in all three field seasons. Swarms have typical recurrence intervals ranging between 10 and 24 days. In addition, we see a substantial decrease in the rate of seismicity detected over the 3 year period that may indicate rift T2 is slowing down. This is confirmed by analysis of rift-propagation rates using satellite imagery, which also shows that the slowdown corresponds to the rift propagating into a suture zone where two ice streams merged. The location of this suture zone is coincident with the easternmost of two thick bands of marine ice under the ice shelf. This leads us to postulate that the slowdown is either related to preexisting fractures and shear bands within the suture zone causing a decrease in the stress concentrated at the tip of the rift or is related to the rift propagating into a thicker band of marine ice. It also suggests that the rift may accelerate once the tip emerges through to the other side of the suture zone, a detail that our current field deployment is designed to detect. Moreover, because the widths of the suture zone and marine ice band are not the same, the timing of any subsequent acceleration in propagation rate will indicate whether the slowdown is caused by the suture zone, marine ice or changes in the background glacio-logical stress.
The dense spacing in the seismic network we deployed enabled us to locate thousands of seismic events originating from the rift, providing the first detailed map of rift-related seismicity. Our locations show that while events cluster around the rift axis, there is also a tail of events up-rift (towards the triple junction), possibly caused by ice blocks falling into the rift from the sides or rearrangement of blocks of ice within the rift. Relative magnitudes are distributed similarly to the Gutenberg–Richter law for earthquakes. However, the temporal behavior of the icequakes is quite different. Unlike earthquakes, there are no foreshock and aftershock sequences. Instead, swarms consist of a series of events of different sizes. This leads us to suggest that rift propagation occurs as a series of micro- and meso-scale cracks ahead of the rift tip. The swarms are the result of coalescence of these smaller cracks once a critical density is reached. This study represents an important first step in understanding the details of the mesoscale physics that control rift propagation over moderate spatio-temporal scales.
Acknowledgements
This work was supported by the US National Science Foundation Office for Polar Programs through grant NSF-0337838, by the Australian Government Antarctic Division and through Australian Antarctic Science grants AAS 2338, 2686 and Australian Research Council grant DP0666733. We thank Raytheon Polar Services and the Australian Government Antarctic Division for logistical support during the fieldwork and T. Sprent and A. Hull for invaluable assistance in the field. Use of the RTD software was provided by the Scripps Institution of Oceanography, under license to Geodetics, Inc. We thank the University Navstar Consortium (UNAVCO) and PASSCAL for loan of the scientific equipment. We also thank T.H. Jacka, B. Kulessa and an anonymous reviewer whose comments substantially improved the clarity of the manuscript.
Appendix A: Event Detection Algorithm
During both seasons 2 and 3 we identified thousands of seismic events with durations lasting from less than 1 s to hundreds of seconds. While hand-picking each of the seismic arrivals would yield more accurate results, the high volume of data necessitated automating the procedure. We developed a variation of the short-term/long-term average (STA/LTA) algorithm specifically designed to detect ice-quakes that we observed with our seismometers (e.g. Reference Baer and KradolferBaer and Kradolfer, 1987). Our implementation of this algorithm consisted of the following steps.
1. Bandpass-filter each component to get a new signal: sn (t), where n ranges from 1 to 3 and indicates the north, east and up components of the seismogram. Most of the energy in the waveforms is in the range 5–50 Hz, but following experiments with a wide variety of bandpass windows we found the best frequency range to be 5–35 Hz.
2. Calculate the envelope function for each component. The envelope function, defined as the analytic signal multiplied by its complex conjugate, forms the positive outline that encloses each wavepacket:
where is the Hilbert transform of the seismogram.3. Calculate the combined envelope function for all channels by taking the mean of the envelope functions for each component:
(As the aperture of our network is only a few kilometers, events that originate within the network will have small separation between P-waves and S-waves. We therefore do not try to distinguish between different phases of the arrival.)4. Use the combined envelope function to calculate the average of the signal over a long interval and the average of the signal over a shorter interval by convolving the combined envelope function with a broad and narrow Gaussian and then taking the ratio of the signal averages. We found that a long-term interval of 1 hour and a shortterm interval of ∼0.25 s gave the best results.
5. When the ratio of the STA/LTA exceeds a predefined threshold, Smin , for an interval longer than a predefined time period, we designate this an event. The time at which the envelope exceeds Smin determines the arrival time. We also record the time at which the STA/LTA exceeds both lower and higher thresholds to determine how emergent or impulsive each waveform is. This also gives an estimate of the uncertainty of our arrival time picks. To prevent detecting the same event twice, we include a latency period of 1.5 s before the algorithm is allowed to ‘trigger’on the next event.
where 
is the Hilbert transform of the seismogram.
(As the aperture of our network is only a few kilometers, events that originate within the network will have small separation between P-waves and S-waves. We therefore do not try to distinguish between different phases of the arrival.)There are a number of tunable parameters in this algorithm. We found that the algorithm was most sensitive to the combination of S min and the short-term averaging interval. We tuned the combination of these two variables to best match the arrival times for 200 events that were picked manually.
Appendix B: Event Association Algorithm
We use a simple approach to associate arrivals detected at multiple stations with common events.
1. Sort all of the arrival times into a list with the earliest detections occurring first, regardless of station.
2. Calculate a matrix of travel times between the ith and jth stations (δij ), assuming a constant velocity. This is an upper bound on the travel time differences between events recorded at two different stations.
3. For each arrival time on the arrival time list we search for arrival times at other stations that fall within
where we had labeled the ith arrival time in the list We associate all arrival times that we find with the master station. These arrival times are then flagged and cannot be associated with other events.4. If we find multiple arrival times at the same station, we decrease our window size until we only have one arrival for each station. For events that are closely spaced, this may result in an incorrect association, but most events are spaced more than several seconds apart. Events separated by <2 s constitute <1% of all detected events.
where we had labeled the ith arrival time in the list 
We associate all arrival times that we find with the master station. These arrival times are then flagged and cannot be associated with other events.Once our catalog of events is constructed, we remove all known earthquakes detected by our network using the US Geological Survey earthquake catalog (http://neic.usgs.gov/neis/epic/epic.html).
Appendix C: Event-location Algorithm
We solve for locations using a three-dimensional grid-location algorithm (see, e.g., Reference ShearerShearer, 1997). To do this we first define a three-dimensional grid of points with horizontal range 10 km, depth 400 m (the ice thickness) and spacing between grid nodes 25 m, centered on the location of the rift tip identified in the field. We then compute travel times from each gridpoint to each of our stations using a seismic velocity ranging from 2000 to 3000 ms−1, in steps of 250ms−1. Using smaller step sizes had a negligible effect on the results. Because we assume a constant velocity, the ray paths are straight lines and the predicted travel time from a source located at Cartesian coordinates (xm, ym, zm ) with an origin time tm to the ith station is the Euclidean distance between source and receiver divided by the velocity, V:
All of our stations are on the surface, therefore zi = 0. The depth of the event, zm , is restricted to occur within the ice thickness, ∼400 m. For each event, the observed arrival times are compared to the predicted travel times at each gridpoint and the best-fitting location is determined from the smallest residual. Following Reference ShearerShearer (1997) we use the L 1 norm (i.e. absolute value of the differences between predicted and observed travel times)
The L 1 norm assigns less weight to outliers than the L 2 norm used in traditional least-squares fitting and is more robust for this application. As pointed out in Reference ShearerShearer (1999), it is not necessary to search for the best-fitting origin time for each gridpoint. The origin time is simply the median or mean of the residuals (for the L 1 or L 2 norm, respectively). We then find the seismic velocity that minimizes the L 1 norm of the residuals for all events at all stations. Although we use a constant velocity, our location algorithm is flexible enough to use a more complicated three-dimensional velocity model. Experimenting with more complicated depth-dependent velocity models did not yield significantly different results for the horizontal locations. However, the depth of each event was highly dependent on the velocity model.
