Spatiotemporal functional modeling of postseismic deformations after the 2011 Tohoku-Oki earthquake

Postseismic deformations continue to occur for a long period after major earthquakes. Temporal changes in postseismic deformations can be approximated using simple functions. Since the 2011 Tohoku-Oki earthquake, operating global navigation satellite system stations have been continuously accumulating a remarkable amount of relevant data. To verify the functional model, we performed statistical data processing on postseismic deformations due to this earthquake and obtained their spatiotemporal distribution. Moreover, we approximated the postseismic deformations over a relatively wide area with a standard deviation of residuals of 1 cm for approximately 10 years using a combined functional model of two logarithmic and one exponential functions; however, the residuals from the functional model exhibited a marked deviation since 2015. Although the pattern of postseismic deformations remained unaltered after the earthquake, a change in the linear deformation occurred from 2015 to date. We reduced the overall standard deviation of the residuals of > 200 stations distributed over > 1000 km to < 0.4 cm in the horizontal component by enhancing the functional model to incorporate this linear deformation. Notably, time constants of the functions were similarly applicable for all stations and components. Furthermore, the spatial distribution of the coefficients of each time constant were nonrandom, and the distribution was spatially smooth, with minute changes in the short wavelengths in space. Thus, it is possible to obtain a gridded model in terms of a spatial function. The spatial distributions of short- and long-period components of the functional model and afterslip and viscoelastic relaxation calculated using the physical model were similar to each other, respectively. Each time function revealed a connotation regarding the physical processes, which provided an understanding of the physical phenomena involved in seismogenesis. The functional model can be used to practical applications, such as discerning small variations and modeling for precise positioning.


Introduction
On March 11, 2011, an earthquake of magnitude (M) 9 hit the Pacific coast near Tohoku, Japan (hereinafter, the 2011 Tohoku-Oki earthquake) and caused coseismic crustal displacement of > 5 m in the horizontal direction and > 1 m in the vertical direction. Even after the earthquake, large postseismic crustal displacements of > 1 m continued to occur, depending on the location (Ozawa et al. 2011;2012). Tobita (2016) implemented simple logarithmic and exponential functions to approximate the time series of the postseismic deformations in this event to accurately predict its future trend over a wide area. The findings suggested that underground changes after a major earthquake are driven by the same physical processes.
Postseismic deformations caused by earthquakes primarily occur because of an afterslip at the plate boundary and viscoelastic relaxation of the upper mantle (e.g., Sun and Wang 2015); however, minor poroelastic rebound (less than several centimeters) may occur (Hu et al. 2014). As understanding the spatiotemporal distribution of postseismic deformations enables accurate estimation of seismogenic processes, such functional models are beneficial for forecasting seismic activity and understanding the physical mechanism of large earthquakes. In addition, they are useful for discerning small crustal deformations caused by other sources by removing large postseismic deformation. The functional model has been utilized by Ozawa et al. (2016), Sakaue et al. (2019), and others to discern crustal deformations other than postseismic deformations of large earthquakes. Tobita (2016) revealed a common time function at all observation stations over a wide area, thereby implying existence of a similar physical phenomenon. Therefore, the spatial distribution of physical parameters (spatial function) must be closely related to the physical phenomena distribution. In particular, the creation of a functional model using statistical processing, that is, an approach independent of a physical model, can be clarified by comparing it with the physical model. In this study, we obtained the spatial function for the first time and compared it with physical modeling results obtained using methods, such as the finite-element method.
Two major objectives are defined in this study. First, we applied our method to the crustal deformation data obtained from the continuous observation stations of the global navigation satellite system (GNSS) on the ground for approximately 10 years after the 2011 Tohoku-Oki earthquake, to verify its effectiveness in describing postseismic deformations and improve its accuracy. Second, we compared the obtained temporal and spatial functions with physical modeling results such as those for afterslip and viscoelastic relaxation to elucidate the physical mechanisms in the seismogenic field. Using our functional model, we could remove the postseismic deformation data from the postearthquake data with high accuracy and efficiently extract smaller nonpostseismic deformation components. The two major purposes of utilizing the spatiotemporal functional model developed in this study are as follows: (1) discerning small variations except those of the postseismic deformations, and (2) modeling for precise positioning such as data represented by the International Terrestrial Reference Frame (Altamimi et al. 2016) and semidynamic datum (Tanaka et al. 2007). In this modeling, it is important to predict time-varying terms of the reference frame in the future. We describe purpose (1) in "Use case of the functional model" section later in this article. However, we do not discuss related specific examples of purpose (2), because the accuracy of the model must be further verified for geodetic use.

Data and methods
We utilized the daily coordinates of the GNSS Earth Observation Network System (GEONET) F3 solution (Nakagawa et al. 2009) operated by the Geospatial Information Authority of Japan as the postseismic deformation data of the 2011 Tohoku-Oki earthquake. Stations were selected so that they are evenly distributed within the area where significant deformation occurred. Stations with large errors were excluded mainly owing to the observation environment, such as multipath of GNSS signal and vegetation. On the other hand, the data for the northern region are required to calculate the slip area that causes postseismic deformation. Therefore, in addition to Tobita's (2016) considerations, we selected stations from Hokkaido, to the northeast of Japan. In other areas, the postseismic deformations were minute. The Fukue GNSS station (station ID: 950462) in Nagasaki Prefecture, located approximately 1500 km west-southwest of the epicenter area, was utilized as the fixed reference station. The original GNSS time series comprises offsets caused primarily by large aftershocks. Therefore, in addition to that for the nine earthquakes considered by Tobita (2016), we rectified the coseismic variations of the following earthquakes: M7.4 November 22, 2016, Fukushimaken-Oki; M6.3 December 28, 2016, Ibarakiken-Hokubu; M6.7 September 6, 2018, Hokkaido-Iburi-Tobu; and M6.7 June 18, 2019, Yamagataken. The correction of the coseismic offset connects the trends before and after each earthquake and did not include the correction for postseismic deformations.
Major causes of postseismic deformations and their relationship with time functions were reported by Tobita (2016). Logarithmic and exponential functions efficiently approximate the observed postseismic deformations. Theoretically, the postseismic surface deformations determined by Maxwell viscoelastic rheology exhibited an exponential, such as relaxation phase (Hetland and Hager 2006). Moreover, postseismic deformations calculated using models with biviscous Burgers rheology are often described by logarithmic functions (Hetland and Hager 2006). The afterslip at the beginning of a postseismic deformation and the deformation related to viscoelastic relaxation can be approximated using a logarithmic function (Perfettini and Avouac 2004;Savage et al. 2005;Hetland and Hager 2006). Tobita (2016) presented seven combinations of the logarithmic (log) and exponential (exp) functions for use: log, exp, log + log, log + exp, exp + exp, log + log + exp, and log + exp + exp. Furthermore, the Akaike information criteria (AIC) value was calculated for each model utilizing the data for the entire period, from March 2011 to December 2020, to quantitatively evaluate the prediction performances. We implemented a mixed model comprising two logarithmic and one exponential functions (log + log + exp) that exhibited the best performance for the 2011 Tohoku-Oki earthquake, as reported by Tobita (2016) and represented by Eq. (1): where D(t) is each component of the postseismic deformation time series at a station; t represents the number of days after the earthquake; ln is the natural logarithm; b, e, and g indicate the relaxation time constants of the logarithmic or exponential functions common to all stations; and V is the steady velocity of each station before 2011. We applied the average velocity between April 1, 1997 and March 31, 2000 as the steady velocity V at each (1) Four stations marked as red squares were used for predicting the time series. Four stations marked as blue squares are used in Figs. 2, 3, and 6, and the stations marked as black squares are the other stations used in this study. Gray area represents land, and a seafloor topographic map is depicted. Gray arrow depicts the horizontal displacement rate observed at the sea bottom after the 2011 Tohoku-Oki earthquake (Honsho et al. 2019). Isodepth contours of the Pacific slab were estimated by Kita et al. (2010), Hasegawa (2006), andNakajima et al. (2009) station. This equation could be utilized to predict the short-and medium-term trends and postseismic deformation changes that vary by location.
The time constants are estimated based on the data of four stations (Minase, Miyako, Yamoto, and Choshi) ( Fig. 1) for two reasons. First, among all the GEONET stations, these four stations possessed high-quality GNSS data considering low noise and seasonal variations. Second, although all four stations subsided during the 2011 Tohoku-Oki earthquake, their vertical postseismic deformations differed (Additional file 1: Fig. S1). The Minase station experienced rapid subsidence after the earthquake. Moreover, the Miyako station subsided rapidly after the earthquake but was reversed to a slow uplift. The Yamoto station uplifted rapidly after the earthquake. Eventually, the Choshi station experienced rapid uplift immediately after the earthquake; however, its vertical movement reduced. All four stations exhibited representative regional crustal deformations; thus, they provided strong constraints for our nonlinear least-squares calculations. We implemented multiple starting points, indicating that the nonlinear least-squares computations start with different initial values. Therefore, verifying the residuals and AIC is necessary while changing the initial values during solution determination.
The coefficients of a, c, d, and f, unique to each station and component, are calculated using the least-squares method. The postseismic deformations include numerous periods, from short-period afterslip components to long-period viscoelastic relaxation components (Sobrero et al. 2020;Sun and Wang 2015). These deformations are expressed as a sum of two logarithms with different time constants and one exponential function (Eq. 1). The nonlinear least-squares method is utilized to determine the time constants b, e, and g, as it is computationally timeconsuming and may provide local solutions. Therefore, owing to the commonality of the time function, after careful and rigorous determination of the time constants at the four observation stations, it is reasonable to determine the coefficients a, c, d, and f at each station using the linear least-squares method, which can determine the solution in a short time. This can be expressed as follows: where time_f(t) is the time function (ln(1 + t/b), ln(1 + t/e), -exp(-t/g), and t) in Eq. (1), and space_f(x, y) is the spatial function (a, c, d, f, and V) at the coordinate (x, y) for each time function. This formula is similar to that used in the principal component analysis (PCA) and independent component analysis (ICA). Munekane (2012) solved the time function statistically using PCA. In these methods, the time functions are separated into orthogonal or independent functions; however, the form of the functions is not given in advance. Our method determines the form of the function in advance as log + log + exp. As a result, it is presumed that the residuals between the observed values and the model values at each station will not represent the minimal possible values; however, as long as the residuals are within an acceptable range, a simplified function can be employed to model the observed values at any station. Furthermore, these methods may provide insight into future predictions.
Consequently, in this study, we verified the method used to limit the time function to logarithmic and exponential functions and discuss the spatial function using 222 stations. Table 1 summarizes the calculated results of the time constants, b, e, and g for the fitting periods of 2.0, 3.9, 5.8, and 8.9 years. These fitting periods are not evenly spaced, because the calculations were performed during data accumulation. In addition, a clear difference in data fitting was observed before and after February 2015 and the fitting up to February 2015 (3.9 years) was distinctively added. Here, for example, a 2.0-year model is one in which we assign 731 days (2.0 years) after the 2011 Tohoku-Oki earthquake as the "fitting period, " and the number of days between the 732nd and after as the "prediction period. " The fitting period started on the day after the earthquake. The time constants b and e tended to increase as the fitting period increased. In contrast, the values of g for 2.0-and 3.9-year models revealed slight changes; however, the values of g for the next 5.8-year model markedly increased. The time constant of 450,000 days for g of the 5.8-and 8.9-year models was > 1200 years; thus, it is unlikely that the time constant can be appropriately calculated with a limited number of years for fitting. Figure 2 illustrates the time series predicted by the functional model with time constants summarized in Table 1 for different fitting periods for the Kawai and Shizugawa stations, which are relatively close to the epicenter of the 2011 Tohoku-Oki earthquake (see Fig. 1), and the Kaminokuni and Nadachi stations, which are relatively distant from the epicenter (see Fig. 1). These plots demonstrate that model parameters estimated from the time series of the four stations can be applicable to other stations. Moreover, Additional file 1: Fig. S1 illustrates the time series for Minase, Miyako, Yamoto, and Choshi. Figure 2 and Additional file 1: Fig. S1 depict the extrapolated predictions up to 2026.5. Furthermore, we evaluated the time series of the short fitting periods of 2.0 years and 3.9 years. At all the stations, the 2.0-year model of the north-south (NS) component revealed a smaller value than the 3.9-year model. Kawai, c Nadachi, and d Shizugawa (Fig. 1). The three components of east-west (EW), north-south (NS), and up-down (UD) were drawn for each station from the day after the 2011 Tohoku-Oki earthquake to 2026.5. Black dots represent the daily observations. Dashed lines depict the time series of predicted values for each fitting period (2.0, 3.9, 5.8, and 8.9 years) Additional file 1: Figure S2 illustrates the distribution of the short-period logarithmic term [the first term in Eq. (1)] for the 2.0-year and 3.9-year models on March 11, 2014, three years after the earthquake, and the difference between them. Although the two curves are consistent, the eastward displacement of the 2.0-year model was smaller near Miyako and south of Yamoto; southward and upward displacements of the 2.0-year model were generally larger than those of the 3.9-year model. NS components of the 2.0-year model were affected by local variations of the fixed reference station Fukue in 2013. Effects of the reference station could be easily detected due to the appearance of same variations at all stations. The effects of the different time constants are described in the Discussion. Notably, in this process, it is difficult to avoid errors in the predictions caused by decreased fitting data. Despite the noise in the fixed reference station, Fig. 2 illustrates that the time series of the 2.0-year and 3.9-year models are parallel in the future, and their long-term trends are similar. This occurs, because the time constant of the long-period exponential function is almost similar for both time series (~ 10 years); Vt, which presumably occurs due to the plate motion, becomes relatively dominant approximately 10 years after the earthquake. Thus, the long-term trend is appropriately included in both models and does not change over time. In addition, although the 2.0-year model is susceptible to noise due to the short period over which the values were obtained, for a prediction period of 10 years, the long-term residual is only several centimeters, which is sufficient for some purposes that can tolerate an error of a few centimeters, because overall postseismic deformation can exceed 1 m.

Time function
Furthermore, the 8.9-year model, which has the longest fitting period, represents the time series with the smallest residuals. Although the 2.0-and 3.9-year models and the 5.8-and 8.9-year models are aligned and yield similar curves, a notable difference is observed between these two groups, thereby indicating that a systematic change occurred between 3.9 and 5.8 years. This is consistent with the fact that the exact value of g cannot be calculated for the 5.8-and 8.9-year models, as summarized in Table 1. Figure 3 illustrates the residuals based on the predictions for the 3.9-year model (until February 2015) to investigate the reason for this change in trend with the fitting period. The residuals of the four stations used to calculate the time constants are depicted in Additional  Figure S1d illustrates no difference between the time series of 5.8 years without the effect of this gap and those of 8.9 years with the effect. Furthermore, in Fig. 3 and Additional file 1: Fig. S3, the predicted values correspond with the observed values within 1 cm until February 2015; however, the residuals accumulated in a certain direction after the M6.9 Sanriku-Oki earthquake on February 17, 2015 (hereinafter, the 2015 Sanriku-Oki earthquake) and the M6.8 Miyagiken-Oki earthquake on May 13, 2015 (hereinafter, the 2015 Miyagiken-Oki earthquake). This change was observed in almost all the stations included in this study (Fig. 1). Thus, since 2015, the residuals increased at a constant rate from Hokkaido to central Japan, thereby indicating that a new phase of postseismic deformation was initiated and the original functional model was ineffective. In contrast, the value of g markedly increased. Therefore, a period longer than 3.9 years should not be adopted as the fitting period to use the functional model in Eq. (1).
Moreover, the east-west (EW) components with the most pronounced changes after mid-2015 reveal that the linear changes in Kaminokuni and Nadachi commenced in the beginning of 2015, whereas very slow changes in Kawai and Shizugawa occurred from February to July 2015. Shizugawa revealed an eastward displacement, caused by the M6.2 Miyagiken-Oki earthquake on April 20, 2020. The common change of approximately 0.5 cm in the NS component in July 2013 was caused by the fixed reference station Fukue. Figure 3 illustrates a linear component appearing after the beginning of 2015 at all the stations and components. A straight line is fitted using the least-squares method for the 5-year period from December 2015 to December 2020 (exactly five years considering the existence of an annual component) to model this new component as a function. Notably, there is a possibility that the new component after 2015 is a function with nonlinear properties; however, none of the changes suggest a function more complex than the linear function. Hence, the new component can be expressed with minimal simple functions. In this study, we propose Eqs. (3) and (4) as a modified functional model: (3) where c' and v are constants for each observation station and component calculated by linear fitting. Here, t 0 denotes February 17, 2015 (t 0 = 1439) and t 1 denotes July 1, 2015 (t 1 = 1573, it is the reading from Fig. 4).  Fig. 1, were less than 1 cm, and were followed by small postseismic deformations. Therefore, it is difficult to separate these from the postseismic deformation of the 2011 Tohoku-Oki earthquake. As a huge direct influence of the local earthquakes ( Fig. 5a 1 ) occurred during this period at the aforementioned stations ( Fig. 4a and c), we select 18 stations ( Fig. 5a) that are not remarkably closer to these earthquakes, and have exhibited larger variations during this period. Figure 4d illustrates the average (stacking) of the 18 stations around northern Honshu (Tohoku region).
Major characteristics of the changes that occurred around 2015 are as follows. First, the EW component was flat before February 2015, and the functional model in Eqs. (1) and (3) was applicable. Second, coseismic deformations associated with the 2015 Sanriku-Oki earthquake, the 2015 Miyagiken-Oki earthquake, and the M6.7 Urakawa-Oki earthquake on January 14, 2016 (hereinafter, the 2016 Urakawa-Oki earthquake) appeared at observation stations near each earthquake. Third, a linear variability with a constant trend was observed between February and July 2015 (pink area in Fig. 4) in a wide area. Fourth, the variability stabilized and changed to a constant rate after January 2016 for all stations. As illustrated in Fig. 4a and c, the second coseismic variation is large in the vicinity of the epicenters; however, in most locations, the third gradual variation from February to July is prominent (Fig. 4). The changes between February and July 2015 include the coseismic and postseismic deformations of the 2015 Sanriku-Oki and 2015 Miyagiken-Oki earthquakes that occurred during this period. The postseismic deformations of these local earthquakes cannot be separated, because they are smaller than several millimeters and their duration remains unclear; thus, the changes associated with these earthquakes are treated as one from (4) 1 In this paper, horizontal distributions such as vertical components in all figures are drawn by grinding and using splines (tension is 0.5) by Smith and Wessel (1990).
February to July. Therefore, the change in displacement velocity observed in July may be apparent; however, a constant velocity component [v in Eq. (4)] was observed after July. The total displacement due to these two fluctuations (the second and the third characteristic changes) is calculated as the difference between the values on July 1, 2015 (t 1 ) in Eq. (4) and the value on February 17, 2015 (t 0 ) in Eqs. (1) and (3). This displacement corresponds to the dashed lines with arrows drawn (Fig. 4) on July 1, 2015. In this study, we refer to the gap as the "gap in 2015. " Although the changes caused by the 2016 Urakawa-Oki earthquake did not occur in 2015, they are closely related to the vt term in Eq. (4), which is described in detail in the Discussion section. In addition, owing to the small value of v at observation stations affected by the 2016 Urakawa-Oki earthquake (Fig. 5b), the error in the gap caused by the different timing of the coseismic deformation is also small. Therefore, we include these in the series of gap in 2015 in our discussion. Thus, the new component of Eq. (4) from 2015 can be described as follows: where gap 2015 is the gap in 2015, and this equation is applicable after t 1 . The variation from 2015 to early 2016 including the occurrence time of three earthquakes cannot be modeled as a simple function. Therefore, for convenience, gap 2015 presumably occurred linearly between February 17 and July 1, 2015 (t 0 ≤ t < t 1 ); Eqs. (3) and (4) are connected so that no step-change is observed in the functional model. Therefore, while using the improved functional model, discontinuities were included from February 2015 to January 2016.
Due to the variability and errors in the daily GNSS observations, the exact date and time of the beginning of the third characteristic change cannot be determined; however, it is estimated to have started around the same time as the 2015 Sanriku-Oki earthquake (February 17, 2015; Fig. 4).
To validate Eqs. (3) and (4), the spatial distribution of the gap in 2015 (gap 2015 ) and v were analyzed ( Fig. 5a and 5b, respectively). The vector distribution was drawn without smoothing the values for each station, and the difference between the values of each component from the neighboring stations was < 1 cm. Depending on the station, the error in the GNSS observation was approximately a few millimeters (Figs. 3 and 4). Figure 5 illustrates a uniform spread of variability from Hokkaido to central Japan. Applying a straight line, vt, to each component of each station may result in over-fitting of the line to local phenomena that are unique to each station. However, both gap 2015 and v were spatially distributed in a smooth manner (Fig. 5) and the systematically occurring phenomena were captured over a large area, which validates the assumption of Eq. (4).  Figure 6 illustrates the time series for each station using Eqs. (3) and (4). The value of the 3.9-year model was used as the overall time constant, and vt, which is a linear variation after 2015, was added. The initial value of each term was drawn as a graph starting from zero, and a constant was added to the exponential term [the fourth term in Eqs. (3) and (4)  Kawai, c Nadachi, and d Shizugawa (Fig. 1). The three components of east-west (EW), north-south (NS), and up-down (UD) are drawn for each station from the day after the 2011 Tohoku-Oki earthquake to 2026.5. Gray dots indicate the daily observations used for the fitting period of the 3.9-year model. Black dots represent the daily observations during the prediction period. Red lines are the predicted values based on Eqs. (3) and (4), and other lines represent specific terms in Eqs. (3) and (4) (3) and (4) was optimal. Figure 6 illustrates the total functional model as well as each term of the model. In Fig. 6, Eqs. (3), (4), and the 3.9-year model are used. Each term is well-defined in the overall function in the horizontal component; however, in the UD component, the function varied depending on the observation station and component, and the exponential term is dominant. Tobita (2016) predicted that although the vertical position of the Yamoto station returned to its pre-earthquake elevation in 2020 (in Additional file 1: Fig. S4, UD zero depicts the elevation just before the 2011 Tohoku-Oki earthquake), the displacement rate of the overall crustal deformation has not returned to the pre-earthquake trend [the Vt term in Eq. (1)] even 10 years after the earthquake. Figure 6 also depicts the extrapolated predictions for the near future, which do not predict a return to the pre-earthquake trend even after several more years.

Spatial function
In this section, we present the spatial distribution of the components of the functional model. Figures 7 and 8 illustrate the distribution of the two logarithmic, exponential, and steady velocity Vt terms three years after the earthquake (March 11, 2014). As determined by Eq. (3), the term vt is zero for this time and changes after July 2015 (Fig. 5b). Compared to the first 3 years, the shortperiod logarithmic, long-period logarithmic, exponential, and Vt terms of the 5-year post-earthquake period were larger by 1.08, 1.22, 1.52, and 1.67 times, respectively; thus, suggesting that the contribution of the longer period terms is larger.
The spatial distribution of all the terms is systematic rather than random. The spatial wavelengths were shorter for smaller relaxation time constants for each term, at the shortest being approximately 100 km. The distribution is spatially smooth and little change was observed in the short wavelengths in space for each component, although it cannot be easily expressed as a mathematical function. Therefore, a gridded spatial function can be obtained by simple (e.g., linear) spatial interpolation, for points that do not require observation.
These spatial distributions are caused by physical phenomena and are separated by time functions with different time constants, indicating that the postseismic deformations are caused by phenomena that are related to time and underground location. For example, the short-period logarithmic term revealed a more complex spatial distribution than that of the other terms, which indicated that the spatial distribution was influenced by a phenomenon occurring in the shallower and narrower areas (Fig. 7a).

Validity of identical time constants
We assumed that the time constants are the same for all stations and all components. The four stations used to calculate the time constants represent the characteristic cases of regional variations in the postseismic deformations. For this reason, the solution of the nonlinear leastsquares method does not converge from a set of arbitrary stations and to obtain a stable solution, the combination of stations must be carefully selected.
In this section, we verify the time constants considering that these constants (Table 1) with various fitting periods of 2.0 and 3.9 years exhibit different values. In our method, first we calculate the time constants [b, e, and g; in the time function in Eq. (2)] using the four stations; second, we calculate the stationwise coefficients [a, c, d, and f; the spatial function in Eq. (2)]. In general, the time period for determining the time constants and that for determining the stationwise coefficients were identical. Here, we examine various combinations of the time periods that are listed in Table 3.
We prepare one functional model time series (same as shown in Fig. 2) in which the spatial function with the fitting period of 2.0 years (up to March 2013) and time constants of the 2.0-year model (Table 1) (a, c, d, and f)  function with the fitting period of 3.9 years and the time constants of the 3.9-year model [S3.9T3.9]) as a reference (Fig. 2). Residuals of these time series are illustrated in Fig. 9 for the same stations as in Fig. 2, and in Additional file 1:   Fig. 9). As a reference value, we use the 3.9-year model (the spatial functions with the fitting period of 3.9 years and the time constants of the 3.9-year model [S3.9T3.9]). b Same as a but for the spatial functions with the fitting period of 3.9 years and the time constants of the 2.0-year model [S3.9T2.0] (red line in Fig. 9). c Histogram of the ratio of the spatial function for the different time constants of the 2.0-year model [S3.9T2.0] to that of the 3.9-year model [S3.9T3.9]. Green color depicts the short-period logarithmic term (a), red color indicates the long-period logarithmic term (d), and blue color represents the exponential term (f) in Eq. (1) (see Additional file 1: Fig. S6 for more information) varies by region, for example, large horizontal changes exist, northwest of Miyako. The residuals in Fig. 10b are markedly higher than those in Fig. 10a; however, the distribution is much smoother spatially than that in Fig. 10a. To understand the reason for this difference, we display histograms (frequency rate) of station numbers for the ratio of each term [the short-and long-period logarithmic terms, and the exponential term in Eq. (1)] between the 2.0-year model [S3.9T2.0] and the 3.9-year model [S3.9T3.9] (the reference model) time constants in Fig. 10c. The short-and long-period logarithmic terms and the exponential term vary by a constant factor of 0.88, 1.05, and 1.16, respectively, as the 3.9-year model [S3.9T3.9] is changed to the 2.0-year model [S3.9T2.0] (see Additional file 1: Fig. S6 for the calculation of the ratio), that is, when the time constants (b, e, and g) vary, the distribution of the spatial function (a, d, and f) is altered at a constant rate for all stations. As the period is shorter for all time constants in the 2.0-year model than in the 3.9-year model (Table 1), the spatial functions corresponding to the time functions with longer periods are larger, because the overall functional model tries to fit the observed values.
Thus, two primary findings are as follows: (1) when the time constants vary, the distribution of the spatial functions is altered at an identical rate at all stations, absorbing the change in the total functional model; and (2) as a result of (1), even if the time constants deviate slightly from the optimal solution, the overall residuals of the functional model can be markedly reduced by changing the spatial functions.
Using the time constants obtained from the subset of the stations that comprehensively include the characteristic changes in the postseismic deformation, the overall functional model can be well-fitted to the time series for each station, although the time constants initially differ for each station. Notably, this does not indicate that the time constants should be the same. Alternatively, this highlights the effectiveness of using the same time constants for each observation station and each component as a means of predicting the postseismic deformation in a simple functional model.
Next, we discuss the slight deviation of the time constants from the optimal solution [the primary finding (2)]. Tobita (2016) reported that a single logarithmic function approximates the sum of other multiple logarithmic functions with different time constants, which explains why the whole can be approximated by a function with a uniform time constant, even if the individual observation stations comprise components with different time constants. In contrast, Hetland and Hager (2006) reported that in layered viscoelastic models, multiple mechanical timescales of the postseismic deformation are present; however, not all of these arise as distinct phases of postseismic relaxation observed at the surface. This may presumably reveal the reason why the whole data can be approximated by one set of uniform time constants.

Comparison between functional and physical models
As the spatiotemporal model reveals time functions with various time constants and its spatial distribution is smoothly spread, in this study, we explore the meaning of the time constant and other factors by comparing them with those of a physical model. Freed et al. (2016) and Suito (2017) constructed 3D viscoelastic relaxation models using the finite element method for the postseismic deformation of the 2011 Tohoku-Oki earthquake. The difference between the viscoelastic relaxation model and the observed data may represent the afterslip. Figure 11a and b exhibit the shortperiod logarithmic term [the first term in Eqs. (3) and (4)] and the afterslip component based on Suito (2017), respectively, as of March 11, 2014. Here, we assume that the observed data include the afterslip, viscoelastic relaxation, and constant velocity (Vt term). Therefore, the afterslip component from the physical model is the observed data minus the Vt term and the viscoelastic relaxation model component.
In Fig. 11a and b, the eastward and subsidence displacements near Miyako and the southeastward and uplift displacements near Choshi are quantitatively in accordance with the functional and physical models. In the central part, between 37° and 40° north latitude, as displayed with dotted line "Volcanos" in 11a and b, a north-northeast-to-south-southwest (NNE-SSW) trending band of relative subsidence (the UD components are smaller than those of surrounding areas) is identified over several hundred kilometers. Takada and Fukushima (2013) highlighted that the volcanic regions deformed and subsided in response to stress changes associated with the 2011 Tohoku-Oki earthquake. The subsidence band in Fig. 11a and b may be due to the inclusion of a postearthquake volcanic subsidence. The physical model is based on the differences in the observed data; therefore, the subsidence band is a manifestation of the original observed data and not that of the modeled data. The two models are consistent in several instances; however, near Yamoto, the functional model reveals a large eastward displacement, whereas the physical model exhibits almost zero horizontal displacement.
The remaining two terms of the functional model, that is, the long-period logarithmic term [the third term in Eqs. (3) and (4)] and the exponential term [the fourth term in Eqs. (3) and (4)], one-by-one, are not related to the physical model (afterslip and/or viscoelastic relaxation) of Suito (2017). The sum of these two and the viscoelastic relaxation physical model are illustrated in Fig. 11c, d, respectively. The horizontal displacement is generally eastward, with a southward displacement in the north and a northward displacement in the south, and the vertical displacement is uplifted along the Pacific coast and subsided along the Sea of Japan coast. In contrast to Fig. 11a and b, the viscoelastic relaxation physical model reveals a large eastward displacement near Yamoto. The differences between Fig. 11a and b (the same as between Fig. 11c and d) are depicted in Fig. 11e, and the difference near Yamoto is large. Moreover, the difference between the two models is likely caused by the shorter time phenomenon in viscoelastic relaxation near Yamoto than in the other areas, and this subsequently affects the short-period logarithmic term of the functional model, which primarily represents the afterslip component ( Fig. 11a and b). This process can be explained as follows: the coseismic slip zone of the main shock of the 2011 Tohoku-Oki earthquake is located close to Yamoto, and the slip zone extends to just near Yamoto (Iinuma et al. 2016;Ozawa et al. 2012;Suito 2017). Therefore, after the main shock, a smaller afterslip  (2017) is observed in this area. This explains the relatively large displacement due to viscoelastic relaxation. According to Suito (2017), the effective viscosity of the mantle wedge, which is located just below the land observation stations, is approximately one-fifth that of the oceanic mantle. Such a low effective viscosity makes ductile deformation in this area more rapid, which leads to a shorter viscoelastic relaxation period; thus, it is certain that this effect occurs near Yamoto. Figure 11 displays a clear difference in the horizontal component between the functional model and the physical model; however, the difference is much smaller in the UD component than the horizontal component, because the pattern of spatial distribution is similar to each model. As the UD components exhibit a large contribution to the exponential term (Fig. 6), the physical model can be improved by increasing the contribution of the long-period term in the UD component. In addition, the effect of viscoelastic relaxation must be considered to infer interseismic vertical displacements (Nishimura 2014).
In this study, we reveal that the viscoelastic relaxation term in the physical model is the sum of the longperiod logarithmic term and the exponential term in the functional model. Moreover, the time constant of the exponential term is > 20 times longer than that of the long-period logarithmic term (Table 1). Therefore, the frequency range of the viscoelastic relaxation is wide, ranging from a few days to thousands of days. The long-period logarithmic and exponential terms are necessary as they effectively capture a spectrum of viscoelastic relaxation times, which occur naturally owing to variation in the viscosity and elastic properties.
As both the models simplify the phenomena, a fundamental difference is that a single physical phenomenon can be divided into multiple time functions, or conversely, multiple physical phenomena can be included in a single time function; hence, the functional and physical models differ. A physical model can be constructed according to the time constant, which may lead to more advanced modeling by incorporating the differences in underground properties.

Slip distribution on the plate surface
Next, using geodetic data inversion to estimate slip at a plate interface, we discuss the extent and location of the movement of the source area. This study has two limitations. First, viscoelastic relaxation is complex in terms of the location and the direction of motion at each location (e.g., Agata et al. 2019). In contrast, slip on a plate surface can be simulated with constraints on the location and direction of the slip. Therefore, we primarily deal with phenomena that can be considered as a slip on the plate surface. The second limitation is that only data from the land side are available, and we are trying to build a model without using data from the ocean side, because the seafloor crustal movement data do not present adequate time series. Hence, offshore slip that is far  Kita et al. (2010), Nakajima and Hasegawa (2006), and Nakajima et al. (2009) from land cannot be estimated from the onshore data alone and the certainty of the model decreases as the distance from land increases (Iinuma et al. 2016).
We estimate the slip distribution on the plate by applying geodetic inversion using the Markov chain Monte Carlo method (Fukuda and Johnson 2008). This method efficiently obtains the desired probability distribution by generating numerous solutions (samples) that follow the posterior probability distribution. The solution does not minimize the residuals, but produces a set of solutions that can explain the observed values within a given error range. We assumed a homogeneous elastic half-space body, constrained the direction of slip to the opposite direction of plate convergence, and used Laplacian smoothing. Additional file 1: Figure S7 illustrates the results from a checkerboard resolution test to evaluate the reliability of the calculation results of the inversion implemented in this study. As mentioned, the certainty of the slip model decreases as the distance from the land increases; however, for plate surfaces deeper than 20 km, this model is extremely useful. Figure 12a depicts the results of the inversion and slip on the plate surface as the driving force for the shortperiod logarithmic term (Fig. 7a) at three years after the earthquake. As discussed above, a quick relaxation in the mantle wedge appears similar to the afterslip; thus, it is modeled together in the same single deformation time scale. The moment magnitude corresponding to this term in the three years is M w 8.3 and it is larger than that of the individual aftershocks that occurred during this period. While calculating the moment magnitudes in this study, we assumed a rigidity of 30 GPa. Figure 12b and c depict the slips on the plate surface estimated from the term gap 2015 in Eq. (5) (Fig. 5a) and term vt in Eq. (4) (Fig. 5b), respectively. The symbols from B to E in the figures indicate the locations of the characteristic variations illustrated in Fig. 12c. Although there are errors due to the aforementioned two restrictions (the viscoelastic relaxation is not considered and there is no offshore observation) in all the three parameters of Fig. 12, these mutual comparisons are meaningful, because the data are analyzed under similar conditions. For example, it is worth discussing that there is no slip at B in Fig. 12a compared to b and c, and that there is a slip spread northeast of C in Fig. 12b.
Although it is possible to calculate the long-period logarithmic term and the exponential term as a slip on the plate surface, we do not discuss the calculations, because they primarily comprise viscoelastic relaxation components, and not afterslip. Figure 12b shows a unique slip from February 2015 to January 2016 apart from the original postseismic deformation that started shortly after the 2011 Tohoku-Oki earthquake, and includes the coseismic and postseismic deformations caused by the three earthquakes with a magnitude of > 6.5. Figure 12b illustrates that these variations can be attributed to two locations: Miyagiken-Oki (C), which is associated with the occurrence of the 2015 Sanriku-Oki and 2015 Miyagiken-Oki earthquakes, and Urakawa-Oki (B) to the north, which corresponds to the 2016 Urakawa-Oki earthquake.

Gap in 2015
The moment magnitude of the slip is M w 7.0 for the Urakawa-Oki (B) area and M w 7.3 for the Sanriku-Oki to Miyagiken-Oki (C) area, indicating that both slip amounts are large and cannot be explained by the three earthquakes and their postseismic deformation alone.
The gap 2015 and the 2015 Sanriku-Oki earthquake of February 17 would have commenced simultaneously (Fig. 4). The focal mechanism of three earthquakes revealed a reverse fault type with a compression axis in the plate motion, and they occurred at the boundary between the Pacific and the continental plates (Earthquake Research Committee 2015a, b;2016). No solid evidence has revealed whether the 2015 Sanriku-Oki earthquake triggered the gap 2015, or the earthquake itself was triggered by some phenomenon; however, Baba et al. (2020) reported that the very low frequency earthquakes (VLFEs) were activated by the 2015 Sanriku-Oki earthquake and an SSE occurred in Sanriku-Oki around February 2015 (Honsho et al. 2019). Therefore, the 2015 Sanriku-Oki earthquake presumably triggered gap 2015 . In either case, based on the location, direction of the slip, and timing of the change, the three earthquakes and the new slip that started in February are related, and it would be appropriate to treat them as a single gap 2015 .

Linear slip after 2015
The temporal changes in the residuals (Fig. 3) based on Eq. (1) between the observed values and functional model occurred systematically and simultaneously over a wide area after 2015 and continued at approximately constant speeds, which are expressed as vt term in Eq. (4). To analyze the cause of this phenomenon, a model of slip rate on the plate surface is illustrated in Fig. 12c, showing a few centimeters of slip per year over larger areas than that of the afterslip. It is uncertain whether the slip in this term is partially caused by viscoelastic relaxation that occurs in the upper mantle, rather than that only on the plate surface; however, the temporal change shown in Fig. 3 proceeds at a constant speed, which is unlikely to occur owing to viscoelastic relaxation, because the change started several years after the 2011 Tohoku-Oki earthquake, and the deformation was observed even in the Urakawa-Oki, located distally from the main shock. Therefore, the constant slip on the plate surface likely continues.
The location of the new slip can be roughly divided into two: Miyagiken-Oki (C) to Fukushimaken-Oki (D) and Urakawa-Oki (B). In the Miyagiken-Oki to Fukushimaken-Oki area, the slip spreads both in the afterslip area (Fig. 12a) and deeper (western) area of gap 2015 and is considered to occur with a constant velocity to add to the afterslip. In contrast, no significant slip has occurred in the Urakawa-Oki area before 2015 but has occurred after 2015. In addition, the total moment magnitudes of these areas from 2015 to 2020 are integrated to be M w 7.4 for B and M w 7.8 for C and D, and they are larger than that of each earthquake that occurred in all individual areas.
Next, the relationship between the new linear deformation vt term after 2015 (Fig. 5a) and steady velocity Vt term that existed before 2011 (Fig. 7d) is discussed.
These v and V have almost opposite directions, and the ratio of their magnitudes is maximum of approximately 0.4. As both present an almost stationary velocity, we can interpret v as a change in velocity V; however, the spatial extent of the distributions of V (Fig. 7d) and v (Fig. 5b) differs, implying that only a small part of the change in the physical mechanism is responsible for v. Moreover, at this point, we could not identify any major event that led to the linear deformation in 2015. Moreover, a huge event likely has a time constant as a secondary after-effect fluctuation, implying that the event may not have a constant velocity as in v. Notably, no major events were noticed in 2015, and a silent shift to a linear deformation occurred.
The underlying cause for this phenomenon, which commenced silently and continued without any major events, remains unknown. However, the linear nature of the phenomenon leads to a hypothesis stating that the site probably returned to its "original" tectonic state long before the 2011 Tohoku-Oki earthquake. This implies that the steady velocity V before the earthquake was an anomaly (Nishimura 2014) that likely recovered after the earthquake and returned to its original state in 2015. Further observations and analyses are required to verify this hypothesis.
Although this slip represents a new "unsteady" linear deformation for the functional model, it has been observed for approximately six years; thus, we may assume that the slip is not a linear deformation but an SSE with an extremely long period (Mallick et al. 2021). Along the Japan Trench, shallow slow earthquakes are frequently observed, whereas long-period SSEs are rarely observed. Furthermore, the slip area is classified as a quasi-stable slip zone with numerous repeating earthquakes (e.g., Obara 2020). Therefore, this slip presumably differs from an SSE owing to its large horizontal extent, long duration, and uniform slip rate. Figure 13 illustrates a series of changes in the location of the slip zone on the plate surface from before 2011 to 2020 based on Fig. 12 and similar analyses. Figure 13a, b, d, and e show two phenomena in the same figure for comparison, one before and one after, in each period. Figure 13a illustrates the slip distribution of the steady velocity Vt term on the plate surface (Eqs. 3 and 4) (> 10 cm/year) before 2011 and the main shock (coseismic) slip of the 2011 Tohoku-Oki earthquake (Ozawa et al. 2012;> 10 m). The slip of the steady velocity Vt term is most likely not caused by the slip on the plate surface but by that on the continental plate (Japanese archipelago side), which is being pushed as the plate boundary interface is locked. This figure depicts a rough distribution of qualitatively strong adherence sites and clarifies the relationship between V and v. The 2011 Tohoku-Oki earthquake occurred on the shallow side of the plate surface with a strong adhesion before 2011 (e.g., Uchida and Matsuzawa 2013). Figure 13b illustrates the coseismic slip of the 2011 Tohoku-Oki earthquake and the slip of the short-period logarithmic term (> 2 m) in three years after the earthquake, primarily representing the afterslip. The afterslip occurred adjacent to the deep side of the coseismic slip of the earthquake. The results obtained for this slip area are consistent with those of Ozawa et al. (2012) and Iinuma et al. (2016); however, notably, this slip distribution includes the viscoelastic relaxation with short period time constants, as observed during comparison with the physical model. Figure 13d illustrates the slip of gap 2015 (> 5 cm), which occurred during the afterslip, immediately following the main shock. Furthermore, the northeastern side of this area indicates the location of the SSEs, VLFEs, and other events or asperities (Baba et al. 2020;Nishikawa et al. 2019). The slip estimated in this study occurred on the northeastern side of the afterslip area partially overlapping it. The slip of gap 2015 is also estimated in the Urakawa-Oki (B) area, which is far from the afterslip area. Considering the time series of the gap 2015 , the VLFEs and SSE occurred near the 2015 Sanriku-Oki earthquake around February 2015 (EQ2 in Fig. 13d) (Baba et al. 2020;Honsho et al. 2019), followed by the occurrence of a linear deformation event from February to July over a wide area, denoted as EQ2 in February, EQ3 (the 2015 Miyagiken-Oki earthquake) in May, and EQ1 (the 2016 Urakawa-Oki earthquake) in January 2016.

Evolution of slip zone on the plate boundary interface
Considering the variations in the Urakawa-Oki (B) area, Fig. 4a illustrates that the slow variation commenced, albeit slightly, in February 2015, although EQ1 occurred in 2016. Thus, EQ1, EQ2, and EQ3 are associated with the series of events, thereby revealing a new observation for the slip areas, as the result cannot be identified without the residual from the precise functional model obtained in this study. Areas B to E represent the locations of the characteristic displacement in e. F represents an area with a small displacement in any of the figures, and A represents the location at which the maximum displacement was estimated to have occurred in the main shock of the 2011 Tohoku-Oki earthquake. Isodepth contours of the Pacific slab were estimated by Kita et al. (2010), Hasegawa (2006), andNakajima et al. (2009) Figure 13e depicts the linear deformation vt term (> 2 cm/year) that occurred simultaneously as gap 2015 and continued to date. The Urakawa-Oki (B) area slightly shifted southwest of the gap 2015 distribution. Moreover, the deformation of the Boso Peninsula offshore (E) area was affected by the SSE in 2018 and did not occur as a linear continuous variation within this period (Ozawa et al. 2019). It is difficult to extract such changes (B, C, and D) without carefully removing the postseismic deformation, and thus, this is the first study to detect this slip area.
Areas B, C, D, and E are nearly identical to the locations in Fig. 13a and e; however, opposite slip direction was observed before 2011 (Fig. 13a) and after 2015 (Fig. 13e). Thus, the slip after 2015 presumably occurred in the direction of releasing strain in areas, where the plate boundary interface was stickier than the surrounding area before 2011. Furthermore, the steady slip (V) decreased by v in B to E (E differs slightly due to the SSE). This reinforces the hypothesis that the steady velocity recovered after the 2011 Tohoku-Oki earthquake and returned to its original state in 2015. Areas C and D were stationary before 2011, revealed postearthquake afterslip, and continued with the new constant velocity slip added after 2015. In addition, B, C, and D presented few earthquakes in common ( Fig. 13c and f ). B is similar to C and D, except that it does not exhibit afterslip immediately after the main shock.
Next, we focus on the areas, where slip is smaller compared to the surroundings in Fig. 13e, such as east of B (Tokachi-Oki), between C and D (Miyagiken-Oki), and F (Aomoriken-Toho-Oki). In these areas, the afterslip of the 2011 Tohoku-Oki earthquake and the linear deformation since 2015 are smaller than the surrounding areas. The 2003 Tokachi-Oki earthquake (Yamanaka and Kikuchi 2003) lies to the east of B; the 2011 Tohoku-Oki earthquake (Fig. 13a, Ozawa et al. 2012) and its afterslip (Fig. 13b) lie between C and D; the 1994 Sanriku-Oki earthquake (Yamanaka and Kikuchi 2004) that has caused asperities is present at F. These locations are unlikely to move except during their unique earthquakes. Tobita (2016) listed six elements of the functional model's prediction limit as follows: (1) uncertainty of the steady velocity V, (2) limitation of function fitting, (3) status change of postseismic deformation mechanisms, (4) variable (inconstant) interseismic velocity, (5) future coseismic deformations, and (6) observation errors.

Error estimation and prediction performance
Based on the observed data (Fig. 6), the effects of (1), (2), and (4) are remarkably smaller, because they have not been detected explicitly in the last 10 years. The coseismic variation in (5) has an effect at stations close to the epicenter of an earthquake with a magnitude > 6, but does not reveal any marked effect at a certain distance from the epicenter. In the case of (6), despite an error due to environmental factors, such as GNSS signal shielding by trees at the observation station, the influence on the model can be minimized using stations with a good observation environment. Recently, the GEONET F5 solution (Muramatsu et al. 2021) has been published in place of the GEONET F3 solution, and a more accurate time series is now available. Although short-period variations such as daily variations have markedly improved, no significant difference is observed for variations over several months, as implemented in this analysis. In addition, for the reference station Fukue, the improvement is likely small even with the F5 solution, because the local variations of Fukue are considered to be an effect of the observation environment near the station. The two factors that deteriorate the accuracy of the functional model found in this study are widespread phenomena that occur in the middle of the time series (from 2015) and the phenomena caused by the fixed reference station. We did not consider seasonal deformations, such as annual variations, for the functional model and fixed station. Although the observed data include annual variations, as illustrated in Figs. 3 and 4, they are markedly small (a few millimeters) such that no correction for annual variations is included. As seasonal deformations are not similar every year (Gualandi et al. 2017), the differences in annual variations can affect the results in a short fitting period. For example, the difference between 2.0-and 3.9-year models could originate from these differences in the annual variations. Therefore, it makes sense to fit the model for each time series to a longer period to reduce the effect of seasonal deformations.
The largest error factor in the residuals includes the linear deformation vt term after 2015, which may represent the aforementioned category (3). This error may not occur due to a change in the mechanism itself; however, the beginning of a new slip phenomenon over a wider area that was entirely unpredicted. In addition, errors related to the fixed reference station are always associated with such observations, and a fundamental solution, for example, using the average of multiple stations as reference, is required.
Therefore, the evaluation of the prediction performance depends on the existence of the linear deformation v. The vt term itself is small; however, as it represents a linear deformation, it accumulates in one direction over time, resulting in an error of a maximum 10 cm in 10 years. Nevertheless, the model in Eq. (4), which considers v, can be modeled with an error of less than a few centimeters in 10 years. The future prediction depends on whether the assumption of Eq. (4) holds true; if a phenomenon, such as v occurs, then the error will increase with time. Therefore, the residuals from the time series at each station in time should be examined in the practical operation of the functional model, and this model should be adjusted as needed. No sudden large shift has been observed in the past 10 years; hence, even if the functional model remains unchanged, it will not pose a practical problem if we use the relative values of time change instead of the absolute values of the functional model (e.g., NS component of 2.0-and 3.9-year models in Fig. 2).

Use case of the functional model
Our functional model can be applied over a wide area at one time and can be computed using only the location and time to predict the future. For example, the uplift continues along the Pacific coast of the Tohoku region; thus, it is extremely challenging to decide the tide level to which the port facilities should be upgraded; however, if the amount and timing of the future uplift can be accurately predicted, the functional model could contribute to such future planning.
In this section, we present two examples of discerning minute crustal deformations. In general, to discern small crustal deformations from the time series, a simple and common approach is to set one fixed station as a reference in the vicinity according to the target location and period and further estimate the trend for each observation station to eliminate the fluctuations over a wider area.
Crustal deformations of the Boso Peninsula offshore SSE that occurred from the end of May to July 2018 (Ozawa et al. 2019) are depicted in Fig. 14 as an example of local variability extraction. In Fig. 14a, wherein a simple difference is considered between the two periods (difference between the 10-day averages of May 20-29 and June 26-July 5, 2018), the postseismic deformations of the 2011 Tohoku-Oki earthquake are scattered, and the pattern is nonuniform; however, in Fig. 14b, where the temporal variation is determined from the time series data after the predictions of the spatiotemporal model were removed, only the SSEs in the Boso Peninsula are clearly extracted, and the extraction power of crustal deformation is incomparable to that of the simple difference between the two periods.
As the widespread phenomenon of the postseismic deformation of a major earthquake appears uniformly within the local range, considering the simple difference between the two periods within the local range is typically sufficient. Therefore, the value of this model is in its ability to extract the phenomena over a wider area. The  Fig. 11a may be one of the examples occurring in a wider area. As an example of discerning minute changes over a wide area, the load deformation due to snow accumulation (Heki 2001) is illustrated in Fig. 15. The residuals of the functional model were calculated from the observed values on and around March 11 of each year after the earthquake (except for 2015), and the average values were drawn, because the deepest snow cover typically occurs in early March (Fig. 15). Creating such a figure for a slight change of 0.1-0.2 cm requires complicated operations, such as modeling the change at each station; however, in this case, the figure can be easily drawn by simply subtracting the functional model.
Considering the daily variations in the GNSS observations, we cannot merely express the accuracy of our functional model; however, if the average value of the GNSS observations is appropriately used, it can be implemented to detect signals of several millimeters over the entire model area; for example, it was the power of the functional model that led to the discovery of the existence of v after 2015. Thus, the model is a powerful tool with various applications for discerning small and gradual changes over a wide area.

Concluding remarks
The functional model for the postseismic deformation of large earthquakes by Tobita (2016) with statistical processing using a model with two logarithmic functions and one exponential function allows us to draw the following conclusions: (1) The functional model can predict the postseismic deformation of the 2011 Tohoku-Oki earthquake with good accuracy even after 10 years over a wide area exceeding 1000 km. (2) The same time constants of the time function are applicable for all observation stations and components; when these time constants change, the distribution of the spatial function for each time function varies at a certain rate for all observation stations, absorbing the change in the total functional model. Therefore, despite the differences in the original time constants for each observation station, the overall functional model can be well fitted to the time variation for each observation station.
(3) From the residuals of this functional model, a linear change began silently in 2015, without being accompanied by major events. The pattern occurred as a form in which a new change had been added as a linear deformation, and continued as of 2020, accumulating at a maximum of approximately 1 cm per year. This phenomenon was observed over the entire area used in this study; furthermore, when the phenomenon was simulated as a slip on the plate surface, in addition to the afterslip area north of Fukushimaken-Oki and its vicinity, a new slip appeared in the Urakawa-Oki area. The areas, where the new slip has sparsely occurred since 2015 coincided with the asperity areas that tend to cause unique earthquakes. (4) The functional model was modified by incorporating the new slip since 2015 in the form of a linear approximation, thereby remarkably improving the overall prediction performance of the functional model. Residuals of the functional model from the observed data were highly effective in discerning the occurrence of other events that were not postseismic deformations; hence, they can be applied in various fields. (5) The spatial distribution of the functional model for each time constant is similar to that of the afterslip and the viscoelastic relaxation by an independently computed physical model; that is, each time function relates to the time-dependent physical processes. (6) The comparison of the functional and physical models revealed that the viscoelastic relaxation near the land had a relatively short time constant, and that the fluctuations caused by the viscoelastic relaxation exhibited a wide range of time constants ranging from a few days to > 1000 days. (7) The evolution of the slip zone on the plate surface revealed that the 2011 Tohoku-Oki earthquake occurred on the shallow side of the pre-2011 sticking zone, whereas the afterslip occurred adjacent to the deep side of the main shock. Subsequently, a temporary slip was observed on the northeast side of the afterslip area and Urakawa-Oki in the first half of 2015, which was followed by assuming the occurrence of a new slip after 2015 that spread to the entire deeper side of the plate surface.

Abbreviations
GEONET: GNSS earth observation network system; GNSS: Global navigation satellite system; SSE: Slow-slip event; VLFE: Very low frequency earthquake.
Additional file 1: Fig. S1. Observed and predicted postseismic deformation using the functional model. Deformation is depicted for the four stations of a Minase, b Miyako, c Yamoto, and d Choshi (Fig. 1). The three components of east-west (EW), north-south (NS), and up-down (UD) are drawn for each station from the day after the earthquake to 2026.5. Black dots represent daily observations. Dashed lines represent the time series of predicted values for each fitting period (2.0, 3.9, 5.8, and 8.9 years). 0-year to the 3.9-year model. Light green line represents the average of the ratios. As the ratios become unstable as the value of the total length becomes smaller, they are sorted by the total length obtained for the 3.9-year model. Therefore, the average calculation is performed sequentially starting from one with the larger length for each component, and is terminated when the standard deviation of the average value becomes markedly larger (shown as the area defined by the light green line). b Same as a for the long-period logarithmic term in Eq.
(1). c Same as a for the exponential term in Eq. (1). Fig. S7. Checkerboard resolution test result of geodetic inversion used in this study. a Synthetic slip on the plate surface. b Inverted slip from the synthetic data of synthetic slip with Gaussian error of 1 cm for horizontal components and 5 cm for vertical components.