Characteristics of postseismic deformation following the 2003 Tokachi-oki earthquake and estimation of the viscoelastic structure in Hokkaido, northern Japan

Postseismic deformation of the 2003 Tokachi-oki earthquake (Mw 8.0) has been observed by GNSS. We analyzed the deformation observed in Hokkaido in the 2nd to the 7th year following the 2003 Tokachi-oki earthquake and examined the effect of two major mechanisms (i.e., afterslip and viscoelastic relaxation) for the observed postseismic deformation by fitting it with a model consisting of afterslip and viscoelastic relaxation. The thickness of the lithosphere, the viscosity of the asthenosphere, and the time decaying constant of afterslip were estimated to be 50 km, 2.0 × 1019 Pa s, and 0.110 year, respectively, which are concordant with those in the Tohoku region estimated in previous studies. The revealed characteristics of postseismic deformation are as follows. At most of the used stations, afterslip played the dominant role in the 2nd year and was still sustained near the coseismic area even in the 7th year. However, the calculated velocity due to viscoelastic relaxation was comparable to that due to afterslip at the stations in northern Hokkaido after the 5th year. Because the calculated velocity due to viscoelastic relaxation was landward near the coseismic slip area, afterslip near the coseismic slip area will be biased to be smaller if viscoelastic relaxation is ignored. A systematic spatial pattern of the residuals considering afterslip only highlights an importance for explaining the observation data. We also examined the effect of viscoelastic relaxation due to afterslip for the parameter estimation and found that it was too small to affect the estimated structure parameters.Graphical Abstract The observed and the calculated postseismic velocity in the 7th year following the 2003 Tokachi-oki earthquake: (left-hand panel) the vector indicates the observed postseismic velocity. (right-hand panel) The vector indicates the calculated velocity of afterslip (blue), viscoelastic relaxation (green), translation (pink) and the sum of them (red). The observed and the calculated postseismic velocity in the 7th year following the 2003 Tokachi-oki earthquake: (left-hand panel) the vector indicates the observed postseismic velocity. (right-hand panel) The vector indicates the calculated velocity of afterslip (blue), viscoelastic relaxation (green), translation (pink) and the sum of them (red).


Introduction
The 2003 Tokachi-oki earthquake (M w 8.0) was an interplate megathrust earthquake that occurred along the Kurile trench on September 26, 2003 (Japan Standard Time). The Japanese GNSS Earth Observation Network System (GEONET), consisting of ~1300 permanent GNSS stations, has detected postseismic deformation from many large earthquakes (e.g., Heki et al. 1997) and revealed that the 2003 Tokachi-oki earthquake was also accompanied by significant postseismic deformation (Ozawa et al. 2004). Postseismic deformation is mainly caused by two mechanisms: afterslip and viscoelastic relaxation (e.g., Scholz 2002). Afterslip is aseismic slip in and around the coseismic slip area in the postseismic period. Estimation of the spatiotemporal distribution of afterslip is important for understanding the healing process of fault strength and the transition from the postseismic period to the interseismic period. On the other hand, viscoelastic relaxation is the time-dependent deformation caused by the flow of the viscoelastic asthenosphere in response to a coseismic stress change (e.g., Thatcher et al. 1980). When an interplate coupling rate in the postseismic period is estimated, viscoelastic relaxation should be corrected or be properly modeled because ignoring its effect can cause a systematic error in the estimation (e.g., Freed et al. 2006).
It is assumed that afterslip plays the dominant role in postseismic deformation and that viscoelastic relaxation has only a secondary role for a few years after the mainshock (e.g., Scholz 2002). However, viscoelastic relaxation cannot be ignored after several years in cases of large earthquakes like the 2003 Tokachi-oki earthquake. Although in a number of other studies (e.g., Miyazaki et al. 2004;Ozawa et al. 2004;Baba et al. 2006) postseismic deformation due to the 2003 Tokachi-oki earthquake was analyzed, all of them considered only the effect of afterslip and ignored the effect of viscoelastic relaxation. Tanaka (2007) analyzed postseismic deformation for the first 2 years and showed that the signal of viscoelastic relaxation could be significantly included in the observed displacement after several months. However, Tanaka (2007) checked similarity between the observed velocity and the calculated one due to viscoelastic relaxation and did not estimate the parameters of viscoelastic structure to fit the observed data. No previous studies have used the data from the 3rd year after the mainshock to estimate parameters of the viscoelastic structure. Although postseismic deformation due to the 1993 Hokkaido Nansei-oki earthquake (M w 7.7) was modeled with viscoelastic relaxation and parameters of viscoelastic structure were estimated by Ueda et al. (2003), they pointed out that the estimated parameters of structure may not explain the viscoelastic behavior of the fore-arc region of Hokkaido enough because the 1993 event occurred in the back-arc region.
Estimation of viscoelastic structure and modeling of postseismic deformation are also important in terms of understanding the earthquake cycle (e.g., Fukahata et al. 2004;Ito and Hashimoto 2004), and they enable us to discuss the earthquake cycle in the Tokachi-oki region where large earthquakes occurred in 1843, 1952(e.g., Abe 1999Yamanaka and Kikuchi 2003) using modern geodetic data.
In this paper, we assumed both viscoelastic relaxation and afterslip as mechanisms of postseismic deformation following the 2003 Tokachi-oki earthquake, and estimated viscoelastic parameters and a time constant of afterslip by using the data from the 2nd to the 7th year after the mainshock (i.e., September 26, 2004-September 26, 2010.

Data analysis
We used the GNSS daily coordinates (so-called the GEONET F3 solution) provided by the Geospatial Information Authority of Japan (GSI) (Nakagawa et al. 2009) at 81 stations in Hokkaido (Fig. 1a). Because the effect of postseismic relaxation following the 1993 Hokkaido Nansei-oki earthquake still continued in 2001 (Ueda et al. 2003;Suwa et al. 2006), and because the eruption of an active volcano (Mt. Usu) occurred in March 31, 2000, we did not use the data in western Hokkaido. As the reference station, we used Station 0151 (Kanita) and considered movement relative to this station. First, we estimated a secular interseismic velocity (Fig. 1b) before the mainshock of the 2003 event and removed it from the observed postseismic displacement. We used 4.5-year-long data from March 1, 1999, to September 1, 2003, to estimate the interseismic velocity because this period was quiet in terms of crustal activity, and long enough to ignore the effect of the annual component, to estimate the linear trend (Blewitt and Lavallée 2002). After the 2003 event, some large earthquakes occurred near Hokkaido and their coseismic displacements were observed by GEONET (Kimura and Miyahara 2013). Therefore, we eliminated their displacements from the observed data. The coseismic displacements of the 2004 Kushiro-oki earthquakes (M w 7.0, M w 6.7) and the 2006 Central Kurile earthquake (M w 7.8) were calculated using the models estimated by Nishimura (2009) and Takahashi and Kasahara (2007), respectively. We eliminated the coseismic displacements of the 2007, 2008, and 2009 Tokachi-oki earthquakes (M w 6.0, M w 6.8, and M w 6.4, respectively) by correcting the offsets in time series. We obtained the amplitude of the offset as the difference between coordinates averaged for 10 days just before and after the earthquakes. We finally modeled the corrected time series as a sum of a piecewise linear term, an annual and a semiannual variation, and an observation error. The postseismic velocity from the 2nd to the 7th year after the mainshock (i.e., September 26, 2004-September 26, 2010 was derived from the piecewise linear term. Section lengths for the piecewise linear terms were 1 and 2 years for the horizontal and the vertical components, respectively. We did not use the data in the first year to estimate the time constant of afterslip and the viscoelastic parameters because velocity and deceleration of immediate postseismic deformation seem to be much faster than the long-term behavior ) and cannot be explained by our simple assumption for the modeling of the postseismic deformation. Coseismic and postseismic deformation of the 2011 Tohokuoki earthquake (M w 9.0) was observed at most stations in Hokkaido. Because it is difficult to distinguish between the transient deformation due to the 2003 and the 2011 events, we used the data until the 7th year from the 2003 Tokachi-oki earthquake. Although the data in the first year were not used to estimate parameters of viscoelastic relaxation, we used the data for 3 months immediately after the earthquake to construct an afterslip model because it is considered to characterize the spatial distribution of afterslip as explained in "Deformation modeling" section.
The main characteristics of the postseismic velocity are as follows. Through the analyzed period, trenchward movement was observed at all the used stations except those located in the easternmost area. Because we subtracted the linear component of the GNSS time series from the data as the steady deformation, trenchward movement indicates that the steady deformation rate before the 2003 Tokachi-oki earthquake had not yet been recovered through the analyzed periods. In the easternmost area, trenchward movement was significant in the 2nd to the 4th years, which could be mainly due to postseismic deformation caused by the 2004 Kushiro-oki earthquakes ( Fig. 2a; Additional file 1: Figure S1a, b). In the 4th year, all the stations moved east by comparison to the stations in other periods (Additional file 1: Figure  S1b), which suggests that postseismic deformation due to the 2006 and 2007 Central Kurile earthquakes was significant. The vertical velocity showed significant uplift along the Pacific coast ( Fig. 2c; Additional file 1: Figure S2a, b), which suggests that afterslip occurred in the down-dip extension of the coseismic slip region.

Deformation modeling
We constructed a postseismic deformation model, which consists of both afterslip and viscoelastic relaxation. We assumed that the surface displacement U(x, t) due to afterslip evolves logarithmically with time, which is represented by the following equation (Marone et al. 1991), where t is years following the mainshock, x is a station, A is a coefficient of amplitude, B is a decaying time constant, and u(x) is surface displacement at the station x calculated from an afterslip model during the first 3 months, which we also estimated using the inversion method of Nishimura (2009) as explained in "Inversion procedure to estimate the afterslip distribution for the first 3 months" section. We estimated the afterslip distribution for the first 2, 3, 4, and 5 months and preferred the displacement calculated from the 3-month model for u(x) (Fig. 3a). This is because the spatial pattern of the calculated surface displacement was the most similar to that of the observed displacement after the second year among the models for 4 periods, and the effect of viscoelastic relaxation was assumed to be small. We checked the effect of viscoelastic relaxation by model calculation using the viscoelastic parameters in the Tohoku region, south of Hokkaido, estimated by Diao et al. (2014) and found that it (<0.3 cm) was small enough to be ignored in the observed postseismic deformation which was up to ~10.0 cm during the first 3 months. It should be noted that we assumed that the spatial distribution of (1), we obtained the following relationship between A and B.
For the modeling of viscoelastic relaxation, we assumed a two-layered structure that consists of an elastic layer overlying a Maxwell viscoelastic half-space. The elastic layer and the viscoelastic half-space correspond to lithosphere and asthenosphere, respectively. The thickness of the elastic layer (H) and the viscosity of the viscoelastic half-space (η) are free parameters to be estimated. We assumed elastic parameters based on Katsumata et al. (2006) and PREM (Dziewonski and Anderson 1981) (Additional file 2: Table S1). We used the coseismic slip model estimated by the Geospatial Information Authority of Japan (2003) and Nishimura (2009) for the 2003 Tokachi-oki and the 2004 Kushiro-oki earthquakes, respectively. We also considered viscoelastic relaxation caused by the 2004 Kushiro-oki earthquakes because coseismic displacement was observed and the stations in eastern Hokkaido suddenly accelerated after these events. Postseismic displacement due to viscoelastic relaxation was calculated using the FORTRAN code named PSGRN/PSCMP (Wang et al. 2006).

Parameter estimation method
We where N is the number of stations, M is the number of periods and components, v obs ij is the observed j-th velocity at the i-th station, v cal ij is the model calculated j-th velocity at the i-th station consisting of afterslip due to the 2003 event, viscoelastic relaxation due to the 2003 and the 2004 events and translation parameters, and σ ij is the GNSS velocity error. σ ij is calculated from the equation proposed by Mao et al. (1999) and Dixon et al. (2000) and includes the effect of white noise and colored (flicker and random-walk) noise. It is important to consider the effect of colored (i.e., time-correlated) noise because Mao et al. (1999) reported that velocity error of the GNSS time series will be underestimated by factors of 5-11 if its effect is ignored. The obtained error σ ij of the horizontal ij and the vertical velocity is about 2.5 and 3.7 mm/year, respectively.

Inversion procedure to estimate the afterslip distribution for the first 3 months
To obtain the term u(x) in the afterslip model expressed as Eq.
(1), we estimated the 3-month afterslip distribution. To obtain the observed displacement for the first 3 months, we didn't use the piecewise linear approximation because the velocity of postseismic deformation decays rapidly and a linear approximation tends to ignore early postseismic displacement and underestimates total one. Instead, we obtained the observed displacement in the first 3 months as the difference of the coordinates between periods immediately and 3 months after the mainshock. When we calculated the displacement, an annual and a semiannual variation estimated in "Data analysis" section was removed from the data. The site coordinate for immediately after the mainshock was calculated by averaging the coordinates for the first 11 days. The coordinate after 3 months was calculated by averaging the coordinates for 21 days centering the date just 3 months after the earthquake.
To estimate the afterslip distribution for the first 3 months, we used the inversion method of linear least squares (Nishimura 2009). Synthetic surface displacement was calculated by the formulation of Okada (1985), assuming a homogeneous elastic half-space. Because the observed vertical component of GNSS data has larger measurements errors, in this case which were about 4 mm and 9 mm for the horizontal and the vertical component, respectively, it was difficult to obtain a model that reproduced both horizontal and vertical displacements simultaneously if the data were weighted just inversely to the square of the observation error. Because the models of previous studies (e.g., Ozawa et al. 2004;Baba et al. 2006) didn't explain the observed vertical component, which means that slip below the land has not been well resolved and discussed, we used the vertical error reduced by multiplying the ratio of the averaged horizontal error to the vertical one and the vertical component was weighted as twice the horizontal component to discuss slip below the land. The observed displacement was inverted to the slip on 240 small rectangular faults, which were placed on the plate interface. The depth of the plate interface followed the model of Miyazaki et al. (2004) (Fig. 1b), which combined the models of the Earthquake Research Committee (2003) and Katsumata et al. (2003). We did not constrain the slip direction, but we constrained the slip outside the modeled region to be zero.
As a priori information, we applied a constraint on the smoothness of the slip distribution, which is indicated by the following term, where s n,m is the amount of slip on the n-th and m-th subfaults along the strike and dip direction. We calculated this term for strike slip and dip slip, respectively. The slip distribution, which minimizes the sum of squared data misfit and the sum of the smoothness terms with a hyperparameter, is the optimum solution. To determine the optimum contribution of a priori information to the solution, the hyper-parameter was determined by minimizing the Akaike's Bayesian Information Criterion (ABIC) (Akaike 1980).
We conducted the checkerboard test to check the resolving power of the used geodetic data for slip on the plate interface. We calculated synthetic surface displacement with the assigned checkerboard-like slip distribution (Additional file 1: Figure S3a, b) and added the synthetic noise with the standard deviation for the 3-month data to the synthetic displacement. The observation error and weighting between the horizontal and vertical components are also same as the inversion of the 3-month data. Slip of 30 cm in the direction of N115°E, opposite to the relative plate motion, was assumed with a ~100-km-wide checkerboard pattern. The value of the assumed slip was determined from the approximate value of the peak of the estimated afterslip (Fig. 3a). The direction of relative plate motion was calculated from a global (4) n m (2s n,m − s n−1,m − s n+1,m ) + (2s n,m − s n,m−1 − s n,m+1 ) 2 plate motion model (e.g., DeMets et al. 1994;Sella et al. 2002).
From the results of resolution test (Additional file 1: Figure S3c, d), slip at the deeper area of the mainshock is considered to be well resolved. This result suggests that slip below the observation stations is well resolved by weighting of the vertical displacement. However, in the most northwestern area of the model region (Additional file 1: Figure S3d) small reverse slip was estimated although no slip was assumed there (Additional file 1: Figure S3b). The results also suggest that slip in the shallower and near the trench region is poorly resolved. Because the input slip has 100 km square or 100 km × 120 km distribution, therefore the recovered slip (Additional file 1: Figure S3c, d) indicates that the slip area larger than 100 or 120 km can be resolved.

Estimated spatial distribution of afterslip
The estimated afterslip distribution explained the observed displacement (Fig. 3b, c), and reproduced surface displacement was used as u(x) in Eq. (1). The estimated afterslip distribution shows that the peak of the slip is located on the up-dip side of the coseismic slip area (Fig. 3a). This characteristic is almost consistent with the estimated afterslip distribution for approximately the first 6 months (Ozawa et al. 2004) and the first year (Baba et al. 2006). A distinct difference in the afterslip distribution between the previous studies and ours is found in the coseismic slip area and its down-dip extension, particularly. The estimated slip in the down-dip extension is necessary to reproduce the observed uplift in the Tokachi Plain, which was not explained by the previous studies, but it is not considered to be totally artificial one according to the resolution test (Additional file 1: Figure S3). In addition to that, although the previous studies estimated little slip in the coseismic slip area surrounded by large slip in the updip and down-dip extension, considerable afterslip in the coseismic slip area was estimated in our model, which is considered to be due to poor resolving power of slip below the offshore region. Indeed, the relative contribution of a penalty function for smoothing [Eq. (4)] to the sum of squared data misfit, which was determined by ABIC, was very large. Therefore, characteristics of our afterslip distribution are attributed to both weighting of datasets and a priori smoothing constraint. An error of the estimated slip was smaller than the estimated slip in the most part of the model region, which supports our inference of the slip characteristics (Additional file 1: Figure S4). Assuming a shear modulus of 44.1 GPa, the seismic moment (M 0 ) of the estimated afterslip would be ~5.36 × 10 20 N m, which is equivalent to M w ~ 7.75.

Estimated parameters for afterslip time constant and viscoelastic structure
The estimated parameters were B = 0.110 year, H = 50 km, and η = 2.0 × 10 19 Pa s. The estimated viscoelastic parameters are concordant with those estimated by previous studies in the Tohoku region (Suito and Hirahara 1999;Diao et al. 2014) where the same Pacific plate subducts. Ueda et al. (2003) proposed that viscosity in the fore-arc region was probably larger than that in the back-arc region based on the studies of the thermal and the seismic velocity structure of northeast Japan. Our estimated viscosity is consistent with their proposal because ours is 5 times larger than their estimated viscosity (4.0 × 10 18 Pa s). The cumulative moment released by afterslip, which is calculated from the estimated value of B, is equivalent to M w ~ 7.95 for the first year and M w ~ 8.12 for 7 years.
The reduced Chi-squared ( Fig. 4; Additional file 1: Figure S5) was obtained by dividing the value of χ 2 of Eq. (3) by the degrees of freedom (711 in this case). The acceptable parameters which are determined by a Chi-squared test with 68.3 % confidence limit are also indicated. From the distribution of the acceptable parameters, it is suggested that our data couldn't constrain an elastic thickness H. The upper bound of the confidence limit for viscosity is larger than 1.0 × 10 20 Pa s, which implies that an afterslip model without viscoelastic relaxation can explain the data within the observation error. However, the estimated viscosity of asthenosphere is also concordant with the value calculated with the experimental constitutive law of olivine (Muto et al. 2013). Therefore, the estimated parameters are valid and the calculated velocity due to viscoelastic relaxation ( Fig. 5; Additional file 1: Figures S6, S7) can be included in the observed velocity.
To check whether viscoelastic relaxation is an ignorable process in the observed postseismic deformation or not, we conducted an additional grid search, in which the search range was from 0.03 to 0.21 year for B and viscoelastic parameters were fixed to H = 50 km, and η = 1.0 × 10 22 Pa s. This model can be regarded as a model of only afterslip because deformation due to viscoelastic relaxation is expected to be small enough to be ignored. The range for B was derived from that of the acceptable parameters in η = 1.0 × 10 20 Pa s (Additional file 1: Figure S5b). The model with B = 0.09 year was obtained as the optimum, and this combination of the model parameters is also acceptable with respect to the best-fitted viscoelastic model according to a Chi-squared test. From the estimated value of B, the cumulative afterslip moment is decreased by about 10 % if viscoelastic relaxation is ignored. It is notable that systematic landward residuals near the coseismic slip and around (142°E, 43.5°N) remained in the horizontal component before the 4th year (Additional file 1: Figure S8a-c). These residuals are similar to the velocity predicted from viscoelastic relaxation ( Fig. 5d; Additional file 1: Figure S6e, f ) and support our assumption that both afterslip and viscoelastic relaxation affects the observed velocity. Landward residuals near the coseismic slip diminished after the 5th year (Additional file 1: Figure S8d-f ) and are probably explained by postseismic deformation due to the 2007,2008,2009 Tokachi-oki earthquakes. We, therefore, suggest that viscoelastic relaxation is necessary to explain the observed postseismic deformation.

Effect of viscoelastic relaxation due to afterslip for the parameter estimation
The calculated seismic moment of afterslip with the estimated parameter B, which is equivalent to M w 7.95, is comparable to that of the mainshock (corresponding to M w 8.0) for the first year. Because it suggests that afterslip can cause significant viscoelastic relaxation and that ignoring it might bias the estimated parameters, we examined how large viscoelastic relaxation due to afterslip affected parameter estimation in the modeling of postseismic deformation. Estimation procedure was the same as the estimation without afterslip viscoelastic relaxation, and only difference was to include velocity due to afterslip viscoelastic relaxation in v cal ij . Calculation procedure of velocity due to afterslip viscoelastic relaxation is summarized in Additional file 3.
The result shows that the optimum viscoelastic parameter was the same in the models both with and without afterslip viscoelastic relaxation. The optimum afterslip  time constant B (0.120 year) is slightly larger than that without afterslip viscoelastic relaxation but located in the confidence limit. Including afterslip viscoelastic relaxation didn't reduce the value of Chi-squared, which suggests that it isn't required by the observed data. Therefore, we can conclude that viscoelastic relaxation due to afterslip doesn't affect the parameter estimation significantly. Horizontal velocity due to afterslip viscoelastic relaxation (Additional file 1: Figure S9) is landward along the coast, and its amplitude, which is smaller than 5 mm/year even in the 7th year, is a few times smaller than the velocity due to viscoelastic relaxation caused by the mainshock.

Characteristics on observed and calculated surface displacement
Although the estimated translation component (Additional file 2: Table S2) can be interpreted as instability of the reference station or large-scale deformation affecting the whole observation network, the horizontal components in the 2nd and the 4th years are larger than those in other periods. That in the 4th year suggests that it mainly represents the postseismic deformation due to the 2006 and 2007 Central Kurile earthquakes. That in the 2nd year might include afterslip due to the 2004 Kushiro-oki earthquakes.
Although our model reproduced the horizontal velocity well ( Fig. 2; Additional file 1: Figure S1), systematic residuals were distributed in some areas. Some of them can be interpreted as the signals of afterslip due to the 2004 Kushiro-oki earthquakes and postseismic deformation due to the 2008 and 2009 Tokachi-oki earthquakes, which were not considered in our model. In the western area of the studied region (around 43°N, 142°E), systematic eastward residuals were distributed throughout the analyzed periods, which suggest that this is partially due to the postseismic deformation of the 1993 Hokkaido Nansei-oki earthquake. For the vertical velocity, our model reproduced the observed uplift along the Pacific coastal region qualitatively but not quantitatively ( Fig. 2; Additional file 1: Figure S2). Systematic residuals larger than the velocity error were distributed in the Tokachi Plain. These uplift residuals suggest that an afterslip distribution since the 2nd year changed from the 3-month distribution and that the afterslip magnitude below the land is larger than that predicted from our logarithmic function model throughout the analyzed periods. Transient uplift after the seventeenth-century Kurile tsunamigenic mega-earthquake, which is suggested to reach almost M9 (Satake et al. 2008;Ioki and Tanioka 2016), was reported by previous studies (e.g., Sawai et al. 2004) and has been considered to be caused by aseismic creep in the downdip extension on the plate interface. Our results also suggest that large afterslip, which occurred in the down-dip region of the coseismic slip region, is a plausible explanation for transient uplift in the postseismic period. These residuals suggest that it is difficult to predict the spatiotemporal evolution of afterslip by our simple model using the early postseismic data (the first 3 months in this case). We believe that our simple assumption that spatial distribution of afterslip in the first 3 months doesn't change is acceptable at least until the first year by comparison with the afterslip distribution estimated by previous studies (e.g., Miyazaki et al. 2004;Ozawa et al. 2004;Baba et al. 2006). However, deviation from our assumption was revealed in the later periods. A model in which spatial distribution of afterslip is evolved with time is required in a further study.
We calculated each contribution of afterslip and viscoelastic relaxation in the modeled postseismic deformation ( Fig. 5; Additional file 1: Figures S6, S7). Because the velocity calculated from viscoelastic relaxation due to the 2004 Kushiro-oki earthquakes was much smaller than that due to the 2003 Tokachi-oki earthquake, it can be ignored. For the horizontal velocity, afterslip played the dominant role until the 4th year except in easternmost Hokkaido and was still sustained near the coseismic slip area even after 7 years ( Fig. 5; Additional file 1: Figure S6). On the other hand, the velocity of viscoelastic relaxation was almost equal to that of afterslip from the 5th year in northern Hokkaido. For the vertical velocity, dominant process in the postseismic deformation was dependent on regions in the 2nd-3rd years; however, after the 4th year, velocity of viscoelastic relaxation was almost equal to that of afterslip in the all regions ( Fig. 5; Additional file 1: Figure S7). A notable characteristic of the calculated velocity of viscoelastic relaxation is that some stations located near the coseismic slip moved landward, which is opposite to the trenchward movement of the calculated velocity of afterslip and the observed velocity. This is theoretically explained by the model calculation, that is, in the coseismic slip area, horizontal displacement of viscoelastic relaxation has a direction opposite to that of the coseismic horizontal displacement (Rundle 1978). In the case of the 2011 Tohoku-oki earthquake, the coseismic slip occurred far from the onshore region (e.g., Ozawa et al. 2011). Because movement toward the trench of the onshore region was predicted by both afterslip and viscoelastic relaxation in the case of the 2011 event (Yamagiwa et al. 2015), it is difficult to distinguish between the contributions of both mechanisms in postseismic deformation. However, in the case of the 2003 Tokachi-oki earthquake, we can conclude that afterslip still continued near the coseismic slip area until at least the 7th year, and its magnitude will be biased to be smaller if we ignore viscoelastic relaxation as we showed in "Estimated parameters for afterslip time constant and viscoelastic structure" section.
The assumed two-layer viscoelastic structure doesn't include the characteristic structure in the subduction zone, namely the subducting slab and the mantle wedge. In previous studies, both layered model (e.g., Ueda et al. 2003;Pollitz et al. 2006;Hoechner et al. 2011;Diao et al. 2014;Yamagiwa et al. 2015) and two-or three-dimensional structure models (e.g., Hu and Wang 2012;Sun et al. 2014) have been used to explain the observed postseismic deformation in the subduction zone around the world and difference between two structure models has been discussed. Here, we discussed the effect of including characteristic structure of subduction zone using models with and without the high-viscosity slab in the theoretical study (Pollitz et al. 2008). Both viscoelastic structure models with and without the slab predict landward motion above the source fault and trenchward motion in a region far from the source fault. The model including the slab predicts larger landward displacement near the coseismic fault and shift of the boundary between landward and trenchward motion to the inland region. The pattern of vertical displacement in the hanging wall is not affected by including slab, that is, subsidence above the down-dip edge of the source fault, uplift in a region with a medium distance from the source, and wide spread subsidence in a far field ( Fig. 5; Additional file 1: Figure S7). However, viscoelastic displacement is reduced in the hanging wall except for the regions above the coseismic fault by introducing the slab with the same asthenospheric viscosity. Therefore, we consider that ignoring a three-dimensional structure did not change our conclusion for an importance of viscoelastic relaxation but bias our estimated viscoelastic parameters. Our estimated viscosity can change to be smaller if high-viscosity slab is introduced. It means that our estimated viscosity is probably larger than that of "real" asthenosphere.

Earthquake hazard inferred from the estimated afterslip model
Interseismic coupling in the Tokachi-oki region was estimated by several previous studies (e.g., Ito et al. 2000;Suwa et al. 2006;Hashimoto et al. 2012), and they suggested high coupling on the plate interface. Comparison between interplate coupling, coseismic, and postseismic slip of the 2003 Tokachi-oki earthquake is useful to discuss a potential of future earthquakes along the Kurile trench. It was reported that coseismic slip released the total accumulated strain due to the interplate coupling if the coupling rate has been constant through the interseismic period since the 1952 event (Miura et al. 2004). Although calculated afterslip of our model for 7 years is ~1 m for the up-dip region of the coseismic slip, which is much smaller than the coseismic slip, high coupling was also found in this region in the interseismic period. Therefore, much of accumulated strain in this region was not released by afterslip. Because coseismic slip of the huge tsunamigenic earthquake in the seventeenth century was estimated in the shallower part of the plate interface that includes the up-dip region of the 2003 event (Satake et al. 2008;Ioki and Tanioka 2016), the accumulated strain due to high interplate coupling there suggests potential of a large tsunamigenic earthquake in the future.
Interseismic coupling was estimated in the east adjacent region (i.e., the Nemuro-oki region) in the previous studies, and its rate seems to be large. Although little afterslip due to the 2003 Tokachi-oki earthquake was estimated in that region, afterslip due to the 2004 Kushiro-oki earthquakes is suggested to occur but might end in the 3rd year ( Fig. 1; Additional file 1: Figure S1e). Therefore, we speculate that accumulated strain in this region could be partially released but a most part of accumulated may still remain. It should be noted that the shallower part of this region also overlapped with the source area of the seventeenth-century giant earthquake.
For the down-dip region of the coseismic slip, especially below the land, there are serious differences of interseismic coupling distribution among the models, that is, the weaker coupling (Suwa et al. 2006) or no coupling (Ito et al. 2000;Hashimoto et al. 2012) was estimated. If we assume significant interseismic coupling there, accumulated strain was partially released by our estimated afterslip that was up to ~0.5 m for 7 years but may still remain. However, plate interface below the land is located deeper than 60 km and its temperature is higher than 300 °C according to the thermal structure model in northeast Japan (e.g., Wada et al. 2015). Accumulated coupling below the land may not lead to the nucleation of earthquake because transition of characteristics on fault friction from velocity weakening to velocity strengthening occurred at about 300 °C and earthquake is not expected to nucleate below a depth of that degree (e.g., Scholz 1998). Aseismic slip may be a key release process of accumulated coupling as inferred in the case of the seventeenth-century Kurile earthquake (e.g., Sawai et al. 2004) if the coupling below the land exists significantly.

Conclusions
We analyzed postseismic deformation of the 2003 Tokachi-oki earthquake observed by GEONET from the 2nd to the 7th year following the mainshock and modeled it with both afterslip and viscoelastic relaxation. We estimated the thickness of lithosphere and the viscosity of the asthenosphere (50 km and 2.0 × 10 19 Pa s, respectively) and found our values concordant with those reported for the Tohoku region, northern Japan. We suggest that viscoelastic relaxation is necessary for explaining the observed postseismic deformation by constructing the only afterslip model and fitting the observation data. Although the horizontal velocity of the calculated postseismic deformation agrees well with the observed one, the vertical velocity was not reproduced quantitatively. This suggests that, over 7 years, the spatial distribution of afterslip changed from the 3-month distribution and that velocity of afterslip that occurred in the down-dip area of the coseismic slip is faster than that predicted from our model. This makes it difficult to predict the spatiotemporal evolution of afterslip using only the early data and the assumptions of the simple model used here. Afterslip played the dominant role in postseismic deformation until the 4th year at most stations and was sustained near the coseismic source area in the 7th year. However, the velocity of viscoelastic relaxation was almost equal to that of afterslip after the 5th year in northern Hokkaido. Because velocity due to viscoelastic relaxation near the coseismic slip was landward, which is opposite to that of afterslip, a few or even several years of afterslip near the coseismic slip will be biased to be smaller if the effect of viscoelastic relaxation is ignored. Although we constructed another model, which additionally includes afterslip viscoelastic relaxation, to examine the effect for the parameter estimation, the estimated viscoelastic parameters didn't change. Velocity due to afterslip viscoelastic relaxation is relatively smaller than that due to other deformation sources. Based on our estimated afterslip distribution, coseismic slip and the interplate coupling distribution of previous studies, we suggest that afterslip following the 2003 Tokachi-oki earthquake released accumulated strain due to interplate coupling in the source region of the seventeenth-century tsunamigenic giant earthquake partially but much of it still remains.