Introduction

The giant Tohoku earthquake on March 11, 2011 occurred off the eastern coast of Honshu in Japan and caused extensive damages and casualties, for the most due to huge tsunami waves that ravaged the coastal areas in the Tohoku region. Location, depth and thrust faulting mechanism indicate that the event originated along a crustal portion of the subduction zone at the boundary between the Pacific and Okhotsk (North American) plates1. As shown by inland GPS observations, in this region the overriding continental plate is dragged down by the subducting Pacific plate and accumulates strain with time2. Sometimes, deformation is rapidly released as co-seismic slip generated by large earthquakes3, which have struck the northeast Honshu island also in the recent past, reaching magnitudes as large as M8.2 in the last century4. These events originated in the deepest part of the seismogenic zone along the descending slab (at 30–50 km depth)4,5, as in the case of the 1978, M7.6, Miyagi-Oki earthquake6.

A magnitude 9.0 event hence sounded as almost unexpected for this zone of the Japan slab7. Several studies have inferred the amount of drag rate as compared to the velocity of the subducting plate (the interplate coupling) along the Japan trench, which provides an indication for the continuous accumulating strain along the slab2,4,8. Results have shown a large slip deficit in the northern part of the Japan arc off the Sendai coast and a large scale heterogeneity with a weaker coupling in the southern part. Such differences also appear in the pattern of interplate seismicity across northeastern (NE) Japan, with the largest historical M8–8.2 earthquakes mostly occurring in the Kuril trench and in the northern part of the Japan trench, north of the Tohoku fault surface.

This different seismic behavior and the mechanical coupling between the northern and the southern parts of the NE Japan slab can be related to the morphology of the slab, whose segmentation may be controlled by the presence of oceanic fracture zones5. In particular, the morphology of the subducting slab is observed to change around 38°N, in correspondence of a landward extending, oceanic fracture zone, which controls the strike and dip angle variation between the southern and the northern parts of the slab (line A in Figure 1a)5.

Figure 1
figure 1

Earthquake source region, stations and slip distribution.

(a) Map of the study region showing the event epicenter from the Japan Meteorological Agency, the main source regions S1, S2 and S3 discussed here, used accelerometer stations (triangles) and the depth to the subducting plate (dashed contours; from4). The red line outlines the fault plane for imaging and the grey lines labeled A and B are fractures dividing the slab into segments8. The background image is the bathymetry with blue colours indicating deeper sea bottom. (b–c) Colour-coded final slip distribution obtained from backprojection in different frequency bands. Frequencies and the respective colour scale are given below each image.

Full size image

Several kinematic models obtained from inversion of geodetic, strong-motion, teleseismic and tsunami data revealed a complex rupture propagation, with a frequency dependent radiation9. While tsunami and geodetic10 as well as low frequency local and teleseismic data indicate an anomalously large slip at shallow depths in the trench vicinity11,12, high frequency radiation around 1 Hz was inferred to be generated in the deeper regions of the plane, beneath the Japanese coast13,14. We here investigated the behavior of the rupture in the intermediate frequency range (0.05–0.4 Hz), which was proven to show some complications when kinematic models are obtained from backprojection of teleseismic data15,16. We apply an improved backprojection technique to local strong-motion data, to reveal an accurate location of radiators and estimate the associated slip amount in this frequency range.

Here we analyzed strong motion records from 246 stations of the K-Net and Kik-Net networks17 located within 450 km from the earthquake epicenter (Figure 1). This dataset provides unprecedented observations of the high frequency radiation in the near-source distance range of a megathrust earthquake. The closest station, MYG011, where a peak ground acceleration of about 1 g has been recorded, clearly displays two main high-frequency transients11,18. These are delayed about 50 seconds from each other and the second wave packet has the same peak amplitude but a larger duration than the first one (Figure 2a). The relatively long time separation and the similar peak amplitudes on near-field accelerograms suggest that the second phase is due to a radiating seismic source rather than to secondary or multi-path arrivals related to wave propagation. The band-pass filtered displacement records confirm the evidence for two main phase arrivals, hereafter referred to as S1 and S2. These two large-amplitude phases can be clearly identified in the whole distance range on the time versus distance seismic section of displacement records (Figure 2b), showing that both phases propagated at a velocity close to the average shear wave speed in the crust (3.7 km/s). Therefore and because of their linear polarization, the phases can be interpreted as S-wave pulses radiated during the initial stage of the Tohoku earthquake rupture from two distinct sources.

Figure 2
figure 2

Waveform data.

(a) Recorded acceleration waveforms at station MYG011 and vertical (UD) component of bandpass filtered displacement record (bottom trace). Green and red bars show theoretical P, S1 and S2 wave arrival times. (b) Bandpass filtered and amplitude-normalized vertical-component displacement records, sorted by distance from the JMA epicenter. The red lines approximate the arrivals of the phases S1 and S2 and the slopes of the red lines indicate an S wave velocity of 3.7 km/s. Time is in seconds after 5:46 UTC on 11 March 2011.

Full size image

With the aim of determining the space-time variation of the earthquake fault slip, we developed and applied a new beamforming and stacking technique to backproject the direct S-wave amplitudes recorded at the Japanese strong-motion networks into the offshore source region (see Methods section and Figure 1). Similar approaches, albeit limited to the identification of strongly radiating spots on the fault, have been previously applied to track the rupture of the giant M9.3 Sumatra earthquake using teleseismic data19, to image the rupture process of local or regional earthquakes with strong-motion records20,21,22 and to locate seismic events by P-wave and S-wave amplitude backprojection23.

Our slip imaging approach combines a spatial weighted estimate on where seismic energy is radiated with proper scaling of the recorded displacement amplitudes to recover the slip rate amplitude at the source as a function of space and time. In the first step, the earthquake source region is subdivided into a grid of small subsources and for any given subsource the delay times are computed to all seismic stations relative to a reference station. Then, the recorded amplitudes at all stations, within a time window bracketing the theoretical times, are summed using a weighted stacking to enhance coherent arrivals and the obtained coherency value is assigned to the current subsource. Repeating this process for all subsources provides the backprojected coherency distribution in the earthquake source region at a given arrival time. Larger values of this function occur in fault regions that are more likely to radiate seismic energy, because the associated wavelet is effectively retrieved at the right time for a large number of stations. After a normalization such that the sum of coherencies of all subsources equals one, this function will be used as a spatial weighting function, assigning a larger value to fault points where most likely high-radiating asperities were located. However, to retrieve the slip rate, the displacement amplitudes need to be corrected for attenuation and radiation pattern coefficients. We therefore probe each subsource and correct the observed amplitudes at all stations before summation by the relative distance and by the radiation pattern coefficient for a constant focal mechanism, here neglecting some anelastic attenuation. This time the stacking weights are avoided, which would otherwise distort the recorded waveforms. This corrected backprojection stack amplitude is a proxy for the slip rate at each subsource location along the fault. The integration over time leads to the estimate of the slip that would have been released at the given location and in the considered time window, if the slip alone was responsible for the observed amplitude averaged among all stations. The final slip map for a given time interval is the product of the latter, amplitude-corrected stack image with the spatial weighting function derived in the first step. The space-time evolution of the slip function can be tracked by repeating the imaging process in a sliding time window along the reference record and correcting by the respective travel times from each subsource (see Methods section).

To track the rupture process of the Tohoku earthquake, we subdivided the subduction plate boundary into a 10 km by 10 km grid of subfaults (Figure 1a) and applied our backprojection imaging to the continuous displacement records filtered in three frequency bands (0.05–0.1 Hz, 0.1–0.2 Hz and 0.2–0.4 Hz). The lowest frequency of 0.05 Hz avoids long-period drifts and ensures that the closest stations are at least one wavelength away from large-amplitude sources, while the upper limit is imposed by a considerable decrease of waveform coherency across the entire array at higher frequencies. Here we used a 20-second analysis time window, covering a full period of the dominant phases and we moved the window further in 5-second increments for the subsequent images. We ran the tracking process from the arrival of the first S-wave onset at the reference station MYG011 until the time when the recorded displacement amplitude becomes comparable to the pre-event level.

Results

Figures 1b, 1c and 1d show the cumulative slip of the entire rupture for the three investigated frequency bands. Most of the slip is released in the 0.05–0.1 Hz band, in which the slip distribution is dominated by three main large-amplitude source regions. The first one is located southwest and closest to the hypocenter, exhibits slip amplitudes of up to 25 m and is responsible for the S1 phase (Figure 2). The second, largest source corresponds to the S2 phase and is located in northeastern, up-dip direction from the hypocenter, nearby the oceanic trench. In the studied frequency band, the S2 source is imaged in a rather small area of about 100 km times 50 km and it shows maximum slip values as large as 50 m, which presumably contributed to the generation of extremely high amplitude tsunami waves. The third strong source, with a slip amount of about 40 m, is located some 200 km southwest, on the deeper portion of the slab and close to the coast. The seismic phase associated with this source region is not as clearly visible as those of the previous two, because it is hidden by the dominant S2 phase for most of the stations along the section (Figure 2b)24.

Going to higher frequency bands, a very similar image is obtained in the intermediate band 0.1–0.2 Hz, whereas the image from the highest-frequency data (0.2–0.4 Hz) shows a more chaotic distribution, with the largest values in the northern half of the rupture area and located in the vicinity of the source regions of largest slip seen in the lower-frequency images. Due to the loss of waveform coherency, higher frequency bands are not interpreted here (see Supplementary Information for data examples). From one frequency band to the consecutive, higher one, our results show a decrease in maximum slip amount in one order of magnitude (50 m, 6.5 m and 0.5 m, respectively). The total moment released on the entire fault plane within each frequency band corresponds to respective moment magnitudes Mw of 8.8, 8.3 and 7.8. In the analyzed frequency range from 0.05 Hz to 0.4 Hz we do not observe a net transition of large-amplitude slip regions toward the deeper portions of the subducting plate as the frequency increases. We still found that within all analyzed frequency bands most of the large-amplitude slip is concentrated in the shallow part of the fault close to the trench, although the amount of radiation of the deeper portion of the fault relatively increases with frequency as compared to the shallow region of the plate15.

In Figure 3a we show the evolution of the rupture as a function of time for the 0.05–0.1 Hz frequency band (see Supplemental Material for images in the higher-frequency bands). The snapshot images reveal that during the first 60 seconds the earthquake rupture developed almost along-dip within a 100 km wide slab stripe. We can observe a rather slow (less than 2 km/s) bilateral rupture in east-southeast and west-northwest directions away from the hypocenter, leading to a nearly simultaneous activation of the source regions of the two phases S1 and S2. The S1 source is activated down-dip the hypocenter, extends in a small stripe of almost 100 km downward toward the coast and lasts for almost 60–75 s, slightly deepening as time goes on. In the same time window, the propagation is dominantly upward along the slab interface, toward the oceanic trench. The phase S2 seems to be associated with a distinct source, produced by the delayed rupture of an about 15 km deep, strong asperity, generating large slip amplitudes at shallow depths. Also S2 source deepens between 45 and 90 s, breaking an adjacent stripe northward the hypocenter.

Figure 3
figure 3

Rupture evolution and moment rate.

(a) Time slices of the rupture evolution for the 0.05–0.10 Hz band in 15-second intervals. Coordinates are in km (UTM grid). (b) Moment rate functions for the three frequency bands obtained from the rupture time slices.

Full size image

The breaking of the S2 asperity near the ocean bottom has subsequently triggered the almost unilateral rupture propagation to the southwest, in a direction parallel to the trench. During this stage the rupture propagated along a distance of about 300 km from S2 at an average velocity of about 3 km/s, corresponding to more than 80% of the crustal shear wave velocity. Finally, the rupture came to arrest at the deeper part of the slab near the third strong asperity (S3) and the total rupture duration was about 160 seconds. The kinematic rupture evolution described here is generally consistent with previous results using either 0.05–0.5 Hz data from a different local seismic array25, or analyzing data from the Japanese strong-motion networks for different dominant wave periods26.

Figure 3b shows the total moment rate as a function of time, computed from the rupture evolution time slices (see Methods section). We estimated that 70% of the seismic moment has been released in a relatively small area of 150 km times 150 km near the hypocenter, covering just about 30% of the total ruptured area. Computing separate moment magnitudes for the three sources S1, S2 and S3 yields Mw 8.1, Mw 8.4 and Mw 8.2, respectively. The peak slip rates are 0.9 m/s at S1, 1.3 m/s at S2 and 1.5 m/s at S3. The rupture duration of each source is about 60 s and we found evidence of repeating slip of the source S1. These results agree with an earlier model based on teleseismic and local seismic and geodetic data12.

Discussion

We developed an earthquake rupture imaging method to obtain the location of high-frequency radiating sources and to retrieve the slip from backprojection of local displacement records. The method generalizes existing approaches and expands the capability of backprojection techniques to accurately reconstruct the slip distribution in space and time during the rupture process of large earthquakes. Additionally, it allows for a fast retrieval of kinematic properties of an earthquake rupture with high resolution images when combining records from stations covering a broad distance–azimuth range.

In this study we imaged the slip for the Tohoku event from near-source, intermediate frequency band-pass filtered, displacement records obtained after the double integration of the acceleration time series recorded by the very dense strong motion K-Net and Kik-Net arrays operating in Japan. The high-resolution images reveal three localized slip patches, which are here referred to as S1, S2 and S3.

According to our results, the giant earthquake initially grew up as twin ruptures mainly developing bilaterally along the slab-dip with a relatively low rupture velocity and being confined within a narrow, 100 km large, slab stripe. During the first 60 seconds of the rupture process, two main sources, S1 and S2, have radiated, near the hypocenter. S1 is located deeper, in the region where the slope of the slab increases, while S2 originates in the shallow part of the slab. Once the breakage of this asperity reached the trench, the rupture propagated in the southwest direction at a higher rupture speed, expanding into the deep, southern part of the NE Japan slab. A third source, S3, radiated close to the southern limit of rupture zone and at a greater depth than S1 and S2. Although S1 and S2 activated simultaneously, they show different properties. The shallower source S2 produced a larger slip (around 50 meters), reaching its maximum amplitude later than the deeper source S1. Our rupture images in different frequency bands show that the S2-to-S1 slip amplitude ratio is larger at lower frequencies that at higher frequencies (Figure 1b,c). This suggests a down-dip migration of the high-frequency component of the slip, consistent with similar observations from teleseismic analyses18. When taking S1 and S2 as single asperities, an omega-square model for their spectral decay predicts a constant amplitude ratio between the two sources as the frequency increases. Since we observe that the relative ratio decreases as a function of the frequency, we argue that the source S2 follows a faster spectral decay in this frequency range.

The low frequency radiation of the shallow asperity S2 lead to an anomalously large slip (up to 50 m) which has been ascribed to a dynamic overshoot occurring in the shallow portion of the fault9. This can be due to dynamic weakening associated with thermal pressurization27 or to a stress increase associated to waves reflected at the free surface28. Such a large slip localized in a very limited region close to the trench was also interpreted as the cause of the large wave pulse clearly observed on two the ocean-bottom pressure gauges TM1 and TM2 off Kamaishi, northern Japan29.

Differently from several studies of the Tohoku event13,14 based on teleseismic and strong motion records, we do not observe the high-frequency radiation prominently generated by deep sources (greater than 30 km depth), closest to the coast. Instead, our analysis supports the hypothesis for the largest slip asperities in the whole investigated frequency range generated along the shallow portion of the slab (smaller than 30 km depth). The relative location of sources S1 and S2 has been obtained here by the backprojection technique which provides the most likely centroid location of the associated wave packets. As compared to the standard picking procedures, backprojection does not assume or need to specify beforehand, which wave wiggle actually belongs to the phase of interest and thus reveals to be a very robust location technique30. In addition, the reliability of the relative location of S1 and S2 was validated through a number of synthetic simulations (see Supplementary Information). We therefore conclude that the data coverage allows for the precise relative location of the S1 and S2, with a possible 10-km smearing effect due to the frequency content, acquisition layout and uncertainties in the velocity model.

The whole rupture process lasted about 160 seconds, with 70% of the moment released by sources S1 and S2 during the first half of the rupture duration, when the rupture remained confined in a 100–150 km wide slab segment in the northern part of the fault. This suggests the presence of strong lateral geometrical and/or mechanical barriers. Difference in radiation between the northern and southern parts of the fault reflects also on the location of large earthquakes along the slab31. In particular, the westward prolongations of oceanic fracture zones (A and B, as denoted in ref. [5] and drawn in Figure 1) delimit a surface along the plate boundary which contains the S1 and S2 source areas. The occurrence of large slip amplitudes at shallow depths might have been a favorable condition to dramatically increase the earthquake size by allowing the rupture to jump across the barriers and propagate into adjacent plate segments. Subducted oceanic fractures may have acted as barriers during the Tohoku event as it was also observed for other megathrust earthquakes, such as the 2004 Peru earthquake32. The slab segment containing S1 and S2 is also characterized by a relatively high P wave velocity, indicating the presence of morphological and compositional changes along the slab surface as compared to the adjacent regions33. Within this zone, subducted oceanic ridges, seamounts or other topographic highs may act as localized asperities which increase the coupling between the two plates33,34,35.

Several previous studies on the Tohoku earthquake demonstrated a strong frequency-dependent source process with short and long periods generated at different locations and times [e.g.15,34]. The intermediate frequency range studied here complements the previous images of this megathrust earthquake and further supports its frequency-dependent characteristics. The location of source S1 is similar to the 1978 and 2005, Miyagi-Oki earthquakes locations6, while the 869 Jogan earthquake36 occurred in the region of S2 source area. Locations and rupture histories of these two main asperities suggest two seismic cycles at the plate boundary, where large earthquakes (M < 8) occur every few decades in the deeper part of the slab, whereas megathrust events (M ≈ 9) originate in the shallow part of slab with a returning period of about a thousand years.

Methods

We applied a baseline correction to the strong-motion acceleration records from 246 K-net and KiK-net stations (Figure 1) located within 450 km from the JMA (Japan Meteorological Agency) epicenter, integrated the time series twice to obtain ground displacement and filtered the displacement traces in three different frequency bands between 0.05 and 0.40 Hz.

The reference station for the backprojection procedure is MYG011 at 155 km distance from the JMA epicenter. To estimate the average shear wave speed of 3.7 km/s in the study area, we searched for the best phase alignment in the displacement records according to delay times computed for the known hypocenter and for various shear wave speeds. We assigned a spatial grid of subfaults with 10 km node spacing following the geometry of the subducting plate in the greater earthquake source region, as outlined in Figure 1. To obtain the space-time evolution of the slip distribution on the fault plane, we developed a new imaging technique which combines standard amplitude backprojection with correction for the distance and radiation pattern. The time window for a single backprojection image is 20 seconds and the reference time increment is 5 seconds. The window length of 20 seconds is based on the period of the dominant phases and it provides a smoother, more stable image than the use of shorter time windows of 5 or 10 seconds.

For computation of the slip rate, the displacement amplitudes appropriately corrected by distance Rij (the geometrical spreading) and radiation pattern Fij37 for a constant focal mechanism are stacked and back-projected for all possible elementary sources located along the fault surface. The resulting slip rate for the i-th grid point and at the time t is obtained as the summation

In the above expression Uj is the observed displacement at the j-th station, tiR the traveltime to the reference station, δtij the time shift between the j-th station and the reference station, Ns the number of stations, ρ = 2900 kg/m3 is the density (mean density of oceanic crust38), vs = 3.7 km/s the S-wave velocity, μ = ρ vs2 = 40 GPa the shear modulus and A the area of the subfault. In the far field condition, this quantity is related to the slip rate released in the subfault. Note that Equation 1 ignores possible anelastic attenuation. However, in the present case with Q = 18039 and the lowest frequency band which yields the largest slip, the amplitude decay due to anelastic attenuation is less than 10% in the distance range considered, which is less than scatter due to site effects. Integration of the above function over time yields the slip Si associated to the i-th subfault, under the hypothesis that such a subfault would be the only responsible for the observed amplitudes on average at all the stations.

However, several subsources could be responsible for the observed amplitudes at the same time on the seismograms. We therefore search for the regions of the fault plane which had more likely radiated the observed amplitudes. With this aim, an objective function for slip source location is defined as the standard weighted back-projected stack amplitude for the i-th grid point and the time t

which corresponds to the 1/n quasi-norm of the displacement at all stations after trace re-alignment based on delay time and cross-correlation, by applying the same procedure as for the slip rate computation.

For nth-root stacking, the nth root of the trace is taken before summation while keeping the sign of each sample and the nth power is computed after summation, again keeping the sign40. Increasing n would enhance the stack of more coherent waveforms, which is indeed a more restrictive criterion that tends to localize the function W on a smaller area of the fault plane. The normalization factor C(t) in Equation (2) ensures that the sum of W over all grid points is 1. According to its definition, the function W represents a spatial weighting function, assigning a larger value to fault points where most likely high amplitudes have been radiated.

We selected a 4th-root stack (n = 4), this choice representing a good compromise between reducing the effect of amplitude differences, enhancing coherent phases and avoiding the introduction of artificial space lengths smaller than the minimum signal wavelength. The final slip map as shown in Figure 1 is at the end computed as the product Wi Si, such a distribution still justifying the observed stacked amplitudes. It is worth to note that the choice of the norm would affect the final slip distribution, as in all kinematic inversions for slip retrieval. However, from n = 2 to n = 4 the total moment changes of about 15% and the magnitude of 0.1, this value being smaller than the error associated with standard estimates.

Each backprojection image belongs to a given arrival time at the reference station MYG011. We can compute the origin time of a given seismic phase arriving at a given reference time at MYG011 by subtracting the travel time of a wave from the identified source location (stack amplitude maximum). This allows to track the rupture propagation in time and space and also accounts for the directivity effect19. Applying this time correction to all backprojection images and all subfaults allows for the reconstruction of the temporal rupture evolution as shown in Figure 3a. Assuming proper amplitude scaling, a spatial integration of each such rupture time slice provides the moment rate at that time, which leads to the moment rate functions shown in Figure 3b. These moment rate functions appear rather smooth due to the finite time window in the backprojection procedure and a secondary smoothing caused by temporal binning during the origin time correction.

To assess the resolution and the reliability of imaged source locations, we computed synthetic S-wave seismograms for two point sources at S1 and S2, having comparable frequency content as the real signals. An application of the same backprojection procedure as for the real data shows resolution as good as one grid node (10 km) in north-south direction, while some amplitude smearing is observed in the east-west direction (see Supplementary Information). Independently from any reasonable imaging velocity, the relative locations of S1 and S2 remain similar, with S2 located closer to the trench than S1. We could also conclude that the approximation of a constant imaging velocity is justified here.