Special Research Articles

Rupture Process of the 2004 Sumatra-Andaman Earthquake

See allHide authors and affiliations

Science  20 May 2005:
Vol. 308, Issue 5725, pp. 1133-1139
DOI: 10.1126/science.1112260


The 26 December 2004 Sumatra-Andaman earthquake initiated slowly, with small slip and a slow rupture speed for the first 40 to 60 seconds. Then the rupture expanded at a speed of about 2.5 kilometers per second toward the north northwest, extending 1200 to 1300 kilometers along the Andaman trough. Peak displacements reached ∼15 meters along a 600-kilometer segment of the plate boundary offshore of northwestern Sumatra and the southern Nicobar islands. Slip was less in the northern 400 to 500 kilometers of the aftershock zone, and at least some slip in that region may have occurred on a time scale beyond the seismic band.

Seismic waves are excited by rapid and varying sliding motions that initiate with a frictional instability. Slip begins as the rupture front spreads across the fault with a velocity usually less than the ambient shear wave speed. Both rupture propagation and local slip history (the temporal variation and total slip at a particular position on a fault) influence the frequency and strength of radiated seismic waves. Different positions on the fault generally have different displacement histories, including variations in the rate and amount of slip. Seismic waves sense these differences, and by using ground motions observed far from the source seismologists can reconstruct the spatial and temporal slip history of faulting.

Several phenomena affect seismic wave excitation during faulting. One is the stress drop at the rupture front. As the rupture front expands, short-period P and S waves are generated from the local stress reduction. For large events, these waves can be used to map the earthquake's rupture expansion. The speed of rupture front propagation, which can be related to the energy partitioning during the faulting process, is an important quantity. The potential energy released during earthquakes is partitioned into seismic radiation, mechanical processes such as creation of fractures, and frictional heat (1). The amount of heat generated by frictional processes during the rupture depends on the absolute stress, total slip, and rupture area. The partitioning of energy between mechanical processes and seismic radiation varies from earthquake to earthquake and provides one method of classifying different faulting processes. Fast ruptures can be associated with a relatively large fraction of seismically radiated energy (1, 2). For many well-studied earthquakes, the rupture speed is 70 to 95% of the shear wave velocity, but important variations have been observed as complex ruptures cross fault-segment boundaries (3). Another important observation is the spatial pattern of slip in large earthquakes. For many shallow earthquakes, slip near the hypocenter is relatively small, indicating to some extent that the earthquake began at a weak region and grew into a much larger event (1). These observations are extracted from analysis of the seismic wave field. The 26 December 2004 Sumatra-Andaman and the 28 March 2005 earthquakes (4) produced the most extensive high-quality broadband seismic data ever recorded for great earthquakes. Here, we exploit signals across a broad bandwidth and every part of the seismic wave field to construct an integrated seismic view of these earthquake ruptures. Our focus is on the first and larger of the two events.

Short-period P-wave directivity. Short-period P-wave radiation (5) for large earthquakes provides direct information about the rupture front propagation. The energy radiated by an expanding rupture front can be observed with the use of the global seismic networks (6) or regional seismic and hydroacoustic arrays (79). One of the simplest measures that can be made is the duration of short-period P-wave radiation from the source region (10, 11). For a long-duration earthquake, a major challenge for P-wave analysis is the interference of later-arriving seismic waves reflected from the surface and discontinuities in the Earth with P waves radiated from later portions of the rupture. Fortunately, most secondary phases involve additional path segments in the highly attenuating upper mantle, and their short-period content is suppressed (12). Applying a high-pass filter can reduce the effects of secondary arrivals. The durations of short-period P waves will be shorter in the direction of rupture propagation and longer in the direction away from the moving source (the rupture front). Data for the Sumatra-Andaman earthquake (Fig. 1) indicate a north-northwest rupture propagation with a speed of about 2.5 km/s and an overall fault length of 1200 to 1300 km, a length consistent with the aftershock distribution (4).

Fig. 1.

High-frequency seismograms (A) and envelopes (B) as a function of azimuth from the epicenter. The signal duration varies smoothly as a function of azimuth, and three relatively energetic amplitude bursts are identified with the gray shading.

The amplitude of the short-period waveforms generated during the rupture also varied about a relatively uniform level. At least three large (from 50 to 150 s, 280 to 340 s, and 450 to 500 s) and several additional seismic energy bursts are apparent (Fig. 1). For a uniform rupture speed, the large bursts are located near latitudes of 4° to 6°N, 8° to 10°N, and 12° to 13.75°N. Ishii et al. (8), with use of the Japanese Hi-net seismic network, identified two large bursts of P-wave energy at about 80 and 300 s after the origin time and a mean rupture speed of about 2.8 km/s. Their estimate of the rupture speed is slightly higher than the one we proposed, but the Hi-net bursts of energy are consistent with the first two features observed in the global network observations. Observations from three hydroacoustic arrays across the Indian Ocean were triangulated to identify the source of T phases excited by the Sumatra-Andaman rupture. The T-phase source propagated northward at a speed of about 2 km/s for the first 600 km of the rupture. A secondary more northern T-phase source began before the first had reached the Nicobar islands (9).

Long-period Rayleigh wave directivity. For earthquakes up to about magnitude 7.5, surface waves with periods greater than 20 s usually sense the entire rupture process, such that the surface wave magnitude, MS, based on 20-s-period surface waves, is indicative of the entire slip event. But for the 2004 earthquake, even 500-s-period surface waves failed to sense the full slip process. Surface-wave phase velocities are close to rupture speeds, which enhances their sensitivity to rupture directivity (azimuthal time delays related to spatial and temporal source finiteness), so comparison of azimuthal patterns in the surface wave amplitudes at different periods provides useful constraints on the overall slip distribution on the fault. Surface waves disperse as they propagate, producing complex waveforms and amplitudes that depend on distance as well as source radiation pattern. For the 2004 event, the faulting geometry inferred from 300- to 500-s-period surface waves is nearly pure thrusting on a shallow plane (4), which produces a predominantly two-lobed azimuthal Rayleigh-wave radiation pattern (modified somewhat by the north-northwest rupture propagation), and a four-lobed Love-wave radiation pattern. We focus on Rayleigh waves because of their simpler pattern.

Rayleigh-wave amplitudes were measured with the use of R1 for shorter periods and R3 for longer periods (which are more robust when measured with longer time windows) (13). To keep the analysis straightforward, we account for effects of dispersion, fault geometry, and distance by making synthetic seismograms for Rayleigh waves. For periods shorter than 800 s, we used synthetic seismograms computed with the spectral-element method (14) for a three-dimensional (3D) Earth model composed of the mantle model s20rts (15), the model Crust 2.0 (16), and the topography from ETOPO5. For periods longer than 800 s, we used synthetics computed by summation of normal modes. The 3D simulations were performed on the Caltech Division of Geological and Planetary Sciences Dell cluster. Initially, the synthetics are for a point source model with no finiteness using the source parameters of the Harvard Centroid Moment Tensor (CMT) solution (4), and we simply measure the amplitude ratio (observed/predicted) and the time shift (observed – predicted) in various pass bands and plot them as a function of azimuth (Fig. 2, top left).

Fig. 2.

Time-domain Rayleigh wave amplitude and phase-shift measurements as a function of azimuth from the epicenter. Observations for a period range of 150 to 1500 s are shown. All periods show the effects of directivity along the rupture direction of 310° to 330°N. (A) The observations are compared with predictions for the Harvard CMT solution for amplitude and time shift. (B) A similar comparison with predictions computed for model III. The smaller scatter for both amplitude and time shift indicates that model III can account for the first-order long-period surface wave observations.

If there were no effects of finiteness, all the amplitude ratios would be unity for all periods. For a uniformly slipping fault with unilateral propagation, one would expect directivity effects to be stronger at shorter periods with larger amplitudes in the direction of rupture. For 150-s and 300-s waves, the amplitudes are enhanced along the rupture direction (around 330° N), as expected. However, the amplitude ratios and the strength of the azimuthal pattern unexpectedly do not decrease much at longer periods, up to at least 1500 s. This behavior cannot be interpreted by using uniform slip along the fault. Modeling of simple propagating ruptures suggests that the observations for periods shorter than 600 s are compatible with north-northwest propagation of a rupture at about 2.5 to 3 km/s for 400 to 600 km from the southern end of the fault, with overall slip averaging about 9 m. This slip accounts for the seismic moment of the Harvard CMT solution (4). But, longer-period surface wave amplitudes are only partly accounted for by this model. The pattern of time shifts (Fig. 2, top) supports this conclusion. At periods of 150 and 300 s, the time shift is relatively small, but, at periods longer than 600 s, the time shift is systematically positive; i.e., observed waves are later than predicted by the CMT. These observations suggest that additional slip extended in either time, space, or both is required to explain the surface-wave data.

Surface-wave moment rate functions. Intermediate-period surface waves (80 to 500 s) were the dominant signals globally recorded from this earthquake (4, 17). Here, we remove the dispersive wave propagation effects by deconvolving the data by point-source normal-mode synthetic seismograms computed for the fault strike and rake of the CMT solution (4), but with a slightly larger dip (11°)(18, 19). We assumed a standard Earth model (20), the USGS origin time and epicenter (4), and a source depth of 15 km. We used two different deconvolution approaches, a frequency-domain method (21) and a time-domain method (22). The results are similar, so we discuss the time-domain results here [see (12) for the frequency-domain estimates]. The result of each deconvolution is a source time function, indicating the time history of Rayleigh-wave radiation at the corresponding azimuth.

Rayleigh-wave source time functions strongly vary with azimuth relative to 330°N [the strike of the Harvard CMT mechanism (4)] (Fig. 3). We thus arranged the data by directivity parameter (23), which transforms features with a sinusoidal azimuthal moveout to a linear moveout. We combined 172 R1 source time functions into 30 average source time functions binned by directivity parameter. The observations are not uniformly distributed with azimuth and reflect the paucity of seismic stations in the Southern Hemisphere. For azimuths close to the direction of rupture (near the bottom of the plot), the observations show a relatively narrow, large-amplitude pulse, and the pulse width increases with decreasing directivity parameter. As the azimuth relative to rupture increases, several smaller amplitude pulses shift systematically to later relative times, the peak amplitude decreases, and the overall duration increases. These pulses are not observable on all the recordings but can be seen clearly on several. Their true temporal relation can be seen on the time functions nearly perpendicular to the rupture. The moveout and pulse shape changes are related to the spatiotemporal slip distribution and rupture propagation as sensed by these waves.

Fig. 3.

Deconvolved Rayleigh wave source time function estimates (black lines). The source functions, obtained by water-level deconvolution, are arranged in order of increasing directivity parameter. The equivalent azimuth relative to the rupture direction (330°N) is shown to the left; the number of observations stacked in each bin is shown to the right. Lines a, b, c, and d identify discrete phases that can be tracked at least across as least several source functions. Predictions from the IRT inversion are shown in gray lines. The fits are best where the numbers of data are large. STFs, source time functions.

These intermediate-period surface-wave time functions can be formally inverted for a 1D model of slip distribution along the fault by using an inverse Radon transform (IRT) (24, 25). This is a first step toward 2D maps of the seismic moment distribution but has the added flexibility that the time dependence of slip does not have to be prescribed. The results (Fig. 4) show that the initial 60 s of the rupture had little moment and that the largest moment occurred about 100 to 200 km northward, about 1 min later. Once the large slip event began, the slip radiating in this bandwidth propagated northward at about 3 km/s. After a few hundred kilometers, the seismic energy in this bandwidth began to decrease steadily, suggesting that most of the rapid slip occurred in the first 700 to 800 km of the rupture. Beyond about 800 km, some slip is imaged, but it is low amplitude (<2 m) and not well resolved. Given that the short-period energy and the longest-period surface wave amplitudes require a rupture propagating farther to the north, we infer that the slip continued on to the central Andaman islands region, about 1300 km north of the hypocenter, but that it ruptured smoothly and radiated relatively little energy in this bandwidth.

Fig. 4.

(Top) Moment-rate density image showing the variation in seismic moment with time and with position along a 1D fault striking 330°N constructed with the use of Rayleigh waves (R1) in the period range from 80 to 400 s. (Bottom) Curve shows the moment-rate time function in units of 1020 N·m/s. The image is the result of 25 local search inversions, stacked to reduce truncated search artifacts. Stacking also smooths the image and moment-rate estimate. Slip radiating seismic energy in the intermediate-period surface wave band (80 to 100 s) clearly extends to about 750 km and was concentrated in the first 250 km. Fits are shown in Fig. 3.

2D analyses. The rupture area outlined by reverse-faulting aftershocks (4) has a substantial width (150 to 200 km), and seismic data are sensitive to variable slip along both the strike and dip directions. We applied three different modeling approaches to estimate the 2D moment distribution of the Sumatra-Andaman earthquake. The use of three methods allows us to identify and to focus on the robust features of the rupture resolved by seismic observations and to assess the sensitivity of those features to necessary choices of fault parameterization, temporal moment distribution parameterization, and seismic data weighting. The results presented here build upon rapid determinations of faulting conducted soon after the event (26).

Each method involves a different parameterization of the source as a function of time, but all three parameterize the slip by using planar faults divided into a grid of subfaults ranging in area from 16 km by 20 km to 50 km by 50 km. Each method includes some constraints on the propagation of the rupture front, although we made efforts to reduce the impact of these on the results. To parameterize the history of sliding on the fault, we used a varying rise time for each subfault or multiple time windows during which slip at the node is allowed to accumulate. In all methods, we restricted slip to occur during a specified time window (up to 180 s long) after the rupture front passes. The slip direction was allowed to vary about the plate motion or the Harvard CMT rake direction. For each approach, we constructed a set of linear or linearized equations relating the moment rate history of each subfault as a function of time to the observed seismograms. Methods used include linear programming to minimize the robust L1 norm for SH waves (27) (model I), a least-squares inversion of group velocity–windowed intermediate-period (80 to 300 s) surface waves combined with regional long-period seismograms (28) (model II), and a search-based method to minimize wavelet coefficients of the observed teleseismic P and SH waves and complete long-period three-component regional seismograms (2931) (model III). The nonlinear dependence on rupture propagation is handled with constraints on the rupture front velocity or by using a grid search specifying a rupture front propagation that is perturbed during the inversion. Waveform fits for each model are included in (12).

The first model (Fig. 5A) was constructed with the use of the method described in (27) and 20 teleseismic SH waveforms, filtered to include periods shorter than 120 s. The rupture surface was parameterized as two fault segments: the first having a strike of 329° and a dip of 8° and the second having a strike of 333° and a dip of 7° (based on the mechanism of the 29 December 2004 MW = 6.0 aftershock). Cells on each fault segment were 40 km by 40 km, and the source time function was divided into 15 12-s time steps. Each cell was allowed to slip in each time step except that the rupture front could not propagate faster than a P wave from the hypocenter. The total moment was constrained to be similar to the Harvard CMT moment of 4.0 × 1022 N·m; however, because the true rupture duration is longer than 180 s, any additional moment after 180 s is forced onto to the final time step, which is then discarded.

Fig. 5.

(A) Fault slip 168s after rupture initiation estimated by using 20 azimuthally distributed teleseismic SH waveforms (Δ ∼ 45° to 85°). The rupture models consists of two faults, the first having a strike of 329° and a dip of 8° and the second having a strike of 333° and a dip of 7° (based on the mechanism of the 29 December 2004 MW = 6.0 aftershock). (B) Slip distribution from method II. The reliance on intermediate-period surface waves and long-period seismograms reduces the detail imaged in the rupture but provides a first-order view of the slip distribution. (C) Slip distribution of finite fault model III using teleseismic body waves (5 to 200 s), intermediate-period three-component regional waves (50 to 500 s), and long-period teleseismic waves (250 to 2000 s). The surface projections of three fault segments are colored on the basis of the slip amplitude. The black thick and thin lines delineate the trench mapped from the ETOPO2 and 50-km iso-depth slab contour. The aftershocks (Ml > 5) downloaded from the National Earthquake Information Center are indicated by black dots. Waveform fits for each model can be found in the electronic supplements. Slip of the 28 March 2005 event is outlined with a dashed line. Area ruptured during the 28 March 2005 event is outlined with a dashed line.

The solution favors at least 10 m of slip near the hypocenter. This amount is consistent with reported uplift of Simeulue Island (32). A second region of large slip, approaching 20 m, is located southeast and west of Great and Little Nicobar islands. A third patch of 5 to 10 m of slip is located near 4°N, but tests show that the data allow more slip in this region. In general, slip is concentrated along the deeper parts of the megathrust in this solution, suggesting that most of the seismic energy in the shorter period SH waves originated from relatively deep on the fault. The model's rupture started slowly near the hypocenter (∼1.3 km/s) but accelerated up to 3.3 km/s toward the Nicobar islands in the north.

The second model (Fig. 5B) was obtained by using a least-squares inversion of regional long-period seismograms in the period range from 100 to 3000 s and regional and teleseismic surface waves in the period range from 80 to 300 s. The surface waves were modeled with use of aspherical Earth model corrections computed for the Harvard phase velocity model (33). The point-source grid spacing was 50 km by 50 km, each node having source time functions with a duration of 40 s. The largest slip predicted by the model was located between about 3° and 6°N, spread over much of the megathrust width, but with larger slip deeper. Slip near the hypocenter was relatively large but decreased quickly in the surrounding 100 km or so. As with the SH waveform results, a second area of larger slip was located northwest of Great and Little Nicobar islands. Slip decreased north of 9° to 10°N, but the model suggest that slip continued to the north into the Andaman islands region, and the total moment was about 6.5 × 1022 N·m, about 1.5 times larger than the Harvard CMT moment and approaching that estimated by using normal modes (34).

The third model (Fig. 5C) was constructed with the use of teleseismic body waves (20 to 200 s), intermediate-period three-component regional seismograms (50 to 500 s), and long-period teleseismic seismograms dominated by R1 and R2 phases (250 to 2000 s). The rupture surface was approximated with the use of three fault segments with strikes approximating the local trench axis. The fault segment dip angles were approximated with the use of seismicity-based slab contours (35), which imply that the dip of the fault segments increases from south to north at 12°, 15°, and 17.5°. Subfault spacing was 20 km along strike and 16 km along dip. The intermediate- and long-period waveforms were computed by using normal mode summation but calibrated by spectral element synthetics computed as described above. A global optimization simulated annealing algorithm was used to estimate the slip amplitude and direction (the rake angle) as well as rupture initiation, rise, and rupture cessation times for each subfault. Rupture initiation times were allowed to vary up to ±150 s from the time that a rupture propagating at 2 km/s would pass the subfault. The 2 km/s average rupture speed was estimated with the use of multiple inversions but is not tightly constrained by the data.

In this model, the accumulated slip across the rupture surface composed of three planar faults (Fig. 5C) lasted for 550 s and produced a total moment of 6.5 × 1022 N·m, which gives a moment magnitude of Mw = 9.1. The model implies that slip was primarily concentrated south of 9.5°N, but slip extended northward into the Andaman Island regions. The area of largest slip is consistent with the surface-wave IRT results as well as model II, and these results all find a decrease in slip (radiating intermediate-band surface waves) along the first 750-km length of the rupture. The region of largest slip extends from about 3°N to about 6°N and includes substantial slip across the entire megathrust width. This is consistent with the large peak in the global surface-wave moment-rate functions. Slip is generally concentrated along the lower half of the megathrust, consistent with the other methods. In contrast with the SH-wave model (Fig. 5A), this model implies that slip near the hypocentral region was relatively low. A second region of strong slip is located west of northern Great Nicobar and Littler Nicobar islands, which matches the SH body-wave results. This model predicts uplift values between 1 to 5 m across a region with dimensions of 900 km by 100 km from the epicenter near 3°N to about 10°N. Uplift is a maximum near the trench between 4°N and 5°N (near-source surface displacement and movies of regional and global seismic velocities predicted by this model are included in the supplemental online material).

Slip maps for the 28 March 2005 (CMT MW = 8.6) event are shown are shown in Fig. 5, B and C. The peak slip (∼5 to 6 m) in both models is located near the hypocenter, and both models include rupture directed primarily to the southeast. The model slip is concentrated between about 20- to 40-km depth, which helps explain the smaller tsunami generated by this event. The substantial outer arc high in this part of the subduction zone also likely played a role in the tsunami generation (36). The model rupture durations were ∼150 s, and both models map the most substantial slip in the region previously identified at the 1861 earthquake rupture (4). Little, if any, slip penetrated into the 1833 rupture zone, which is a possible site for the next large event in the region.

Discussion. A number of general features are apparent from the directivity observations and fault rupture models. Like the 1960 Chile, the 1964 Alaska, and the 1952 Kamchatka earthquakes, the 2004 Sumatra-Andaman event ruptured largely unilaterally. Little aftershock activity penetrated south of the epicenter until the 28 March 2005 earthquake. Although there are differences in the rupture models described above, the first-order attributes of the rapid-slip faulting are well established. The moment rate functions (the combined effects of slip and rupture area expansion as a function of time) show that the fault sliding began relatively slowly but grew rapidly after about 40 to 60 s as large amounts of slip occurred off the west coast of Sumatra between 3°N and 4°N. The rapid increase in moment rate and the accompanying burst of short-period energy (Fig. 6) suggest the possible failure of a relatively strong section of the megathrust at that time. Slip amplitudes in the region are also the largest anywhere on the fault, approaching 15 m offshore of Sumatra. Rupture and slip continued to the north, but, after about 180 s, the moment rate decreased gradually to relatively low levels by about 450 to 600 s. The slip models obtained from inversions of body and surface waves (models 2 and 3) include gradually decreasing slip extending to 13° to 14°N. These models match the low-order normal-mode amplitudes to within about 10% (34). The lower panels in Fig. 2 show the predictions of model III (Fig. 5C) on the long-period Rayleigh wave directivity observations. The spread of the measurements in the top figures, both amplitude ratio and time shift, is greatly reduced, the amplitude ratios are relatively constant and close to unity, and the time shift is close to zero regardless of the period and azimuth (12). We conclude that the tapered slip between the Nicobar and Andaman islands is responsible for the observed azimuthal patterns of amplitude ratio and time shifts of Rayleigh waves.

Fig. 6.

Moment rate functions from each of the four rupture imaging methods in our analyses. From top to bottom, surface-wave IRT imaging and finite-fault inversion results using methods I, II, and III. All are presented on the same amplitude scale, and the seismic moment is listed above and to the left of each signal. For method I, the moment shown is the value reached 168 s after rupture commencement. Peak moment rates are in the range of 4 × 1020 N·m/s. Models constructed with the use of body waves are generally higher frequency. The surfacewave models recover only the smoother components of the rupture. The numbers at the bottom identify the apparent times of high-frequency energy bursts.

This faulting model suggests a relationship between megathrust coupling and rupture velocity and/or slip rate: The results indicate that the fault was well-coupled in the south, somewhat less coupled in the central portion, and weakly coupled in the north of the rupture zone. The subducting slab dip angle, age, and plate motion obliquity all increase from the southern (Sumatra) segment to northern (Andaman) segments of the rupture (4), perhaps contributing to reduction of interplate coupling as a function of distance northward. The reduction of slip just north of the Great Nicobar Island coincides with a northward rotation of the trench, and the rupture terminated in a region where the trench is parallel with the interplate motion (or even extensional) (4).

Although our models explain seismological data ranging from body waves to the gravest normal mode (period of 54 min) satisfactorily, the slip in the models to the north of 8°Nistoo small to explain global positioning system (GPS) displacements observed in the Nicobar Island (1 to 2 m vertical and 5 m horizontal) and the Andaman Island (1 to 2 m vertical and 3 m horizontal) (37). If we are to explain the deformation of the islands with the megathrust fault model estimated in slip inversions, we must increase the slip in the section north of 8°N by a factor of 2 to 3 (fig. S1). However, adding rapid slip of this magnitude considerably reduces the fit to the normal-mode amplitudes. Thus, most of this additional slip was probably slow and occurred at a time scale beyond the seismic band. More detailed analyses of tsunami, normal mode, and GPS data will be required to resolve the time scale of this additional slip.

Supporting Online Material


Materials and Methods

Figs. S1 to S13

Table S1

Movies S1 to S3

References and Notes

View Abstract

Stay Connected to Science

Navigate This Article