Rupture Reactivation during the 2011 M_w 9.0 Tohoku Earthquake: Dynamic Rupture and Ground‐Motion Simulations

The 11 March 2011 M_w 9 Tohoku earthquake occurred in the subduction zone between the Pacific and North American plates in northeastern Honshu, Japan. This giant event was recorded by a vast Global Positioning System (GPS) and seismic network, providing a unique opportunity to study the complex rupture process of a megathrust event (e.g., Ide et al., 2011; Simons et al., 2011; Tajima et al., 2013). Visual inspection of the strong ground motion recorded by the K‐NET and KiK‐net networks and the high‐rate 1 Hz GPS displacements recorded by the GNSS Earth Observation Network System (GEONET) network along the Japanese coast clearly show two groups of waves (S1 and S2 in Fig. 1) arriving about 40 s apart in both the Miyagi and Iwate regions. A third group of waves (S3 in Fig. 1) is observed further south in the Ibaraki and Chiba regions. These observations, together with kinematic source inversion models and backprojection source imaging results (e.g., Ide et al., 2011; Ishii, 2011; Koketsu et al., 2011; Lee et al., 2011; Meng et al., 2011; Simons et al., 2011; Suzuki et al., 2011; Roten et al., 2012; Yagi et al., 2012), indicate that the Tohoku earthquake featured complex rupture patterns involving multiple rupture fronts. These studies collectively reveal three remarkable features: (1) large shallow slip exceeding 40 m near the trench; (2) depth‐dependent frequency content of seismic radiation with stronger high‐frequency radiation in the deeper region and larger low‐frequency slip at shallow depth; and (3) repeated slip in the region near and up‐dip of the hypocenter.

Within the framework of dynamic rupture models, a number of studies have proposed possible mechanisms to explain these unusual rupture features. Although shallow megathrusts generally display velocity‐strengthening friction due to the presence of sediments in the fault zone (e.g., Marone and Scholz, 1988; Oleskevich et al., 1999), which promotes stable sliding, it is now recognized that sufficiently energetic dynamic ruptures can break through frictionally stable regions and low‐stress regions (Kaneko et al., 2010). The possibility of large slip near the trench can be further enhanced by dynamic stresses in the subduction wedge associated with seafloor effects (Kozdon and Dunham, 2013; Huang et al., 2014) and by dynamic weakening mechanisms operating at large slip or slip velocity in materials that are velocity strengthening at low speed (Mitsui, Kato, et al., 2012; Noda and Lapusta, 2013; Cubas et al., 2015). Enhanced high‐frequency radiation at depth has been attributed to the presence of deep asperities on the megathrust (Mitsui, Iio, and Fukahata, 2012; Huang et al., 2013, 2014) and its deprivation at shallow depth to thermal pressurization (e.g., Noda and Lapusta, 2013) and inelastic failure in the wedge (Ma and Hirakawa, 2013).

Repeated slip during the Tohoku earthquake has been mostly characterized as a rupture front backpropagating down‐dip after the initial up‐dip front reached the trench (Ide et al., 2011). This phenomenon is usually found in dynamic rupture simulations of dipping faults and is caused by free surface effects (e.g., Dalguer et al., 2001a,b; Huang et al., 2014). Some observational studies conclude instead that repeated slip occurred by repeated nucleation of rupture near the initial hypocenter. The presence of two sequential rupture fronts starting from the same hypocentral area during the Tohoku earthquake, a phenomenon denoted here as “rupture reactivation,” was first reported by Lee et al. (2011) in a kinematic source inversion study based on teleseismic, strong ground motion, and GPS data. The feature is illustrated in Figure 2a (see also Gabriel et al., 2012), which shows the spatiotemporal evolution of slip velocity along a cross section of the fault oriented along‐dip and passing through the hypocenter. Two rupture fronts emerge from the hypocenter area (black lines): the first one propagates mainly toward the trench, and the second one starts 40–50 s after the earthquake origin time and propagates bilaterally. Lee et al. (2011) parametrized the slip velocity functions in their source inversion by multiple time windows. Rupture reactivation is also apparent in another source inversion model derived from teleseismic, strong motion, and geodetic data (Koketsu et al., 2011), in backprojection source imaging results based on data from the Metropolitan Seismic Observation Network (MeSO‐net) in the Tokyo metropolitan area (Honda et al., 2011), and in a source model inferred from regional ground‐motion data by an iterative deconvolution and stacking technique (Zhang et al., 2014). A common methodological aspect of these studies is their allowance for multiple rupture fronts by either an unrestricted parameterization or a nonparametric description of the source.

Two mechanical models have been proposed to explain rupture reactivation: rupture in a heterogeneous initial stress field (Goto et al., 2012) and renucleation by a growing stress concentration left behind an initial self‐similar slip pulse (Gabriel et al., 2012). Here, we develop a new rupture reactivation model based on a double‐slip‐weakening friction law similar to that introduced by Kanamori and Heaton (2000). These authors proposed that melting or fluid pressurization induced by frictional heating can reduce fault strength when slip exceeds a certain critical slip distance (D_r). They argued that, when superimposed to the conventional (isothermal) slip‐weakening process with shorter critical slip distance (D_c), these thermally activated weakening mechanisms lead to a slip‐dependent friction model with two sequential strength drops (Fig. 3).

The main goal of this work is to show that the double‐slip‐weakening friction in the form proposed by Kanamori and Heaton (2000) (Fig. 3) offers a plausible model for the rupture process of the 2011 Tohoku earthquake, including large shallow slip, rupture reactivation, and large rupture extent. First, we present evidence supporting this friction model, obtained from seismological observations, laboratory experiments, and theoretical considerations. Then, we extend the dynamic rupture model developed by Galvez et al. (2014; hereafter referred to as G14) by modifying the slip‐weakening friction model to account for rupture reactivation. Finally, we show that this model produces (1) a rupture pattern consistent with a kinematic source inversion model that features rupture reactivation and (2) ground‐motion patterns along the Japanese coast consistent with the observations.

The stress changes inferred from kinematic models provide constraints on frictional strength evolution during rupture (e.g., Bouchon, 1997; Ide and Takeo, 1997; Dalguer et al., 2002). We compute the dynamic fault stress changes implied by the Lee et al. (2011) kinematic source model by applying its spatiotemporal distribution of slip velocity as a boundary condition along the fault in a spectral element seismic‐wave propagation simulation done with the SPECFEM3D code. Figure 2b, 2c, and 2d, respectively, shows the along‐dip slip velocity, slip, and computed shear stress change around the hypocenter as a function of time. The stress change features two sequential drops, correlated in time with two peaks of slip velocity. Figure 2e shows the stress change as a function of slip. The similarity with the proposed slip‐weakening friction model shown in Figure 3 is remarkable. In particular, it shows that the critical slip for the onset of the second weakening is D_r≈20  m in the hypocentral region.

Additional evidence of secondary frictional weakening is found in laboratory rock experiments. Dramatic weakening is often observed in rock friction experiments at sliding velocities representative of seismic slip (∼1  m/s) and has been attributed to thermomechanical and physicochemical processes, including flash heating, melting, and gelification, decarbonation, and dehydration reactions (e.g., Di Toro et al., 2012). Han et al. (2007) performed experiments on calcite at seismic slip rates and found that decarbonation induces pronounced fault weakening, which they attributed to flash heating on ultrafine particles produced by thermal decomposition. The high‐velocity frictional experiments of O’Hara et al. (2006) on bituminous coal gouge, proposed as an analog for hydrous fault zones, showed two major strength drops (see fig. 2 in O’Hara et al., 2006). The second weakening occurred after 20–25 m of slip and was attributed to gouge compaction and ensuing pressurization of gas released by thermal decomposition.

Theoretical studies provide further support for a double‐weakening friction model enabled by thermochemical decomposition. We first note that thermal pressurization alone produces gradual weakening (Rice, 2006), in contrast to the conceptual model by Kanamori and Heaton (2000), and that melting temperatures are not reached during earthquake slip in rupture models that include thermal pressurization (Garagash, 2012). Sulem and Famin (2009) developed a numerical model of fault‐strength evolution, including shear heating, thermal pressurization, and the mechanical effects of calcite volume loss, heat consumption, and CO₂ production. In their simulations at prescribed slip velocity, a second strength drop occurs due to additional pressurization induced by rapid decarbonation when a critical temperature of about 700°C is reached. Strengthening is subsequently caused by the increase of porosity and permeability due to the solid decomposition. Brantut et al. (2010) developed analytical models for strength evolution of hydrous fault zones and showed that dehydration reactions can cause a significant second weakening phase if the reaction rate is comparable or faster than the thermal pressurization process, which is favored by a thick slip zone. Further theoretical analysis of the problem was presented by Platt, Brantut, and Rice (2015) and Platt, Viesca, and Garagash (2015).

On the basis of estimates of peak temperature achieved on a fault undergoing thermal pressurization, Garagash (2012) argued that thermal decomposition could be triggered at depths larger than ∼15  km, or ∼6  km in a damaged fault zone. Theoretical work indicates that the occurrence and significance of weakening by thermal decomposition depends not only on fault‐zone material properties, but also on background stress. Platt, Viesca, and Garagash (2015) showed that weakening by thermal decomposition becomes significant if it overcomes the strengthening effect of hydraulic diffusion; this occurs if the temperature exceeds a critical value, which they found to be higher if the background shear stress is lower. Moreover, in dynamic rupture simulations with thermal pressurization, the temperature rise is higher if background stress is larger (Noda et al., 2009). This implies that, comparing fault areas of identical material properties and initial temperature, thermal decomposition weakening is favored in areas of higher stress, which can operate as asperities (large slip areas) during an earthquake, and is less expected in their surrounding, lower‐stress areas.

Field evidences that mineral decomposition reactions may occur during earthquakes have been found in cores from the Chelungpu fault that hosted the 1999 Chi‐Chi, Taiwan earthquake (Hirono et al., 2008, 2015; Hamada et al., 2009) and from the Nojima fault that hosted the 1995 Kobe, Japan, earthquake (Famin et al., 2008). In particular, subduction fault gouges rich in clay minerals or phyllosilicates such as chlorite, talc, and serpentinite can host dehydration reactions at temperatures reached during earthquake slip.

The seismological, experimental, and theoretical arguments presented in the previous discussion support the possibility of double weakening, in particular induced by thermal decomposition. However, the purpose of the present work is not to simulate or to assess the viability of a specific physical mechanism in detail. This would require the implementation of coupled thermohydrochemical processes in a 3D rupture‐dynamics simulation code. It would also require efforts to constrain the relevant fault‐zone parameters and to evaluate the effect of their (often large) uncertainties. Instead, our goal is to demonstrate conceptually, through simulations based on a phenomenological slip‐weakening friction law, that the assumption of double weakening is consistent with basic rupture and ground‐motion patterns.

We extend the model developed by G14 to account for rupture reactivation. G14 developed an asperity model to simulate the dynamic rupture process of the 2011 Tohoku earthquake, in particular to reproduce qualitatively the observed depth‐dependent frequency content of seismic radiation and large slip close to the trench. The model included a nonplanar megathrust fault geometry adapted from Simons et al. (2011). The unstructured mesh‐generation software CUBIT was utilized to create a spectral element mesh accounting for these geometrical complexities. Our simulations employ the unstructured spectral element code SPECFEM3D (Peter et al., 2011), in which G14 implemented dynamic rupture capabilities to enable accurate simulations of dynamic rupture and seismic‐wave propagation. G14 demonstrated the stability and accuracy of SPECFEM3D for dynamic rupture simulations in the subduction wedge geometry of Tohoku. In particular, results of spatial resolution tests presented in appendix A of G14 show that the rupture process is accurately resolved with an average grid size of 500 m on the fault, a value we also adopt in the present model (i.e., fourth‐order spectral elements of size 2000 m). In that test, decreasing the average grid size to 125 m (fourth‐order spectral elements of size 500 m) provided similar results.

The G14 model is composed of two large asperities (1 and 2) and a collection of small circular asperities in deeper regions, as shown in Figure 4. The location of the two main asperities corresponds approximately to the large shallow slip inferred from low‐frequency (LF) data through kinematic source inversions (e.g., Lee et al., 2011; Suzuki et al., 2011; Yagi and Fukahata, 2011; Yoshida, 2011; Yue and Lay, 2011; Wei et al., 2012). The small asperities are placed in the general areas of high‐frequency (HF) radiation identified by backprojection techniques (e.g., Ishii, 2011; Meng et al., 2011; Simons et al., 2011; Roten et al., 2012; Yagi et al., 2012). Figure 4a compares the location of the asperities with the final slip model of Lee et al. (2011), denoted with contour lines, and with the HF radiation points (green circles) obtained by Meng et al. (2011). The largest asperity (asperity 1) was defined as an ellipse that delimits the region of large slip inferred by Lee et al. (2011). The depth dependency of LF and HF radiation implied by this rupture model was presented by G14, and other studies based on asperity models were developed by Aochi and Ide (2011, 2014), Huang et al. (2013, 2014), and Ide and Aochi (2013). The shallow zone of megathrusts is characterized by velocity strengthening friction due to the presence of sediments in the fault zone (e.g., Marone and Scholz, 1988; Oleskevich et al., 1999). To approximately mimic this mechanism, we impose a negative stress drop at shallow depth (e.g., Dalguer et al., 2008; Pitarka et al., 2009).

Here, we extend the G14 model to account for rupture reactivation by incorporating the double‐slip‐weakening friction shown in Figure 3. The rupture in the largest asperity is governed by our proposed double‐slip‐weakening friction. This asperity breaks with two sequential stress drops, each of 4.5 MPa, resulting in a total stress drop of 9 MPa. In the rest of the fault, where initial stresses are lower, rupture is governed by the standard linear slip‐weakening friction model (Andrews, 1976). The distributions of (1) nominal stress drop, defined as (initial along‐dip shear stress) – (dynamic friction coefficient) × (initial normal stress), (2) nominal strength excess, defined as (static friction coefficient) × (initial normal stress) – (initial along‐dip shear stress), and (3) critical slip distance (shown in Figs. 4b, 4c, and 5, respectively) were obtained by a systematic tuning of the dynamic parameters (initial stresses, friction coefficients, and D_c) to qualitatively fit the kinematic source model of Lee et al. (2011) and the observed ground‐motion pattern shown in Figure 1. The blue and red regions of Figure 4b and 4c represent the artificial nucleation patch, where the static coefficient is slightly reduced so that the initial yield strength, defined as (static friction coefficient) × (initial normal stress), is slightly below the initial along‐dip shear stress. We adopt the 1D layered velocity structure of Fukuyama et al. (1998) (Table 1).

Our dynamic rupture model generates a rupture reactivation process consistent with the one observed in the kinematic source model of Lee et al. (2011). Snapshots of the simulated slip rate are shown in Figure 6. An initial rupture front propagates up‐dip, with subshear speed, and reaches the trench at t≈35  s (time is relative to the rupture onset), breaking a small portion of it. Simultaneously, the rupture also propagates down‐dip, although slowly, and starts to break the closest deep small asperities. Hereafter, we refer to this first episode as phase 1. At t≈40  s, a new episode initiates with an energetic rupture reactivation around the hypocenter (phase 2). To better illustrate this, in Figure 7 we show the spatiotemporal evolution of slip velocity along the longitudinal transect passing through the hypocenter (indicated as section PQ in the first snapshot of Fig. 6). Similar to the kinematic source model of Lee et al. (2011) (Fig. 2a), two sets of rupture fronts emerge from the hypocenter area, corresponding to phases 1 and 2 (black lines). During phase 2, the new rupture front propagates both up‐dip and down‐dip. The upward front propagates with supershear speed, reaches the trench at t≈50  s, and then propagates bilaterally along strike. The downward rupture propagates with subshear speed, penetrates the deepest region, and generates high peak‐slip velocities when it breaks the small asperities. These sharp peak‐slip velocities are expected to radiate at the high frequencies imaged by backprojection studies but remain unresolved in lower‐frequency kinematic source models (Fig. 2a). At t≈80  s, the deep rupture front starts propagating southward, then travels slowly a few kilometers until it reaches the second big asperity (asperity 2) at t≈100  s. Then a third episode initiates, leading to the complete rupture of asperity 2 (phase 3). The whole rupture stops in the southern part at t≈160  s.

The final slip distributions resulting from our dynamic model and from the kinematic model of Lee et al. (2011) are compared in Figure 8. The locations of the maximum slip area are consistent in both models.

Inspection of the simulated ground‐velocity wavefield and comparison with the distinct rupture phases (Figs. 6 and 7) reveal the correspondence between the rupture process and the prominent wave groups arriving to the Japanese coast. The seismic waves are well resolved by our simulation up to 0.2 Hz. Ⓔ The movie, available in the electronic supplement to this article, shows in detail the evolution of the east–west ground velocity, filtered in the 0.01–0.05 Hz frequency range to highlight the main groups of waves. Four selected snapshots of ground velocity are shown in Figure 9a, along with simulated seismograms (Fig. 9b) at the locations of the K‐NET and KiK‐net seismic stations distributed along the Japanese coast (indicated by white squares in Fig. 9a). Three prominent wave groups are identified. The first group arrives at the Japanese coast between 50 and 65 s and corresponds to the waves radiated from phase 1 of the rupture. The second group of waves reaches the coast between 90 and 100 s and is attributed to radiation from the rupture reactivation during phase 2 of the rupture. The third group of waves hits the coast of the Ibaraki region at about 120 s and is associated with the rupture of the southern asperity during phase 3.

These three groups of waves resemble the S1, S2, and S3 wavefronts identified in the observed ground motions (Fig. 1). The left and middle panels of Figure 10 show, respectively, simulated and observed ground velocities at the selected K‐NET and KiK‐net stations. Their similarity is remarkable, suggesting that our dynamic rupture model with rupture reactivation explains the ground‐motion pattern recorded in the Tohoku earthquake. As a reference, the right panel of Figure 10 shows the ground velocities generated by the dynamic model of G14, which does not include rupture reactivation. This model does not reproduce clearly separated wavefronts S1 and S2, but generates a single S2 front with larger amplitudes and longer duration than those observed, because the entire stress drop (9 MPa) in asperity 1 is released at once.

In Figure 11, we compare the three components of simulated and observed ground velocities at three representative KiK‐net stations located onshore of the north (IWTH20), middle (MYGH08), and south (IBR002) portions of the rupture of the model with reactivation and model G14 without reactivation. The seismograms are low‐pass filtered with a frequency cutoff of 0.2 Hz. In general, synthetics of our model with reactivation follow the waveform patterns of observations better than those of the G14 model. The best agreement is achieved at the middle station (MYGH08) and the worse at the southern station (IBR002). Waveform shapes and amplitudes in the horizontal components agree better than in the vertical component.

We also compare the ground displacements simulated by our model with rupture reactivation (Fig. 12a) and by model G14 without reactivation (Fig. 12b) to the observed continuous GPS data recorded along the Japanese coast (east–west component). Only our model with reactivation shows good consistency with the observed displacements, in shape and amplitude. In particular, it reproduces the three main groups of waves that are also observed in the GPS data (Fig. 1). As shown in Figure 12a, our model with rupture reactivation better explains the displacements at stations located in the middle latitude of the rupture (Miyagi and Fukushima areas) where the two groups of waves S1 and S2 are radiated (according to our model) from the repeated rupture of asperity 1.

These comparisons suggest that the two episodes of stress drop leading to rupture reactivation are needed to reproduce the observed ground‐motion pattern.

In the framework of our dynamic model, we investigate the role of rupture reactivation and deep small asperities in enabling the Tohoku earthquake to reach its final magnitude, M_w 9. For this purpose, results of three dynamic rupture models are compared (Fig. 13a). Model A is the dynamic rupture model with rupture reactivation in asperity 1 (two stress drops of 4.5 MPa each) and multiple small asperities at greater depth. Model B is identical to model A but without rupture reactivation (asperity 1 has a single stress drop of 4.5 MPa). Model C is identical to model A but with the small asperities removed.

Model B produces an M_w 8.6 earthquake, which breaks only the small asperities located just below asperity 1 and leaves the second large asperity unbroken. Model C produces an M_w 8.75 event, again without breaking asperity 2. Figure 13b shows the final slip distribution of each model. Maximum slip values are around 50, 20, and 40 m, respectively, for models A, B, and C. The slip of models B and C spans an area confined to the vicinity of asperity 1. Figure 13c shows the seismic moment rates of each model. Our preferred model, model A, generates three pulses in the seismic moment rate, corresponding to the three rupture phases described previously (see Figs. 6 and 7) and to the three wave groups identified in Figures 1 and 9. Model B generates only the first pulse (phase 1), with amplitude slightly lower than the one produced by model A. Model C generates only the first and second pulses (phases 1 and 2) and the latter, caused by rupture reactivation, is smaller (about one‐third) than in model A.

These model comparisons suggest that both the lack of rupture reactivation (model B) and the absence of small asperities (model C) cause the down‐dip rupture to die out earlier, limiting the rupture area and final magnitude of the earthquake. The series of deep small asperities serves as a bridge to connect the rupture of asperities 1 and 2, which enhances the accumulation of slip close to the trench and enables rupture growth up to an M_w 9.0 event. In the context of our study, both rupture reactivation and deep small asperities are necessary to generate a megathrust event of a size consistent with observations.

We developed a 3D dynamic model of the rupture process of the 2011 Tohoku earthquake that includes the rupture reactivation phenomenon (repeated nucleation of rupture in the hypocentral area) previously inferred from seismological observations. Our model features a slip‐weakening friction with two sequential strength drops, supported by the evolution of dynamic stresses derived from a finite source inversion model and by laboratory experiments involving thermochemical weakening mechanisms associated with mineral decomposition reactions.

Three prominent groups of seismic waves observed in ground motions recorded along the Japanese coast are generated in our model by three notable rupture phases. The first two originate from the main asperity, which contains the hypocenter and experiences rupture reactivation. The third rupture phase involves a second large asperity located in the southern portion of the rupture. The ground velocity and displacement waveforms recorded by the K‐NET and KiK‐net strong‐motion networks and the GEONET GPS network are generally reproduced by our numerical simulation. These results show that a friction model with rupture reactivation can explain the complex dynamic rupture process of the 2011 Tohoku earthquake and the observed ground‐motion patterns near the source.

Our dynamic model has two main ingredients, rupture reactivation of the main asperity and the presence of deep small asperities. The latter were previously found to explain the enhanced high‐frequency radiation in the deeper portions of the Tohoku earthquake rupture zone, as revealed by teleseismic backprojection source imaging. Here, we found that both ingredients may have played a crucial role in enabling the large size of this megathrust rupture. Our modeling results suggest that the small asperities serve as a bridge to connect the rupture of the hypocentral asperity to that of the southern asperity; however, without rupture reactivation in the first asperity, this bridge is not activated.

The additional strength drop caused by decomposition reactions triggered by shear heating can significantly enhance the final rupture size and rupture complexity of large earthquakes.

The strong ground motion data used in this work are from the K‐NET and KiK‐net database (http://www.kyoshin.bosai.go.jp, last accessed November 2013) and the Global Positioning System (GPS) records are from the GNSS Earth Observation Network System (GEONET) of Japan, operated by the Geospatial Information Authority of Japan (GSI; http://www.gsi.go.jp/ENGLISH/index.html, last accessed November 2013).

This study was supported by the Quantitative Estimation of Earth’s Seismic Sources and Structure (QUEST) project funded by the Seventh Framework Programme of the European Commission, the Advanced Simulation of Coupled Earthquake and Tsunami Events (ASCETE) Project funded by the Volkswagen Foundation within the program “New Conceptual Approaches to Modeling and Simulation of Complex Systems,” and National Science Foundation (NSF) CAREER Award EAR‐1151926. We thank the Swiss project “High‐Rate GNSS for Seismology,” funded by the Swiss National Foundation, for providing us with the Global Positioning System (GPS) records of the GNSS Earth Observation Network System (GEONET). Simulations were done at the Swiss National Supercomputing Center (CSCS), under the production projects “Development of Dynamic Rupture Models to Study the Physics of Earthquakes and Near‐Source Ground Motion” and “Development of a Database of Physics‐Based Synthetic Earthquakes for Ground Motion Prediction.”