Rupture Process of the 2011 Tohoku-Oki Earthquake Based upon Joint Source Inversion of Teleseismic and GPS Data

This study investigated 18 broadband teleseismic records and 451 near field GPS coseismic deformation data to determine the spatial and temporal slip distribution of the 2011 Tohoku-Oki earthquake (M 9.0). The results show a large triangular shaped slip zone with several asperities. The largest asperity centered above the hypocenter at about 5 30 km depth. A secondary large asperity was found in the deeper subduction zone beneath the hypocenter. The average slip on the fault is ~15 m and the maximum displacement on the biggest asperity is > 30 m. The temporal rupture process shows that the slip nucleated near the hypocenter at the beginning, and then ruptured to the shallow fault plane forming the largest asperity. The slip developed in the deeper subduction zone in the second stage. Finally, the rupture propagated toward the north and south of the fault along the Japan Trench. The source time function shows three segments of energy releases with two large peaks related to the development of the asperities. The overall rupture process is ~180 seconds. This source model coincides well with the aftershock distribution and provides a first-order information on the source complexity of the earthquake which is crucial for further studies.


InTRODucTIOn
On 11 March 2011, a giant earthquake struck Japan with magnitude of M 9.0.The earthquake reported by the Japan Meteorological Agency (JMA) and US Geological Survey (USGS) both show that the hypocenter of the mainshock was located off the Pacific coast of northeast Honshu (Fig. 1).The initial report from the JMA suggested that the earthquake was a magnitude of 7.9 which was then modified to 8.9 after an hour into the event.Further analysis of the seismic data resulted in an upgrade to a magnitude of 9.0.According to the USGS W-phase moment tensor inversion result, faulting during the earthquake has been interpreted as reverse faulting on a low-angle plane along the Japan Trench.
The studies of the rupture process of the Tohoku-Oki earthquake have been carried out intensively (e.g., Ide et al. 2011;Ozawa et al. 2011;Simons et al. 2011).However, in most of the studies only one kind of data or few data sets were used to invert the source model, thus the results might not be able to explain other coseismic observations.For example, using only GPS coseismic displacement data (i.e., Ozawa et al. 2011) cannot provide temporal rupture process information in the inversion and might not be able to explain the waveform characteristics observed at regional seismic stations.In order to obtain first-order information concerning the origin of this earthquake, especially how the rupture develops into the megathrust event, this study performs a joint source inversion using teleseismic body wave and GPS coseismic deformation data.This large earthquake was well recorded by the teleseismic stations worldwide providing a good opportunity to arrive at a quick determination of fault rupture behavior.The teleseismic data have good data quality and azimuthal coverage to the source.Furthermore, the near field GPS coseismic deformation data also provide a good constraint on the total slip pattern of this giant event.

Data
The 18 teleseismic stations used in this study are shown in Fig. 2 and their station parameters are listed in Table 1.The stations are chosen for the epicentral distances between 30° to 100° to avoid the complex earth structure.The data are digital recordings obtained from the Incorporated Research Institutions for Seismology (IRIS), and the stations chosen for the study provide a good azimuthal coverage from the source.Instrument responses were removed from the original recordings.The P wave arrivals were picked in the raw data first, and then retained 10 seconds before and 190 seconds after the P arrivals.A band-pass filter was applied on the displacement waveforms between 0.01 to 0.2 Hz.The sampling rate used in inversion is 2 points per second.For such a large event which usually slips over a large area, an earthquake cannot be treated as a point source, so it is not easy to rotate the original E and N component data into the R and T component to obtain a precise SH waveform.Thus, this study used the P wave only as shown in Fig. 2, the SH wave has not been considered in the inversion.
The GPS coseismic displacements of the event are shown in Fig. 3.This data set was computed from the difference between the averages of 2 day GPS site positions before and after the main shock; notably, this data set might be contaminated by large aftershocks or post seismic deformation.In total, 451 sites of three component GPS coseismic displacements were compiled by the Taiwan Earthquake Research Center (TEC) GPS Lab (http://gps.earth.sinica.edu.tw).All original GEONET RINEX data was provided to the TEC GPS Lab by the Geospatial Information Authority (GSI) of Japan.The GPS data shows that most of the coseismic deformation was mainly toward the east of Honshu.The large horizontal displacements are moved eastward and focused toward the offshore area near the mainshock epicenter with a maximum value of about 4.2 m; vertical displacements are relatively small which show a maximum downward motion of about 0.5 m along eastern coast and mostly uplift motion (less than 0.2 m) in western Honshu.

Inversion Method
Finite fault inversion problems are usually formulated in a linear form, Ax = b where A is a Green's functions matrix, b is the observed data vector and x is the solution vector of slip on each subfault.A misfit function, defined as (Ax -b) 2 / b 2 , is used to evaluate the quality of a solution.An important approach to the inversion problem is the introduction of multiple-time windows resulting in a better spatial and temporal resolution of slip (Hartzell and Heaton 1983).However, an increase in the number of time windows would lead to a large expansion of matrix A making the solution of this problem very costly in terms of computer time.Lee et al (2006) develop a parallel Non-Negative Least-Squares (Parallel NNLS) inversion technique which decomposes matrix A into different computing nodes and solves for the vector x for each time window, and thus promoting program performance.Further information about this technique and its application to finite fault inversion for the rupture process of large earthquakes can be found in Lee et al. (2006Lee et al. ( , 2008) ) and Konstantinou et al. (2009).

Inversion Settings
The fault plane solution determined from the USGS W-phase moment tensor inversion (strike 193°, dip 14° and rake 81°) is considered in this study.The hypocenter on the fault is given by the JMA earthquake report (142.8°E,38.05°N, 24 km in depth).With the combination of the JMA earthquake location and USGS W-phase focal mechanism, the shallowest part of the fault plane can fit better with the surface features of the Japan Trench.The fault plane is composed by 32 × 12 subfaults each having a dimension of 20 × 20 km 2 covering a sufficiently broad area of length and down-dip width of 640 and 240 km respectively (see Fig. 1).
In matrix A, the teleseismic Green's functions are calculated by using the approach developed by Kikuchi and Kanamori (1982).The near source structure and the receiver side structure both use the global 1D PREM model (Dziewonski and Anderson 1981).The source and receiver functions are computed using Haskell's propagator matrix following Bouchon (1976) and Haskell (1960).Once Green's functions are computed, the seismic synthetic waveforms are filtered between 0.01 -0.2 Hz in the same way as the observed data in teleseismic Green's functions.For the geodetic Green's functions, the analytic expression provided by Okada (1985) is used to calculate the horizontal and vertical static displacement for surface deformation resulting from a uniform slip over each subfault.The inversion is based on the full time-space inversion approach using the Parallel NNLS technique (Lee et al. 2006).There are 60 time windows set in the inversion wherein each subfault is allowed to slip in any of the 4 s time windows following the passage of the rupture front while each window has an overlap of 2 s.Thus, each subfault can slip by flexible source time function (by any combination of 4 s time window) within a time period of 122 seconds after the rupture front passes through.A larger initial rupture velocity of 3.5 km s -1 is given in the inversion in order to capture the possible supershear rupture behavior.By the use of 60 time windows, the variation in the rupture velocity can vary from zero to over shear wave velocity.The final rupture velocity at different part of the fault plane is determined based on the inversion result.Finally, some stability constraints are imposed such as minimization of the seismic moment to the value derived in this study, damping at the edge of the parameterized fault (except the shallowest side) and smoothing on the slip between adjacent subfaults in order to avoid unrealistic slip distributions.In order to make the two data sets and inversion constraints contribute equally in the inversion, a normalized weight for each data set is given: Observed data value of the teleseismic waveform is in the unit of micro-meter and the GPS coseismic displacement is in the unit of centimeter.The number of data for the teleseismic waveform is 7200 and that for the GPS data is 1353.By following this rule, the normalized weightings for teleseismic and GPS data are 0.015 and 0.5 respectively.This study considers 60 time windows in the inversion; the temporal resolution can reach a high resolution in every 2 second during the rupture.The inversions require 60 CPU cores and 39.2 gigabytes of memory which are carried out on IES HPC cluster.The average wall-clock time in one inversion is 4.5 hours.

Spatial Slip Distribution
The slip distribution determined from teleseismic and GPS joint source inversion is shown in Fig. 1.Inversion results show that the rupture occurred on a large triangular shaped slip zone.It has an apex located in the deep western part of the fault plane near Sendai.The slip was concentrated in several areas, causing one large asperity and several secondary asperities on the fault.The largest asperity developed near the epicenter between 5 -30 km in depth (Asperity I) and a second asperity located in the deeper subduction zone beneath the hypocenter at about 45 km depth (Asperity II).The other two asperities occurred in the north and south in relation to the hypocenter and both were centered at about 10 km in depth (Asperity III and IV).The slips of these two asperities overlapped with Asperity I and were recognized as they developed during the temporal rupture process which is discussed in the next section.The slip of the largest asperity was > 30 m with a predominately reverse motion and covered a large area of about 200 × 120 km 2 .The slips of asperities II, III and IV were smaller which are about 10 -20 m.Rupture behaviors in these areas were a little different; they were mainly right-lateral, left-lateral and predominately thrust slip respectively.A close relationship between the fault slip projection onto the ground surface and the aftershocks distribution can be seen in Fig. 1.Most of the aftershocks occurred around the large slip areas, and no big aftershocks (M F 6) were located on these asperities.
The comparisons between inverted teleseismic synthetic waveform and observed data are shown in Fig. 2. Peak amplitudes and later phases can fit sufficiently well, especially for the stations on the northwest side where waveforms are complex and have larger amplitudes.Conversely, stations located in the southeast display smaller amplitude waveforms, but their characteristics are generally comparable.The teleseismic waveform misfit in the joint inversion is 0.18.Synthetic and observed GPS coseismic displacements are shown in Fig. 3.The large horizontal displacements that focus to the epicenter are predicted from this model and its amplitude as well as azimuth are recovered by the inversion.For the vertical displacement, the fits are also good in most areas; however, the fits are poor for a few GPS sites located in northern Honshu.In general, the synthetic GPS coseismic displacement can explain the pattern of observed GPS displacement.The GPS data misfit in the joint inversion result is 0.03.The total misfit of the joint inversion is 0.14.

Temporal Rupture Process
The snapshots of the rupture process are shown in Fig. 4a.Large apparent slips (about 5 m) were observed the first twenty seconds at the middle part of the fault plane close to the hypocenter.Although they were relatively small compared to the final rupture, these slips might indicate the nucleation of the earthquake.At about 40 second, a large thrust movement quickly occurred above the hypocenter.This area continued to slip in the shallow part for more than 40 seconds and merged with the initial slip near the epicenter, forming an extremely large slip area and becoming the first formed asperity on the fault plane (Asperity I).After about 80 second, a rupture developed in the deeper part of the fault and formed the second asperity (Asperity II).Although these two asperities centered at different depths, their slips were largely over lapped which formed the main slip zone offshore of the Tohoku.After about 100 second, the shallow part of the fault seemed to have significant movement again.Then the rupture extended bilaterally to both north and south along the Japan Trench, forming Asperity III in the north and then at about 150 second forming the last asperity (Asperity IV) in the south.
Figure 4b is the inverted moment rate function.The moment release during the earthquake is very complex; compared to the rupture snapshot, three energy release segments can be found.The first lasted for the first 40 seconds (T1) with a relatively small peak amplitude at ~30 second, then the second release (T2) ran from 40 to 95 second with a peak value at approximately 70 second.It looks to be quieting between 95 -100 seconds, and then the third moment release occurred between 100 to 180 seconds (T3).From the temporal slip inversion result, the Tohoku-Oki earthquake appears to have been caused by a chain of slip events that are related to the development of asperities.By comparing the moment rate function with the rupture snapshots, we can infer that the time period T1 was related to the occurrence of the rupture nucleation.The second energy release time (T2) included the growth of the biggest asperity (Asperity I) at the shallow part above the hypocenter.Almost at the same period, Asperity II formed at the deeper part of the fault beneath the hypocenter.The slip on the fault plane was quiet at ~100 second, and then the shallow portion of the fault seemed to slip again and extend along the Japan Trench in the time period T3.The overall rupture duration time was about 180 seconds for a total seismic moment of 0.345 × 10 30 dyne-cm, which is equivalent to an earthquake of M w 8.96.
The variation of the rupture speed can be approximately obtained by lining the rupture front (a 5-m contour of accumulated slip) of the snapshots as shown in Fig. 4a.At the beginning, the rupture velocity near the hypocenter was initially ~1.7 km s -1 , smaller to common observation of  about 85% of the shear wave velocity.But it increased to ~2.0 km s -1 at the shallow part (Asperity I) between about 40 -80 seconds.After that, the rupture velocity rapidly decreased to 0.0 -0.5 km s -1 .This could be due to the fact that the rupture developed in the deeper part of the fault below the hypocenter between 80 to 100 seconds.After 140 second, the rupture front began to move in the shallow subduction zone along the Japan Trench with a rupture velocity of approximately 1.0 -2.0 km s -1 .
This giant earthquake was well recorded by dense regional seismic networks (e.g., K-NET, KiK-net and F-net) which can provide a further constraint on the details in the temporal rupture process and rupture velocity variations.A joint source inversion analysis with the constraint of regional waveform data is ongoing.

cOncLuSIOn
This joint source inversion analysis uses teleseismic body wave with constraint of near field GPS coseismic deformation data.The results provide first-order information of the rupture process of the 2011 Tohoku-Oki earthquake.The spatial slip distribution reveals a large triangular shaped slip zone off the Pacific coast of Tohoku.A large anomalous asperity with a dimension of about 200 × 120 km 2 centered at 5 -30 km in depth above the hypocenter is observed.The slip in this area is over 30 m.A large secondary asperity is found in the deeper subduction zone centered at depth of 45 km beneath the hypocenter.The spatial slip distribution pattern coincides well with the aftershock distribution.The temporal slip process shows that the rupture nucleated near the hypocenter at the beginning, and then ruptured to the shallow part forming an anomalously large asperity.The slip developed in the deep subduction zone in the second stage.Finally, the rupture propagated to the shallow fault plane along the Japan Trench.Three time periods of energy releases relating to the development of asperities are found.The overall rupture duration time is about 180 seconds for a total seismic moment of 0.345 × 10 30 dyne-cm that is equivalent to an earthquake of M w 8.96.Such a giant earthquake will surely be the focus of in depth study for the foreseeable future.

Fig. 1 .
Fig. 1.Map view of the spatial slip distribution of the Tohoku-Oki earthquake.The slip values are shown as color scale indicated at the bottom.The vectors indicate slip directions.The red star shows the hypocenter reported by the JMA.Aftershocks occurred within 2 months after the main shock are shown by solid circles.The beach ball shows the focal mechanism derived from USGS W-phase moment tensor inversion.

Fig. 2 .
Fig. 2. Comparison between teleseismic observations (black lines) and synthetic waveforms (red lines).All the waveforms are vertical component in displacement type.The maximum observed ground displacement is shown at the end of each waveform.

Fig. 4 .
Fig. 4. (a) Snapshots of the rupture process.The open blue squares show the asperities that occurred at different times during the rupture.The black lines that link the rupture front (a 5-m slip contour) between each snapshot are used to identify the variation of the rupture speed.Red numbers show the estimated rupture speed value.(b) The moment rate function.Three segments of moment release time periods are separated by dotted lines.They are labeled by T1, T2 and T3 respectively.