Omori-like slow decay (p < 1) of postseismic displacement rates following the 2011 Tohoku megathrust earthquake

We infer the postseismic deformation mechanisms following the 2011 Tohoku megathrust earthquake via an analysis of onshore geodetic observations. We focus on the temporal decay characteristics of postseismic deformation using continuous time-series data at time scales that span many orders of magnitude by means of high-rate GNSS data. Our analysis indicates Omori-like power-law decay of the horizontal ground displacement rates, with p-value (0.69) that is significantly smaller than that of the aftershock occurrence (~ 1). This slow decay implies a (non-Maxwellian) viscoelastic relaxation mechanism other than afterslip since immediately after the mainshock, which is inferred using only onshore geodetic data. Spatial distribution of the Omori parameters implies that the postseismic deformation will continue over 100 years in a down-dip area of the northern part of the mainshock fault. The decay characteristics of vertical displacement rates are also almost Omori-like, but data deviation from the fitting line several 100 days after the mainshock might reflect the change of the dominant mechanism of the postseismic deformation. This multi-time scale geodetic approach will provide important constraints for future viscoelastic models of Earth’s interior.


Introduction
Postseismic deformation following large earthquakes has been quantitatively investigated using geodetic observations (Segall 2010). Rock rheology and friction control postseismic deformation via two primary mechanisms, aseismic transient afterslip along the mainshock fault (Heki et al. 1997) and viscoelastic relaxation of the asthenosphere (Nur and Mavko 1974), with both mechanisms often coexisting during deformation (Thatcher 1983;Wang et al. 2012). However, the ability to differentiate between these two mechanisms has been difficult due to limited observations and their similar spatial deformation patterns especially in horizontal components (e.g., Barker 1976), which complicates studies investigating the frictional and/or rheological characteristics of Earth's interior. This also presents a challenge for seismic hazard mitigation since afterslip can accelerate the rupture of adjacent earthquake source regions (Matsuzawa et al. 2004).
A recent continental earthquake study (Ingleby and Wright 2017) demonstrated that afterslip (or localized power-law creep in a shear zone) is the dominant postseismic deformation mechanism. They focused on nearfault geodetic observations at time scales that spanned over five orders of magnitude, and found that the postseismic velocities were inversely proportional to time since the earthquake (following the Omori law for aftershocks (Omori 1894), with a p-value of ~ 1), with logarithmic increase in postseismic displacement over time. This deformation characteristic was consistent with a rate-and-state-dependent friction law (Marone et al. 1991) and a hypothesis that aftershocks are triggered by afterslip (Perfettini and Avouac 2004). Linear relations between aftershock number and afterslip following the 2005 Nias earthquake and the 2010 Maule earthquake have been reported in the subduction megathrusts (Hsu et al. 2006;Lange et al. 2014). For the 2011 Tohoku earthquake, Ozawa et al. (2012) suggested a good correlation between the aftershock number and the moment of the afterslip based on daily data during 7 months after the mainshock (but they did not propose a linear relation of them).
For subduction zone interplate earthquakes, viscoelastic relaxation as well as afterslip has been identified as a notable mechanism (Wang et al. 2012). Seafloor geodetic observations of the 2011 Tohoku megathrust earthquake (moment magnitude (Mw) 9.0) provided clear evidence for viscoelastic relaxation as the dominant deformation mechanism (Sun et al. 2014;Watanabe et al. 2014). However, diverse postseismic deformation models that incorporate afterslip and viscoelastic relaxation have been proposed (Diao et al. 2014;Yamagiwa et al. 2015;Hu et al. 2016;Muto et al. 2016;Iinuma et al. 2016;Freed et al. 2017;Sobolev and Muldashev 2017;Suito 2017;Agata et al. 2019), with the model results being strongly dependent on their structures, constitutive laws, and viscoelasticity parameters. Some studies assumed layered structures of parameters (Diao et al. 2014;Yamagiwa et al. 2015), some reduced to two-dimensional problems (Muto et al. 2016;Sobolev and Muldashev 2017), but the others considered three-dimensional parameter distribution, with various constitutive laws such as the Maxwell viscosity, the bi-viscous Burgers rheology, or power-law rheology. We suggest that an effective method to constrain the complex deformation processes following a megathrust earthquake is via the analysis of continuous (onshore) time-series data at time scales that span many orders of magnitude, following the aforementioned continental earthquake study (Ingleby and Wright 2017).
Here we clarify the temporal evolution of postseismic deformation following the 2011 Tohoku Earthquake. We accomplish this analysis using continuous (onshore) time-series data at time scales that span many orders of magnitude immediately after the mainshock. We compare the ground displacements with aftershock occurrence to infer the dominant deformation mechanism.

GNSS data
We used GNSS (Global Navigation Satellite System) data from the GEONET array, which are provided by GSI (Geospatial Information Authority of Japan), to study the postseismic deformation of the Mw 9.0 Tohoku megathrust earthquake on March 11, 2011, at 14:46 JST (Japan Standard Time).
We obtained ground displacements via analysis of the 30-s RINEX files ("high-rate GNSS") obtained between 9:00 JST on March 10, 2011, and 9:00 JST on March 18, 2011, using the GSILIB software (GSI 2015). We obtained the kinematic time-series displacements at the GEONET stations relative to a fixed station, Otaru 1 (0013; 141.031° E, 43.178° N). 300 GEONET stations in northeastern Japan are used in this analysis. The GNSS settings are as follows: (i) ephemeris from the IGS final orbit, with a 15-min interval; (ii) elevation mask of 15°; (iii) solid Earth tide correction based on the IERS model (McCarthy and Petit 2003); (iv) ionospheric delay correction via estimation of the slant total electron content; (v) tropospheric delay correction via estimation of the horizontal gradient of the zenith total delay; (vi) phase center variation corrections for both the satellites and receivers; and (vii) P1-C1 differential code bias correction.
We also used the daily ground displacement solutions provided by GSI from March 18, 2011, to December 31, 2018(Nakagawa et al. 2009). Therefore, our analysis employs two types of GNSS time-series data at time scales that span many orders of magnitude. Figure 1 shows an example of the time-series data.

Postseismic deformation following the Tohoku Earthquake
We calculated the ground displacement rates following the mainshock at 14:46 JST on March 11, 2011, for both the N-S and E-W (horizontal) components of the GNSS time-series data. We calculated the displacement rates also for the vertical component, but they contained much noise due to the satellite geometry problem (e.g., Berber et al. 2012). Therefore, in the followings, we mainly use N-S and E-W (horizontal) components for analyses. In "Discussion" section, we will provide analyzed results for the vertical component. The time frames after the mainshock for the calculation are as follows: 0-6 h (hour), 6-21 h, 21-33 h, 33-81 h, 81-153 h, 7-21 d (day), 21 d to 2 mo (month), 2-6 mo, 6-11 mo, 11-23 mo, 23-45 mo, 45-69 mo, and 69-93 mo.
We had to extract the deformation signal from the 30-s data since they contain high noise levels. We first removed the coseismic offsets of two large aftershocks that occurred near the coast on March 11, 2011 ), from the time-series data. One was a Mw 7.4 earthquake at 15:08 JST, which occurred along the northern extension of the mainshock fault, and the other was a Mw 7.7 earthquake at 15:15 JST, which occurred along the southern extension of the mainshock fault. Next, using a moving-average method (Cleveland et al. 1990), we removed a periodic diurnal component and an irregular component from the time-series data. We then smoothed the cleaned time-series data using Friedman's super smoother (Friedman 1984) to Page 3 of 10 Morikami and Mitsui Earth, Planets and Space (2020) 72:37 obtain regression lines of each time-series, with a start time of 14:51 JST on March 11, 2011. We then calculated the ground displacement rates. Lastly, we excluded the noisy observation points for the displacement rates, where we define the outliers as those > 1.5 times the interquartile range from the third quartile value during each observation period. Friedman's super smoother was also applied to the daily ground displacement data every 0.1 years after the mainshock, for the calculation of the displacement rates. Finally, we removed interseismic (pre-Tohoku) displacement rates from the above estimation in order to obtain true signals of the postseismic deformation. We identified a representative period from December 1996 to December 2002, because transient deformation signals after 2003 in northeastern Japan were reported (Suito et al. 2011;Ozawa et al. 2012;Heki and Mitsui 2013). Although postseismic deformation of the 1994 Sanriku-oki earthquake (Mw 7.7) might have continued in the selected representative period (Ozawa et al. 2007), we determined it as a better choice than the following period after 2003. Then we obtained 6-year deformation amount as difference between medians of the first and last months, and calculated the interseismic displacement rates. We can use the time-series data at 141 GEONET stations ( Fig. 1) in this study.
We expect the ground displacement rates to decrease with time since the mainshock. The temporal evolution of the horizontal ground displacement rates for every GNSS station is shown in Fig. 2. The decay in ground displacement rate resembles the modified Omori law for aftershocks (Utsu 1957): where d(t) is the ground displacement rates for time t since the mainshock, and a d , c d , and p d are constants. We estimate the parameters via a least-squares fit (all values are assumed to be non-negative). (1)

Seismicity decay after the Tohoku Earthquake
We extracted the aftershocks from the JMA earthquake catalog that occurred between March 11, 2011, and December 31, 2017 (this date is different from that of the postseismic deformation), with focal depths shallower than 100 km and epicenters located within the map in Fig. 1, for a comparison with the postseismic ground displacement rates. The decay in aftershock occurrence follows the modified Omori law as: where s(t) is the seismicity rate for time t , and a s , c s , and p s are constants. We estimate the parameters via a leastsquares fit (all values are assumed to be non-negative). We determined a lower magnitude limit of 4.0 based on the size-frequency distribution of the above data and the Gutenberg-Richter law (Gutenberg and Richter 1944). In order to discuss the aftershock decay without effects of inland earthquakes, we limit the analysis area near the source fault of the mainshock as shown by Fig. 3.

Postseismic deformation following the Tohoku Earthquake
For the "Omori-like" law of the postseismic displacement rates (Fig. 2), we estimated the p-value ( p d ) as 0.69(±0.01) via a least-squares fit (all values are (2) s(t) = a s (t + c s ) p s , assumed to be non-negative) in logarithmic space, with a d = 2.1(±0.1) × 10 −7 m s −1 and c d = 0.28(±0.07) d, where ± mean 95% confidence intervals. A reference of convergence time of postseismic deformation, t c , where we define t c as time of d(t) = 1 cm/year (a typical value of tectonic deformation in the research area, referring to Ozawa et al. (2012)), is about 30 years. If we change 1 cm/ year to 1 mm/year, t c is about 860 years. We found that the p-value ( p d ) is significantly smaller than that for continental earthquakes (Ingleby and Wright 2017), which indicates the slow decay of postseismic deformation following the 2011 Tohoku Earthquake. The integration of d over time indicates that the displacement since the mainshock is not logarithmic, which is in contrast to a p-value of ~ 1. Because a p-value of ~ 1 (logarithmic increase in displacement) was related to frictional afterslip (Marone et al. 1991), our result implies additional mechanisms of the postseismic deformation. Figure 4 shows examples of the Omori-like fitting for the horizontal displacement rates at each GNSS station. Additional file 1: Table S1 compiles all of the fitted parameters. We calculate the spatial distribution of the reference of convergence time of postseismic deformation, t c , as time of d(t) = 1 cm/year (a typical value of tectonic deformation). Figure 5 illustrates the t c distribution (see Additional file 1: Table S1 for details). t c varied by two orders of magnitude in the studied region. We found that the postseismic deformation will continue over 100 years especially in a down-dip area of the northern part of the mainshock fault. By contrast, the postseismic deformation to the southwest of the mainshock fault will soon end.

Seismicity decay after the Tohoku Earthquake
We investigated the modified Omori law of the aftershocks for comparison with the postseismic deformation (Fig. 6). We estimated the p-value ( p s ) as 0.95(±0.02) via a least-squares fit (all values are assumed to be non-negative) in logarithmic space, with a s = 520(±59) d −1 and c s = 0.32(±0.07) d.
We found that the p-value of the aftershocks ( p s = 0.95) is ~ 1, which is considerably larger than the p-value of the ground displacement ( p d = 0.69). This difference ( p d < p s ) is significant because it indicates that the postseismic ground deformation decayed more slowly than the aftershock occurrence. For example, the seismicity rates several years after the mainshock are almost one order of magnitude lower than those calculated using the p-value derived from the ground displacement rates, as shown in Fig. 6.
Our analysis clearly highlights that the ground displacement rates following the 2011 Tohoku megathrust earthquake were not inversely proportional to time, as observed for continental earthquakes (Ingleby and Wright 2017), with a different aftershock decay compared to that of the Nias megathrust earthquake (Hsu et al. 2006). An exhaustive literature review indicates that such deformation characteristics ( p d < p s ) have not been reported to-date.
One might consider that p d does not necessarily need to equal p s in terms of rock friction, particularly since sophisticated earthquake models that employ a rate-and-state-dependent friction (Dieterich 1994) to trigger earthquakes via afterslip (Helmstetter and Shaw 2009) have predicted cases with different p d and Fig. 4 Evolution of the horizontal ground displacement rates after the 2011 Tohoku Earthquake at representative GNSS stations (Fig. 1). The red dashed lines represent the best-fit based on Eq. (1) Page 6 of 10 Morikami and Mitsui Earth, Planets and Space (2020) 72:37 p s values. However, these studies suggest that p d ≥ p s , making p d < p s impossible for p d = 1. The model prediction ( p d ≥ p s ) is consistent with the well-known concept that "each earthquake generates additional earthquakes" (e.g., Kagan and Knopoff 1981;Ogata 1988). Therefore, our result of p d < p s indicates additional postseismic deformation mechanisms that do not effectively trigger aftershocks.

Relation between the Omori-like decay and deformation mechanisms
Maxwellian viscoelastic models have been used in many studies to explain postseismic deformation at subduction zones. For instance, Suito (2017) developed a three-dimensional heterogeneous model incorporating the subducting Pacific Plate and the mantle wedge. Figure 7 shows comparison of Model 4 in Suito (2017) (final model with a weak mantle wedge and thin layer beneath the subducting slab, of which viscosities are on the order of 10 18 Pa s) with our results of the horizontal displacement rates at the same stations. We confirm that the Maxwellian viscoelastic relaxation can well explain our results for a later span (since one hundred d after the mainshock). However, as investigated by Ingleby and Wright (2017), the Maxwellian viscoelastic relaxation does not behave as the Omori-like law and cannot explain the early span of our results. Some mechanisms play key roles. Afterslip has been frequently raised as the most dominant mechanism (e.g., Diao et al. 2014).
Our results of a slower postseismic deformation decay ( p d = 0.69) than the aftershock occurrence indicate a non-negligible contribution of additional mechanisms other than afterslip (with a p-value of ~ 1) immediately after the mainshock. This implies that the interpretation of early postseismic deformation as afterslip in the literature (e.g., Mitsui and Heki 2013) is inappropriate, at least in subduction zone environments. Non-Maxwellian (bi-viscous Burgers rheology or power-law rheology) (2) with p s = 0.95 , whereas the thin red dashed line shows a (failed) prediction with the same p-value as the case for the horizontal ground displacement rates ( p d ) Fig. 7 Comparison between the evolution of the horizontal ground displacement rates (Fig. 2) and a Maxwellian viscoelastic simulation (Suito 2017) after the 2011 Tohoku Earthquake. The green points represent the simulation results of Model 4 (final model) in Suito (2017) Page 7 of 10 Morikami and Mitsui Earth, Planets and Space (2020) 72:37 viscoelastic relaxation of the asthenosphere is one plausible mechanism. Another possible mechanism of poroelastic rebound with fluid flow may not cause large deformation far from the source fault, compared with viscoelastic relaxation (Hu et al. 2014).
Although previous studies have highlighted the importance of viscoelastic relaxation, we note that the estimated viscoelastic parameters are heavily dependent on the assumed parameter structures and constitutive laws. For example, some studies of the 2011 Tohoku Earthquake have estimated the viscosity of the asthenosphere at ~ 10 19 Pa s (Diao et al. 2014;Yamagiwa et al. 2015;Hu et al. 2016), whereas other studies have estimated it at ~ 10 18 Pa s (Freed et al. 2017;Suito 2017). The existence of a subducting slab and localized weaker zones have been suggested (Hu et al. 2014(Hu et al. , 2016Muto et al. 2016;Freed et al. 2017;Suito 2017). Furthermore, a stressdependent and non-constant viscosity due to power-law dislocation creep (Freed and Bürgmann 2004) has been implemented in recent studies (Sobolev and Muldashev 2017;Agata et al. 2019), where a sophisticated study (Sobolev and Muldashev 2017) simulated ground displacement rates by varying the model parameters in a two-dimensional body. This behavior is intrinsic to giant earthquakes which mostly can occur only in subduction zones, because a giant earthquake has a huge effect on the viscosity due to the power-law of dislocation creep.
Decay characteristics of postseismic deformation must depend on viscoelastic parameters, structures, and constitutive laws, according to the previous numerical experiments of viscoelasticity (e.g., Figures 8-9 in Sobolev and Muldashev 2017;Figure 3 in Ingleby and Wright 2017). Although multiple mechanisms may work during the postseismic deformation as many studies exhibited, the evolution of the horizontal displacement rates seemed to obey a single ). This point is interesting for a future issue of the mechanisms of postseismic deformation (for example, do we consider "afterslip" as a separate mechanism from viscoelastic relaxation?). Analysis of the p-value distribution using continuous time-series data at time scales that span many orders of magnitude may, therefore, be an important constraint in guiding future viscoelastic models of Earth's interior.

Omori-like decay of vertical displacement rates
In order to further discuss the mechanisms of the postseismic deformation, we examined the postseismic decay of the displacement rates in the vertical component although the vertical component of the GNSS time-series contained much noise (e.g., Berber et al. 2012) especially for the 30-s sampling high-rate data. On the basis of the same methods as the horizontal components, we calculate the vertical displacement rates after the 2011 Tohoku earthquake. Figure 8a shows the calculation results. The evolution of the vertical displacement rates can be roughly fitted by the Omori-like law using the same p d and c d as those for the horizontal components (with a d = 2.7(±0.2) × 10 −8 m s −1 ), however, the fitting seems worse for the initial (until around 1e+01 d) and last (since several hundred d) spans in Fig. 8a. The initial span consists of the 30-s sampling high-rate data, of which much noise may lead to artificially large displacement rates. By contrast, during the last span, the deviation of the data from the fitting line possibly reflects change of the dominant mechanism of the postseismic deformation in mid-2013, as implied by a recurrent neural network method for the horizontal components (Yamaga and Mitsui 2019). An enlarged figure (Fig. 8b) with median values of the displacement rates for the time frames emphasizes this fact. Moreover, the data during the last span also slightly deviate from the Maxwellian viscoelastic relaxation simulation (Suito 2017), which implies another mechanism after mid-2013.

Spatial distribution of estimated parameters
The spatial p-value distribution provides a hint for discussing postseismic deformation. Figure 9 shows the estimated p-value ( p d ) distribution (see Additional file 1: Table S1 for details). The p-values to the southwest of the mainshock fault are generally larger than the value for every stations (0.69), approaching 1. This implies that the additional deformation mechanism for the early span, which is likely (non-Maxwellian) viscoelastic relaxation of the asthenosphere, is disturbed and afterslip (with logarithmic increase in displacement) dominates the postseismic deformation in this southwestern region. The existence of another subducting slab (Philippine Sea Plate) beneath this region may be related to this mutation, as one numerical study showed that the depth extent of the subducting slab strongly influences viscoelastic relaxation (Freed et al. 2017). Another possible origin of the larger p-values in the southwestern region is difference in source depths of the postseismic deformation, which may lead to different rheological properties, as Iinuma et al. (2016) estimated that afterslip (from April 23, 2011, to December 10, 2011) occurred at the shallower plate interface in the southern part than the northern part. By contrast, the p-values around the Echigo Plain (due west of the centroid moment tensor of the 2011 Tohoku Earthquake and along the western coast of Japan) tend to be smaller than the value for every stations (0.69), although the Echigo Plain is far from the source fault of the 2011 Tohoku Earthquake. It might be related to a possible reverse fault with stable sliding beneath the Echigo Plain ; Meneses-Gutierrez and Sagiya 2016).

Conclusions
We observed Omori-like decay of the horizontal ground displacement rates following the 2011 Tohoku Earthquake, with a p-value (0.69) that is significantly smaller than that of the observed aftershock occurrence, which is ~ 1. These observations may reflect a contribution of a (non-Maxwellian) viscoelastic relaxation mechanism since immediately after the mainshock, because (rateand-state) frictional afterslip should be accompanied by a p-value of ~ 1. The spatial p-value distribution at each observation station also supported this conclusion. The decay characteristics of the vertical ground displacement rates were also almost Omori-like, but the data seemed deviated from the fitting line several 100 days after the mainshock which might reflect the change of the dominant mechanism of the postseismic deformation. Analysis of the p-value distribution using continuous time-series data at time scales that span many orders of magnitude would be an important constraint in guiding future viscoelastic models of Earth's interior.
Additional file 1: Table S1. Results of the Omori-like fitting for the horizontal displacement rates at each GNSS station.

Abbreviations
GNSS: Global Navigation Satellite System; GSI: Geospatial Information Authority of Japan; JMA: Japan Meteorological Agency; JST: Japan Standard Time.