Analysis of the role of branching angle in the dynamic rupture process on a 3‐D branching fault system

The fault branching phenomenon, which may heavily influence the patterns of rupture propagation in fault systems, is one of the geometric complexities of fault systems that is widely observed in nature. In this study, we investigate the effect of the branching angle on the rupture inclination and the interaction between branch planes in two‐fork branching fault systems by numerical simulation and theoretical analysis based on Mohr’s circle. A friction law dependent on normal stress is used, and special attention is paid to studying how ruptures on the upper and lower branch planes affect the stress and rupture on each other separately. The results show that the two branch planes affect each other in different patterns and that the intensity of the effect changes with the branching angle. The rupture of the lower branch plane has a negative effect on the rupture of the upper branch plane in the case of a small branching angle but has almost no negative effect in the case of a large branching angle. The rupture of the upper branch plane, however, suppresses the rupture of the lower branch plane regardless of whether the branching angle is large or small.


Introduction
Branching faults are widely observed in nature, and many natural earthquakes involve them, such as the 2008 M s 8.0 Wenchuan earthquake (Hao KX et al., 2009), the 1979 M w 6.4 Imperial Valley earthquake (Niazi, 1986), the 1992 M w 7.3 Landers earthquake (Sieh et al., 1993) and the 2002 M s 7.9 Denali earthquake (Eberhart-Phillips et al., 2003). Because fault branching is one of the complexities responsible for the complex patterns of earthquake source ruptures, it is important to study the dynamic rupture behavior of branching faults for earthquake hazard assessment and disaster prevention.
Since the 1990s, seismologists have made it technically possible to consider fault branching in the research on earthquake source dynamics. Because the boundary integral equation method (BIEM) includes the feature of reducing the dimensions of the considered problem by one, the BIEM is advantageous in dealing with complex fault geometry. Perhaps the simplest case is a static antiplane crack on branching faults, which was first considered by Tada and Yamashita (1997). Kame and Yamashita (1999) then sim-ulated the spontaneous bifurcation of a crack tip. However, neither study specifically focused on the geometry of branching faults. Poliakov et al. (2002) theoretically analyzed the stress field around a dynamic crack tip and predicted how the prestress state would affect the most favored direction for dynamic branching of a rupture. Kame et al. (2003) conducted simulations of dynamic rupture processes on branching faults and studied the effects of the prestress state, branching angle, and rupture velocity on the propensity to rupture. Bhat et al. (2007) studied the effect of branch length and predicted the interaction between a rupture on the main fault and on the branch plane. All the above studies used 2-D fault models. Aochi et al. (2000a) constructed a nonsingular BIEM in a 3-D elastic medium, which was used to study the effect of branching angles on the interactions between planes of branching faults (Aochi et al., 2000b) and the effect of normal stress on the rupture selectivity on branching faults (Aochi et al., 2002). Recently, Ando and Kaneko (2018) studied the effect of complex 3-D fault geometry on rupture propagation and termination during the 2016 M w 7.9 Kaikoura earthquake, which occurred on branching faults.
Other simulation methods, such as the finite element method (FEM), have also been used to study the rupture dynamics of branching faults. Oglesby et al. (2003) simulated the 1999 M w 7.1 Hector Mine earthquake on a complex fault system with a branch-ing fault to the north. DeDontney et al. (2012) investigated the earthquake rupture behavior on branching faults in an elastic-plastic medium. Xu DD et al. (2016) computed the dynamic ruptures along branching faults in the Longmen Mountain fault and analyzed the effects of the prestress state and fault geometry on rupture path selection. A limitation of the FEM is that different artificial treatments on the branching junctions in branching faults may result in different simulation outcomes (DeDontney et al., 2012).
Our study aimed to determine the effect of branching angles on the rupture propagation on branching faults and provide a physical explanation. This topic has been studied previously. Kame et al. (2003) pointed out the interactions between two branch planes, but they did not clarify the mechanism of the interactions. Aochi et al. (2000b) and Dong S and Zhang HM (2019) studied the stress interaction between the two branch planes, but the friction law they used was independent of normal stress, and they did not separate the stress effect generated by each part of the fault system. In this research, we simulated the dynamic propagation of a rupture on a 3-D two-fork branching fault model more realistically by using the BIEM with the slip-weakening Coulomb friction law (Kame et al., 2003), which is dependent on normal stress. First, the inclination of different angles of branch planes to rupture was analyzed explicitly by using a Mohr-Coulomb diagram. The stress accumulations produced by the slip of one of the branch planes on the other were then separated to analyze how the interactions between the two branch planes were affected by the branching angle.

N x
We used the BIEM to simulate dynamic rupture propagation. By considering an unstructured mesh of a branching fault system containing triangular elements, the discretized boundary integral equation can be written as (Qian F et al., 2019) is the stress of element at time p, is the initial stress of element m, is the stress kernel based on the triangular elements (Feng X and Zhang HM, 2017), which represents the stress increment of element at time generated by the unit slip of element at time , and is the slip rate of element at time . Equation (1) implies that the stress change of an element is influenced by the previous slip history of all elements, and the contribution of the slip of each element can clearly be calculated. Therefore, by using the BIEM, we can explicitly analyze the interaction between any two segments of a fault system.
As is usually adopted in the simulation of spontaneous dynamic ruptures, the slip-weakening friction law (Ida, 1972) was introduced in this study. For any given element on the fault system, it begins to rupture as soon as the shear stress on which exceeds the peak strength (yield stress) , and the shear stress decreases linearly with the increase in its slip distance until it reaches the residual stress level , as shown in Figure 1. Variable is the slipweakening distance related to the rock property and is set at 0.57 m in this study. We introduced the Coulomb criteria into the slip-τ p τ r μ s μ d weakening friction law (Kame et al., 2003), which defines and of an element proportional to its normal stress, and the scale factors were the static and dynamic friction coefficients, and , respectively. By combining the slip-weakening friction law and the BIEM, we could determine the slip rate of each element at each time step, based on which the corresponding slip and stress could also be obtained.

Δt
In this research, we considered a branching fault system in a 3-D unbounded homogeneous isotropic elastic medium, as shown in Figure 2a. Previous studies have shown that when a fault system is located at a depth of more than 1 km below the surface, the influence from the surface is very small and can be neglected (Zhang HM and Chen XF, 2006). The entire branching fault system consisted of three fault planes: a 30km × 10km primary plane a and two 15km × 10km branch planes b and c. We used the Cartesian coordinate system and took point A as the origin of the coordinate, with the x and y axes being parallel to OA and OO', respectively. The angle between planes a and b and that between planes a and c are denoted as and , both of which are counted counterclockwise. The fault model was discretized by unstructured triangular elements. Elements I, II, and III lie on the midline of the fault system and are next to the intersection, as shown in Figure 2b and 2c. The minimum of the inscribed circle radiuses of all the triangular elements was 249 m, and the time step of our simulation was 0.044 s. The density of the medium was 3,000 kg/m 3 , and the P-wave and S-wave velocities in the medium were 5.6 km/s and 3.2 km/s, respectively. The entire fault system was loaded with a uniform initial stress field. It was assumed that the shear and normal stresses on planes a, b, and c had no y component, and neither did the slip on each plane. As is shown in Figure 2a, a round nucleation zone with a radius of 2.5 km was located at the center of plane a, in which the shear stress was slightly greater than . The rupture began from the nucleation zone and propagated outward spontaneously. Given the initial stress components , and , the normal and shear stresses on each plane could be determined, and the spatiotemporal history of the rupture process in order of time could then be calculated by using the BIEM and the friction law. Throughout this study, MPa, and MPa. The initial shear stress in the nucleation zone was 1.2 to generate an Figure 1. Slip-weakening Coulomb friction law (Kame et al., 2003), where is the normal stress, is the slip-weakening distance, and and are static and dynamic friction coefficients, respectively. The peak strength ( ), the residual stress ( ), and the shear stress ( at any amount of slip ( is proportional to the normal compressive stress ( ).
initial crack. Because we were primarily focusing on the effect of , on the rupture patterns of planes b and c and the interaction between them in this study, all the values of geometric and physical quantities except and were fixed in our simulations.

Numerical Results
In our simulations, was always greater than , which always placed plane b above plane c. Therefore, the included angle of planes b and c, , were always positive. The rupture processes on two series of faults were computed: series 1, in which remained constant, and α 2 = −10°, −20°, −25°, −30°, −40°, and series 2, in which remained constant, and . Analysis of the numerical results will be helpful in understanding the effect of branching angles on how a rupture on the upper or lower branch planes separately influences rupture propagation on the other branch plane.  (Figure 3a), and the suppression becomes stronger as the absolute value of becomes smaller, i.e., becomes smaller. As shown in Figure 5a, the rupture propagation velocities on plane b in different cases can be compared in the space-time diagram of the rupture front. The corresponding result on plane b in the folding fault model is also shown for reference. Note that the rupture propagation velocity on plane b when is almost the same as that of the folding fault. In other cases, this velocity becomes smaller with a smaller .
In Table 1, we list the times when elements I, II, and III begin to crack in different cases of series 1, which indicate the times when the rupture reaches the intersection and begins on planes b and c, respectively. In cases and plane c cracks before the rupture front on plane a reaches the intersection, whereas plane b cracks after the rupture on plane a reaches the intersection. In cases and , both planes b and c crack after the rupture on plane a reaches the intersection. Generally, the larger the branching angle is, the sooner plane b cracks, and the later plane c cracks. Figure 4a shows that the rupture can propagate to and throughout plane c without an effect on plane b. The rupture on plane c dominates in cases with and , whereas the rupture on plane b dominates in cases with and . When the rupture propagates almost equally on planes b and c. The slip rate on plane c decreases slightly from to and violently when and . In Figure 5b, the rupture front propagates at almost the same velocity on plane c in cases with and and at a slightly slower velocity when . In contrast, in the case of , the crack on plane c obviously begins late and stops quickly.  Table 2 shows the times when the rupture reaches the intersection and begins on planes b and c. In cases with and , plane c cracks before the rupture on plane a reaches the intersection, after which plane b cracks. When , planes c and b both crack before the rupture on plane a arrives at the intersection, and plane c cracks sooner than plane b. Different from this scenario, when , plane b cracks first, and then the rupture on plane a arrives at the intersection and plane c cracks at almost the same time.

Effect of Branching Angles on the Chronological Order of Cracks
As shown in Tables 1 and 2, the chronological order at which elements I, II, and III begin to crack varies with the values of and , which can be illustrated by the Mohr-Coulomb diagram shown in Figure 6. Given the initial stress field, we can draw a Mohr's circle, on which any plane with a certain direction in this stress field can be donated as a point on the circle. We introduce the slip-weakening Coulomb friction law into the Mohr's circle to evaluate the inclination of cracking for different planes. A point close to the peak strength line marked by in Figure 6 on the Mohr's circle corresponds to a plane that cracks easily. A plane is not permitted to crack if its initial shear stress is less than the residual stress level marked by in Figure 6 (Aochi et al., 2000a). The direction of the plane that cracks the easiest can be defined as the optimal plane, as shown below: Φ opt −15.5 • where is measured from the maximum principal stress direction. We can then estimate the preference for rupture propagation by measuring the angle between a particular plane and the optimal plane. In this study, the angle measured from plane a to the optimal plane was (counted counterclockwise). In series 1, because plane a is closer to the optimal plane than is plane b, the crack cannot begin on plane b unless the rupture has reached the intersection. Plane c, with and is closer to the optimal plane than is plane a. Thus, element III cracks before element I does. In series 2, element III cracks before element I does when because plane c is closer to the optimal plane than is plane a.
The Mohr-Coulomb diagram can explain only the preference for rupture propagation before either of the two branch planes cracks because a branch plane will produce a stress effect on the other once it cracks. Take and as an example. Plane c is supposed to crack before the rupture arrives at the intersection because it is closer to the optimal plane than is plane a. In fact, however, it cracks after the rupture reaches the intersection because plane b cracks first and thereafter affects the rupture on plane c. Because element II is the first element on plane b to crack, the stress change on plane b indicates its potential for rupture growth. In series 1, we investigated the normal and shear stress change on element II produced by the rupture on planes a and c, denoted as and separately, by changing the sum in Equation (1) as follows where elements numbered from to are all the elements on plane a and those numbered from to are all the elements on plane c. The slip rate on plane a is almost independent of planes b and c, so according to Equation (3), its effect on the stress field of plane b barely changes with . We show this effect alone in Figure 7a to indicate that the rupture on plane a promotes the cracking of plane b. The effects of plane c alone on element II with different are shown in Figure 7b-7f.

Effect of the Lower Branch on the Upper One
In the case of (Figure 7b), the shear stress decreases rapidly while the normal stress decreases slowly on element II, quickly moving the stress on element II away from the peak strength line on the Mohr-Coulomb diagram. The stress variation on element II produced by plane c is almost opposite that produced by plane a. Thus, if plane c cracks prior to plane b, it rapidly produces a negative effect on the rupture of plane b, which is in accordance with the fact that plane b cracks both very late and very weakly when . When and (Figure 7c and 7d), the normal stress decreases more quickly and the shear stress decreases more slowly on element II. The stress on element II still moves away from the peak strength line, but at a slower speed. When (Figure 7e), plane c cracks before plane b, but the effect of the rupture on plane c slowly moves the stress on element II away from the peak strength line so that the rupture can begin on plane b shortly after plane c cracks and propagates throughout plane b. Finally, when , the rup-    Figure 6. Mohr-Coulomb diagram for planes with different dips in (a) series 1 and (b) series 2. is the angle between plane a and a branch plane b or c. The black filled hexagrams represent plane a, and the black filled pentagrams represent the optimal plane. Lines and are the peak strength and residual stress line, with the static and dynamic friction coefficients and , respectively. ture of plane c is weak, and the trend of its effect on element II (Figure 7f) indicates that the rupture on plane b is scarcely impeded.
Generally, the rupture on plane c (i.e., the lower branch) reduces both the normal stress and the shear stress on element II. As increases, the shear stress decreases more slowly and the normal stress decreases more quickly, which lessens the suppression on the rupture of plane b.

Effect of the Upper Branch on the Lower One
In series 2, we calculated the stress variation on element III generated by planes a and b, as shown in Figure 8. Because the rupture on plane a is barely affected by the ruptures on planes b and c, its effect on the stress on element III changes little in different cases in series 2. Figure 8a shows that the rupture on plane a causes the stress on element III to approach the peak strength line quickly, which promotes its cracking.
When and (Figure 8b-8d), the rupture on plane b α 1 moves the stress on element III away from the peak strength line, but this effect is very weak because of the suppressed rupture on plane b. To investigate specifically how the rupture on plane b affects the stress on element III without the influence of plane a resulting from different in these three cases, we assumed that the direction of plane a equaled that of plane b and that plane c did not produce any stress effect on plane b so that the rupture on plane a would propagate to plane b unimpeded. We then calculated the effect of plane b on the stress on element c, as shown in Figure 9. The only difference between Figure 9 and Figure 8b-8d is that the ruptures on plane b in Figure 9 are stronger than those in Figure 8b-8f.

Δα Δα
Figure 8e-8f and Figure 9 show that the rupture on plane b increases the normal stress and decreases the shear stress on element III. When the value of is large, the normal stress increases more rapidly and the shear stress decreases more slowly, whereas the normal stress increases more slowly and the shear stress decreases more rapidly when becomes smaller. Within the range of investigated in this study, the stress on element III moves away from the peak strength line no matter whether the value of is large or small, and the medium value of causes the stress on element III away from the peak strength line to be the fastest. When (Figure 8f), plane b cracks very early and begins to produce a negative effect on the rupture of plane c, making it difficult for plane c to crack, which is in accordance with the fact that plane c hardly cracks when .

Conclusions
In this paper, we applied the BIEM to 3-D branching fault systems to simulate the dynamic rupture process on different fault models with different branch angles in two series. Under a given uniform initial stress state, the direction of a plane determines the normal and shear stress on that plane, which affects its tendency to crack. We calculated the direction of the optimal plane and used the Mohr-Coulomb diagram to illustrate that a plane cracks easily if it is close to the optimal plane on the Mohr-Coulomb diagram and is difficult to crack if it is far from the optimal plane.
In a branching fault system with two branch planes, the stress change produced as one of the branch planes ruptures will promote or suppress the rupture on the other branch plane. We investigated this effect by separating the stress change on a branch fault produced by the other branch plane from the total stress change, as summarized in Table 3. After the lower branch plane cracks, it causes both the normal and shear stress on the upper branch plane to decrease. When the angle between the two branch planes is small, it decreases the shear stress rapidly and decreases the normal stress slowly on the upper branch plane, which suppresses the rupture on the plane. In contrast, when the angle becomes larger, the shear stress decreases more slowly and the normal stress decreases more rapidly; thus, the suppression of the rupture on the upper branch plane is weakened. The rupture on the upper branch plane causes the normal stress to increase and the shear stress to decrease on the lower branch plane. The normal stress increases more slowly and the shear stress decreases more rapidly when the angle between the two branch planes becomes larger. Regardless of whether the angle between the two branch planes is large or small, the rupture on the lower branch plane is suppressed by the rupture on the upper branch plane.
. The ordinate and abscissa of a dot represent the shear stress and normal stress on element III, respectively, and the color represents the time of the dot. The rupture on plane a is assumed to propagate to plane b unimpeded.