Various slips at plate interfaces in subduction zones have been documented, and the role of aseismic slips in strain release and stress re-distribution has been revealed. Prior to the M9 2011 Tohoku earthquake in Japan, various small-scale aseismic slip events which possibly affected the occurrence of the 2011 Tohoku earthquake had been reported. However, the small-scale aseismic slip behaviors after the M9 2011 Tohoku earthquake have not been well-clarified because the dominant postseismic deformation of the 2011 Tohoku earthquake has hidden minor aseismic slips. In this study, we performed trajectory models for the GNSS time-series after eliminating the postseismic deformation of the 2011 Tohoku earthquake, and we extracted afterslips followed by M7 class earthquakes, and an ongoing long-term aseismic slip event (L-ASE) since late 2019. The spatial extent of these aseismic slips was also clarified, and the 2019 L-ASE was estimated to be in the down-dip area of the seismogenic zone. As a similar L-ASE was also found before the 2011 Tohoku earthquake, a comparison of these L-ASEs suggests that the repeated L-ASEs promoted strain accumulation in the shallow plate interface, including the M9 rupture area, and may be markers preceding M7 class seismicity. Because L-ASEs are key phenomena in assessing sequential seismic behaviors in the Tohoku subduction zone, monitoring the ongoing L-ASE and investigating its physical interaction with other slip behaviors is important.

Geodetic fault slip inversions have been generally performed by employing a least squares method with a spatial smoothing constraint. However, this conventional method has various problems: difficulty in strictly estimating non-negative solutions, assumption that unknowns follow the Gaussian distributions, unsuitability for expressing spatially non-uniform slip distributions, and high calculation cost for optimizing many hyper-parameters. Here, we have developed a trans-dimensional geodetic slip inversion method using the reversible-jump Markov chain Monte Carlo (rj-MCMC) technique to overcome the problems. Because sub-fault locations were parameterized by the Voronoi partition and were optimized in our approach, we can estimate a slip distribution without the spatial smoothing constraint. Moreover, we introduced scaling factors for observational errors. We applied the method to the synthetic data and the actual geodetic observational data associated with the 2011 Tohoku-oki earthquake and found that the method successfully reproduced the target slip distributions including a spatially non-uniform slip distribution. The method provided posterior probability distributions with the unknowns, which can express a non-Gaussian distribution such as large slip with low probability. The estimated scaling factors properly adjusted the initial observational errors and provided a reasonable slip distribution. Additionally, we found that checkerboard resolution tests were useful to consider sensitivity of the observational data for performing the rj-MCMC method. It is concluded that the developed method is a powerful technique to solve the problems of the conventional inversion method and to flexibly express fault-slip distributions considering the complicated uncertainties.