Slip model for the 2011 Mw 9.0 Sendai (Japan) earthquake and its Mw 7.9 aftershock derived from GPS data

Using GPS-measured coseismic and post-seismic displacements for the 8 h following the Mw 9.0 Sendai earthquake of March 11, 2011, coseismic and post-seismic fault slip models were developed based on a layered crustal model. The geodetic moment magnitude of the main shock was measured as approximately Mw 8.98. The slip exhibits clear reverse characteristics, with a maximum near the hypocenter, and a magnitude of about 23.3 m. Some strike-slip behavior may occur on the two sides of the peak rupture zone. Almost 90% of the seismic moments released by the main shock occurred at depths less than 40 km. The energy released by the fault slip in the 8 h following the main shock is approximately equal to an earthquake of Mw 8.13. With a maximum of ∼1.5 m, the post-seismic slip was concentrated in the southwestern part of the coseismic rupture fault, which agrees well with the location and behavior of the Mw 7.9 aftershock. This implies that the post-seismic deformation in the 8 h after the main shock was mainly induced by the Mw 7.9 aftershock. In addition, a post-seismic slip of 0.2–0.4 m was observed at the down-dip extension of the coseismic rupture, which may have been caused by the effect of after-slip during this period.

At UTC 05 : 46 on March 11, 2011, a M w 9.0 earthquake ruptured along the Japan Trench, 130 km east of Sendai off the northeastern coast of Honshu. About 30 min later, a M w 7.9 aftershock occurred 250 km southwest of the epicenter of the main shock. The earthquakes and subsequent consequences, including tsunamis and the nuclear crisis at Fukushima, have been disastrous for residents.
Along the northeastern part of Japan, the Pacific Plate subducts beneath the North American Plate at a velocity of ~9 cm/a [1]. Many large earthquakes have occurred along this plate boundary, including, most recently, the M w 9.0 Sendai earthquake and its M w 7.9 aftershock of March 2011 ( Figure 1). In contrast with continental earthquakes, the structural genesis of earthquakes occurring within subduction zones is well understood. During the interseismic period, the Pacific Plate subducts westward at a stable rate, although some asperities on the plate interface may be "locked" during this period resulting in the accumulation of energy. When the accumulated energy exceeds the strength of the rock, transient rupture and earthquake occur. Based on ground observations between 1996 and 2000, Hashimoto et al. [2] investigated the slip deficit along the Kuril-Japan trench and found that the plate interface near Sendai was almost completely "locked" with a slip-deficit rate of ~9 cm/a.
After the Sendai earthquake, a number of rupture models for this earthquake were constructed using teleseismic waveform data. However, significant differences are apparent between the various results (http://supersites.earthobservations.org/sendai.phd). Near-field geodetic observations provide additional constraints on slip models [3]. An accurate slip model is critically important to constrain estimates of the Coulomb stress change on surrounding faults, and to assess resultant seismic hazards and analyze possible postseismic mechanisms. In this paper, we employ a constrained Figure 1 Tectonic setting of the M w 9.0 Sendai earthquake. The two red beachballs show focal mechanisms and epicenters of the main shock and the M w 7.9 aftershock (http://earthquake.usgs.gov/earthquakes/seqs/events/usc0001xgp/). The purple rectangle is the surface projection of the fault rupture plane. The green stars and gray ellipses indicate epicenters and source regions, respectively, for large interplate earthquakes (M w ≥7.5) that occurred in the past century [2]. The white dots are GPS stations. In the top left panel: AM, the Amurian Plate; PH, the Philippine Plate; PA, the Pacific Plate; NA, the North American Plate. least squares method based on a layered earth model using near-field GPS observations. This enables the derivation of a coseismic slip distribution for the M w 9.0 Sendai earthquake and the post-seismic slip distribution in the 8 h following the main shock.

Data
After the earthquake, the ARIA group at Jet Propulsion Laboratory (JPL) (California, American) and California Institute of Technology (California, American) processed the original GEONET RINEX GPS data provided by the Geospatial Information Authority of Japan (Tokyo, Japan) (ftp://sideshow.jpl.nasa.gov/pub/usrs/ARIA). The coseismic displacements and crustal displacements in the 8 h following the main shock provided an important database for slip model inversion. GIPSY-OASIS software (developed by JPL) and JPL's rapid orbit data were used in the data processing; the coordinates of each station were calculated every 5 min. The coseismic displacements of the main shock were obtained using the differences in station coordinates between UTC 5 : 40 and UTC 5 : 55, whereas the post-seismic displacements 8 h after the main shock were calculated using the station coordinates at UTC 5 : 55 and UTC 14 : 00. Because the data times used in the calculation of coseismic displacements is very close to the earthquake rupture time, the estimated coseismic displacements are mainly induced by the main shock without the effects of after-slip and aftershocks. However, the post-seismic displacements in the 8 h following the main event are complicated because this information may include contributions from both strong aftershocks and aseismic after-slip.
As shown in Figure 2, with a maximum horizontal displacement of 5 m, clear eastward crustal movements of northeastern Japan were observed after the M w 9.0 earthquake. In addition, significant vertical displacements were also found near the coastline. The largest horizontal displacement was close to 0.5 m in the 8 h after the main shock. However, compared with coseismic displacements, the area with the most significant post-seismic displacements is closest to the M w 7.9 aftershock, ~250 km southwest of the main shock ( Figure 3). We used both horizontal and vertical displacements in the coseismic slip model inversion, but only horizontal displacements were used in the post-seismic slip  model inversion due to the poor signal-to-noise ratio of vertical displacements.

Methods and parameterization
The method proposed by Wang et al. [4] was used in slip model inversions of the Sendai earthquake. When the fault geometry is fixed, the relationship between surface observations and the fault slip can be expressed by where G is the Green's function, which here is a function of the parameterization of fault geometry, such as strike, dip, length, depth and location of the fault. The Green's function can be calculated using dislocation theory for an elastic halfspace or a layered earth model; b is the slip on sub-faults, including both strike and dip components; y denotes the ground observation of GPS displacements; and ε donates the error, which includes observational error and theoretical errors induced by model simplification.
To obtain a slip model with high resolution, the fault plane is generally discretized into many uniform rectangular sub-faults. However, when the number of unknown slip parameters on the sub-faults is larger than that of the number of observations or the equations show a strong correlation, the solution becomes unstable. To obtain a stable and reliable result, a priori constraints are introduced, including smoothing the slip and stress on the fault plane. Thus, solving eq. (1) becomes a constrained least square problem, and the objective function can be expressed as [4]: where τ is the shear stress drop caused by the distributed slip on the whole fault plane, H is the finite difference approximation of the Laplacian operator, and α is the smoothing factor, which controls the trade-off between model roughness and data misfit. Before running the inversion, we first determine the geometric parameters of the fault rupture plane. Because of the absence of surface rupture information, we fix the location using the aftershock distribution, hypocenter, and slab contours [5]. The north end of the fault was fixed at 40°N based on the aftershock distribution in the initial stage. The length and width of the fault plane were set to be 500 and 200 km, respectively. In addition, when inverting the fault slip for the 8 h after the main shock, the southwest end of the fault was prolonged by 100 km from the aftershock distribution.
We selected a mean value and fixed the strike at 200°, although strike angles of the ruptured fault inferred from teleseismic inversions show slight differences with this value. The dip is inferred from seismic studies to increase from the trench towards the land [6]. The smoothing factors were estimated from the trade-off curves between model roughness and data misfit (Figure 4). The layered crustal model from CRUST 2.0 [7] was used to calculate the Green's function. Note that we used a layered crustal model on the hanging wall of the fault (Table 1), as all stations are on the hanging wall.

Results and discussion
Based on surface GPS observations, the layered earth model, and the fault geometry discussed in section 2, we inverted for the coseismic slip and the slip model for the 8 h following the main shock. The results are shown in Figure 5. The geodetic moment magnitude of the earthquake is M w 8.98, which is close to the M w 9.0 inferred from long period W Phase and the M w 9.1 from the Global CMT. This demonstrates that the magnitude of the Sendai earthquake is indeed around M w 9.0, and argues that this earthquake is one of the biggest instrumentally detected earthquakes in history.
The earthquake is dominantly pure thrust slip with a maximum slip of 23.3 m at a depth of 25 km near the epicenter. About 90% of the energy is radiated at depths less than 40 km, as a result of the relatively small dip angle of the fault. In addition, some strike-slip slips can also be found along the two sides of the slip concentration. The movement of the fault is, however, dominated by thrust slips, which are even obvious in areas with large slips.
We let the rake angle vary freely from 0° to 180° during the inversion process, and no other constraints were added to the rake angle. The appearance of small strike-slip components may be caused by two factors. (1) The subduction of the Pacific Plate may not be completely perpendicular to the strike of the fault and some strike-slip strain may accumulate at the plate interface. (2) Strike-slip lateral movement may have resulted in response to blockage of the subducting plate. Other studies have also observed some strike-slip components in their slip models (e.g. http://earthquake.usgs. gov/earthquakes/eqinthenews/2011/usc0001xgp/neic_c0001 xgp_cmt.php, http://tectonics.caltech. edu/slip_history/2011_ taiheiyo-oki/index.html, http://www. geol.tsukuba.ac.jp/~yagiy/EQ/Tohoku/). Thus, the existence of the strike-slip component of motion and its proportion require further investigation.
Our slip model fits the observed data well (Figures 2 and  3). For the main shock, the root mean square residuals in the N-S, W-E and vertical directions are 2.5, 3.7 and 4.8 cm, respectively. In the slip model for the 8 hours after the main shock, the root mean square residuals are 1.4 and 1.3 cm in the N-S and W-E directions. In general, the root mean square residual is close to or better than the observational error.
The fault slip accumulated during the 8 h after the main shock is shown in Figure 5(b). The total energy released by the fault slip in this period is equivalent to a M w 8.13 earthquake. Based on our model, the maximum slip is up to 1.5 m and shows a pure thrust-slip characteristic. Although the slip concentration is located mainly southwest of the main slip area, slips of 0.2-0.4 m can also be found at the downdip extension of the main rupture ( Figure 5(b)).
Besides the fault slip, the post-seismic viscoelastic relaxation of the lower crust and the upper mantle can also induce significant post-seismic deformation [8]. To investigate the possibility of a viscoelastic effect, we applied a simple layered model in which both the lower crust and the upper mantle are set to be viscoelastic media. The coseismic  slip model inferred from this work was applied as the coseismic stress loading. The code PSGRN/PSCMP [9] was used to calculate the post-seismic displacement at each GPS station that was induced by 8 h of viscoelastic relaxation.
Numerical results show that post-seismic displacements generated by viscoelastic relaxation in this time period are less than 1 mm. Thus, essentially all of these post-seismic deformations are induced by fault slip and not viscoelastic relaxations.
Although the effect of viscoelastic relaxation is negligible, two other mechanisms may be responsible for observed post-seismic deformation: strong aftershocks and aseismic after-slip on the fault. It is normally difficult to distinguish the proportion of the two mechanisms. However, the aftershocks that occurred during the 8 h after the main shock provide the opportunity to analyze each contribution. Figure 6(b) shows the spatial distribution of aftershocks  and the slip distribution. During the initial 8 h, only two aftershocks with magnitudes larger than 7 occurred. One was located 200 km east of the main shock with the focal mechanism of a normal earthquake. To analyze the effect on the fault slip inversion caused by this aftershock, we first built a unified slip model according to the focal mechanism and magnitude of this M w 7.6 aftershock to forward model the surface displacements. Then we removed the effect of the aftershock from the observed post-seismic displacement, and repeated the inversion for the fault slip. The result indicates that this aftershock has only a very small role in the post-seismic deformation field. The induced maximum discrepancy is only 0.02 m between the corrected slip model and the original one, far less than the slip value (~1.5 m). The other M w 7.9 aftershock, which occurred 250 km southwest of the main shock, shows a strong correlation with the inverted slip concentration both in the spatial location and focal mechanism. Moreover, the moment released by the aftershock is close to that released by the slip concentration near the rupture of the aftershock. Based on these factors, we deduce that the post-seismic dislocation in the southwest part of the main shock is mainly caused by the M w 7.9 aftershock, whereas the slips at the down-dip extension of the coseismic rupture are mostly due to aseismic after-slip.
There is a good complementary relationship between the aftershock, after-slip and coseismic slip distributions. The after-slip and aftershocks mainly occurred on the surrounding areas of the major rupture zone of the main shock. This is similar to the relationship between the coseismic slip and aftershock distributions of the Morgan Hill earthquake observed by Bakun et al. [10]. A possible reason for this phenomenon is that the coseismic slip increases the stress field in the surrounding areas of the main rupture zone, and aftershocks and the after-slip are consequences of the gradual release of the accumulated stress field [11].
Frequent aftershocks produce interference between seismic waves, which makes it difficult to invert for the rupture process using seismological observations. However, although static geodetic observations cannot distinguish contributions of aftershocks and after-slips, they do provide independent static slip models, which gives important evidence to illuminate the coseismic and post-seismic slip distributions. With further analysis of displacement time-series records, if deformation caused by the strong aftershocks can be isolated, we can obtain the slip distribution of strong aftershocks as well as the after-slip.
The simplification of our fault model is noteworthy. In this model, only a unified 500 km straight fault plane is adopted and the local effects of a curved fault are neglected. To assess the effect of the simplification on the inversion, we constructed a curved fault plane roughly based on the trend of the trench. Following the same procedures for the inversion mentioned above, we obtained slip models for both of the coseismic and post-seismic 8 h periods. The results of the unified fault and the curved fault models are quite similar (Figure 7); the only significant difference is that the slip value of the curved model is slightly larger than that inferred from the unified model. In the curved fault model, the maximum of the coseismic slip and that in the 8 h following the main shock, are ~26 and ~1.7 m, respectively, whereas in the straight model they are ~23.3 and ~1.5 m. Thus, the difference between the two models is about 10%. Meanwhile, the similar slip values in the two fault models and the slip pattern show a high correlation and the lateral slip component is also observed in the curved model.
Because of the larger distances (> 200 km) between observation stations and fault patches, the resolution of the slip distribution in the upper part of the fault is poorer than that in the lower part, as seen in our checkerboard resolution tests. Moreover, there are trade-offs among the fault shape, smoothing factors and slip distributions. For these reasons, taking the straight fault as a first order simplification is reasonable.

Conclusions
Constrained by surface-based GPS observations of the M w 9.0 Sendai, Japan earthquake, we inverted for a coseismic slip model of the main earthquake and the fault slip in the 8 h after the main shock. The estimated geodetic moment magnitude is M w 8.98, close to the results inferred from seismic studies. With a maximum slip of 23.3 m, coseismic slips mainly occurred near the hypocenter and show a clear reverse behavior. However, slips with a strike-slip component were also found along the two sides of the slip concentration. The energy released by the fault slip in the 8 h after the main shock corresponds to an earthquake of M w 8.13. Having a slip maximum of ~1.5 m, the slips in this time period are mostly concentrated to the southwest of the main rupture, which shows clear correlation with the location and behavior of the M w 7.9 aftershock. In addition, with a magnitude of 0.2-0.4 m, these slips can also be found near the lower part of the main rupture, which may mostly be due to aseismic after-slip.