Research ArticleGEOPHYSICS

Along-trench variation in seafloor displacements after the 2011 Tohoku earthquake

See allHide authors and affiliations

Science Advances  19 Jul 2017:
Vol. 3, no. 7, e1700113
DOI: 10.1126/sciadv.1700113


The 2011 Tohoku-oki earthquake was the largest earthquake ever observed with seafloor geodetic techniques in and around its source region. Large crustal deformation associated with both the coseismic rupture and the rapid postseismic deformation has been reported. However, these observations are insufficient to describe the postseismic deformation processes occurring around the broad rupture area. We report the first results of seafloor Global Positioning System and acoustic ranging (GPS-A) observations based on repeated campaign surveys conducted over nearly 4 years using the extended GPS-A network deployed along the Japan Trench in September 2012. The observed postseismic displacement rates (DRs) show evident spatial variation along the trench: (i) distinct landward DRs in the large coseismic slip area [primary rupture area (PRA)], evidencing the predominance of viscoelastic relaxation; (ii) remarkable trenchward DRs in the south of the PRA, indicating rapid afterslip; and (iii) slight trenchward DRs in the north of the PRA. These features provide great insights into constructing a more complete model of viscoelastic relaxation, and they also indicate spatial variation of afterslip and fault locking along the plate interface with clear spatial resolution, providing invaluable information for the improvement of seismic hazard assessment.


Postseismic deformation following a great subduction earthquake is controlled by multiple physical processes, of which viscoelastic relaxation (VR), afterslip (AS), and fault locking (FL) are the most important (1). VR and AS are transient responses to the impulsive coseismic stress change, and they are governed by the rheological structure of the earth interior (2) and frictional properties on the plate interface (3), respectively. FL, also representing fault nature, causes strain accumulation due to the subduction of an incoming plate.

The 2011 Tohoku-oki earthquake (hereinafter, the Tohoku earthquake) shows large crustal deformation not only in the coseismic period (46) but also in the postseismic period (4, 79). The previous reports on the postseismic deformation after the Tohoku earthquake have shown complex patterns of the deformation field. Extensive trenchward movement was observed by onshore Global Positioning System (GPS) networks (4, 7), whereas the directions of motion varied considerably in the offshore area. Trenchward motion was also observed in the southern offshore area, but significant landward motion, that is, completely opposite to that observed onshore, was observed in the primary rupture area (PRA) (Fig. 1) by offshore GPS-A [combined geodetic technique of GPS and acoustic ranging (10)] measurements (8).

Fig. 1 Observed postseismic DRs after the Tohoku earthquake through the repeated GPS-A observations.

Black and green vectors show the DRs with 1σ error ellipses relative to the North American plate estimated by this study and by the Japan Coast Guard (8), respectively. The observation periods are from September 2012 to May 2016 for black vectors and from September 2012 to January 2014 for green vectors. Gray vectors in the onshore region show the DRs relative to the North American (NA) plate from September 2012 to November 2015 at the GNSS Earth Observation Network (GEONET) sites. Yellow and red contours represent 20- and 50-m contours of the coseismic slip distribution (40), indicating the PRA and very large slip area (VLSA), respectively.

The observed deformation pattern in the early postseismic period has been modeled in terms of VR and local AS (1115). It has been established that VR was an essential process in generating the landward movement in the PRA, whereas AS in the downdip extension of the PRA and below the southern offshore area was required to explain the trenchward movement, suggesting strong spatial variation in the postseismic deformation processes. However, the available offshore GPS-A observation sites are located mostly around the PRA, and coverage was insufficient to consider the entire scope of postseismic deformation across the broad rupture area of the Tohoku earthquake (13, 14).

To broaden our scope, we constructed a dense and wide GPS-A observation network along the Japan trench (fig. S1) in September 2012 (16); we conducted six campaigns at most sites until May 2016. On the basis of these new observations, we have already reported faster postseismic displacement rate (DR) at a site on the incoming Pacific plate than the steady plate convergence rate, and we have shown that the internal deformation due to VR accounts for the observed rapid DR (9). This study elucidates the postseismic DR field along the landward area of the trench based on the new GPS-A data, with the objective of analyzing it in terms of the elementary processes controlling postseismic deformation (that is, VR, AS, and FL).


The obtained postseismic DR field shows remarkable spatial variation (Fig. 1). Because the trench-normal component is dominant in the observed DRs, we show the profile of the trench-normal component of the DRs along the trench in Fig. 2. The DRs are consistently directed landward among the stations above the PRA, whereas prominent trenchward movement is found to the south of the PRA. To the north of the PRA, trenchward DRs are also observed, but their magnitudes are smaller than those in the southern area. We compared the observed DRs with those expected based on the existing VR model (VR-DRs) (11) in Fig. 2 to examine how VR contributes to the observed postseismic deformation. VR-DRs can approximately explain the observed landward DRs, although some discrepancies are evident. Outside the PRA, where the contribution of VR to deformation is small, the residuals between the observed DRs and VR-DRs are expected to reflect the contributions of FL and/or AS. In the following discussion, we divide the observation region into five segments (Fig. 2) on the basis of the magnitudes and signs of the residuals to identify the most relevant factor.

Fig. 2 Comparison of the observed DRs and the estimates by the VR model and FL.

(A) Black vectors are the same as those in Fig. 1. Blue, red, and green vectors show VR-DRs, the modified VR-DRs, and VR-DRs with contributions of FL assuming full-locking conditions, respectively. (B) DRs of 14 sites [shown as gray squares in (A)] close to the trench projected on line AA′ [shown in (A)] to investigate the spatial characteristics near the trench. Horizontal axis indicates the trench-normal component (113°N) of the DRs. Vertical axis indicates the distance in kilometers along the line AA′ shown in (A). The observed DRs in the trench-normal component are represented as black squares with 1σ errors. Colored curves are smoothed cubic spline fitting of the modeled DRs: VR-DRs (blue), the modified VR-DRs (red dashed), and VR-DRs with contributions of FL (green dashed). Horizontal orange lines show northern and southern limits of the PRA.


The largest discrepancy between the observed DRs and VR-DRs can be found in the southernmost segment (Segment-5), where the observed trenchward DRs are obviously inconsistent with the motion expected from the VR model. The observed DR requires a deformation process generating trenchward movement. However, AS is only the process among the three possible fundamental postseismic processes causing significant trenchward motions near the trench; therefore, the observed DRs must be manifestations of the occurrence of extensive AS in the segment. Our results show clear spatial variation of the DRs within this segment; the DR is largest at G17, and it decreases toward the south or southwest, as evidenced at G19 and G20. The observed spatial variation is interpreted as the result of the heterogeneous distribution of AS, that is, large AS occurs at the shallow plate interface, but it is somewhat diminished in the deeper portion of the megathrust and in the southernmost part of the study area. Although previous studies on postseismic deformation (1215) and on small repeating earthquakes (17) have already highlighted the occurrence of AS in this segment (Fig. 3), its spatial distribution has not been well constrained. Our results are direct evidence of the occurrence of AS and of its spatial variation.

Fig. 3 Comparison of the observed postseismic DRs with postseismic slip estimated by previous studies.

Black and green vectors and yellow and red contours are the same as those in Fig. 1. Black rectangle shows the fault zone of the 1896 Meiji-Sanriku earthquake (Eq.) (18). Blue contours show the AS distribution (12). The cumulative postseismic slip for the 23 April 2011 to 10 December 2011 period derived by small repeating earthquakes (17) is shown as colored patches. Gray patches indicate no aseismic slip during the period. Black curve denotes the northeastern limit of the Philippine Sea plate (42).

In the northernmost segment (Segment-1), the observed DRs show small trenchward movement of approximately <5 cm/year, comparable with the VR contribution, suggesting that the contributions of AS and FL are insignificant. However, an earlier study on small repeating earthquakes (17) has suggested the occurrence of AS on the western side of the segment, whereas the segment includes the rupture area of the 1896 Meiji-Sanriku earthquake [moment magnitude (Mw), 8.0] (18), where FL on the shallow plate boundary could account for the occurrence of the earthquake (Fig. 3). Therefore, it is important to examine whether AS and/or FL contributes to the observed DRs. Here, we estimated the DRs by FL with various distributions (see Materials and Methods and fig. S3) and found that the observed small DRs do not support strong FL on the shallow fault. Similarly, forward slip, causing DRs with opposite directions to the DRs by FL, would not explain the observations when significant forward slip is assumed on the shallow interface. The DRs reported here are insensitive to slip on the deep fault, and they do not contradict the existence of AS represented by the increased activities of small repeating earthquakes. Therefore, it is probable that DRs with different components, that is, trenchward due to AS and landward due to FL, largely cancel each other, resulting in insignificantly small DRs. However, this conjecture requires broader coverage of data on DRs to elucidate their variations in the trench-normal direction.

The landward DRs observed in Segment-2, Segment-3, and Segment-4 in the PRA and its surroundings are generally larger (approximately >10 cm/year) than the VR-DRs. This prompted us to examine the reasons for the additional landward deformation required to fill the gap between the observed and modeled DRs. One possible reason is the difference in viscosities of the model and the actual medium, and another is the contribution of FL. Here, we show two end-membered models: the VR model assuming lower viscosities (modified VR model; see Materials and Methods) that maximally explains our observations and a model that includes deformation that is caused by VR and complete locking uniformly on the plate boundary fault (VR model with full locking; see Materials and Methods). The modified VR model could reduce the discrepancies in Segment-2 and Segment-4 (Fig. 2), as could the VR model with full locking. Comparison of our observations and these models suggests that adjustment of the rheological model for VR estimation and/or introduction of FL should be considered in the postseismic deformation modeling. However, these two new models increase the residuals in Segment-3, which we will discuss later.

Although lower viscosity or the presence of FL appears to account for the large landward DRs in the PRA, the landward components of the DRs observed in the northern half of Segment-2 remain systematically large, and it is suspected that the viscosity of the assumed model is also too large. However, a model with even lower viscosity would be inconsistent with the observed DRs without significant trench-parallel components because the VR-DRs in the area have significant trench-parallel components (Fig. 2). Therefore, we have to consider different factors to enhance the landward components of DRs in the area. Here, we suggest that the PRA would be expanded more to the north than assumed in the previous VR modeling. Several slip inversions using tsunami data (19, 20) have yielded coseismic slip models with larger spatial extents than those based on geodetic modeling. As shown in fig. S2, the landward components of the VR-DRs are increased by including additional coseismic slip patches to the north of the assumed PRA (that is, around Segment-2), whereas the trench-parallel components are not. However, the idea suggested here must be verified quantitatively using both coseismic and postseismic deformation data with a realistic rheology structure because the tsunami source in the north might not be relevant to the faulting (21) and because the spatial extent of the coseismic rupture remains under dispute.

The relatively smaller DRs observed in Segment-3 suggest various heterogeneities, which have not been considered thus far. The heterogeneity in the rheological structure would be responsible for the unevenness of the observed VR, although short-wavelength variations of viscosity in the mantle are unlikely. It would be more reasonable to interpret the variation of landward DRs as the result of the spatial variation of FL strength. Here, we make a calculation using a simple model of variable FL to fit for the observations, in which the plate interface in the PRA is completely locked except for the VLSA (Fig. 1) where zero coupling is assumed; this practically assumes temporal delay of coupling recovery in VLSA. The calculated DRs roughly reproduce the observed DRs (fig. S4). Although the present model is overly simple, the spatial variation of the observed DRs suggests that introducing heterogeneous FL in the PRA would help explain the observations, which imply complexities of fault relocking after the mainshock.

As explained, the results of the GPS-A observations presented in this study have great potential for constraining the stick/slip condition along the plate boundary fault, as well as the rheological behavior. Therefore, it is necessary to construct a comprehensive postseismic deformation model that is able to explain the observed features to further assess the future seismic hazards in the subduction zone. Furthermore, the temporal variation of the deformation field must be important for further quantitative investigation of the various fundamental postseismic processes, and we will continue to conduct the GPS-A surveys in the area.


GPS-A seafloor geodetic observations and analysis

GPS-A observation is a unique method for measuring seafloor horizontal displacement in global geodetic coordinates—that was developed conceptually by the Scripps Institution of Oceanography (10, 22). The GPS-A observation system is composed of kinematic GPS positioning on a sea surface platform and acoustic ranging between the sea surface platform and a seafloor acoustic unit. Here, research vessels with GPS antennas and a hull- or side-mounted acoustic transducer were used as sea surface platforms. The seafloor acoustic unit at each GPS-A site consists of three to six seafloor transponders forming a square or triangular array. The detailed measurement settings are the same as those in our previous study (9).

Absolute positioning of the GPS antenna on the vessel was estimated at 0.5-s intervals using the long-baseline kinematic GPS analysis by IT (Interferometric Translocation) software (23). The reference onshore GPS station (AOBL) for the baseline analysis was installed at Aobayama campus of Tohoku University, which is located in the city of Sendai in Japan (fig. S1). The daily coordinates of the reference site in the ITRF2008 reference frame (24) were estimated by precise point positioning analysis (25) using the GPS data processing package GIPSY-OASIS version 6.1.2 (26). We determined the absolute position of the vessel’s transducer using predetermined antenna-transducer geometry and the vessel’s attitude as measured by a GPS gyro composed of three GPS antennas and/or a dynamic motion sensor.

The acoustic ranging system measures round trip travel times between the vessel’s transducer and the seafloor transponders. We used an acoustic signal with a duration of 25 ms using a 10-kHz carrier wave, whose phase was modulated by the seventh-order M-sequence. We detected precise travel times by correlating the transmitted and received signals at a sampling rate of 100 kHz.

To estimate precise seafloor positions, we conducted two types of survey at each site: One was a moving survey to determine the initial positions of the individual transponders for configuring the array geometry, and the other was a point survey at the center of the array for array positioning (27, 28). In both surveys, a stratified sound-speed structure was assumed, and temporal variation of its vertically averaged quantity was estimated for each ping. Using point survey data, we estimated the horizontal array position for each ping relative to the predetermined position, based on the assumption that the array geometry was fixed during the observation period. Note that the array positioning using data of point surveys specializes in determining horizontal movement but not the vertical one because of its intrinsic trade-off with sound speed, whereas some of the previous GPS-A studies successfully resolved vertical movement using substantial data of moving surveys (8, 29).

It is known that the estimated array position in a campaign would fluctuate over a 0.5- to 3.0-hour time scale because of the oscillation of the lateral sound-speed structure caused by internal gravity waves (9, 22, 27). To reduce the effect of this fluctuation, we performed a point survey for about 6 to 12 hours at each site, which amounted to hundreds or thousands of pings. Then, we took their average as the final position and their 1σ SD as the error ellipse for the campaign. The estimated array positions and the total numbers of pings in each campaign are listed in table S1.

Aftershock correction

Because the time series of the observed array positions could include coseismic steps caused by large aftershocks of the Tohoku earthquake, we removed their contribution to extract pure contributions due to AS, VR, and FL. Evaluating the effects of the aftershock at each site, we found only two successive aftershocks (Mw 7.1 and Mw 7.2) that occurred on 7 December 2012 (fig. S1), which caused non-negligible displacement at several sites. We calculated the horizontal coseismic steps of these aftershocks in a uniform elastic half-space (30) for point sources of Global Centroid Moment Tensor (GCMT) solutions. The contributions of the two sources were merged and corrected for when the total steps exceeded 1 cm. The calculated coseismic displacements for GPS-A sites are listed in table S2. The maximum coseismic step due to the aftershocks was found to be 8 cm in the west-northwest direction at G12. The time series before and after the correction are shown in fig. S5.

We also investigated the responses of VR due to the large aftershocks on 7 December 2012 using VISCO1D software (31). We constructed a viscoelastic earth model composed of three layers: an elastic layer with a thickness of 50 km and two Maxwell viscoelastic medium representing the asthenosphere and upper mantle, with thicknesses of 170 and 430 km, respectively. The elasticity parameters for the medium were set on the basis of the Preliminary reference Earth model (32). Viscosity in the upper mantle layer was defined as 1.0 × 1020 Pa∙s based on past VR studies in the Tohoku-oki region (33), and the viscosity in the asthenosphere mantle layer was determined to be 1.5 × 1019 Pa∙s by an optimization mentioned later. We used the fault parameters of the aftershocks estimated from the GCMT solutions using a scaling law (34) for the VR calculation. The total displacement caused by the VR of the aftershocks from the time of occurrence of the aftershocks to May 2016 for each site was about 1 cm at most, and therefore, it could be ignored in this study.

DR calculation

Because the onshore and offshore geodetic observations (7, 8) reported that the rapid decay of the postseismic deformation was limited to approximately within 1 year after the mainshock, we assumed constant postseismic DRs when considering the accuracy of our GPS-A observations. The postseismic DRs in the northern and eastern components at each site were calculated from the time series of array positions corrected for the aftershocks. We performed a weighted M-estimation method (that is, a robust linear regression technique) using Tukey’s biweight function when the time series had five or more observation campaigns. Because such a robust linear regression method is unsuitable for handling small numbers of data, we performed a simple weighted linear regression for the cases with fewer than five campaigns. This robust fitting method has been used in various GPS-A studies by other research groups (35, 36). Because the obtained postseismic DRs are in ITRF2008, we converted them relative to the North American plate on the basis of the MORVEL model (37) by subtracting the North American plate motion at each site (derived by the no-net-rotation MORVEL56 model) (38) from the estimated DRs. The postseismic DRs relative to the North American plate are listed in table S3.

We also calculated the DRs from the time series data in a former study provided by the Japan Coast Guard (8) using an M-estimation method (Figs. 1 and 3). Because the observation period of that study was from April 2011 to January 2014, we used only those data from after September 2012 to fit our data period. Furthermore, the VR-DRs (Fig. 2) were calculated from the time series data predicted by the VR model (11) during our observation period for each site (blue lines in fig. S5) using simple temporal averaging. The DRs of onshore GEONET sites (Fig. 1) were calculated from the F3 coordinate solutions for the same observation period provided by the Geospatial Information Authority of Japan (39) using a classical linear fitting method. These rates were also transformed into the North American plate reference as outlined above.

Modification of the VR model

The rheological structure used in the VR model (11) is optimized to explain the onshore and offshore geodetic observations. However, it was constructed considering only a few offshore GPS-A sites located little far from the trench. Therefore, the viscosities in the rheological structure might not be constrained well to explain our observation sites located near the trench because the spatial change in the rheological structure normal to the trench may not be small. Then, we tried to adjust the viscosities of the VR model only from our offshore observation results and to test whether our results could be explained only by the contributions of VR, or not. Because the VR-DRs in a structure with lower (higher) viscosity generally show higher (lower) DRs, we modified the VR-DRs by simple amplification. To minimize the misfits between the observed and modified VR-DRs, we estimated a multiplication factor of 1.32, indicating slightly lower viscosities. This operation practically lowered viscosities both in the mantle wedge and in the oceanic mantle from the original VR model. We used the operation for the sake of simplicity, although it is desirable to adjust only the oceanic mantle viscosities because the significant landward motion above PRA was strongly controlled by the oceanic mantle viscosities (13, 15). We show the modified VR-DRs in Fig. 2.

Investigation of the additional coseismic slip patches

To compensate the misfits between the observed DRs and the VR-DRs in Segment-2, we investigated the contributions of additional coseismic slip patches at the shallow portion of the plate interface. The coseismic slip modeling study (40), used in the VR model, indicated that the coseismic GPS and GPS-A observation data allow the additional coseismic slip of approximately 20 m in this segment. Thus, we laid out the additional coseismic slip patches with a slip of 20 m and a rake of 81° (fig. S2, A and C) and calculated their viscoelastic responses using VISCO1D software (31) and the Earth model mentioned above. Although it is desirable to calculate the viscoelastic responses based on the same rheological structure used in the VR model (11), we assumed the layered VR structure model for the simplicity, which can calculate the viscoelastic responses similar to the VR model by optimizing the viscosity of the asthenosphere mantle layer by following these steps: (i) calculate the viscoelastic responses to the coseismic model of the Tohoku earthquake (40) at the GPS-A sites varying the viscosity in an increment of 0.5 × 1019 Pa∙s and calculate the DRs similar with the calculation of the VR-DRs mentioned above and (ii) search the optimal viscosity by minimizing the VR-DRs and the DRs calculated by our layered VR structure models. As a result, the asthenosphere viscosity that can reproduce the viscoelastic responses of the Tohoku earthquake similar to the VR model is determined to be 1.5 × 1019 Pa∙s (fig. S2B). Then, we calculated the viscoelastic responses of the additional coseismic slip patches, and they caused the moderate but further landward motions (<1.5 cm/year) near the trench that close the misfits between the observation and the calculation shown in Fig. 2. The motions are a little small to completely compensate the misfits. However, because other coseismic slip modeling study (20) that used both the geodetic and the tsunami data allowed the coseismic slip to be up to 30 m in this segment, the coseismic slip amounts could be varied depending on models. If the additional coseismic slip patches with a slip of 30 m, the additional landward motions would be amplified to ~2 cm/year. Furthermore, as we pointed out in the Discussion, the viscosities used in the VR model might need to be revised, which can magnify the contributions of VR. Thus, while the amounts of deformation due to VR caused by introducing the additional patches are still controversial, we have shown that the additional coseismic slip certainly improves the fit of the model to the observations.

Investigation of the FL contributions

We estimated the DRs caused by FL on the entire or local areas of the plate interface assuming full-locking conditions. We divided the plate interface into small fault patches (10 km × 10 km), whose geometries were set on the basis of the Slab1.0 model (41). Then, we calculated the seafloor displacements caused by the slip deficit rate of 8.3 cm/year (293°N) (37) for each small patch in a uniform elastic half-space (30) as FL. This model uses an end-membered assumption that maximally considers the contributions of FL in the whole-observation regions; thus, note that it does not represent actual effects of FL. The VR-DRs with this FL model, assuming full-locking conditions across the entire plate interface (fig. S6), are shown in Fig. 2. Furthermore, we also examined the VR-DRs for the case of local FL, that is, the shallow and deep portion of Segment-1 (fig. S3), and the PRA except the VLSA (fig. S4).


Supplementary material for this article is available at

fig. S1. GPS-A observation site distribution.

fig. S2. Additional coseismic slip patches for VR calculation and postseismic DRs calculated by the layered VR structure model.

fig. S3. Comparison of the observed DRs and VR-DRs with local FL in Segment-1.

fig. S4. Comparison of the observed DRs and VR-DRs with local FL in PRA.

fig. S5. Time series of the postseismic horizontal displacements at GPS-A observation sites.

fig. S6. Fault patches for calculating the FL effects assuming full-locking conditions.

table S1. Estimated array positions at GPS-A observation sites.

table S2. Coseismic steps caused by the aftershocks on 7 December 2012.

table S3. Estimated postseismic DRs relative to the North American plate.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank O. L. Colombo at the NASA Goddard Space Flight Center for providing us with the kinematic GPS software “IT.” The figures were generated using Generic Mapping Tools software. A part of the observation data was obtained using the Research Vessel Shinsei Maru, supported by the Collaborative Research Program of Atmosphere and Ocean Research Institute, University of Tokyo. Funding: This research was supported by the Ministry of Education, Culture, Sports, Science and Technology, Japan, in the project for “Development of GPS/Acoustic Technique” and by the Japan Science and Technology Agency with the Council for Science, Technology and Innovation and the Cross-ministerial Strategic Innovation Promotion Program “Enhancement of societal resiliency against natural disasters.” Author contributions: F.T. participated in the GPS-A observation cruises and the positioning analysis and wrote the manuscript. M.K., Y.O., T.I., and R.H. participated in the GPS-A observation cruises and the discussion of the results. T.I. analyzed the aftershocks effect. Y.O. participated in obtaining the onshore GPS observation data. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors. We use the calculation results of VR following the 2011 Tohoku earthquake, which were provided by T. Sun, University of Victoria. The daily coordinates of the GEONET sites were provided by the Geospatial Information Authority of Japan.

Stay Connected to Science Advances

Navigate This Article