Monitoring of the spatio-temporal change in the interplate coupling at northeastern Japan subduction zone based on the spatial gradients of surface velocity field

Monitoring of the spatio-temporal change in the interplate coupling at northeastern Japan... Summary A monitoring method to grasp the spatio-temporal change in the interplate coupling in a subduction zone based on the spatial gradients of surface displacement rate fields is proposed. I estimated the spatio-temporal change in the interplate coupling along the plate boundary in northeastern (NE) Japan by applying the proposed method to the surface displacement rates based on global positioning system observations. The gradient of the surface velocities is calculated in each swath configured along the direction normal to the Japan Trench for time windows such as 0.5, 1, 2, 3 and 5 yr being shifted by one week during the period of 1997–2016. The gradient of the horizontal velocities is negative and has a large magnitude when the interplate coupling at the shallow part (less than approximately 50 km in depth) beneath the profile is strong, and the sign of the gradient of the vertical velocity is sensitive to the existence of the coupling at the deep part (greater than approximately 50 km in depth). The trench-parallel variation of the spatial gradients of a displacement rate field clearly corresponds to the trench-parallel variation of the amplitude of the interplate coupling on the plate interface, as well as the rupture areas of previous interplate earthquakes. Temporal changes in the trench-parallel variation of the spatial gradient of the displacement rate correspond to the strengthening or weakening of the interplate coupling. We can monitor the temporal change in the interplate coupling state by calculating the spatial gradients of the surface displacement rate field to some extent without performing inversion analyses with applying certain constraint conditions that sometimes cause over- and/or underestimation at areas of limited spatial resolution far from the observation network. The results of the calculation confirm known interplate events in the NE Japan subduction zone, such as the post-seismic slip of the 2003 M8.0 Tokachi-oki and 2005 M7.2 Miyagi-oki earthquakes and the recovery of the interplate coupling around the rupture area of the 1994 M7.6 Sanriku-Haruka-oki earthquake. The results also indicate the semi-periodic occurrence of slow slip events and the expansion of the area of slow slip events before the 2011 Tohoku-oki earthquake (M9.0) approaching the hypocentre of the Tohoku-oki earthquake. Satellite geodesy, Seismic cycle, Transient deformation, Time-series analysis 1 INTRODUCTION Many large (M > 8) megathrust earthquakes occur in plate subduction zones where oceanic plates are subducting beneath continental plates. Such earthquakes rupture plate interface faults to release the strain energy due to the coupling between the continental and subducting plates during the interseismic period. Coseismic fault slip and interseismic fault coupling generate crustal deformation as well as quasi-static fault slip and transient processes such as visco-elastic relaxation and poro-elastic rebound during post-seismic periods. In order to consider the crustal deformation due to the interplate coupling during the interseismic period, Savage (1983) introduced the ‘backslip model’, in which the interplate coupling is expressed by two overlapping virtual components, namely, the stable sliding between the continental and subducting plates and the normal faulting at the areas where full or partial locking between the plates actually occurs. Plate subduction with stable sliding is regarded as producing no crustal deformation, with the exception of steady plate motion, while the normal faulting on the plate interface fault yields deformation in both plates. Following the proposal of the backslip model (Savage 1983), a number of studies have applied this idea to the estimation of the distribution of interplate coupling patches in subduction zones based on surface displacement fields (e.g. Hyndman & Wang 1995; Norabuena et al.1998; Lundgren et al.1999; Wallace et al.2004). A number of studies have examined the distribution of interplate coupling between the subducting Pacific and overriding continental plates in the northeast (NE) Japan subduction zone (Fig. 1) based on the surface displacement rate field by means of geodetic observations [mainly global positioning system (GPS) observations] (e.g. El-Fiky & Kato 1999; Mazzotti et al.2000; Ito et al.2000; Nishimura et al.2004; Suwa et al.2006; Hashimoto et al.2009; Loveless & Meade 2010). In such studies, the degree of coupling, which has generally been expressed in terms of the backslip (or slip deficit) rate at the plate interface, is estimated by applying linear inversion methods while imposing some constraint conditions, such as the smoothness of the spatial distribution of the backslip, the non-negativity of the backslip, Dirichlet boundary condition and the a priori distribution of the backslip (e.g. Yabuki & Matsu’ura 1992; Mosegaard & Tarantola 1995; Matsu’ura et al.2007; Meade & Loveless 2009). These constraint conditions prevent the solution from becoming unstable, because it is difficult to resolve the backslip distribution at the plate interface far off the coast based only on terrestrial observations. Figure 1. View largeDownload slide Map of the northeast Japan subduction zone. GPS stations are represented by circles, and the blue circles indicate stations used in this study. Note that this map shows the GPS stations for which the available displacement time-series are longer than one year during the entire analysis period. Examples of time-series analysis is presented for two stations, Oshika and Erimo, in Fig. 2. The brown dashed contours indicate the depth of the subducting plate interface (Nakajima & Hasegawa 2006; Nakajima et al.2009; Kita et al.2010). The grey line is the trace of Japan Trench. The intervals of the thin and thick contours are 10 and 50 km, respectively. The black thin lines indicate the triangular subfault used in the synthetic test. The colour contours indicate the rupture area of past large interplate earthquakes. Three events (M7.5, M7.3, and M7.4) in 1938 off Fukushima and Ibaraki Prefectures [Murotani (2003), green, purple and navy], four events off Miyagi Prefecture in 1936 [Yamanaka & Kikuchi (2004), M7.5, maroon], 1978 [Yamanaka & Kikuchi (2004), M7.4, magenta], 1981 [Yamanaka & Kikuchi (2004), M7.1, lime] and 2003 [Yamanaka (2003), M6.8, cyan], the Tokachi-oki earthquakes in 1968 [Yamanaka & Kikuchi (2004), M7.9, blue] and 2003 [Yamanaka & Kikuchi (2003), M8.0, teal], the Nemuro-oki earthquake in 1973 [Yamanaka (2005), M7.4, red] and the 2011 Off the Pacific Coast of Tohoku Earthquake [Iinuma et al. (2012), M9.0, black] are plotted. The names of prefectures used in the text are also shown. Figure 1. View largeDownload slide Map of the northeast Japan subduction zone. GPS stations are represented by circles, and the blue circles indicate stations used in this study. Note that this map shows the GPS stations for which the available displacement time-series are longer than one year during the entire analysis period. Examples of time-series analysis is presented for two stations, Oshika and Erimo, in Fig. 2. The brown dashed contours indicate the depth of the subducting plate interface (Nakajima & Hasegawa 2006; Nakajima et al.2009; Kita et al.2010). The grey line is the trace of Japan Trench. The intervals of the thin and thick contours are 10 and 50 km, respectively. The black thin lines indicate the triangular subfault used in the synthetic test. The colour contours indicate the rupture area of past large interplate earthquakes. Three events (M7.5, M7.3, and M7.4) in 1938 off Fukushima and Ibaraki Prefectures [Murotani (2003), green, purple and navy], four events off Miyagi Prefecture in 1936 [Yamanaka & Kikuchi (2004), M7.5, maroon], 1978 [Yamanaka & Kikuchi (2004), M7.4, magenta], 1981 [Yamanaka & Kikuchi (2004), M7.1, lime] and 2003 [Yamanaka (2003), M6.8, cyan], the Tokachi-oki earthquakes in 1968 [Yamanaka & Kikuchi (2004), M7.9, blue] and 2003 [Yamanaka & Kikuchi (2003), M8.0, teal], the Nemuro-oki earthquake in 1973 [Yamanaka (2005), M7.4, red] and the 2011 Off the Pacific Coast of Tohoku Earthquake [Iinuma et al. (2012), M9.0, black] are plotted. The names of prefectures used in the text are also shown. However, the interplate coupling distributions estimated by previous studies are not unique. The differences between the results may be caused by the differences of the analysis period, the geometrical model of the plate interface, the weight of the vertical displacement and the applied constraint conditions. In particular, the constraint conditions produce artefacts at the shallow plate interface, where the spatial resolution is poor in order to compensate for the errors of observation, the plate geometry and the Green’s function (because of inelastic, non-linear and/or anisotropic effects). Moreover, if we adopt the non-negativity constraint, it messes up the estimated backslip distribution when aseismic slip having a rate larger than the plate convergence rate occurs at some portion of the plate interface, because it requires negative backslip. Therefore, it is not easy to distinguish whether the interplate coupling state temporally changes or not by comparing the results of the inversion analyses. Methods to monitor the changes in the interplate coupling on the subduction plate interface based on the displacement rate field itself may shed light on the temporal variation of the coupling state. If the method can be automatically applied to the real-time geodetic data, it will play an important role in developing the operational monitoring system of the interplate coupling that is essential to predict megathrust earthquakes. Meanwhile, some studies have applied surface displacement rate field based on GPS observations to detect transient deformations in and around plate boundary zones (e.g. Ohtani et al.2010; Holt & Shcherbenko 2013; Crowell et al.2016). Ohtani et al. (2010) and Holt & Shcherbenko (2013) focused to the horizontal strain rate field to find transient deformation events and determine their locations, while Crowell et al. (2016) analysed 3-D displacement time-series at station by station to detect episodic tremor and slip events automatically. However, Ohtani et al. (2010) and Holt & Shcherbenko (2013) did not involve vertical displacement rates in their analysis even though Ohtani et al. (2010) mentioned the extensibility of their method to include the vertical component. In this study, I propose a method for monitoring the spatio-temporal variation of the interplate coupling based on calculating the spatial gradients of the horizontal and vertical surface displacement rate fields. Most of the anomalies detected in the displacement rate gradients correspond to previously reported slow slip events on the plate interface, and the remaining events are unidentified events. It is not possible to quantitatively estimate changes in interplate coupling by applying the monitoring method, although it is possible to grasp the spatio-temporal variation of the change in the interplate coupling in a qualitative manner and to constrain the trench-parallel widths of the perturbed areas and the durations of the events to some extent. 2 DATA AND ANALYSIS 2.1 Displacement time-series and displacement rates In this subsection, I explain the procedure of data processing for the displacement time-series in order to deduce the displacement rate fields. This curve fitting procedure does not discriminate the tectonic deformation signals from non-tectonic ones, but remove random noise due to the monument instabilities, meteorological effects and so on (e.g. Agnew 1992; Williams et al.2004) and steps due to earthquakes and station maintenance. Note that the effect of the spatially correlated noise due to the errors of the reference frame and the satellites’ orbits are reduced by taking spatial gradient of the displacement rate field as described in Section 2.2. Daily site coordinate time-series data obtained from the GPS observation network of the Geospatial Information Authority of Japan (F3 solution, Nakagawa et al.2009) were used to estimate the displacement rate at each GPS station (Fig. 1). The displacement rate is extracted from the site coordinate time-series by means of least-squares fitting for a mathematical function that consists of the long-term linear trend, annual and biannual trigonometric curves, and steps due to earthquakes and station maintenance, such as antenna and receiver replacement. The displacement time-series at the ith site is then expressed as:   \begin{eqnarray} \boldsymbol {u}^i ( t ) & = & \boldsymbol {u}_0^i + \sum _{j=1}^{J-1} \boldsymbol {a}_j^i \left( T_j -T_{j-1} \right) + \boldsymbol {a}_J^i \left( t - T_{J-1} \right)+ \boldsymbol {b}^i \cos \left( \frac{2 \pi t}{\tau } \right) +\boldsymbol {c}^i \sin \left( \frac{2 \pi t}{\tau } \right) + \boldsymbol {d}^i \cos \left( \frac{4 \pi t}{\tau } \right) + \boldsymbol {e}^i \sin \left( \frac{4 \pi t}{\tau } \right)\nonumber \\ &&+\, \sum _{k=1}^K \boldsymbol {f}_k^i H \left( t_k - t \right) \qquad \left( T_{J-1} < t \le T_J \right)\!, \end{eqnarray} (1)where $$\boldsymbol {u}(t)$$ is a displacement vector composed of east–west, north–south and up–down components at date t; Tj is the ending date of the jth sub-time-window in which the long-term linear trend is a constant, where T0 is the starting date of the time-window (refer to the following paragraphs for the time-window and sub-time-windows); tk is the date at which the earthquakes that were large enough to be considered occurred (see Appendix A) or station maintenance was performed; τ (=365.2422) is the number of days in one year; H(t) is the Heaviside step function; and $$\boldsymbol {a}, \boldsymbol {b}, \boldsymbol {c}, \boldsymbol {d}, \boldsymbol {e}, \boldsymbol {f}$$ and $$\boldsymbol {u}_0$$ are the coefficient vectors of the corresponding terms and a constant offset (initial value) vector. Because F3 solutions are usually updated on every Monday and are available since 1996 March 22, 5-yr time-windows ending every Monday from 1997 March 24 to 2016 December 26 were assigned. Before 2001 March 26, all daily coordinate data were applied although the time-window became shorter than 5 yr. For the long-term linear trends, two sets of sub-time-windows are configured to divide the 5-yr time-window into (i) five successive 1-yr sub-time-windows and (ii) two successive 4.5- and 0.5-yr sub-time-windows. The displacement values at the junction of the sub-time-windows are continuous. Then, the linear trend estimated for the latest 1- and 0.5-yr sub-time-windows were regarded as the average displacement rate during the last 1- and 0.5-yr sub-time-windows at the ending date of each 5-yr time window. Note that displacement time-series data during the first four 1- and 4.5-yr sub-time-windows are used for the estimation of annual and biannual trigonometric curves and steps but are not directly used for the average displacement rates at the last sub-time-windows. The time-window is able to be shortened from 5 yr to an arbitrary length as long as we can estimate the annual and biannual variations and steps robustly. It may be possible to optimize the lengths of the time-window and the sub-time-windows in an objective manner based on the observed time-series data, but it is beyond the scope of this paper. The average displacement rate during (iii) the entire 5-yr time-window was also estimated with no sub-time-window and while performing least-squares fitting in the same manner. Fig. 2(a) shows an example of the fitting. In this example, (A) when the five 1-yr sub-time-windows are configured, T0, T1, T2, T3, T4 and T5 in eq. (1) are assigned as 2006 March 07, 2007 March 07, 2008 March 07, 2009 March 07, 2010 March 07 and 2011 March 07, respectively, (B) when 4.5- and 0.5-yr sub-time-windows are configured, T0, T1 and T2 in eq. (1) are assigned as 2006 March 07, 2010 September 07 and 2011 March 07, (C) when no sub-time-window is configured, T0 and T1 in eq. (1) are assigned as 2006 March 07 and 2011 March 07, respectively. Then, t1, t2, t3, t4, t5, t6 and t7, which are the dates of moderate earthquakes and are assigned as 2006 November 15, 2007 January 13, 2008 May 08, June 14, July 19 and 25 and 2010 March 14, respectively, are independent of the sub-time-window setting. Finally, $$\boldsymbol {a}_5, \boldsymbol {a}_2$$ and $$\boldsymbol {a}_1$$ are regarded as the average displacement rates for the last 1- and 0.5-yr sub-time-windows and 5-yr time-window on 2011 March 07, just before the 2011 Tohoku-oki earthquake (M 9) on 2011 March 11. Figure 2. View largeDownload slide Examples of the results of the time-series analysis. (a) Case in which the 5-yr time-window from 2006 March 07 to 2011 March 07 is considered. The EW (top), NS (middle) and UD (bottom) components of the displacement time-series at site 0550 at Oshika are shown. The red dots indicate raw data. The green, blue and orange dots indicate the synthesized displacements obtained using a 5-yr time-window with the offsets (−5 mm for the EW and NS components and −10 mm for the UD component) five 1-yr sub-time-windows, and 4.5- and 0.5-yr sub-time-windows to estimate the linear trends with the offsets (5 mm for the EW and NS components and 10 mm for the UD component). Annual and semi-annual trigonometric curves and steps due to the earthquakes, for which the dates, magnitudes and locations, are exhibited, are also included in the synthesized displacements. Estimated linear trends for time-windows and sub-time-windows are indicted in corresponding colours. (b) Case in which the 6-yr time-window from 1997 September 22 to 2003 September 22 is considered. The EW (top), NS (middle) and UD (bottom) components of the displacement time-series at site 0019 at Erimo are shown. The red dots indicate raw data. The green and blue dots indicate the synthesized displacements using two 3-yr sub-time-windows with the offset (10 mm for all components) and three 2-yr sub-time-windows to estimate the linear trends. Annual and semi-annual trigonometric curves and steps due to earthquakes and site maintenance, for which the dates, magnitudes and locations are exhibited, are also included in the synthesized displacements. Estimated linear trends for sub-time-windows are indicted in corresponding colours. Figure 2. View largeDownload slide Examples of the results of the time-series analysis. (a) Case in which the 5-yr time-window from 2006 March 07 to 2011 March 07 is considered. The EW (top), NS (middle) and UD (bottom) components of the displacement time-series at site 0550 at Oshika are shown. The red dots indicate raw data. The green, blue and orange dots indicate the synthesized displacements obtained using a 5-yr time-window with the offsets (−5 mm for the EW and NS components and −10 mm for the UD component) five 1-yr sub-time-windows, and 4.5- and 0.5-yr sub-time-windows to estimate the linear trends with the offsets (5 mm for the EW and NS components and 10 mm for the UD component). Annual and semi-annual trigonometric curves and steps due to the earthquakes, for which the dates, magnitudes and locations, are exhibited, are also included in the synthesized displacements. Estimated linear trends for time-windows and sub-time-windows are indicted in corresponding colours. (b) Case in which the 6-yr time-window from 1997 September 22 to 2003 September 22 is considered. The EW (top), NS (middle) and UD (bottom) components of the displacement time-series at site 0019 at Erimo are shown. The red dots indicate raw data. The green and blue dots indicate the synthesized displacements using two 3-yr sub-time-windows with the offset (10 mm for all components) and three 2-yr sub-time-windows to estimate the linear trends. Annual and semi-annual trigonometric curves and steps due to earthquakes and site maintenance, for which the dates, magnitudes and locations are exhibited, are also included in the synthesized displacements. Estimated linear trends for sub-time-windows are indicted in corresponding colours. Similarly, I assigned 6-yr time-windows ending every Monday from 1997 March 24 to 2016 December 26 to calculate the displacement rates for the last 2- and 3-yr sub-time-windows at the ending dates. Each 6-yr time-window is divided into (D) two successive 3-yr sub-time-windows and (E) three successive 2-yr sub-time-windows, and the displacement values at the junction of the sub-time-windows are forced to be continuous. The linear trends are estimated by least-squares fitting, and the trends for the latest sub-time-windows are regarded as the average displacement rates during the last 3- and 2-yr sub-time-windows at the ending date of each 6-yr time-window. Fig. 2(b) shows an example of the fitting for 6-yr time-window. In this example, (D) when two 3-yr sub-time-windows are configured, T0, T1 and T2 in eq. (1) are assigned as 1997, September 22, 2000 September 22 and 2003 September 22, respectively, (E) when three 2-yr sub-time-windows are configured, T0, T1, T2 and T3 in eq. (1) are assigned as 1997 September 22, 1999 September 22, 2001 September 22 and 2011 September 22, respectively. Moreover, t1, t2 and t2, which are the date of an M 7.0 earthquake on 2000 January 28 and the dates of antenna replacements on 2001 December 07 and 2002 November 20, are independent of the sub-time-window setting. Finally, $$\boldsymbol {a}_2$$ and $$\boldsymbol {a}_3$$ are regarded as the average displacement rates for the last 3- and 2-yr sub-time-windows on 2003 September 22, just before the 2003 Tokachi-oki (M 8) earthquake on 2003 September 23. Applying these time-series fitting procedures to coordinate time-series of all GPS stations, we can obtain weekly surface displacement rate fields for the last 0.5-, 1-, 2- and 3-yr sub-time-windows, and 5-yr time-window. Naturally, the displacement rate fields for the longer windows are more robust than those for the shorter windows. Figs 3 and 4, and Fig. S1 in the Supporting Information show the examples of estimated displacement rate fields. The horizontal displacement rate field for the last 0.5-yr sub-time-window (Fig. S1, Supporting Information) shows many outlier vectors in directions different from the general trend. Meanwhile, the vertical displacement rate fields for the last 0.5- and 1-yr sub-time-windows include short-wavelength spatial variations with large amplitudes. For instance, the displacement rate at a site at about 40 km in the distance from E140° in the profiles shown in Fig. S1 in the Supporting Information has difference from that at the neighbourhood sites as large as 1σ error (∼10 mm yr−1), while there are no such outlier in the vertical displacement rate profiles for the last 2- and 3-yr sub-time-windows. Therefore, the displacement rate fields are not robustly estimated for the last 0.5-yr sub-time-window with respect to the horizontal component and for the last 0.5- and 1-yr sub-time-windows with respect to the vertical component. Figure 3. View largeDownload slide Examples of the displacement rate profiles and spatial gradients for the Tohoku region. (a) Case in which the centreline of the swath runs through at N38.7°, E140.0° is considered. The displacement rates at sites for which the distance to the swath centre does not exceed 30 km (red vectors) are chosen in order to calculate the spatial gradient. The horizontal component of the displacement rate field for a 5-yr time-window from 2006 March 07 to 2011 March 07 is shown. The bottom panel shows the profile of the horizontal displacement rates along the direction of the swath (N105°E) with respect to the distance from the trench. (b) The vertical components for the case described in (a) are plotted. Figure 3. View largeDownload slide Examples of the displacement rate profiles and spatial gradients for the Tohoku region. (a) Case in which the centreline of the swath runs through at N38.7°, E140.0° is considered. The displacement rates at sites for which the distance to the swath centre does not exceed 30 km (red vectors) are chosen in order to calculate the spatial gradient. The horizontal component of the displacement rate field for a 5-yr time-window from 2006 March 07 to 2011 March 07 is shown. The bottom panel shows the profile of the horizontal displacement rates along the direction of the swath (N105°E) with respect to the distance from the trench. (b) The vertical components for the case described in (a) are plotted. Figure 4. View largeDownload slide Examples of the displacement rate profiles and spatial gradients for the Hokkaido region. (a) Case in which the centreline of the swath runs through N43.5°, E140.0° is considered. The displacement rates at sites for which the distance to the swath centre does not exceed 20 km (red vectors) are chosen in order to calculate the spatial gradient. The horizontal component of the displacement rate field for the 2-yr sub-time-window from 2001 September 22 to 2003 September 22 is shown. The bottom panel shows the profile of the horizontal displacement rates along the direction of the swath (N120°E) with respect to the distance from the trench. (b) The vertical components for the case described in (a) are plotted. Figure 4. View largeDownload slide Examples of the displacement rate profiles and spatial gradients for the Hokkaido region. (a) Case in which the centreline of the swath runs through N43.5°, E140.0° is considered. The displacement rates at sites for which the distance to the swath centre does not exceed 20 km (red vectors) are chosen in order to calculate the spatial gradient. The horizontal component of the displacement rate field for the 2-yr sub-time-window from 2001 September 22 to 2003 September 22 is shown. The bottom panel shows the profile of the horizontal displacement rates along the direction of the swath (N120°E) with respect to the distance from the trench. (b) The vertical components for the case described in (a) are plotted. 2.2 Spatial gradient of the displacement rate field Next, the spatial gradients of the displacement rate fields were calculated for each swath. Swaths were configured along the direction near the plate convergence, namely, N105°E for the Tohoku region and N120°E for the Hokkaido region. The centrelines of the swaths were defined as passing the reference points configured at every 0$${^{\circ}_{.}}$$1 of latitude from N34.5 to N42.0 for the Tohoku region and from N43.0 to N46.5 for the Hokkaido region along meridian E140°. The width of the swath was configured as 40, 60, or 80 km. Spatial gradients of the surface horizontal and vertical displacement rates are estimated with linear regression for the displacement rates of GPS stations inside the swaths (See Appendix B for the details of the calculation and error estimation of the displacement rate gradient). The horizontal displacement rates are projected into the direction of the swaths using the east–west and north–south components to calculate the spatial gradients. The displacement rates estimated from time-series shorter than one year are excluded. Stations at longitudes of less than E139.5° are also excluded. An example of displacement rate field and spatial gradients for the last 5-yr time-window on 2011 March 07 is exhibited in Fig. 3. Fig. 3 also shows the displacement rate profile in a swath of which centreline passes through N38.7°, E140.0°. GPS stations for which the distance to the swath centreline is less than 30 km were included in the calculation of the spatial gradient (the width of the swath is as 60 km). Since the period is just before the 2011 Tohoku-oki earthquake, interplate coupling off Miyagi prefecture was strong, as reflected by large negative value of the horizontal displacement rate gradient. The Pacific coast around Miyagi Prefecture was subsiding during the interseismic period (e.g. Kato & Tsumura 1979; El-Fiky & Kato 1999; Aoki & Scholz 2003; Nishimura et al.2004), and the vertical displacement rate gradient in this profile has a large negative value. Fig. 4 shows the example for the last 2-yr sub-time-window on 2003 September 22. Almost the same features are observed in the profile for which the centreline passes though N43.5°, E140.0° just before the 2003 Tokachi-oki earthquake (Figs 4a and b). However, since, in this case, the width of the swath was 40 km, whereas it was 60 km in the previous case, the number of sites in the swath were fewer. Moreover, as the displacement rates were estimated for the last 2-yr sub-time-window, the error for each displacement rate was larger than in the previous case. Therefore, the spatial gradient of the displacement rates was generally less robust than in the case where the swath was wide and/or the displacement rates were estimated for a longer sub-time-window. This means that there is some trade-off between the robustness of the spatial gradient of the displacement rate field and the resolution of the spatial and temporal change in the spatial gradients. By plotting spatial gradients in all swaths, we can easily grasp the spatial change in the displacement rate gradients along the trench-parallel direction. Fig. 5 shows examples of the distributions of the spatial gradient of the displacement rate fields. The width of the swaths and the length of sub-time-window in which the displacement rates were estimated correspond to the cases shown in Figs 3 and 4. Trench-parallel variations of the displacement rate gradients are clearly observed in each panel, which may reflect the spatial variation in the degree of coupling on the subducting plate interface. For example, Fig. 5(a) shows opposite sign in vertical components between Tohoku and Hokkaido, while the horizontal components show negative spatial gradient in both regions. This difference may be caused by the post-seismic deformation associated with the 2003 Tokachi-oki earthquake, and indicates that the afterslip on the deep plate interface continues longer than that at the shallow portion. Meanwhile, the spatial variation of vertical displacement rate gradient in Tohoku district shown in Fig. 5(b) is more heterogeneous than that shown in Fig. 5(a). Some of the vertical displacement rate gradients shown in Fig. 5(b) are not significantly larger than the estimation errors, but the positive gradients around N39° must be caused by the post-seismic deformation associated with an M7.1 earthquake off Miyagi Prefecture on 2003 May 26 (Okada & Hasegawa 2003). Figure 5. View largeDownload slide Examples of the spatial (trench-parallel) variations of the displacement rate gradients. The horizontal and vertical components are represented by red and blue bars, the origins of which correspond to the intersection points of the swath centrelines and the trench axis. The lengths of the bars are proportional to the absolute values of the displacement rate gradients. The scales are shown in the figure. The bars toward the southeast (seaward) indicate that the signs of the gradients are positive, whereas the bars toward the northwest (landward) indicate that the signs are negative. The thicknesses of the colours of bars correspond to the estimation errors. (a) The spatial gradients of the displacement rate field for the 5-yr time-window from 2006 March 07 to 2011 March 07 is shown together with the horizontal component of the displacement rates. [These vectors are identical to those in Fig. 3(a).] As in Figs 3(a) and (b), the width of each swath is 60 km. The spatial gradients of the displacement rate field for the 2-yr sub-time-window from 2001 September 22 to 2003 September 22 is presented along with the vertical component of the displacement rates. [These vectors are identical to those shown in Fig. 4(b).] As in Figs 4(a) and (b), the width of each swath is 40 km. Figure 5. View largeDownload slide Examples of the spatial (trench-parallel) variations of the displacement rate gradients. The horizontal and vertical components are represented by red and blue bars, the origins of which correspond to the intersection points of the swath centrelines and the trench axis. The lengths of the bars are proportional to the absolute values of the displacement rate gradients. The scales are shown in the figure. The bars toward the southeast (seaward) indicate that the signs of the gradients are positive, whereas the bars toward the northwest (landward) indicate that the signs are negative. The thicknesses of the colours of bars correspond to the estimation errors. (a) The spatial gradients of the displacement rate field for the 5-yr time-window from 2006 March 07 to 2011 March 07 is shown together with the horizontal component of the displacement rates. [These vectors are identical to those in Fig. 3(a).] As in Figs 3(a) and (b), the width of each swath is 60 km. The spatial gradients of the displacement rate field for the 2-yr sub-time-window from 2001 September 22 to 2003 September 22 is presented along with the vertical component of the displacement rates. [These vectors are identical to those shown in Fig. 4(b).] As in Figs 4(a) and (b), the width of each swath is 40 km. Fig. 6 shows the distribution of the average estimated error. The figure clearly indicates that the horizontal displacement rate gradients for sub-time-windows not shorter than one year have estimation errors of less than 1 and 2 mm yr−1 100 km−1 in the Hokkaido region and the Tohoku region north to N38.5° (see Fig. S3, Supporting Information for the 1- and 3-yr sub-time-windows). Since the distances between the coastline and longitude line of E139.5° are small in the Tohoku region south to N38.5°, the estimation errors of the horizontal displacement rates are sometimes larger than 2 mm yr−1 100 km−1 for the swaths for which the width is less than 80 km. With respect to the vertical displacement component, the estimation errors are less than 2 and 5 mm yr−1 100 km−1 in the Hokkaido and Tohoku regions, respectively, for sub-time-windows longer than one year. For most swaths having a width larger than 40 km, the errors are less than 2 mm yr−1 100 km−1 in the area north to N38.5° in the Tohoku region. Figure 6. View largeDownload slide Trench-parallel distribution of the average estimation error of the displacement rate gradients. The upper panel shows the vertical component, whereas the lower panel shows the horizontal component. All cases for the width of the swaths are plotted using various styles of lines (see the legends in the panels). The blue and red lines correspond with 0.5- and 2-yr sub-time-windows. The black lines show the cases for 5-yr time-window. Figure 6. View largeDownload slide Trench-parallel distribution of the average estimation error of the displacement rate gradients. The upper panel shows the vertical component, whereas the lower panel shows the horizontal component. All cases for the width of the swaths are plotted using various styles of lines (see the legends in the panels). The blue and red lines correspond with 0.5- and 2-yr sub-time-windows. The black lines show the cases for 5-yr time-window. Next, I performed synthetic tests in order to understand the relationship between the displacement gradients and the interplate coupling distribution and to examine the detectability while taking the estimation errors into account. 3 INTERPLATE COUPLING AND DISPLACEMENT RATE GRADIENTS 3.1 Synthetic tests In order to confirm the relationship between the spatial gradients of the surface displacement rate field and the backslip distribution on the plate interface, I carried out synthetic tests by applying forward calculations based on dislocation models while assuming various backslip distributions. I used the plate interface model devised by Nakajima & Hasegawa (2006), Nakajima et al. (2009) and Kita et al. (2010). The plate interface fault shallower than 100 km in depth is approximated using triangular subfaults, the edges of which are approximately 20 km in length, to calculate the surface displacements due to the dislocation at each subfault using the computer code provided by Meade (2007). These subfaults are divided according to depth into ten 10-km ranges. For instance, the first range is composed of triangular subfaults shallower that 10 km and the jth range include depths of (j − 1) × 10 to j × 10 km (e.g. the 10th range is composed of subfaults at 90–100 km in depth). The surface displacement rate field due to normal-fault type slip on the subfaults that belong to each range with a slip rate of 4 cm yr−1 (corresponding to approximately half the plate convergence rate) or 8 cm yr−1 (approximately equal to the plate convergence rate) was calculated. As a result, approximately 50 per cent, or full locking, occurs on the plate interface. Fig. 7 shows the calculated displacement rate fields for four of the (310 − 1) cases. The total number of cases is obtained from zero, half and full locking for 10 ranges, excluding the case in which zero locking occurs in every range. I chose the cases in which full locking at 30–50 km in depth (Fig. 7a), half locking at shallower than 30 km in depth, and full locking at 30–50 km in depth (Fig. 7b), full locking at shallower than 50 km in depth (Fig. 7c), and half locking at shallower than 30 or 50–90 km in depth, and full locking at 30–50 km in depth (Fig. 7d) are assumed as typical cases. Figure 7. View largeDownload slide Spatial distributions of the displacement rate gradients using the synthetic displacement rate fields calculated from the assumed interplate coupling distributions. (a) Full coupling (8 cm yr−1) on the plate interface at 30–50 km in depth is assumed (green). (b) Half coupling (4 cm yr−1) on the plate interface shallower than 30 km (turquoise) is added to the case shown in (a). (c) Full coupling (8 cm yr−1) on the plate interface shallower than 50 km is assumed. (d) Half coupling (4 cm yr−1) on the plate interface at 50 to 90 km in depth (turquoise) is added to the cases shown in (b). The horizontal displacement rate field is shown in (a) and (b) along with the displacement rate gradient distribution, while vertical displacement rates are plotted in (c) and (d). Figure 7. View largeDownload slide Spatial distributions of the displacement rate gradients using the synthetic displacement rate fields calculated from the assumed interplate coupling distributions. (a) Full coupling (8 cm yr−1) on the plate interface at 30–50 km in depth is assumed (green). (b) Half coupling (4 cm yr−1) on the plate interface shallower than 30 km (turquoise) is added to the case shown in (a). (c) Full coupling (8 cm yr−1) on the plate interface shallower than 50 km is assumed. (d) Half coupling (4 cm yr−1) on the plate interface at 50 to 90 km in depth (turquoise) is added to the cases shown in (b). The horizontal displacement rate field is shown in (a) and (b) along with the displacement rate gradient distribution, while vertical displacement rates are plotted in (c) and (d). The results indicate that full locking on the plate interface at intermediate depth (∼30 to ∼50 km, Fig. 7a) produces significant negative spatial gradients of the horizontal displacement rates, namely, ∼ −6 and ∼ −4 mm yr−1 100 km−1 in Tohoku and Hokkaido regions, respectively, while large negative spatial gradients of the horizontal displacement rates, namely, ∼ −11 and ∼ −7 mm yr−1 100 km−1 in the Tohoku and Hokkaido regions, respectively, are produced by full locking on the plate interface at shallow and intermediate depths (less than ∼50 km, Fig. 7c). Taking into account the fact that the estimation errors for the spatial gradient of the horizontal displacement rate are less than 2 mm yr−1 100 km−1 in most of the analysis regions (Fig. 6), such a drastic change in the interplate coupling state can be monitored by means of the spatial gradients of the horizontal displacement rates. On the other hand, small (0 to half or half to full) changes in the interplate coupling at the shallow (less than 30 km) plate interface may not be well detected by monitoring the spatial gradients of the horizontal displacement rates, because interplate coupling produces a change in the horizontal displacement rates that may be as large as the estimation errors (Figs 7b and 6). The difference between the horizontal displacement rate fields shown in Figs 7(b) and (d) indicates that the change in the coupling state on the deep (greater than 50 km) plate also causes an insignificant change in the spatial gradients of the horizontal displacement rates. The spatial gradients of the vertical displacement rates have almost no dependency on the change in the interplate coupling state on the intermediate and shallow portions (Figs 7a–c). The difference between the vertical displacement rate fields shown in Figs 7(a) and (c) produces the difference between the spatial gradients of the vertical displacement rates smaller than 2 mm yr−1 100 km−1, which is not significantly larger than the estimation error (Fig. 6). Therefore, we can deduce no information about the change in the degree of coupling on the plate interface at intermediate and shallow depths by monitoring the spatial gradients of the vertical displacement rates. However, the spatial gradients of the vertical displacement rates are prominently affected by the interplate coupling at the deep portion. The sign of the spatial gradients of the vertical displacement rates become negative when the interplate coupling occurs at the deep plate interface, namely, when the plate interface is located beneath land, that is, beneath the GPS observation network (e.g. Fig. 7d). This inversion of the spatial gradient of the vertical displacement rate is natural considering that the hinge line is located above the downdip limit of the dip-slip fault (e.g. Segall 2010). Next, I performed another set of the synthetic tests to examine the constrainablitiy of trench-parallel expansion of the area in which the interplate coupling changes based on the change in the trench-parallel distribution of the displacement rate gradients. The surface displacement rate fields were calculated by modifying the interplate coupling distribution shown in Fig. 7(c), namely, full locking on the plate interface shallower than 50 km in depth. Figs 8(a) and (c) show the examples of the cases in which the coupling in certain areas on the plate interface shallower than 30 km in depth is weakened to 4 cm yr−1. The trench-parallel expansion of a weak coupling zone was configured according to the distance from the centreline of a swath as shown in Fig. 8. The interplate coupling at an area within 40 km from the swath centreline that runs through at N39.0°, E140.0° was weakened in the case shown in Fig. 8(a), while the case in which the swath centreline runs through at N44.0°, E140.0° is exhibited in Fig. 8(c). The displacement rate gradient distribution was calculated and the distribution of the difference of the displacement rate gradient was taken by subtracting that due to the assumed interplate coupling distributions with the weak coupling zone (e.g. Figs 8a and c) from that due to the assumed interplate coupling distributions without the weak coupling zone (Fig. 7c). Changing the centreline of the weak coupling zone as passing the reference points configured at every 1° of latitude from N35 to N41.0 and from N44.0 to N46.0 along meridian E140°, the distributions of the difference were obtained as shown in Figs 9(a) and (c). Figure 8. View largeDownload slide Spatial distributions of the displacement rate gradients using the synthetic displacement rate fields calculated from the assumed interplate coupling distributions with perturbations given by half coupling zones to the coupling distribution shown in Fig. 7(c). (a) The coupling on the plate interface shallower than 30 km in depth and within 40 km from the swath centreline that runs through at N39.0°, E140.0° is weakened to 4 cm yr−1. (b) The half coupling at 50–90 km in depth within 40 km from the swath centreline that runs through at N39.0°, E140.0° is added. (c) The coupling on the plate interface shallower than 30 km in depth and within 40 km from the swath centreline that runs through at N44.0°, E140.0° is weakened to 4 cm yr−1. (d) The half coupling at 50–90 km in depth within 40 km from the swath centreline that runs through at N44.0°, E140.0° is added. The horizontal displacement rate field is shown in (a) and (c) along with the displacement rate gradient distribution, while vertical displacement rates are plotted in (b) and (d). Figure 8. View largeDownload slide Spatial distributions of the displacement rate gradients using the synthetic displacement rate fields calculated from the assumed interplate coupling distributions with perturbations given by half coupling zones to the coupling distribution shown in Fig. 7(c). (a) The coupling on the plate interface shallower than 30 km in depth and within 40 km from the swath centreline that runs through at N39.0°, E140.0° is weakened to 4 cm yr−1. (b) The half coupling at 50–90 km in depth within 40 km from the swath centreline that runs through at N39.0°, E140.0° is added. (c) The coupling on the plate interface shallower than 30 km in depth and within 40 km from the swath centreline that runs through at N44.0°, E140.0° is weakened to 4 cm yr−1. (d) The half coupling at 50–90 km in depth within 40 km from the swath centreline that runs through at N44.0°, E140.0° is added. The horizontal displacement rate field is shown in (a) and (c) along with the displacement rate gradient distribution, while vertical displacement rates are plotted in (b) and (d). Figure 9. View largeDownload slide Spatial distributions of the difference between the displacement rate gradients using the synthetic displacement rate fields calculated from assumed interplate coupling distribution shown in Fig. 7(c) and from perturbed coupling distribution (examples are shown in Fig. 8. (a) and (c) The backslip on the plate interface shallower than 30 km and within 40 km from each swath centreline is weakened to 4 cm yr−1. The horizontal axis indicates the latitude of the crossing point of the centreline of the swath and meridian E140°. The legend of each coloured line indicates the latitude of the crossing point of the centreline of the area where the interplate coupling is perturbed and meridian E140°. (b) and (d) Half coupling (4 cm yr−1) on the plate interface at 50–90 km in depth is added. The differences of the horizontal displacement rate gradient are plotted in (a) and (b), while the vertical components are shown in (c) and (d). Figure 9. View largeDownload slide Spatial distributions of the difference between the displacement rate gradients using the synthetic displacement rate fields calculated from assumed interplate coupling distribution shown in Fig. 7(c) and from perturbed coupling distribution (examples are shown in Fig. 8. (a) and (c) The backslip on the plate interface shallower than 30 km and within 40 km from each swath centreline is weakened to 4 cm yr−1. The horizontal axis indicates the latitude of the crossing point of the centreline of the swath and meridian E140°. The legend of each coloured line indicates the latitude of the crossing point of the centreline of the area where the interplate coupling is perturbed and meridian E140°. (b) and (d) Half coupling (4 cm yr−1) on the plate interface at 50–90 km in depth is added. The differences of the horizontal displacement rate gradient are plotted in (a) and (b), while the vertical components are shown in (c) and (d). Fig. 9(a) shows that the difference of the horizontal displacement gradient mostly distributes around the latitudes at which the centreline of the weak coupling zone is configured, even though the amplitude of the difference depends on the location of the weak coupling zone. The amplitude of the difference may be not large enough to be detected taking the estimation error of the horizontal displacement rate gradient that is discussed in the previous section into account. However, because the change in horizontal displacement rate gradient is almost proportional to the change in the coupling on the shallow plate interface, we can detect the change in the interplate coupling at the shallow portion between full coupling and no coupling based on the monitoring of the displacement rate gradient. Because the number of GPS stations are small, the weak coupling zones near the tips of the land do not change the displacement gradient distribution large enough to be detected. However, it is possible to constrain the trench-parallel expansion of an area where the interplate coupling on the shallow plate interface changes to some extent based on the monitoring of the horizontal displacement rate gradient. Note that the change in the coupling state at the shallow plate interface causes no significant difference of the vertical displacement gradient as Fig. 9(c) shows. Fig. S4 in the Supporting Information shows other cases in which the widths of the swath (40, 60 and 80 km) and the width of the week coupling zone (40, 80, 120, 160 and 200 km) are examined. The distribution of the difference sometimes becomes bimodal when the horizontal displacement rate gradient in narrow (40 km in width) swath is calculated for the displacement field synthesized by assuming wide (≥80 km) weak coupling zone. In such case, we should utilize wider swath to calculate the displacement rate gradient, and it is possible to estimate the true centreline of the weak coupling zone. With respect to the change in the coupling state on the plate interface at deep portion, the cases in which the weak coupling zone at 50–90 km in depth is added were examined. Figs 8(b) and (d) show the examples of the assumed interplate coupling distributions. The distributions of the difference of the displacement rate gradients were calculated in the same manner as the tests with respect to the weak coupling zone in the shallow portion (Figs 9b and d). It shows that the changes in the vertical displacement rate gradient are concentrated around added deep coupling zones, while the distribution of the difference with respect to the horizontal component do not well correspond with the change in the interplate coupling. Fig. S5 in the Supporting Information shows other cases as well as Fig. S4 in the Supporting Information. The relationship between the widths of the swath and the weak coupling zone in the deep portion has the same tendency as the relationship between the widths of the swath and the weak coupling zone in the shallow portion as described above. In summary, we can detect the temporal change in the coupling state on the plate interface offshore and onshore in the NE Japan subduction zone to a certain degree by monitoring the temporal changes in the spatial gradients of the horizontal and vertical displacement rates, respectively. And it is possible to constrain the location of the area where the interplate coupling change occurs to some extent. Of course, it is difficult to distinguish the ranges at which the change in the degree of coupling actually occurs among depth ranges of 10–50 km. However, low spatial resolution is inevitable in the trench-normal direction, even if we apply geodetic inversion analysis, as long as only onshore GPS data are used. 4 SPATIO-TEMPORAL CHANGE IN INTERPLATE COUPLING 4.1 Spatio-temporal variation in the velocity gradient Fig. 10 shows the spatio-temporal variations of the horizontal and vertical displacement rate gradients for the last 1-yr sub-time-windows and 5-yr time-windows in the 60-km-wide swaths. (See Fig. S6 in the Supporting Information for the cases of other sub-time-window lengths and swath widths.) After the 2011 Tohoku-oki earthquake, the displacement rate field in the Tohoku region is strongly influenced by the visco-elastic response (e.g. Watanabe et al.2014; Sun et al.2014; Freed et al.2017), and it is difficult to deduce the information with respect to the interplate coupling based on the spatio-temporal changes in the displacement rate gradients. Figure 10. View largeDownload slide Spatio-temporal variations of the displacement rate gradients. Results for (a) the 1-yr sub-time-window and (b) the 5-yr time-window. The horizontal axis indicates the latitudes of the crossing point of the centreline of the swath and meridian E140°. The panels on the left-hand side correspond to the Kanto to Tohoku regions, and the panels on the right-hand side correspond to the Hokkaido region. The vertical axis indicates the final dates of the sub-time-windows and time-windows. The results for the 1-yr sub-time-window and the 5-yr time-window include 1032 Mondays when the sub-time-windows end from 1997 March 24 to 2016 December 26 and 823 Mondays when the time-windows end from 2001 March 26 to 2016 December 26, respectively, but the half of them (every other Monday) are plotted in each panel. The top and bottom panels show the vertical and horizontal displacement rate gradients, respectively. Note that different colour scales are used for the horizontal and vertical components. Figure 10. View largeDownload slide Spatio-temporal variations of the displacement rate gradients. Results for (a) the 1-yr sub-time-window and (b) the 5-yr time-window. The horizontal axis indicates the latitudes of the crossing point of the centreline of the swath and meridian E140°. The panels on the left-hand side correspond to the Kanto to Tohoku regions, and the panels on the right-hand side correspond to the Hokkaido region. The vertical axis indicates the final dates of the sub-time-windows and time-windows. The results for the 1-yr sub-time-window and the 5-yr time-window include 1032 Mondays when the sub-time-windows end from 1997 March 24 to 2016 December 26 and 823 Mondays when the time-windows end from 2001 March 26 to 2016 December 26, respectively, but the half of them (every other Monday) are plotted in each panel. The top and bottom panels show the vertical and horizontal displacement rate gradients, respectively. Note that different colour scales are used for the horizontal and vertical components. The spatio-temporal change in the horizontal displacement rate gradients for 5-yr time windows (Fig. 10) clearly indicates the transient change in the interplate coupling that were reported in previous studies. For instance, The decrease in the negative gradients at N36.5° to N38.5° in 2008–2011 (before the Tohoku-oki earthquake) corresponding to the weak interplate coupling at the shallow and intermediate plate interfaces, as described in the previous section (Ozawa et al.2012; Mavrommatis et al.2014; Yokota & Koketsu 2015). The increase in the negative gradients around N41.0° before the 2011 M 9.0 earthquake representing the recovery process of the interplate coupling at and around the rupture area of the 1994 Sanriku-Haruka-oki earthquake (Mw7.7, Ozawa et al.2007). The dissipation of large negative gradients at N43.5° to N45.0° around 2006 corresponding to the afterslip of the 2003 Tokachi-oki earthquake (M 8.0, Miyazaki et al.2004; Murakami et al.2006). The effect of the afterslip due to the 2005 Miyagi-oki earthquake (M 7.2 on 16 August 2005) can be observed, but the magnitude of the change in the displacement gradient is comparatively smaller than in the cases of the above large and long-term events. Considering the vertical component, negative displacement gradients of up to −5 ± 2 mm yr−1 100 km−1 are consistently observed in the Tohoku region, N37.5° to N40.5°, before 2011, and of up to −3 ± 2 mm yr−1 100 km−1 are consistently observed in the Hokkaido region, N44.2° to N46.2°, before 2005. This implies that the interplate coupling at deep plate interfaces was strong before megathrust earthquakes, such as the 2011 Tohoku-oki and 2003 Tokachi-oki earthquakes. Since the displacement rate fields that are estimated for 5-yr time-windows do not contain localized deformation as shown in Fig. 3, the spatial gradients are robustly estimated, but temporal changes in shorter periods cannot be monitored. The results shown in Figs 10(a) and S6 indicate that there are significant changes in the horizontal displacement rate gradients for periods shorter than five years. This indicates that the interplate coupling on the plate interface at shallow and/or intermediate depths may change depending over the course of a couple of years. 4.2 Periodic change in interplate coupling The spatio-temporal variation in the horizontal displacement rate gradients for a 1-yr sub-time-window (Fig. 10a) reveals that the decrease in the negative gradients at N36.5° to N38.5° did not occur monotonically. The decreases occurred around 1997, 2001, 2005 and 2009, and the northern limit of the fluctuating profiles advanced northward, namely, the decrease in the displacement rate gradient occurred around N37.0° in 1997, around N37.4° in 2001, around N37.6° in 2005 and around N38.5°, in 2009. Note that the very small (≤ −2.5 mm yr−1 100 km−1, coloured blue) displacement rate gradients estimated during 2006 between N38.0° and N39.0° correspond to the afterslip of the 2005 Miyagi-oki earthquake (M7.2) (Miura et al.2006). Fig. 11 shows a time-series of the horizontal displacement rate gradients for 3- and 5-yr sub-time-windows. The quasi-periodic fluctuation events described above are also observed in these time-series in Fig. 11(a). For instance, the increase (decrease of the negative value) of the displacement rate gradient occurred for the profiles between N36.2° and N37.4° during the period of 2001–2002, between N36.8° and N37.8° during the period of 2005–2006 and between N36.8° and N38.4° during the period of 2009–2010. Note that the periods of fluctuation are shifted approximately one year later here, as compared to the discussion of the previous paragraph, because the length of sub-time-window (three years) is longer than one year. On the other hand, we cannot clearly observe the fluctuations during the periods of 2001–2002 and 2005–2006 based on the time-series of the displacement rate gradients of 5-yr time-windows shown in Fig. 11(b). This indicates that the duration of the change in the interplate coupling at shallow and intermediate depths off Fukushima and Ibaraki Prefectures is shorter than five years. Figure 11. View largeDownload slide Time-series of the horizontal displacement rate gradients for the 40-km-wide swaths at N36.0° to N39.0° before the 2011 Tohoku-oki earthquake. The results for the (a) 3-yr and (b) 5-yr time-windows are shown. The data for the swath are plotted with colours corresponding to the profile latitude at E140.0°, as shown at the right-hand side of each panel. The latitude numbers used to label the horizontal lines correspond to the zero displacement rate gradients for each profile. The grey squares indicate fluctuations mentioned in the text. Figure 11. View largeDownload slide Time-series of the horizontal displacement rate gradients for the 40-km-wide swaths at N36.0° to N39.0° before the 2011 Tohoku-oki earthquake. The results for the (a) 3-yr and (b) 5-yr time-windows are shown. The data for the swath are plotted with colours corresponding to the profile latitude at E140.0°, as shown at the right-hand side of each panel. The latitude numbers used to label the horizontal lines correspond to the zero displacement rate gradients for each profile. The grey squares indicate fluctuations mentioned in the text. 5 DISCUSSION As described in Section 2.1, the displacement rate field deduced from the time-series analysis performed in this study includes both tectonic and non-tectonic deformations. It is quite difficult to distinguish and eliminate non-tectonic signals due to causes such as ocean tide loading, atmospheric mass loading and continental water loading from the estimated displacement rate field in a strict manner. Nevertheless, the deformation due to the seasonal water variations are able to be removed by including annual and biannual trigonometric curves in the time-series fitting to some extent, and non-seasonal non-tectonic deformation with a short time-scale should be disappeared by taking the long-term trends. Comparing between the results of different sub-time-windows, swath widths and components may help us to distinguish long-term non-tectonic deformation from the tectonic deformation. Recently, Uchida et al. (2016) revealed the periodic changes in the interplate coupling in the NE Japan subduction zone based primarily on the activity of small repeating earthquakes. They found that slow slip on the plate interface had occurred repeatedly at intervals of 1–6 yr, depending on the location. They also used the horizontal displacement gradient time-series of the 1-yr sub-time-window to estimate the recurrence interval and concluded that the recurrence intervals based on the seismic and geodetic data analyses are mutually corresponding. Based on the small repeating earthquakes, Uchida et al. (2016) estimated the periods of slow slip events that occur off Fukushima Prefecture (corresponding to the swaths between N37° and N38° in this study) at shallow and intermediate depths as 3–5 yr (Uchida et al.2016, Fig. 4). Taking this into account, the dissipation of the interplate coupling in the area south of the main rupture of the 2011 Tohoku-oki earthquake did not progress monotonically, but rather episodically, while expanding northward. This implies that the process of earthquake generation on the subducting plate interface is strongly affected by the background periodic changes in the interplate coupling, and such changes can be monitored not only based on the small repeating earthquakes but also on temporal changes in the displacement rate gradients. The displacement rate distribution can, in principle, be estimated for a time-window shorter than 0.5 yr. However, as shown in Fig. 6 and Fig. S3 in the Supporting Information, since the horizontal displacement rate gradients estimated for a 0.5-yr sub-time-window and vertical displacement rate gradients estimated for 0.5- and 1-yr sub-time-windows have large estimation errors. And as shown in Fig. 10 and Fig. S6 in the Supporting Information large perturbations in a very short period of time, the displacement rates for sub-time-windows shorter than 0.5 yr are not robust. On the other hand, with respect to the width of the swath, 40 km is large enough to robustly estimate displacement rate gradients for most profiles, except those near the ends of Honshu and Hokkaido Islands (Fig. 6 and Fig. S3, Supporting Information). The width of the swaths corresponds to the monitoring resolution of the change in the interplate coupling along the trench-parallel direction. Although the most appropriate width must depend on the density of the GPS observation stations, we can apply the method to other subduction zones where continuous GPS observation networks are sufficiently dense, such as Cascadia and Hikurangi. The density and expansion of the observation network in both trench-normal and trench-parallel directions determine the availability of the proposed method to other subduction zones. It is acceptable to use campaign sites, if we permit that the temporal resolution becomes lower. The graphical expression introduced in this study may make it easy to grasp the spatio-temporal change in the interplate coupling. We cannot estimate the detailed distribution of the area where the interplate coupling changes based only on the change in the displacement rate gradient, especially along the trench-normal direction. However, the recent development of seafloor observation networks equipped with ocean-bottom pressure gauges (OBPGs) such as DONET (Kaneda et al.2015; Kawaguchi et al.2015) and S-net (Uehira et al.2012; Kanazawa 2013) may lead to further applications of the monitoring method based on the displacement rate gradients. We may infer the relative vertical displacement rates in the seafloor observation networks from the OBPG records, even though the sensor drift of OBPGs makes the estimation of the absolute vertical displacement at each OBPG site difficult. 6 CONCLUSIONS A general framework of monitoring method for spatio-temporal change in the interplate coupling based on the spatial gradients of surface displacement rate fields is proposed and applied to NE Japan subduction zone. The results show that we can detect known interplate events in the NE Japan subduction zone, such as the post-seismic slip of the 2003 M8.0 Tokachi-oki and 2005 M7.2 Miyagi-oki earthquakes and the recovery of the interplate coupling around the rupture area of the 1994 M7.6 Sanriku-Haruka-oki earthquake. The results also imply the semi-periodic occurrence of slow slip events and the expansion of the area of slow slip events before the 2011 Tohoku-oki earthquake (M9.0) approaching the hypocentre of the Tohoku-oki earthquake. The proposed method can be applied to other subduction zones and to the operational monitoring after detailed examinations for temporal and spatial estimation error distributions and optimization of the parameters that depend on the site density of the observation network and geometry between the subducting plate interface and the land. Acknowledgements Discussions with Prof. Akira Hasegawa, Prof. Toru Matsuzawa, Prof. Ryota Hino and Dr. Naoki Uchida were valuable in the development of the monitoring method introduced in this paper. I thank two anonymous reviewers and the editor, Prof. Duncan Agnew, for their constructive comments that greatly improved this manuscript. I used Generic Mapping Tools devised by Wessel & Smith (1998) to generate the figures in this paper. This research was partly supported by JSPS KAKENHI grant number JP17H05422 and the project ‘Evaluation and disaster prevention research for the coming Tokai, Tonankai and Nankai earthquakes’ of the Ministry of Education, Culture, Sports, Science and Technology of Japan. REFERENCES Agnew D.C., 1992. The time-domain behavior of power-law noises, Geophys. Res. Lett. , 19, 333– 336. https://doi.org/10.1029/91GL02832 Google Scholar CrossRef Search ADS   Aoki Y., Scholz C.H., 2003. Vertical deformation of the Japanese islands, 1996–1999, J. geophys. Res. , 108, 2257, doi:10.1029/2002JB002129. Crowell B.W., Bock Y., Liu Z., 2016. Single-station automated detection of transient deformation in GPS time series with the relative strength index: a case study of Cascadian slow slip, J. geophys. Res. , 121, 9077– 9094. https://doi.org/10.1002/2016JB013542 Google Scholar CrossRef Search ADS   El-Fiky G.S., Kato T., 1999. Interplate coupling in the Tohoku district, Japan, deduced from geodetic data inversion, J. geophys. Res. , 104, 20,361– 20,377. https://doi.org/10.1029/1999JB900202 Google Scholar CrossRef Search ADS   Freed A.M., Hashima A., Becker T.W., Okaya D.A. Sato H., Hatanaka T., 2017. Resolving depth-dependent subduction zone viscosity and afterslip from postseismic displacements following the 2011 Tohoku-oki, Japan earthquake, Earth planet. Sci. Lett. , 459, 279– 290. https://doi.org/10.1016/j.epsl.2016.11.040 Google Scholar CrossRef Search ADS   Hashimoto C., Noda A., Sagiya T., Matsu’ura M., 2009. Interplate seismogenic zones along the Kuril-Japan trench inferred from GPS data inversion, Nat. Geosci. , 2, 141– 144. https://doi.org/10.1038/ngeo421 Google Scholar CrossRef Search ADS   Holt W.E., Shcherbenko G., 2013. Toward a continuous monitoring of the horizontal displacement gradient tensor field in Southern California using cGPS observations from Plate Boundary Observatory (PBO), Seismol. Res. Lett. , 84, 455– 467. https://doi.org/10.1785/0220130004 Google Scholar CrossRef Search ADS   Hyndman R.D., Wang K., 1995. The rupture zone of Cascadia great earthquakes from current deformation and the thermal regime, J. geophys. Res. , 100, 22 133– 22 154. https://doi.org/10.1029/95JB01970 Google Scholar CrossRef Search ADS   Iinuma T. et al.  , 2012. Coseismic slip distribution of the 2011 off the Pacific Coast of Tohoku Earthquake (M9.0) refined by means of seafloor geodetic data, J. geophys. Res. , 117, B07409, doi:10.1029/2012JB009186. https://doi.org/10.1029/2012JB009186 Google Scholar CrossRef Search ADS   Ito T., Yoshioka S., Miyazaki S., 2000. Interplate coupling in northeast Japan deduced from inversion analysis of GPS data, Earth planet. Sci. Lett. , 176, 117– 130. https://doi.org/10.1016/S0012-821X(99)00316-7 Google Scholar CrossRef Search ADS   Kanazawa T., 2013. Japan trench earthquake and tsunami monitoring network of cable-linked 150 ocean bottom observatories and its impact to Earth disaster science, in 2013 IEEE International Underwater Technology Symposium (UT) , pp. 1– 5, IEEE, Tokyo. Kaneda Y. et al.  , 2015. Development and application of an advanced ocean floor network system for megathrust earthquakes and tsunamis, in Seafloor Observatories , pp. 643– 662, eds Favali P. et al.  , Springer Praxis Books. Google Scholar CrossRef Search ADS   Kato T., Tsumura K., 1979. Vertical land movement in Japan as deduced from tidal record (1951–1978) (in Japanese), Bull. Earthq. Res. Inst., Univ. Tokyo , 54, 559– 628. Kawaguchi K., Kaneko S., Nishida T., $ Komine T., 2015. Construction of the DONET real-time seafloor observatory for earthquakes and tsunami monitoring, in Seafloor Observatories , pp. 211– 228, eds Favali P. et al.  , Springer Praxis Books. Google Scholar CrossRef Search ADS   Kikuchi M., 2001. Fault model: Earthquake source processes (in Japanese), in Encyclopedia of Earthquakes , 2nd edn, pp. 259– 284, eds Utsu T., Shima E., Yoshii T., Yamashina K., Asakura Publishing, Tokyo. Kita S., Okada T., Hasegawa A., Nakajima J., Matsuzawa T., 2010, Anomalous deepening of a seismic belt in the upper-plane of the double seismic zone in the Pacific slab beneath the Hokkaido corner: possible evidence for thermal shielding caused by subducted forearc crust materials, Earth planet. Sci. Lett. , 290, 415– 426. https://doi.org/10.1016/j.epsl.2009.12.038 Google Scholar CrossRef Search ADS   Loveless J.P., Meade B.J., 2010. Geodetic imaging of plate motions, slip rates, and partitioning of deformation in Japan, J. geophys. Res. , 115, B02410, doi:10.1029/2008JB006248. https://doi.org/10.1029/2008JB006248 Google Scholar CrossRef Search ADS   Lundgren P., Protti M., Donnellan A., Heflin M., Hernandez E., Jefferson D., 1999. Seismic cycle and plate margin deformation in Costa Rica: GPS observations from 1994 to 1997, J. geophys. Res. , 104, 28 915– 28 926. https://doi.org/10.1029/1999JB900283 Google Scholar CrossRef Search ADS   Matsu’ura M., Noda A., Fukahata Y., 2007. Geodetic data inversion based on Bayesian formulation with direct and indirect prior information, Geophys. J. Int. , 171, 1342– 1351. https://doi.org/10.1111/j.1365-246X.2007.03578.x Google Scholar CrossRef Search ADS   Mazzotti S., Le Pichon X., Henry P., Miyazaki S., 2000. Full interseismic locking of the Nankai and Japan-west Kurile subduction zones: an analysis of uniform elastic strain accumulation in Japan constrained by permanent GPS, J. geophys. Res. , 105, 13 159–13 177. https://doi.org/10.1029/2000JB900060 Mavrommatis A.P., Segall P., Johnson K.M., 2014. A decadal-scale deformation transient prior to the 2011 Mw 9.0 Tohoku-oki earthquake, Geophys. Res. Lett. , 41, 4486– 4494. https://doi.org/10.1002/2014GL060139 Google Scholar CrossRef Search ADS   Meade B.J., 2007, Algorithms for the calculation of exact displacements, strains, and stresses for triangular dislocation elements in a uniform elastic half space, Comput. Geosci. , 33, 1064– 1075. https://doi.org/10.1016/j.cageo.2006.12.003 Google Scholar CrossRef Search ADS   Meade B.J., Loveless J.P., 2009. Block modeling with connected fault-network geometries and a linear elastic coupling estimator in spherical coordinates, Bull. seism. Soc. Am. , 99, 3124– 3139. https://doi.org/10.1785/0120090088 Google Scholar CrossRef Search ADS   Miura S., Iinuma T., Yui S., Uchida N., Sato T., Tachibana K., Hasegawa A., 2006 Co- and post-seismic slip associated with the 2005 Miyagi-oki earthquake (M7.2) as inferred from GPS data, Earth Planets Space , 58, 1567– 1572. https://doi.org/10.1186/BF03352662 Google Scholar CrossRef Search ADS   Miyazaki S., Segall P., Fukuda J., Kato T., 2004. Space time distribution of afterslip following the 2003 Tokachi-oki earthquake: implications for variations in fault zone frictional properties, Geophys. Res. Lett. , 31, L06623, doi:10.1029/2003GL019410. Mosegaard K., Tarantola A., 1995. Monte Carlo sampling of solutions to inverse problems, J. geophys. Res. , 100, 12 431– 12 447. https://doi.org/10.1029/94JB03097 Google Scholar CrossRef Search ADS   Murakami M., Suito H., Ozawa S., Kaidzu M., 2006. Earthquake triggering by migrating slow slip initiated by M8 earthquake along Kuril Trench, Japan, Geophys. Res. Lett. , 33, L09306, doi:10.1029/2006GL025967. https://doi.org/10.1029/2006GL025967 Google Scholar CrossRef Search ADS   Murotani S., 2003. Rupture processes of large Fukushima-Oki earthquakes in 1938, Master’s thesis , University of Tokyo, Tokyo. Nakagawa H.et al.  , 2009. Development and validation of GEONET new analysis strategy (Version 4) (in Japanese), J. Geogr. Surv. Inst. , 118, 1– 8. Nakajima J., Hasegawa A., 2006. Anomalous low-velocity zone and linear alignment of seismicity along it in the subducted Pacific slab beneath Kanto, Japan: reactivation of subducted fracture zone?, Geophys. Res. Lett. , 33, L16309, doi:10.1029/2006GL026773. https://doi.org/10.1029/2006GL026773 Google Scholar CrossRef Search ADS   Nakajima J., Hirose F., Hasegawa A., 2009. Seismotectonics beneath the Tokyo metropolitan area, Japan: effect of slab-slab contact and overlap on seismicity, J. geophys. Res. , 114, B08309, doi:10.1029/2008JB006101. https://doi.org/10.1029/2008JB006101 Google Scholar CrossRef Search ADS   Nishimura T., Hirasawa T., Miyazaki S., Sagiya T., Tada T., Miura S., Tanaka K., 2004. Temporal change of interplate coupling in northeastern Japan during 1995–2002 estimated from continuous GPS observations, Geophys. J. Int. , 157, 901– 916. https://doi.org/10.1111/j.1365-246X.2004.02159.x Google Scholar CrossRef Search ADS   Norabuena E., Leffler-Griffin L., Mao A., Dixon T., Stein S., Sacks I.S., Ocola L., Ellis M., 1998. Space geodetic observations of Nazca-South America convergence across the Central Andes, Science , 279, 358– 362. https://doi.org/10.1126/science.279.5349.358 Google Scholar CrossRef Search ADS PubMed  Ohtani R., McGuire J.J., Segall P., 2010. Network strain filter: a new tool for monitoring and detecting transient deformation signals in GPS arrays, J. geophys. Res. , 115, B12418, doi:10.1029/2010JB007442. https://doi.org/10.1029/2010JB007442 Google Scholar CrossRef Search ADS   Okada Y., 1992. Internal deformation due to shear and tensile faults in a half-space, Bull. seism. Soc. Am. , 82, 1018– 1040. Okada T., Hasegawa A., 2003. The M7.1 May 26, 2003 off-shore Miyagi Prefecture Earthquake in northeast Japan: source process and aftershock distribution of an intra-slab event, Earth Planets Space , 55, 731– 739. https://doi.org/10.1186/BF03352482 Google Scholar CrossRef Search ADS   Ozawa S., Suito H., Nishimura T., Tobita M., Munekane H., 2007. Possibility of recovery of slip deficit rate between the North American plate and the Pacific plate off Sanriku, northeast Japan, Geophys. Res. Lett. , 34, L20308, doi:10.1029/2007GL030477. https://doi.org/10.1029/2007GL030477 Google Scholar CrossRef Search ADS   Ozawa S., Nishimura T., Munekane H., Suito H., Kobayashi T., Tobita M., Imakiire T., 2012. Preceding, coseismic, and postseismic slips of the 2011 Tohoku earthquake, Japan, J. geophys. Res. , 117, B07404, doi:10.1029/2011JB009120. https://doi.org/10.1029/2011JB009120 Google Scholar CrossRef Search ADS   Savage J.C., 1983. A dislocation model of strain accumulation and release at a subduction zone, J. geophys. Res. , 88, 4984– 4996. https://doi.org/10.1029/JB088iB06p04984 Google Scholar CrossRef Search ADS   Segall P., 2010. Earthquake and Volcano Deformation , pp. 51– 91, Princeton University Press. Sun T.et al.  , 2014. Prevalence of viscoelastic relaxation after the 2011 Tohoku-oki earthquake, Nature , 514, 84– 87. https://doi.org/10.1038/nature13778 Google Scholar CrossRef Search ADS PubMed  Suwa Y., Miura S., Hasegawa A., Sato T., Tachibana K., 2006. Interplate coupling beneath NE Japan inferred from three-dimensional displacement field, J. geophys. Res. , 111, B04402, doi:10.1029/2004JB003203. Uchida N., Iinuma T., Nadeau R.M., Bürgmann R., Hino R, 2016. Periodic slow slip triggers megathrust zone earthquakes in northeastern Japan, Science , 351, 488– 492. https://doi.org/10.1126/science.aad3108 Google Scholar CrossRef Search ADS PubMed  Uehira K. et al.  , 2012. Ocean bottom seismic and tsunami network along the Japan trench, in Abstract OS41C-1736 presented at 2012 Fall Meeting , AGU, San Francisco, California, 3–7 December. Wallace L.M., Beavan J., McCaffrey R., Darby D., 2004. Subduction zone coupling and tectonic block rotations in the North Island, New Zealand, J. geophys. Res. , 109, B12406, doi:10.1029/2004JB003241. https://doi.org/10.1029/2004JB003241 Google Scholar CrossRef Search ADS   Watanabe S., Sato M., Fujita M., Ishikawa T., Yokota Y., Ujihara N., Asada A., 2014. Evidence of viscoelastic deformation following the 2011 Tohoku-Oki earthquake revealed from seafloor geodetic observation, Geophys. Res. Lett. , 41, 5789– 5796. https://doi.org/10.1002/2014GL061134 Google Scholar CrossRef Search ADS   Wessel P., Smith W.H. F., 1998. New, improved version of Generic Mapping Tools released, EOS, Trans. Am. geophys. Un. , 79, 579, doi:10.1029/98EO00426. https://doi.org/10.1029/98EO00426 Google Scholar CrossRef Search ADS   Williams S.D.P., Bock Y., Fang P., Jamason P., Nikolaidis R.M., Prawirodirdjo L., Miller M., Johnson D.J., 2004. Error analysis of continuous GPS position time series, J. geophys. Res. , 109, B03412, doi:10.1029/2003JB002741. Yabuki T., Matsu’ura M., 1992. Geodetic data inversion using a Bayesian information criterion for spatial distribution of fault slip, Geophys. J. Int. , 109, 363– 375. https://doi.org/10.1111/j.1365-246X.1992.tb00102.x Google Scholar CrossRef Search ADS   Yamanaka Y., 2003. Off Fukushima-ken earthquake, 31 October 2003, EIC Seismological Notes , 141. Available at http://www.eri.u-tokyo.ac.jp/sanchu/Seismo_Note/EIC_News/031031.html, last accessed 31 May 2017. Yamanaka Y., Kikuchi M., 2003. Source process of the recurrent Tokachi-oki earthquake on September 26, 2003, inferred from teleseismic body waves, Earth Planets Space , 55, e21– e24. https://doi.org/10.1186/BF03352479 Google Scholar CrossRef Search ADS   Yamanaka Y., Kikuchi M., 2004. Asperity map along the subduction zone in northeastern Japan inferred from regional seismic data, J. geophys. Res. , 109, B07307, doi:10.1029/2003JB002683. https://doi.org/10.1029/2003JB002683 Google Scholar CrossRef Search ADS   Yamanaka Y., 2005. Asperity Map along the subduction zone in Hokkaido region (in Japanese), in Abstract A055 Presented at 2005 Fall Meeting , Seism. Soc. Jpn., Sapporo, Japan, 19–21 October. Yokota Y., Koketsu K., 2015. A very long-term transient event preceding the 2011 Tohoku earthquake, Nature Commun. , 6, 5934, doi:10.1038/ncomms6934. https://doi.org/10.1038/ncomms6934 Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1. Examples of the displacement rate field for 0.5-year, 1-year, 2-year and 3-year sub-timewindows from 7 March 2011 (see Figure 3 for 5-year time-window). The displacement rate profiles in the swath for which the centerline runs through N38.7°, E140.0° are also shown. Figure S2. All profiles of an example displacement rate field for a 5-year time-window from 7 March 2006 to 7 March 2011. The width of the swaths is configured to be 60 km. The red vectors indicate the horizontal displacement rates at the sites in the swaths. The bottom panel shows the profile of the horizontal and vertical displacement rates along the direction of the swath with respect to the distance from the trench. Figure S3. Trench-parallel distribution of the average estimation error of the displacement rate gradients. The upper panel shows the vertical component, whereas the lower panel shows the horizontal component. All cases for the width of the swaths are plotted using various styles of lines (see the legends in the panels). Orange and green lines correspond with 1-year and 3-year sub -timewindows. Figure S4. Spatial distributions of the difference between displacement rate gradients using the synthetic displacement rate fields calculated from the assumed interplate coupling distribution shown in Figure 7(c) and the interplate coupling distribution with perturbation on plate interface shallower than 30 km in depth. The widths of the swaths (D) and perturbed regions (W) are implemented in the panels. Figure S5. Spatial distributions of the difference between displacement rate gradients using the synthetic displacement rate fields calculated from the assumed interplate coupling distribution shown in Figure 7(c) and the interplate coupling distribution with perturbation on plate interface at 50 to 90 km in depth. The widths of the swaths (D) and perturbed regions (W) are implemented in the panels. Figure S6. Spatio-temporal variations of the displacement rate gradients for all cases. The horizontal axis indicates the latitudes of the crossing point of the centerline of the profile region and meridian E140°. The panels on the left-hand side correspond to the Kanto to Tohoku regions, and the panels on the right-hand side correspond with the Hokkaido region. The vertical axis indicates the final dates of the sub-time-windows and the time-windows. The top and bottom panels represent the vertical and horizontal displacement rate gradients. Note that different color scales are used for horizontal and vertical components. Figure S7. Time series of the horizontal displacement rate gradients at N36.0° to N39.0° before the 2011 Tohoku-oki earthquake for all cases. The data for each profile region are plotted with a color corresponding to the profile latitude at E140.0°, as shown at the right-hand side of each panel. The horizontal line identified by the latitude number corresponds with the zero displacement rate gradient for each profile. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. APPENDIX A: SELECTION OF EARTHQUAKES IN THE DISPLACEMENT TIME-SERIES ANALYSIS The earthquake catalogue of the Japan Meteorological Agency for the period of 1996–2016 was used to configure the dates for the step functions in the time-series analysis. First, M earthquakes of magnitude greater than 4.0 were listed. Next, dates for which the horizontal or vertical displacement was greater than the threshold were sought for each site, namely, $$\sqrt{{dE^i_k}^2 + {dN^i_k}^2 } > 2$$ mm or $$|dU^i_k| > 5$$ mm. Here, $$dE^i_k$$, $$dN^i_k$$ and $$dU^i_k$$ are the east–west, north–south and up–down components of the displacement at the ith site due to the earthquake(s) that occurred on tk. These displacement components are calculated as:   \begin{eqnarray} dE^i_k = \sum _{m=1}^M \delta _{kD(m)} de^i_{m} \mbox{,} \end{eqnarray} (A1)  \begin{eqnarray} dN^i_k = \sum _{m=1}^M \delta _{kD(m)} dn^i_{m} \quad \mbox{and} \end{eqnarray} (A2)  \begin{eqnarray} dU^i_k = \sum _{m=1}^M \delta _{kD(m)} du^i_{m} \mbox{,} \end{eqnarray} (A3)where δkj and D(m) are the Kronecker Delta and the occurrence date of the mth earthquake, respectively, $$de^i_m$$, $$dn^i_m$$ and $$du^i_m$$ are the east–west, north–south and up–down components, respectively, of the displacement at the ith site due to the mth earthquake. The displacement at a GPS station due to an earthquake was calculated as follows. Estimating the area of the earthquake source fault based on its magnitude by applying the equation of the stress drop, Δσ, the earthquake moment, Mo and the fault area, S, derived from the dislocation theory under the assumption that the fault length is twice the width, i.e., $$\Delta \sigma = 2.5 M_o S^{-\frac{3}{2}}$$ (e.g. Kikuchi 2001). Δσ = 10 MPa was assumed, and Mo was calculated from its magnitude by regarding it to be the same as the moment magnitude. Calculating the fault length, L, from estimated fault area, S, as $$L=\sqrt{2S}$$. Estimating the fault slip amount from Mo and S by assuming the rigidity to be 40 GPa. Calculating the displacement at the ith site by means of Okada’s (1992) method. The following two fault mechanisms are assumed: Pure strike-slip faulting on a vertical fault along the direction rotated by 45° from the line between the epicentre of the mth earthquake and the ith GPS site. Pure thrust faulting on a dipping (30°) fault, the strike of which is perpendicular to the line between the epicentre of the mth earthquake and the ith GPS site. $$de^i_m$$, $$dn^i_m$$ and $$du^i_m$$ were obtained by taking the larger of the absolute values between the two cases for each component. The procedure may be improved by introducing centroid moment tensor (CMT) solutions from a catalogue, such as the Global CMT catalogue, but displacement steps at each GPS site should be estimated from the displacement time-series data. Therefore, the estimated displacement rates do not change if the dates for step functions are not changed. APPENDIX B: CALCULATION OF THE DISPLACEMENT RATE GRADIENT The spatial gradients of the displacement rate field are basically calculated by applying the weighted least-squares linear regression. Thus, the spatial gradients of horizontal and vertical displacement rates, gh and gv, are calculated as:   \begin{eqnarray} g_h &=& \frac{\displaystyle \left[ \sum _{n=1}^N \frac{1}{\left( {\delta a_h}_n\right)^2} \right] \left[ \sum _{n=1}^N \frac{x_n {a_h}_n}{ \left({\delta a_h}_n \right) ^2} \right] - \left[ \sum _{n=1}^N \frac{x_n}{\left( {\delta a_h}_n\right)^2} \right] \left[ \sum _{n=1}^N \frac{{a_h}_n}{ \left({\delta a_h}_n \right) ^2} \right] }{\displaystyle \left[ \sum _{n=1}^N \frac{1}{\left( {\delta a_h}_n\right)^2} \right] \left[ \sum _{n=1}^N \frac{x_n^2}{ \left({\delta a_h}_n \right) ^2} \right] - \left[ \sum _{n=1}^N \frac{x_n}{\left( {\delta a_h}_n\right)^2} \right]^2}, \quad \mbox{and} \end{eqnarray} (B1)  \begin{eqnarray} g_v &=& \frac{\displaystyle \left[ \sum _{n=1}^N \frac{1}{\left( {\delta a_v}_n\right)^2} \right] \left[ \sum _{n=1}^N \frac{x_n {a_v}_n}{ \left({\delta a_v}_n \right) ^2} \right] - \left[ \sum _{n=1}^N \frac{x_n}{\left( {\delta a_v}_n\right)^2} \right] \left[ \sum _{n=1}^N \frac{{a_v}_n}{ \left({\delta a_v}_n \right) ^2} \right] }{\displaystyle \left[ \sum _{n=1}^N \frac{1}{\left( {\delta a_v}_n\right)^2} \right] \left[ \sum _{n=1}^N \frac{x_n^2}{ \left({\delta a_v}_n \right) ^2} \right] - \left[ \sum _{n=1}^N \frac{x_n}{\left( {\delta a_v}_n\right)^2} \right]^2}, \end{eqnarray} (B2)where ahn, avn, δahn, δavn and xn are the horizontal and vertical displacement rates at the nth site in a swath, their estimation errors, and the distance of the nth site in the profile direction from the crossing point of the profile centreline and meridian E140°, respectively, and N is the number of sites in the swath. In practice, the displacement rate gradients and their errors, εh and εv, respectively, are estimated by different procedures depending on N, as follows: When N = 2, gh and gv are calculated using eqs (B1) and (B2). Moreover, εh and εv are calculated as $$\displaystyle \varepsilon _h = \sqrt{{\delta a_h}_1^2 + {\delta a_h}_1^2}$$, and $$\varepsilon _v = \sqrt{{\delta a_v}_1^2 + {\delta a_v}_1^2}$$. When 2 < N ≤ 10, gh and gv are calculated by taking the averages of ghn and gvn, where ghn and gvn (here, n = 1, …, N) are the spatial gradients calculated from the displacement rates at the sites in the swath, excluding the nth site. Moreover, εh and εv are defined as the standard deviations of ghn and gvn, respectively. When 10 < N ≤ 20, gh and gv are calculated by taking the averages of ghmn and gvmn, where ghmn and gvmn (here, n = 1, …, N − 1 and m = n + 1, …, N) are the spatial gradients calculated from the displacement rates at the sites in the swath, excluding the mth and nth sites. Here, εh and εv are the standard deviations of ghmn and gvmn, respectively. When 20 < N, gh and gv are calculated by taking the averages of ghlmn and gvlmn, respectively, where ghlmn and gvlmn (here, n = 1, …, N − 2, m = n + 1, …, N − 1 and l = m + 1, …, N) are the spatial gradients calculated from the displacement rates at the sites in the swath, excluding the lth, mth and nth sites. Here, εh and εv are defined as the standard deviations of ghlmn and gvlmn, respectively. Since displacement rates at GPS sites sometimes have systematic errors larger than the estimation errors based on the least-squares fitting of the displacement time-series, I adopted the above procedure to estimate the spatial gradients and their errors. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

Monitoring of the spatio-temporal change in the interplate coupling at northeastern Japan subduction zone based on the spatial gradients of surface velocity field

Loading next page...
 
/lp/ou_press/monitoring-of-the-spatio-temporal-change-in-the-interplate-coupling-at-2cz9J4wgv3
Publisher
The Royal Astronomical Society
Copyright
© The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx527
Publisher site
See Article on Publisher Site

Abstract

Summary A monitoring method to grasp the spatio-temporal change in the interplate coupling in a subduction zone based on the spatial gradients of surface displacement rate fields is proposed. I estimated the spatio-temporal change in the interplate coupling along the plate boundary in northeastern (NE) Japan by applying the proposed method to the surface displacement rates based on global positioning system observations. The gradient of the surface velocities is calculated in each swath configured along the direction normal to the Japan Trench for time windows such as 0.5, 1, 2, 3 and 5 yr being shifted by one week during the period of 1997–2016. The gradient of the horizontal velocities is negative and has a large magnitude when the interplate coupling at the shallow part (less than approximately 50 km in depth) beneath the profile is strong, and the sign of the gradient of the vertical velocity is sensitive to the existence of the coupling at the deep part (greater than approximately 50 km in depth). The trench-parallel variation of the spatial gradients of a displacement rate field clearly corresponds to the trench-parallel variation of the amplitude of the interplate coupling on the plate interface, as well as the rupture areas of previous interplate earthquakes. Temporal changes in the trench-parallel variation of the spatial gradient of the displacement rate correspond to the strengthening or weakening of the interplate coupling. We can monitor the temporal change in the interplate coupling state by calculating the spatial gradients of the surface displacement rate field to some extent without performing inversion analyses with applying certain constraint conditions that sometimes cause over- and/or underestimation at areas of limited spatial resolution far from the observation network. The results of the calculation confirm known interplate events in the NE Japan subduction zone, such as the post-seismic slip of the 2003 M8.0 Tokachi-oki and 2005 M7.2 Miyagi-oki earthquakes and the recovery of the interplate coupling around the rupture area of the 1994 M7.6 Sanriku-Haruka-oki earthquake. The results also indicate the semi-periodic occurrence of slow slip events and the expansion of the area of slow slip events before the 2011 Tohoku-oki earthquake (M9.0) approaching the hypocentre of the Tohoku-oki earthquake. Satellite geodesy, Seismic cycle, Transient deformation, Time-series analysis 1 INTRODUCTION Many large (M > 8) megathrust earthquakes occur in plate subduction zones where oceanic plates are subducting beneath continental plates. Such earthquakes rupture plate interface faults to release the strain energy due to the coupling between the continental and subducting plates during the interseismic period. Coseismic fault slip and interseismic fault coupling generate crustal deformation as well as quasi-static fault slip and transient processes such as visco-elastic relaxation and poro-elastic rebound during post-seismic periods. In order to consider the crustal deformation due to the interplate coupling during the interseismic period, Savage (1983) introduced the ‘backslip model’, in which the interplate coupling is expressed by two overlapping virtual components, namely, the stable sliding between the continental and subducting plates and the normal faulting at the areas where full or partial locking between the plates actually occurs. Plate subduction with stable sliding is regarded as producing no crustal deformation, with the exception of steady plate motion, while the normal faulting on the plate interface fault yields deformation in both plates. Following the proposal of the backslip model (Savage 1983), a number of studies have applied this idea to the estimation of the distribution of interplate coupling patches in subduction zones based on surface displacement fields (e.g. Hyndman & Wang 1995; Norabuena et al.1998; Lundgren et al.1999; Wallace et al.2004). A number of studies have examined the distribution of interplate coupling between the subducting Pacific and overriding continental plates in the northeast (NE) Japan subduction zone (Fig. 1) based on the surface displacement rate field by means of geodetic observations [mainly global positioning system (GPS) observations] (e.g. El-Fiky & Kato 1999; Mazzotti et al.2000; Ito et al.2000; Nishimura et al.2004; Suwa et al.2006; Hashimoto et al.2009; Loveless & Meade 2010). In such studies, the degree of coupling, which has generally been expressed in terms of the backslip (or slip deficit) rate at the plate interface, is estimated by applying linear inversion methods while imposing some constraint conditions, such as the smoothness of the spatial distribution of the backslip, the non-negativity of the backslip, Dirichlet boundary condition and the a priori distribution of the backslip (e.g. Yabuki & Matsu’ura 1992; Mosegaard & Tarantola 1995; Matsu’ura et al.2007; Meade & Loveless 2009). These constraint conditions prevent the solution from becoming unstable, because it is difficult to resolve the backslip distribution at the plate interface far off the coast based only on terrestrial observations. Figure 1. View largeDownload slide Map of the northeast Japan subduction zone. GPS stations are represented by circles, and the blue circles indicate stations used in this study. Note that this map shows the GPS stations for which the available displacement time-series are longer than one year during the entire analysis period. Examples of time-series analysis is presented for two stations, Oshika and Erimo, in Fig. 2. The brown dashed contours indicate the depth of the subducting plate interface (Nakajima & Hasegawa 2006; Nakajima et al.2009; Kita et al.2010). The grey line is the trace of Japan Trench. The intervals of the thin and thick contours are 10 and 50 km, respectively. The black thin lines indicate the triangular subfault used in the synthetic test. The colour contours indicate the rupture area of past large interplate earthquakes. Three events (M7.5, M7.3, and M7.4) in 1938 off Fukushima and Ibaraki Prefectures [Murotani (2003), green, purple and navy], four events off Miyagi Prefecture in 1936 [Yamanaka & Kikuchi (2004), M7.5, maroon], 1978 [Yamanaka & Kikuchi (2004), M7.4, magenta], 1981 [Yamanaka & Kikuchi (2004), M7.1, lime] and 2003 [Yamanaka (2003), M6.8, cyan], the Tokachi-oki earthquakes in 1968 [Yamanaka & Kikuchi (2004), M7.9, blue] and 2003 [Yamanaka & Kikuchi (2003), M8.0, teal], the Nemuro-oki earthquake in 1973 [Yamanaka (2005), M7.4, red] and the 2011 Off the Pacific Coast of Tohoku Earthquake [Iinuma et al. (2012), M9.0, black] are plotted. The names of prefectures used in the text are also shown. Figure 1. View largeDownload slide Map of the northeast Japan subduction zone. GPS stations are represented by circles, and the blue circles indicate stations used in this study. Note that this map shows the GPS stations for which the available displacement time-series are longer than one year during the entire analysis period. Examples of time-series analysis is presented for two stations, Oshika and Erimo, in Fig. 2. The brown dashed contours indicate the depth of the subducting plate interface (Nakajima & Hasegawa 2006; Nakajima et al.2009; Kita et al.2010). The grey line is the trace of Japan Trench. The intervals of the thin and thick contours are 10 and 50 km, respectively. The black thin lines indicate the triangular subfault used in the synthetic test. The colour contours indicate the rupture area of past large interplate earthquakes. Three events (M7.5, M7.3, and M7.4) in 1938 off Fukushima and Ibaraki Prefectures [Murotani (2003), green, purple and navy], four events off Miyagi Prefecture in 1936 [Yamanaka & Kikuchi (2004), M7.5, maroon], 1978 [Yamanaka & Kikuchi (2004), M7.4, magenta], 1981 [Yamanaka & Kikuchi (2004), M7.1, lime] and 2003 [Yamanaka (2003), M6.8, cyan], the Tokachi-oki earthquakes in 1968 [Yamanaka & Kikuchi (2004), M7.9, blue] and 2003 [Yamanaka & Kikuchi (2003), M8.0, teal], the Nemuro-oki earthquake in 1973 [Yamanaka (2005), M7.4, red] and the 2011 Off the Pacific Coast of Tohoku Earthquake [Iinuma et al. (2012), M9.0, black] are plotted. The names of prefectures used in the text are also shown. However, the interplate coupling distributions estimated by previous studies are not unique. The differences between the results may be caused by the differences of the analysis period, the geometrical model of the plate interface, the weight of the vertical displacement and the applied constraint conditions. In particular, the constraint conditions produce artefacts at the shallow plate interface, where the spatial resolution is poor in order to compensate for the errors of observation, the plate geometry and the Green’s function (because of inelastic, non-linear and/or anisotropic effects). Moreover, if we adopt the non-negativity constraint, it messes up the estimated backslip distribution when aseismic slip having a rate larger than the plate convergence rate occurs at some portion of the plate interface, because it requires negative backslip. Therefore, it is not easy to distinguish whether the interplate coupling state temporally changes or not by comparing the results of the inversion analyses. Methods to monitor the changes in the interplate coupling on the subduction plate interface based on the displacement rate field itself may shed light on the temporal variation of the coupling state. If the method can be automatically applied to the real-time geodetic data, it will play an important role in developing the operational monitoring system of the interplate coupling that is essential to predict megathrust earthquakes. Meanwhile, some studies have applied surface displacement rate field based on GPS observations to detect transient deformations in and around plate boundary zones (e.g. Ohtani et al.2010; Holt & Shcherbenko 2013; Crowell et al.2016). Ohtani et al. (2010) and Holt & Shcherbenko (2013) focused to the horizontal strain rate field to find transient deformation events and determine their locations, while Crowell et al. (2016) analysed 3-D displacement time-series at station by station to detect episodic tremor and slip events automatically. However, Ohtani et al. (2010) and Holt & Shcherbenko (2013) did not involve vertical displacement rates in their analysis even though Ohtani et al. (2010) mentioned the extensibility of their method to include the vertical component. In this study, I propose a method for monitoring the spatio-temporal variation of the interplate coupling based on calculating the spatial gradients of the horizontal and vertical surface displacement rate fields. Most of the anomalies detected in the displacement rate gradients correspond to previously reported slow slip events on the plate interface, and the remaining events are unidentified events. It is not possible to quantitatively estimate changes in interplate coupling by applying the monitoring method, although it is possible to grasp the spatio-temporal variation of the change in the interplate coupling in a qualitative manner and to constrain the trench-parallel widths of the perturbed areas and the durations of the events to some extent. 2 DATA AND ANALYSIS 2.1 Displacement time-series and displacement rates In this subsection, I explain the procedure of data processing for the displacement time-series in order to deduce the displacement rate fields. This curve fitting procedure does not discriminate the tectonic deformation signals from non-tectonic ones, but remove random noise due to the monument instabilities, meteorological effects and so on (e.g. Agnew 1992; Williams et al.2004) and steps due to earthquakes and station maintenance. Note that the effect of the spatially correlated noise due to the errors of the reference frame and the satellites’ orbits are reduced by taking spatial gradient of the displacement rate field as described in Section 2.2. Daily site coordinate time-series data obtained from the GPS observation network of the Geospatial Information Authority of Japan (F3 solution, Nakagawa et al.2009) were used to estimate the displacement rate at each GPS station (Fig. 1). The displacement rate is extracted from the site coordinate time-series by means of least-squares fitting for a mathematical function that consists of the long-term linear trend, annual and biannual trigonometric curves, and steps due to earthquakes and station maintenance, such as antenna and receiver replacement. The displacement time-series at the ith site is then expressed as:   \begin{eqnarray} \boldsymbol {u}^i ( t ) & = & \boldsymbol {u}_0^i + \sum _{j=1}^{J-1} \boldsymbol {a}_j^i \left( T_j -T_{j-1} \right) + \boldsymbol {a}_J^i \left( t - T_{J-1} \right)+ \boldsymbol {b}^i \cos \left( \frac{2 \pi t}{\tau } \right) +\boldsymbol {c}^i \sin \left( \frac{2 \pi t}{\tau } \right) + \boldsymbol {d}^i \cos \left( \frac{4 \pi t}{\tau } \right) + \boldsymbol {e}^i \sin \left( \frac{4 \pi t}{\tau } \right)\nonumber \\ &&+\, \sum _{k=1}^K \boldsymbol {f}_k^i H \left( t_k - t \right) \qquad \left( T_{J-1} < t \le T_J \right)\!, \end{eqnarray} (1)where $$\boldsymbol {u}(t)$$ is a displacement vector composed of east–west, north–south and up–down components at date t; Tj is the ending date of the jth sub-time-window in which the long-term linear trend is a constant, where T0 is the starting date of the time-window (refer to the following paragraphs for the time-window and sub-time-windows); tk is the date at which the earthquakes that were large enough to be considered occurred (see Appendix A) or station maintenance was performed; τ (=365.2422) is the number of days in one year; H(t) is the Heaviside step function; and $$\boldsymbol {a}, \boldsymbol {b}, \boldsymbol {c}, \boldsymbol {d}, \boldsymbol {e}, \boldsymbol {f}$$ and $$\boldsymbol {u}_0$$ are the coefficient vectors of the corresponding terms and a constant offset (initial value) vector. Because F3 solutions are usually updated on every Monday and are available since 1996 March 22, 5-yr time-windows ending every Monday from 1997 March 24 to 2016 December 26 were assigned. Before 2001 March 26, all daily coordinate data were applied although the time-window became shorter than 5 yr. For the long-term linear trends, two sets of sub-time-windows are configured to divide the 5-yr time-window into (i) five successive 1-yr sub-time-windows and (ii) two successive 4.5- and 0.5-yr sub-time-windows. The displacement values at the junction of the sub-time-windows are continuous. Then, the linear trend estimated for the latest 1- and 0.5-yr sub-time-windows were regarded as the average displacement rate during the last 1- and 0.5-yr sub-time-windows at the ending date of each 5-yr time window. Note that displacement time-series data during the first four 1- and 4.5-yr sub-time-windows are used for the estimation of annual and biannual trigonometric curves and steps but are not directly used for the average displacement rates at the last sub-time-windows. The time-window is able to be shortened from 5 yr to an arbitrary length as long as we can estimate the annual and biannual variations and steps robustly. It may be possible to optimize the lengths of the time-window and the sub-time-windows in an objective manner based on the observed time-series data, but it is beyond the scope of this paper. The average displacement rate during (iii) the entire 5-yr time-window was also estimated with no sub-time-window and while performing least-squares fitting in the same manner. Fig. 2(a) shows an example of the fitting. In this example, (A) when the five 1-yr sub-time-windows are configured, T0, T1, T2, T3, T4 and T5 in eq. (1) are assigned as 2006 March 07, 2007 March 07, 2008 March 07, 2009 March 07, 2010 March 07 and 2011 March 07, respectively, (B) when 4.5- and 0.5-yr sub-time-windows are configured, T0, T1 and T2 in eq. (1) are assigned as 2006 March 07, 2010 September 07 and 2011 March 07, (C) when no sub-time-window is configured, T0 and T1 in eq. (1) are assigned as 2006 March 07 and 2011 March 07, respectively. Then, t1, t2, t3, t4, t5, t6 and t7, which are the dates of moderate earthquakes and are assigned as 2006 November 15, 2007 January 13, 2008 May 08, June 14, July 19 and 25 and 2010 March 14, respectively, are independent of the sub-time-window setting. Finally, $$\boldsymbol {a}_5, \boldsymbol {a}_2$$ and $$\boldsymbol {a}_1$$ are regarded as the average displacement rates for the last 1- and 0.5-yr sub-time-windows and 5-yr time-window on 2011 March 07, just before the 2011 Tohoku-oki earthquake (M 9) on 2011 March 11. Figure 2. View largeDownload slide Examples of the results of the time-series analysis. (a) Case in which the 5-yr time-window from 2006 March 07 to 2011 March 07 is considered. The EW (top), NS (middle) and UD (bottom) components of the displacement time-series at site 0550 at Oshika are shown. The red dots indicate raw data. The green, blue and orange dots indicate the synthesized displacements obtained using a 5-yr time-window with the offsets (−5 mm for the EW and NS components and −10 mm for the UD component) five 1-yr sub-time-windows, and 4.5- and 0.5-yr sub-time-windows to estimate the linear trends with the offsets (5 mm for the EW and NS components and 10 mm for the UD component). Annual and semi-annual trigonometric curves and steps due to the earthquakes, for which the dates, magnitudes and locations, are exhibited, are also included in the synthesized displacements. Estimated linear trends for time-windows and sub-time-windows are indicted in corresponding colours. (b) Case in which the 6-yr time-window from 1997 September 22 to 2003 September 22 is considered. The EW (top), NS (middle) and UD (bottom) components of the displacement time-series at site 0019 at Erimo are shown. The red dots indicate raw data. The green and blue dots indicate the synthesized displacements using two 3-yr sub-time-windows with the offset (10 mm for all components) and three 2-yr sub-time-windows to estimate the linear trends. Annual and semi-annual trigonometric curves and steps due to earthquakes and site maintenance, for which the dates, magnitudes and locations are exhibited, are also included in the synthesized displacements. Estimated linear trends for sub-time-windows are indicted in corresponding colours. Figure 2. View largeDownload slide Examples of the results of the time-series analysis. (a) Case in which the 5-yr time-window from 2006 March 07 to 2011 March 07 is considered. The EW (top), NS (middle) and UD (bottom) components of the displacement time-series at site 0550 at Oshika are shown. The red dots indicate raw data. The green, blue and orange dots indicate the synthesized displacements obtained using a 5-yr time-window with the offsets (−5 mm for the EW and NS components and −10 mm for the UD component) five 1-yr sub-time-windows, and 4.5- and 0.5-yr sub-time-windows to estimate the linear trends with the offsets (5 mm for the EW and NS components and 10 mm for the UD component). Annual and semi-annual trigonometric curves and steps due to the earthquakes, for which the dates, magnitudes and locations, are exhibited, are also included in the synthesized displacements. Estimated linear trends for time-windows and sub-time-windows are indicted in corresponding colours. (b) Case in which the 6-yr time-window from 1997 September 22 to 2003 September 22 is considered. The EW (top), NS (middle) and UD (bottom) components of the displacement time-series at site 0019 at Erimo are shown. The red dots indicate raw data. The green and blue dots indicate the synthesized displacements using two 3-yr sub-time-windows with the offset (10 mm for all components) and three 2-yr sub-time-windows to estimate the linear trends. Annual and semi-annual trigonometric curves and steps due to earthquakes and site maintenance, for which the dates, magnitudes and locations are exhibited, are also included in the synthesized displacements. Estimated linear trends for sub-time-windows are indicted in corresponding colours. Similarly, I assigned 6-yr time-windows ending every Monday from 1997 March 24 to 2016 December 26 to calculate the displacement rates for the last 2- and 3-yr sub-time-windows at the ending dates. Each 6-yr time-window is divided into (D) two successive 3-yr sub-time-windows and (E) three successive 2-yr sub-time-windows, and the displacement values at the junction of the sub-time-windows are forced to be continuous. The linear trends are estimated by least-squares fitting, and the trends for the latest sub-time-windows are regarded as the average displacement rates during the last 3- and 2-yr sub-time-windows at the ending date of each 6-yr time-window. Fig. 2(b) shows an example of the fitting for 6-yr time-window. In this example, (D) when two 3-yr sub-time-windows are configured, T0, T1 and T2 in eq. (1) are assigned as 1997, September 22, 2000 September 22 and 2003 September 22, respectively, (E) when three 2-yr sub-time-windows are configured, T0, T1, T2 and T3 in eq. (1) are assigned as 1997 September 22, 1999 September 22, 2001 September 22 and 2011 September 22, respectively. Moreover, t1, t2 and t2, which are the date of an M 7.0 earthquake on 2000 January 28 and the dates of antenna replacements on 2001 December 07 and 2002 November 20, are independent of the sub-time-window setting. Finally, $$\boldsymbol {a}_2$$ and $$\boldsymbol {a}_3$$ are regarded as the average displacement rates for the last 3- and 2-yr sub-time-windows on 2003 September 22, just before the 2003 Tokachi-oki (M 8) earthquake on 2003 September 23. Applying these time-series fitting procedures to coordinate time-series of all GPS stations, we can obtain weekly surface displacement rate fields for the last 0.5-, 1-, 2- and 3-yr sub-time-windows, and 5-yr time-window. Naturally, the displacement rate fields for the longer windows are more robust than those for the shorter windows. Figs 3 and 4, and Fig. S1 in the Supporting Information show the examples of estimated displacement rate fields. The horizontal displacement rate field for the last 0.5-yr sub-time-window (Fig. S1, Supporting Information) shows many outlier vectors in directions different from the general trend. Meanwhile, the vertical displacement rate fields for the last 0.5- and 1-yr sub-time-windows include short-wavelength spatial variations with large amplitudes. For instance, the displacement rate at a site at about 40 km in the distance from E140° in the profiles shown in Fig. S1 in the Supporting Information has difference from that at the neighbourhood sites as large as 1σ error (∼10 mm yr−1), while there are no such outlier in the vertical displacement rate profiles for the last 2- and 3-yr sub-time-windows. Therefore, the displacement rate fields are not robustly estimated for the last 0.5-yr sub-time-window with respect to the horizontal component and for the last 0.5- and 1-yr sub-time-windows with respect to the vertical component. Figure 3. View largeDownload slide Examples of the displacement rate profiles and spatial gradients for the Tohoku region. (a) Case in which the centreline of the swath runs through at N38.7°, E140.0° is considered. The displacement rates at sites for which the distance to the swath centre does not exceed 30 km (red vectors) are chosen in order to calculate the spatial gradient. The horizontal component of the displacement rate field for a 5-yr time-window from 2006 March 07 to 2011 March 07 is shown. The bottom panel shows the profile of the horizontal displacement rates along the direction of the swath (N105°E) with respect to the distance from the trench. (b) The vertical components for the case described in (a) are plotted. Figure 3. View largeDownload slide Examples of the displacement rate profiles and spatial gradients for the Tohoku region. (a) Case in which the centreline of the swath runs through at N38.7°, E140.0° is considered. The displacement rates at sites for which the distance to the swath centre does not exceed 30 km (red vectors) are chosen in order to calculate the spatial gradient. The horizontal component of the displacement rate field for a 5-yr time-window from 2006 March 07 to 2011 March 07 is shown. The bottom panel shows the profile of the horizontal displacement rates along the direction of the swath (N105°E) with respect to the distance from the trench. (b) The vertical components for the case described in (a) are plotted. Figure 4. View largeDownload slide Examples of the displacement rate profiles and spatial gradients for the Hokkaido region. (a) Case in which the centreline of the swath runs through N43.5°, E140.0° is considered. The displacement rates at sites for which the distance to the swath centre does not exceed 20 km (red vectors) are chosen in order to calculate the spatial gradient. The horizontal component of the displacement rate field for the 2-yr sub-time-window from 2001 September 22 to 2003 September 22 is shown. The bottom panel shows the profile of the horizontal displacement rates along the direction of the swath (N120°E) with respect to the distance from the trench. (b) The vertical components for the case described in (a) are plotted. Figure 4. View largeDownload slide Examples of the displacement rate profiles and spatial gradients for the Hokkaido region. (a) Case in which the centreline of the swath runs through N43.5°, E140.0° is considered. The displacement rates at sites for which the distance to the swath centre does not exceed 20 km (red vectors) are chosen in order to calculate the spatial gradient. The horizontal component of the displacement rate field for the 2-yr sub-time-window from 2001 September 22 to 2003 September 22 is shown. The bottom panel shows the profile of the horizontal displacement rates along the direction of the swath (N120°E) with respect to the distance from the trench. (b) The vertical components for the case described in (a) are plotted. 2.2 Spatial gradient of the displacement rate field Next, the spatial gradients of the displacement rate fields were calculated for each swath. Swaths were configured along the direction near the plate convergence, namely, N105°E for the Tohoku region and N120°E for the Hokkaido region. The centrelines of the swaths were defined as passing the reference points configured at every 0$${^{\circ}_{.}}$$1 of latitude from N34.5 to N42.0 for the Tohoku region and from N43.0 to N46.5 for the Hokkaido region along meridian E140°. The width of the swath was configured as 40, 60, or 80 km. Spatial gradients of the surface horizontal and vertical displacement rates are estimated with linear regression for the displacement rates of GPS stations inside the swaths (See Appendix B for the details of the calculation and error estimation of the displacement rate gradient). The horizontal displacement rates are projected into the direction of the swaths using the east–west and north–south components to calculate the spatial gradients. The displacement rates estimated from time-series shorter than one year are excluded. Stations at longitudes of less than E139.5° are also excluded. An example of displacement rate field and spatial gradients for the last 5-yr time-window on 2011 March 07 is exhibited in Fig. 3. Fig. 3 also shows the displacement rate profile in a swath of which centreline passes through N38.7°, E140.0°. GPS stations for which the distance to the swath centreline is less than 30 km were included in the calculation of the spatial gradient (the width of the swath is as 60 km). Since the period is just before the 2011 Tohoku-oki earthquake, interplate coupling off Miyagi prefecture was strong, as reflected by large negative value of the horizontal displacement rate gradient. The Pacific coast around Miyagi Prefecture was subsiding during the interseismic period (e.g. Kato & Tsumura 1979; El-Fiky & Kato 1999; Aoki & Scholz 2003; Nishimura et al.2004), and the vertical displacement rate gradient in this profile has a large negative value. Fig. 4 shows the example for the last 2-yr sub-time-window on 2003 September 22. Almost the same features are observed in the profile for which the centreline passes though N43.5°, E140.0° just before the 2003 Tokachi-oki earthquake (Figs 4a and b). However, since, in this case, the width of the swath was 40 km, whereas it was 60 km in the previous case, the number of sites in the swath were fewer. Moreover, as the displacement rates were estimated for the last 2-yr sub-time-window, the error for each displacement rate was larger than in the previous case. Therefore, the spatial gradient of the displacement rates was generally less robust than in the case where the swath was wide and/or the displacement rates were estimated for a longer sub-time-window. This means that there is some trade-off between the robustness of the spatial gradient of the displacement rate field and the resolution of the spatial and temporal change in the spatial gradients. By plotting spatial gradients in all swaths, we can easily grasp the spatial change in the displacement rate gradients along the trench-parallel direction. Fig. 5 shows examples of the distributions of the spatial gradient of the displacement rate fields. The width of the swaths and the length of sub-time-window in which the displacement rates were estimated correspond to the cases shown in Figs 3 and 4. Trench-parallel variations of the displacement rate gradients are clearly observed in each panel, which may reflect the spatial variation in the degree of coupling on the subducting plate interface. For example, Fig. 5(a) shows opposite sign in vertical components between Tohoku and Hokkaido, while the horizontal components show negative spatial gradient in both regions. This difference may be caused by the post-seismic deformation associated with the 2003 Tokachi-oki earthquake, and indicates that the afterslip on the deep plate interface continues longer than that at the shallow portion. Meanwhile, the spatial variation of vertical displacement rate gradient in Tohoku district shown in Fig. 5(b) is more heterogeneous than that shown in Fig. 5(a). Some of the vertical displacement rate gradients shown in Fig. 5(b) are not significantly larger than the estimation errors, but the positive gradients around N39° must be caused by the post-seismic deformation associated with an M7.1 earthquake off Miyagi Prefecture on 2003 May 26 (Okada & Hasegawa 2003). Figure 5. View largeDownload slide Examples of the spatial (trench-parallel) variations of the displacement rate gradients. The horizontal and vertical components are represented by red and blue bars, the origins of which correspond to the intersection points of the swath centrelines and the trench axis. The lengths of the bars are proportional to the absolute values of the displacement rate gradients. The scales are shown in the figure. The bars toward the southeast (seaward) indicate that the signs of the gradients are positive, whereas the bars toward the northwest (landward) indicate that the signs are negative. The thicknesses of the colours of bars correspond to the estimation errors. (a) The spatial gradients of the displacement rate field for the 5-yr time-window from 2006 March 07 to 2011 March 07 is shown together with the horizontal component of the displacement rates. [These vectors are identical to those in Fig. 3(a).] As in Figs 3(a) and (b), the width of each swath is 60 km. The spatial gradients of the displacement rate field for the 2-yr sub-time-window from 2001 September 22 to 2003 September 22 is presented along with the vertical component of the displacement rates. [These vectors are identical to those shown in Fig. 4(b).] As in Figs 4(a) and (b), the width of each swath is 40 km. Figure 5. View largeDownload slide Examples of the spatial (trench-parallel) variations of the displacement rate gradients. The horizontal and vertical components are represented by red and blue bars, the origins of which correspond to the intersection points of the swath centrelines and the trench axis. The lengths of the bars are proportional to the absolute values of the displacement rate gradients. The scales are shown in the figure. The bars toward the southeast (seaward) indicate that the signs of the gradients are positive, whereas the bars toward the northwest (landward) indicate that the signs are negative. The thicknesses of the colours of bars correspond to the estimation errors. (a) The spatial gradients of the displacement rate field for the 5-yr time-window from 2006 March 07 to 2011 March 07 is shown together with the horizontal component of the displacement rates. [These vectors are identical to those in Fig. 3(a).] As in Figs 3(a) and (b), the width of each swath is 60 km. The spatial gradients of the displacement rate field for the 2-yr sub-time-window from 2001 September 22 to 2003 September 22 is presented along with the vertical component of the displacement rates. [These vectors are identical to those shown in Fig. 4(b).] As in Figs 4(a) and (b), the width of each swath is 40 km. Fig. 6 shows the distribution of the average estimated error. The figure clearly indicates that the horizontal displacement rate gradients for sub-time-windows not shorter than one year have estimation errors of less than 1 and 2 mm yr−1 100 km−1 in the Hokkaido region and the Tohoku region north to N38.5° (see Fig. S3, Supporting Information for the 1- and 3-yr sub-time-windows). Since the distances between the coastline and longitude line of E139.5° are small in the Tohoku region south to N38.5°, the estimation errors of the horizontal displacement rates are sometimes larger than 2 mm yr−1 100 km−1 for the swaths for which the width is less than 80 km. With respect to the vertical displacement component, the estimation errors are less than 2 and 5 mm yr−1 100 km−1 in the Hokkaido and Tohoku regions, respectively, for sub-time-windows longer than one year. For most swaths having a width larger than 40 km, the errors are less than 2 mm yr−1 100 km−1 in the area north to N38.5° in the Tohoku region. Figure 6. View largeDownload slide Trench-parallel distribution of the average estimation error of the displacement rate gradients. The upper panel shows the vertical component, whereas the lower panel shows the horizontal component. All cases for the width of the swaths are plotted using various styles of lines (see the legends in the panels). The blue and red lines correspond with 0.5- and 2-yr sub-time-windows. The black lines show the cases for 5-yr time-window. Figure 6. View largeDownload slide Trench-parallel distribution of the average estimation error of the displacement rate gradients. The upper panel shows the vertical component, whereas the lower panel shows the horizontal component. All cases for the width of the swaths are plotted using various styles of lines (see the legends in the panels). The blue and red lines correspond with 0.5- and 2-yr sub-time-windows. The black lines show the cases for 5-yr time-window. Next, I performed synthetic tests in order to understand the relationship between the displacement gradients and the interplate coupling distribution and to examine the detectability while taking the estimation errors into account. 3 INTERPLATE COUPLING AND DISPLACEMENT RATE GRADIENTS 3.1 Synthetic tests In order to confirm the relationship between the spatial gradients of the surface displacement rate field and the backslip distribution on the plate interface, I carried out synthetic tests by applying forward calculations based on dislocation models while assuming various backslip distributions. I used the plate interface model devised by Nakajima & Hasegawa (2006), Nakajima et al. (2009) and Kita et al. (2010). The plate interface fault shallower than 100 km in depth is approximated using triangular subfaults, the edges of which are approximately 20 km in length, to calculate the surface displacements due to the dislocation at each subfault using the computer code provided by Meade (2007). These subfaults are divided according to depth into ten 10-km ranges. For instance, the first range is composed of triangular subfaults shallower that 10 km and the jth range include depths of (j − 1) × 10 to j × 10 km (e.g. the 10th range is composed of subfaults at 90–100 km in depth). The surface displacement rate field due to normal-fault type slip on the subfaults that belong to each range with a slip rate of 4 cm yr−1 (corresponding to approximately half the plate convergence rate) or 8 cm yr−1 (approximately equal to the plate convergence rate) was calculated. As a result, approximately 50 per cent, or full locking, occurs on the plate interface. Fig. 7 shows the calculated displacement rate fields for four of the (310 − 1) cases. The total number of cases is obtained from zero, half and full locking for 10 ranges, excluding the case in which zero locking occurs in every range. I chose the cases in which full locking at 30–50 km in depth (Fig. 7a), half locking at shallower than 30 km in depth, and full locking at 30–50 km in depth (Fig. 7b), full locking at shallower than 50 km in depth (Fig. 7c), and half locking at shallower than 30 or 50–90 km in depth, and full locking at 30–50 km in depth (Fig. 7d) are assumed as typical cases. Figure 7. View largeDownload slide Spatial distributions of the displacement rate gradients using the synthetic displacement rate fields calculated from the assumed interplate coupling distributions. (a) Full coupling (8 cm yr−1) on the plate interface at 30–50 km in depth is assumed (green). (b) Half coupling (4 cm yr−1) on the plate interface shallower than 30 km (turquoise) is added to the case shown in (a). (c) Full coupling (8 cm yr−1) on the plate interface shallower than 50 km is assumed. (d) Half coupling (4 cm yr−1) on the plate interface at 50 to 90 km in depth (turquoise) is added to the cases shown in (b). The horizontal displacement rate field is shown in (a) and (b) along with the displacement rate gradient distribution, while vertical displacement rates are plotted in (c) and (d). Figure 7. View largeDownload slide Spatial distributions of the displacement rate gradients using the synthetic displacement rate fields calculated from the assumed interplate coupling distributions. (a) Full coupling (8 cm yr−1) on the plate interface at 30–50 km in depth is assumed (green). (b) Half coupling (4 cm yr−1) on the plate interface shallower than 30 km (turquoise) is added to the case shown in (a). (c) Full coupling (8 cm yr−1) on the plate interface shallower than 50 km is assumed. (d) Half coupling (4 cm yr−1) on the plate interface at 50 to 90 km in depth (turquoise) is added to the cases shown in (b). The horizontal displacement rate field is shown in (a) and (b) along with the displacement rate gradient distribution, while vertical displacement rates are plotted in (c) and (d). The results indicate that full locking on the plate interface at intermediate depth (∼30 to ∼50 km, Fig. 7a) produces significant negative spatial gradients of the horizontal displacement rates, namely, ∼ −6 and ∼ −4 mm yr−1 100 km−1 in Tohoku and Hokkaido regions, respectively, while large negative spatial gradients of the horizontal displacement rates, namely, ∼ −11 and ∼ −7 mm yr−1 100 km−1 in the Tohoku and Hokkaido regions, respectively, are produced by full locking on the plate interface at shallow and intermediate depths (less than ∼50 km, Fig. 7c). Taking into account the fact that the estimation errors for the spatial gradient of the horizontal displacement rate are less than 2 mm yr−1 100 km−1 in most of the analysis regions (Fig. 6), such a drastic change in the interplate coupling state can be monitored by means of the spatial gradients of the horizontal displacement rates. On the other hand, small (0 to half or half to full) changes in the interplate coupling at the shallow (less than 30 km) plate interface may not be well detected by monitoring the spatial gradients of the horizontal displacement rates, because interplate coupling produces a change in the horizontal displacement rates that may be as large as the estimation errors (Figs 7b and 6). The difference between the horizontal displacement rate fields shown in Figs 7(b) and (d) indicates that the change in the coupling state on the deep (greater than 50 km) plate also causes an insignificant change in the spatial gradients of the horizontal displacement rates. The spatial gradients of the vertical displacement rates have almost no dependency on the change in the interplate coupling state on the intermediate and shallow portions (Figs 7a–c). The difference between the vertical displacement rate fields shown in Figs 7(a) and (c) produces the difference between the spatial gradients of the vertical displacement rates smaller than 2 mm yr−1 100 km−1, which is not significantly larger than the estimation error (Fig. 6). Therefore, we can deduce no information about the change in the degree of coupling on the plate interface at intermediate and shallow depths by monitoring the spatial gradients of the vertical displacement rates. However, the spatial gradients of the vertical displacement rates are prominently affected by the interplate coupling at the deep portion. The sign of the spatial gradients of the vertical displacement rates become negative when the interplate coupling occurs at the deep plate interface, namely, when the plate interface is located beneath land, that is, beneath the GPS observation network (e.g. Fig. 7d). This inversion of the spatial gradient of the vertical displacement rate is natural considering that the hinge line is located above the downdip limit of the dip-slip fault (e.g. Segall 2010). Next, I performed another set of the synthetic tests to examine the constrainablitiy of trench-parallel expansion of the area in which the interplate coupling changes based on the change in the trench-parallel distribution of the displacement rate gradients. The surface displacement rate fields were calculated by modifying the interplate coupling distribution shown in Fig. 7(c), namely, full locking on the plate interface shallower than 50 km in depth. Figs 8(a) and (c) show the examples of the cases in which the coupling in certain areas on the plate interface shallower than 30 km in depth is weakened to 4 cm yr−1. The trench-parallel expansion of a weak coupling zone was configured according to the distance from the centreline of a swath as shown in Fig. 8. The interplate coupling at an area within 40 km from the swath centreline that runs through at N39.0°, E140.0° was weakened in the case shown in Fig. 8(a), while the case in which the swath centreline runs through at N44.0°, E140.0° is exhibited in Fig. 8(c). The displacement rate gradient distribution was calculated and the distribution of the difference of the displacement rate gradient was taken by subtracting that due to the assumed interplate coupling distributions with the weak coupling zone (e.g. Figs 8a and c) from that due to the assumed interplate coupling distributions without the weak coupling zone (Fig. 7c). Changing the centreline of the weak coupling zone as passing the reference points configured at every 1° of latitude from N35 to N41.0 and from N44.0 to N46.0 along meridian E140°, the distributions of the difference were obtained as shown in Figs 9(a) and (c). Figure 8. View largeDownload slide Spatial distributions of the displacement rate gradients using the synthetic displacement rate fields calculated from the assumed interplate coupling distributions with perturbations given by half coupling zones to the coupling distribution shown in Fig. 7(c). (a) The coupling on the plate interface shallower than 30 km in depth and within 40 km from the swath centreline that runs through at N39.0°, E140.0° is weakened to 4 cm yr−1. (b) The half coupling at 50–90 km in depth within 40 km from the swath centreline that runs through at N39.0°, E140.0° is added. (c) The coupling on the plate interface shallower than 30 km in depth and within 40 km from the swath centreline that runs through at N44.0°, E140.0° is weakened to 4 cm yr−1. (d) The half coupling at 50–90 km in depth within 40 km from the swath centreline that runs through at N44.0°, E140.0° is added. The horizontal displacement rate field is shown in (a) and (c) along with the displacement rate gradient distribution, while vertical displacement rates are plotted in (b) and (d). Figure 8. View largeDownload slide Spatial distributions of the displacement rate gradients using the synthetic displacement rate fields calculated from the assumed interplate coupling distributions with perturbations given by half coupling zones to the coupling distribution shown in Fig. 7(c). (a) The coupling on the plate interface shallower than 30 km in depth and within 40 km from the swath centreline that runs through at N39.0°, E140.0° is weakened to 4 cm yr−1. (b) The half coupling at 50–90 km in depth within 40 km from the swath centreline that runs through at N39.0°, E140.0° is added. (c) The coupling on the plate interface shallower than 30 km in depth and within 40 km from the swath centreline that runs through at N44.0°, E140.0° is weakened to 4 cm yr−1. (d) The half coupling at 50–90 km in depth within 40 km from the swath centreline that runs through at N44.0°, E140.0° is added. The horizontal displacement rate field is shown in (a) and (c) along with the displacement rate gradient distribution, while vertical displacement rates are plotted in (b) and (d). Figure 9. View largeDownload slide Spatial distributions of the difference between the displacement rate gradients using the synthetic displacement rate fields calculated from assumed interplate coupling distribution shown in Fig. 7(c) and from perturbed coupling distribution (examples are shown in Fig. 8. (a) and (c) The backslip on the plate interface shallower than 30 km and within 40 km from each swath centreline is weakened to 4 cm yr−1. The horizontal axis indicates the latitude of the crossing point of the centreline of the swath and meridian E140°. The legend of each coloured line indicates the latitude of the crossing point of the centreline of the area where the interplate coupling is perturbed and meridian E140°. (b) and (d) Half coupling (4 cm yr−1) on the plate interface at 50–90 km in depth is added. The differences of the horizontal displacement rate gradient are plotted in (a) and (b), while the vertical components are shown in (c) and (d). Figure 9. View largeDownload slide Spatial distributions of the difference between the displacement rate gradients using the synthetic displacement rate fields calculated from assumed interplate coupling distribution shown in Fig. 7(c) and from perturbed coupling distribution (examples are shown in Fig. 8. (a) and (c) The backslip on the plate interface shallower than 30 km and within 40 km from each swath centreline is weakened to 4 cm yr−1. The horizontal axis indicates the latitude of the crossing point of the centreline of the swath and meridian E140°. The legend of each coloured line indicates the latitude of the crossing point of the centreline of the area where the interplate coupling is perturbed and meridian E140°. (b) and (d) Half coupling (4 cm yr−1) on the plate interface at 50–90 km in depth is added. The differences of the horizontal displacement rate gradient are plotted in (a) and (b), while the vertical components are shown in (c) and (d). Fig. 9(a) shows that the difference of the horizontal displacement gradient mostly distributes around the latitudes at which the centreline of the weak coupling zone is configured, even though the amplitude of the difference depends on the location of the weak coupling zone. The amplitude of the difference may be not large enough to be detected taking the estimation error of the horizontal displacement rate gradient that is discussed in the previous section into account. However, because the change in horizontal displacement rate gradient is almost proportional to the change in the coupling on the shallow plate interface, we can detect the change in the interplate coupling at the shallow portion between full coupling and no coupling based on the monitoring of the displacement rate gradient. Because the number of GPS stations are small, the weak coupling zones near the tips of the land do not change the displacement gradient distribution large enough to be detected. However, it is possible to constrain the trench-parallel expansion of an area where the interplate coupling on the shallow plate interface changes to some extent based on the monitoring of the horizontal displacement rate gradient. Note that the change in the coupling state at the shallow plate interface causes no significant difference of the vertical displacement gradient as Fig. 9(c) shows. Fig. S4 in the Supporting Information shows other cases in which the widths of the swath (40, 60 and 80 km) and the width of the week coupling zone (40, 80, 120, 160 and 200 km) are examined. The distribution of the difference sometimes becomes bimodal when the horizontal displacement rate gradient in narrow (40 km in width) swath is calculated for the displacement field synthesized by assuming wide (≥80 km) weak coupling zone. In such case, we should utilize wider swath to calculate the displacement rate gradient, and it is possible to estimate the true centreline of the weak coupling zone. With respect to the change in the coupling state on the plate interface at deep portion, the cases in which the weak coupling zone at 50–90 km in depth is added were examined. Figs 8(b) and (d) show the examples of the assumed interplate coupling distributions. The distributions of the difference of the displacement rate gradients were calculated in the same manner as the tests with respect to the weak coupling zone in the shallow portion (Figs 9b and d). It shows that the changes in the vertical displacement rate gradient are concentrated around added deep coupling zones, while the distribution of the difference with respect to the horizontal component do not well correspond with the change in the interplate coupling. Fig. S5 in the Supporting Information shows other cases as well as Fig. S4 in the Supporting Information. The relationship between the widths of the swath and the weak coupling zone in the deep portion has the same tendency as the relationship between the widths of the swath and the weak coupling zone in the shallow portion as described above. In summary, we can detect the temporal change in the coupling state on the plate interface offshore and onshore in the NE Japan subduction zone to a certain degree by monitoring the temporal changes in the spatial gradients of the horizontal and vertical displacement rates, respectively. And it is possible to constrain the location of the area where the interplate coupling change occurs to some extent. Of course, it is difficult to distinguish the ranges at which the change in the degree of coupling actually occurs among depth ranges of 10–50 km. However, low spatial resolution is inevitable in the trench-normal direction, even if we apply geodetic inversion analysis, as long as only onshore GPS data are used. 4 SPATIO-TEMPORAL CHANGE IN INTERPLATE COUPLING 4.1 Spatio-temporal variation in the velocity gradient Fig. 10 shows the spatio-temporal variations of the horizontal and vertical displacement rate gradients for the last 1-yr sub-time-windows and 5-yr time-windows in the 60-km-wide swaths. (See Fig. S6 in the Supporting Information for the cases of other sub-time-window lengths and swath widths.) After the 2011 Tohoku-oki earthquake, the displacement rate field in the Tohoku region is strongly influenced by the visco-elastic response (e.g. Watanabe et al.2014; Sun et al.2014; Freed et al.2017), and it is difficult to deduce the information with respect to the interplate coupling based on the spatio-temporal changes in the displacement rate gradients. Figure 10. View largeDownload slide Spatio-temporal variations of the displacement rate gradients. Results for (a) the 1-yr sub-time-window and (b) the 5-yr time-window. The horizontal axis indicates the latitudes of the crossing point of the centreline of the swath and meridian E140°. The panels on the left-hand side correspond to the Kanto to Tohoku regions, and the panels on the right-hand side correspond to the Hokkaido region. The vertical axis indicates the final dates of the sub-time-windows and time-windows. The results for the 1-yr sub-time-window and the 5-yr time-window include 1032 Mondays when the sub-time-windows end from 1997 March 24 to 2016 December 26 and 823 Mondays when the time-windows end from 2001 March 26 to 2016 December 26, respectively, but the half of them (every other Monday) are plotted in each panel. The top and bottom panels show the vertical and horizontal displacement rate gradients, respectively. Note that different colour scales are used for the horizontal and vertical components. Figure 10. View largeDownload slide Spatio-temporal variations of the displacement rate gradients. Results for (a) the 1-yr sub-time-window and (b) the 5-yr time-window. The horizontal axis indicates the latitudes of the crossing point of the centreline of the swath and meridian E140°. The panels on the left-hand side correspond to the Kanto to Tohoku regions, and the panels on the right-hand side correspond to the Hokkaido region. The vertical axis indicates the final dates of the sub-time-windows and time-windows. The results for the 1-yr sub-time-window and the 5-yr time-window include 1032 Mondays when the sub-time-windows end from 1997 March 24 to 2016 December 26 and 823 Mondays when the time-windows end from 2001 March 26 to 2016 December 26, respectively, but the half of them (every other Monday) are plotted in each panel. The top and bottom panels show the vertical and horizontal displacement rate gradients, respectively. Note that different colour scales are used for the horizontal and vertical components. The spatio-temporal change in the horizontal displacement rate gradients for 5-yr time windows (Fig. 10) clearly indicates the transient change in the interplate coupling that were reported in previous studies. For instance, The decrease in the negative gradients at N36.5° to N38.5° in 2008–2011 (before the Tohoku-oki earthquake) corresponding to the weak interplate coupling at the shallow and intermediate plate interfaces, as described in the previous section (Ozawa et al.2012; Mavrommatis et al.2014; Yokota & Koketsu 2015). The increase in the negative gradients around N41.0° before the 2011 M 9.0 earthquake representing the recovery process of the interplate coupling at and around the rupture area of the 1994 Sanriku-Haruka-oki earthquake (Mw7.7, Ozawa et al.2007). The dissipation of large negative gradients at N43.5° to N45.0° around 2006 corresponding to the afterslip of the 2003 Tokachi-oki earthquake (M 8.0, Miyazaki et al.2004; Murakami et al.2006). The effect of the afterslip due to the 2005 Miyagi-oki earthquake (M 7.2 on 16 August 2005) can be observed, but the magnitude of the change in the displacement gradient is comparatively smaller than in the cases of the above large and long-term events. Considering the vertical component, negative displacement gradients of up to −5 ± 2 mm yr−1 100 km−1 are consistently observed in the Tohoku region, N37.5° to N40.5°, before 2011, and of up to −3 ± 2 mm yr−1 100 km−1 are consistently observed in the Hokkaido region, N44.2° to N46.2°, before 2005. This implies that the interplate coupling at deep plate interfaces was strong before megathrust earthquakes, such as the 2011 Tohoku-oki and 2003 Tokachi-oki earthquakes. Since the displacement rate fields that are estimated for 5-yr time-windows do not contain localized deformation as shown in Fig. 3, the spatial gradients are robustly estimated, but temporal changes in shorter periods cannot be monitored. The results shown in Figs 10(a) and S6 indicate that there are significant changes in the horizontal displacement rate gradients for periods shorter than five years. This indicates that the interplate coupling on the plate interface at shallow and/or intermediate depths may change depending over the course of a couple of years. 4.2 Periodic change in interplate coupling The spatio-temporal variation in the horizontal displacement rate gradients for a 1-yr sub-time-window (Fig. 10a) reveals that the decrease in the negative gradients at N36.5° to N38.5° did not occur monotonically. The decreases occurred around 1997, 2001, 2005 and 2009, and the northern limit of the fluctuating profiles advanced northward, namely, the decrease in the displacement rate gradient occurred around N37.0° in 1997, around N37.4° in 2001, around N37.6° in 2005 and around N38.5°, in 2009. Note that the very small (≤ −2.5 mm yr−1 100 km−1, coloured blue) displacement rate gradients estimated during 2006 between N38.0° and N39.0° correspond to the afterslip of the 2005 Miyagi-oki earthquake (M7.2) (Miura et al.2006). Fig. 11 shows a time-series of the horizontal displacement rate gradients for 3- and 5-yr sub-time-windows. The quasi-periodic fluctuation events described above are also observed in these time-series in Fig. 11(a). For instance, the increase (decrease of the negative value) of the displacement rate gradient occurred for the profiles between N36.2° and N37.4° during the period of 2001–2002, between N36.8° and N37.8° during the period of 2005–2006 and between N36.8° and N38.4° during the period of 2009–2010. Note that the periods of fluctuation are shifted approximately one year later here, as compared to the discussion of the previous paragraph, because the length of sub-time-window (three years) is longer than one year. On the other hand, we cannot clearly observe the fluctuations during the periods of 2001–2002 and 2005–2006 based on the time-series of the displacement rate gradients of 5-yr time-windows shown in Fig. 11(b). This indicates that the duration of the change in the interplate coupling at shallow and intermediate depths off Fukushima and Ibaraki Prefectures is shorter than five years. Figure 11. View largeDownload slide Time-series of the horizontal displacement rate gradients for the 40-km-wide swaths at N36.0° to N39.0° before the 2011 Tohoku-oki earthquake. The results for the (a) 3-yr and (b) 5-yr time-windows are shown. The data for the swath are plotted with colours corresponding to the profile latitude at E140.0°, as shown at the right-hand side of each panel. The latitude numbers used to label the horizontal lines correspond to the zero displacement rate gradients for each profile. The grey squares indicate fluctuations mentioned in the text. Figure 11. View largeDownload slide Time-series of the horizontal displacement rate gradients for the 40-km-wide swaths at N36.0° to N39.0° before the 2011 Tohoku-oki earthquake. The results for the (a) 3-yr and (b) 5-yr time-windows are shown. The data for the swath are plotted with colours corresponding to the profile latitude at E140.0°, as shown at the right-hand side of each panel. The latitude numbers used to label the horizontal lines correspond to the zero displacement rate gradients for each profile. The grey squares indicate fluctuations mentioned in the text. 5 DISCUSSION As described in Section 2.1, the displacement rate field deduced from the time-series analysis performed in this study includes both tectonic and non-tectonic deformations. It is quite difficult to distinguish and eliminate non-tectonic signals due to causes such as ocean tide loading, atmospheric mass loading and continental water loading from the estimated displacement rate field in a strict manner. Nevertheless, the deformation due to the seasonal water variations are able to be removed by including annual and biannual trigonometric curves in the time-series fitting to some extent, and non-seasonal non-tectonic deformation with a short time-scale should be disappeared by taking the long-term trends. Comparing between the results of different sub-time-windows, swath widths and components may help us to distinguish long-term non-tectonic deformation from the tectonic deformation. Recently, Uchida et al. (2016) revealed the periodic changes in the interplate coupling in the NE Japan subduction zone based primarily on the activity of small repeating earthquakes. They found that slow slip on the plate interface had occurred repeatedly at intervals of 1–6 yr, depending on the location. They also used the horizontal displacement gradient time-series of the 1-yr sub-time-window to estimate the recurrence interval and concluded that the recurrence intervals based on the seismic and geodetic data analyses are mutually corresponding. Based on the small repeating earthquakes, Uchida et al. (2016) estimated the periods of slow slip events that occur off Fukushima Prefecture (corresponding to the swaths between N37° and N38° in this study) at shallow and intermediate depths as 3–5 yr (Uchida et al.2016, Fig. 4). Taking this into account, the dissipation of the interplate coupling in the area south of the main rupture of the 2011 Tohoku-oki earthquake did not progress monotonically, but rather episodically, while expanding northward. This implies that the process of earthquake generation on the subducting plate interface is strongly affected by the background periodic changes in the interplate coupling, and such changes can be monitored not only based on the small repeating earthquakes but also on temporal changes in the displacement rate gradients. The displacement rate distribution can, in principle, be estimated for a time-window shorter than 0.5 yr. However, as shown in Fig. 6 and Fig. S3 in the Supporting Information, since the horizontal displacement rate gradients estimated for a 0.5-yr sub-time-window and vertical displacement rate gradients estimated for 0.5- and 1-yr sub-time-windows have large estimation errors. And as shown in Fig. 10 and Fig. S6 in the Supporting Information large perturbations in a very short period of time, the displacement rates for sub-time-windows shorter than 0.5 yr are not robust. On the other hand, with respect to the width of the swath, 40 km is large enough to robustly estimate displacement rate gradients for most profiles, except those near the ends of Honshu and Hokkaido Islands (Fig. 6 and Fig. S3, Supporting Information). The width of the swaths corresponds to the monitoring resolution of the change in the interplate coupling along the trench-parallel direction. Although the most appropriate width must depend on the density of the GPS observation stations, we can apply the method to other subduction zones where continuous GPS observation networks are sufficiently dense, such as Cascadia and Hikurangi. The density and expansion of the observation network in both trench-normal and trench-parallel directions determine the availability of the proposed method to other subduction zones. It is acceptable to use campaign sites, if we permit that the temporal resolution becomes lower. The graphical expression introduced in this study may make it easy to grasp the spatio-temporal change in the interplate coupling. We cannot estimate the detailed distribution of the area where the interplate coupling changes based only on the change in the displacement rate gradient, especially along the trench-normal direction. However, the recent development of seafloor observation networks equipped with ocean-bottom pressure gauges (OBPGs) such as DONET (Kaneda et al.2015; Kawaguchi et al.2015) and S-net (Uehira et al.2012; Kanazawa 2013) may lead to further applications of the monitoring method based on the displacement rate gradients. We may infer the relative vertical displacement rates in the seafloor observation networks from the OBPG records, even though the sensor drift of OBPGs makes the estimation of the absolute vertical displacement at each OBPG site difficult. 6 CONCLUSIONS A general framework of monitoring method for spatio-temporal change in the interplate coupling based on the spatial gradients of surface displacement rate fields is proposed and applied to NE Japan subduction zone. The results show that we can detect known interplate events in the NE Japan subduction zone, such as the post-seismic slip of the 2003 M8.0 Tokachi-oki and 2005 M7.2 Miyagi-oki earthquakes and the recovery of the interplate coupling around the rupture area of the 1994 M7.6 Sanriku-Haruka-oki earthquake. The results also imply the semi-periodic occurrence of slow slip events and the expansion of the area of slow slip events before the 2011 Tohoku-oki earthquake (M9.0) approaching the hypocentre of the Tohoku-oki earthquake. The proposed method can be applied to other subduction zones and to the operational monitoring after detailed examinations for temporal and spatial estimation error distributions and optimization of the parameters that depend on the site density of the observation network and geometry between the subducting plate interface and the land. Acknowledgements Discussions with Prof. Akira Hasegawa, Prof. Toru Matsuzawa, Prof. Ryota Hino and Dr. Naoki Uchida were valuable in the development of the monitoring method introduced in this paper. I thank two anonymous reviewers and the editor, Prof. Duncan Agnew, for their constructive comments that greatly improved this manuscript. I used Generic Mapping Tools devised by Wessel & Smith (1998) to generate the figures in this paper. This research was partly supported by JSPS KAKENHI grant number JP17H05422 and the project ‘Evaluation and disaster prevention research for the coming Tokai, Tonankai and Nankai earthquakes’ of the Ministry of Education, Culture, Sports, Science and Technology of Japan. REFERENCES Agnew D.C., 1992. The time-domain behavior of power-law noises, Geophys. Res. Lett. , 19, 333– 336. https://doi.org/10.1029/91GL02832 Google Scholar CrossRef Search ADS   Aoki Y., Scholz C.H., 2003. Vertical deformation of the Japanese islands, 1996–1999, J. geophys. Res. , 108, 2257, doi:10.1029/2002JB002129. Crowell B.W., Bock Y., Liu Z., 2016. Single-station automated detection of transient deformation in GPS time series with the relative strength index: a case study of Cascadian slow slip, J. geophys. Res. , 121, 9077– 9094. https://doi.org/10.1002/2016JB013542 Google Scholar CrossRef Search ADS   El-Fiky G.S., Kato T., 1999. Interplate coupling in the Tohoku district, Japan, deduced from geodetic data inversion, J. geophys. Res. , 104, 20,361– 20,377. https://doi.org/10.1029/1999JB900202 Google Scholar CrossRef Search ADS   Freed A.M., Hashima A., Becker T.W., Okaya D.A. Sato H., Hatanaka T., 2017. Resolving depth-dependent subduction zone viscosity and afterslip from postseismic displacements following the 2011 Tohoku-oki, Japan earthquake, Earth planet. Sci. Lett. , 459, 279– 290. https://doi.org/10.1016/j.epsl.2016.11.040 Google Scholar CrossRef Search ADS   Hashimoto C., Noda A., Sagiya T., Matsu’ura M., 2009. Interplate seismogenic zones along the Kuril-Japan trench inferred from GPS data inversion, Nat. Geosci. , 2, 141– 144. https://doi.org/10.1038/ngeo421 Google Scholar CrossRef Search ADS   Holt W.E., Shcherbenko G., 2013. Toward a continuous monitoring of the horizontal displacement gradient tensor field in Southern California using cGPS observations from Plate Boundary Observatory (PBO), Seismol. Res. Lett. , 84, 455– 467. https://doi.org/10.1785/0220130004 Google Scholar CrossRef Search ADS   Hyndman R.D., Wang K., 1995. The rupture zone of Cascadia great earthquakes from current deformation and the thermal regime, J. geophys. Res. , 100, 22 133– 22 154. https://doi.org/10.1029/95JB01970 Google Scholar CrossRef Search ADS   Iinuma T. et al.  , 2012. Coseismic slip distribution of the 2011 off the Pacific Coast of Tohoku Earthquake (M9.0) refined by means of seafloor geodetic data, J. geophys. Res. , 117, B07409, doi:10.1029/2012JB009186. https://doi.org/10.1029/2012JB009186 Google Scholar CrossRef Search ADS   Ito T., Yoshioka S., Miyazaki S., 2000. Interplate coupling in northeast Japan deduced from inversion analysis of GPS data, Earth planet. Sci. Lett. , 176, 117– 130. https://doi.org/10.1016/S0012-821X(99)00316-7 Google Scholar CrossRef Search ADS   Kanazawa T., 2013. Japan trench earthquake and tsunami monitoring network of cable-linked 150 ocean bottom observatories and its impact to Earth disaster science, in 2013 IEEE International Underwater Technology Symposium (UT) , pp. 1– 5, IEEE, Tokyo. Kaneda Y. et al.  , 2015. Development and application of an advanced ocean floor network system for megathrust earthquakes and tsunamis, in Seafloor Observatories , pp. 643– 662, eds Favali P. et al.  , Springer Praxis Books. Google Scholar CrossRef Search ADS   Kato T., Tsumura K., 1979. Vertical land movement in Japan as deduced from tidal record (1951–1978) (in Japanese), Bull. Earthq. Res. Inst., Univ. Tokyo , 54, 559– 628. Kawaguchi K., Kaneko S., Nishida T., $ Komine T., 2015. Construction of the DONET real-time seafloor observatory for earthquakes and tsunami monitoring, in Seafloor Observatories , pp. 211– 228, eds Favali P. et al.  , Springer Praxis Books. Google Scholar CrossRef Search ADS   Kikuchi M., 2001. Fault model: Earthquake source processes (in Japanese), in Encyclopedia of Earthquakes , 2nd edn, pp. 259– 284, eds Utsu T., Shima E., Yoshii T., Yamashina K., Asakura Publishing, Tokyo. Kita S., Okada T., Hasegawa A., Nakajima J., Matsuzawa T., 2010, Anomalous deepening of a seismic belt in the upper-plane of the double seismic zone in the Pacific slab beneath the Hokkaido corner: possible evidence for thermal shielding caused by subducted forearc crust materials, Earth planet. Sci. Lett. , 290, 415– 426. https://doi.org/10.1016/j.epsl.2009.12.038 Google Scholar CrossRef Search ADS   Loveless J.P., Meade B.J., 2010. Geodetic imaging of plate motions, slip rates, and partitioning of deformation in Japan, J. geophys. Res. , 115, B02410, doi:10.1029/2008JB006248. https://doi.org/10.1029/2008JB006248 Google Scholar CrossRef Search ADS   Lundgren P., Protti M., Donnellan A., Heflin M., Hernandez E., Jefferson D., 1999. Seismic cycle and plate margin deformation in Costa Rica: GPS observations from 1994 to 1997, J. geophys. Res. , 104, 28 915– 28 926. https://doi.org/10.1029/1999JB900283 Google Scholar CrossRef Search ADS   Matsu’ura M., Noda A., Fukahata Y., 2007. Geodetic data inversion based on Bayesian formulation with direct and indirect prior information, Geophys. J. Int. , 171, 1342– 1351. https://doi.org/10.1111/j.1365-246X.2007.03578.x Google Scholar CrossRef Search ADS   Mazzotti S., Le Pichon X., Henry P., Miyazaki S., 2000. Full interseismic locking of the Nankai and Japan-west Kurile subduction zones: an analysis of uniform elastic strain accumulation in Japan constrained by permanent GPS, J. geophys. Res. , 105, 13 159–13 177. https://doi.org/10.1029/2000JB900060 Mavrommatis A.P., Segall P., Johnson K.M., 2014. A decadal-scale deformation transient prior to the 2011 Mw 9.0 Tohoku-oki earthquake, Geophys. Res. Lett. , 41, 4486– 4494. https://doi.org/10.1002/2014GL060139 Google Scholar CrossRef Search ADS   Meade B.J., 2007, Algorithms for the calculation of exact displacements, strains, and stresses for triangular dislocation elements in a uniform elastic half space, Comput. Geosci. , 33, 1064– 1075. https://doi.org/10.1016/j.cageo.2006.12.003 Google Scholar CrossRef Search ADS   Meade B.J., Loveless J.P., 2009. Block modeling with connected fault-network geometries and a linear elastic coupling estimator in spherical coordinates, Bull. seism. Soc. Am. , 99, 3124– 3139. https://doi.org/10.1785/0120090088 Google Scholar CrossRef Search ADS   Miura S., Iinuma T., Yui S., Uchida N., Sato T., Tachibana K., Hasegawa A., 2006 Co- and post-seismic slip associated with the 2005 Miyagi-oki earthquake (M7.2) as inferred from GPS data, Earth Planets Space , 58, 1567– 1572. https://doi.org/10.1186/BF03352662 Google Scholar CrossRef Search ADS   Miyazaki S., Segall P., Fukuda J., Kato T., 2004. Space time distribution of afterslip following the 2003 Tokachi-oki earthquake: implications for variations in fault zone frictional properties, Geophys. Res. Lett. , 31, L06623, doi:10.1029/2003GL019410. Mosegaard K., Tarantola A., 1995. Monte Carlo sampling of solutions to inverse problems, J. geophys. Res. , 100, 12 431– 12 447. https://doi.org/10.1029/94JB03097 Google Scholar CrossRef Search ADS   Murakami M., Suito H., Ozawa S., Kaidzu M., 2006. Earthquake triggering by migrating slow slip initiated by M8 earthquake along Kuril Trench, Japan, Geophys. Res. Lett. , 33, L09306, doi:10.1029/2006GL025967. https://doi.org/10.1029/2006GL025967 Google Scholar CrossRef Search ADS   Murotani S., 2003. Rupture processes of large Fukushima-Oki earthquakes in 1938, Master’s thesis , University of Tokyo, Tokyo. Nakagawa H.et al.  , 2009. Development and validation of GEONET new analysis strategy (Version 4) (in Japanese), J. Geogr. Surv. Inst. , 118, 1– 8. Nakajima J., Hasegawa A., 2006. Anomalous low-velocity zone and linear alignment of seismicity along it in the subducted Pacific slab beneath Kanto, Japan: reactivation of subducted fracture zone?, Geophys. Res. Lett. , 33, L16309, doi:10.1029/2006GL026773. https://doi.org/10.1029/2006GL026773 Google Scholar CrossRef Search ADS   Nakajima J., Hirose F., Hasegawa A., 2009. Seismotectonics beneath the Tokyo metropolitan area, Japan: effect of slab-slab contact and overlap on seismicity, J. geophys. Res. , 114, B08309, doi:10.1029/2008JB006101. https://doi.org/10.1029/2008JB006101 Google Scholar CrossRef Search ADS   Nishimura T., Hirasawa T., Miyazaki S., Sagiya T., Tada T., Miura S., Tanaka K., 2004. Temporal change of interplate coupling in northeastern Japan during 1995–2002 estimated from continuous GPS observations, Geophys. J. Int. , 157, 901– 916. https://doi.org/10.1111/j.1365-246X.2004.02159.x Google Scholar CrossRef Search ADS   Norabuena E., Leffler-Griffin L., Mao A., Dixon T., Stein S., Sacks I.S., Ocola L., Ellis M., 1998. Space geodetic observations of Nazca-South America convergence across the Central Andes, Science , 279, 358– 362. https://doi.org/10.1126/science.279.5349.358 Google Scholar CrossRef Search ADS PubMed  Ohtani R., McGuire J.J., Segall P., 2010. Network strain filter: a new tool for monitoring and detecting transient deformation signals in GPS arrays, J. geophys. Res. , 115, B12418, doi:10.1029/2010JB007442. https://doi.org/10.1029/2010JB007442 Google Scholar CrossRef Search ADS   Okada Y., 1992. Internal deformation due to shear and tensile faults in a half-space, Bull. seism. Soc. Am. , 82, 1018– 1040. Okada T., Hasegawa A., 2003. The M7.1 May 26, 2003 off-shore Miyagi Prefecture Earthquake in northeast Japan: source process and aftershock distribution of an intra-slab event, Earth Planets Space , 55, 731– 739. https://doi.org/10.1186/BF03352482 Google Scholar CrossRef Search ADS   Ozawa S., Suito H., Nishimura T., Tobita M., Munekane H., 2007. Possibility of recovery of slip deficit rate between the North American plate and the Pacific plate off Sanriku, northeast Japan, Geophys. Res. Lett. , 34, L20308, doi:10.1029/2007GL030477. https://doi.org/10.1029/2007GL030477 Google Scholar CrossRef Search ADS   Ozawa S., Nishimura T., Munekane H., Suito H., Kobayashi T., Tobita M., Imakiire T., 2012. Preceding, coseismic, and postseismic slips of the 2011 Tohoku earthquake, Japan, J. geophys. Res. , 117, B07404, doi:10.1029/2011JB009120. https://doi.org/10.1029/2011JB009120 Google Scholar CrossRef Search ADS   Savage J.C., 1983. A dislocation model of strain accumulation and release at a subduction zone, J. geophys. Res. , 88, 4984– 4996. https://doi.org/10.1029/JB088iB06p04984 Google Scholar CrossRef Search ADS   Segall P., 2010. Earthquake and Volcano Deformation , pp. 51– 91, Princeton University Press. Sun T.et al.  , 2014. Prevalence of viscoelastic relaxation after the 2011 Tohoku-oki earthquake, Nature , 514, 84– 87. https://doi.org/10.1038/nature13778 Google Scholar CrossRef Search ADS PubMed  Suwa Y., Miura S., Hasegawa A., Sato T., Tachibana K., 2006. Interplate coupling beneath NE Japan inferred from three-dimensional displacement field, J. geophys. Res. , 111, B04402, doi:10.1029/2004JB003203. Uchida N., Iinuma T., Nadeau R.M., Bürgmann R., Hino R, 2016. Periodic slow slip triggers megathrust zone earthquakes in northeastern Japan, Science , 351, 488– 492. https://doi.org/10.1126/science.aad3108 Google Scholar CrossRef Search ADS PubMed  Uehira K. et al.  , 2012. Ocean bottom seismic and tsunami network along the Japan trench, in Abstract OS41C-1736 presented at 2012 Fall Meeting , AGU, San Francisco, California, 3–7 December. Wallace L.M., Beavan J., McCaffrey R., Darby D., 2004. Subduction zone coupling and tectonic block rotations in the North Island, New Zealand, J. geophys. Res. , 109, B12406, doi:10.1029/2004JB003241. https://doi.org/10.1029/2004JB003241 Google Scholar CrossRef Search ADS   Watanabe S., Sato M., Fujita M., Ishikawa T., Yokota Y., Ujihara N., Asada A., 2014. Evidence of viscoelastic deformation following the 2011 Tohoku-Oki earthquake revealed from seafloor geodetic observation, Geophys. Res. Lett. , 41, 5789– 5796. https://doi.org/10.1002/2014GL061134 Google Scholar CrossRef Search ADS   Wessel P., Smith W.H. F., 1998. New, improved version of Generic Mapping Tools released, EOS, Trans. Am. geophys. Un. , 79, 579, doi:10.1029/98EO00426. https://doi.org/10.1029/98EO00426 Google Scholar CrossRef Search ADS   Williams S.D.P., Bock Y., Fang P., Jamason P., Nikolaidis R.M., Prawirodirdjo L., Miller M., Johnson D.J., 2004. Error analysis of continuous GPS position time series, J. geophys. Res. , 109, B03412, doi:10.1029/2003JB002741. Yabuki T., Matsu’ura M., 1992. Geodetic data inversion using a Bayesian information criterion for spatial distribution of fault slip, Geophys. J. Int. , 109, 363– 375. https://doi.org/10.1111/j.1365-246X.1992.tb00102.x Google Scholar CrossRef Search ADS   Yamanaka Y., 2003. Off Fukushima-ken earthquake, 31 October 2003, EIC Seismological Notes , 141. Available at http://www.eri.u-tokyo.ac.jp/sanchu/Seismo_Note/EIC_News/031031.html, last accessed 31 May 2017. Yamanaka Y., Kikuchi M., 2003. Source process of the recurrent Tokachi-oki earthquake on September 26, 2003, inferred from teleseismic body waves, Earth Planets Space , 55, e21– e24. https://doi.org/10.1186/BF03352479 Google Scholar CrossRef Search ADS   Yamanaka Y., Kikuchi M., 2004. Asperity map along the subduction zone in northeastern Japan inferred from regional seismic data, J. geophys. Res. , 109, B07307, doi:10.1029/2003JB002683. https://doi.org/10.1029/2003JB002683 Google Scholar CrossRef Search ADS   Yamanaka Y., 2005. Asperity Map along the subduction zone in Hokkaido region (in Japanese), in Abstract A055 Presented at 2005 Fall Meeting , Seism. Soc. Jpn., Sapporo, Japan, 19–21 October. Yokota Y., Koketsu K., 2015. A very long-term transient event preceding the 2011 Tohoku earthquake, Nature Commun. , 6, 5934, doi:10.1038/ncomms6934. https://doi.org/10.1038/ncomms6934 Google Scholar CrossRef Search ADS   SUPPORTING INFORMATION Supplementary data are available at GJI online. Figure S1. Examples of the displacement rate field for 0.5-year, 1-year, 2-year and 3-year sub-timewindows from 7 March 2011 (see Figure 3 for 5-year time-window). The displacement rate profiles in the swath for which the centerline runs through N38.7°, E140.0° are also shown. Figure S2. All profiles of an example displacement rate field for a 5-year time-window from 7 March 2006 to 7 March 2011. The width of the swaths is configured to be 60 km. The red vectors indicate the horizontal displacement rates at the sites in the swaths. The bottom panel shows the profile of the horizontal and vertical displacement rates along the direction of the swath with respect to the distance from the trench. Figure S3. Trench-parallel distribution of the average estimation error of the displacement rate gradients. The upper panel shows the vertical component, whereas the lower panel shows the horizontal component. All cases for the width of the swaths are plotted using various styles of lines (see the legends in the panels). Orange and green lines correspond with 1-year and 3-year sub -timewindows. Figure S4. Spatial distributions of the difference between displacement rate gradients using the synthetic displacement rate fields calculated from the assumed interplate coupling distribution shown in Figure 7(c) and the interplate coupling distribution with perturbation on plate interface shallower than 30 km in depth. The widths of the swaths (D) and perturbed regions (W) are implemented in the panels. Figure S5. Spatial distributions of the difference between displacement rate gradients using the synthetic displacement rate fields calculated from the assumed interplate coupling distribution shown in Figure 7(c) and the interplate coupling distribution with perturbation on plate interface at 50 to 90 km in depth. The widths of the swaths (D) and perturbed regions (W) are implemented in the panels. Figure S6. Spatio-temporal variations of the displacement rate gradients for all cases. The horizontal axis indicates the latitudes of the crossing point of the centerline of the profile region and meridian E140°. The panels on the left-hand side correspond to the Kanto to Tohoku regions, and the panels on the right-hand side correspond with the Hokkaido region. The vertical axis indicates the final dates of the sub-time-windows and the time-windows. The top and bottom panels represent the vertical and horizontal displacement rate gradients. Note that different color scales are used for horizontal and vertical components. Figure S7. Time series of the horizontal displacement rate gradients at N36.0° to N39.0° before the 2011 Tohoku-oki earthquake for all cases. The data for each profile region are plotted with a color corresponding to the profile latitude at E140.0°, as shown at the right-hand side of each panel. The horizontal line identified by the latitude number corresponds with the zero displacement rate gradient for each profile. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. APPENDIX A: SELECTION OF EARTHQUAKES IN THE DISPLACEMENT TIME-SERIES ANALYSIS The earthquake catalogue of the Japan Meteorological Agency for the period of 1996–2016 was used to configure the dates for the step functions in the time-series analysis. First, M earthquakes of magnitude greater than 4.0 were listed. Next, dates for which the horizontal or vertical displacement was greater than the threshold were sought for each site, namely, $$\sqrt{{dE^i_k}^2 + {dN^i_k}^2 } > 2$$ mm or $$|dU^i_k| > 5$$ mm. Here, $$dE^i_k$$, $$dN^i_k$$ and $$dU^i_k$$ are the east–west, north–south and up–down components of the displacement at the ith site due to the earthquake(s) that occurred on tk. These displacement components are calculated as:   \begin{eqnarray} dE^i_k = \sum _{m=1}^M \delta _{kD(m)} de^i_{m} \mbox{,} \end{eqnarray} (A1)  \begin{eqnarray} dN^i_k = \sum _{m=1}^M \delta _{kD(m)} dn^i_{m} \quad \mbox{and} \end{eqnarray} (A2)  \begin{eqnarray} dU^i_k = \sum _{m=1}^M \delta _{kD(m)} du^i_{m} \mbox{,} \end{eqnarray} (A3)where δkj and D(m) are the Kronecker Delta and the occurrence date of the mth earthquake, respectively, $$de^i_m$$, $$dn^i_m$$ and $$du^i_m$$ are the east–west, north–south and up–down components, respectively, of the displacement at the ith site due to the mth earthquake. The displacement at a GPS station due to an earthquake was calculated as follows. Estimating the area of the earthquake source fault based on its magnitude by applying the equation of the stress drop, Δσ, the earthquake moment, Mo and the fault area, S, derived from the dislocation theory under the assumption that the fault length is twice the width, i.e., $$\Delta \sigma = 2.5 M_o S^{-\frac{3}{2}}$$ (e.g. Kikuchi 2001). Δσ = 10 MPa was assumed, and Mo was calculated from its magnitude by regarding it to be the same as the moment magnitude. Calculating the fault length, L, from estimated fault area, S, as $$L=\sqrt{2S}$$. Estimating the fault slip amount from Mo and S by assuming the rigidity to be 40 GPa. Calculating the displacement at the ith site by means of Okada’s (1992) method. The following two fault mechanisms are assumed: Pure strike-slip faulting on a vertical fault along the direction rotated by 45° from the line between the epicentre of the mth earthquake and the ith GPS site. Pure thrust faulting on a dipping (30°) fault, the strike of which is perpendicular to the line between the epicentre of the mth earthquake and the ith GPS site. $$de^i_m$$, $$dn^i_m$$ and $$du^i_m$$ were obtained by taking the larger of the absolute values between the two cases for each component. The procedure may be improved by introducing centroid moment tensor (CMT) solutions from a catalogue, such as the Global CMT catalogue, but displacement steps at each GPS site should be estimated from the displacement time-series data. Therefore, the estimated displacement rates do not change if the dates for step functions are not changed. APPENDIX B: CALCULATION OF THE DISPLACEMENT RATE GRADIENT The spatial gradients of the displacement rate field are basically calculated by applying the weighted least-squares linear regression. Thus, the spatial gradients of horizontal and vertical displacement rates, gh and gv, are calculated as:   \begin{eqnarray} g_h &=& \frac{\displaystyle \left[ \sum _{n=1}^N \frac{1}{\left( {\delta a_h}_n\right)^2} \right] \left[ \sum _{n=1}^N \frac{x_n {a_h}_n}{ \left({\delta a_h}_n \right) ^2} \right] - \left[ \sum _{n=1}^N \frac{x_n}{\left( {\delta a_h}_n\right)^2} \right] \left[ \sum _{n=1}^N \frac{{a_h}_n}{ \left({\delta a_h}_n \right) ^2} \right] }{\displaystyle \left[ \sum _{n=1}^N \frac{1}{\left( {\delta a_h}_n\right)^2} \right] \left[ \sum _{n=1}^N \frac{x_n^2}{ \left({\delta a_h}_n \right) ^2} \right] - \left[ \sum _{n=1}^N \frac{x_n}{\left( {\delta a_h}_n\right)^2} \right]^2}, \quad \mbox{and} \end{eqnarray} (B1)  \begin{eqnarray} g_v &=& \frac{\displaystyle \left[ \sum _{n=1}^N \frac{1}{\left( {\delta a_v}_n\right)^2} \right] \left[ \sum _{n=1}^N \frac{x_n {a_v}_n}{ \left({\delta a_v}_n \right) ^2} \right] - \left[ \sum _{n=1}^N \frac{x_n}{\left( {\delta a_v}_n\right)^2} \right] \left[ \sum _{n=1}^N \frac{{a_v}_n}{ \left({\delta a_v}_n \right) ^2} \right] }{\displaystyle \left[ \sum _{n=1}^N \frac{1}{\left( {\delta a_v}_n\right)^2} \right] \left[ \sum _{n=1}^N \frac{x_n^2}{ \left({\delta a_v}_n \right) ^2} \right] - \left[ \sum _{n=1}^N \frac{x_n}{\left( {\delta a_v}_n\right)^2} \right]^2}, \end{eqnarray} (B2)where ahn, avn, δahn, δavn and xn are the horizontal and vertical displacement rates at the nth site in a swath, their estimation errors, and the distance of the nth site in the profile direction from the crossing point of the profile centreline and meridian E140°, respectively, and N is the number of sites in the swath. In practice, the displacement rate gradients and their errors, εh and εv, respectively, are estimated by different procedures depending on N, as follows: When N = 2, gh and gv are calculated using eqs (B1) and (B2). Moreover, εh and εv are calculated as $$\displaystyle \varepsilon _h = \sqrt{{\delta a_h}_1^2 + {\delta a_h}_1^2}$$, and $$\varepsilon _v = \sqrt{{\delta a_v}_1^2 + {\delta a_v}_1^2}$$. When 2 < N ≤ 10, gh and gv are calculated by taking the averages of ghn and gvn, where ghn and gvn (here, n = 1, …, N) are the spatial gradients calculated from the displacement rates at the sites in the swath, excluding the nth site. Moreover, εh and εv are defined as the standard deviations of ghn and gvn, respectively. When 10 < N ≤ 20, gh and gv are calculated by taking the averages of ghmn and gvmn, where ghmn and gvmn (here, n = 1, …, N − 1 and m = n + 1, …, N) are the spatial gradients calculated from the displacement rates at the sites in the swath, excluding the mth and nth sites. Here, εh and εv are the standard deviations of ghmn and gvmn, respectively. When 20 < N, gh and gv are calculated by taking the averages of ghlmn and gvlmn, respectively, where ghlmn and gvlmn (here, n = 1, …, N − 2, m = n + 1, …, N − 1 and l = m + 1, …, N) are the spatial gradients calculated from the displacement rates at the sites in the swath, excluding the lth, mth and nth sites. Here, εh and εv are defined as the standard deviations of ghlmn and gvlmn, respectively. Since displacement rates at GPS sites sometimes have systematic errors larger than the estimation errors based on the least-squares fitting of the displacement time-series, I adopted the above procedure to estimate the spatial gradients and their errors. © The Author(s) 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.

Journal

Geophysical Journal InternationalOxford University Press

Published: Apr 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off