## Abstract

Modeling of postseismic deformation following great earthquakes has revealed the viscous structure of the mantle and the frictional properties of the fault interface. However, for giant megathrust events, viscoelastic flow and afterslip mechanically interplay with each other during the postseismic period. We explore the role of afterslip and viscoelastic relaxation and their interaction in the aftermath of the 2011 *M*_{w} (moment magnitude) 9.0 Tohoku earthquake based on a detailed model analysis of the postseismic deformation with laterally varying, experimentally constrained, rock rheology. Mechanical coupling between viscoelastic relaxation and afterslip notably modifies both the afterslip distribution and surface deformation. Thus, we highlight the importance of addressing mechanical coupling for long-term studies of postseismic relaxation, especially in the context of the geodynamics of the Japan trench across the seismic cycle.

## INTRODUCTION

The postseismic deformation due to great and giant earthquakes provides useful constraints on the rheological properties of rocks in subduction zones (*1*–*7*). Laboratory rock deformation experiments illustrate a stress dependency to the viscosity of rocks (*8*, *9*), in addition to temperature, water content, and pressure dependencies. The transient and nonlinear viscous response to coseismic stress change implies a substantial reduction in the effective viscosity during the postseismic period (*3*, *4*, *10*–*16*). The coseismic stress changes along the plate boundary fault also drive afterslip along the plate boundary fault (*6*, *17*–*21*), which is thought to occur on velocity-strengthening regions that surround velocity-weakening asperities, with a link to aftershock production (*22*–*25*). The early stages of postseismic deformation result from a combination of transient nonlinear deformation of mantle rocks (*3*, *4*, *11*, *16*) and stress-driven afterslip on the fault (*6*, *17*–*21*). Thus, numerical models of quasi-static deformation may provide estimates of the fault extent and viscosity structure in the subsurface.

After the 2011 *M*_{w} (moment magnitude) 9.0 Tohoku earthquake, seismogeodetic observations and modeling have illuminated the short-term dynamics of the northeastern (NE) Japan island arc subduction zone (*1*, *6*, *16*, *26*–*33*). Often, these models have assumed no change in the rheological structure during the postseismic period and combine stress-driven relaxation based on a Maxwell or Burgers rheology with a kinematic inversion of afterslip to explain the postseismic displacement fields (*1*, *27*–*29*, *31*, *33*), implicitly assuming no significant stress interactions between fault slip and viscoelastic flow. However, as afterslip may accrue a substantial geodetic moment, mostly aseismic, and as viscoelastic flow redistributes stress during the postseismic period, these two processes may be interacting (i.e., be mechanically coupled) as they diffuse the sudden stress changes imposed by the seismic cycle. The goal of this study is to assess the importance of the coupling between viscoelastic flow and afterslip and whether it constitutes a critical aspect of the postseismic relaxation.

The deformation of crust and mantle rocks, constrained by laboratory experiments on aggregates or single crystals, is accommodated by several micromechanisms that include diffusion creep, dislocation creep, dislocation glide, and Peierls creep (*8*, *9*, *34*, *35*). During the steady state, i.e., at constant stress or strain rate, a constitutive relationship predicts the deformation rate based on stress and other in situ physical conditions. During periods of rapid change of stress or strain rate, other deformation processes take place within the rock, leading to a short phase of hardening that converges to the steady-state response, usually referred to as transient creep (*34*). In contrast to the steady state, the micromechanics of transient creep are poorly known and different constitutive relationships have been proposed (*4*, *36*, *37*), but plastic anisotropy within single crystals may be an important factor (*38*). Because of the rapid change imposed by the seismic cycle, transient creep may be operating during postseismic relaxation and is sometimes approximated with the linear, bi-viscous Burgers rheology (*39*, *40*). However, Freed *et al.* (*11*, *12*) and Masuti *et al.* (*4*) argue that both transient creep and power-law flow may be needed to explain the geodetic data set that they considered such that transient creep constitutes a pathway to the steady state. Following these considerations, we combine both the power-law rheology at the steady state with a power-law rheology for transient creep. We model viscoelastic relaxation with a Burgers spring-slider assembly with the following nonlinear steady-state flow law in the Maxwell element*A*_{M} is the pre-exponential factor, σ is the norm of the deviatoric stress tensor, *n* is the stress exponent, _{M} is the steady-state strain rate, *Q* is the activation energy, *P* is the pressure, *V* is the activation volume, *R* is the universal gas constant, *T* is the temperature, and *C*_{OH} and *r* are the water content and its exponent. The strain rate of the transient phase is represented by a deforming Kelvin element (*4*)*A*_{K} is the pre-exponential factor, *G*_{K} is the strain-hardening coefficient, and ε_{K} is the cumulative transient strain. We assume that the transient phase has the same sensitivity to temperature and stress as in the steady state (*34*). At the steady state before the mainshock, we assume that transient creep does not operate because of the long recurrence time between large earthquakes. Accordingly, we set 2*G*_{K}ε_{K} = σ_{0} as the initial condition, where σ_{0} is a background stress of the order of 1 kPa.

Stress-driven afterslip is modeled with a velocity-strengthening friction law (*17*–*21*)*41*, *42*) used to simulate all phases of the seismic cycle (*43*–*47*), where Δτ is the stress change after the mainshock, *a* − *b* is the steady-state friction parameter, σ_{n} is the effective normal stress on the fault accounting for the effect of pore pressure, and *v* represents the velocity of afterslip. Equation 3 features a nonlinear response to large stress and regularization at vanishing stress compatible with linearization of the constitutive relationship (*48*). The numerical simulations are then performed using the integral method within the two-dimensional in-plane strain approximation using a combination of surface and volume cross-sectional elements (*14*, *49*, *50*).

## RESULTS

Many coastal areas in NE Japan experienced coseismic subsidence and are now undergoing postseismic uplift (Fig. 1), although they have not yet recovered to their pre-earthquake levels (*51*, *52*). Figure 1A shows the observed 5-year postseismic displacement field (*32*), Fig. 1B shows the total displacements along the transect, and Fig. 1C shows the time series of displacements for the GPS stations identified in Fig. 1A. Notably, the forearc is still uplifting (e.g., Oshika station in Fig. 1C), while the subsidence along the volcanic arc and in the backarc (e.g., Naruko and Toyoura stations) has apparently ceased within the observation period (Fig. 1C).

We develop a rheology model to explain the postseismic deformation of the Tohoku earthquake over 5 years (Fig. 1), taking into account two dominant deformation mechanisms of the postseismic deformation: viscoelastic relaxation and afterslip. By comparison to the observed displacement patterns and time series, the constitutive parameters for viscoelastic flow and fault friction are determined as follows (see also Table 1). We created a thermal structure extended from (*53*) that includes several key characteristics of the subducting zone including a cold subducting plate and corner flow (Fig. 1D). The thermal boundary condition beneath the backarc was set on the basis of the potential temperature of 1400°C at 180-km depth in mantle wedge (*53*). A preliminary sensitivity study showed that the surface deformation is not greatly sensitive to the thermal structure below 200-km depth, within reasonable bounds on the adiabatic thermal gradient. We used the laboratory-derived flow law of olivine (*8*) to calculate the steady-state viscosity with the following parameters: *Q* = 430 kJ/mol, *r* = 1.2, *n* = 3, *V* = 13.5 cm^{3}/mol, and the pre-exponential factor *A*_{M} = 10^{0.56} s^{−1} MPa^{−n}. A recent study suggests a different sensitivity to water content of transient creep and steady-state creep, if they represent activation of different slip systems in olivine crystals (*38*). However, surface geodetic observations alone cannot distinguish the roles of water content and other prefactors independently, and lacking further constraints from laboratory experiments, we assume the same water sensitivity and prescribe *A*_{K} = 10 *A*_{M}, compatible with a large reduction of effective viscosity during transient creep (*11*, *39*). As the postseismic response of a nonlinear material depends on the background stress, we assumed the presence of a horizontal compressive prestress compatible with an average uniaxial strain rate across the island arc of 10^{−16} to 10^{−15}/s.

To explain the observations in detail, we assume heterogeneity in water content in the mantle wedge and the oceanic mantle. However, given the epistemic uncertainties related to the rheology of transient creep and the specific focus of the study to assess the importance of the mechanical coupling between afterslip and viscoelastic flow in the years that follow a giant earthquake, we do not place much weight on the inferred lateral variations of water content. However, the average water contents of oceanic mantle and mantle wedges are estimated to be ~260 and ~3100 parts per million (ppm) H/Si, respectively, compatible with dehydration of the oceanic mantle at spreading centers and hydration of the mantle wedge by subducted sediments and metamorphic reactions downdip from the slab. The frictional parameter of stress-driven afterslip was constrained in particular by the surface vertical motion. Data are best explained by a frictional parameter of (*a − b*)σ ~ 0.8 MPa and *V*_{0} = 1.9 × 10^{−9} m/s throughout the downdip extension of the main rupture region.

The combination of a power-law rheology for distributed deformation with velocity-strengthening stress-driven afterslip on the plate interface explains the total observed two-dimensional horizontal (Fig. 2A) and vertical (Fig. 2B) displacement fields well. The horizontal displacement is dominated by viscoelastic relaxation. However, the vertical deformation is more sensitive to rheological heterogeneities. In some regions, the surface displacements caused by afterslip and viscoelastic relaxation may have the opposite sense, in particular the vertical motion in the forearc region (Fig. 2B). In those regions, deeper viscoelastic relaxation results in subsidence, while shallower afterslip drives uplift in the 5-year period. In general, these results are compatible with previous studies (*6*, *29*).

Figure 2 (C to H) illustrates the GPS time series of selected locations across the NE Japan arc (site locations in Fig. 1A). With a similar spatial variation as the observed displacement, viscoelastic relaxation dominates the horizontal component of the time series at almost all locations. On the other hand, opposite contributions of both mechanisms can be seen in the time series of vertical deformation: Afterslip drives steady uplift, and viscoelastic relaxation results in subsidence at the forearc (Oshika station in Fig. 2D). In particular, the contribution of afterslip is large in the forearc region close to the source area and becomes smaller with distance from the source. Additionally, for both components of displacement, the time dependence is different for the two mechanisms: Viscoelastic deformation is initially rapid because of the transient and nonlinear response but decays with time according to temporal changes in effective viscosity driven by the coseismic stress change (Fig. 3, A and B, and figs. S4 and S5).

Figure 4 shows the afterslip distribution at 5 years after the coseismic event. The afterslip mainly occurs at the region directly downdip of the main rupture area, where coseismic stress changes are large (depth of ~30 km; figs. S1 and S2). We find a peak afterslip after 5 years of ~2.5 m and a slip rate (average rate of ~0.5 m/year in 5 years; Fig. 4A), consistent with some previous estimates from the early postseismic observations (*28*, *29*). The second sharp peak in afterslip at the 180-km distance from the trench is likely driven by large gradients in coseismic slip associated with local stress concentration.

## DISCUSSION

Although the afterslip rate gradually decreases with time in the 5 years of the postseismic period (Fig. 4B), it causes persistent uplift in the coseismically subsided coastal area (Fig. 2). The lack of deep afterslip in the postseismic deformation of the Tohoku earthquake had previously been inferred (*33*) from 5 years of GPS horizontal motion. The difference in results between these analyses may be due to the continuation of stress-driven afterslip embedded in a material deforming with a power-law rheology. Strong viscosity reduction and mechanical coupling causes continuous loading of the plate boundary fault driving persistent afterslip directly downdip of the main rupture (Fig. 3A and figs. S4 and S5).

Equipped with a reference model that fits 5 years of geodetic observations, we may now discuss the importance of the mechanical coupling between afterslip and viscoelastic flow. Here, the coupling allows two different interactions: viscoelastic flow triggered by the occurrence of afterslip (afterslip to viscoelastic relaxation = AV coupling) and persistent afterslip influenced by viscoelastic flow (= VA coupling). The details of the coupling depend on the geometry of the viscoelastic materials and the megathrust, but for the NE Japan subduction zone, we find that neglecting the mechanical interaction between viscoelastic flow and afterslip on the plate boundary fault introduces a variance of up to 25% (Fig. 4A) on the amplitude of afterslip below the coast (around 70- to 80-km depth). The difference in afterslip rate by the coupling reaches up to half of the steady-state subduction rate (Fig. 4C) and can alter the surface displacements (dashed lines in Fig. 2 and fig. S3). Only allowing changes in tractions on the plate boundary fault due to viscoelastic flow (VA coupling) decreases the afterslip amount at around 70- to 100-km depth. This is close to the region where the viscosity of the mantle wedge is decreased by the coseismic stress change (e.g., Fig. 3, A and B, and figs. S4 and S5), indicating that viscoelastic flow in this region effectively relaxes the stress on the plate boundary fault. Conversely, the viscoelastic flow loads the fault at around 50- to 80-km depth, driving persistent afterslip (Fig. 4). Notably, only allowing the stress on viscoelastic materials caused by afterslip (AV coupling) results in more vertical deformation (i.e., uplift) in the coastal regions (dotted line in fig. S3B). This is due to the stress on the plate boundary that can only be relaxed by afterslip, not being damped by the global viscoelastic relaxation (VA coupling) occurring in mantle wedges and oceanic mantle, and causing a stress-driven afterslip directly downdip of the main rupture.

The effect of mechanical coupling is most important in the region directly downdip of the main rupture in the forearc mantle wedge (Fig. 5), where the effective viscosity was greatly reduced in the transient phase by the large stress perturbation from the mainshock (Fig. 3, A and B). The importance of coupling grows as the postseismic transient develops, and up to 15% of viscous strain in mantle wedge over 5 years is caused by the afterslip via this coupling (Fig. 3, C and D). It is smaller in the early phase of the postseismic deformation and may, in some cases, be neglected during this period (time series in Fig. 2 and fig. S3). However, the coupling becomes more pronounced with time. Eventually, the coupling drives 20 to 30% of the total strain at the downdip region of the main rupture (Fig. 5). Substantial differences in the postseismic uplift on the order of ~10 cm might occur at the coastal regions over 5 years (fig. S3B). Thus, we highlight the importance of addressing mechanical coupling for long-term studies of postseismic relaxation, in particular when combining kinematic inversion with forward models for viscoelastic relaxation, where estimates of the recovery of plate coupling and postseismic uplift may be significantly modified. Incorporating the constitutive rock behavior constrained by mineral physics and including the mechanical coupling between localized and distributed deformation may be key to better predict the future recovery of the coseismically subsided coastal area through postseismic uplift and quantify the seismic hazards of inland faults with the progress of the postseismic deformation.

## MATERIALS AND METHODS

### Numerical model

We developed a two-dimensional model to explain the postseismic deformation of the Tohoku earthquake over 5 years (Fig. 1). The simulation was performed using the integral method in the condition of in-plane strain (*14*, *49*, *50*). We calculated the stress interactions among the embedded surface and volume elements with analytic solutions for both displacements and distributed strain caused by viscoelastic deformation of a finite cross-sectional area. The evolution of stress, slip rate, strain rate, and the frictional state parameter was then solved numerically throughout the postseismic period using the integral method.

The simplification of the model introduced by the two-dimensional approximation allows us to explore a realistic rheology for the NE Japan subduction zone compatible with mineral physics and to explain the wealth of geodetic observations. Following a previous study (*29*), we construct a two-dimensional transect (black line in Fig. 1A) normal to the trench axis around the Miyagi-Yamagata area in the central part of NE Japan, where the peak coseismic slip was observed (*54*). The coseismic slip distribution calculated along the transect based on a three-dimensional coseismic slip distribution by (*54*) in a previous study (*29*) was used. The model covers 800 km to the west and 200 km to the east from the trench with 400-km depth. We considered slip on the plate boundary fault (patch width of 2 km) and distributed viscoelastic deformation of the surrounding inland lower crust, mantle wedge, and oceanic mantle (rectangular areas of 5 km × 5 km).

The steady-state viscosity (Fig. 1E) was estimated from the laboratory-derived flow laws of olivine (*8*). We included spatial variations in viscosity due to the temperature-dependent rheology obtained from a subduction zone thermal flow model (Fig. 1D) (*53*). The thermal model from (*53*) only covered 600 km west from the trench, so we extended it to 800 km by solving a one-dimensional half-space cooling model. Thermal structure beneath the oceanic lithosphere was created according to (*55*) and combined into the thermal model of the island arc by (*53*). The steady-state viscosity (Fig. 1E) was estimated using the thermal structure with heterogeneous water distribution. The viscosity structure also includes notable effects from the weak volcanic front and strong forearc (cold nose) based on the results of previous studies (*1*, *29*, *31*). We searched the following optimum ranges of water contents to explain the observations: average water contents of 260 and 3100 ppm H/Si for oceanic mantle and mantle wedge, respectively. The former is consistent with a previous estimate of the viscosity structure of the oceanic mantle from the early postseismic deformation of the 2011 Tohoku earthquake (*29*).

We modeled afterslip using the regularized rate-strengthening approximation of the rate-and-state–dependent friction law (*17*–*21*). Afterslip is not allowed in fault regions of negative coseismic stress change. Plate boundary was modeled by shallow coseismic slip region and deep afterslip region. We explored the values of reference slip velocity *V*_{0} and (*a − b*)σ values in Eq. 3 that best explain the geodetic observations in combination with viscoelastic relaxation under the stable subduction of the Pacific plate (8 cm/year). In particular, afterslip location at the downdip region of the main rupture can be constrained by comparison to the vertical motion at the forearc (e.g., Oshika station).

### Geodetic observation

For the postseismic deformation of the 2011 Tohoku earthquake, we incorporated observations from continuous GPS stations on land deployed by the Geospatial Information Authority of Japan (GEONET, gray arrows in Fig. 1A) and Tohoku University as well as seafloor observations deployed by Japan Coast Guard and Tohoku University (white arrows) (*14*). In the 6 years of the postseismic period, dense and wide GPS-A observation network was constructed along the Japan trench in September 2012 and six campaign surveys were conducted at most site until May 2016 (*32*). Accumulated postseismic displacements were calculated for the period. The displacements for both the on-land GPS and seafloor GPS-A data were estimated relative to the North American plate. We did not correct the effect of interplate coupling during the period because we could not obtain the pre-earthquake deformation rate for the seafloor observations that was installed after the earthquake.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/9/eaaw1164/DC1

Fig. S1. Coseismic stress change calculated by the coseismic slip model of (*20*).

Fig. S2. Stress change caused by cumulative afterslip at 6 years after the 2011 Tohoku earthquake in a fully mechanical coupling model.

Fig. S3. Effect of mechanical coupling on the GPS time series.

Fig. S4. Temporal change in effective viscosities of the Kelvin element within the Burgers rheology.

Fig. S5. Temporal change in effective viscosities of Maxwell medium in Burgers rheology.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We acknowledge S. Horiuchi for help on the creation of thermal structure. We also acknowledge R. Bürgmann and an anonymous reviewer for many helpful suggestions that improved the manuscript.

**Funding:**This study was funded by a JSPS KAKENHI grant numbers JP15K21755 and JP26109007 and by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of Japan, under its Earthquake and Volcano Hazards Observation and Research Program. This research was supported, in part, by the National Research Foundation (NRF) of Singapore under the NRF Fellowship scheme (National Research Fellow Award NRF-NRFF2013-04) and by the Earth Observatory of Singapore under the Research Centres of Excellence initiative.

**Author contributions:**J.M., J.D.P.M., and S.B. designed the study and wrote the manuscript with contributions from all other coauthors. J.M. conducted numerical simulation. T.I. and Y.O. prepared GPS data. H.I. conducted the numerical model on the thermal structure.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).