Source process with heterogeneous rupture velocity for the 2011 Tohoku-Oki earthquake based on 1-Hz GPS data

A rupture model with varying rupture front expansion velocity for the March 11, 2011, Tohoku-Oki earthquake was obtained by the joint inversion of high-rate Global Positioning System (GPS) data and ocean bottom GPS/acoustic (OB-GPS) data. The inverted rupture velocity with a complex distribution gradually increases near the hypocenter and shows rapid rupture expansion at the shallowest part of the fault. The entire rupture process, which lasted 160 s, can be divided into three energy release stages, based on the moment rate function. The preferred slip model, showing a compatible relationship with aftershocks, has a primary asperity concentrated from the hypocenter to the trench and a small asperity located on the southern fault. Source time functions for subfaults and temporal rupture images suggest that repeated slips occurred in the primary rupture, which is consistent with that from seismic waveforms. Our estimated maximum slip and total seismic moment are ~65 m and 4.2 × 1022 Nm (Mw 9.0), respectively. The large slip, stress drop, and rupture velocity are all concentrated at shallow depths, which indicates that the shallow part of the fault radiated high-frequency as well as low-frequency seismic waves.Graphical abstract Figure (a) shows the preferred cumulative slip distribution with a large slip concentrated in the shallow part of the fault near the trench, northeast of the hypocenter. The estimated maximum slip is ~65 m, and the seismic moment estimated with the variable rigidities based on the layered velocity structure used for Green’s function computation is ~4.2 × 1022 Nm (Mw 9.0). The region with slip of at least 5 m spans 150 km down-dip from the top boundary of the fault, corresponding to the depth from 10 km to 41 km, and extends along the strike for approximately 300 km. Figure (b) shows a heterogeneous rupture velocity distribution ranging from the minimum of 1.3 km/s to the maximum of 3.5 km/s on the fault plane. The rupture velocity within the vicinity of hypocenter shows a regular pattern, gradually increasing from 1.0 km/s to 2.0 km/s with increasing distance from the hypocenter. The rupture at the shallowest part of the fault significantly speeded up, with a velocity of larger than 3.0 km/s, consistent with that obtained by Lee et al. (2011). As for the southern part of the fault, 120 km away from the hypocenter, the rupture velocity is much smaller than in other regions, ranging from 1.0 km/s to 2.0 km/s. It indicates that the rupture front firstly arrived at the shallow part of the fault where large slip occurred, and then spread to the surrounding area. Figure (a) shows the preferred cumulative slip distribution with a large slip concentrated in the shallow part of the fault near the trench, northeast of the hypocenter. The estimated maximum slip is ~65 m, and the seismic moment estimated with the variable rigidities based on the layered velocity structure used for Green’s function computation is ~4.2 × 1022 Nm (Mw 9.0). The region with slip of at least 5 m spans 150 km down-dip from the top boundary of the fault, corresponding to the depth from 10 km to 41 km, and extends along the strike for approximately 300 km. Figure (b) shows a heterogeneous rupture velocity distribution ranging from the minimum of 1.3 km/s to the maximum of 3.5 km/s on the fault plane. The rupture velocity within the vicinity of hypocenter shows a regular pattern, gradually increasing from 1.0 km/s to 2.0 km/s with increasing distance from the hypocenter. The rupture at the shallowest part of the fault significantly speeded up, with a velocity of larger than 3.0 km/s, consistent with that obtained by Lee et al. (2011). As for the southern part of the fault, 120 km away from the hypocenter, the rupture velocity is much smaller than in other regions, ranging from 1.0 km/s to 2.0 km/s. It indicates that the rupture front firstly arrived at the shallow part of the fault where large slip occurred, and then spread to the surrounding area.


Introduction
Seismological and geodetic observations have been widely used to constrain seismic source parameters in modern seismology. Seismic waveforms including teleseismic waves and strong-motion records at different frequencies have been commonly used to estimate the finite fault rupture process model (Kikuchi and Kanamori 1982;Yoshida et al. 1996;Chi et al. 2001), whereas GPS measurements are used mostly to obtain the coseismic static slip distribution (Yabuki and Matsu'ura 1992). Joint inversion of seismic and geodetic data can integrate the kinematic features of seismic data and the static features of geodetic data (Delouis et al. 2002); however, the weighting between different datasets is always an issue.
Both seismic and geodetic tools for measuring ground motion during earthquakes have inherent advantages and limitations. Far-field seismic waveforms recorded by seismometers provide a smoothed estimate of the moment release rate and constrain the rupture direction but have less information on the near-field terms that reflect the detailed rupture history and static slip on the fault. Nearfield strong-motion acceleration waveforms are sensitive to the rupture history of earthquakes and the ground motion during earthquakes; however, they are prone to be affected by saturation for large earthquakes (Bilich et al. 2008), coseismic ground tilts, and errors owing to analog-to-digital conversion (Boore 2003). As a result, displacement waveforms obtained by double integration of strong-motion data often exhibit large drifts, making the accurate recovery of static displacements difficult (Boore et al. 2002;Boore and Bommer 2005). In contrast, high-rate GPS (hr-GPS) waveforms can directly capture both instantaneous and static displacements without saturation and large drift Bock et al. 2004).
The rapid development of a high-frequency sampling technique for GPS enables the rupture process inversion Open Access *Correspondence: wangzhen@eri.u-tokyo.ac.jp 1 Earthquake Research Institute, University of Tokyo, 1-1-1 Yayoi, Bunkyo-ku, Tokyo 113-0032, Japan Full list of author information is available at the end of the article for a finite fault by using hr-GPS waveforms. One primary advantage of rupture inversion using hr-GPS data is that it essentially fits both the dynamic and static features to avoid inadequate weighting in different data types as long as the Green's functions, including both dynamic and static responses, are computed. Thus far, relatively few works exist on rupture process models solely inverted from hr-GPS (Miyazaki et al. 2004;Yokota et al. 2009;, or on joint inversion combining hr-GPS with other datasets such as seismic waveforms  or static GPS displacement (Vigny et al. 2011). The available studies have shown the ability of hr-GPS data to infer the rupture process of mediumand large-sized earthquakes. However, most of the previous works that model the rupture process using hr-GPS data assumed a priori rupture velocity and solved a linear system for slip history during the coseismic period. If an appropriate rupture velocity and a subfault source time function are selected that are long enough to span the actual rupture time, the rupture velocity will have negligible influence on the final cumulative slip . However, Lee et al. (2011) showed that heterogeneous rupture velocity could influence the rupture process by the inversion of seismic waveforms. Konca et al. (2007) showed that non-uniform rupture velocities from 1.5 to 2.5 km/s are required to fit the seismic data for the 2005 Mw 8.6 Nias-Simeulue earthquake. In the present study, to better understand the rupture propagation process, we conduct rupture process inversions that allow for heterogeneous rupture velocity.
During the 2011 Mw 9.0 Tohoku-Oki earthquake, abundant continuous GPS stations in the Japan Islands enabled the inversion for slip history and rupture velocity using hr-GPS waveforms. Using hr-GPS data in single inversion  or joint inversion with other available data Bletery et al. 2014) for the rupture process has been performed with a constant rupture velocity to obtain a consistent main slip asperity with a simple rupture feature up-dip of the hypocenter to the trench. A rupture model with freely varying rupture velocity for each subfault has been presented by Minson et al. (2014) from joint inversion of geodetic data and tsunami data using a single-time-window method. In the present study, we use hr-GPS data in northeastern Japan and ocean bottom GPS/acoustic (OB-GPS) data to invert for the rupture process of the 2011 Mw 9.0 Tohoku-Oki earthquake. The analysis uses the inherently dynamic and static features of hr-GPS data and considers the inhomogeneous rupture velocity.
In the discussion, the preferred slip model obtained in present study is compared with the corresponding results of previous studies on the 2011 Tohoku earthquake. Our slip model shows that large slip is concentrated in the shallow part of the fault extending from the hypocenter to the trench, which is consistent with several previous studies. Furthermore, the distributions of slip, rupture velocity, and stress drop are compared to discuss the point from which the high-frequency or low-frequency seismic waves are possibly radiated and the corresponding radiation efficiency.

Observations
The Geospatial Information Authority of Japan (GSI) via Nippon GPS Data Service Company (NGDS) provided 2-h 1-Hz GPS data that span the 2011 Mw 9.0 Tohoku-Oki earthquake from 414 hr-GPS stations of the Global Navigation Satellite System (GNSS) Earth Observation Network System (GEONET) over northeastern Japan. GPS Solutions processed the raw GPS data and verified that the dataset was postprocessed in "real-time mode" with real-time orbits and clock corrections from the Veripos/Apex global service using RTNet software. We used the estimated GPS ground motion time series provided by GPS Solutions. The estimated maximum horizontal displacement exceeded 5 m, and the maximum vertical displacement was larger than 1 m. The average standard deviations of hr-GPS recordings during the 300 s prior to the mainshock were 0.017, 0.019, and 0.038 m for the east, north, and vertical components, respectively. The estimated data error of vertical displacement was larger than twice that for horizontal displacement.
From the 414 stations, we selected 53 stations distributed uniformly in northeastern Japan to invert for the rupture process of the 2011 Tohoku-Oki earthquake ( Fig. 1). To include near-field GPS waveforms, we selected several coastal stations located in northeast Honshu that were not used by . We extracted 700-s-long time windows from the hr-GPS data, including 300-s-and 400-s-long time windows before and after the earthquake origin time of 05:46:18 UTC March 11, 2011, respectively, as reported by the Japan Meteorological Agency (JMA). To reduce the offset before the mainshock, we subtracted the average value of the 300-s-long data before the origin time from the 700-s-long data. Delouis et al. (2010) reported that hr-GPS data low-pass-filtered at 0.03 Hz are appropriate for recovering the main characteristic of the rupture process for a mega-earthquake (Mw 8.8), and  demonstrated that very similar inverted results are found by using low-pass corner frequencies of 0.025-0.05 Hz. We thus applied a 0.05-Hz fourth-order low-pass Butterworth filter to the hr-GPS data and extracted waveforms of 300-s-long time windows after the origin time for the inversion. In the present study, we assumed that each station is uncorrelated; thus, the weight matrix for the hr-GPS data reduced to a diagonal matrix. The diagonal elements were set to the inverse of the standard deviation of the hr-GPS data, which means the weight of the horizontal component was about twice that of the vertical component. We used this weighting matrix for hr-GPS data in the inversion. Yagi and Fukahata (2011) proposed a new strategy to address the problem of uncertainty of Green's function in waveform inversion. Their method uses the corresponding data covariance with significant off-diagonal components, which is different from that assumed in the present study. Even though they reported that the modeling errors are of essential importance in waveform inversion analyses, such analysis is beyond the scope of the present study.
For the 2011 Tohoku-Oki earthquake, all available hr-GPS stations are located west of the rupture area and lack the constraints on the fault rupture near the Japan Trench. Seafloor geodetic measurements using the GPS/ acoustic technique above the rupture area reported by Sato et al. (2011) and Kido et al. (2011) with observation errors up to ~0.5 m can fill the gap. The peaks of horizontal and vertical displacements near the hypocenter were ~31 and ~3.9 m, respectively. Because the last preseismic and first postseismic GPS/acoustic measurement were conducted several months before the mainshock and about one month after the mainshock, respectively, the observed displacements include coseismic displacement of the mainshock and postseismic movement during the first month as well as coseismic displacement of foreshocks and aftershocks. However, the total displacement other than that of the coseismic signal by the mainshock has been reported as not larger than 1 m (Sato et al. 2011). Considering the signal-to-noise ratio, the OB-GPS data still provide valuable constraints on the static coseismic behavior near the Japan Trench. We thus used the coseismic static displacements at six stations shown in Fig. 1 with equal weighting in the joint inversion together with hr-GPS data. The relative weight between OB-GPS data and hr-GPS data is analyzed in "Results" section.

Fault parameters and method
We used a planar fault with a length of 600 km along strike and a width of 240 km along the dip with a uniform strike of 201° and a uniform dip angle of 11° to approximate the subduction plate interface following Zhou et al. (2014) (Fig. 1). Previous studies have demonstrated that the assumption of planar fault geometry does not cause significant modeling effects through inversions of static geodetic displacement (Zhou et al. 2014;Lee et al. 2011). We thus employed a single planar fault in this study. To obtain the seismic slip history, the fault plane was divided into 20 × 8 subfaults with dimensions of 30 km × 30 km. The hypocenter (142.8°E, 38.05°N, 24 km depth) determined by JMA is located at the center of the 11th subfault along the strike and at the center of the 6th subfault along the dip referring to the northwesternmost subfault ( Fig. 1).
We used a multi-time-window inversion method (Olson and Apsel 1982;Lee et al. 2006) to estimate the rupture process of the 2011 Tohoku-Oki earthquake. The moment rate function for each subfault was expressed by a linear combination of 20 symmetric triangular basis functions with expansion coefficients. The duration time of each triangle was 10 s, and the shift time of adjacent triangles was 5 s. Suzuki et al. (2011) indicated that the 78-s slip duration for each subfault is sufficient for stable inversion analysis of the 2011 Tohoku-Oki earthquake, whereas the 105-s slip duration in the present study is long enough to span the entire rupture process. The expansion coefficients for the triangular basis functions for all subfaults are unknown parameters in the inversion. To allow for heterogeneous rupture velocity, we assumed the hypocenter to be the initial rupture location. Moreover, we employed a sufficiently large a priori rupture velocity of 3.5 km/s that approaches the shear wave velocity following the "zero rupture delay time" approach of Lee et al. (2006). After the inversion with the a priori rupture velocity, we defined a rupture delay time for each subfault as the time of cumulative slip exceeding a threshold following Lee et al. (2006), which was measured from the initial rupture at the hypocenter. We then calculated the rupture velocity for each subfault as the shortest onfault distance from the hypocenter to the subfault divided by the rupture delay time. By definition, this estimated rupture velocity cannot exceed the a priori value.
To directly invert for the rupture model from the complete ground motion recordings, it is necessary to use Green's functions including both dynamic and static displacement responses. We used the frequency-wavenumber (FK) synthetic seismogram package of fk3.1 developed by Zhu to generate the Green's functions (Zhu and Rivera 2002). We used a layered velocity structure shown in Additional file 2: Table S1, which is based on the local velocity structure of Wu et al. (2008) and the Preliminary Reference Earth Model (PREM; Dziewonski and Anderson 1981) to calculate the Green's functions. In forward modeling, to compute the synthetic waveforms accurately, we divided each of the 30 km × 30 km subfaults into 25 finer subfaults with dimensions of 6 km × 6 km, and we obtained the Green's function as the mean of displacements owing to unit slip on the 25 finer subfaults. A 0.05-Hz fourth-order low-pass Butterworth filter, which was applied to the 1-Hz hr-GPS data, was also employed to filter the Green's functions to ensure common spectral content.
We imposed a spatial smoothing constraint on the triangular basis functions that span the same time period for adjacent subfaults to penalize unrealistically rough slip distributions. Boundary conditions along the edge of the fault can affect the slip inversion, particularly for shallow subfaults near the trench, owing to limitation of the spatial resolution of on-land GPS; however, the choice of smoothing operator, whether gradient or Laplacian, does not severely affect the slip inversion (Zhou et al. 2014). We thus employed a discrete Laplacian regularization with a free slip boundary condition along all the boundaries to smooth the fault slip. A trade-off curve between data misfit and slip roughness was used to determine the optimal value for the relative weight placed on fitting the data versus agreement with the smoothing constraint. In addition, to avoid the unphysical slip direction, we employed a non-negative constraint on slip to limit the rake range between 45° and 135°. We minimized the following objective function by using a non-negative least squares inversion (Lawson and Hanson 1995) where G is the coefficient matrix representing the linear theoretical relationship between the vector of observation data, d, and the vector of slip on the fault surface, s ; Σ d is the covariance matrix of observation errors; and L is the smoothing matrix obtained from Laplacian operator, and is the smoothing parameter that determines the relative weight of the smoothing constraint against the data fit.

Cumulative slip distribution
As described in "Fault parameters and method" section, we employed the spatial smoothing constraint and the non-negative constraint as a priori information to penalize unrealistically rough slip distribution. The smoothing parameter directly influences the slip distribution and slip amplitude. We used a trade-off curve between the weighted residual norm, Σ −1 / 2 d (d − Gs) , and the roughness norm of the slip, Ls , to determine the optimal smoothing parameter. To avoid changing the scale of the axes to modify the trade-off curve, we used a normalized trade-off curve to select the optimal smoothing parameter. A list of smoothing parameters from 10 −6 to 10 6 was tested. When the smoothing parameter is 10 −6 (10 6 ), 2 is equal to 10 −12 (10 12 ). According to Eq. (1), the result depends almost only on the observed data or the smooth constraint, respectively. A broad range of smoothing parameters was tested to support enough information to plot the trade-off curve and clearly show the location of the bend in the trade-off curve. We subtracted the minimum values from the weighted residual norm and the roughness norm, and we normalized them by their maximum values, respectively. In Fig. 2, the optimal value of 15.0 marked by point B was obtained by measuring the minimum distance from the normalized trade-off curve to the origin of the figure. The smoothing parameter that minimizes the Akaike's Bayesian Information Criterion (ABIC; Akaike 1980) was also plotted with point A. Point C was selected to show the comparison of points A and B.
(1) Figure 3 shows a comparison of three slip distributions inverted from three smoothing parameters corresponding to points A, B, and C of Fig. 2. The ABIC-based slip distribution (point A) was concentrated in a small region with slip sharply increasing from 20 m to a maximum of 80 m. The slip model corresponding to the smoothing parameter of point C extended in the entire fault plane with the slip smoothly increasing from almost 0 m to a maximum of 24 m. Although the ABIC-based method has been widely used in the source process inversion (Ide et al. 1996), it is actually incorrect for the case of employing non-negative constraints on slip, which is proved by the normalized trade-off curve as shown above. Fukuda and Johnson (2008) theoretically demonstrated this incorrectness and suggested that the ABIC-based method in general may significantly underor over-smooth the solution under the non-negative constraint.
We then examined the resolution of these datasets by using checkerboard tests. Two types of input model are shown in Fig. 4a, b, covering up to 3 × 3 nodes and 5 × 3 nodes. We assigned 0 and 10 m for slips at white and black grids, respectively. To create a set of synthetic data, we added observation errors generated from Gaussian distributions, with zero mean and standard deviation of 4 cm for horizontal hr-GPS, 2 cm for vertical hr-GPS, and 0.1 m for OB-GPS. The synthetic data were then inverted by using the same inversion strategy as that performed for the actual data. Figure 4c-f shows the slip distributions from the inversion of hr-GPS synthetic datasets generated by using the prior checkerboard slip model, with and without the aid of OB-GPS synthetic data. As shown in Fig. 4c, d, the inversion of hr-GPS data can almost perfectly retrieve the assumed slip model with good resolution in the deep part of the fault. The joint inversion of hr-GPS and OB-GPS data, shown in Fig. 4e, f, indicates that the OB-GPS data slightly improved the resolution of the fault plane in the shallowest central part. According to these figures, the obtained slip distribution using the actual data can recover the slip distribution fairly well with some cautions of 20-30% error at the shallower part of the faults. Figure 5a shows the preferred cumulative slip distribution with a large slip concentrated in the shallow part of the fault near the trench northeast of the hypocenter. The estimated maximum slip is ~65 m, and the seismic moment estimated with the variable rigidities based on the layered velocity structure used for Green's function computation is ~4.2 × 10 22 Nm (Mw 9.0). The region with slip of at least 5 m spans 150 km down-dip of the top boundary of the fault, corresponding to a depth of 10-41 km, and extends along the strike for approximately 300 km. The slip model shows that the deepest part of the fault did not slip significantly during the mainshock. Our preferred slip pattern is in good agreement with several previous results, such as the joint inversion of teleseismic body waves and static geodetic displacement (Kubo and Kakehi 2013), the inversion of tsunami waveform data (Satake et al. 2013), the inversion of strong-motion data (Suzuki et al. 2011), and the inversion of high-rate GPS data .
In Fig. 6, the preferred coseismic slip model is compared with the spatial distribution of the aftershocks. Most of the aftershocks within 24 h following the mainshock were located around the rupture area with slip larger than 30 m. The large coseismic rupture zone shows an excellent complementary relationship with the aftershock distribution. This feature is also consistent with the previous studies of other interplate earthquakes such as the 2012 Nicoya, Costa Rica, earthquake (Mw ~ 7.6;) and the 2005 Nias-Simeulue earthquake (Mw ~ 8.7;Hsu et al. 2006).

Rupture velocity
As described in "Fault parameters and method" section, the rupture delay time for each subfault can be determined from the time at which the accumulated slip exceeds a threshold, which is regarded as the time required to pass the rupture front. Lee et al. (2006) used a threshold of 0.1 m for the 1999 Chichi, Taiwan, earthquake (Mw 7.6). A higher accumulated slip of 0.5 m was employed as the threshold in this study by considering the poorer temporal resolution of the triangular basis functions with duration of 10 s and relatively large inversion uncertainty for this mega-thrust event. The rupture velocity was estimated as the ratio of the shortest distance between the hypocenter and each subfault to the rupture delay time at the corresponding subfault. Figure 5b shows a heterogeneous rupture velocity distribution ranging than 3.0 km/s, which is consistent with that obtained by Lee et al. (2011). The rupture velocity of the southern part of the fault 120 km from the hypocenter was significantly lower than that in other regions, ranging from 1.0 to 2.0 km/s. This indicates that the rupture front first arrived at the shallow part of the fault where large slip occurred and then spread to the surrounding area.

Temporal rupture process
The rupture process is shown in Fig. 7 by slip images obtained in 5-s increments. Similar to Lee et al. (2011), during the first 25 s of rupture, slip of about 5 m was limited in the region of the hypocenter. From 25 to 40 s, the rupture continuously propagated northeastward to the shallow part of the fault from the hypocenter, whereas the region around the hypocenter gradually became quiet. After a short break of about 40 s, the slip reoccurred at the hypocenter and propagated to the surrounding region and to the shallow part of the fault near the trench. Between 40 and 50 s, when the second rupture at the hypocenter began, the slip in the northeastern shallow part of the fault caused by the first rupture did not cease but continued and merged with the second rupture.
The significantly large slip occurred at the shallow part of the fault near the trench during 50-90 s. In addition to the shallow slip, the smaller slip located down-dip of the hypocenter occurred between 70 and 85 s. Figure 7 clearly shows that the second rupture occurred between 45 and 100 s and that a large seismic moment was released during the same period. After 100 s, slip only in the southern part of the fault continued until the end of the entire rupture, whereas the region around the hypocenter and the northeastern part of the fault was calm. Figure 8 shows the moment rate function for the entire fault and that for the individual subfaults. According to Figs. 7 and 8, we divided the entire rupture process into three time periods as defined in Fig. 8a. The percentages of seismic moment release during Periods I, II, and III were 9, 81, and 10%, respectively (Fig. 8a). This indicates that a major part of the total seismic moment was released during Period II, as shown in Fig. 8.
According to the rupture time periods on each subfault, we marked the rupture regions corresponding to Periods I, II, and III by a red rectangle (Region I), green rectangle (Region II), and blue rectangle (Region III), respectively, in  (Regions II and  III). At 40 s from the onset of rupture, the small moment rate in Region I indicates a transition from the first event (Period I) to the second event (Period II). In Period I, ~9% of the entire seismic moment had been released. During Period II, the Region I rupture repeated, and the moment rate rapidly accelerated in Regions I and II, releasing a large amount of moment. During Period III, the slip was concentrated in the southern part of the fault (Region III) and continued until the end of the mainshock with lower rupture velocities of 1.5-2.0 km/s.

Data fit
Eight hr-GPS stations, including four near-field and four far-field stations, were selected to show the comparison of the observed and synthetic recordings (Fig. 9). The comparisons of the observed and synthetic waveforms for all 53 hr-GPS stations are shown in Additional file 1: Figure S1. The synthetic waveforms reproduced the observation recordings fairly well. We also defined the misfit as �d obs −d syn� 2 �d obs � 2 , where d obs is the observed data vector and d syn is the synthetic waveform vector, to quantitatively evaluate the quality of the solution. Except for station 2001, which had the largest misfits for all the three components, the observed horizontal displacements were fit by the synthetic waveforms significantly better than the observed vertical displacements with misfit values of 0.0016 and 0.1403, respectively. The better fitting for horizontal displacements was somewhat expected because the weighting for horizontal displacements is twice that for vertical displacements. Obvious shifts between the observed and synthetic waveforms in the latter portions of many vertical waveforms, which are not found in horizontal waveforms (Yokota et al. 2009), resulted mainly in worse fitting for vertical displacement than that for horizontal displacement. Additional file 1: Figure S1 shows that the synthetic recordings could not trace the complex waveforms in the middle portions of several vertical waveforms, which can be attributed to the poor temporal resolution by modeling the moment rate functions with the triangular basis functions with a duration time of 10 s. In this study, static displacements of six OB-GPS stations were used in the joint inversion together with hr-GPS data. Issues of how to weight different types of datasets always occur in joint inversions. Our preferred weighting was obtained by performing a few tests to evaluate the displacement fits for OB-GPS data. Four sets of relative weights of 1:20, 1:50, 1:100, and 1:120 between hr-GPS data and OB-GPS data were tested. Additional file 1: Figure S2 shows a comparison of slip distribution inverted from these four different relative weights between hr-GPS data and OB-GPS data. Although an increase in relative weight from 1:20 to 1:50 increased the maximum slip and made the slip distribution more Aftershocks (Mw > 6) are shown by significantly larger circles compact, increasing the relative weight from 1:50 to 1:100 did not significantly change the slip distribution. Slip distribution inverted from the relative weights of 1:100 and 1:120 was almost identical. It is worth noting that although it is obvious that the fit to the OB-GPS data improved by increasing the weight of OB-GPS data, the final slip distribution was not be strongly affected by the relative weight in the range of 1:50-1:120. In addition, considering the possible contamination by postseismic deformation or aftershocks in the data of OB-GPS, we selected the relative weight of 1:100 between hr-GPS data and OB-GPS data in this study. Figure 10 shows that the predicted OB-GPS displacements effectively reproduced the observations for the relative weight of 1:100. Furthermore, we independently inverted the hr-GPS recordings without OB-GPS data to clarify the influence of OB-GPS data to the maximum slip. The inversion result showed a maximum slip of ~50 m, which is significantly less than that of our preferred rupture model. The OB-GPS data significantly increased the maximum slip. In addition, we compared the maximum slip inverted from different relative weights between hr-GPS data and OB-GPS data to determine the degree to which the OB-GPS data affects the maximum slip. The maximum slips were about 60, 64, 65, and 66 m for relative weights of 1:20, 1:50, 1:100, and 1:120 between hr-GPS data and OB-GPS data, respectively. Changing the weight of OB-GPS data from 20 to 120 slightly increased the maximum slip from 60 to 65 m.

Discussion
Our preferred rupture model had a seismic moment of 4.2 × 10 22 Nm, whereas other published models have those of 3.4 × 10 22 -5.8 × 10 22 Nm (Tajima et al. 2013). The maximum model slip was ~65 m, which is slightly larger than the 25-60 m in most previous results based on seismic waveform, geodetic, and tsunami data (Ozawa et al. 2011;Simons et al. 2011;Koketsu et al. 2011;Yokota et al. 2011;Saito et al. 2011;Suzuki et al. 2011;Lee et al. 2011;Zhou et al. 2014) and approaches that of 64 m obtained by Bletery et al. (2014). The main rupture area spans ~400 km along the strike and ~200 km along the dip. Under the free boundary condition, which allows the slip to reach the free surface, our preferred model showed that the greatest slip occurred up-dip of the hypocenter near the trench. This result is in agreement with that of other works (Lee et al. 2011;Suzuki et al. 2011;Bletery et al. 2014), although the amplitudes are inconsistent.
Our moment rate function characterized with three prominent slip events is consistent with that reported in previous research (Lee et al. 2011;Suzuki et al. 2011). The repeating slip behavior in Region I is consistent with the inversion result of solely strong-motion data by Suzuki et al. (2011) and the joint inversion of teleseismic, strongmotion, and geodetic data by Lee et al. (2011).  inverted for a rupture process assuming a constant rupture velocity from hr-GPS data. Although their result showed only one main slip event on the shallow part of the fault, the cumulative slip distribution was similar to that of our result. This indicates that if a priori homogeneous rupture velocity is selected and the a priori source time function for each subfault is sufficient for spanning the actual rupture time, the cumulative slip distribution is stable, but the rupture process is strongly affected by the rupture velocity. In addition, our results of the moment rate function and the slip images all showed that somewhat small slip at the deepest depths (sixth and seventh rows along dip in Fig. 8b) occurred only at the second rupture stage. These results are in agreement with Suzuki et al. (2011) but differ significantly from Ide et al. (2011). This discrepancy may be attributed to the inversion method. Suzuki et al. (2011) used a method of multiple time windows similar that used in the present study; however, Ide et al. (2011) implemented an empirical Green's function in the source process.
It is undisputed that a shallow region with large slip radiates very low-frequency seismic waves (Ide et al. 2011;Suzuki et al. 2011;Lee et al. 2011). However, the location in which high-frequency seismic waves are radiated is debatable. The vicinity of the edges of coseismic slip in which the rupture velocity changes abruptly emits highfrequency radiation (Madariaga 1977). Ide et al. (2011) considered that the shallow depths have a quiet rupture and that high-frequency seismic waves are radiated from deep depths. Simons et al. (2011) suggested that relatively high-frequency seismic waves are generated from the edge of coseismic slip at the deepest depths. The results of back-projection imaging show that the short-period energy is radiated down-dip from the hypocenter as well (Koper et al. 2011;Wang and Mori 2011). Dynamic source simulation suggests that the shallowest and deepest edges of coseismic slip all radiate high-frequency seismic waves (Duan 2012). By using kinematic inversion, Lee et al. (2011) demonstrated that deep depths as well as the shallow parts of a fault radiate high-frequency seismic waves.
A large decrease in stress was concentrated in the shallow region between the trench and the hypocenter obs. syn.

Fig. 10
Comparison between observed and synthetic static coseismic displacements, shown as red and black arrows, respectively, for hr-GPS stations, with average displacement of 270-300 s, and OB-GPS stations. The weight of OB-GPS data relative to hr-GPS data is 100. The left-hand figure shows a comparison of horizontal displacements with legends of 10 m for OB-GPS stations and 2.5 m for hr-GPS stations; that on the right shows a comparison of vertical displacements with legends of 2 m for OB-GPS stations and 0.5 m for hr-GPS stations (Fig. 11). The large slip exceeding 50 m, high rupture velocity, and slip repetition up-dip of the hypocenter inferred from our hr-GPS inversion indeed indicate that the shallow fault experienced a dramatic slip process and released a large amount of energy. The large stress drop (Fig. 11) and rupture velocity illustrate that the shallow part radiated high-frequency seismic waves. However, the longer duration and wider slip region relative to the deep depths demonstrate that the low-frequency radiations originated from shallow depths. According to the final rupture velocity distribution and temporal rupture process, part of the deep area had rupture velocities approaching 3.0 km/s, shorter durations, and abrupt velocity changes, which are characteristics of the highfrequency radiation in this area.
According to the crack model for a thrust earthquake, the relationship between radiation efficiency, η R , and rupture velocity, V R , can be considered as (Kanamori and Brodsky 2004) where c R and β are the Rayleigh wave and S wave velocities, respectively. We demonstrated that the shallow part of the fault with larger rupture velocity has higher radiation efficiency. Thus, the high-frequency seismic waves radiated from shallow (deep) depths up-dip (down-dip) of the hypocenter, and the shallow part had higher radiation efficiency.

Conclusions
Based on the assumption of multiple time windows with a large initial rupture velocity, we inverted hr-GPS data for a rupture model with varying rupture velocity for the 2011 Tohoku earthquake. By comparison of a slip model and rupture velocity distribution, the large slip areas were coupled with areas of large rupture expansion speed. Similar to the result from seismic waveforms, which are not limited by temporal resolution, the temporal rupture process and the total moment rate function showed three energy release stages during the event. The slip propagated from the hypocenter to the shallow part during the first energy release stage. In the second stage, the rupture in Region I that experienced a repeating slip released an extremely large amount of energy. The shallow part ruptured through the overall rupture duration time, forming a large rupture asperity with a maximum slip of about 65 m. The patterns of slip, stress drop, and rupture velocity demonstrate that high-frequency seismic waves radiated from both shallow and deep depths, whereas low-frequency waves radiated from the shallow fault. Therefore, varying the rupture velocity allowed us to obtain a more appropriate rupture model than that when using an assumed constant rupture velocity. Higher temporal resolution is preferred to determine the kinematic rupture process in detail by using hr-GPS data.

Additional files
Additional file 1: Figure S1. Comparison of observed (red line) and synthetic (black line) 300-s three-component hr-GPS ground motions observations for all 53 hr-GPS stations used in the inversion, where the EW, NS, and UD components are east, north, and vertical, respectively. Hr-GPS station numbers and peak absolute displacement amplitude in meters are shown on each waveform. Figure S2. (left) Slip models inverted from hr-GPS and OB-GPS data. Those inverted from relative weights of 1:20, 1:50, 1:100, and 1:120 between hr-GPS data and OB-GPS data are shown from top to bottom, respectively. (middle) and (right) The corresponding Ob-GPS data fits for horizontal and vertical displacements, respectively. The observations and synthetics are marked by red and black arrows, respectively Additional file 2: Table S1. Velocity structure used for Green's functions.

Fig. 11
Static stress decrease distribution for the 2011 Tohoku earthquake calculated by using a fast algorithm following Ripperger and Mai (2004)