Improvement on spatial resolution of a coseismic slip distribution using postseismic geodetic data through a viscoelastic inversion

Obvious crustal deformation is observed during a postseismic period as well as a coseismic period associated with a large earthquake. Major mechanisms of transient postseismic deformation are known as afterslip and viscoelastic relaxation. Since the viscoelastic relaxation occurs as a response to a coseismic slip, postseismic deformation provides information on coseismic deformation through the viscoelastic response. However, most previous studies have not thoroughly utilized postseismic geodetic observational data for revealing coseismic slip behaviors. In this study, we developed a slip inversion method that simultaneously estimates coseismic slip and postseismic slip distributions from coseismic and postseismic geodetic observational data using viscoelastic Green’s function (viscoelastic inversion method). We investigated the performance of the viscoelastic inversion method via two synthetic tests: one assumed a strike–slip event along an inland fault, while the other assumed a dip–slip event along a plate interface in a subduction zone. Both synthetic tests demonstrated that when extensive postseismic observational data were given, the viscoelastic inversion method provided a superior spatial resolution of coseismic slip distributions compared to conventional elastic inversion distributions. We also applied the viscoelastic inversion method to co- and post-seismic deformations associated with the 2011 Tohoku-oki earthquake. The seafloor geodetic observational network of the off-Tohoku region has been widely extended after the occurrence of the mainshock. Using this extended seafloor geodetic observational data, we successfully improved the spatial resolution of the coseismic slip distribution through the viscoelastic inversion method. Furthermore, using the seafloor observational data during the postseismic period, our inversion method enables us to obtain high spatial resolution of the coseismic slip in the offshore area and a reasonable coseismic slip distribution even if seafloor observational data during the coseismic period are unavailable. These results clarify the importance of deploying a geodetic observational network even after large coseismic events to assess past coseismic slip behaviors by considering the viscoelasticity of the Earth.


Introduction
Associated with the occurrence of large earthquakes, characteristic crustal deformation is measured during both co-and post-seismic periods. Coseismic deformation results from a coseismic fault slip, whereas postseismic deformation is governed by multiple physical mechanisms: afterslip, viscoelastic relaxation, and fault locking (e.g., Wang et al. 2012). These postseismic mechanisms provide important information on the frictional properties of faults (e.g., Miyazaki et al. 2004) and on an underground rheological structure. Therefore, geodetic observational networks have been often urgently extended after large coseismic events to precisely detect postseismic deformations; examples include an extension of the offshore GNSS-Acoustic observational network after the 2011 Tohoku-oki earthquake, Japan (Kido et al. 2015;Tomita et al. 2017) and an extension of the land-based GNSS network close to inland faults of the 2007 Chuetsu-oki earthquake, Niigata, Japan (e.g., Ohta et al. 2008;Iinuma et al. 2008). Postseismic displacements inherently contain information about coseismic slips via a viscoelastic Earth because viscoelastic relaxations are caused by stress perturbations associated with coseismic ruptures. Therefore, we can constrain a coseismic slip distribution by installing or extending a geodetic network after the occurrence of an earthquake, even if the geodetic observational network is spatially insufficient during the coseismic period.
Most previous studies that modeled coseismic slip distributions for large earthquakes based on observations have utilized only coseismic observational data assuming an elastic Earth. Additionally, most studies subtracted viscoelastic responses to coseismic slip models estimated only using coseismic observational data from original postseismic geodetic data, and then modeled postseismic slip distributions from the subtracted postseismic geodetic data assuming an elastic Earth (e.g., Diao et al. 2013;Lubis et al. 2013;Iinuma et al. 2016;Freed et al. 2017). Herein, we refer to this inversion method assuming an elastic Earth as the "elastic inversion" method; the elastic inversion method discards the information on coseismic slips contained in the postseismic geodetic data. Meanwhile, some recent studies have simultaneously estimated contributions of afterslip and viscoelastic relaxation using a stress kernel, which is independent of a coseismic slip model (e.g., Tsang et al. 2016;Moore et al. 2017;Qiu et al. 2018). This approach is quite useful and reasonable, but it also cannot model coseismic behaviors by utilizing postseismic deformation data. Additional previous studies have attempted to simultaneously model co-and post-seismic deformation considering a viscoelastic Earth by introducing viscoelastic Green's function (herein, we call this inversion method as the "viscoelastic inversion" method) (e.g., Ito and Hashimoto 2004;Hoechner et al. 2011;Yamagiwa et al. 2015;Ito et al. 2018). Slip inversions using the viscoelastic Green's function have been often performed for modeling interseismic deformations (e.g., Fukahata et al. 2004;Hashimoto et al. 2009;Li et al. 2015), whereas the viscoelastic inversion method, estimating both the co-and post-seismic slip distributions, has not been popularly performed. Furthermore, improvements on the spatial resolution of coseismic slip distributions stemming from considering postseismic geodetic observation data have not been thoroughly investigated using the viscoelastic inversion method. Since the estimation of fault-slip distributions with better spatial resolution is crucial to investigate the process of strain accumulation and release and to evaluate the temporal change in the seismic potential along the faults during a seismic cycle, it is important to assess potential performance of the viscoelastic inversion method.
In this study, we develop a viscoelastic inversion method that simultaneously inverts co-and post-seismic geodetic observational data. We apply this method to synthetic observational data and investigate its performance, especially for an improvement on the spatial resolution of a coseismic slip distribution by introducing postseismic observation data. Then, we also apply our method to actual geodetic observation data associated with the 2011 Tohoku-oki earthquake.

Observation equation
We introduce an observation equation that links coand post-seismic displacements with co-and post-seismic slip distributions via viscoelastic Green's function. Although our approach is similar to Yamagiwa et al. (2015) and Ito and Hashimoto (2004), our method has a few notable differences; these include optimization procedures of the viscosity, the thickness of an elastic layer, and relative weight between co-and post-seismic observation data, as pointed out in this section. Viscoelastic responses due to a fault slip can be modeled by assuming linear fluid rheology. Besides, although the poroelastic rebound is known as one of the postseismic relaxation processes (e.g., Lubis et al. 2013;Hu et al. 2014;Gunawan et al. 2019), it is ignored in this study. This is because its contributions to the long-term postseismic deformation are relatively small compared with other processes (e.g., Lubis et al. 2013) and because the main aim of this study is the assessment of the utility of the viscoelastic inversion. However, note that poroelastic rebound can be also involved by calculating poroelastic Green functions as same with viscoelastic relaxation, because the poroelastic rebound is also linked to the coseismic stress changes.
The viscoelastic surface displacement u i for ith component caused by a slip motion a j for the jth direction on a plate interface Σ can be written with the following equation (Fukahata et al. 2004): where G ij (x, t; ξ, τ ) is the viscoelastic Green's function and indicates viscoelastic displacement at a point x at a time t on the surface, due to a unit step slip at point ξ at time τ on the plate interface, and t 0 is the timing of a mainshock. Note that contributions from pre-seismic slip behaviors were ignored in this study ( a j (ξ, t) = 0 if t < t 0 ). The time derivative form of the (1) u i (x, t) = slip motion a j can be decomposed into slip motions due to the coseismic slip a c j and to the postseismic slip a p j as follows: where H represents a Heaviside step function. Note that the postseismic slip is defined as a combination of the afterslip (positive postseismic slip) and the fault locking (negative postseismic slip) in this study. According to this decomposition, we can rewrite Eq. (1) as follows: where the first and second terms on the right side of Eq. (3) represent a viscoelastic response to the coseismic slip and a viscoelastic response to the postseismic slip, respectively. We then decompose the viscoelastic response into an elastic (instantaneous) response and a viscous (time-delayed) response in terms of reformatting the Green's function as follows: where G e ij indicates the elastic Green's function and G v ij represents the viscous Green's function. We then can rewrite the first term on the right side of Eq. (3) as follows: where the first and second terms on the right side of Eq. (5) represent the surface displacement during the coseismic period due to an elastic response to the coseismic slip and the surface displacement during the postseismic period due to a viscous response to the coseismic slip, respectively. Note that the surface displacement during the postseismic period does not include the coseismic step and that the step function in Eq. (4) is considered to be 1 in Eq. (5), because we assume t ≥ t 0 here. Moreover, we can also rewrite the second term on the right-hand of Eq. (3) as follows: with: The first and second terms on the right side of Eq. (6) represent surface displacements during the postseismic period due to elastic and viscous responses to the cumulative postseismic slip caused from the time t 0 to the current time t , respectively. Combining Eqs. (5) and (6), Eq. (3) is rewritten as follows: in which the second term on the right side represents the viscous response to the postseismic slip and uniquely includes the time derivative form of postseismic slip and a temporal integral.
Here, an observation equation for observed coseismic displacement d c j is represented via u i as follows: where e c i indicates an error in the observed coseismic displacement. If we divide the postseismic period into K time windows, an observation equation for observed cumulative postseismic displacement d p k,i during the kth time window ( t k−1 ≤ t < t k ) can be written as follows: Tomita et al. Earth, Planets and Space (2020) 72:84 with: where e p k,i (x) indicates an observational error of the postseismic displacement. In this equation, the temporal integration is approximated by a simple summation of the time windows. Note that the third term in the equation is not required for the postseismic displacements during the first postseismic time window. Combining Eqs. (9) and (10), an observation equation that expresses both the co-and the post-seismic deformation is obtained in the vector form as follows:

Viscoelastic inversion method
To estimate the co-and the post-seismic slip distributions based on the observation equation, we provided two types of a priori information: smoothness of the coseismic slip distribution and smoothness of the postseismic slip distribution. Moreover, we investigated the relative weights between the co-and the post-seismic observational data as hyper-parameters (e.g., Funning et al. 2014), and we also investigated the appropriate viscosity in a lower viscoelastic layer and the thickness of an elastic layer (elastic thickness) as hyperparameters contained in the coefficient matrix of the viscoelastic Green's function when assuming a twolayered viscoelastic structure. Optimal values of these (10) hyper-parameters were determined by minimizing Akaike's Bayesian Information Criterion (ABIC, Akaike 1980;Yabuki and Matsu'ura 1992;Fukahata and Wright 2008). Although various previous studies included temporal smoothing of transient fault slips in the framework of ABIC (e.g., Fukahata et al. 2004;Yoshioka et al. 2015), we omit the temporal smoothing to deal with dramatic spatio-temporal evolution of the postseismic slip in this study. Properties of the co-and the post-seismic observational data were different in terms of magnitudes of displacements and of the preprocessing methods, although both data were obtained through geodetic observations. Thus, we determined their relative weight using ABIC, similar to determining the relative weights among different types of observational data (e.g., Funning et al. 2014). Assuming the error vectors, e c and e p k , follow Gaussian distributions with zero means, and covariances of σ 2 c E c and σ 2 p k E p k , respectively, an error matrix can be written as follows: with: where γ 2 k is a hyper-parameter representing the relative weight. A stochastic model for the data vector d can be then described as follows: where G(η, H) indicates the coefficient matrix of viscoelastic Green's function assuming a viscosity of η and the elastic thickness of H . N is the number of data (the dimensions of the vector d ). Although we calculated the viscoelastic Green's function assuming a two-layered viscoelastic half-space, as denoted later in this paper, the viscosity and elastic thickness are non-linear unknown parameters in the inversion problem. Thus, we optimized these parameters through the framework of Fukahata and Wright (2008), which enables us to determine the unknown parameter contained in the coefficient matrix using ABIC.
We impose spatial smoothing on the co-and post-seismic slip distributions using the Laplacian operator (e.g., Jónsson et al. 2002;Zhou et al. 2014) L as follows: where L c and L p k are the Laplacian operators for spatial smoothing on the co-and post-seismic slip distributions, respectively. Assuming that error vectors associated with constraints on the spatial smoothing for the co-and the post-seismic slip distributions follow Gaussian distributions with zero means and covariances of ρ 2 c I and ρ 2 p k I , the prior constraints are represented in the form of a probability density function (e.g., Yabuki and Matsu'ura 1992;Fukahata et al. 2004;Funning et al. 2014) as follows: where M is the number of parameters for slips (the dimensions of the vector a ). Based on Bayes' theorem (Bayes 1763), using Eqs. (16) and (19), the posterior probability density function is described as follows: with: where α 2 and β 2 k represent hyper-parameters adjusting the spatial smoothing for the co-and post-seismic slip distributions, which are defined as σ 2 c /ρ 2 c and σ 2 c /ρ 2 p k , respectively, and c is a normalizing factor independent of the model parameters and the hyper-parameters. The ABIC value is obtained by the following equation: where C ′ indicates a constant. Representing the optimum hyper-parameters that minimize ABIC as α 2 , β 2 k , γ 2 k , η and Ĥ , the optimal model parameters are obtained as: The covariance matrix, C, and the resolution matrix, R, are represented as:

General approach of synthetic tests
We performed two synthetic tests to clarify the utility of the viscoelastic inversion. In the first test, we (24) assumed a reverse faulting earthquake and its afterslip in a subduction zone, and in the second test, we assumed a shallow left-lateral strike-slip earthquake and its afterslip in an inland region. In both cases, we assigned a single time window in the postseismic period. Thus, the observation equation (Eq. 13) is simply rewritten as follows: For the synthetic tests, we assumed that â p expressed the first 1-year cumulative postseismic slip and that d p expressed the first 1-year cumulative postseismic displacements. We then generated co-and post-seismic displacements through Eq. (26) assuming synthetic slip distributions during the co-and the post-seismic periods. The viscoelastic Green functions were calculated assuming a two-layered viscoelastic half-space with an upper elastic layer and a lower viscoelastic layer Matsu'ura 2005, 2006). Although a rheology model considering a rapid temporal decay of the viscoelastic response (such as the Burgers rheology) has been often adopted (e.g., Pollitz 2003;Sun et al. 2014) for modeling the Earth rheology, this study assumed the linear Maxwell rheology to simply assess utility of the viscoelastic inversion method as a first step. Rigidity, density, and Poisson's ratio were set at 40 GPa, 2800 kg/m 3 , and 0.25 in the elastic layer, and at 67 GPa, 3300 kg/m 3 , and 0.28 in the viscoelastic layer, respectively. For the elastic thickness and viscosity in the viscoelastic layer, we assigned different values for the subduction zone and the inland cases. Finally, we estimated slip distributions based on Eq. (23). In the inversion, we searched for the optimal elastic thickness and viscosity based on ABIC (Eq. 22) using viscoelastic Green's function with various combinations of elastic thickness and viscosity values. In this study, to demonstrate the spatial resolutions of slip distributions, resolution values at sub-faults were calculated from the diagonal elements of the resolution matrix (Eq. 25). We estimated slip amounts for two orthogonal directions to express spatial variation of the slip direction. Thus, we defined a value of spatial resolution at each sub-fault as follows (e.g., Yoshioka and Matsuoka 2013): where R 1 and R 2 indicate the diagonal elements of the resolution matrix, which belong to the same sub-fault.
Moreover, we also estimated slip distributions by a conventional elastic inversion for comparison. Through the conventional elastic inversion, we first estimated a coseismic slip distribution from the coseismic displacements alone using the elastic Green's function. Then, we calculated viscoelastic responses on the surface due to the estimated coseismic slip and subtracted them from the postseismic displacements. Finally, we estimated a postseismic slip distribution from the subtracted postseismic displacements using the elastic Green's function. The details of the elastic inversion approach are explained in Additional file 1: Text S1. Orange triangles indicate synthetic GNSS sites for both co-and post-seismic period, while blue triangles indicate those only for post-seismic period. b A schematic image of cross section of rheological structure (the two-layered viscoelastic structure) Tomita et al. Earth, Planets and Space (2020) 72:84 Case study of strike-slip along an inland fault

Analysis conditions
We assumed a left-lateral strike-slip fault with a dip of 90° and with a strike of 0° from the north as shown by a green line in Fig. 1. The fault length is 60 km, and the fault domain is set from the surface to the depth of 20 km. We randomly distributed synthetic GNSS observation sites around the fault. For the coseismic period, 40 sites were widely distributed within a range of 100 km from the fault (orange triangles in Fig. 1a). In addition to these sites, 30 sites densely distributed within the range of 25 km from the fault were added for the postseismic period (blue triangles in Fig. 1a). Synthetic co-and post-seismic slip distributions are shown in Fig. 2a, b, respectively. We provided two separated peaks of coseismic rupture along the fault with maximum slips of 750 cm, while we provided one peak of afterslip along the fault with a maximum slip of 90 cm. The moment magnitudes (M w ) of the co-and post-seismic slips are 7.16 and 6.35 (seismic moments of 6.99 × 10 19 Nm and 4.21 × 10 18 Nm), respectively. Rakes of the coseismic slip and of the afterslip were fixed to 0° (pure left-lateral strike-slip). Co-and postseismic displacements were calculated by Eq. (26) using viscoelastic Green's function assuming the thickness of the elastic layer of 25 km and the viscosity of 3.0 × 10 18 Pa s. This viscosity is often found as optimal value in inland regions from geodetic modeling studies (e.g., Ohzono et al. 2012). We added Gaussian noises to the synthetic displacements: 2 cm, 4 cm, 0. In both the viscoelastic inversion and the elastic inversion, we restricted rakes of the estimated co-and post-seismic slip to fall within a range from − 45° to 45°; this was achieved by estimating slip in the direction of rake = − 45° and in the direction of rake = 45° with non-negative constraints. The non-negative constraint was applied using a calculation program for the bounded variables least squares problem (Lawson and Hanson 1995). Moreover, we provided zero slip constraint on the edges of the fault except the top edge for avoiding rank deficiency of the smoothing matrix (e.g., Iinuma 2009).

Figure 2c
, d shows the estimated co-and post-seismic slip distributions from the elastic inversion, respectively (elastic model). The estimation error and resolution for the elastic model are shown in Figs. 3a, b and 4a, b, respectively. The synthetic and calculated displacements are shown in Additional file 1: Figure S1, while those due to viscoelastic relaxation and postseismic slip in the postseismic period are shown in Additional file 1: Figure  S2. Although the magnitude of the estimated coseismic slip (M w 7.16, the seismic moment of 6.97 × 10 19 Nm) is comparable to that of the synthetic slip (M w 7.16), the estimated coseismic slip distribution is obviously oversmoothed. Therefore, the estimated slip distribution failed to reproduce the input slip distribution. Because the observation sites during the coseismic period are sparse, the spatial resolution is generally low (Fig. 4a). Consequently, a strong smoothing constraint was required to obtain stable solutions. Besides, the estimated postseismic slip distribution is also over-smoothed obviously although the observational sites are added in the postseismic period. The magnitude of the estimated postseismic slip (M w 6.31, the seismic moment of 3.68 × 10 18 Nm) is underestimated (synthetically, M w 6.35). Since the estimated coseismic slip distribution provided slightly incorrect contributions of the postseismic viscoelastic responses (Additional file 1: Figure S2a), the postseismic displacements subtracting the viscoelastic responses did not properly express the contributions of the postseismic slip. Consequently, the estimated postseismic slip distribution is distorted, although the spatial resolution during the postseismic period is improved compared to the distribution during the coseismic period (Fig. 4b).
Figure  Figure S4. The viscoelastic model demonstrated two distinguishable peaks of coseismic slip and deep coseismic slips, which cannot be reproduced by the elastic model, and the magnitude of the coseismic slip was estimated to be M w 7.21 (the seismic moment of 8.27 × 10 19 Nm). This improvement in the viscoelastic model is because the viscoelastic inversion improved the spatial resolution of the coseismic slip distribution by introducing postseismic observational data via the viscoelastic Green's function, especially deeper part of the fault zone as pointed out later. Since the coseismic slip distribution was improved in the viscoelastic model, its viscoelastic responses were accurately calculated (Additional file 1: Figure S4a). Thus, the postseismic slip distribution was also accurately estimated by utilizing the viscoelastic inversion (Fig. 2f ), although its maximum slip amount was slightly underestimated. Unlike the elastic model, the magnitude of the estimated postseismic slip (M w 6.35, the seismic moment of 4.16 × 10 18 Nm) is comparable to that of the synthetic slip (M w 6.35).
Through the viscoelastic inversion, we determined the optimal elastic thickness and the optimal viscosity in the viscoelastic layer. Additional file 1: Figure S5 shows the difference of ABIC values depending on these two parameters. The parameters tend to show a trade-off relationship, but the minimum ABIC value successfully shows the assumed elastic thickness and the assumed viscosity. Additional file 1: Figure S6 shows the viscoelastic model without rake constraints. In this inversion condition, the estimated slip distributions for both the co-and post-seismic periods are over-smoothed to some extent. Thus, the rake constraints are useful for obtaining plausible slip distributions.
The estimation error of the coseismic slip in the viscoelastic model (Fig. 3c) is generally smaller than that in the elastic model (Fig. 3a), while the estimation error of the postseismic slip in the viscoelastic model (Fig. 3d) is larger than that in the elastic model (Fig. 3b). Conversely, resolution of the coseismic slip in the viscoelastic model is lower than that in the elastic model except the deep portion of the fault zone (Fig. 4e), while resolution of postseismic slip in the viscoelastic model is higher than that in the elastic model (Fig. 4f ). This relationship is basically controlled by the weight of the spatial smoothing constraints; stronger and weaker spatial smoothing constraints were applied to the co-and the post-seismic slip distributions in the viscoelastic model, respectively. Although the strong smoothing constraint was applied to the coseismic slip distribution in the viscoelastic model (note that strong smoothing constraint degraded spatial resolution), slightly high spatial resolution was obtained in the deep portion of the fault zone (Fig. 4e). This indicated that the viscoelastic responses observed by the postseismic geodetic observational sites improved the spatial resolution of the coseismic slip in the viscoelastic model.
In order to further demonstrate the utility of the viscoelastic inversion, we investigated improvement of the spatial resolution of the coseismic slip assuming various observational site distributions. Figure 5 shows the difference of spatial resolution between the elastic model and the viscoelastic model assuming the same smoothing and weighting parameters for both the elastic and the viscoelastic models. All site distributions improve spatial resolution in the deep fault zone by considering the viscoelastic responses, even if no observational site is added in the postseismic period (Fig. 5a). Especially, the postseismic observational sites near the fault trace on the ground surface are useful for improving the deep spatial resolution (Fig. 5b, d). However, it is difficult to improve spatial resolution of the shallow fault zone regardless of site distributions because viscoelastic responses to the shallow coseismic rupture are generally smaller than those of the deep coseismic rupture. Such investigation on difference of the spatial resolution is useful for determining adequate extension of the geodetic observational network.

Case study of dip-slip along a plate interface in a subduction zone Analysis conditions
We assumed a dip-slip fault with a dip of 15° imitating a plate interface in a subduction zone as shown in Fig. 6. The fault domain was 500 km (along strike) × 310 km (along dip), and its edge reached to the surface (i.e., trench). For the coseismic period, we randomly distributed 150 synthetic sites within a range of 200-400 km from the trench, which imitated onshore GNSS sites. In addition to these sites, 30 sites densely distributed within a range of 0-150 km from the trench were added for the postseismic period, which imitated offshore GNSS-Acoustic (GNSS-A) sites. This geodetic site distribution is similar to the observational sites in the subduction zones around Japan (Nankai and Tohoku regions). Especially, the additional GNSS-A sites in the postseismic period imitated those of Tohoku region after the 2011 Tohoku-oki earthquake (Kido et al. 2015).
Synthetic co-and post-seismic slip distributions are shown in Fig. 7a, c, respectively. We provided three separated peaks of coseismic rupture along the plate interface. Their spatial extents, depths, and maximum slip amounts are different, and the peak at an intermediate depth represents the main rupture, with a maximum slip amount of 42 m denoting a large subduction earthquake. The moment's magnitude of the synthetic coseismic slip is M w 8.85 (the seismic moment of 2.34 × 10 22 Nm). We provided two separated peaks of afterslip along the plate interface with a maximum slip of 1.6 m. The moment magnitude of the synthetic postseismic slip is M w 7.87 (the seismic moment of 7.91 × 10 20 Nm). The rakes of the coseismic slip and of the afterslip were fixed to 90° (pure dip-slip). Co-and post-seismic displacements were calculated by using viscoelastic Green's function assuming a viscosity of 1.5 × 10 19 Pa s and a thickness of the elastic layer of 50 km following Eq. (26). Note that we added Gaussian noises to the synthetic displacements. For GNSS sites (over 200 km far from the trench), these observational noises are 2.0 cm, 4.0 cm, 0.5 cm/year, and 1.0 cm/year in 1σ for the coseismic horizontal, the coseismic vertical, the postseismic horizontal, the postseismic vertical components, respectively. For the postseismic GNSS-A sites, they are 2.0 cm/year and 4.0 cm/year in 1σ for the horizontal and the vertical components, respectively.
In the same manner as the strike-slip case, we restricted rakes of the estimated co-and post-seismic slip. For the coseismic slip, the rake was restricted within the range from 45° to 135° by estimating slip in the direction of rake = 45° and in the direction of rake = 135° with non-negative constraints. For the postseismic slip, we only estimated slips in the direction of rake = 90°, and  we restricted slips over − 8.3 cm/year. The negative postseismic slip represents the contributions of the interplate locking, so we restricted the negative postseismic slip within a plate convergence rate; in this study, we assumed the plate convergence rate of the subduction zone in the off-Tohoku region, Japan. The plate convergence rate is calculated as the plate motion rate of the Pacific plate relative to the North American plate by the MORVEL model (DeMets et al. 2010). The non-negative constraint and the restriction of the negative postseismic slip were performed by the program of Lawson and Hanson (1995) as with the strike-slip case. Additionally, we also provided zero slip constraint on the edges of the fault except the trench edge. Results Figure 7b, f shows the estimated co-and post-seismic slip distributions from the elastic inversion, respectively (elastic model). The estimation error and resolution for the elastic model are shown in Figs. 8a, c, and 9a, d, respectively. The synthetic and calculated displacements are shown in Additional file 1: Figure S7, while those due to viscoelastic relaxation and postseismic slip are shown in Additional file 1: Figure S8. Since the spatial resolution near the trench is obviously low (Fig. 9a), the estimated coseismic slip distribution cannot reproduce the shallowest peak. Furthermore, the spatial resolution could not adequately decompose the main peak at the intermediate depth and the deepest peak; therefore, the estimated coseismic slip distribution was slightly over-smoothed. Note that the magnitude of the estimated coseismic slip is M w 8.80 (the seismic moment of 1.99 × 10 22 Nm) is slightly underestimated compared with that of the input slip (M w 8.85). In contrast, the estimated postseismic slip distribution successfully reproduced the assumed slip distribution although a negative postseismic slip partially occurred. The magnitude of the estimated positive postseismic slip is M w 7.90 (the seismic moment of 9.13 × 10 20 Nm). This is because the postseismic GNSS-A sites sufficiently resolve the peaks of the postseismic slip, even those located near the trench (Fig. 9d). Moreover, since the viscoelastic responses near the peaks of the postseismic slip are small (Additional file 1: Figure S8a, b), the estimated postseismic slip distribution was not distorted, unlike the strike-slip case. respectively (viscoelastic model). The synthetic and calculated displacements are shown in Additional file 1: Figure S9, while those due to viscoelastic relaxation and postseismic slip are shown in Additional file 1: Figure  S10. The spatial resolution of the coseismic slip distribution in the viscoelastic model improved compared with that of the elastic model (Fig. 9a-c); especially regions where the postseismic observation sites are located demonstrated improvement of the spatial resolution on the coseismic slip. The viscoelastic model provides higher maximal coseismic slip than that estimated by the elastic model, which is suitable to reproduce the input coseismic slip distribution. However, it is difficult to reproduce the shallowest coseismic slip patch even by the viscoelastic model, and then the magnitude of the estimated coseismic slip was M w 8.81 (the seismic moment of 2.05 × 10 22 Nm), which is comparable to that of the input slip (M w = 8.85). This is because the spatial resolution near the trench is still low (Fig. 9b, c) and because a trade-off relationship between viscoelastic responses and the negative postseismic slip (interplate coupling effect). The large coseismic rupture in the shallow portion of the plate interface provided substantial motion to the down-dip direction on the surface in the postseismic period because of viscoelastic relaxation. However, since the negative postseismic slip can also produce motion to the down-dip direction, the contributions of the negative postseismic slip and the viscoelastic responses potentially show a trade-off relationship. Thus, if we somehow remove interplate coupling effects from the postseismic observational data, the viscoelastic model constraining negative postseismic slip (constrained viscoelastic model) can well reproduce the input coseismic slip distribution even including the shallowest patch (Fig. 7d). For the postseismic slip, the viscoelastic model successfully resolved the assumed slip distribution as well as the elastic model ( Fig. 7e-g). The magnitude of the estimated positive postseismic slip was M w 7.91 (the seismic moment of 9.23 × 10 20 Nm), which is almost same as that obtained from the elastic model. Moreover, both of the elastic and the viscoelastic models show similar the estimation error and the spatial resolution on the postseismic slip.
The optimal elastic thickness and the optimal viscosity in the viscoelastic layer are determined by the minimum ABIC framework as shown in Additional file 1: Figure  S11. We successfully identified the assumed values of these parameters through the framework.
We investigated the slip distributions by the viscoelastic inversion without the rake constraint on the coseismic slip and the constraint on the negative postseismic slip (Additional file 1: Figure S12). As a result, we obtained an over-smoothed distribution for the coseismic slip. Furthermore, we obtained an unnatural negative postseismic slip; the postseismic slip of approximately − 30 cm/year was estimated where the maximum peak of the coseismic rupture was assigned. This unnatural negative postseismic slip was caused by the trade-off relationship between the viscoelastic responses and the negative postseismic slip as indicated above. Therefore, the constraint on the negative postseismic slip is an essential analysis condition for a viscoelastic inversion when estimating co-and postseismic slip distributions in a subduction zone.
In the same manner as the strike-slip case, we investigated improvement of the spatial resolution of the coseismic slip assuming various observational site distributions. Figure 10 shows the difference of spatial resolution between the elastic model and the viscoelastic model assuming the same smoothing and weighting parameters for both the elastic and the viscoelastic models. Additional observational sites in the postseismic period basically increase the spatial resolution of the coseismic slip for sub-faults just below the sites (Fig. 10b-e). However, no matter how close to the trench the additional observational sites on the upper plate are, spatial resolution near the trench in the coseismic plate is hardly improved (Fig. 10d). To improve the spatial resolution near the trench, the additional observational sites off the trench are effective (Fig. 10f ). Moreover, it is effective in the improvement of the spatial resolution to deploy multiple observational sites in the along-dip direction (Fig. 10b, c).

Data and analysis conditions
In this study, we applied a viscoelastic inversion to the co-and post-seismic deformations associated with the 2011 Tohoku-oki earthquake. To investigate the utility of the viscoelastic inversion, the method can be applied to an event in which a geodetic observation network is greatly extended during the postseismic period. The 2011 Tohoku-oki earthquake is among the largest recorded megathrust earthquakes with a magnitude of 9.0, and its coseismic deformation was recorded not only by onshore geodetic observations (e.g., onshore GNSS observational sites shown as circles in Additional file 1: Figure S13), but also by seafloor geodetic observations such as ocean bottom pressure (OBP) gauges (shown as inversed triangles in Additional file 1: Figure S13) and seafloor GNSS-Acoustic (GNSS-A) observations (shown as triangles and diamonds in Additional file 1: Figure S13). However, these seafloor geodetic observations were too sparse to monitor deformation in the entire off-Tohoku region. After the occurrence of the 2011 Tohoku-oki earthquake, 20 GNSS-A sites were newly installed in September 2012 (Kido et al. 2015) to monitor postseismic deformation (shown as squares in Additional file 1: Figure S13), and the overall features of the postseismic deformation along the Japan trench were successfully detected (Tomita et al. 2017). Thus, using the extended GNSS-A observation data during the postseismic period, we investigated the utility of the viscoelastic inversion method to the actual observational data. Although the viscoelastic inversion method has been performed for modeling the co-and post-seismic slip distributions (Yamagiwa et al. 2015;Ito et al. 2018), these studies did not include the extensive GNSS-A observation data during the postseismic period.
For the coseismic observational data, we employed three types of geodetic observational results: onshore GNSS observations, OBPs (Ito et al. 2011b;Maeda et al. 2011), and GNSS-A observations (Sato et al. 2011;Kido et al. 2011). The details of these observational data and their pre-processing are shown in Additional file 1: Text S2.
For the postseismic observational data, we employed two types of geodetic observational results: onshore GNSS observations and GNSS-A observations (Tomita et al. 2017;Yokota et al. 2018). Since the GNSS-A observations were conducted by campaign-style surveys, the temporal resolution of the GNSS-A observational data during the postseismic period was low, unlike the onshore GNSS data. Thus, we assigned a one time-window during the postseismic period in the inversion analysis. Accordingly, the observation equation of this analysis b e c a f d Fig. 10 Improvement of the spatial resolution assuming various postseismic site distributions for the dip-slip case. Each panel shows difference of the spatial resolution on the coseismic slip same as with Fig. 9c assuming different site distributions. Gray triangles indicate the observational sites during both the co-and the post-seismic periods, while inverse gray triangles indicate the observational sites only during the postseismic period is expressed by Eq. (26) in the same manner as the synthetic tests. The time window was set to be compatible with the GNSS-A observational data obtained by the extended GNSS-A observational network; therefore, its period was from September 2012 to September 2016. We calculated cumulative postseismic displacements of the onshore GNSS and offshore GNSS-A observations during the time window. Note that the postseismic displacements were obtained relative to the North American plate (DeMets et al. 2010). We then estimated the cumulative postseismic slip distribution through the inversion analysis. The details of these observational data and their pre-processing are shown in Additional file 1: Text S2.
To perform a viscoelastic inversion, the viscoelastic Green's function was calculated in a two-layered viscoelastic structure Matsu'ura 2005, 2006) assuming the same rigidity, density, and Poisson's ratio with the synthetic tests. Since we did not employ early postseismic observational data, we utilized the linear Maxwell rheology without considering rapid temporal decay of the viscoelastic responses and contribution of the poroelastic rebound. We employed a compiled model from various studies (Nakajima and Hasegawa 2006;Nakajima et al. 2009;Kita et al. 2010;Hirose et al. 2008) for the geometry of the plate interface. We calculated the viscoelastic response at each point of 5-km grids on the plate interface. Then, we calculated a viscoelastic response on each sub-fault at spacing of 20 km as a viscoelastic Green's function by summing all the viscoelastic responses of 5-km grids within the sub-fault (there were 482 sub-faults in total). This methodology was used to reflect effects of the detailed (5 km grids) plate interface geometry on the viscoelastic Green's function. In calculating the viscoelastic responses for the seafloor geodetic sites, we considered seafloor depths by biasing the depths of the sub-faults along the plate interface in the same manner as Iinuma et al. (2012). Note that we fixed an elastic thickness of 50 km in the off-Tohoku region in this same manner as previous studies in this region (Diao et al. 2013;Yamagiwa et al. 2015) to reduce calculation time. Note that we determined the optimal viscosity by the ABIC framework in the same manner as the synthetic tests. Like the dip-slip case of the synthetic tests, the rake of the coseismic slip was constrained between 45° and 135°, and the negative postseismic slip was restricted to over − 8.3 cm/year. Zero slip constraints were assigned to the edge of the fault zone, except the trench.

Coseismic slip distribution estimated by elastic inversion
We first show the coseismic slip distribution, which was estimated by applying the elastic inversion (the elastic model) in Fig. 11a. Note that the masked patches in the slip distribution show a low-resolution area [resolution < 0.05, calculated by Eq. (27) as the same criteria with Yokota et al. (2016)]. Detailed estimation results for slip directions, estimation errors, resolution, and misfits between the observed surface displacements and the calculated surface displacements, are shown in Fig. 12a, b, e, and f, respectively. The estimated seismic moment in the elastic model is 4.18 × 10 22 Nm, which is equivalent to M w 9.01. In the elastic model, the primary rupture area (PRA), which is defined as a large coseismic slip area with a slip greater than 20 m, is concentrated in the off-Miyagi region (Central region of Tohoku-oki, Additional file 1: Figure S13) where the spatial resolution is high (Figs. 11a  and 12e). The estimated coseismic slip distribution is similar to those in previous models that have used for both on-shore and off-shore geodetic data (e.g., Ito et al. 2011a;Ozawa et al. 2012;Sato et al. 2013;Silverii et al. 2014). However, some previous studies have indicated that the large coseismic slip was concentrated further near the trench; such results were dependent on boundary conditions and/or on the weights of data (e.g., Iinuma et al. 2012;Zhou et al. 2014;Hashima et al. 2016). to be 1.5 × 10 19 Pa s by minimizing ABIC (Additional file 1: Figure S14). Detailed estimation results for the coseismic slip are shown in Fig. 12c, d, g, h, and those for the postseismic slip are shown in Fig. 13. The seismic moment of the coseismic slip is 4.46 × 10 22 Nm, which is equivalent to M w 9.03, while the cumulative seismic moment of postseismic slips during the data period is 3.29 × 10 21 Nm, which is equivalent to M w 8.28.

Co-and post-seismic slip distributions estimated by viscoelastic inversion
PRA in the viscoelastic model is wider than that in the elastic model, and it extends to the southern region (Fig. 11a, b). The extended PRA spatially corresponds to areas where the resolution of the coseismic slip in the viscoelastic model is improved from that in the elastic model, which is due to the introduction of the postseismic observational data (Figs. 11a, b and 12h). As large coseismic ruptures generally generate postseismic landward motions above the rupture area via viscoelastic relaxation (e.g., Sun et al. 2014), the extension of PRA in the viscoelastic model explains the postseismic landward displacements in the offshore region, which can be prominently observed within the new GNSS-A network (Fig. 13c). The coseismic RMS misfits in the elastic model for offshore and onshore data are 45.96 cm and 2.99 cm, respectively, while those in the viscoelastic model are 29.44 cm and 2.90 cm, respectively. The misfits of the coseismic displacements in the viscoelastic model are almost the same as those in the elastic model (pattern of misfits: Fig. 12a  . Black vectors shown in the columns of "All", "Postseismic Slip", and "Viscoelastic Relaxation" represent the observed postseismic displacements. Black vectors shown in the columns of "Residuals" represent the misfits between the observation and the calculation. Blue vectors shown in the columns of "All" represent the calculated postseismic displacements. Blue vectors shown in the columns of "Postseismic Slip" and "Viscoelastic Relaxation" represent contributions of postseismic slip and viscoelastic relaxation of the calculated postseismic displacements, respectively. Red and blue bars shown in the columns of "All", "Postseismic Slip", and "Viscoelastic Relaxation" represent the observed postseismic displacements in uplift and subsidence directions, respectively. Red and blue bars shown in the column of "Residuals" represent the misfits. Orange and cyan bars represent the calculated postseismic displacements in uplift and subsidence directions, respectively, and they show contributions of the labeled postseismic deformation process. Colors of sub-faults in a, d, and e show postseismic slip, 1σ file 1: Figure S14). The postseismic RMS misfits in the viscoelastic model are 2.25 cm/year and 1.10 cm/year for the offshore and onshore data, respectively. The viscoelastic model indicates positive postseismic slip (afterslip) areas in the downdip region of the PRA and in the off-Ibaraki and off-Fukushima regions (Southern regions of Tohoku-oki, Additional file 1: Figure S13) including near the trench (Fig. 11c). In addition, a negative postseismic slip (fault locking) area is estimated in the off-Miyagi region overlapping with the PRA (Fig. 11b, c).
As shown above, we successfully improved the spatial resolution of the coseismic slip distribution using the postseismic GNSS-A observational data through the viscoelastic inversion. However, the estimated slip distributions likely have systematic errors, especially in the coastal area of the postseismic slip distribution, due to the assumed viscoelastic rheology. The influence of the viscoelastic rheology is discussed in the next section.

Influence of viscoelastic structure
The optimal viscosity was determined to be 1.5 × 10 19 Pa s through the viscoelastic inversion. This viscosity is similar to that estimated by Diao et al. (2013) (2.0 × 10 19 Pa s) assuming a layered viscoelastic structure, but it is generally larger than asthenosphere viscosities employed in other studies assuming a realistic 3-dimensional viscoelastic structure (e.g., Sun et al. 2014;Hu et al. 2016;Suito 2017;Freed et al. 2017). This difference in the viscosity can be explained as follows: (1) No data from the early postseismic period was used in this study. Because the viscoelastic responses during the early postseismic period are generally large due to transient creep (e.g., Wang et al. 2012), the apparent viscosity in the late period appears to be relatively large.
(2) A two-layered viscoelastic half-space was employed to calculate the viscoelastic responses in this study; as the low viscosity media extended deeper in the twolayered viscoelastic half-space, this simple rheological structure provided much more viscoelastic deformation than the 3-dimensional viscoelastic models considering an oceanic slab and different viscosities in the oceanic and continental asthenospheric mantles. Therefore, the apparent viscosity in the two-layered viscoelastic structure also appears to be high. If we assumed much lower viscosities, a physically unnatural situation would take place with a substantial afterslip having occurred in the PRA; this also supports the suitability of our optimal viscosity method.
The viscoelastic responses obtained in this study accurately reproduce the postseismic landward motions in the offshore area, especially prominent near the trench (Fig. 13c), which should be mainly explained by viscoelastic relaxation (e.g., Watanabe et al. 2014;Suito 2017). Thus, the obtained viscosity is plausible at least for the offshore area. However, the layered a b c Fig. 14 Co-and post-seismic slip models without using coseismic seafloor observations. a-c The co-and post-seismic slip models estimated by the elastic and viscoelastic inversions without using seafloor geodetic observational data during the coseismic period. These panels are drawn in the same manner as with Fig 11, except the gray patches expressing low spatial resolution in a and b. Gray-scaled patches in panels a and b indicate degree of low spatial resolution (resolution < 0.05) viscoelastic structure provides the high viscosity, even under the onshore area where viscosity is generally low relative to the offshore area (e.g., Muto et al. 2013;Suito 2017). Thus, the viscoelastic responses in the onshore region might be underestimated in this study. This effect, which is due to the layered viscoelastic structure, might produce overestimations of the afterslip near the coastline as discussed by Wang et al. (2018). Therefore, to obtain more reliable slip distributions, viscoelastic Green's function should be calculated assuming a realistic 3-dimensional viscoelastic structure as similar to Ito et al. (2018). Because this study aims to investigate improvements on the spatial resolution of a coseismic slip distribution using postseismic geodetic observations from the viscoelastic inversion method, we will introduce a 3-dimensional viscoelastic structure in a future study.

Importance of extension of postseismic geodetic observational network
We finally emphasize the effectiveness of using postseismic geodetic data through the viscoelastic inversion method to provide constraints on the coseismic slip distribution. Reinforcement using geodetic observations near a mainshock fault, even after the occurrence of an earthquake, is valuable if the earthquake is large enough to induce considerable viscoelastic deformation. In our demonstration of the usefulness of postseismic geodetic data, we estimated the co-and post-seismic slip distributions of the 2011 Tohoku-oki earthquake without using coseismic offshore data (Fig. 14). This experiment assumes a situation whereby a large earthquake occurred in a region where spatially poor geodetic observations existed during the coseismic period. Using the elastic inversion, the estimated coseismic slip was broadly distributed along the trench, but it did not reach into the trench (Fig. 14a); this result is in accordance with previous source models derived by onshore geodetic data alone (e.g., Miyazaki et al. 2011;Iinuma et al. 2011). In contrast, the viscoelastic inversion method provided superior performance in revealing the important characteristics of the Tohoku earthquake, such as the concentration of the coseismic rupture in the off-Miyagi region and the coseismic rupture near the trench (Fig. 14b).
Recoveries of coseismic slip distributions through the viscoelastic inversion method will likely also be applicable to other large earthquakes where an extensive viscoelastic relaxation occurred, such as in the Sumatra or Chile subduction zones. Furthermore, extensive postseismic seafloor geodetic observations provide valuable data for determining whether large subduction earthquakes are accompanied by a shallow coseismic slip; this information cannot be easily constrained by coseismic observational data obtained in a far field setting.

Conclusion
Focusing on the link between a coseismic slip distribution and postseismic deformation via viscoelastic relaxation, we developed an inversion approach that simultaneously estimates co-and post-seismic slip distributions from co-and post-seismic geodetic observational data using the viscoelastic Green's function. We applied this viscoelastic inversion to two cases of synthetic observational data: a strike-slip event in an inland region and a dip-slip event in a subduction zone. The synthetic tests demonstrated the viscoelastic inversion method's considerable improvement of the spatial resolution of coseismic slip distributions when the geodetic observational network was extended during the postseismic period. Moreover, we also applied the viscoelastic inversion method to the actual observational data associated with the 2011 Tohoku-oki earthquake. The seafloor GNSS-A observational network, which was extended after the occurrence of the mainshock, enabled us to improve the spatial resolution of the coseismic slip distribution through the viscoelastic inversion method. The estimation results showed characteristic slip distributions such as a large coseismic rupture near the trench in the central region of Tohoku-oki and a large afterslip near the trench in the southern region of Tohoku-oki. However, since we employed a simple two-layered viscoelastic structure to calculate viscoelastic Green's function, the estimated slip distributions might have large systematic errors, especially for the postseismic slip distribution near the coastline. Introducing a 3-dimensional viscoelastic structure is therefore important for further precise modeling, and it will be our next topic of research. Finally, a geodetic observational network (especially a seafloor geodetic network for an event in a subduction zone) should be installed or extended, even after a massive seismic event, to improve the spatial resolution of a coseismic slip distribution and to assess coseismic rupture behaviors through the viscoelastic inversion approach.
Additional file 1: Text S1. Elastic inversion method. Text S2. Pre-processing for observational data. Figure S1. Co-and post-seismic displacements for the strike-slip case in the elastic model. Figure S2. Contributions of viscoelastic relaxation and postseismic slip for the strike-slip case in the elastic model. Figure S3. Co-and post-seismic displacements for the strike-slip case in the viscoelastic model. Figure S4. Contributions of viscoelastic relaxation and postseismic slip for the strike-slip case in the elastic model. Figure S5. ABIC values assuming various combination of elastic thickness and viscosity for the strike-slip case. Figure S6. Co-and post-seismic slip distributions in the viscoelastic inversion without rake constraints for the strike-slip case. Figure S7. Co-and post-seismic displacements for the dip-slip case in the elastic model. Figure S8.
Contributions of viscoelastic relaxation and postseismic slip for the dip-slip case in the elastic model. Figure S9. Co-and post-seismic displacements for the dip-slip case in the viscoelastic model. Figure S10. Contributions of viscoelastic relaxation and postseismic slip for the dip-slip case in the viscoelastic model. Figure S11. ABIC values assuming various combination of elastic thickness and viscosity for the dip-slip case. Figure S12.
Co-and post-seismic slip distributions in the viscoelastic inversion without rake and negative postseismic slip constraints for the dip-slip case. Figure  S13. Geodetic observation site distribution off Tohoku region. Figure S14. RMS (root mean square) misfits and ABIC values assuming various viscosities. Figure S15. Examples of timeseries of postseismic displacements at onshore GNSS sites. Table S1. Postseismic displacement rate at GNSS-A observation sites.