Secular Crustal Deformation Characteristics Prior to 1 the 2011 Tohoku-Oki Earthquake of Mw 9 . 1 2 Detected from IGS Array , 2002-2011 3


 In order to reveal the secular deformation evolution and strain accumulation progress along subduction zone, we focus on high-precision Global Positioning System (GPS) data processing, GPS time series analysis and regional spatiotemporal filtering. The GPS position time series is modeled as a combination of the long-term plate movements, short-term noise and the frequency-dependent variations and tectonic deformation information. Common Mode Errors (CME) is removed and the deformation signals is then extracted from the residual time series by Principal Component Analysis (PCA). Combing the filtered displacement array and the dynamic regional strain array derived from GPS baseline, a spatiotemporal tectonic evolution map is exhibited obviously. We process IGS time series over the 2002-2011 period in Japan island and reveal the dynamic evolution of crustal deformation along subduction zone. The results show that the regional reference frame transformation by constraining one site location and one stable GPS baseline direction can further improve the Signal-to-Noise Ratio (SNR) of GPS observations. A integrated analysis combining the baseline length with azimuth change can show the crustal motion feature more comprehensively. The observed behaviors agree well with the simulation experiment results of the rock rupture in the laboratory. We divide the pre-seismic deformation into four stages, the stable linear strain accumulation stage, the formation stage of the local locked area, local decoupling stage, and strain release stage. We also give an appropriate dynamical model for the interplate coupling to explain the observed deformation characteristics. Research result provides some of the observational new evidence for inter-seismic deformation anomaly detection and the medium- to longer-term seismic hazard assessments.


Introduction 40
Japanese islands are located in the convergence of the Pacific plate, the north 41 American plate, the Philippine plate and the Eurasian plate, where plate tectonic 42 activities is strong, volcanoes and earthquakes is frequent. Different from 43 intra-continental plate earthquakes, the Pacific plate is under the Japan trench 44 (Fig.1). The Pacific plate is subducting westwards beneath the Eurasian plate at a 45 rate of about 80-90mm/yr and the Philippine Sea Plate is subducting beneath 46 47 convergence between the Pacific and Okhotsk(or north American) plates forms 48 an EW compressional stress field constraining the tectonic activity of the 49 region (Ohzono et al., 2013). Both shallow and deep source strong earthquakes in 50 this subduction zone are very active. The seismic slip rate in this area and along the 51 adjacent subduction zone to the south is about 1/4 of the plate convergence rate, 52 which has an important implication for the long-term seismic hazard in this 53 area (Kanamori et al., 2006). Contemporary Global Positioning System (GPS) 54 measurements have revealed a significant strain concentration zone in Tohoku 55 region in northeastern Japan (Miura et al., 2004). According to the historical 56 earthquake catalog from Global Centroid Moment Tensors(GCMT), a total of 6 57 earthquakes of over Mw 6.9 have frequently occurred in the region during

2.1.Modeling GPS Position Time Series 122
The original GPS observation array show significant fluctuation mainly due to 123 non-tectonic factors such as frame errors and un-modeled GPS processing errors, 124 seasonal effect, periodic change, or CME (Dong et al., 2006). In addition, the 125 coordinate time sequence is usually discontinuous or has large errors because of 126 bad observation environment or instrument obstacle. For these purposes, it has 127 essential importance to preprocess the data of GPS coordinate time sequence for 128 extracting tectonic deformation signals, such as eliminating the outfield value, 129 recovering the missing data in the sequence, and improving SNR by 130 spatiotemporal filtering. In this paper, the missing data is firstly interpolated by 131 cubic splines function, the time series with abnormally large residuals over 2 root 132 mean square (RMS) (>2.5, >2.5, and >6.5 mm in north, east, and up directions, 133 respectively) are regarded as outliers. The time series is modeled for removing the 134 trends and seasonal effects. CME is then removed and the deformation signal is 135 extracted from the residual time series by Principal Component Analysis (PCA). 136 The time series of daily GPS site coordinate are fitted as a combination of the 137 annual and semi-annual periodic terms, the initial position and a long-term trend 138 term and the step terms to obtain the residual.

2.2.Spatiotemporal Filtering and Signal Extraction 153
PCA approach is known as the empirical orthogonal function analysis, which 154 has been widely applied in geodesy, both to filter common mode noise in GPS 155 networks (Dong et al., 2006) and to detect weak regional tectonic deformation 156 signals by improve SNR in GPS coordinate time series (Ji and Herring, 2013;157 Kositsky and Avouac, 2010;Xu et al., 2015). For GPS residual time series with n 158 stations and m days, the matrix X(m×n) is constructed, and each column denotes 159 the residual value for N,E,U components for each station at all epochs, each rows 160 denotes N,E,U components of all stations at a given epoch. The elements at i-th 161 row, j-th column of the covariance matrix B is defined as, 162 where  is diagonal matrix with n nonzero eigenvalues. V is orthonormal matrix 167 with n eigenvectors.
X can expressed as, where k a , k v is the temporal variations and spatial response of k-th principal 172 component of X, respectively. The first few PCs usually represent the biggest 173 contributors to the variance of the network residual time series, related to the 174 common source time function. Therefore, CME can be expressed as, 175 where p is the first p PCs. 177

2.3.Dynamic Strain Array Calculation 178
Strain rate tensor is independent of CME and the reference frame, and have the 179 unique advantage for the detection of the crustal deformation. We calculate the 180 strain time series based on the high-precision GPS baseline vectors (see Fig.6). 181 Suppose that the plane baseline vector before deformation is ( , )

3.1.GPS Data Processing Results 214
This earthquake is recorded by a available modern GPS observation network 215 (http://geodesy.unr.edu/gsrm.php), unfortunately, most sites whose available time less 216 than 3 years prior to the earthquake in source area. The time scale is extremely lack 217 relative to the earthquake preparation scale, thus not enough to investigate the 218 pre-seismic crustal motion and evolution. Therefore, GPS observation data used in 219 this study are mainly from IGS network with 7 continuous stations observed from orbits. Dual-frequency phase data are used to compute double-difference 224 ionosphere-free residuals, which are then inverted through a least squares 225 procedure for solving the unknown parameters including single-day positions, 226 satellite orbital elements, phase ambiguity integer numbers, tropospheric delay 227 residuals, etc. Uncertainty of 10 mm for the carrier beat phase data,30 s of sampling 228 intervals, relax mode for the satellite orbit, the IERS2003 model for earth gravity 229 field, solid tide and pole tide, the FES2004 mode for the global ocean tide, the GMF 230 for troposphere mapping functions, and the parameter of zenith delay is estimated 231 every 2h. 232 completely. AIRA, the furthest site from the epicenter with about 1300 km, the 236 coseismic displacement is close to zero, with less influence from the earthquake, 237 whereas the stations MIZU and KSMV, about 140km and 280km from the epicenter, 238 have more significant coseismic effect, with -1.25m, 2.35m, -0.11m and 0.06m,0.71m, 239 -0.26m in N,E,U direction, respectively. Comparing with MIZU, the vertical 240 displacement of KSMV is more significant than the horizontal displacement. 241

267
For the clarity of the internal deformation, we transform the station 268 coordinates from ITRF frame into local reference frame by similarity

300
The secular displacement rate is mainly caused by plate movement. In order to 301 illustrate the internal tectonic movement feature of this region, we calculate the 302 velocity vectors with respect to the stable Eurasian plate based on the Eurler 303 parameters (55.25°N,98.3°W,0.256°/Ma) from Xu and Wu (2014), as shown in Fig.8. 304 The Japanese Islands show a dominated westward motion in the region with a 305 average velocity of 20mm/yr, due to the subduction of Pacific plate and Philippine 306 Sea plate. Comparing with the horizontal motion, the vertical motion is not 307 significant, with about 1mm/yr. 308

3.3.Filtering Results of the Residual Time Series 312
The residual time series in Fig.7 still show a significant spatial correlation 313 errors among different stations, so the SNR of the GPS residual time series need to 314 be further improved by PCA spatiotemporal filtering. Fig.9 shows the first 5 315 eigenvalues normalized by the highest one in N, E, U components. The first mode 316 have significant eigenvalue and other modes are less than 50% of that of the mode 317 1 signal for the N, E, U components.  Fig.10, Fig.11 and 322 Fig.12, respectively. For all three components, their PC1 in Fig.10, Fig.11 and 323 Although there is no generally accepted consensus on calculation CME at a large 326 scale, for a small GPS network with a spatial scale of 5• × 5•, the mode is considered 327 as CME where most sites (more than 50%) stations have significant normalized 328 responses with larger than 25% and the corresponding eigenvalues exceed 1% of 329 the summation of all eigenvalues (Dong et al., 2006). According to this criterion, we 330 regard the PC1 is regarded as CME, and focus on other modes PCs for tectonic 331 deformation analysis. 332 20 of 30 The temporal variations and spatial response of other modes PC2,3 in the E, N 333 and U component are shown in Fig.10, Fig.11 and Fig.12.right. For all three 334 components, their temporal variations show the a long-period turning motion 335 process, from westward to eastward, from northward to southward, from upward 336 to downward, respectively. The spatial distribution have the opposite sign in the 337 east and west sides of network. The station KSMV has the most significant spatial 338 responses in three components, which are explained that the station is located on 339 the boundary of ocean-continent, where has the strong impact from subduction. of Japan island has strong impact from subduction zone. The opposite spatial 353 distribution in the east and west sides of Japan island probably be explained that 354 the strain accumulation will release and turn to move in a oppose direction when 355 the deformation reaches the limit. 356  The arrows and solid green curve are defined as in Fig.10.  finally. According to the theory of rock mechanics and deformation, the whole 388 observed behaviors are described as stable linear motion and instable nolinear 389 motion. In first stage, the deformation show a trend of linear increase, referred to 390 as a elastic deformation. The deformation will restore quickly, once the external 391 force is unloaded. In the second stage, the deformation starts to deviate the linear 392 movement trend, both of the elastic and inelastic deformation coexist. Contrary to 393 the first stage, the deformation will not restore completely if the external force is 394 unloaded. In the stage, the strain of different section on subduciton zone begin to 395 appear the features of the differentiation change. The displacement and the 396 dilation strain rate of part of areas decrease gradually, close to the locked status. 397

3.4.Regional Strain Time Series 357
After a short period of locking, the part locked areas will decouple and appear the 398 reverse motion trend. Meantime, the shear strain transfer to the remaining locked 399 region and become more and more concentrated, resulting in its acceleration 400 accumulation on the whole region, until the occurrence of the earthquake. Our 401 findings are consistent with the simulation experiment results of the rock rupture 402 in the laboratory (Ma, 2016). 403 According to the above spatio-temporal evolution process of crustal motion we 404 have observed and the previous research (Chen et al., 2013;Hoshiba, 2006), we give 405 an model for the interplate coupling to explain the present three-dimensional 406 deformation, as shown in Fig.14. The evolution process is divided into the 407 following four stages. 408

411
(1) The stable linear strain accumulation stage( Fig.14.a). The preliminary linear 412 movement can be explained as the relatively stable motion between two plates in 413 the subduction process of the ocean-continental plate at the beginning of the 414 earthquake preparation. Because the subducting ocean plate is getting stuck, the 415 overriding continental plate will be squeezed and dominated by compression 416 deformation, the leading edges is dragged down, while the hinterland areas in 417 near-field surface is bulged upward. As a result, the strain near the boundary of the 418 plate starts to accumulate, and the seismogenic region is formed initially. 419 (2)The formation stage of the local locked area (Fig.14.b). With the continuing 420 of the relative motion between plates in far-field, the frictional stress level gradually 421 increases due to the heterogeneity of the ocean-continental subduction belt contact 422 surface. The ability to resist the deformation gradually reduces, in this case, the 423 deformation deviates from the originally linear trend. The local seismogenesis 424 region presents a deformation characteristic of gradual attenuation from deep to 425 shallow. The westward and upward displacement of the region gradually slows 426 down and reaches the locked state, whereas the remaining regions still keep on 427 relative movement under the action of tectonic stress, so the stress begin to transfer 428 gradually to the locked area and become more and more concentrated. 429 (3) The local decoupling stage (Fig.14.c).The strain accumulation and the stress 430 level in the locked region continuingly increases and reaches the limit finally. The 431 seismological element is decoupled to become the slip area, the part strain energy 432 begins to release slowly, and the surface deformation then appears the reverse 433 motion. In this stage, there is usually pre-slip and increased small earthquake 434 activity. However, the shear strain energy still keep continuing accumulation 435 during locked status. The remaining major locked area get the further concentration 436 of strain and redistribution of stress, resulting in the accelerated accumulation of 437 the shear strain and increasingly high stress concentration degree. 438 (4) The great rupture and strain release stage (occurrence of great earthquake) 439 ( Fig.14.d). With more and more concentrated stress and strain accumulation in the 440 remaining locked area, the frictional resistance between plates is overcome 441 completely. As a result, the widely distributed locked zones become instable and 442 appear the instantaneous mutation and a large strain release, the elastic rebound 443 happen on the leading edge of the jammed overriding plate, where the obduction 444 and upheaval cause the tsunami phenomenon, while the bulge behind the leading 445 edge is collapsed. 446 In the different stages, the deformation characteristics and the time experienced 447 are different. After four stages of evolution, the whole seismic process is completed. 448

Conclusions 449
The SNR in coordinate time series for regional GPS networks can be further 450 improved by spatiotemporal filtering from the residual GPS observation array. 451 Integrating filtered displacement time series, regional strain, baseline change and 452 processing of similarity transformation, the significant pre-seismic deformation 453 evolution around a period is detected prior to 2011 Tohoku-Oki earthquake of Mw 454 9.1. The displacement and the dilation strain rate of part regions decrease 455 gradually, close to the locked status, then reverse motion, instead, the shear strain 456 is more and more concentrated, resulting in the acceleration accumulation. The 457 tectonic evolution for earthquake occurrence along the subduction zone can be 458

Authors' contributions 479
Keke Xu designed the study, performed the methodology research, and drafted the manuscript. Rong He 480 performed data processing and analysis. Kezhao Li performed the validation and assessment of GPS data.

Competing interests 490
The authors declare that they have no competing interests.