**2. Model and method**

#### **2.1 Hierarchical asperity model for M9 earthquakes**

#### **2.1.1 Concept**

Aseismic sliding occurs during the M9 earthquake cycle if the frictional property in the M9 source area is apparently stable enough, in other words, if the nucleation size is large enough. As a result, the aseismic sliding can load the smaller unstable locked patches within the M9 source area. If the nucleation size and fracture energy of the patches are much smaller than they are in the rest of the source area, unstable slip occurs on the patches but does not propagate to their outside. This is a possible mechanism that can account for M<9 events within the source area of an M9 event. This mechanism is similar to the foreshock model (Matsu'ura et al., 1992) although the space-time scale is much larger here. The events on the unstable patches can occur repeatedly until the accumulated strain energy is sufficient to cause an M9 event in the entire seismogenic zone, as an unstable slip on a patch triggers the rupture of the area including large fracture energy. Note that similar multi-scale heterogeneity in fracture energy is assumed for the 2011 Tohoku earthquake (Aochi & Ide, 2011) although they focus on the dynamic rupture process.

In the following discussion, the unstable patches and the M9 source area, including the unstable patches, are called regular asperities and a hyper asperity, respectively. The source area outside the regular asperities is called a conditional asperity, because in this area both aseismic sliding and seismic slip occur, depending on the stress and fault strength conditions as shown later.

#### **2.1.2 Friction law**

2 Will-be-set-by-IN-TECH

earthquake (Mw=8.8∼9.0), M∼7 earthquakes occurred in 1904 also in 1993 within the source area of the 1952 earthquake (Johnson & Satake, 1999). After 29 years of the 1957 Alutian earthquake (Mw=8.6), M=7.7 earthquake occurred within the source area (Johnson et al., 1994). Furethermore, for the 1906 Colombia-Ecuador earthquake (Mw=8.8), three Mw=7.7∼7.9 earthquakes occurred in 1942, 1958 and 1979 (Kanamori & McNally, 1982). Among the above M9 earthquakes, it is known that some of them have occurred repeatedly as in the southern Chile. For example, in the Kamchatka an M9 earthquake also occurred in 1737 (Johnson & Satake, 1999). For the 2011 Tohoku earthquake, identical tsunami deposite distribution has been found in Sendai area, northeast Japan (Minoura et al., 2001). The estimated recurrence time interval is several hundreds to a thousand years. Hence, we assume that the M9 earthquake occurrence is the fundamental rupture mode in each subduction zone. The key question then becomes: How can M=7 ∼8 earthquakes occur within the source area of an M9 event during the long-term seismic cycle of M9 earthquakes? Based on the concept of a hierarchical asperity model that was applied to the M3 sequence within an M5 asperity off Kamaishi along the Japan trench (Hori and Miyazaki, 2010), we speculate that such events

In this chapter, first we will introduce the hierarchical asperity model for M9 earthquakes and explain how to represent them as numerical models. Then, we will describe the simulation results, part of which is identical to our short paper (Hori & Miyazaki, 2011). In discussion, we will introduce the fault strength that is introduced by Nakatani (Nakatani, 2001), and discuss the generation mechanism of M9 recurrence with M=7∼8 earthquakes during the M9 interseismic period. Finally, we will compare our results with observation data, especially for

Aseismic sliding occurs during the M9 earthquake cycle if the frictional property in the M9 source area is apparently stable enough, in other words, if the nucleation size is large enough. As a result, the aseismic sliding can load the smaller unstable locked patches within the M9 source area. If the nucleation size and fracture energy of the patches are much smaller than they are in the rest of the source area, unstable slip occurs on the patches but does not propagate to their outside. This is a possible mechanism that can account for M<9 events within the source area of an M9 event. This mechanism is similar to the foreshock model (Matsu'ura et al., 1992) although the space-time scale is much larger here. The events on the unstable patches can occur repeatedly until the accumulated strain energy is sufficient to cause an M9 event in the entire seismogenic zone, as an unstable slip on a patch triggers the rupture of the area including large fracture energy. Note that similar multi-scale heterogeneity in fracture energy is assumed for the 2011 Tohoku earthquake (Aochi & Ide, 2011) although

In the following discussion, the unstable patches and the M9 source area, including the unstable patches, are called regular asperities and a hyper asperity, respectively. The source area outside the regular asperities is called a conditional asperity, because in this area both aseismic sliding and seismic slip occur, depending on the stress and fault strength conditions

could be explained.

the 2011 Tohoku earthquake.

**2.1 Hierarchical asperity model for M9 earthquakes**

they focus on the dynamic rupture process.

**2. Model and method**

**2.1.1 Concept**

as shown later.

In the present study we employ a laboratory-derived rate- and state-dependent friction law (Rice, 1993) that has been succeeded in modeling slip histories over seismic cycles:

$$
\mu = \mu\_\* + a \ln \left(\frac{V}{V\_\*}\right) + b \ln \left(\frac{V\_\* \theta}{d\_\varepsilon}\right),
\tag{1}
$$

where *μ* is a coefficient of friction, *V* and *θ* are slip velocity and state variable on a fault surface, *μ*<sup>∗</sup> and *V*<sup>∗</sup> are references of frictional coefficients and slip velocity, respectively, *a* is an increment of frictional coefficients for a velocity step (direct effect), *b* is a transient decrement of frictional coefficients for the fault to slip over the evolution length *dc*.

In the framework of the rate- and state-dependent friction law, nucleation size *L*∞ and fracture energy *Gc* are given by

$$L\_{\infty} = \frac{c}{\pi} \left(\frac{b}{b-a}\right)^2 \frac{G}{b\sigma} d\_{c\prime} \tag{2}$$

$$G\_c = \frac{1}{2} \left( \log \frac{\theta\_{\rm in}}{\theta\_{\rm out}} \right)^2 b \sigma d\_{c\prime} \tag{3}$$

respectively (Rubin & Ampuero, 2005), where *G*, *σ*, *θin* and *θout* are shear modulus, effective normal stress, state variables just outside and inside the crack tip, respectively. *c* is a geometrical factor that is about 2 for a circular crack. Clearly the nucleation size and fracture energy increase with a state evolution length *dc* in a similar way as by a slip weakening distance *Dc*. Bizzarri & Cocco, for example, suggested from fully dynamic rupture simulations that *dc* is propotional to *Dc* (Bizzarri & Cocco, 2003). Furthermore, it is suggested that *dc* and *Dc* depends on observation scales, explaining its discrepancy between laboratory and seismological observations (Perfettini et al., 2003; Shibazaki & Matsu'ura, 1998). A multi-scale patch distribution of scale-dependent *Dc* is considered to investigate multi-scale heterogeneous dynamic rupture propagation within a continuum media (Ide & Aochi, 2005). They demonstrated that such heterogenities are able to produce realistic features of earthquake size distribution and dynamic rupture patterns on a fault within the continuum class of models.

#### **2.1.3 Representation of a hyper asperity**

Based on those previous studies, we assume that a state evolution length *dc* is a dominant factor to determine the size of an earthquake. Then we introduce a heterogeneous distribution of *dc* in a single asperity (hyper asperity). Areas of smaller *dc* in the hyper asperity represent regular asperities and would have small fracture energy *Gc*. Hence they would produce smaller size of earthquakes. Rest of the same asperity of larger *dc* represents conditional asperity and would have larger *Gc*. This part is persistent to be locked for a longer period and would induce a failure of the entire hyper asperity to generate large size of earthquakes. It should be noted that the size of regular asperities is another factor to determine the style of fault slips. Earthquakes occur if the size is larger than the nucleation length, but otherwise slow slip events occur [e.g., (Kato, 2003)].

On the other hand, we assume that (*a* − *b*)*σeff* (*a* − *b* < 0) would be uniform over the entire hyper asperity because apparent stress drop is nearly constant over a wide range of earthquake size (Ide & Beroza, 2001). In the following sections we set up a numerical model for a hierarchical asperity model and perform numerical simulations. Then we demonstrate that our model is a candidate to account for the earthquake sequence of M=7∼8 within M∼9 source area.

#### **2.2 Governing equations**

In our model, we devide a flat fault, which represents a subducting plate, into many triangular fault cells and calculate variables for each cell. The governing equations are expressed as following. The first equation is a quasi-dynamic equilibrium of shear stress

$$0 = \sum\_{j} K\_{\bar{j}\bar{j}} (u\_{\bar{j}} - V\_{pl}t) - \mu\_{\bar{i}} \sigma\_{\bar{i}} - \frac{G}{2\mathcal{B}} \frac{du\_{\bar{i}}}{dt} \,' \tag{4}$$

where *uj* is a fault slip of *j*-th fault cell, *Kij* is a shear stress responce of the *i*-th cell caused by a unit dislocation of the *j*-th cell, *σ<sup>i</sup>* is an effective normal stress of the *i*-th cell, *μ<sup>i</sup>* is a coefficient of friction of the *i*-th cell, *Vpl* is a plate convergence velocity that is assumed constant over the entire model region, *β* is a shear wave velocity and *G* is a rigidity. To calculate the shear stress response at the center of each cell, which is triangular, stress changes calculated with angular dislocations are summed (Comninou & Dundurs, 1975). The last term on the right hand side is a radiation damping term to approximate energy radiation as elastic wave for high velocity slip (Rice, 1993).

The second equation is a laboratory-derived rate- and state-dependent friction law as (1),

$$
\mu\_i = \mu\_\* + a\_i \ln \left(\frac{V\_i}{V\_\*}\right) + b\_i \ln \left(\frac{V\_\* \theta\_i}{d\_{ci}}\right) \,\tag{5}
$$

where *Vi* = *dui dt* is a slip velocity of the *i*-th cell, *θ<sup>i</sup>* is a state variable for the *i*-th cell, *μ*<sup>∗</sup> is a reference for frictonal coefficients, *V*<sup>∗</sup> is a reference velocity and is set to *Vpl* in this study, and *ai*, *bi*, *dci* are friction parameters for the *i*-th cell. The state variable *θ<sup>i</sup>* follows a state evolution law. Among several versions of evolution laws, we employ a composite law

$$\frac{d\theta\_i}{dt} = \exp\left(-\frac{V\_i}{V\_c}\right) - \left(\frac{\theta\_i V\_i}{d\_{ci}}\right) \ln\left(\frac{\theta\_i V\_i}{d\_{ci}}\right) \tag{6}$$

where *Vc* is a cut-off velocity and is 1.0 <sup>×</sup> <sup>10</sup>−<sup>8</sup> m/s (Kato & Tullis, 2001). This evolution law approximates slip law for *V* >> *Vc* and aging law for *V* << *Vc*.

From the above equations, we derive set of ordinary differential equations for slip velocity *Vi* and state variable *θi*. The governing equations are solved with adaptive time step 4th-order Runge-Kutta algorithm (Press et al., 1996). For initial condition, slip velocity *Vi* and *θ<sup>i</sup>* are assumed to be uniform, and set to be 0.9*Vpl* and *dci*/*Vi*, respectively.

#### **2.3 Model setup in subduction zone**

We show spatial distributions of frictional parameters *A* − *B* (*A* = *aσ*, *B* = *bσ*) and *dc* in Figure 1. From rock friction experiments, *A* − *B* primarily depends on temperature and hence depends on depths, and *dc* increases where temperature is high under wet condition (Blanpied et al., 2001). For simplicity, *A* is assumed constant on the entire fault. Parameters *B* and *dc*

are, hence, assumed to depend on depths, representing their temperature dependence. *B* decreases and *dc* increases around the shallower and deeper end of seismogenic zone (Hillers et al., 2006). *dc* is large enough in conditional asperity so as to become conditionally stable (Boatwright & Cocco, 1996). The large value of *dc* may be consistent with the idea that plate boundary is a fairly matured fault composed by thick shear zone and *dc* is larger for the thicker shear zone width in direct shear experiments of gouge layer (Marone, 1998). Although existence of thick shear zone is not directly observed in subduction zone, low velocity layer is found on plate boundary in Japan trench area (Miura et al., 2003; Takahashi et al., 2000). So called subduction channel (Vannucchi et al., 2008) is also a candidate of substance for a large *dc* fault.

As shown in Figure 1, the entire seismogenic zone is modeled as one large hyper asperity. Of course, this is an over-simplification, and it results in simple recurrence of the entire seismogenic zone rupture for an M9 event, a virtually homogeneous slip distribution in the strike direction, a positive stress drop in the entire M9 source area, no M<9 aftershock occurrence in and just around the source area, and so on. To introduce such realistic complexity associated with an M9 event, heterogeneity both in *A* − *B* and *dc* should be considered.

In the seismogenic zones, where *A* − *B* is negative, we set two regular asperities. *dc* is one order smaller in regular asperities than the conditional asperity. The model is a simple and conceptual one and is not intended for any specific subduction zone. The distribution and size of the regular asperities significantly affect the entire rupture patterns and the recurrence time intervals. However, qualitative characteristics, such as M<9 event occurrence within the M9 source area shown below, can be reproduced in many cases that have a different size and distribution of regular asperities. In such cases, *dc* in the regular asperity should be one order or more smaller than it is in the conditional asperity. This order difference in *dc* on a fault is consistent with the discreteness in distribution function for *Gc* (Fukao & Furumoto, 1985; Ide & Aochi, 2005).

### **3. Results**

4 Will-be-set-by-IN-TECH

earthquake size (Ide & Beroza, 2001). In the following sections we set up a numerical model for a hierarchical asperity model and perform numerical simulations. Then we demonstrate that our model is a candidate to account for the earthquake sequence of M=7∼8 within M∼9

In our model, we devide a flat fault, which represents a subducting plate, into many triangular fault cells and calculate variables for each cell. The governing equations are expressed as

*Kij*(*uj* <sup>−</sup> *Vplt*) <sup>−</sup> *<sup>μ</sup>iσ<sup>i</sup>* <sup>−</sup> *<sup>G</sup>*

where *uj* is a fault slip of *j*-th fault cell, *Kij* is a shear stress responce of the *i*-th cell caused by a unit dislocation of the *j*-th cell, *σ<sup>i</sup>* is an effective normal stress of the *i*-th cell, *μ<sup>i</sup>* is a coefficient of friction of the *i*-th cell, *Vpl* is a plate convergence velocity that is assumed constant over the entire model region, *β* is a shear wave velocity and *G* is a rigidity. To calculate the shear stress response at the center of each cell, which is triangular, stress changes calculated with angular dislocations are summed (Comninou & Dundurs, 1975). The last term on the right hand side is a radiation damping term to approximate energy radiation as elastic wave for high velocity

The second equation is a laboratory-derived rate- and state-dependent friction law as (1),

 *Vi V*∗ 

reference for frictonal coefficients, *V*<sup>∗</sup> is a reference velocity and is set to *Vpl* in this study, and *ai*, *bi*, *dci* are friction parameters for the *i*-th cell. The state variable *θ<sup>i</sup>* follows a state evolution

where *Vc* is a cut-off velocity and is 1.0 <sup>×</sup> <sup>10</sup>−<sup>8</sup> m/s (Kato & Tullis, 2001). This evolution law

From the above equations, we derive set of ordinary differential equations for slip velocity *Vi* and state variable *θi*. The governing equations are solved with adaptive time step 4th-order Runge-Kutta algorithm (Press et al., 1996). For initial condition, slip velocity *Vi* and *θ<sup>i</sup>* are

We show spatial distributions of frictional parameters *A* − *B* (*A* = *aσ*, *B* = *bσ*) and *dc* in Figure 1. From rock friction experiments, *A* − *B* primarily depends on temperature and hence depends on depths, and *dc* increases where temperature is high under wet condition (Blanpied et al., 2001). For simplicity, *A* is assumed constant on the entire fault. Parameters *B* and *dc*

+ *bi* ln

*dt* is a slip velocity of the *i*-th cell, *θ<sup>i</sup>* is a state variable for the *i*-th cell, *μ*<sup>∗</sup> is a

 ln *θiVi dci*

 *θiVi dci*

 *<sup>V</sup>*∗*θ<sup>i</sup> dci*

*μ<sup>i</sup>* = *μ*<sup>∗</sup> + *ai* ln

*dθ<sup>i</sup> dt* <sup>=</sup> exp

**2.3 Model setup in subduction zone**

approximates slip law for *V* >> *Vc* and aging law for *V* << *Vc*.

assumed to be uniform, and set to be 0.9*Vpl* and *dci*/*Vi*, respectively.

law. Among several versions of evolution laws, we employ a composite law

 − *Vi Vc* − 2*β dui*

*dt* , (4)

, (5)

, (6)

following. The first equation is a quasi-dynamic equilibrium of shear stress

0 = ∑ *j*

source area.

slip (Rice, 1993).

where *Vi* = *dui*

**2.2 Governing equations**

#### **3.1 Slip history and coseismic slip distribution**

In order to see total earthquake cycle behavior, slip history at some sampling points is shown first. The points are on the same depth as indicated in Figure 1. Within regular asperities, earthquakes occur repeatedly (Figure 2). The slip amount is several times larger in the entire area rupture (giant earthquakes, EQ0, 5) than in the cases of rupture for each regular asperity (regular earthquakes, EQ1-4). The distribution of coseismic slip, whose velocity is higher than 0.01 m/s, is shown in Figure 3. Each moment magnitude is also shown. The rupture areas of EQ1-4 roughly correspond to the two regular asperities. Including both asperities, EQ5 ruputres almost the entire seismogenic zone.

The recurrence time interval is much longer after the giant earthquake than after the regular earthquakes (Figure 2). This recurrence pattern variation is similar to the one in time-predictable model (Shimazaki & Nakata, 1980). Large and small slips correspond to ruptures of a hyper asperity and a regular asperity, respectively (Figure 3). The hyper asperity ruptures with about 1200 years recurrence interval and its moment magnitude is 9.0. During this hyper-cycle, only in the latter stage of the cycle, two regular asperities within the hyper asperity rupture independently. The recurrence interval for the regular cycles is about 200 ∼ 300 years. The superposition of the hyper- and regular-cycles results in the time-predictable recurrence pattern in our model.

On the other hand, in the conditional asperity (*x*=100, 300, 500 in Figure 2), nearly constant aseismic slip and afterslip appear, as well as large coseismic slip and a locked state, depending upon the stage in the M9 earthquake cycle. Such slip pattern variation can also be reproduced in the velocity-strengthening area with constant *dc* (Kaneko et al., 2010). Our model and their model are end-members of the heterogeneity in *A* − *B* and *dc*. We believe that the actual subduction plate boundary has the properties between them.

During a giant earthquake coseismic slip occurs and then locked after long period more than 500 years. After the occurrence of EQ1 or 2, slow sliding appears whose slip velocity is 1/3 ∼ 2/3 of the plate convergence rate. While after EQ3 and/or 4, afterslip occurs and then slow steady slip occurs. Note that interseismic slow sliding speed is almost constant for more than 100 years (Figure 2).

#### **3.2 Space-time evolutions of slip velocity and shear stress**

Snap shots of spatio-temporal variation in slip velocity and stress are shown in Figure 4. The entire area of the hyper asperity is locked after a giant earthquake. A snapshot of slip velocity for 47 years after the giant earthquake is shown in Figure 4a, and the timing is shown in Figure 2. During the interseismic period, the locked area shrinks and slow sliding area spreads. As shown in Figure 4b, stress becomes higher along the edge of the locked area. When the slow sliding area reaches the bottom of a regular asperity within the hyper asperity, rupture starts in the regular asperity (Figure 4b, EQ1). The rupture propagates within the regular asperity but decelerates at its outside (Figure 4c). After the event, low velocity slip like afterslip occurs around the asperity, especially in its deeper extent (Figure 4d). Another regular earthquake (Figure 3b, EQ2) occurs at the other regular asperity (Figure 4e). Then both asperities and their surrounding area are locked again (Figure 4f) but the total locked area is smaller than the stage after the giant earthquake (Figure 4a). Similar size of regular earthquakes (Figure 3c, d, EQ3, 4) and following afterslip occurs after stable slip reaches again the regular asperity (Figure 4g, h, i, j). In this case, afterslip more widely spreads than after the former regular earthquakes (see Figure 4d and i). The locked area becomes further shurunk (Figure 4k). When stable slip reaches the regular asperity at the third time, stress level of around the regular asperity is fairly high (Figure 4l). Thus rupture does not decelerate outside the regular asperity and becomes a giant earthquake (Figure 4 m, Figure 3e, EQ5). After the giant earthquake, the shallower part and in and around the regular asperities are locked first, and afterslip occurs in the surrounding area (Figure 4 n). More than a year after EQ5, whole the seismogenic zone is locked and afterslip continues in the deeper extent (Figure 4 o). Then the entire area is locked again like Figure 4a. As shown in Figure 4 b, g, l, rupture initiates near the deeper edge of a regular asperity for all the cases.

#### **4. Discussion**

#### **4.1 Mechanism of slip variation**

#### **4.1.1 Strength in rate-state friction**

To consider the mechanism of space-time variation in slip pattern as shown above, we use the fault strength introduced by Nakatani (Nakatani, 2001). First, we introduce an alternative state variable.

6 Will-be-set-by-IN-TECH

300 years. The superposition of the hyper- and regular-cycles results in the time-predictable

On the other hand, in the conditional asperity (*x*=100, 300, 500 in Figure 2), nearly constant aseismic slip and afterslip appear, as well as large coseismic slip and a locked state, depending upon the stage in the M9 earthquake cycle. Such slip pattern variation can also be reproduced in the velocity-strengthening area with constant *dc* (Kaneko et al., 2010). Our model and their model are end-members of the heterogeneity in *A* − *B* and *dc*. We believe that the actual

During a giant earthquake coseismic slip occurs and then locked after long period more than 500 years. After the occurrence of EQ1 or 2, slow sliding appears whose slip velocity is 1/3 ∼ 2/3 of the plate convergence rate. While after EQ3 and/or 4, afterslip occurs and then slow steady slip occurs. Note that interseismic slow sliding speed is almost constant for more than

Snap shots of spatio-temporal variation in slip velocity and stress are shown in Figure 4. The entire area of the hyper asperity is locked after a giant earthquake. A snapshot of slip velocity for 47 years after the giant earthquake is shown in Figure 4a, and the timing is shown in Figure 2. During the interseismic period, the locked area shrinks and slow sliding area spreads. As shown in Figure 4b, stress becomes higher along the edge of the locked area. When the slow sliding area reaches the bottom of a regular asperity within the hyper asperity, rupture starts in the regular asperity (Figure 4b, EQ1). The rupture propagates within the regular asperity but decelerates at its outside (Figure 4c). After the event, low velocity slip like afterslip occurs around the asperity, especially in its deeper extent (Figure 4d). Another regular earthquake (Figure 3b, EQ2) occurs at the other regular asperity (Figure 4e). Then both asperities and their surrounding area are locked again (Figure 4f) but the total locked area is smaller than the stage after the giant earthquake (Figure 4a). Similar size of regular earthquakes (Figure 3c, d, EQ3, 4) and following afterslip occurs after stable slip reaches again the regular asperity (Figure 4g, h, i, j). In this case, afterslip more widely spreads than after the former regular earthquakes (see Figure 4d and i). The locked area becomes further shurunk (Figure 4k). When stable slip reaches the regular asperity at the third time, stress level of around the regular asperity is fairly high (Figure 4l). Thus rupture does not decelerate outside the regular asperity and becomes a giant earthquake (Figure 4 m, Figure 3e, EQ5). After the giant earthquake, the shallower part and in and around the regular asperities are locked first, and afterslip occurs in the surrounding area (Figure 4 n). More than a year after EQ5, whole the seismogenic zone is locked and afterslip continues in the deeper extent (Figure 4 o). Then the entire area is locked again like Figure 4a. As shown in Figure 4 b, g, l, rupture initiates

To consider the mechanism of space-time variation in slip pattern as shown above, we use the fault strength introduced by Nakatani (Nakatani, 2001). First, we introduce an alternative

recurrence pattern in our model.

100 years (Figure 2).

**4. Discussion**

**4.1 Mechanism of slip variation 4.1.1 Strength in rate-state friction**

subduction plate boundary has the properties between them.

**3.2 Space-time evolutions of slip velocity and shear stress**

near the deeper edge of a regular asperity for all the cases.

$$
\Theta\_{\dot{l}} = b\_{\dot{l}} \ln \left( \frac{V\_{\ast} \theta\_{\dot{l}}}{d\_{c\dot{l}}} \right) . \tag{7}
$$

With this state variable, the constitutive equation (5) becomes

$$\frac{\tau\_{\dot{l}}}{\sigma\_{\dot{l}}} = (\mu\_\* + \Theta\_{\dot{l}}) + a\_{\dot{l}} \ln \left( \frac{V\_{\dot{l}}}{V\_\*} \right) \tag{8}$$

and the state evolution equation (6) becomes

$$\frac{d\Theta\_{\rm i}}{dt} = \frac{b\_{\rm i}}{(d\_{\rm ci}/V\_{\ast})} \exp\left(-\frac{\Theta\_{\rm i}}{b\_{\rm i}}\right) \exp\left(-\frac{V\_{\rm i}}{V\_{\rm c}}\right) - \frac{V\_{\rm i}}{d\_{\rm ci}} \left\{\Theta\_{\rm i} - b\_{\rm i} \ln\left(\frac{V\_{\ast}}{V\_{\rm i}}\right)\right\}.\tag{9}$$

Here, the equation (8) can be rewritten as follows

$$V\_{\bar{l}} = V\_\* \exp\left[\frac{\tau\_{\bar{l}} - (\mu\_\* + \Theta\_{\bar{l}})\sigma\_{\bar{l}}}{a\_{\bar{l}}\sigma\_{\bar{l}}}\right].\tag{10}$$

If we define a variable (*μ*<sup>∗</sup> + Θ*i*)*σ<sup>i</sup>* as strength of the fault surface, equation (10) gives a slip velocity *Vi* of a fault surface whose strength is (*μ*<sup>∗</sup> + Θ*i*)*σ<sup>i</sup>* when a shear stress *τ<sup>i</sup>* is applied as shown in Figure 5a (Nakatani, 2001). Hence the state evolution equation (9) describes the evolution of the fault strength. The first term of the right hand side represents the healing process of strength which is proportional to log(*t*) (Figure 5b), and the second term represents the decrease of strength in terms of slip (i.e., slip-weakening) as shown in Figure 5c. It should be noted that the strength weakening of the rate- and state-friction law is not rate-weakening but slip-weakening. Rate dependence appears in the steady-state strength <sup>Θ</sup>*SSi* <sup>=</sup> *bi* ln *<sup>V</sup>*<sup>∗</sup> *Vi* 

#### **4.1.2 Mechanism of large variation in coseismic slip amount**

The amount of coseismic slip is determined by how the slip is stopped or decelerated. As described in the fault constitutive law (10), slip velocity depends on the difference between stress and strength relation. Here, the stress is the elastic stress ∑*<sup>j</sup> Kij*(*uj* − *Vplt*) minus radiation damping *μiσ<sup>i</sup>* (4). We compare the variation in the elastic stress and strength near the edge of the right asperity in Figure 1 between the cases of giant and regular earthquakes (EQ5 and EQ1, 3). As shown in Figure 6a, both the stress and the strength decrease in 8-10 MPa during less than 2 meter slip. And the amount of the stress drop is almost the same for regular earthquakes and the giant earthquake. Thus the difference in slip amount is not due to the difference in the stress drop. Although the strength decreases rapidly during about 2m-slip in both cases, the stress decreases more gently in the giant earthquake than in the regular one. This lower weakening rate of the stress keeps the slip velocity of the giant earthquake higher than that of the regular earthquake. Thus the coseismic slip amount of the former is larger than the latter.

What controls such differences in the weakening behavior? Weakening rate of the strength at a point depends mainly on the own slip rate (9), while the stress weakening rate is affected not only own slip history but also the slip distribution around it (4). Slip weakening rate of the stress should be low when slip occurs in wider area around a monitoring point because in equation (4) *Kii* < 0 and *Kij* > 0 (*i* �= *j*). Thus, the difference in stress weakening rate comes from the extent of the slip area. The slip distribution just before regular earthquakes (Figure 4b and g) shows that asperities and their surrounding area, especially above 15km in depth, are almost locked. Since the stress is much lower than the strength in such a locked area (equation 10), large stress increase is necessary for rupture propagation. Until such an area starts to slip, the stress weakening rate in the already sliding area should be high. Thus, the slip is decelerated in the asperity before enough slip occurs in the surrounding area. This is the reason why slip amount is small in regular earthquakes.

On the other hand, slip distribution just before the giant earthquake shows that only the small portion of the entire seismogenic zones including the asperities is locked, and slow sliding occurs around the asperities with around half of the plate convergence rate (Figure 4l). In such sliding area, the stress is close to the strength. Thus, slip around the asperities can be easily triggerd with low stress increase and slip can occur in wider area.

#### **4.1.3 Mechanism for variation in slip velocity among earthquake cycles**

At the middle point between the regular asperities, not only seismic slip but also slow sliding and afterslip occurred (Figure 2). During the slow sliding, the strength is slightly higher than the stress and the weakening is almost negligible (Figure 6b). After the regular earthquakes (EQ3, 4), the stress increases and afterslip occurs here. During the afterslip, strength weakening occurs with low weakening rate because of the slower slip velocity. Because the strength does not significantly reduce, the slip is decelerated. Hence both the amount of stress drop and slip amount of the afterslip are significantly smaller than those of the giant earthquake. During the giant earthquake, typical slip weakening occurs as similar to the one at the regular asperity (Figure 6a). The amount of the stress drop is also similar. Only the difference is the weakening rate. Because this point has much larger *dc* than in the regular asperity, the slip weakening distance is about 12 m, which is 6 times more than the regular asperity. As explained in 4.1.2, because wide area slips in the giant earthquake, weakening rate of the stress is low enough and slip can continue there.

#### **4.1.4 A mechanism of time predictable behavior**

As shown in Figure 4, the area around the asperities is locked after the giant earthquakes but slowly sliding after the regular earthquakes. This causes the difference in stressing rate around the hypocenter. After the giant earthquake, stressing rate becomes low because the wide area is locked (Figure 7). On the other hand, the stressing rate is high after the regular earthuqkes because the surrounding area is sliding (Figure 7). Therefore, such inherent stressing rate variation during the cycle of the giant earthquake causes the variation in recurrence interval depending on the earthquake size. This mechanism of time predictable behavior is completely different with the original time predictable model in which the stressing rate and the peak strength are constant and the stress drop varies depending on the earthquake size (Shimazaki & Nakata, 1980). It should be noted that the time predictable behavior here does not necessarily appear in general. If regular asperities exist in the deeper portion of the seismogenic zone, they should be ruptured earlier because of the slow sliding propagation from the bottom of the seismogenic zone.

#### **4.2 Comparison with observation**

8 Will-be-set-by-IN-TECH

(Figure 4b and g) shows that asperities and their surrounding area, especially above 15km in depth, are almost locked. Since the stress is much lower than the strength in such a locked area (equation 10), large stress increase is necessary for rupture propagation. Until such an area starts to slip, the stress weakening rate in the already sliding area should be high. Thus, the slip is decelerated in the asperity before enough slip occurs in the surrounding area. This

On the other hand, slip distribution just before the giant earthquake shows that only the small portion of the entire seismogenic zones including the asperities is locked, and slow sliding occurs around the asperities with around half of the plate convergence rate (Figure 4l). In such sliding area, the stress is close to the strength. Thus, slip around the asperities can be

At the middle point between the regular asperities, not only seismic slip but also slow sliding and afterslip occurred (Figure 2). During the slow sliding, the strength is slightly higher than the stress and the weakening is almost negligible (Figure 6b). After the regular earthquakes (EQ3, 4), the stress increases and afterslip occurs here. During the afterslip, strength weakening occurs with low weakening rate because of the slower slip velocity. Because the strength does not significantly reduce, the slip is decelerated. Hence both the amount of stress drop and slip amount of the afterslip are significantly smaller than those of the giant earthquake. During the giant earthquake, typical slip weakening occurs as similar to the one at the regular asperity (Figure 6a). The amount of the stress drop is also similar. Only the difference is the weakening rate. Because this point has much larger *dc* than in the regular asperity, the slip weakening distance is about 12 m, which is 6 times more than the regular asperity. As explained in 4.1.2, because wide area slips in the giant earthquake, weakening

As shown in Figure 4, the area around the asperities is locked after the giant earthquakes but slowly sliding after the regular earthquakes. This causes the difference in stressing rate around the hypocenter. After the giant earthquake, stressing rate becomes low because the wide area is locked (Figure 7). On the other hand, the stressing rate is high after the regular earthuqkes because the surrounding area is sliding (Figure 7). Therefore, such inherent stressing rate variation during the cycle of the giant earthquake causes the variation in recurrence interval depending on the earthquake size. This mechanism of time predictable behavior is completely different with the original time predictable model in which the stressing rate and the peak strength are constant and the stress drop varies depending on the earthquake size (Shimazaki & Nakata, 1980). It should be noted that the time predictable behavior here does not necessarily appear in general. If regular asperities exist in the deeper portion of the seismogenic zone, they should be ruptured earlier because of the slow sliding

is the reason why slip amount is small in regular earthquakes.

rate of the stress is low enough and slip can continue there.

**4.1.4 A mechanism of time predictable behavior**

propagation from the bottom of the seismogenic zone.

easily triggerd with low stress increase and slip can occur in wider area.

**4.1.3 Mechanism for variation in slip velocity among earthquake cycles**

#### **4.2.1 M**<**9 earthquakes in M9 source area**

As shown in Figure 3, M<9 earthquake occurrence in M9 source area is reproduced. The slip amount of M7 events is nearly one order smaller than that of M9 event at the regular asperites (Figure 3). This is consistent with the slip amount difference between the 1978 off Miyagi earthquake and the 2011 Tohoku earthquake as mentioned in Introduction.

In our simulation, M7 earthquakes occur only in the later half of the hyper cycle. However, they occur decades after the M9 earthquakes in Kamchatka, Alutian and Colombia-Ecuador. For Kamchatka, the time interval between M<9 and M9 earthquakes are 41 years, which is around 20% of the recurrence interval of 215 years. The M<9 earthquake occurrence timing discrepancy after M9 earthquakes between the data and our simulation results can be explained as in 4.1.4. Actually, the above mentioned M<9 earthquake sources located deeper part of the M9 source areas. This indicates that the 1978 off Miyagi type earthquake also can occur in the early stage of M9 cycle because its asperity is located at the deeper edge of the 2011 Tohoku earthquake (Kato & Yoshida, 2011). It should be noted that such M<9 earthquakes are not usual aftershocks just after M9 earthquakes but occurred after more than a few decades.

#### **4.2.2 Interseismic period**

As shown in Figure 2, aseismic sliding with the velocity of 30 ∼ 60 % of the plate convergence rate occurs before the M9 earthquakes and seismic slip occurs during the M9 earthquakes. This result seems to be consistent with the sliding behavior in the southern half of the 2011 Tohoku earthquake. In the area, aseismic sliding occurred with the velocity of less than 50% of the plate convergence rate (Uchida et al., 2009). In this area, coseismic slip also occurs in the 2011 Tohoku earthquake but the slip amount is relatively small (Iinuma et al., 2011). Furthermore, in the aseismically sliding area, afterslip after M=6∼7 earthquakes extended significantly wider than the coseismic slip area (Suito et al., 2011). Such large area afterslip can be found in the later stage of the M9 earthquake as shown in Figure 4i. Hence, such afterslip occurrence is possibly an indicator of the high stress level before the greater event.

#### **4.2.3 Preseismic slip and rupture initiation**

It should be noted here that the wide area afterslip mentioned above decelerates and locked area appears again (Figure 4j). In other words, the afterslip does not accelerate to produce the M9 earthquake as a preslip. As shown in Figure 4l, the initiation of the M9 rupture is triggered by the rupture of a regular asperity, and the nucleation process is similar to the one for the M7 earthquakes in our simulation (Figure 4b and g). This indicates that large preslip is not necessary observed. Actually, preseismic acceleration of the slip near the hypocenter of the 2011 Tohoku earthquake has not been found (Hirose, 2011). On the other hand, it is important to examine whether the initial rupture of the 2011 Tohoku earthquake corresponds to an asperity identified before in this area or not.

#### **4.2.4 Coseismic period**

One of the characteristics of the 2011 Tohoku earthquake is the several tens meter slip along the Japan trench axis (Maeda et al., 2011). As shown in Figure 3, the large slip in the shallower part of the fault is reproduced. Since the fault plane cut the free surface, there is no loading source shallower than the fault plane. This causes the accumulation of the slip deficit easier than the deeper portion of the seismogenic zone which is loaded due to the sliding of the deeper extent. Additionally, since the state evolution length *dc* is large, slip cannot easily be triggered by M7 earthquakes. Furthermore, when the seismic slip occurs in the M9 earthquake, free surface enhance the sliding. These may explain the observed large slip near the trench axis.

Since our simulation is quasi-dynamic, we cannot estimate rupture propagation velocity quantitatively. However, large *dc* in the conditional asperity results in the lower weakening rate and hence slower rupture propagation than the ordinal earthquakes. This is consistent with the slow rupture velocity estimation for the 2011 Tohoku earthquake (Wang & Mori, 2011). Although the weakening rate during the coseismic rupture has not been estimated yet, our model predicts that higher and lower weakening rate in the regular and conditional asperities, respectively (Figure 6). This is another verification of our hierarchical asperity model and friction law used here.

#### **4.2.5 Postseismic period**

Furthermore, our model predicts the postseismic transient behavior which can be tested by the observation during the coming years after the 2011 Tohoku earthquake. As shown in Figure 4n, the asperities of the past M=7∼8 events will be locked first, and afterslip occurs in the surrounding area. Then, after a year, the entire seismogenic zone will be locked and afterslip will continue only in the deeper extent. Such spatio-temporal variation in the locked zone after the 2011 Tohoku earthquake can be examined by repeating earthquake distribution and more directly by seafloor deformation observed using ocean bottom pressure gauges off Miyagi. Long term afterslip in the deeper portion will cause uplift along the coast where subsidence occurred seismically. This will be observed by GPS onland for decades.

#### **4.2.6 Other subduction zones**

The most important implication from our results and observation is the wide area afterslip observed in the later stage of the M9 earthquake cycle as shown in 4.2.2. Such wide area afterslips comparing with coseismic slip areas are observed in Hyuga-nada, southwest Japan (Yagi et al., 2001). Although no interplate coupling or even forward slip has been obtained in the southern part of this area, a M=7.7 earthquake occurred in 1662, causing strong motion and a tsunami in wide areas of Hyuga (Hatori, 1969). The area off Hyuga may now be in the later stage of a cycle of M∼8 earthquake such as in 1662, while the Japanese Government gives only 10% chance of such an earthquake occurring in this area (Headquarters for Earthquake Research Promotion, 2004). Furthermore, this area may be ruptured with great earthqaukes along the Nankai trough.

Another example of wide area afterslip is observed in Tokachi-Nemuro area, northeast Japan, after the 2003 Tokachi-oki earthquake of Mw∼8 (Miyazaki et al., 2004). In 1600's, significantly larger tsunami than the one caused by 2003 earthquake attacked the coast of Tokachi-Nemuro area (Nanayama et al., 2003). The estimated source area includes the 2003 coseismic slip area, the 1973 Nemuro coseismic slip area (Mw=7.8), and the wide afterslip area. Note that most of the recurrence interval estimated from the tsunami deposits in this area is less than 500 years (Sawai et al., 2009).

of the fault is reproduced. Since the fault plane cut the free surface, there is no loading source shallower than the fault plane. This causes the accumulation of the slip deficit easier than the deeper portion of the seismogenic zone which is loaded due to the sliding of the deeper extent. Additionally, since the state evolution length *dc* is large, slip cannot easily be triggered by M7 earthquakes. Furthermore, when the seismic slip occurs in the M9 earthquake, free surface enhance the sliding. These may explain the observed large slip near the trench axis. Since our simulation is quasi-dynamic, we cannot estimate rupture propagation velocity quantitatively. However, large *dc* in the conditional asperity results in the lower weakening rate and hence slower rupture propagation than the ordinal earthquakes. This is consistent with the slow rupture velocity estimation for the 2011 Tohoku earthquake (Wang & Mori, 2011). Although the weakening rate during the coseismic rupture has not been estimated yet, our model predicts that higher and lower weakening rate in the regular and conditional asperities, respectively (Figure 6). This is another verification of our hierarchical asperity

Furthermore, our model predicts the postseismic transient behavior which can be tested by the observation during the coming years after the 2011 Tohoku earthquake. As shown in Figure 4n, the asperities of the past M=7∼8 events will be locked first, and afterslip occurs in the surrounding area. Then, after a year, the entire seismogenic zone will be locked and afterslip will continue only in the deeper extent. Such spatio-temporal variation in the locked zone after the 2011 Tohoku earthquake can be examined by repeating earthquake distribution and more directly by seafloor deformation observed using ocean bottom pressure gauges off Miyagi. Long term afterslip in the deeper portion will cause uplift along the coast where

The most important implication from our results and observation is the wide area afterslip observed in the later stage of the M9 earthquake cycle as shown in 4.2.2. Such wide area afterslips comparing with coseismic slip areas are observed in Hyuga-nada, southwest Japan (Yagi et al., 2001). Although no interplate coupling or even forward slip has been obtained in the southern part of this area, a M=7.7 earthquake occurred in 1662, causing strong motion and a tsunami in wide areas of Hyuga (Hatori, 1969). The area off Hyuga may now be in the later stage of a cycle of M∼8 earthquake such as in 1662, while the Japanese Government gives only 10% chance of such an earthquake occurring in this area (Headquarters for Earthquake Research Promotion, 2004). Furthermore, this area may be ruptured with great earthqaukes

Another example of wide area afterslip is observed in Tokachi-Nemuro area, northeast Japan, after the 2003 Tokachi-oki earthquake of Mw∼8 (Miyazaki et al., 2004). In 1600's, significantly larger tsunami than the one caused by 2003 earthquake attacked the coast of Tokachi-Nemuro area (Nanayama et al., 2003). The estimated source area includes the 2003 coseismic slip area, the 1973 Nemuro coseismic slip area (Mw=7.8), and the wide afterslip area. Note that most of the recurrence interval estimated from the tsunami deposits in this area is less than 500 years

subsidence occurred seismically. This will be observed by GPS onland for decades.

model and friction law used here.

**4.2.5 Postseismic period**

**4.2.6 Other subduction zones**

along the Nankai trough.

(Sawai et al., 2009).

Fig. 1. Spatial distributions of frictional parameters (a) *A* − *B* and (b) *dc* on the model fault. Crosses indicate the cumulative slip shown in Figure 2.

Fig. 2. Temporal variation of slip at the points shown in Figure 1. *x* and *d* are strike and depth in Figure 1. *Vpl* is the plate convergence rate (0.05m/yr). Arrows indicate the timing of images shown in Figure 4.

Fig. 2. Temporal variation of slip at the points shown in Figure 1. *x* and *d* are strike and depth in Figure 1. *Vpl* is the plate convergence rate (0.05m/yr). Arrows indicate the timing

of images shown in Figure 4.

Fig. 3. Coseismic slip distribution and the moment magnitude for each events.

Fig. 4. Spatio-temporal variations in the stress and slip velocity on the fault. Blue, white, yellow-green and red colors for the slip velocity correspond to locking, slip with plate convergence rate, aseismic slip (faster than the plate convergence rate) and seismic slip, respectively. Numerals attached to the arrows indicate the time intervals between two successive snapshots.

Fig. 4. Spatio-temporal variations in the stress and slip velocity on the fault. Blue, white, yellow-green and red colors for the slip velocity correspond to locking, slip with plate convergence rate, aseismic slip (faster than the plate convergence rate) and seismic slip, respectively. Numerals attached to the arrows indicate the time intervals between two

successive snapshots.

Fig. 5. (a) Relation between stress, strength and slip velocity in the rate- and state-dependent friction law. (b) Time dependent variation in strength. (c) Slip dependent variation in strength.

Fig. 6. Solid, dashed and thick dashed lines indicate stress, strength and slip velocity variation, respectively. (a) in a regular asperity. (b) at a point between the two asperities.

Fig. 7. Stress and strength variation in time at the hypocenter. Solid and dashed lines indicate stress and strength variation, respectively. Coseismic period is excluded.
