Abstract

Large tectonic earthquakes lead to significant deformations in the months and years thereafter. These so-called post-seismic deformations include contributions mainly from afterslip and viscoelastic relaxation, quantification of their relative influence is of importance for understanding the evolution of post-seismic crustal stress, strain and aftershocks. Here, we investigate the post-seismic deformation processes following the 2011 Mw 9.0 Tohoku earthquake using surface displacement data as observed by the onshore global positioning system network in the first ∼1.5 yr following the main shock. We explore two different inversion modelling strategies: (i) we simulate pure afterslip and (ii) we simulate the combined effect of afterslip and viscoelastic relaxation. By assuming that the afterslip is solely responsible for the observed post-seismic deformation, we find most afterslip activities to be located close to the downdip area of the coseismic rupture at 20–80 km depth with a maximum cumulative slip of ∼3.8 m and a seismic moment of 2.3 × 1022 Nm, equivalent in moment to an Mw 8.84 earthquake. By assuming a combination of afterslip and viscoelastic components, the best data fit is found for an afterslip portion that is spatially consistent with the pure afterslip model, but reveals a decreased seismic moment of 2.1 × 1022 Nm, or Mw 8.82. In addition, the combined model suggests an effective thickness of the elastic crust of ∼50 km overlying an asthenosphere with a Maxwell viscosity of 2 × 1019 Pa s. Temporal analysis of our model inversions suggests that the rate of afterslip rapidly decreases with time, consistent with the state- and rate-strengthening frictional law. The spatial pattern of afterslip coincides with the locations of aftershocks, and also with the area of coseismically increased Coulomb failure stress (CFS). Only a small part of the coseismically increased CFS was released by the afterslip in 564 d after the event. The effect of the viscoelastic relaxation within this initial stage only plays a secondary role, but it shows an increasing tendency, that is, the contribution of viscoelastic relaxation increases with time. Further geodetic observations are needed for a robust quantification of the role of the viscoelastic relaxation in the post-seismic deformation.

1 INTRODUCTION

At subduction zones, the seismic cycle is depicted by time-dependent surface deformation due to thrust faulting, involving an elastic and a viscoelastic compartment (Thatcher & Rundle 1979). Associated with an earthquake, in specific, the Earth's upper layer behaves foremost elastically (instantaneously), whereas the deeper regions of the lower crust and upper mantle behave viscoelastic (time-dependent). The rapid development of spaceborne geodetic techniques, such as global positioning system (GPS) and interferometric synthetic aperture radar (InSAR), has brought valuable data sets for studying such seismic cycles at different scales and stages. This has led to construct coseismic slip models with high spatial detail (e.g. Shen et al. 2009; Ozawa et al. 2011; Simons et al. 2011; Diao et al. 2012), and post-seismic models that allow analysing transient stress and strain changes occurring over the months to decades following an earthquake, relevant for identifying regions of fault slip deficits and assessment of regional seismic hazards (e.g. Diao et al. 2011a; Ozawa et al. 2012).

Post-seismic deformation processes following large earthquakes include (i) aseismic afterslip surrounding the coseismic slip area, (ii) viscoelastic relaxation of the lower crust or upper mantle, (iii) poroelastic rebound of the medium near the fault plane and (iv) combined effects that include two or more of those processes (Segall & Davis 1997; Jónsson et al. 2003; Hsu et al. 2006; Pollitz et al. 2006; Freed 2007; Bürgmann & Dresen 2008; Bruhat et al. 2011; Diao et al. 2011a; Hoechner et al. 2012; Wang et al. 2012). The location and amount of afterslip at a given fault and/or the rheological structure under the seismogenic zone, which are indispensable for fault evolution research and regional geodynamic modelling, can be inferred from model simulations of post-seismic deformation data (Bürgmann & Dresen 2008).

It is usually assumed that poroelastic rebound can only affect post-seismic deformation within tens of kilometres (km) near fault ruptures, whereas afterslip and viscoelastic relaxation are considered the two main mechanisms for observed long-term and wide-range post-seismic deformations. However, distinguishing the difference between deep afterslip and viscoelastic relaxation by using limited data is difficult, because of the similarities in their induced deformations (Bürgmann & Dresen 2008). For this purpose, further observations to increase the spatial and temporal data coverage are required. In this work, we make the attempt to model post-seismic deformation occurrence for pure afterslip and the combination of afterslip and viscoelastic relaxation.

The Tohoku earthquake occurred off the coast of northeastern Japan, at a site where the Pacific Plate subducts beneath the North American Plate and the Philippine Plate at a rate of ∼9 cm yr−1 (DeMets 1992). Several large earthquakes of magnitude 7.0–8.0 have occurred in the 20th century along this plate boundary (Fig. 1). Onshore GPS observations collected in the years before the Tohoku earthquake already allowed to investigate the slip deficit along megathrust fault, allowing to identify that the plate interface near Sendai is almost completely ‘locked’ (Hashimoto et al. 2009; Loveless & Meade 2010). On 2011 March 11, the Mw 9.0 Tohoku earthquake occurred at the segment of the trench 230 km east of Sendai. This earthquake and associated secondary disasters, including tsunami and a nuclear leakage, caused massive damage to human life and infrastructure. Also aftershocks may be damaging for the months to years following the main shock, locations and parameters of which can be better assessed by analysing the post-seismic stress and strain.

Tectonic setting of the Mw 9.0 Tohoku earthquake. The red beachballs show the focal mechanisms of the main shock and the Mw 7.9 aftershock. The purple dashed line is the surface projection of the fault plane used for afterslip inversion. The green stars and grey ellipses indicate epicentres and source regions, for large interplate earthquakes (≥Mw 7.5) that occurred in the past century, respectively. The white dots are GPS stations. Black dot in the top left panel is the reference site of the post-seismic displacements. AM, the Amurian Plate; PH, the Philippine Plate; PA, the Pacific Plate; NA, the North American Plate.
Figure 1.

Tectonic setting of the Mw 9.0 Tohoku earthquake. The red beachballs show the focal mechanisms of the main shock and the Mw 7.9 aftershock. The purple dashed line is the surface projection of the fault plane used for afterslip inversion. The green stars and grey ellipses indicate epicentres and source regions, for large interplate earthquakes (≥Mw 7.5) that occurred in the past century, respectively. The white dots are GPS stations. Black dot in the top left panel is the reference site of the post-seismic displacements. AM, the Amurian Plate; PH, the Philippine Plate; PA, the Pacific Plate; NA, the North American Plate.

The Tohoku earthquake caused a coseismic slip at the megathrust fault that has been constructed based on onshore and offshore GPS and acoustic data. The inverted coseismic slip, which has a slip maximum of ∼50 m, was found to be very shallow, concentrated at the updip area of the hypocentre and close to the trench (Simons et al. 2011; Diao et al. 2012; Ozawa et al. 2012). This earthquake caused a crustal displacement that severely changed the stress state not only on the rupture fault, but also on the surrounding media leading to post-seismic deformation (Ozawa et al. 2012). Recent studies attribute the initial post-seismic deformation following the Tohoku earthquake to the aseismic fault slip (afterslip), whereas other mechanisms such as viscoelastic relaxation are neglected (Ozawa et al. 2012). As viscoelastic relaxation, however, mainly affects the long-term and far-field post-seismic deformation (Lorenzo-Martín et al. 2006; Wang et al. 2009; Bruhat et al. 2011; Diao et al. 2011a), its discounting may change the afterslip pattern, lead to an overestimation of the afterslip and consequently to an incorrect source process interpretation.

In this study, we investigate the time-dependent post-seismic deformation in the ∼1.5 yr (564 d) following the Tohoku earthquake by using the afterslip model and the combined model of afterslip and viscoelastic relaxation.

2 DATA AND OBSERVATIONS

The GPS network in Japan is of high quality and very dense, now providing valuable data sets for studying post-seismic processes associated with the 2011 Tohoku earthquake. We make use of the GPS displacement data recorded in the first 564 d after the 2011 Tohoku earthquake, which are divided into eight periods (0–7, 8–20, 21–48, 49–107, 108–200, 201–289, 290–381 and 382–564 d) and recently released by the Geospatial Information Authority of Japan. The GPS data are processed by using the BERNESE 5.0 software and the strategy described by Ozawa et al. (2011). Because the data format, its technical specifications, precision and accuracy, and the GPS network itself have been detailed earlier, herein we provide only a brief summary. In addition, the post-seismic displacements used in this study are relative to a reference site (34.77°N, 131.92°E), which is ∼1400 km away from the earthquake epicentre. This reference site was selected because it is stable and far enough from the rupture fault, which should not significantly affected by the earthquake and other local deformations.

In the GPS data, we expect overlapping deformation processes. Given that post-seismic deformation and interseismic signals are captured by post-seismic GPS observations, we first subtract interseismic signals from the observations. For this we use the interseismic velocity fields solved from GPS observations before the event (Hatanaka et al. 2003; Loveless & Meade 2010). Assuming that the interseismic velocity field has remained constant, we can remove the interseismic signals from the post-seismic displacement records. In addition, the displacements induced by strong aftershocks (M > 7) are modelled and removed from the observation in each period using the method of Diao et al. (2011b).

This leaves us with a corrected post-seismic displacement data set from eight periods, still likely containing the signals of overlapping post-seismic processes. The observed post-seismic displacements vectors are found most pronounced between latitude 35°N and 41°N, and are oriented generally perpendicular to the trench. In other words, the spatial pattern of post-seismic displacement resembles that of the coseismic displacement field. To allow comparison of the time-dependent characteristics of the post-seismic displacement, we transfer the displacement in each period to the daily crustal velocity.

GPS velocities show a decreasing displacement rate in time (see ‘observation’ in Fig. 2). The maximum velocity is approximately 2 cm d−1 in the initial period of (0–7 d) and decreases to 0.1 cm d−1 at the last period of (382–564 d). The maximum of the accumulative horizontal post-seismic displacement within the 564-d period is 1.05 m, which is almost 20 per cent of the maximum coseismic displacement recorded onshore. The following sections will address the processes that may explain this significant displacement.

Time-dependent post-seismic deformations following the 2011 Mw 9.0 Toholu earthquake, in which the interseismic deformation has been removed. The first column is the observed GPS displacement, the second column is the prediction of the afterslip model (Model 1) and the last column is the residual of the data fitting.
Figure 2.

Time-dependent post-seismic deformations following the 2011 Mw 9.0 Toholu earthquake, in which the interseismic deformation has been removed. The first column is the observed GPS displacement, the second column is the prediction of the afterslip model (Model 1) and the last column is the residual of the data fitting.

3 MODELLING METHODS AND RESULTS

We consider different post-seismic processes in modelling simulations. All make use of an inversion routine, which aim is to minimize the misfit between the model and data. In specific, we explore two different inversion modelling strategies: (i) we simulate pure afterslip only and (ii) we simulate the combined model of afterslip and viscoelastic relaxation. In the following we detail the used methods, the estimated results and a resolution test that was performed to better judge and challenge the validity of the afterslip model and the combined model.

3.1 Afterslip model (Model 1)

In this first model, we simulate post-seismic dislocation along the megathrust fault only. The assumption that all post-seismic deformations are induced by afterslip, which is similar to the assumption made by Ozawa et al. (2012).

Because afterslip is known to occur often in areas surrounding the rupture of a main shock (e.g. Wang et al. 2009; Ozawa et al. 2012), we herein extend the coseismic fault surface proposed by Diao et al. (2012) in all strike and dip directions. For the inversion of the afterslip, we accordingly construct a fault plane with a spatial scale of 1000 km × 500 km (Fig. 1), with a curvature that approximates the curvature of the trench. The fault dip increases with depth from 6° to 45° according to the slab isodepth (Miura et al. 2005). The fault geometry is subdivided into rectangular fault patches with a regular size of 30 km × 30 km to obtain a distribution of afterslip.

The inversion code, SDM, developed by Wang et al. (2013) is used to derive the afterslip distribution. The name SDM stands for steepest descent method, an iterative algorithm used for the constrained least-squares optimization. The physical constraint is introduced to get a smooth slip model, which is realized through a roughness term to be minimized with the misfit to data. In difference to most previous inversion codes, the smoothing can be optionally applied to stress-drop on the whole fault. Actually, a constant or fairly smooth stress drop within slip asperities is usually assumed in the seismology. The objective function can be written as
(1)
where τ is the shear stress drop caused by the distributed slip on the whole fault plane, H is the finite difference approximation of the Laplacian operator multiplied by a weighting factor relating to the slip amplitude and α is the smoothing factor. b is the slip on each patch to be solved, y is the observed GPS displacements and G is the Green's function. The weighting factor is imported to the Laplacian operator to resolve the slip concentrations effectively. We choose the smoothing factor based on the trade-off curve between the stress roughness of the model and the fitting residual. Green's functions are calculated by using a layered earth model, which is explained in the subsequent section.

Fig. 3(a) shows the estimated cumulative afterslip 564 d following the main shock. The maximum afterslip of 3.8 m is mainly distributed at the downdip area of the main rupture (between 20 and 80 km), which slips generally perpendicular to the trench. The seismic moment released by the afterslip reaches up to 2.3 × 1022 Nm, which is approximately 60 per cent of the moment released by the main shock and corresponds to an Mw 8.84 earthquake. The distribution of the aftershock 564 d after the main shock is generally consistent with the distribution of the afterslip.

The fault slip rates in each period are calculated to infer the time-dependent characteristics of the afterslip. The maximum slip rate decreases rapidly from 11.0 cm d−1 in the initial period of (0–7 d) to 0.5 cm d−1 in the last period of (382–564 d), but the spatial distribution of the afterslip remains stable in each period (Fig. 4). The afterslip model can match the observations at a high level (Fig. 2). The horizontal root mean square (rms) residuals in each period are less than 5 mm, which is close to the uncertainty of the GPS observation. However, the rms residuals of the vertical component are generally close to 1 cm, which is significantly larger than that of the horizontal component.

Distribution of afterslip in 564 d following the main shock. (a) is the result derived from Model 1 and (b) is the result derived from Model 2. The black contours are the coseismic slip with an interval of 5 m. The grey contours are the depth of the fault plain with an interval of 20 km. The grey dots are aftershocks in this period with M ≥ 5.0 (aftershock catalogue is obtained from the USGS).
Figure 3.

Distribution of afterslip in 564 d following the main shock. (a) is the result derived from Model 1 and (b) is the result derived from Model 2. The black contours are the coseismic slip with an interval of 5 m. The grey contours are the depth of the fault plain with an interval of 20 km. The grey dots are aftershocks in this period with M ≥ 5.0 (aftershock catalogue is obtained from the USGS).

Time-dependent afterslip inverted from the afterslip model. Black star denotes the epicentre of the main shock. The grey contours are the coseismic slip with an interval of 5 m. The afterslip rate was drawn with a non-linear colour scale.
Figure 4.

Time-dependent afterslip inverted from the afterslip model. Black star denotes the epicentre of the main shock. The grey contours are the coseismic slip with an interval of 5 m. The afterslip rate was drawn with a non-linear colour scale.

3.2 Combined model (Model 2)

Although the afterslip model (Model 1) can explain the observation well, we are not certain if the afterslip is the only mechanism that controls the post-seismic process because a large number of free parameters (fault patches) are used in this model.

Viscoelastic relaxation in the medium surrounding a ruptured fault has found to be evident following other megathrust earthquakes (e.g. Lorenzo-Martín et al. 2006; Pollitz et al. 2008; Hoechner et al. 2012; Hu & Wang 2012). Several studies have shown that post-seismic displacement can be explained more reasonably based on combined models of afterslip and viscoelastic relaxation (Sheu & Shieh 2004; Bruhat et al. 2011; Diao et al. 2011a; Hoechner et al. 2012; Hu & Wang 2012). Therefore, we construct a combined model by considering the effects of afterslip and viscoelastic relaxation simultaneously (Model 2).

Our viscoelastic earth model includes four layers (Fig. 5): an elastic upper crust with thickness E as a free parameter and three viscoelastic layers beneath the elastic crust, representing the asthenosphere, upper and lower mantle, for all of which the Maxwell rheology is used. The asthenosphere viscosity is defined as another free parameter η, but the upper- and lower-mantle viscosities are fixed at 1 × 1020 and 1 × 1021 Pa s, respectively. The crust structure is collected from the Crust 2.0 model (Mooney et al. 1998), whereas the elasticity parameters for the medium under the crust are adopted from the Preliminary Reference Earth Model (PREM; Dziewonski & Anderson 1981). This layered earth model is considered reasonable for modelling the viscoelastic relaxation induced by large subduction zone earthquakes (Pollitz et al. 2006; Han et al. 2008).

Layered earth model and medium parameters used in our study.
Figure 5.

Layered earth model and medium parameters used in our study.

The modelling of the viscoelastic relaxation is carried out by using the PSGRN/PSCMP code developed by Wang et al. (2006) based on the viscoelastic-gravitational dislocation theory. The coseismic slip distribution is adopted from Diao et al. (2012). The grid search method is used to find the optimal values of the two parameters, E (crust thickness) and η (asthenosphere viscosity). Considering that the slip of the main shock mostly occurs at depths less than 35 km, we set the lower boundary of the thickness of the elastic layer at 35 km.

The grid search is performed in two steps. (1) For each given combination of E and η (gridpoint), the viscoelastic part of the post-seismic deformation is calculated and removed from the observed data. (2) The afterslip in each observation period is inverted from the residual data from (1), and the final rms misfit between the observation and the combined model is calculated by
(2)
where Obs is the GPS displacement data in both horizontal and vertical direction; Calvis is the displacement calculated from the viscoelastic model, which is a non-linear function of the elastic thickness (E) and viscosity (η); Calafterslip is the displacement contributed by the afterslip; i is the index for the post-seismic periods and j is the index for the GPS stations. We consider the secondary viscoelastic relaxation caused by afterslip when modelling the viscoelastic relaxation effect in each period, that is, the coseismic slip and the accumulated afterslip obtained in Model 1 are introduced as the deriving source of the viscoelastic relaxation in following period.

The grid search results in the combined model are shown in Fig. 6(a). The optimal values of the elastic thickness and viscosity are 50 km and 2 × 1019 Pa s, respectively. The afterslip in the initial 564 d after the main shock obtained in Model 2 reveals similar slip patterns to that of Model 1 (Fig. 3b). However, the afterslip magnitude in Model 2 is slightly decreased compared with that of Model 1 because of the effect of viscoelastic relaxation. The seismic moment released by the afterslip in Model 2 is 2.1 × 1022 Nm, which is 93 per cent of the afterslip moment estimated by Model 1 and corresponds to an Mw 8.82 earthquake. The slip rate also shows rapid decay with time (Fig. 7), which is similar to that of Model 1.

Misfit variations in searching for the optimal viscosity and thickness of the elastic layer. (a) is the searching result for Model 2; (b) is the searching result for the input combined model.
Figure 6.

Misfit variations in searching for the optimal viscosity and thickness of the elastic layer. (a) is the searching result for Model 2; (b) is the searching result for the input combined model.

Time-dependent afterslip inverted from Model 2. Black star denotes the epicentre of the main shock. The grey contours are the coseismic slip with an interval of 5 m. The afterslip rate was draw with a non-linear colour scale.
Figure 7.

Time-dependent afterslip inverted from Model 2. Black star denotes the epicentre of the main shock. The grey contours are the coseismic slip with an interval of 5 m. The afterslip rate was draw with a non-linear colour scale.

The predicted and residual displacements obtained by using the optimal values of E = 50 km and η = 2 × 1019 Pa s are shown in Fig. 8. The combined model can explain the observation equally well as the ‘pure’ afterslip model. Post-seismic displacements are mainly induced by the afterslip in the initial periods. However, the contribution of the viscoelastic relaxation can be found increases with time and becomes larger in later periods (Fig. 8).

The first column is the same as that in the Fig. 2, the second column are the contribution of the fault afterslip in combined model (Model 2), the third column are the contribution of the viscoelastic relaxation in the combined model (Model 2) and the last column is the residual of the model fitting.
Figure 8.

The first column is the same as that in the Fig. 2, the second column are the contribution of the fault afterslip in combined model (Model 2), the third column are the contribution of the viscoelastic relaxation in the combined model (Model 2) and the last column is the residual of the model fitting.

3.3 Model resolution test

The checkerboard test indicates that the input slip model can generally be recovered with a spatial resolution of about 150 km by onshore GPS observations (Fig. 9b). In comparison, the slip close to the trench is less resolved than that in depth. However, as mentioned in previous paper (Diao et al. 2012), the resolution of the fault patches close to the trench can be significantly increased by adding the seafloor observation, which perhaps could be obtained in the future.

Checker board test of the data resolution of the data. (a) Input slip distribution; (b) Recovered slip distribution based on assumptions in Model 1 and (c) Recovered slip distribution from the combined model.
Figure 9.

Checker board test of the data resolution of the data. (a) Input slip distribution; (b) Recovered slip distribution based on assumptions in Model 1 and (c) Recovered slip distribution from the combined model.

To validate the reliability of the combined model, we first construct a combined model that includes fault slip and viscoelastic relaxation and generate the synthetic GPS data based on this model. We then invert the synthetic data using the method mentioned in Section 3.2 to test whether the input slip distribution and parameters of the viscoelastic model can be recovered. In the input model, we set the thickness of the elastic layer at 50 km and the viscosity of the asthenosphere layer at 5 × 1018 Pa s and added a slip model shown in Fig. 9(a). The grid-searched optimal values are 50 km and 4 × 1018 Pa s (Fig. 6b), respectively, very close to the input values. The afterslip can also be recovered well simultaneously (Fig. 9c). This test indicates that the setting and results of the combined model are reliable. In addition, the misfits show larger gradient at low viscosities, thus the lower bound to the apparent viscosity is stronger than the upper bound.

4 DISCUSSION

Determination of the relevance of the different processes responsible for post-seismic deformation following large earthquakes has remained challenging. Several studies consider the afterslip to be the governing mechanism of post-seismic deformations, and neglect the effect of viscoelastic relaxation in the initial years following megathrust earthquakes (e.g. Hsu et al. 2006; Ozawa et al. 2012). However, other studies indicate that post-seismic displacements can be similarly well simulated by using viscoelastic relaxation processes only (e.g. Pollitz et al. 2006; Kogan et al. 2011). In this study, we add evidence that both, afterslip and viscoelastic relaxation, may jointly occur after larger earthquakes, but play different roles in different periods. In the spatial domain, post-seismic deformation predicted by afterslip generally show short-wavelength character compared with that predicted by viscoelastic relaxation. Thus, afterslip and viscoelastic relaxation may control the near-field and far-field post-seismic deformation, respectively. Recent studies (Wang et al. 2009; Bruhat et al. 2011; Diao et al. 2011a; Wang et al. 2012), including this one, consider the combined effect of afterslip and viscoelastic relaxation when analysing post-seismic displacements, offering a reasonable way for further post-seismic research.

4.1 Temporal and spatial evolution of fault afterslip

In the case of the Tohoku earthquake, the afterslip is mainly distributed in the downdip areas of the hypocentre and remains basically stationary throughout the 564-d period, whereas the coseismic slip is mostly located at the updip areas of the hypocentre (Fig. 3). The overall distribution of afterslip obtained in this paper show a general consistence with the spatial distribution of afterslip inverted by Ozawa et al. (2012), who used the post-seismic displacement data within the first 7 months after the earthquake.

Partial overlapping between the coseismic slip and the post-seismic afterslip can be found near the hypocentre, where two significant earthquakes with M > 7 occurred in the 20th century (Fig. 1). Moreover, a clear afterslip cluster is found around the rupture of the Mw 7.9 aftershock (141°E, 35.5°N). In addition, we find that the afterslip rate decays rapidly with time. Our analysis reveals that this decay is robust, being independent on the viscoelastic relaxation effect (Figs 4 and 7). The moment released by the aftershocks (M > 5) in this time span is 3.1 × 1020 Nm, just ∼1.5 per cent of the moment released by the afterslip in Model 2. Thus, the afterslip mainly occurs aseismically.

Post-seismic displacements in both horizontal and vertical directions are used to investigate the post-seismic process of the Tohoku earthquake but the vertical component was set with relatively small weight in afterslip inversions due to the large uncertainty (Ozawa et al. 2012). In case of the Tohoku earthquake, we found that the vertical components only play a second role compared with the horizontal components. As the afterslip model estimated by using both horizontal and vertical displacements does not show clear difference compared with the afterslip model estimated just using the horizontal component.

4.2 Regional rheological structure of the Japan subduction zone

The optimal viscosity of the asthenosphere obtained from the combined model is 2 × 1019 Pa s, which resembles the steady-state viscosity of 1019 Pa s obtained from the post-seismic studies of other megathrust earthquakes (e.g. Hoechner et al. 2012; Hu & Wang 2012). The estimated elastic thickness of 50 km is consistent with the depth range where 80 per cent of the aftershocks have occurred and the major coseismic rupture has been located. It should be noted that we considered the secondary viscoelastic relaxation caused by afterslip in the inversion (we call it the secondary viscoelastic relaxation as the viscoelastic relaxation caused by the coseismic slip plays a dominant role). This enlarges the driving source and results in a larger viscosity compared with the viscosity of 8 × 1018 Pa obtained by neglecting the secondary effect.

Although the viscosity of 2 × 1019 Pa s obtained in this paper in dimension agrees well with other viscosity estimates near the northeast Japan arc (Suito & Hirahara 1999; Muto et al. 2013), it is larger than values proposed by Kogan et al. (2011). In the latter, a Maxwell viscosity of 5–10 × 1017 Pa s is proposed based on a relatively small number of GPS observations following the 2006–2007 Kuril earthquakes. The viscosity obtained in this paper is also larger than that obtained by Ohzono et al. (2012), who proposed 2.4–4.8 × 1018 Pa s by using 2-yr continuous GPS data following the 2008 Iwate–Miyagi Nairiku earthquake. One may speculate that a reason for their small viscosity is that they model the post-seismic displacements using the viscoelastic relaxation only, whereas the fault afterslip effect was not considered. Another reason might be related to the timescale of the considered viscoelastic relaxation process, which depends on the earthquake magnitude as proposed by Wang et al. (2012). We therefore emphasize that viscoelastic relaxation is a long-term process, and a reliable recognition of viscosity and its transient characters require more observations.

To consider other modelling complexities, we made tests using numerical codes. In specific we test how the crustal deformation can be affected by inserting an elastic slab in the asthenosphere. To this aim a set of 2-D finite element models were constructed using the Pylith software (Aagaard et al. 2013). The post-seismic displacement induced by the viscoelastic relaxation is found to be clearly reduced by an elastic slab inserted into the layered model (Fig. 10a). These results are also in consistence with the findings of Pollitz et al. (2008) and Tanaka et al. (2009). However, the effect of inserting the slab is negligible when calculating the displacements induced by fault slip (Fig. 10b). Thus, by inserting the slab cannot affect the dominant role of the afterslip in the short-term post-seismic deformation process. Nevertheless, neglecting the effect of the slab may cause the viscosity to be slightly overestimated. Therefore, we suggest that more rigorous methods, such as Finite Element Modeling (FEM) method including 3-D rheological structure of the Japan subduction zone, may further improve our understanding of the long-term viscoelastic relaxation processes to come.

Sketch maps of the 2-D FEM models (left-hand side) and comparisons of the surface displacements predicted by models A and B (right-hand side). (a) is the comparison of displacements predicted by 2-yr viscoelastic relaxation of models A and B and (b) is the comparison of displacements predicted by fault slip of models A and B. The 2-D coseismic slip distribution was used as the source.
Figure 10.

Sketch maps of the 2-D FEM models (left-hand side) and comparisons of the surface displacements predicted by models A and B (right-hand side). (a) is the comparison of displacements predicted by 2-yr viscoelastic relaxation of models A and B and (b) is the comparison of displacements predicted by fault slip of models A and B. The 2-D coseismic slip distribution was used as the source.

4.3 Seismic moment released by afterslip and viscoelastic relaxation

Our model results suggest an afterslip moment of 2.3 × 1022 Nm that has accumulated in 564 d after the main shock, which is about 60 per cent of the moment released by the coseismic rupture and corresponds to an Mw 8.84 earthquake. The afterslip moment in the initial 7 months obtained by Ozawa et al. (2012) is 18 per cent of coseismic moment with a half-space model, which is clearly smaller compared with our result by 43 per cent in the same period obtained with a layered Earth structure. Time-dependent seismic moments released by the fault afterslip in the afterslip model and the combined model are shown in Fig. 11(a). In the initial periods (<0.2 yr), the contribution of the viscoelastic relaxation is relatively small compared with that of the afterslip (<5 per cent). However, the effect of the viscoelastic relaxation increases with time in the following periods. In the last period (i.e. 382–564 d), the moments induced by the viscoelastic relaxation can reach up to 11 per cent of the total moment. Moreover, the relative contributions of the afterslip and viscoelastic relaxation can also be inferred from the displacements predicted by the afterslip and viscoelastic relaxation in the combined model (Fig. 8). The time-dependent seismic moment in the combined model suggests that viscoelastic relaxation can only be neglected in the initial period after the earthquake and may need to be considered in the following periods.

Estimated time-dependent moment released by afterslip and viscoelastic relaxation following the Tohoku earthquake. In (a) the secondary relaxation effect induced by afterslip is considered, while in (b) this effect is neglected. Black and grey dots are seismic moments derived from the ‘pure’ afterslip model (Model 1) and afterslip model in the combined model (Model 2). The blue colour shows that contribution of afterslip in the combined model, whereas the orange colour is the difference between the moment that released by afterslip in Models 1 and 2, which is considered to be the moment released by viscoelastic relaxation in the combined model. The black and orange numbers are the percentage of the moment that released by afterslip and viscoelastic relaxation in the combined model for each period.
Figure 11.

Estimated time-dependent moment released by afterslip and viscoelastic relaxation following the Tohoku earthquake. In (a) the secondary relaxation effect induced by afterslip is considered, while in (b) this effect is neglected. Black and grey dots are seismic moments derived from the ‘pure’ afterslip model (Model 1) and afterslip model in the combined model (Model 2). The blue colour shows that contribution of afterslip in the combined model, whereas the orange colour is the difference between the moment that released by afterslip in Models 1 and 2, which is considered to be the moment released by viscoelastic relaxation in the combined model. The black and orange numbers are the percentage of the moment that released by afterslip and viscoelastic relaxation in the combined model for each period.

In addition, it should be noted that the secondary viscoelastic effect induced by afterslip is considered in our combined model. If this secondary effect was neglected, the inverted Maxwell viscosity and elastic thickness are 8 × 1018 Pa s and 40 km, and the contribution of viscoelastic relaxation will be significantly overestimated (Fig. 11b). Therefore, it is necessary to consider the secondary relaxation effect induced by afterslip when analysing long-term post-seismic processes.

4.4 Coulomb failure stress (CFS) changes induced by coseismic and post-seismic fault slip

CFS change induced by coseismic and post-seismic slip may directly relate to further seismic activities around the rupture area (Toda et al. 2011), thus investigation of the CFS change induced by fault slip of the Tohoku earthquake might be relevant for further hazard assessment along the Japan trench. Fig. 12 shows the CFS change induced by the coseismic rupture, afterslip and combined effect of them. The CFS changes are calculated by ΔCFS = Δτ − μ′Δσn, where ΔCFS is the CFS change, Δτ is the maximum shear stress change on fault patch, μ′ is the effective friction coefficient including the effect of pore-fluid change and Δσn is the normal stress change. A wide range of μ′ is considered, but the obtained CFS change does not change significantly. We use μ′ = 0.5 in this paper.

Coulomb stress change (colour) on the plate interface induced by the coseismic rupture (a), fault afterslip (b) and combined effect of them (c). The grey contours are afterslip in Model 2 (0.5 m intervals). The grey dots are aftershocks in this period with M ≥ 5.0 (aftershock catalogue is obtained from the USGS). Two black stars denote the epicentre of the main shock and the Mw 7.9 aftershock. The four green squares are selected fault patches.
Figure 12.

Coulomb stress change (colour) on the plate interface induced by the coseismic rupture (a), fault afterslip (b) and combined effect of them (c). The grey contours are afterslip in Model 2 (0.5 m intervals). The grey dots are aftershocks in this period with M ≥ 5.0 (aftershock catalogue is obtained from the USGS). Two black stars denote the epicentre of the main shock and the Mw 7.9 aftershock. The four green squares are selected fault patches.

The maximum CFS drop of ∼10 MPa is found within the major coseismic slip asperity. Considerable increases of CFS appear in areas surrounding the slip asperities and show a good spatial correlation with the aftershock and afterslip activities (Fig. 12a). The increased CFS at downdip area of the coseismic rupture was partly unloaded by afterslip in 564���d after the event (Fig. 12b). However, a large part of the CFS increments persists, which was (and/or in the future will be) released by long-term afterslip and viscoelastic relaxation or transient rupture (Fig. 12c). The CFS on shallow patches at the two ends of the main rupture were increased by co- and post-seismic fault slip, indicating that the seismic hazard in the two areas has been enlarged by the Tohoku event. However, the CFS was just slightly changed near Tokyo at the southern part of the fault.

To investigate the frictional behaviour of the fault, the slip rate and CFS change as functions of time are shown in Fig. 13 for four selected locations (Fig. 12a). The CFS on points A, C and D are clearly increased by the coseismic rupture, which is partially released by time-dependent afterslip, whereas the CFS on point B is significantly decreased by coseismic slip and is slightly increased by followed afterslip (Fig. 13b). In addition, the CFS decrease almost linearly with the decrease of logarithm of slip rate on points A, C and D, which show general consistence with the behaviour of state- and rate-strengthening frictional law as indicated by Miyazaki et al. (2004) and Hsu et al. (2006). However, point B shows different behaviour that could not be explained by the state- and rate-strengthening frictional law. The slopes of d(C F S)/d log(slip - rate) = Aσn are obtained by using the least-square fittings on points A, C and D, where σnis the effective normal stress, A is a positive rheological parameter representing the fault frictional behaviour. For depths of 48, 21 and 70 km, the effective normal stress would be 1290, 560 and 1900 MPa, and the derived values of A are 2.4 × 10−4, 1.6 × 10−4 and 0.4 × 10−4, which are consistent with previous findings (Miyazaki et al. 2004; Hsu et al. 2006). We also find that the A value on deep patch (point D) is clearly smaller than that of shallow patches, which perhaps indicates the spatial discrepancy of the fault slip behaviour.

Slip rate history (a), Coulomb stress change evolution (b) and CFS change in log (slip rate) phase space on selected fault patch (Fig. 12a). The lines in (c) was obtained by fitting the observation (black dots) using the least-square method.
Figure 13.

Slip rate history (a), Coulomb stress change evolution (b) and CFS change in log (slip rate) phase space on selected fault patch (Fig. 12a). The lines in (c) was obtained by fitting the observation (black dots) using the least-square method.

5 CONCLUSIONS

We use the time-dependent post-seismic displacements in the initial ∼1.5 yr (564 d) following the Tohoku earthquake recorded by the onshore GPS stations to investigate the involved physical mechanism. The following conclusions are drawn:

  1. Both the afterslip and combined models can well explain the GPS observations and indicate significant afterslip in the downdip area of the coseismic rupture. The seismic moments released by the aseismic slip in the two models are 2.3 × 1022 and 2.1 × 1022 Nm, which are equivalent to Mw 8.84 and 8.82 earthquakes, respectively. The effect of the viscoelastic relaxation within this short term only plays a secondary role, but it shows an increasing tendency.

  2. The afterslip and aftershock are mostly distributed in the areas with increased Coulomb stress, thus indicating that the afterslip and aftershocks are strongly related to the rupture of the main shock. Afterslip in 564 d releases only a small part of the Coulomb stress increasing downdip of the coseismic rupture area.

  3. The optimal viscosity of the asthenosphere in the combined model is found to be 2 × 1019 Pa s by using a grid search method. This finding agrees well with the viscosity of 1019 Pa s obtained from previous post-seismic studies of large thrust earthquakes.

  4. In both the afterslip and combined models, the slip rate on patches downdip of the coseismic rupture decreases rapidly with time and obeys the state- and rate-strengthening frictional law.

We thank Dr T. Nishimura, who provided the GPS data used in this paper and gave valuable suggestions for improvement of the initial manuscript. We thank the Editor Dr Kosuke Heki and two anonymous reviewers for their insightful comments. This work is supported by the National Natural Science Foundation of China (Grant Nos. 41304017, 41074052, 41174086, 41021003), Key Development Program of the Institute of Geodesy and Geophysics, CAS (Y309411041) the International Partnership Program for Creative Research Teams of Chinese Academy and Science/State Administration of Foreign Experts Affairs of Chinese Academy of Sciences (CAS/SAFEA) (Grant No. KZZD-EW-TZ-05), the Open Research Fund Program of the Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, China (Grant No. 12-01-02). The figures were generated by using the open-source Generic Mapping Tools software (Wessel & Smith 1991).

REFERENCES

Aagaard
B.
Kientz
S.
Knepley
M.
Strand
L.
Williams
C.
PyLith User Manual
2013
 
Available at: (last accessed 11 October 2013) Version 1.9.0, Computational Infrastructure of Geodynamics. geodynamics.org/cig/software/pylith/pylith_manual-1.9.0.pdf/view
Bruhat
L.
Barbot
S.
Avouac
J.-P.
Evidence for postseismic deformation of the lower crust following the 2004 Mw 6.0 Parkfield earthquake
J. geophys. Res.
2011
, vol. 
116
 pg. 
B08401
  
doi:10.1029/2010JB008073
Bürgmann
R.
Dresen
G.
Rheology of the lower crust and upper mantle: evidence from rock mechanics, geodesy, and field observations
Ann. Rev. Earth planet. Sci.
2008
, vol. 
36
 (pg. 
531
-
567
)
DeMets
C.
Oblique convergence and deformation along the Kuril and Japan trenches
J. geophys. Res.
1992
, vol. 
97
 
B12
(pg. 
17 615
-
17 625
)
Diao
F.
Xiong
X.
Wang
R.
Mechanism of transient post-seismic deformation following the 2001 Mw 7.8 Kunlun (China) earthquake
Pure appl. Geophys.
2011a
, vol. 
168
 (pg. 
767
-
779
)
Diao
F.
Xiong
X.
Ni
S.
Zheng
Y.
Ge
C.
Slip model for the 2011 Mw 9.0 Sendai (Japan) earthquake and its Mw 7.9 aftershock derived from GPS data
Chin. Sci. Bull.
2011b
, vol. 
56
 (pg. 
1999
-
2005
)
Diao
F.
Xiong
X.
Zheng
Y.
Static slip model of the Mw 9.0 Tohoku (Japan) earthquake: results from joint inversion of terrestrial GPS data and seafloor GPS/acoustic data
Chin. Sci. Bull.
2012
, vol. 
57
 
16
(pg. 
1990
-
1997
)
Dziewonski
A.M.
Anderson
D.L.
Preliminary Reference Earth Model
Phys. Earth planet. Inter.
1981
, vol. 
25
 (pg. 
297
-
356
)
Freed
A.M.
Afterslip (and only afterslip) following the 2004 Parkfield, California earthquake
Geophys. Res. Lett.
2007
, vol. 
34
 
6
pg. 
L06312
  
doi:10.1029/2006GL029155
Han
S.-C.
Sauber
J.
Luthcke
S.B.
Ji
C.
Pollitz
F.F.
Implications of postseismic gravity change following the great 2004 Sumatra-Andaman earthquake from the regional harmonic analysis of GRACE intersatellite tracking data
J. geophys. Res.
2008
, vol. 
113
 pg. 
B11413
  
doi:10.1029/2008JB005705
Hashimoto
C.
Noda
A.
Sagiya
T.
Matsu'ura
M.
Interplate seismogenic zones along the Kuril–Japan trench inferred from GPS data inversion
Nat. Geosci.
2009
, vol. 
2
 (pg. 
141
-
144
)
Hatanaka
Y.
Iizuka
T.
Sawada
M.
Yamagiwa
A.
Kikuta
Y.
Johnson
J.M.
Rocken
C.
Improvement of the analysis strategy of GEONET
Bull. Geogr. Surv. Inst.
2003
, vol. 
49
 (pg. 
11
-
37
)
Hoechner
A.
Sobolev
S.V.
Einarsson
I.
Wang
R.
Investigation on afterslip and steady state and transient rheology based on postseismic deformation and geoid change caused by the Sumatra 2004 earthquake
Geochem. Geophys. Geosyst.
2012
, vol. 
12
 pg. 
Q07010
  
doi:10.1029/2010GC003450
Hsu
Y.J.
et al. 
Frictional afterslip following the Mw 8.7, 2005 Nias-Simeulue earthquake, Sumatra
Science
2006
, vol. 
312
 (pg. 
1921
-
1926
)
Hu
Y.
Wang
K.
Spherical-Earth finite element model of short-term postseismic deformation following the 2004 Sumatra earthquake
J. geophys. Res.
2012
, vol. 
117
 pg. 
B05404
  
doi:10.1029/2012JB009153
Jónsson
S.
Segall
P.
Pedersen
R.
Björnsson
G.
Post-earthquake ground movements correlated to poro-pressure transients
Nature
2003
, vol. 
424
 (pg. 
179
-
183
)
Kogan
M.G.
Vasilenko
N.F.
Frolov
D.I.
Freymueller
J.T.
Steblov
G.M.
Levin
B.W.
Prytkov
A.S.
The mechanism of postseismic deformation triggered by the 2006–2007 great Kuril earthquakes
Geophys. Res. Lett.
2011
, vol. 
38
 pg. 
L06304
  
doi:10.1029/2011GL046855
Lorenzo-Martín
F.
Roth
F.
Wang
R.
Inversion for rheological parameters from postseismic surface deformation associated to the 1960 Valdivia earthquake, Chile
Geophys. J. Int.
2006
, vol. 
164
 
6
(pg. 
75
-
87
)
Loveless
J.P.
Meade
B.J.
Geodetic imaging of plate motions, slip rates, and partitioning of deformation in Japan
J. geophys. Res.
2010
, vol. 
115
 pg. 
B02410
  
doi:10.1029/2008JB006248
Miura
S.
Takahashi
N.
Nakanishi
A.
Tsuru
T.
Kodaira
S.
Kaneda
Y.
Structural characteristics off Miyagi forearc region, the Japan Trench seismogenic zone, deduced from a wide-angle and refraction study
Tectonophysics
2005
, vol. 
407
 (pg. 
165
-
188
)
Miyazaki
S.
Segall
P.
Fukuda
J.
Kato
T.
Space time distribution of afterslip following the 2003 Tokachi-oki earthquake: implications for variations in fault zone frictional properties
Geophys. Res. Lett.
2004
, vol. 
31
 pg. 
L06623
  
doi:10.1029/2003GL019410
Mooney
W.D.
Laske
G.
Masters
T.G.
CRUST 5.1: a global crustal model at 5° × 5°
J. geophys. Res.
1998
, vol. 
103
 (pg. 
727
-
747
)
Muto
J.
Shibazaki
B.
Ito
Y.
Iinuma
T.
Ohzono
M.
Matsumoto
T.
Okada
T.
Two-dimensional viscosity structure of the northeastern Japan islands arc-trench system
Geophys. Res. Lett.
2013
, vol. 
40
  
doi:10.1002/grl.50906
Ohzono
M.
Ohta
Y.
Iinuma
T.
Miura
S.
Muto
J.
Geodetic evidence of viscoelastic relaxation after the 2008 Iwate-Miyagi Nairiku earthquake
Earth Planets Space
2012
, vol. 
64
 
9
(pg. 
759
-
764
)
Ozawa
S.
Nishinura
T.
Suito
H.
Kobayashi
T.
Tobita
M.
Iakiire
T.
Coseismic and postseismic slip of the 2011 magnitude-9 Tohoku-Oki earthquake
Nature
2011
, vol. 
475
 (pg. 
373
-
376
)
Ozawa
S.
Nishimura
T.
Munekane
H.
Suito
H.
Kobayashi
T.
Tobita
M.
Imakiire
T.
Preceding, coseismic, and postseismic slips of the 2011 Tohoku earthquake, Japan
J. geophys. Res.
2012
, vol. 
117
 pg. 
B07404
  
doi:10.1029/2011JB009120
Pollitz
F.
Banerjee
P.
Grijalva
K.
Nagarajan
B.
Bürgmann
R.
Effect of 3-D viscoelastic structure on post-seismic relaxation from the 2004 M = 9.2 Sumatra earthquake
Geophys. J. Int.
2008
, vol. 
173
 
1
(pg. 
189
-
204
)
Pollitz
F.F.
Bürgmann
R.
Banerjee
P.
Post-seismic relaxation following the great 2004 Sumatra-Andaman earthquake on a compressible self-gravitating Earth
Geophys. J. Int.
2006
, vol. 
167
 (pg. 
397
-
420
)
Segall
P.
Davis
J.L.
GPS applications for geodynamics and earthquake studies
Ann. Rev. Earth planet. Sci.
1997
, vol. 
25
 (pg. 
301
-
336
)
Shen
Z.K.
et al. 
Slip maxima at fault junctions and rupturing of barriers during the 2008 Wenchuan earthquake
Nat. Geosci.
2009
, vol. 
2
 
10
(pg. 
718
-
724
)
Sheu
S.Y.
Shieh
C.F.
Viscoelastic-afterslip concurrence: a possible mechanism in the early post-seismic deformation of the Mw 7.6, Chi-Chi (Taiwan) earthquake
Geophys. J. Int.
2004
, vol. 
159
 (pg. 
1112
-
1124
)
Simons
M.
et al. 
The 2011 magnitude 9.0 Tohoku-Oki Earthquake: mosaicking the megathrust from seconds to centuries
Science
2011
, vol. 
332
 (pg. 
1421
-
1425
)
Suito
H.
Hirahara
K.
Simulation of postseismic deformation caused by the 1896 Riku-u earthquake, northeast Japan: re-evaluation of the viscosity in the upper mantle
Geophys. Res. Lett.
1999
, vol. 
26
 (pg. 
2561
-
2564
)
Tanaka
Y.
Klemann
V.
Fleming
K.
Martinec
Z.
Spectral finite element approach to postseismic deformation in a viscoelastic self-gravitating spherical Earth
Geophys. J. Int.
2009
, vol. 
176
 (pg. 
715
-
739
)
Thatcher
W.
Rundle
J.B.
A model for the earthquake cycle in underthrust zones
J. geophys. Res.
1979
, vol. 
84
 
B10
(pg. 
5540
-
5556
)
Toda
S.
Lin
J.
Stein
S.
Using the 2011 Mw 9.0 off the Pacific coast of Tohoku earthquake to test the Coulomb stress triggering hypothesis and to calculate faults brought closer to failure
Earth Planets Space
2011
, vol. 
63
 
7
(pg. 
725
-
730
)
Wang
K.
Hu
Y.
He
J.
Deformation cycles of subduction earthquakes in a viscoelastic Earth
Nature
2012
, vol. 
484
 (pg. 
327
-
332
)
Wang
L.
Wang
R.
Roth
F.
Enescu
F.
Hainzl
B.S.
Ergintav
S.
After-slip and viscoelastic relaxation following the 1999 M7.4 Ízmit earthquake, from GPS measurements
Geophys J. Int.
2009
, vol. 
178
 (pg. 
1220
-
1237
)
Wang
R.
Lorenzo-Martín
F.
Roth
F.
PSGRN/PSCMP—a new code for calculating co- and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory
Comput. Geosci.
2006
, vol. 
32
 (pg. 
527
-
541
)
Wang
R.
Diao
F.
Hoechner
A.
SDM—a geodetic inversion code incorporating with layered crust structure and curved fault geometry
Proceedings of the EGU General Assembly 2013
2013
, vol. 
15
  
EGU2013-2411-1
Wessel
P.
Smith
W.H.F.
Free software helps map and display data
EOS, Trans. Am. geophys. Un.
1991
, vol. 
72
 (pg. 
445
-
446
)