The effect of thermal pressurization on dynamic fault branching

. Finally, we considered the case where there is a free surface above the branch fault system. In this case, the rupture can propagate along both faults because of interaction between the free surface and the branch fault, in addition to TP on the main fault.


I N T RO D U C T I O N
One important issue with regard to future large earthquakes is the assessment of strong ground motions and tsunamis. To accurately evaluate seismic and tsunami hazards, we must have a fundamental understanding of the physics of the dynamic rupture propagation along faults with complex geometry because, as geological and seismic observations have indicated, large earthquakes generally occur on complex fault systems including jogs, step-overs, bends and branches rather than a single planar fault (e.g. Sibson 1986;Wesnousky 1988;Sieh et al. 1993). Here, we focus on fault branching. Branch faults have been observed in shallow inland strikeslip earthquakes (e.g. Archuleta 1984;Aochi & Fukuyama 2002). Recently, branch faults have also been observed along subduction interfaces at shallow depths by seismic reflection surveys, such as in the Nankai subduction zone off southwest Japan (e.g. Park et al. 2002;Moore et al. 2007), in the Japan Trench (Tsuji et al. 2013), and in the north Ecuador-south Colombia oceanic margin (Collot et al. 2008).
Rupture propagation on branch fault systems has been studied both theoretically and numerically. Theoretical analyses suggested that rupture velocity and pre-stress states control rupture paths (Poliakov et al. 2002). These features have also been investigated using numerical simulations. Kame et al. (2003) performed many numerical simulations of 2-D mode II ruptures in an infinite elastic medium. They showed that rupture propagation paths are controlled by the inclination of the maximum principal stress, the branching angle, and incoming rupture velocity at the branching intersection. In addition, Bhat et al. (2007) suggested that the finite length of the branching fault could affect the rupture dynamics. Numerical models of earthquake ruptures including branch faults have mainly focused on strike-slip events (Aochi & Fukuyama 2002;Oglesby et al. 2003;Bhat et al. 2004;Fukuyama & Mikumo 2006). To model a megathrust earthquake rupture, Tamura & Ide (2011) introduced material heterogeneity and a free surface. They considered a spontaneously propagating mode II crack on a bi-material interface and a branching fault, suggesting that material contrast plays an important role in dynamic fault branching.
In subduction zones, the existence of fluid has been suggested (e.g. Hasegawa et al. 2009), and fluid at the plate interface is thought to play an important role during rupture propagation. As one possible effect of the fluid's presence, Sibson (1973Sibson ( , 1977 proposed thermal pressurization (TP), an increase in pore fluid pressure due to frictional heating in wet conditions, that results in a decrease in effective normal stress. Field observations and rock experiments have suggested that TP can be in effect during earthquakes. Ishikawa et al. (2008) analysed the compositions of major and trace elements and isotope ratios of borehole core samples from the Chelungpu fault, Taiwan, and suggested the coseismic presence of high-temperature fluids (>350 • C). They mentioned that TP may have been in effect during the M w 7.6 1999 Chi-Chi earthquake because the fault's hydraulic diffusivity was estimated to be low (Doan et al. 2006). Ujiie et al. (2007) studied the ancient subduction thrust faults in the Shimanto accretionary complex of southwest Japan and observed that the ultracataclasite layer indicates that granular material was injected into the adjacent wall rocks. They suggested that the injected material represents fluidization during a seismic slip and that TP appears to be, of the various fluidization mechanisms, the most consistent with their observations. Ujiie & Tsutsumi (2010) conducted high-velocity friction experiments on clay-rich fault gouge obtained from the megasplay fault zone in the Nankai subduction zone. They observed a rapid increase in the gouge's temperature during rapid slip weakening under wet conditions, which implies that rapid weakening results from TP.
Numerical simulations suggested that TP can increase breakdown stress drop (e.g. Andrews 2002;Bizzarri & Cocco 2006b), peak slip rate and rupture velocity (Bizzarri & Cocco 2006a;Bizzarri & Spudich 2008;Urata et al. 2008). As a result of these effects, TP also affects rupture processes when geometric complexities, such as step-overs, are taken into consideration (Urata et al. 2012). TP can also modify fault behaviour at longer timescales. Noda & Lapusta (2013) proposed a model applicable for the M w 9.0 2011 Tohoku-Oki earthquake in which a creeping unlocked part of the fault can become an unstable seismically slipping patch due to TP's effects. Their model for the time of the earthquake produced higher slip but a smoother slip-rate history in the part of the fault where TP was active, making the model qualitatively consistent with seismic observations. It should be noted that they considered a single fault plane without geometrical complexity and did not consider the free surface. As far as we know, the geometry of shallow thrust faults has never been considered in dynamic rupture simulations with TP, although such geometry is common in subduction zones (e.g. Park et al. 2002;Moore et al. 2007;Collot et al. 2008;Tsuji et al. 2013).
As mentioned above, TP could affect dynamic fault branching phenomena strongly. However, previous studies on branching issues have not taken TP into account. The purpose of this study is to explore the effects of TP on dynamic fault branching. We also examine the effects of the free surface and TP simultaneously. We perform 2-D numerical simulations of spontaneous ruptures with TP, based on the branch fault model of Kame et al. (2003). The principal stresses, the inclination of the maximum principal stress with respect to the main fault, and the branching angle are the same as those of Kame et al. (2003). We also assume that the maximum principal stress axis is horizontal in cases with free surfaces (i.e. the maximum stress axis is parallel to the surface). We investigate whether the rupture propagates along the main and/or branch faults to demonstrate how TP alters dynamic fault branching. Finally, we simulate a case where there is a free surface above the branch fault system.

Model
We consider a 2-D mode II (in-plane) rupture that propagates along a planar main fault and encounters an intersection with a pre-existing branch fault (Fig. 1). The branch fault makes an angle θ with the main fault. The fault system is in an infinite, homogenous, isotropic and linear elastic medium, subjected to uniform external stresses σ max and σ min (compression is taken as positive). The fault planes are on the xz-plane. The direction of σ max is parallel to the x-axis, and it makes an angle ψ with the main fault plane.

Numerical procedure for dynamic ruptures
The numerical algorithm for dynamic ruptures is based on the 2-D boundary integral equation method (BIEM) with integration kernels described by eq. (A1) of Tada & Madariaga (2001). These kernels link both in-plane shear slip (mode II) and opening slip (mode I) finite dislocation to stress tensor change in the medium surrounding the slipping crack, as functions of time. Please note that in this study, opening slip is set to zero for faults, and is computed for free surface elements as vertical displacement (see Appendix A). Our fault model consists of a collection of individual fault elements (any geometry consisting of a combination of small planar fault elements can be considered). Each of them has its own amount of dislocation at a given time step and thus its own contribution to the stress change in the surrounding medium, thanks to the integration kernels. The spontaneous rupture is obtained by solving the loading stress applied to each fault element at each time step together with the prescribed frictional behaviour as a function of its slip history (constitutive law).
A rupture is initiated by artificially adding shear stress in a small patch on the main fault ( Fig. 1). It then proceeds spontaneously, governed by a slip-weakening law with the Coulomb failure criteria where τ is the shear stress, μ s is the static coefficient of friction, σ eff n is the effective normal stress, μ d is the dynamic coefficient of friction, D c is the characteristic displacement and u is the slip. The friction coefficients and D c are assumed to be uniform and the same on the two faults.

Calculation of pore pressure change due to TP
We now consider the effect of frictional heating. Once slip is initiated and the slip rate increases, the pore fluid pressure p increases because of frictional heating. Then, σ eff n decreases according to the relation σ eff n = σ n − p, where σ n is the normal stress. We assume that the pore pressure change p at location x on the fault surface due to TP obeys the following equation (Bizzarri & Cocco 2006a), where the dimensionless parameter γ is α f β f c, w is the shear zone thickness, χ is thermal diffusivity, hydraulic diffusivity ω is k ηβ f φ, k is permeability, φ is porosity, v is slip velocity, ε is an arbitrarily small positive real number and erf ( ) is the error function defined below.
The definitions of α f , β f , c and η are listed in Table 1. The assumptions in eq.
(2) and their validity are set forth in Bizzarri & Cocco (2006a) and Urata et al. (2008). ω and φ are set to be constant. We assumed that fluid migration between the main and branch faults does not affect each other. We discuss the validity of this assumption in Section 4. Note that the yield stress values τ s = μ s σ eff n are not affected by TP because the pore pressure is constant ( p = p 0 ) before the rupture front arrives.

Validation tests
Before computing the dynamic fault branching with TP, we tested our numerical code. For dynamic rupture propagation with TP on one fault plane in an infinite medium, we compare the results with those of the 3-D finite-difference method (FDM) used by Urata et al. (2008). The initial stresses, μ s , μ d and D c are assumed uniform in the comparison. In the 3-D FDM calculations, we set the fault width along the y-direction long enough to stay in a pure 2-D regime and measured the slip evolution at y = 0, which we compared with that measured using the 2-D BIEM (Fig. 3a). Figs 3(b) and (c) show the propagation of the rupture front and the evolution of the slip velocity at the centre of the fault for the BIEM (red curves) and the FDM (black curves). In the simulations, μ s , μ d and D c are 0.74, 0.5 and 0.25 m, respectively. The initial effective normal stress σ eff n0 = σ n0 − p 0 and the initial shear stress τ 0 are 102 and 66 MPa, respectively. The FDM grid spacing ( x, y and z) and time interval are 12.5 m and 0.001 s, respectively. The BIEM grid spacing x 1 and time interval are 12.5 m and 0.0009 s, respectively. The density ρ, P-wave velocity v p and S-wave velocity v s are 2670 kg m -3 , 6.000 km s -1 and 3.464 km s -1 , respectively. The TP parameters are shown in Table 1. As shown in Figs 3(b) and (c), the red and black curves overlap; the solutions obtained by the BIEM and the FDM are almost identical. This comparison validates our calculation of dynamic rupture on a planar fault with TP.
To compare the distribution of off-fault stress changes, we first compute the stresses on the plane parallel to the fault (Fig. 3a). Fig. 4 shows the changes in the shear and normal stresses obtained  (c), because TP is absent. The off-fault stress changes calculated by the two methods are quite similar (Fig. 4). We confirmed that our numerical code properly provides dynamic stress changes on an off-fault plane parallel to a fault plane.
We also tested the off-fault static stress change by comparing our results with those computed by the Okada's static deformation code (Okada 1992). To compute the off-fault static stress change tensor using the 2-D BIEM, we calculated the slip evolution on a fault in full space and measured the stress change after a sufficiently long time had passed. TP is not in effect in this computation. Since we only needed the final stress change, the initial stress on the fault is set higher than the static shear strength. The rupture initiated simultaneously all over the fault and slip evolved according to the slip-weakening law (eq. 1, Fig. 2). The fault length was 50 m. σ eff n0 , τ 0 , μ s , μ d and D c were 100 MPa, 37 MPa, 0.36, 0.1 and 0.008 m, respectively. ρ, v p and v s are the same as those used in Fig. 3. The grid spacing and the time interval were 12.5 m and 0.0009 s, respectively. In Okada's computation, which deals with finite rectangular sources in 3-D half-space, we put the fault width along the y-direction centred at y = 0. The fault width was set long enough to stay in a pure 2-D regime and the fault was located deep enough not to be affected by the free surface. The static stress change tensor was evaluated along the centre of the fault length (i.e. y = 0). We then compared Okada's static stress change tensors and our BIEM computation and confirmed that they were identical. This experiment certifies that our off-fault stress change tensors are correctly computed.
Finally, we compared our dynamic fault branching behaviour in the absence of TP with the results of Kame et al. (2003). We examined cases where the angles ψ between σ max and the main fault are 13 • and 25 • with branching angles θ of +15 • and +30 • . In our computations, the rupture velocities when the rupture front reached the branching point were around 0.6v s and 0.8v s . The final rupture traces of these above eight cases are the same as figs. 12 and 13 of Kame et al. (2003), except for the case with an incoming rupture velocity of 0.8v s , ψ of 13 • and θ of +30 • . In Kame et al. (2003), the rupture for those conditions propagates on both the main and branch faults. In our simulations, however, the rupture propagates only on the main fault, similar to the case with an incoming rupture velocity of 0.6v s . If we increase the incoming rupture velocity to 0.9v s , the rupture propagates on both the main and branch faults. We noted that the conditions that cause a rupture to propagate on both faults are very delicate, and might depend on the rupture initiation. Thus, our rupture propagation paths are slightly different from those of Kame et al. (2003).
As described above, we confirmed that our numerical code works properly.

Parameters for simulations
We considered the four cases where the dip angle of the main fault ψ is either 13 • or 25 • with a branching angle θ of either 15 • or 30 • (Fig. 1, Table 2). Fig. 5 shows the pre-stress state we used for all four cases, which is the same as that used by Kame et al. (2003). We assume that σ eff n0 and τ 0 on the main fault are 100 and 24 MPa, respectively, for all cases. The values of σ eff n0 and τ 0 on the branch faults depend on both ψ and θ (see Fig. 5 for details). The rupture Table 2. Simulation parameters for fault geometry. Note: ψ, θ , L, L 1 and L 2 are defined in Fig. 1.  initiation is the same in all cases. The length of the nucleation patch is 1.0 km and the initial shear stress inside the nucleation patch is τ s + 2.0 MPa. TP is not in effect inside the nucleation patch regardless of whether or not it is in effect on the main fault. The distance from the nucleation patch to the branching intersection is 1.9 km. The rupture velocity when the rupture front reaches the branching point is about 0.6v s for the cases where TP is not in effect.
Tables 1 and 3 list the parameters used in this study. Ujiie et al. (2007) observed that ultracataclasite in an ancient subduction thrust is mostly a few centimetres thick, so the value of w in our simulations (1 cm) is reasonable. Permeability k and porosity φ for the sediments of four different subduction zones ranges from 10 −19 to 10 −15 m 2 and from 0.2 to 0.9, respectively, at effective stress of up to 12.4 MPa (Gamage et al. 2011). We set k at 5 × 10 −17 m 2 and φ at 0.025, which is smaller than the measured values of φ. Our choices are still reasonable because k and φ are found to decrease dramatically with increasing effective stress (e.g. Tanikawa & Shimamoto 2009;Ikari & Saffer 2012), which could be larger in seismogenic conditions than 12.4 MPa of the measurements. If we use a larger value of φ, TP would be more effective because the hydraulic diffusivity ω = k ηβ f φ is proportional to k φ. From the values in Table 1, ω is approximately 10 -2 m 2 s -1 , much larger than the thermal diffusivity χ = 10 −6 m 2 s −1 . Therefore, ω dominates TP's effect (eq. 2).
We assumed μ s to be 0.6, consistent with laboratory experiments using fault gouge from the megasplay fault zone in the Nankai subduction zone (Ujiie & Tsutsumi 2010). The S value in the absence of TP, as defined by Andrews (1976), is set at 3.0 on the main fault. D c is 0.64 m, consistent with seismic observations for strike-slip earthquakes (Mikumo et al. 2003;Fukuyama & Mikumo 2007).

Low main fault dip angle (ψ = 13 • )
First, we will describe the dynamic rupture branching behaviour calculated for a main fault dip angle ψ of 13 • and a branch angle θ of 15 • (Model 1 in Table 2). Fig. 6(a) shows the evolution of slip velocity on the main and branch faults when TP is not in effect on either fault. The rupture initiated on the main fault propagates bilaterally along the main fault. After the rupture encounters the intersection, it terminates. A rupture on the branch fault, however, is triggered and propagates successfully. Thus, the rupture propagates along the branch fault. TP can alter the rupture propagation paths. Fig. 6(b) shows the time history of the slip velocity when TP is in effect only on the main fault. TP increases the slip velocity and the rupture velocity on the main fault, as described by previous numerical studies (e.g. Bizzarri & Cocco 2006a,b;Bizzarri & Spudich 2008;Urata et al. 2008Urata et al. , 2012. In this case, the rupture on the main fault continues to propagate after encountering the branch point. A rupture on the branch is triggered, but it immediately dies out. Due to TP on the main fault, the rupture propagates along the main fault. However, when TP is in effect on both the main and branch faults, the rupture propagates along the branch (Fig. 7a).
The rupture propagation paths depend on the branch angle θ when TP is not in effect on either fault ( Fig. 7; Kame et al. 2003). The rupture propagates along the branch fault in Model 1 (θ = 15 • ; Figs 6a and 7a), but along the main fault in Model 2 (θ = 30 • ; Fig. 7b). The rupture propagation paths, however, do not depend on the branch angle θ when TP is in effect. Ruptures propagate along the main fault in cases with TP on the main fault, and along the branch fault when TP is in effect on both faults (Fig. 7). This is observed in both Models 1 (θ = 15 • ) and 2 (θ = 30 • ). Thus, dynamic rupture processes are affected more strongly by TP than by θ .
Why does TP alter rupture propagation paths? Rupture propagation paths are controlled by the pre-stress state when TP is not in effect on either fault. Fig. 8 shows the distribution of stress   (Model 1). The white colour corresponds to the initial value on the main fault (0.24). The red and blue colours indicate τ/σ eff n values higher and lower than 0.24, respectively. During slips, τ/σ eff n obeys the slip-weakening law (Fig. 2). The black curves indicate the evolution of the rupture front. The dark blue regions at distances of 1.9-2.3 km on the main fault, which appear after the rupture front arrives, show that the slip arrests and that τ/σ eff n no longer obeys the slip-weakening law. The blue lines show the branch intersection. ratio τ/σ eff n in Model 1 (θ = 15 • ). The stress ratio increases on both faults owing to slip on the main fault before the rupture encounters the branch intersection, and ruptures on both faults are triggered. The rupture triggered on the branch fault is more likely to propagate because the pre-stress state on the branch fault is closer to the Coulomb failure criterion than on the main fault in Model 1 (θ = 15 • ; Fig. 5). Once the rupture propagates, the stress shadow (Yamashita & Umeda 1994) of the branch fault, which corresponds to the blue regions in Fig. 8 at distances of 2.5-3.5 km on the main fault, hinders rupture propagation along the main fault. Although in Model 2 (θ = 30 • ; Fig. 5), the pre-stress state on the branch is a little closer to the Coulomb failure criterion than that on the main fault, the fracture energy G c = 1 2 (μ s − μ d )σ eff n0 D c is much higher (solid and dashed lines in Fig. 5), and thus the rupture propagates along the main fault. When TP is in effect on the main fault, the increase in pore pressure (left-hand panel in Fig. 9a) rapidly weakens the shear resistance, promoting rupture on the main fault. A large, quick stress drop promotes the stress shadow effect (right-hand panel in Fig. 9a),  Fig. 8. Dark blue regions at distances of 1.9-2.1 km on the main fault in the right-hand panel of (b) and at distances of 0-0.2 km on the branch fault in the right-hand panel of (c) indicate that the slip arrests and that τ /σ eff n no longer obeys the slip-weakening law.
suppressing rupture propagation on the branch fault. When TP is in effect on both faults, the rupture propagation paths are controlled by the TP's effectiveness. TP on the branch fault is more effective and the pore pressure on the branch fault increases much more than on the main fault (left-hand panel in Fig. 9b). The effective TP on the branch fault results from high initial normal and shear stresses (Fig. 5). The high initial stresses increase frictional heating, and the high initial shear stress can increase the stress drop if TP is in effect. The TP on the branch fault encourages rupture propagation on the branch fault and enhances the stress shadow effect (righthand panel in Fig. 9b). As a result, the rupture propagates along the branch fault.

High main fault dip angle (ψ = 25 • )
The initial stresses on the branch faults are different at ψ values of 13 • and 25 • , although the initial stresses on the main fault are the same (Fig. 5). Therefore, the dynamic ruptures are the same before the rupture front reaches the branching intersection, but the dynamic fault branching is different. When ψ is 25 • , the rupture propagates along the main fault for both Models 3 (θ = 15 • ) and 4 (θ = 30 • ), whether or not TP is in effect. The pre-stress state on the main fault is closer to the Coulomb failure criterion than on the branch faults, and the fracture energy G c of the main fault is smaller than that of the branch faults (Fig. 5), so the main fault is more likely to break. In addition, TP on the branch fault is less effective than when ψ is 13 • . Fig. 9(c) shows that the increase in the pore pressure on the branch fault is small and that the rupture propagates along the main fault even when TP is in effect on both faults.

Free surface effect
We considered the case where a free surface is located 20 m above the faults' upper edge. Appendix A details the numerical method we employed. It should be noted that full-space and half-space computations are the same until the reflected waves from the free surface arrive at the branch fault system. However, once the reflected waves from the free surface arrive, they change the stress field and thus alter the slips on the faults. We examined cases where the dip angle of the main fault ψ is 13 • with a branching angle θ of 15 • , with TP in effect and not in effect on the main fault. Please note that TP is not in effect on the branch.
When TP is not in effect on either of the faults, the behaviour of the dynamic rupture branching is similar to the case without the free surface: the rupture propagates along the branch fault and the stress shadow effect hinders the rupture from propagating along the main fault. In contrast, when TP is in effect on the main fault, the free surface alters the rupture propagation path. Fig. 10(a) shows the slip velocity on both faults. First, the rupture propagates along the main fault. This is also observed when there is no free surface (Fig. 6b). During rupture propagation along the main fault, seismic waves reflected at the free surface arrive at both faults (Fig. 10b). Then, the reflected waves trigger rupture near the centre of the branch fault, and the rupture propagates bilaterally (Fig. 10a). Slip on the branch fault suppresses slip on the main fault (the white regions shown after 1.5 s at distances of 2.0-3.5 km on the main fault in Fig. 10a). Therefore, the rupture finally propagates along both faults because of the free surface and TP on the main fault.
It should be noted that the asymmetric pattern of slip velocity distribution in Fig. 10(a) results from the inclination of the main fault with respect to the free surface. The figure shows that the traveltimes of reflected S waves emitted from the nucleation patch (Fig. 10b) are well correlated with the spatiotemporal variation of slip velocity (Fig. 10a). The arrival of the reflected S waves changed the stress field and promoted slips behind the crack tip. It should also be kept in mind that the slip velocity on the faults is also affected by the reflected waves emitted from the propagating crack.

Implications for the Nankai subduction zone
Large interplate earthquakes have occurred repeatedly along the Nankai subduction zone off the coast of southwest Japan (Ando 1975), accompanied by disastrous tsunamis, most recently in 1944 (the Tonankai earthquake) and 1946 (the Nankai earthquake). Tsunami data suggested possible coseismic slips on splay faults during both earthquakes (Cummins & Kaneda 2000;Baba et al. 2006), although the resolution of seismic and tsunami waveforms was not fine enough to distinguish whether slip was dominant on the splay fault or on the plate interface. Seismic reflection data along the Nankai subduction zone suggested that there is a megasplay fault that runs continuously from the plate interface to the surface and cuts the old thrust faults within the frontal accretionary prism (Moore et al. 2007). From this, Moore et al. (2007) speculated that the megasplay fault was active during the 1944 Tonankai earthquake. Recently, the Integrated Ocean Drilling Program (IODP) Nankai Trough Seismogenic Zone Experiment (NanTroSEIZE; Tobin & Kinoshita 2006) drilled and cored several holes in the shallow portions of the megasplay and frontal plate boundary faults. Vitrinite reflectance geothermometry revealed that both fault zones experienced localized temperatures of more than 380 • C, implying that coseismic slip would have propagated at least once to the updip end of the megasplay fault and to the toe of the accretionary wedge (Sakaguchi et al. 2011). Tanikawa et al. (2012) measured the permeability of the core materials from the megathrust and frontal thrust of the Nankai subduction zone before and after sliding friction tests. They showed that the permeability of the splay fault materials under wet conditions is higher than that of the frontal thrust materials before the friction test, but the permeability of both fault materials decreases to around 10 −19 m 2 at a confining pressure of 3.5 MPa. These measurements imply that TP could be more effective on the frontal thrust than on the splay fault when the splay fault has not experienced a large slip, and that TP could be in effect on both frontal and thrust faults when the splay fault has experienced slips. Combining these implications with our results for cases with a small main fault dip angle (ψ = 13 • ; Fig. 7), we can envision the following scenarios: the rupture may propagate along the frontal thrust in the former case, but may propagate along the splay fault in the latter case. These implications are consistent with previous studies that suggested that the both splay and frontal thrusts would have slipped coseismically at least once (Sakaguchi et al. 2011) and that seismic slip could have occurred on the splay fault during the most recent event (Moore et al. 2007).

Validity of our results
The element size x 1 and the time step t in our simulations (Table 3) are sufficiently small so as not to change our results. The effects of x 1 and t are described in detail in Appendix B.
We assumed that fluid migration between the main and branch faults does not affect each other. This assumption is reasonable in our simulations. Since the hydraulic diffusivity ω is ω ≈ 10 −2 m 2 s −1 (Section 2.5), the pore pressure diffusion length 2 √ ωt is 0.2 m for t = 1 s. This is quite small compared to the element dimension, so the effect of fluid migration can indeed be ignored.
To investigate whether or not melting would occur during a rupture, we calculated the temperature increase due to frictional heating. We examined Model 1 in full-space without a free surface. We assumed that the temperature change at location x on the fault surface obeys the following equation (Bizzarri & Cocco 2006a): in the same way as pore pressure change (eq. 2). In both cases, where TP is in effect on only the main fault and where it is in effect on both faults, the temperature increase does not exceed 900 • C during the calculated time period (Figs 9a and b), except in and near the nucleation patch. If the melting temperature is approximately 1000 • C and the initial temperature is lower than 100 • C, which corresponds to the temperature 3 to 4 km deep, melting will not occur except in and near the nucleation patch. The temperature increase is high in the nucleation patch because we assumed that TP is not in effect there so that the nucleation processes are the same in all cases. TP lowers the temperature increase, preventing melting, as suggested by Sibson (1973).
In half-space cases, the rupture processes may differ depending on the depth of the fault intersection point because the time until the reflected waves trigger the rupture on the branch fault varies with depth. The results shown in Section 3.3 might not change qualitatively despite this variability.
In this study, we assumed uniform external stresses, while Kitajima & Saffer (2012) suggested that effective stress increases with distance from the trench axis by combining P-wave velocities obtained from geophysical surveys with empirical relationships between P-wave velocity, porosity, and effective mean stress defined by laboratory deformation tests on drill core samples of oceanic sediment. We also ignored the hydraulic fracturing that might be caused by increased pore pressure due to TP that was suggested by Ujiie et al. (2007). We assumed constant porosity, but dilatancy, which is an increase in pore space caused by fault slips, can affect dynamic ruptures (e.g. Bizzarri 2012 and references cited therein). These effects could be included for more realistic simulations in future work. We assumed the linear slip-weakening law (eq. 1). The slip-weakening law represents physically reasonable features in dynamic ruptures, along with other fault governing models (reviewed by Bizzarri 2011). We focused on 2-D mode II (in-plane) ruptures in this study. Our results would be still reasonable if the fault width along the trough axis is small in 3-D geometry. If the fault width is long, the mode III (anti-plane) component should be taken into account.

C O N C L U S I O N S
We performed a series of numerical simulations of dynamic fault branching including TP, based on the branch fault model of Kame et al. (2003). We revealed that TP can alter rupture propagation paths when the dip angle of the main fault is small. The rupture propagation paths depend on the branching angle when TP is not in effect on either of the faults. In contrast, the rupture propagates along the main fault in cases where TP is in effect on the main fault, but along the branch when TP is in effect on both faults. These features occur regardless of the branching angle. Thus, dynamic rupture processes are controlled more strongly by TP than by the branching angle. Finally, we considered the case when there is a free surface above the branch fault system. In this case, the rupture can propagate along both faults because of the combination of the free surface and TP on the main fault.

A C K N O W L E D G E M E N T S
This work was supported by the JST J-Rapid/ANR FLASH-Japan DYNTOHOKU project (ANR-11-JAPN-0009), the NIED project for Development of Monitoring and Forecasting Technology for Crustal Activity, and the NIED joint project with JAMSTEC on High Frequency Source Modelling. Two anonymous reviewers' comments were quite valuable in improving our manuscript. Professors Sánchez-Sesma and Iturrarán-Viveros kindly provided their computation code for the Garvin solution. Generic Mapping Tool (GMT version 4.5.6, Wessel & Smith 1998) software was used to draw some of the figures.

A P P E N D I X A : N U M E R I C A L M E T H O D F O R T H E F R E E S U R FA C E
The free surface was implemented in the code in the same way Hok & Fukuyama (2011) performed their 3-D BIEM. We used the kernels for full space derived by Tada & Madariaga (2001). In order to simulate a half-space, an interface composed of slipping crack elements was added at the location of the free surface, with a traction-free boundary condition as a constitutive relation. In practice, on the free surface elements, instead of solving eq. (1) for the fault elements, the following condition is solved: where τ is the shear stress and σ eff n is the effective normal stress. For each element, these equations are solved together with the BIEM kernels, in which the instantaneous slip velocity and the past time step loading can be separated. This leads to correct estimations of the surface motion as a response to the previous time step loading. On the fault elements, stress and slip velocity are estimated step by step. It should be noted that on the free surface elements, since the stresses are given, the free parameters are displacements in three directions.
To validate the 2-D implementation, we computed the Garvin's problem as described by Sánchez-Sesma & Iturrarán-Viveros (2006). This 2-D problem consists of placing an explosive line source 100 m below the surface to generate planar P and Rayleigh waves propagating at the surface. As we cannot generate a strictly explosive source, we use two perpendicular opening cracks; each is 0.5 m long. Fig. A1 shows synthetic seismograms for four surface receivers. The model parameters are the same as those used in fig.  (2006): a triangular slip source with duration of 0.25 s, v p of 2 km s -1 , and v s of 1 km s -1 . We set the density and the maximum of slip source time function at 1000 kg/m 3 and 0.5 m, respectively. The length of the surface elements and the time interval are 10 m and 0.0025 s, respectively, for the BIEM. As shown in Fig. A1(b), the waveforms computed by the BIEM and the analytic solutions by Sánchez-Sesma & Iturrarán-Viveros (2006) are quite similar. The wave amplitude for BIEM is slightly smaller than for the analytic solutions. This difference could be due to spatiotemporal discretization in BIEM. This comparison validates our computation with the free surface.

A P P E N D I X B : A C C U R A C Y O F O U R C O M P U TAT I O N S
The element size x 1 and the time step t in our simulations (Table 3) are sufficiently small so as not to change our results. To test the influence of x 1 and t on our results, we computed cases with TP using the same values of x 1 and t as those in the cases without TP (four times longer). We did not include the free surface in these calculations. We confirmed that the rupture propagation paths for the fine grid (Fig. 7) and coarse grid cases were all quite similar. The rupture velocity, however, was slightly larger with smaller x 1 and t values (Fig. B1) because the slip-weakening curves and stress concentration near the rupture front were reproduced more precisely in the fine grid case. We decided to show the results of our computations including TP with smaller x 1 and t in this paper because the slip velocity showed significantly less oscillations (Fig. B1b) than when using larger x 1 and t values. Actually, the values of x 1 and t in our simulations are smaller than those used in previous studies dealing with 3-D spontaneous ruptures (Bizzarri & Cocco 2006a;Urata et al. 2008).