Function Model Based on Nonlinear Transient Rheology of Rocks: An Analysis of Decadal GNSS Time Series After the 2011 Tohoku‐oki Earthquake

The postseismic Global Navigation Satellite System (GNSS) time series after the 2011 Tohoku‐oki earthquake is often empirically explained by one or more logarithmic and exponential decay functions. However, a function model with these decay functions struggles to illuminate various deformation mechanisms. Here, we propose a new function model incorporating laboratory‐derived constitutive laws (power‐law flow and velocity‐strengthening friction) utilized in mechanical postseismic models of the last decade. We demonstrated that our function model accurately predicts the decade‐long time series of inland GNSS stations, including coastal areas where significant postseismic uplift continues now. Moreover, it decomposes them into displacements due to viscoelastic relaxation and afterslip, similar to results produced by previous stress‐dependent postseismic models of the 2011 Tohoku‐oki earthquake. This physics‐based function model is an effective tool to ease the process of forecasting long‐term GNSS time series by dividing them into two dominant mechanisms at the plate interface and the surrounding viscoelastic mantle.

2 of 13 ple, a logarithmic function empirically explains the decay of AS (Freed et al., 2007;C. Marone, 1998;Marone et al., 1991;Scholz, 1998) as well as the biviscous deformation (Hetland & Hager, 2006), whereas exponential functions represent Maxwellian viscous deformation (Hetland & Hager, 2006). Using an empirical fit, previous studies (Fujiwara et al., 2022;Tobita, 2016) proposed several combinations of logarithmic and exponential functions to explain the postseismic GNSS time series after the 2011 Tohoku-oki earthquake. However, these empirical function models often lack evidence and persuasiveness of the physical mechanisms.
Rapid deformations after a large earthquake reveal the coexistence of two phases: transient (short-term) and steady-state (long-term) deformation. Power-law rheology (Freed & Bürgmann, 2004) and transient rheology (i.e., Burgers rheology, Pollitz, 2003) were proposed to explain the biphasic deformation. Several postseismic models combined transient creep and power-law flow to explain the postseismic geodetic observations of the 1999 Hector Mine earthquake (Freed et al., 2010(Freed et al., , 2012, 2004 Sumatra earthquake (Masuti et al., 2016), and 2011 Tohoku-oki earthquake (Agata et al., 2019;Dhar et al., 2022Muto et al., 2019). Therefore, similar to Tobita (2016) and Fujiwara et al. (2022), the search for an appropriate function model analytically equivalent to the nonlinear stress-dependent rheology may lead to a further understanding of the physical mechanisms of postseismic deformation.
Here, we propose a physics-based analytical function for fitting and predicting the decadal GNSS time series. We built our decay functions using the laboratory-derived constitutive laws associated with the postseismic deformation of the 2011 Tohoku-oki earthquake. Utilizing our function model, the decomposition of GNSS displacements into the displacements caused by VR and AS can be achieved. The analytical functions can also facilitate the successful inversion of postseismic deformation to infer various source mechanisms in the Japan trench and island arc of NE Japan.

GNSS Observations
We used the daily coordinates of F5 solutions (Muramatsu et al., 2021) provided by the GNSS Earth Observation Network System (GEONET), a dense and widespread geodetic network operated by the Geospatial Information Authority of Japan (GSI). Among the observations at ∼1,300 permanent GNSS stations referenced in ITRF2014, we selected relatively high-quality data from 110 GNSS stations for this study (Figure 1). We analyzed the GNSS time series from 12 March 2011-11 March 2021 after removing the coseismic offsets of large earthquakes (Table  S1 in Supporting Information S1) to maintain the continuous trend of the time series. The common changes and systematic noises were removed by referencing the original time series with the observations at Fukue station (station ID: 950462 in Figure 1) (e.g., Fujiwara et al., 2022;Tobita, 2016). We utilized the east-west, northsouth, and up-down displacements of six GNSS stations: 940028 (Miyako), 960549 (Yamoto), 940042 (Hitachi), 950193 (Minase), 950151 (Kanita), and 950233 (Ryoutsu2) (blue squares in Figure 1) as a training data set for our main modeling. These stations were selected because they display high-quality data, wider variation in displacement trends, and broader coverage across NE Japan (Tobita, 2016). The data set of these GNSS stations (such as Miyako, Minase, Yamoto) helps obtain a stable solution for our function model, similar to Tobita (2016). We deployed the remaining 104 GNSS stations (red circles in Figure 1) to evaluate the performance of the proposed function model.

Formulation of Function Model
We built the function model for postseismic deformation given that surface deformations occur by viscoelastic relaxation and afterslip. Our function model includes the time-dependent behavior of these two mechanisms, without involving any structures of the subduction zone. This modeling strategy is motivated by several powerlaw models that qualitatively explains the observed postseismic deformation after the Tohoku-oki earthquake, regardless of differences in their modeling structures (Dhar et al., 2023 and references therein). The effectiveness of power-law flow was also advocated for other earthquakes (e.g., Freed & Bürgmann, 2004). The displacement D at each GNSS station in any direction can be modeled as a linear combination of responses due to viscoelastic strain ( V ) in the mantle and continuous afterslip ( ) at the plate interface, where is the time after the mainshock. The VR , AS , and 3 are scalar constants with different values for each station and each directional component. The steady velocity was assigned to the average linear velocity measured from 1997 to 2000 for each GNSS station (Fujiwara et al., 2022;Tobita, 2016). We adopted the semi-analytical expression of viscoelastic strain due to power-law creep and fault slip due to rate-strengthening friction from previous studies (Barbot & Fialko, 2010;Barbot et al., 2009;Tang et al., 2019), as follows: where VR and AS represent coseismic stress perturbation on mantle and plate interface respectively; 1 , , and 0 indicate relaxation timescale, shear modulus of the mantle, aspect ratio of mantle geometry, stress-exponent, the effective stiffness of slip area, direct-effect parameter, normal stress, and the reference velocity for the AS, respectively. In postseismic deformation, VR in the mantle shows a short-term transient creep accompanied by long-term steady-state deformation (e.g., Wang et al., 2012). This biphasic deformation is often represented by Burgers rheology model where total strain is the sum of strains in two spring-dashpot systems (Maxwell and Kelvin elements) connected in series. Therefore, Equation 1 becomes, where 1V and 2V indicate the strains due to the steady-state and transient creep, respectively. By combining Equations 2-4, we obtained the following: where 1 and 2 represent the relaxation timescales for steady-state and transient creep of the viscoelastic mantle, respectively. The 2 indicates strain-hardening coefficient. The above equation can be rewritten as, In the above equations, the is a dimensionless parameter related to the frictional parameter of the plate interface. The scalar constants 1 and 2 are time-invariant. The is set to 3 (Agata et al., 2019;Muto et al., 2019;Qiu et al., 2018;Tang et al., 2020). Because the value of is poorly known in laboratory experiments, we adopted = 1 (i.e., 1 = 2 ) following the previous postseismic models (Agata et al., 2019;Barbot, 2020;Dhar et al., 2022;Fukuda & Johnson, 2021;Hu et al., 2014;Hu, Bürgmann, Uchida, et al., 2016;Muto et al., 2019;Shi et al., 2020). We adopted 0 = 0.08 m/yr same as the stable subduction rate (Apel et al., 2006;Sella et al., 2002). Previous studies indicated that ∼90 km downdip width of the plate interface hosts significant AS in the postseismic period (Dhar et al., 2022;Iinuma et al., 2016;Muto et al., 2019). Thus, we set the linear dimension of slip area ≈ 90 km and the shear modulus = 30 GPa ( * ≈ ∕ ; Barbot et al., 2009) to yield 2 0 * = 0.0533 MPa/yr. Because the average stress perturbation by the Tohoku-oki rupture is in the order of ∼1 MPa (Iinuma et al., 2012), we set AS at 1 MPa, which also allows us to directly compare to the estimates of in previous postseismic models. Therefore, Equation 6 becomes, We used Equation 8 as a regression model with three coefficients ( 1 , 2 , and 3) and two decay functions ( VR for VR; AS for AS) in the form of, and

Training of Function Model
We assigned the first 3.9 years of the GNSS time series (12 March  the regression coefficients are calculated using the derivative-free simplex search method (Lagarias et al., 1998;Montgomery et al., 2015). We evaluated the prediction performance of our function model by calculating the residuals (cumulative observed minus modeled displacements) for 10 years after the mainshock. Cumulative displacements in the observed time series were obtained using a 14-day moving average which reduces the adverse effect of high time-series noises (particularly in the vertical components) on the stability of model solution. We deployed a Bayesian search involving the parameters 1 , 2 , and in the range of 0.001-10 years, 0.001-10 years, and 0.1-2, respectively, which covers the values estimated by the previous function models (Fujiwara et al., 2022;Tobita, 2016) and numerous postseismic models (Fukuda & Johnson, 2021 and references therein). To optimize the model, we used the data set of six training stations (see Section 2.1), the likelihood function following Fukuda and Johnson (2021) and the Gaussian Process Upper Confidence Bound (GP-UCB) for sampling the parameter space (Cox & John, 1992;Srinivas et al., 2012). The details of Bayesian search and GP-UCB methods are provided in Text S1 in Supporting Information S1.

Performance of Function Model
Figure 2 shows our function model's fitting and prediction performance at the six GNSS stations assigned for training in the main modeling. The fitting and prediction performances were reasonably good (an average residual of ∼1 cm for 10 years) for both the horizontal and vertical time series. The residuals may increase slightly with time (∼0.1 cm/year) due to the increasing model uncertainties (see residual time series in Figure S1 in Supporting Information S1). The decomposition of the calculated displacements into the VR (= 1 VR) and AS (= 2 AS) components demonstrates that the VR dominates the horizontal time series (left panel, Figure 2). For vertical deformation, the VR contribution varies from station to station (right panel, Figure 2). For several stations along the Pacific coast (e.g., Yamoto 960549 in Figure 2b), the VR initiates the postseismic uplift (manifestation of the cold nose; Luo & Wang, 2021), but a large AS supersedes in the first year onwards. The temporal evolution of VR and AS components is broadly consistent with that of stress-dependent models (Dhar et al., 2022;Fukuda & Johnson, 2021;Muto et al., 2019). Although the minimum 2-year fitting period is required to avoid the influence of seasonal variation on the model parameters (see Text S2 and Figures S2-S4 in Supporting Information S1 for details), we chose a 3.9-year fitting period to keep the residuals minimum as suggested by Fujiwara et al. (2022).
In fact, all fitting periods longer than 2 years provide almost same fitting parameters, indicating the robustness of our analysis ( Figure S3 in Supporting Information S1).
We further evaluated the 10-year prediction by our function model at 104 GNSS stations (red dots in Figure 1). Figures 3a and 3b illustrate the residual displacements of all stations (averages of ∼0.62 cm, ∼1.40 cm, and ∼0.95 cm for north-south, east-west, and up-down components, respectively) estimated by our best-fit model (see also Table S2 in Supporting Information S1). These residuals are well explained by the uncertainties in the prediction model, as illustrated by our time-series forecast in five representative stations, assigned for evolutions (red line + orange shade in Figures 3c-3g). Similarly, the consistency between the observed and modeled time series of all stations is demonstrated in Figures S5-S20 in Supporting Information S1.

Relaxation Time Scales and Transient Rheology
The best-fit values of 1, 2, and were estimated to be ∼0.9754 years, ∼0.0428 years, ∼0.9013, respectively (see Figure S21 and Table S3 in Supporting Information S1 for details). The best-fit parameters did not change significantly for different sets of training stations, indicating our model's robustness ( Figure S22 in Supporting Information S1). When varies, the estimates of 1 and 2 are affected, although the changes in the , modeled time series, and modeled VR and AS displacements remain subtle ( Figure S23 in Supporting Information S1). Similarly, the choice of AS only affects the estimates of the nonlinear parameters, although it imposes insignificant effect on the modeled displacements ( Figure S24 in Supporting Information S1). Our estimated pair of relaxation timescales ( 1 and 2 ) is at the lower end of the value range (0.0038-10 years), reported previously by logarithmic-exponential functions (Freed et al., 2017;Fujiwara et al., 2022;Hu et al., 2014;Tobita, 2016) values. The relaxation timescale can be defined as the ratio of viscosity to shear modulus, regardless of linear or nonlinear rheology (e.g., Freed et al., 2012;Fukahata & Matsu'ura, 2006). Therefore, the relaxation timescales of ∼0.0428 and ∼0.9754 years yield the viscosities of ∼4 × 10 16 and ∼9 × 10 17 Pa·s, respectively, for dashpots with a shear modulus of 30 GPa. Our estimated transient viscosity in the order of ∼10 16 Pa·s and the increase of magnitude by an order for steady-state viscosity are consistent with several postseismic models (Hu, Bürgmann, Uchida, et al., 2016;Moore et al., 2017;Muto et al., 2019;Qiu et al., 2018;Tang et al., 2019).

Viscoelastic Relaxation and Afterslip
Figures 4a and 4b illustrate the decomposition of 3-year (appropriate for transient periods) GNSS observations into the VR and AS components. Both the VR and AS impose deformations in the same direction at most stations for horizontal motion (Figure 4a) but show different behavior for vertical motion (Figure 4b). The relative contributions of VR and AS are illustrated by the ratio of VR to AS where VR/AS>1 and VR/AS<1 indicate the VR 9 of 13 and AS dominance, respectively (Figures 4c and 4d). The VR dominates the horizontal deformation in Tohoku region (Figure 4c) as reported by many studies (Agata et al., 2019;Dhar et al., 2022;Freed et al., 2017;Fukuda & Johnson, 2021;Iinuma et al., 2016;Muto et al., 2019;Suito, 2017;Sun et al., 2014). Such VR dominance is controlled by the average viscosity of the mantle (Luo & Wang, 2022). Using high-rate GNSS time series, Morikami and Mitsui (2020) also reported similar power-law decay of the horizontal motion after the Tohoku-oki earthquake implying the domination of non-Maxwellian VR.
Particularly for the coast of Ibaraki and northern Iwate (near station 940028), VR dominance in the horizontal motion disagrees with a few postseismic models that employ kinematic inversion of AS (Freed et al., 2017;Iinuma et al., 2016). However, the relative contributions of VR and AS in these areas are yet to be resolved by stress-dependent models. With no mechanical framework, the function model includes only weak interactions between the VR and AS components, which is indicated by small correlation coefficients (<0.5) among 1, 2, and (Table S4 in Supporting Information S1). Changes in the value of and AS impose insignificant changes in the correlation coefficients (Figures S4c and S5c in Supporting Information S1). Several postseismic studies reported no or weak interaction between VR and AS (Fukuda & Johnson, 2021;Tian et al., 2020Tian et al., , 2021. Such interactions may affect the region where the VR and AS components are similar in direction and magnitude. For example, our results predict an uplift farther inland from the coast (Figure 4d) by the AS component, in contrast to the results of Freed et al. (2017) and Fukuda and Johnson (2021).
Whereas the combined logarithmic-exponential functions explain the observed time series empirically, our decay functions tie themselves to the governing physics of VR and AS. Consequently, our conclusion on VR dominance in horizontal motion and AS dominance in coastal uplift agrees well with the recently developed viscoelastic models (Agata et al., 2019;Dhar et al., 2022;Freed et al., 2017;Fukuda & Johnson, 2021;Hu, Bürgmann, Uchida, et al., 2016;Iinuma et al., 2016;Muto et al., 2019;Sun et al., 2014;Wang et al., 2018) but contradict the analysis of Fujiwara et al. (2022) and Tobita (2016). Previous studies illuminated the dominance of power-law flow at least in the short-term (i.e., transient phase, Qiu et al., 2018;Tang et al., 2019) and both short-and longterm postseismic deformation (i.e., steady-state phase, Agata et al., 2019;Dhar et al., 2022;Muto et al., 2019). This nonlinear response owes to the dislocation creep in the mantle under the high-stress condition after the mainshock. The power-law decay function has a long tail (initial rapid with subsequent long-lasting deformation); thus, it explains the postseismic geodetic time series more accurately in both short-term and long-term deformation (Agata et al., 2019;Dhar et al., 2022;Freed et al., 2010Freed et al., , 2012Muto et al., 2019;Weiss et al., 2019). The logarithmic function representing an AS component shows asymptotic divergences in smaller slip rates (Barbot et al., 2009); thus, our can be a better substitution in predicting long-term evolution of the AS component (Barbot et al., 2009;Tang et al., 2019Tang et al., , 2020. Consequently, our function model produces subtle residuals even in 10 years of postseismic observations (∼0.1 cm/year; Figure S25a in Supporting Information S1) (c.f. ∼1 cm/yr for logarithmic-exponential function; Fujiwara et al., 2022). These residuals may vary negligibly (<1 cm) when the observed time series are smoothed with the moving average of different time windows (Figures S25b-S25c in Supporting Information S1). Notably, the contribution of VR and AS is easily defined in each time series as the VR and AS functions were built on the physics of distinct deformation mechanisms. Moreover, an efficient function model can be constructed by at least a 2-year time series (or 3.9 years for our best-fitted results) of 10 of 13 continuous postseismic observations. For example, our strategies are applicable to the postseismic deformation of the 1999 Chi-Chi earthquake (14-year observations, Tang et al., 2019), 2007 Bengkulu earthquake (5-year observations, Luo & Wang, 2021), 2010 Maule earthquake (6-year observations, Weiss et al., 2019;Peña et al., 2019), and 2012 Indian Ocean earthquake (3-year observations, . Because we can see the VR and AS components in GNSS time series of those earthquakes, our function model may illuminate the complex interplay of these mechanisms in other subduction zones and postseismic deformations of future earthquakes.

Conclusions
We propose a new function model (Equation 8) based on the laboratory-derived constitutive laws of rock deformation for the postseismic period of the 2011 Tohoku-oki earthquake. Our function model accurately predicts the decade-long GNSS time series following the 2011 Tohoku-oki earthquake; thus, it may also help forecast longterm GNSS time series. We ensured the precision and robustness of the function model by deploying Bayesian optimization. Moreover, it demonstrates the relative contributions of VR and AS with greater ease. Our results are supported by previous observations on the repeating earthquakes, high-rate GNSS time series, and recently developed stress-dependent postseismic models. Our proposed model may also illuminate governing mechanisms in the postseismic deformations of other megathrust earthquakes.