Fourteen‐Year Acceleration Along the Japan Trench

An acceleration of the background seismicity and a shortening of the slow slip events recurrence intervals in the Boso peninsula (Japan) suggest a slow unlocking of the Philippine Sea‐North America (PHS‐NAM) subduction interface from 1990 to 2011. Motivated by these previous observations, we used GPS (Global Positioning System) time series to study the 1997–2011 evolution of interface locking offshore Honshu with a specific focus on the Kanto region. We newly processed the GPS data in double difference and analyzed them with a trajectory model that accounts for seismic and aseismic variations, and that includes an inter‐seismic acceleration term. We inverted the surface acceleration obtained, on both the Pacific‐North America (PAC‐NAM) and the PHS‐NAM interfaces. The inverted slip rate changes over time compare well with previous studies: we infer slip deceleration between 39° – 41° N and acceleration between 37° – 38° N (PAC‐NAM), with a maximum amplitude of 1.75 mm/ yr2 corresponding to an equivalent geodetic locking change of 0.3 over the 14 years. More notably, our analysis reveals a novel and robust slip acceleration South of 36° N that we interpret as an unlocking of the PAC‐NAM interface. It is located noticeably far from the 2011 Tohoku earthquake rupture and is therefore unlikely connected to it. The slip acceleration inverted from the surface displacements is comparable to the changes observed in the seismicity rate. Our results further demonstrate that inter‐seismic slip rate can significantly evolve over years to decades, and suggest a simple relationship between the seismicity and the slip on the subduction interface.


Introduction
A common assumption regarding the seismic cycle is to consider the inter-seismic strain rate as being constant over time (Savage & Thatcher, 1992). Recently, long-term changes in slip rate, or interface locking, have been inferred or suggested in the context of subduction zones. Long-term slow slip events (L-SSEs) with a duration of a few years have been documented in Alaska (duration of 2-9 years) (Li et al., 2016;Rousset et al., 2019) and in various regions of Japan (e.g., Tokai district, Kanto region, Kii peninsula, and Bungo channel, with duration from 1 to 5 years) (Hirose et al., 1999;Kobayashi, 2014;Kobayashi & Tsuyuki, 2019;Miyazaki et al., 2003;Ochi & Kato, 2013;Ozawa et al., 2001Ozawa et al., , 2013Ozawa, Suito, Imakiire, & Murakmi, 2007;Tanaka & Yabe, 2017;Yagi & Kikuchi, 2003;Yoshioka et al., 2015). Variations at even longer time scales (decades) have also been inferred. In Sumatra, Prawirodirdjo et al. (2010) infer a decrease of locking between the 1990s and 2010 in the Batu and Enggano islands. Based on coral observations that Abstract An acceleration of the background seismicity and a shortening of the slow slip events recurrence intervals in the Boso peninsula (Japan) suggest a slow unlocking of the Philippine Sea-North America (PHS-NAM) subduction interface from 1990 to 2011. Motivated by these previous observations, we used GPS (Global Positioning System) time series to study the 1997-2011 evolution of interface locking offshore Honshu with a specific focus on the Kanto region. We newly processed the GPS data in double difference and analyzed them with a trajectory model that accounts for seismic and aseismic variations, and that includes an inter-seismic acceleration term. We inverted the surface acceleration obtained, on both the Pacific-North America (PAC-NAM) and the PHS-NAM interfaces. The inverted slip rate changes over time compare well with previous studies: we infer slip deceleration between  39 E - 41 E N and acceleration between  37 E - 38 E N (PAC-NAM), with a maximum amplitude of 1.75 mm/ 2 yr E corresponding to an equivalent geodetic locking change of 0.3 over the 14 years. More notably, our analysis reveals a novel and robust slip acceleration South of  36 E N that we interpret as an unlocking of the PAC-NAM interface. It is located noticeably far from the 2011 Tohoku earthquake rupture and is therefore unlikely connected to it. The slip acceleration inverted from the surface displacements is comparable to the changes observed in the seismicity rate. Our results further demonstrate that inter-seismic slip rate can significantly evolve over years to decades, and suggest a simple relationship between the seismicity and the slip on the subduction interface.
Plain Language Summary Subduction zones, where an oceanic plate dives under a continental plate, host both seismic and aseismic slip events. In Japan, along the Sagami Trough where the Philippine Sea (PHS) plate subducts beneath the North America (NAM) plate, the increase of the number of earthquakes between 1990 and 2011 suggests that the PHS plate slips at a higher rate under the NAM plate. Along the Japan Trench, North of the Sagami Trough, the Pacific (PAC) plate subducts beneath the NAM plate; changes in the slip rate of the PAC plate have already been inferred from geodetic measurements by previous studies. To further investigate those changes in slip rate at both interfaces, we processed and analyzed GPS (Global Positioning System) time series from 1997 to 2011, and compared them with changes in earthquake rates. From those time series, we determine the slip rate of the continental plate. On the PAC-NAM interface, the slip either decelerates or accelerates depending on the location. Most notably, our analysis reveals a novel slip acceleration near the Boso peninsula that we interpret as an acceleration of the PAC plate's slip beneath the NAM plate, proving that plates can subduct at non-stationary velocities.
We here present a study of inter-seismic slip rate changes at the Japan Trench-North America (NAM) plate and Pacific (PAC) plate interface-and the Sagami Through-NAM plate and Philippine Sea (PHS) plate interface-subduction zones (Figure 1), from January 1, 1997 to February 6, 2011. As large scale variations in geodetic movements can be caused by reference frame issues, a first goal of this study is to provide a full re-analysis of the Japanese GPS data that is completely independent of the F3 solution (Nakagawa et al., 2009) on which are based previous studies of decadal slip rate variations (Mavrommatis et al., 2014;Yokota & Koketsu, 2015). The second goal is to complement the observations of Mavrommatis et al. (2014) and Yokota and Koketsu (2015) by extending the investigation of the surface displacements to the Kanto region ( Figure 1). This region is known to have undergone a strong acceleration of the background seismicity rate on both PHS and PAC plates (Marsan et al., 2017;Reverso et al., 2016), as well as a shortening of the Boso SSEs recurrence times from 1996 to 2011 on the PHS plate (Fukuda, 2018;Gardonio et al., 2018;Hirose et al., 2012). These observations suggest that plate unlocking, involving the three tectonic plates that control surface deformation in Kanto (the NAM, PAC, and PHS plates), has been going on for at least a decade before the 2011 w E M 9.0 Tohoku earthquake. Given the distance between the Kanto region and potential mainshock candidates that could possibly explain the northern Honshu changes in slip rate (namely, the 2003 w E M 8.0 Tokachi earthquake, see Figure 1), any variation in local slip rate in Kanto could not be conceivably related to unmodeled post-seismic processes of such mainshocks. The Kanto region is therefore a natural candidate to study possible decadal slip rate variations. The third goal of this study is to compare the evolution of the slip on the subduction interface derived from surface displacements to the variations of the seismicity rate, and to propose a simple relationship that relates both observables.
In more detail, this study complements the analyses of Mavrommatis et al. (2014) and Yokota and Koketsu (2015) by (a) performing a re-analysis of the GPS data, (b) including the Kanto region, and (c) testing the sensitivity of the inverted changes in slip rate relative to the inclusion of slip on the PHS plate, of different a priori distributions, including a seismicity-based prior, and of GPS vertical components. After summarizing the tectonic setting of the Honshu island (Section 2), we describe our processing and analysis of GPS data (Section 3) as well as the inversion of slip rate distribution (Section 4). We then analyze the mean surface velocity and the acceleration fields, and infer the inter-seismic locking, the locking difference between 2011

Data Processing
We processed the 1,421 stations of the GNSS (Global Navigation Satellite System) Earth Observation Network System in Japan (GEONET) and 44 IGS (International GNSS Service, https://www.igs.org/) sites worldwide, from January 1, 1997 to February 6, 2011 using only GPS (Global Positioning System) satellite system and following a double-difference approach, using the GAMIT/GLOBK software suite (Herring et al., 2015(Herring et al., , 2018a(Herring et al., , 2018b. As in Herring et al. (2016), we assemble our stations into sub-networks for the daily processing. We reduce 24-hr measurement sessions to daily estimates of station position, choosing the ionosphere-free combination and fixing the ambiguities to integer values. We use precise orbit positions  Note. 1998 Iwate volcanic unrest from Miura et al. (2000); 2000 Miyakejima volcanic unrest from Cattania et al. (2017), Nakada et al. (2005) and Uhira et al. (2005); Boso slow slip events (SSEs) from Fukuda (2018) and Gardonio et al. (2018); Tohoku SSE from Y. Ito et al. (2013). a Computed from the total moment release:   19 0 3.6 10 E M N m (Cattania et al., 2017). b Estimation based on the observed time series displacements.
We first proceed to a local estimate for the jumps affecting the time series. The 400-data point window allows to have a good estimate of the jumps without polluting the estimation with long-term effects. To do so, we go through the time series and take the first jump (antenna change or co-seismic) or slow deformation (hereafter considered as a "jump" also); if there is no jump in the whole time series, the program proceeds to the next step. We consider first a 400-data point window centered on the jump date ( s E t ). If another jump happens between the first jump and the end of the window, we extend the window with an extra 200 data points after the second jump. We keep extending the window this way until no more jump stays in the last 200 data points of the window. Whenever fewer than the required number of data points are available (at the beginning or end of the analyzed period), we take as many data points as possible. To illustrate this process, Figure 2c shows a window with three jumps (the 2000 Miyakejima volcanic unrest and two antenna changes). Once the window length is fixed, we fit Equation 4 within that window. As given in Equation 4, an inter-seismic velocity E v and post-seismic transients are also modeled for the window, to get the most accurate estimate of the jumps as possible. In case there is a post-seismic transient in the window, we force the transient to have the same sign as the associated co-seismic jump (   0 i j E m c for and E j corresponding to the same earthquake), individually for each component (NS, EW, or vertical). Moreover, in this window, we take R E T = 10 days because after testing different values we observe a better co-seismic jump estimation with a smaller R E T . The value of the co-seismic estimation is important because the estimation of the post-seismic displacement (in the next step) can vary according to the estimation of the co-seismic jump. Then, we go to the next jump not included in the current window and start again, until it reaches the last jump. Figure 2 illustrates how the fitting and modeling perform for some selected time windows and stations. Once we have the a E b , j E c and s E d optimized for all the jumps of the time series, we remove these jumps to obtain which is thus the time series corrected for co-seismic offsets, antenna jumps and slow transients (Table 2).
In the second step, we compute i E m by optimizing the values of the post-seismic transients using the time series from 2 years of data before the earthquake to 2 years of data after and fitting: We take 2-year data before the earthquake to estimate the seasonal and semi-seasonal terms within the window (or less than 2 years of data depending on data availability), which results in a better estimation of the post-seismic transient. We also take 2 years of data (or less if not available) after the earthquake. If there is another post-seismic transient within the 2 years of data, as for the 400-data jump window, we extend the window to 2-year data after the transient, until there is no post-seismic transient during the last 2 years of data. Then, we iterate this step to all subsequent w E M  E 6.8 earthquakes influencing the station and subtract the modeled post-seismic transients from 1 ( ) In the third step, we model x t step2 ( ) as: and s , and k E c by optimizing the seasonal and inter-seismic phenomena using the whole time series. Then we remove these contributions to get the final residuals of our raw time series.
10.1029/2020JB021226 8 of 24 In Equation 1, we added an acceleration term ( E a ) that is usually absent from trajectory models (Bevis & Brown, 2014). This term is motivated by the observation that the residuals of linear trajectory models (  0 E a ) very often display a decadal curvature (visible in the residuals in Figure 3b) and by previous finding of Heki and Mitsui (2013), Mavrommatis et al. (2014), and Yokota and Koketsu (2015) who have first identified such surface motion acceleration. We add a quadratic term, representative of the acceleration, to test if the inter-seismic velocity could have undergone significant changes in Honshu over the 1997-2011 period. The results of the quadratic model are shown in Figures 3c and 3d. To properly estimate the acceleration, E a , we require the time series to have at least 8 years and a half (70%) of data between 1997 and 2011.

Statistical Significance of the Acceleration
We test the significance of adding a quadratic term using synthetic data. To generate synthetic time series, we first compute the Fourier Transform of the quadratic model residuals and randomize the phase. Finally, the inverse Fourier Transform yields a synthetic time series in which colored noise follows the spectral characteristics of the original residual time series. For each station individually, we generate 100 synthetic

Inversion Model of the Loading Rate and Its Acceleration on the Subduction Interface
To determine the average locking and the slip acceleration on the subduction interface, we perform separate least-square slip inversions (Tarantola & Valette, 1982) of the surface velocity and the acceleration fields, respectively, using a modified version of Kositsky and Avouac (2010)'s software package, including the regularization of Radiguet et al. (2011). We only use horizontal displacement time series at this stage; the addition of the vertical displacements will be discussed later (Section 5.2). The forward model is  E d Gm , where E d is the observed data (surface velocity or acceleration field), E G , the transfer function matrix computed using Okada (1985), and E m , the model on the fault (slip rate or its acceleration, depending on the analysis). The best model ( E m ) is determined using the misfit function: where 0 E m is the prior model, m E C is the covariance matrix for the model parameters (detailed in Section 4.4), d E C is the covariance matrix for the data (velocity or acceleration field, detailed in Section 4.3), and denotes the transpose operation. Minimizing ( ) E S m of Equation 7, we obtain the model function:

Plate Geometry and Imposed Rake
As the Kanto region ( Figure 1) is characterized by a double subduction, we run two tests for each inversion (velocity and acceleration fields): (a) considering only the Pacific-North America (PAC-NAM) subduction, and (b) considering both the Philippine Sea-North America (PHS-NAM) and the PAC-NAM subductions. We discretize the faults into triangular sub-fault patches of  E 15 km size. For the PAC plate interface, we use the Kamchatka-Kuril Island-Japan region of the Slab 2 model (Hayes et al., 2018). We keep only the part from the Sagami Trough (   34.2 E N, see Figure 1) to the North of Honshu island (   41 E N), and from the trench down to 90 km depth. For the PHS plate interface, we use Ishida (1992)'s model. For both interfaces, we impose a fixed rake angle for each sub-fault.
To determine this angle, we compute the PAC and PHS plate velocity vectors relative to the NAM plate at each sub-fault using the Euler poles from Nishimura et al. (2007), given in Table 1. Then, we determine the rake, which will be different for each sub-fault, by projecting the velocity vector direction on the sub-fault surface.

Prior Models
Several prior models have been tested. As default models, we consider (a) a fully locked model for the coupling, hence slip rates equal zero, and (b) a zero acceleration model for the slip acceleration. If the prior is not specified, this default model was applied.
We additionally test a fully unlocked prior model, hence slip rates equal to the convergence rate: this prior brings only small differences in the well-restored area of the PAC plate ( Figure S1 in Supporting Information S1). We decided to use the fully locked model as our default model rather than the fully unlocked one because the locking in the not well-restored area is much more likely to be locked. We notice when comparing Figure 4; Figures S2 and S3 in Supporting Information S1, that while changing the PHS plate prior model does not affect the inverted PAC interface locking, the PAC plate prior model impacts the inverted PHS interface locking. Because of this sensitivity to the PAC prior, we must keep in mind that the results for the PHS plate are to be taken with caution.
For the slip acceleration, we test a prior model based on seismicity acceleration to test whether slip and background seismicity are linked. The computation of the seismicity acceleration and the slip acceleration prior model, as well as the inversion results are detailed later in Section 6.

Data Covariance Matrix
The data covariance matrix ( d E C ) is a diagonal matrix whose dimension is two times the number of stations, since there are two horizontal channels per station. The diagonal values are determined by the data uncertainties for the North and East components: where i is the index of the station, sta E n is the number of stations, i E errN , the uncertainty of the data i according to the North component, and i E errE , the uncertainty of the data i according to the East component.
As velocity uncertainties, we take the RMS of the quadratic model and divide it by the number of data to take into account the data gap. The 5% and 95% quantiles of the horizontal velocity uncertainties are found at 0.14 and 0.27 mm/yr, respectively, with a mean of 0.19 mm/yr. As acceleration uncertainties, we take the standard deviation of the acceleration computed from the synthetic time series (see Section 3.3). The 5% and 95% quantiles of the horizontal acceleration uncertainties are found at 0.010 and 0.126 mm/ 2 yr E , respectively, with a mean of 0.050 mm/ 2 yr E .

Model Covariance Matrix
The model covariance matrix ( while  2 E measures the quality of the fit, and is defined as: For the velocity field inversion, we take  m E = 0.53 10 E to limit the slip rate maximum to the plate velocity. Using a fully unlocked prior model for the inversion, the L-curve in Figure S4 in Supporting Information S1 We therefore smooth the inverted slip acceleration so that it can reach 1.75 mm/ 2 yr E at maximum, hence  E L  E 1.75 mm/ 2 yr E . Given this constraint, the L-curve ( Figure S5 in Supporting Information S1, obtained using a zero-acceleration prior model) gives  m E = 0.23 10 E for the acceleration inversion. We finally notice that the  2 E values are large (  2 E = 307 and  2 E = 16.3 for the slip rate and the acceleration, respectively), which is due to the relatively small uncertainties we obtain. The latter are likely under-estimated since they only account for estimation, not model, errors. We however use them anyway as they allow to weight (in a relative sense) the inversion; absolute values of  E are therefore of little use here.

Slip Restitution
For each inversion, we evaluate the resolution by assessing the ability of each sub-fault to resolve a unit slip. We compute the surface displacements caused by a unit slip located on one sub-fault, and invert these displacements using the parameterization previously explained, to obtain a field of slip values. The associated resolution matrix from Tarantola and Valette (1982) is: The resolution matrix diagonal values range from 0 if the slip is not resolved, to one if it is fully resolved. Rather than the resolution itself, we look at the restitution index of the sub-faults corresponding to the sum of E R along each row, ranging from 0 to 1.5 ( Figure S6 in Supporting Information S1 for the slip rate and Figure S7 in Supporting Information S1 for the slip acceleration). We display with gridded meshes the subfaults with less than 90% of the slip restored.

Inter-Seismic Velocity Field and Locking
The spatial pattern of the horizontal inter-seismic surface displacements is given by the velocity terms E v of the linear trajectory model (Equation 1 when imposing E a = 0; see Section 3.2; this corresponds to choosing an averaged displacement rate over the 1997-2011 interval. A similar analysis is shown in Figure S8 in Supporting Information S1 when considering either the 1997 or the 2011 displacement rates). To compare our velocity field (Figure 4, black arrows) to the one of Mavrommatis et al. (2014), we transpose our velocity field from ITRF2014 to the NAM reference frame. We compute the velocity of the NAM plate at each station site using the absolute plate rotation poles of the NAM plate in the ITRF2014 given by Altamimi et al. (2017), and we subtract it from our velocity field. The subtraction gives us the velocity field in the NAM reference frame. Figure 4 shows the locking obtained when performing the inversion described in Section 4, with a fully locked prior model, jointly on both the PAC-NAM and the PHS-NAM interfaces. We also show in Figures S9 and S10 in Supporting Information S1 the locking obtained on PAC and PHS for smaller and larger  m E values, hence allowing for a smoother or rougher locking. Moreover, the locking obtained with a fully unlocked prior model is shown in Figure S1a in Supporting Information S1.  Table S2 in Supporting Information S1.

Acceleration Field
Our acceleration field (Figure 6, black and gray arrows) exhibits a landward acceleration North of Honshu with the amplitude decreasing when going West, and a trenchward acceleration South of 38.  5 E N. The first order of our acceleration pattern and amplitude are compatible with the ones observed by Mavrommatis et al. (2014) and Yokota and Koketsu (2015). Nevertheless, a notable difference is found in the Kanto region, where Yokota and Koketsu (2015) show a very limited acceleration compared to the rest of Honshu. Instead, we find an eastward acceleration with an amplitude of the same order as the north Honshu acceleration. According to our statistical test of the acceleration significance (see Section 3.3) most acceleration values in the Kanto region are significant ( Figure 6, black arrows). This acceleration suggests that some plate unlocking took place between 1997 and 2011. We conducted some refined analysis and sensitivity tests, as presented in Section 5.3, to track the origin of this difference but all our tests still show an acceleration in the Kanto region. Ultimately, we cannot reproduce exactly the same analysis as in Mavrommatis et al. (2014) and Yokota and Koketsu (2015) and therefore we cannot determine specifically which difference in our analysis could explain the difference in results.
The vertical accelerations ( Figure S11b in Supporting Information S1, outer circles) do not exhibit a spatially coherent pattern, likely because of the noise (that includes all non-tectonic phenomena) contaminating the acceleration signal.

Locking Change
The acceleration observed in the surface deformation implies a change in locking that can be simply evaluated. To do so, we estimate the locking both in 1997 and 2011 using the displacement rates modeled at all stations at the dates January 1, 1997 and February 6, 2011. Namely, we extract the modeled velocity at the first date, which is by definition the linear term of the quadratic model ( E v in Equation 1), and compute the associated locking ( Figure S8a in Supporting Information S1). Adding the acceleration term, we then compute the velocity at the second date as predicted by the trajectory model and invert for the associated locking ( Figure S8b in Supporting Information S1). We finally compute the difference in locking Δ E L ( Figure 5) and find that the maximum in absolute value is Δ E L = 0.3 as mentioned in Section 4.4.

Slip Acceleration
We perform a direct estimate of the accelerated slip on the subduction interface. First, we invert the slip acceleration field on the PAC-NAM subduction interface only. Figure 6 shows both accelerated (red) and decelerated (blue) slip rates. The slip acceleration corresponds to an increased forward slip on the subduction interface, or equivalently a decrease in locking. As explained in Section 4.4, a change in locking Δ E L = 0.3 over 14 years corresponds to a slip acceleration slip E a = 1.75 mm/yr 2 . Then, Figures 5 and 6 display similar relative variations, up to a simple proportionality factor imposed by the convergence rate and the duration of observation (Equation 13).
The residuals between the observed and modeled surface acceleration ( Figure S12 in Supporting Information S1) show little spatial coherence at the large scale, implying that the model does well in capturing the large scale pattern, even though the residuals are almost always significantly larger than the uncertainties as discussed in Section 4.4.
A remarkable result is the acceleration area South of 36.  5 E N, consistent with the observed surface acceleration mentioned earlier. This acceleration remains clearly visible even while varying the  m E parameter ( Figure S13 in Supporting Information S1). Because the presence of the PHS plate cannot be ignored for the stations in Kanto, we perform a second inversion using both the PAC and PHS plates ( Figure S14 in Supporting Information S1). Very limited changes are observed, with an increase of the slip acceleration on the PAC plate limited to 0.6 mm/ 2 yr E at maximum, and a deceleration of the order of 0.3 mm/ 2 yr E on the PHS plate. Since a trade-off exists between slip acceleration/deceleration on either plate (see Section 4.2), we cannot resolve with accuracy the proportion of change in slip rate hosted by each plate. Our analysis however clearly demonstrates that there exists a significant increase in slip rate between 1997 and 2011 South of 36.  5 E N, that is compatible with the observed acceleration of seismicity on the PAC-NAM interface as shown in Marsan et al. (2017).
We also evaluate the impact of adding the vertical displacement to the horizontal ones when inverting for the slip acceleration. Comparing Figure 6 and Figure S11 in Supporting Information S1, we observed that the vertical acceleration has very little effect on the inverted slip acceleration. Figure 6. Slip acceleration of the Pacific (PAC) plate with a zero acceleration prior model. The black and gray arrows depict the observed acceleration field: we distinguish the stations which have a signal-to-noise ratio greater than 3 (black), and those which fail to meet this criterion (gray, see Section 3.3). The green arrows represent the acceleration field predicted by the model. The amount of slip acceleration on the PAC slip interface is shown with the blue to red color scale. The gray grid shows the sub-faults for which the slip is poorly restored (see Section 4.5). Other elements are described in Figure 1.
As we present nine figures (within the main figures and supplementary figures in Supporting Information S1) based on PAC-NAM and/or PHS-NAM slip acceleration, we summarize all the model parameters associated with each figure in Table S3 in Supporting Information S1.

Refined Analyses and Sensitivity Tests
We here provide further analyses to compare our results with those of Mavrommatis et al. (2014), Tanaka and Yabe (2017), and Yokota and Koketsu (2015). We notably explore the effects on our results of the reprocessing of the GPS data (Section 5.3.1) and of the inclusion of slow deformation transients in the trajectory model (Section 5.3.2).

This Study's GAMIT/GLOBK Solution Versus F3 Solution
A significant difference with the data presented in Mavrommatis et al. (2014), Tanaka and Yabe (2017), and Yokota and Koketsu (2015) is the GPS solution used. Here, we have done a complete reprocessing of the GPS data, while previous studies used the GEONET F3 solution published by the Geospatial Information Authority of Japan (Nakagawa et al., 2009;Tsuji et al., 2017). Both solutions have been obtained with different processing strategies. Main differences include the use of different softwares (GAMIT/GLOBK vs. Bernese), and the projection in a different version of the International Terrestrial Reference Frame (ITRF2014 vs. ITRF2005) (Altamimi et al., 2007(Altamimi et al., , 2017. To compare our solution with the F3 solution, we use the Helmert parameters given by the ITRF website (http://itrf.ensg.ign.fr) to rotate the F3 solution into the ITRF2014 (more details in Text S1 and Table S1 in Supporting Information S1). Then, we apply the analysis described in Section 3 to the converted F3 solution, by again fitting the trajectory model of Equation 1, and then inverting the acceleration field on the PAC plate (Section 4). The obtained acceleration field, shown in Figure S15 in Supporting Information S1, displays similar first-order features than in Figure 6 and in Mavrommatis et al. (2014) and in Yokota and Koketsu (2015): North of Honshu, the acceleration field is pointing toward the East. The F3 solution implies higher southward acceleration amplitudes (  E 0.5 mm/yr) than our solution as shown by Figure S15a in Supporting Information S1.
The reason for this difference is not clear. It is likely not related to the difference in reference frame, since we have estimated that the transformation from ITRF2005 to ITRF2014 generates average changes in acceleration rate of the order of  E 0.0001 mm/ 2 yr E in Japan. Instead, the F3 solution is clearly more noisy than our solution, notably before 2000, which may lead to different results (see Figure S16 in Supporting Information S1). Indeed, the scatter in 1997 in the F3 solution is almost three times the scatter in 1997 in the GAMIT solution, and even in 2010, the scatter in the F3 solution is two times the one in the GAMIT solution.

Comparison With Previous Studies: Influence of Post-Seismic Modeling and Slow Slip Events
The other differences we find with respect to previous studies (Mavrommatis et al., 2014(Mavrommatis et al., , 2015Yokota & Koketsu, 2015) seem to be related to choices made during the time series analysis.
At the latitudes 36.
N, we find a weaker eastward acceleration than previous studies, both for the GAMIT and F3 data sets (Figure 6; Figures S15a and S16 in Supporting Information S1). This is probably explained by the consideration of different approaches for the post-seismic modeling in the time series analysis. Mavrommatis et al. (2014) perform a physical model for their post-seismic relaxation, while we compute a best fit that is independent from one station to the other. This will naturally lead to the fit of any signal visible in the position time series, and to the reduction of the westward acceleration found in the area, hence of the inverted accelerated slip.
Looking more specifically at Kanto, we observe a significant acceleration with both the F3 and GAMIT solutions (black arrow, Figure 6 and Figure S15 in Supporting Information S1) that were not reported by Mavrommatis et al. (2014) nor by Yokota and Koketsu (2015).
Contrary to Mavrommatis et al. (2014), Mavrommatis et al. (2015), and Yokota and Koketsu (2015), we here consider several transient deformation events (SSEs in Boso and offshore Honshu, volcanic unrests, see list in Table 2), and correct for those events in the time series modeling (see Section 3.2): we therefore compute the inter-SSE velocities and associated locking. To evaluate the impact of removing the SSEs, we additionally compute the inter-seismic velocities and the associated locking without correcting for the SSEs in the trajectory model (Figures S17a and S18a in Supporting Information S1). As expected, the changes between both models are mainly located around the Boso SSEs area (Figures S17b and S18b in Supporting Information S1). However, the changes extend to the North due to the Tohoku SSE. The changes are tiny and involve an area that is larger than the strict extent of the SSEs because of the inversion smoothing and because the SSEs can be seen in a wide range of stations. On the PAC plate interface, the inter-SSE model is more locked than the inter-seismic model by  E 0.02 under the Boso SSEs area ( Figure S17b in Supporting Information S1). For the PHS plate, the locking is on average larger in the inter-SSE model by  E 0.05 ( Figure  S18b in Supporting Information S1). These results can be explained by the fact that the SSEs contribute to the release of some of the accumulated slip.
More importantly, we find that not considering the Boso SSEs nor the Tohoku SSE during the computation of the trajectory model increases even further the slip acceleration on the PAC plate in the wider Kanto region ( Figure S19 in Supporting Information S1). (2017) Tanaka and Yabe (2017). The acceleration we infer in the Kanto region includes the slip area inferred by Tanaka and Yabe (2017) for both SSEs. We showed in this section that if SSEs were not modeled, a stronger acceleration would result from our analysis. Thus our results are compatible with the observations of Tanaka and Yabe (2017), in the sense that our long-term acceleration could equally be interpreted as long-duration SSEs with a total accommodated slip which rate would accelerate at the 14-year scale. Moreover, the background seismicity used by Tanaka and Yabe (2017) to support their SSEs observations show an acceleration coherent with our seismic study (Section 6).

Seismicity Acceleration and Slip Acceleration
The confirmation of acceleration slip allows us to relate it to the concomitant acceleration of seismicity. To do so, and as with the GPS time series, aftershock sequences must be removed from the data, as they would otherwise dominate the activity. This is traditionally done by using declustering methods. This approach was followed by Marsan et al. (2017), who studied how the rate of background earthquakes changed along the Japan Trench from 1990 to 2011. They made use of two distinct declustering methods, so to only keep acceleration patterns that remain robust with regard to the specific assumptions underlying each of these methods. As a result, locations with incoherent acceleration values obtained with these two methods were discarded from the analysis. This becomes a problem for our analysis as large parts of the PAC slip interface are discarded, and cannot then be used to construct a seismic prior, unless applying specific smoothing constraints. We therefore re-analyze the seismicity, by (a) keeping all the 1997-2011 (up to March) earthquakes within 20 km of the PAC slip interface; (b) using the same regular grid as with our slip inversions; (c) selecting at each grid node (center of the triangle patches) the earthquakes within a radius of 100 km (hence an earthquake can be used at several neighboring nodes) and extending this range if necessary until finding 500 earthquakes; and (d) model all the obtained time series with an earthquake rate: the space-varying seismicity acceleration parameter E a with our inverted slip rate changes. We measure the seismicity acceleration with: Figure 7a shows the spatial variations of the seismicity acceleration and Figure S20 in Supporting Information S1 shows the temporal evolution of the seismicity at 10 patches.
We convert this seismicity acceleration into the a priori expected slip acceleration over the same period, defining a new prior model as shown in Figure 7c. This conversion is based on the assumption that the slip rate ( dx dt slip / ) and the background seismicity rate (  E ) are proportional (the proportionality coefficient is assumed constant in time, but is allowed to vary spatially). This factor would likely be controlled by physical parameters, including the size distribution and density of asperities, the resulting seismic coupling at the local scale, and so on. In our modeling, the modeled slip rate is dx dt v a t t We here use 1 E t = 1997 so that 1 E t = R E t (same reference time than for Equation 1), and 2 E t = 2011.1. We thus get that the slip acceleration when comparing slip rate at those 1 E t and 2 E t is: where E L is the locking coefficient and PAC E v = 76 mm/yr. As we only examine changes over time, this approach allows to discard the unknown proportionality coefficient between the slip and seismicity rates.
The slip acceleration prior obtained using our analysis of the seismicity is shown in Figure 7c, and show good correlation with the already inverted slip acceleration of Figure 6, with the notable exception of the downdip 36.  5 E - 39 E N latitude range, for which seismicity acceleration is very clear and significant, while our slip inversion finds an overall deceleration. Seismicity in the co-seismic zone of the 2011 Tohoku earthquake is not strong during the 1997-2011 (up to March) period, as can be expected given the high degree of inter-seismic locking there, and the seismicity acceleration is very mild compared to the inverted slip rate change. Mavrommatis et al. (2015) also studied the acceleration in the seismicity but using repeating earthquake time series. Our seismic prior displays the same large-scale features, although the strong downdip acceleration we find is observed more updip in their analysis. Also, and contrary to them, we find a strong acceleration in the Kanto region and offshore Kanto. This particular feature is discussed in more detail in Section 7.
We show the inverted slip acceleration field using this seismic prior in Figure 7d. The locking used to compute the prior values in Figure 7b is the one obtained with a fully locked prior model; Figure S21 in Supporting Information S1 shows the results using the locking as given by a fully unlocked prior model. The accelerated slip obtained using the seismic prior shows similar features to the one shown in Figure 6. A noticeable change is the more focused acceleration in the Boso area, which now appears separated from the acceleration patch in the updip  36 It must be argued here that our inversions of the slip rate acceleration assumed a fully elastic medium. Visco-elasticity from the lower crust downwards could affect this inversion, especially as the time scale of our study (14 years) is comparable to typical Maxwell characteristic times in subduction modeling (Wang et al., 2012). A full investigation of visco-elastic models is however beyond the scope of the present study.
It has been argued by Hasegawa and Yoshida (2015) that the downdip acceleration (between  37 E and  38 E N) must have contributed to further stressing of the main asperity that eventually failed in 2011. Whether this acceleration is a common phenomenon in this region or cannot be answered with certainty. Hardebeck (2012) discusses the near complete release of the differential stress at the time of the 2011 shock. The released seismic moment only amounts to a 0.23-0.77 proportion of the moment that accumulated in the 1,142 years since the last equivalent megathrust earthquake in this region, the 869 CE Jogan earthquake that triggered a tsunami with comparable characteristics to those of the 2011 Tohoku event (Minoura et al., 2001(Minoura et al., , 2013Sugawara et al., 2012). Indeed, even accounting for large post-seismic relaxation lasting several decades, as seen after the 1960 Chile earthquake (Melnick et al., 2018;Wang et al., 2007), as well as incomplete plate coupling, a w E M 9 earthquake would need to occur every 300-700 years to match the average modeled rate of moment accumulation (Ozawa et al., 2012). The stressing rate (hence the slip rate on the interface) is thus likely to vary during the extended inter-seismic period. The short-term (i.e., over 14 years) average coupling could then be an overestimation of the longer time-averaged coupling. The fact that we infer a 14-year long decrease of coupling (even after removing co-seismic and post-seismic effects of the intermediate size earthquakes that have occurred in the region during this time period) is coherent with the fact that there must exist significant variations of coupling during the earthquake cycle in this region.
A novel feature of our analysis is the broad slip acceleration in the  35 E to  36 E N area, although weaker (by a factor 2) than for the  37 E - 38 E N region. This was not identified by previous studies, although this is a robust feature of our inversions (Figures 6 and 7). Even though some trade-off exists between changes in slip rate along the PAC-NAM and the PHS-NAM interfaces, our analysis suggests that the slip acceleration along the former dominates ( Figure 6 and Figure S14 in Supporting Information S1). Slip acceleration along the PHS-NAM interface is however independently suggested by the acceleration of seismicity documented in Reverso et al. (2016) offshore Boso, as well as the shortening of the Boso SSE recurrence interval that took place between 1996 and 2014, and which can be seen as the signature of an increasing loading in this area (Fukuda, 2018;Hirose et al., 2012;Ozawa, 2014). as suggested for the PAC-NAM interface by Nishikawa and Ide (2018) in the case of the 2008 w E M 7.0 Ibaraki earthquake. Intermittent slow slip therefore seems significant in this region, and accelerated seismicity could thus imply an acceleration of the controlling aseismic forcing rate. Given the mechanical coupling (i.e., contact) between the PAC and PHS plates offshore Boso, this interpretation is coherent with (a) the strong acceleration in background rate observed by Marsan et al. (2017) over the 1990-2011 period, (b) the remarkable acceleration in seismicity rate in the Boso area for the same period (Reverso et al., 2016), as well as (c) the shortening of Boso inter-SSE periods of times (Ozawa, 2014). All these observations point to an unlocking of the PAC-NAM interface there, estimated to be about 0.2-0.3 as shown in Figure 5. We find that the slip rate decelerates along the PHS-NAM interface ( Figure S14b in Supporting Information S1), in contradiction with the seismicity acceleration and SSE interval shortening. The smoothing performed by the inversion procedure implies a nearly uniform slip rate change over the modeled PHS plate, although the geometrical complexity of this plate (especially compared to PAC) could be expected to result in a more heterogeneous slip rate change; also, the inversion trade-off between the two PAC and PHS plates could hide more restricted areas of acceleration that would be compatible with the above observations of Ozawa (2014) and Reverso et al. (2016).
This inferred accelerated slip in the Kanto area is noticeably far from the 2011 rupture and remains therefore to be explained. It adds up to previous observations in other plate convergence contexts (Materna et al., 2019;Meltzner et al., 2015;Prawirodirdjo et al., 2010) that the slip rate, and thus the plate locking, could evolve significantly over years to decades outside the widely observed post-seismic slip following any mainshock.

Data Availability Statement
Raw time series data are accessible on the Observatoire des Sciences de l'Univers de Grenoble website associated to the DOI: https://doi.org/10.17178/GNSS.products.Japan.