Rupture Dynamics of Cascading Earthquakes in a Multiscale Fracture Network

Fault‐damage zones comprise multiscale fracture networks that may slip dynamically and interact with the main fault during earthquake rupture. Using 3D dynamic rupture simulations and scale‐dependent fracture energy, we examine dynamic interactions of more than 800 intersecting multiscale fractures surrounding a listric fault, emulating a major listric fault and its damage zone. We investigate 10 distinct orientations of maximum horizontal stress, probing the conditions necessary for sustained slip within the fracture network or activating the main fault. Additionally, we assess the feasibility of nucleating dynamic rupture earthquake cascades from a distant fracture and investigate the sensitivity of fracture network cascading rupture to the effective normal stress level. We model either pure cascades or main fault rupture with limited off‐fault slip. We find that cascading ruptures within the fracture network are dynamically feasible under certain conditions, including: (a) the fracture energy scales with fracture and fault size, (b) favorable relative pre‐stress of fractures within the ambient stress field, and (c) close proximity of fractures. We find that cascading rupture within the fracture network discourages rupture on the main fault. Our simulations suggest that fractures with favorable relative pre‐stress, embedded within a fault damage zone, may lead to cascading earthquake rupture that shadows main fault slip. We find that such off‐fault events may reach moment magnitudes up to Mw ≈ 5.5, comparable to magnitudes that can be otherwise hosted by the main fault. Our findings offer insights into physical processes governing cascading earthquake dynamic rupture within multiscale fracture networks.

forming a fault core surrounded by a damage zone that includes subsidiary faults, highly fractured material, and distributed macro-fractures (Faulkner et al., 2010). Chester and Logan (1986) proposed a fault-zone model based on field observations of the Punchbowl fault, Southern California, with a damage zone generated by a combination of co-seismic events (Chester et al., 1993;Mitchell & Faulkner, 2009 and aseismic/quasi-static processes (Faulkner et al., 2011;Griffith et al., 2012). The damage zone comprises brittle, heterogeneous, anisotropic, discontinuous materials, and multiscale fractures oriented in different strike and dip directions (Schulz & Evans, 2000;Faulkner et al., 2010). These fracture networks may have formed under different tectonic stress conditions in the past. The ensemble of multiscale fractures of different ages and orientations generates strong spatial variations in mechanical properties (Ostermeijer et al., 2020).
The properties of the fault-damage zone play an essential role in fault mechanics and earthquake rupture dynamics, which in turn generates co-seismic off-fault damage and geometrical complexity (Andrews, 2005;Dunham et al., 2011;Gabriel et al., 2013;Okubo et al., 2019;Cappa et al., 2014). However, dynamic source characteristics of earthquake rupture within a fracture network are largely unexplored.
Several fault segments may rupture in a single large earthquake, like during the 1992 M w 7.3 Landers, the 2012 M w 8.6 off-shore Sumatra, the 2016 M w 7.8 Kaikōura, or the 2019 M w 6.4, and M w 7.1 Ridgecrest earthquakes. For large magnitude multi-segment earthquakes (M w > 7), dynamic rupture branching or "jumping" across fault segments has been inferred from observations or in simulations (e.g., Hauksson et al. (1993); Harris and Day (1993); Oglesby (2008); Meng et al. (2012); Wollherr et al. (2019); Ross et al. (2019); Taufiqurrahman et al. (2023)). The characteristics of such cascading multi-segment rupture include relatively closely spaced segments (Harris & Day, 1993), some fault segments being optimally oriented and/or critically stressed within the ambient stress field (e.g., Ulrich et al. (2019)). For example, rupture nucleation on a favorably oriented fault segment may create a cascading earthquake, while nucleation on a less favorably oriented segment may lead to early rupture arrest (e.g., Oglesby and Mai (2012) ;Kyriakopoulos et al. (2019); Lozos and Harris (2020)). In this process, co-seismic stress transfer plays a pivotal role, altering the imminent stress conditions of adjacent fault segments by decreasing or increasing local stresses.
Correspondingly, further rupture propagation may be promoted, or impeded by pre-stress fault conditions, depending on fault orientation within the ambient stress field.
In this study, we explore smaller-scale cascading earthquake ruptures, for example, in fracture networks of a fault damage zone or forming a geo-reservoir. Generally, a fracture network comprises more than one fracture population; each population of fractures is then characterized by its size distribution and its dominant strike direction, which depends on the stress orientations and tectonic conditions at the time of fracture formation (Zoback & Kohli, 2019). The prescribed fracture size, density, strike, and dip orientation can be described statistically (Dershowitz et al., 2019) and constrained from observations (Mitchell & Faulkner, 2009;Savage & Brodsky, 2011).
In particular, we investigate the potential for cascading earthquakes on multiscale fracture networks with a very large number of fractures. Instead of considering only several fault segments (as in the above- Figure 1. Schematic view of multiscale fault structure and its associated damage zone that accommodates continued deformation in strike-slip faulting (modified after Tchalenko (1970)).
mentioned studies for M w > 7), we include in our simulations over 800 intersected fractures of different sizes surrounding a listric main fault. We examine how dynamic earthquake rupture cascades may occur in such a fracture network if a subset of fractures is favorably oriented. We also explore mechanisms by which a rupture cascade starting on a fracture may promote or suppress rupture on the main fault. Pre-stress conditions may or may not lead to a sustained dynamic rupture cascade within the fracture network, or to the activation of the main fault. Therefore, we investigate dynamic rupture processes within a complex fracture network under varying pre-stress conditions, rupture nucleation points, and levels of fluid overpressure. To facilitate comparisons with observations, we then analyze earthquake source parameters (i.e., magnitude, centroid moment tensor, stress drop, rupture speed) and general aspects of seismic wave radiation.

Model Setup
We model dynamic earthquake rupture scenarios with a principal listric fault and many off-fault fractures. Our model represents a typical fault structure in a sedimentary basin near a subsurface reservoir (Figure 2a; e.g., Gibbs (1984); Hardman and Booth (1991); Ward et al. (2016)), as we are also interested in examining general aspects of induced earthquakes in geo-reservoirs. In deep sedimentary basins, listric faults often extend to basement depth and form traps for natural resources near geological reservoirs (e.g., Withjack et al. (2002); Onajite (2013); Dixon et al. (2019)).
Dynamic rupture simulations are governed by initial stress conditions, including fault strength, fault and fracture geometry, subsurface material properties, the nucleation procedure, and the friction law. In this study, we vary the direction of maximum horizontal stress and constrain other parameters from observations. In the following, we explain how we construct the fracture network, considering both field observations and statistical estimates. We summarize the applied material and fault properties and how the initial stress and fault strength conditions are specified.

Fracture Network Construction
Numerous studies have characterized the fault damage zone focusing on its width and fracture distribution, as these are critical for fluid transport, naturally in the context of fault-valving effects and fault healing, but also in industrial applications like oil & gas exploration, geothermal production, or waste-water injection (Rice, 1992;Faulkner et al., 2010Faulkner et al., , 2011. Numerical simulations and field observations suggest that the formation of off-fault fractures is affected by factors like fault geometry and co-seismic displacement, tectonic environment and ambient stress state (Faulkner et al., 2011;Okubo et al., 2019;Wu et al., 2019;Gabriel et al., 2021;Sainoki et al., 2021). A key property of off-fault fractures is that their density decreases with increasing distance from the fault core, typically following a power-law scaling that depends on long-term stress evolution, rock type, and fault maturity (Sainoki et al., 2021).
Using a statistical approach, we generate two approximately conjugate fracture families to achieve a high degree of fracture connectivity (henceforth referred to as fracture family 1 and family 2). The conceptual geometry is illustrated in Figure 2a. The fracture network is generated using the commercial software FRACMAN (Dershowitz et al., 2019), whereby we consider four key quantitative properties: fracture density, fracture size distribution, fracture orientation distribution, and fracture shape. We discuss these properties and our assumptions in the following section, but state already here that the resulting fracture network comprises in total 854 fracture surfaces, with 423 and 431 fracture surfaces in family 1 and 2, respectively. 78% of all fractures are connected to more than one other fracture, 3% fractures are connected to only one other fracture, and 19% are unconnected.

Fracture Density
Field observations suggest a power law decay of fracture density with the distance from the fault core (Mitchell & Faulkner, 2009;Savage & Brodsky, 2011). They use the number of fractures per unit length (P 10 ) to quantify fracture density. However, P 10 is inefficient for describing a volumetric fracture density distribution because it only specifies fracture distribution along one line. To overcome this limitation, we adopt fracture density based on P 32 , defined as the ratio between the area of all fractures and the rock volume. However, observing P 32 in nature is difficult. We thus constrain P 32 using an inferred P 10 value based on a multi-dimensional intensity metric (Wang, 2005). The particular relation between P 32 and P 10 is implemented in FRACMAN, which we apply to generate the fracture network. To validate the P 32constrained fracture network, we measure P 10 of the modeled fracture network using an average number from linear transects. The such defined average fracture density P 10 in our model is well approximated by -5-manuscript submitted to JGR: Solid Earth (1) The exponential decay follows m = 0.8 and the fault-specific constant is C = 2.5, whereby d is the fracture distance (Savage & Brodsky, 2011). The nearest neighbor method is used to populate the volumetric fault zone with fractures. As the fracture-network generation is statistically performed, small fluctuations in the final fracture density are expected, hence the lower slope at distances d < 200 m (Figure 2b), indicating that the FRACMAN-generated number of smaller fractures is slightly lower than given by Eq.1.

Fracture Size Distribution and Main Fault Geometry
Fracture lengths within a fault damage zone range from micrometer-to kilometer-scales for mature faults (Tchalenko, 1970). The maximum fracture length is limited by the width of the fault damage zone (Mitchell & Faulkner, 2009). Various factors influence fracture size, including fault length, fault roughness, in-situ stress conditions, maximum fault displacement, and distance from the fault core (Perrin et al., 2016;Sainoki et al., 2021). In studies of natural fracture networks, a power-law relation is commonly used to characterize the fracture-size distribution within assumed maximum and minimum fracture length (Panza et al., 2018). Our chosen fracture size distribution follows a power-law with fracture lengths between 500 m and 100 m (Figures 2c and S1). A similar distribution is described by Lavoine et al. (2019) who develop a theoretical approximation based on natural fracture networks in geological reservoirs considering a power law exponential decay of 2 which we use as well ( Figure S1). The maximum fracture size is determined by our assumed damage zone thickness of 1 km, the shortest fracture length of 100 m is due to computational considerations.
Our main fault extends 8 km along-strike and 4 km along-dip. For reference, we define the strike direction as N270E with respect to the Cartesian system in Figure (Thingbaijam et al. (2017)).

Fracture Orientation
We prescribe fracture orientations based on field observations and numerical models. For example, the observed orientation of large-scale fractures and conjugate faults in the damage zone of the Caleta Coloso strike-slip fault occur around 25 • − 30 • with respect to the main fault (Mitchell & Faulkner, 2009). 2D dynamic rupture simulations reveal that small-scale and meso-scale discrete off-fault fractures spontaneously group into two sets of conjugate fracture orientations separated by ∼ 100 • (Dalguer et al., 2003;Ando & Yamashita, 2007;Okubo et al., 2019;Gabriel et al., 2021).
Our model comprises two fracture populations. The first one, fracture family 1, has an average strike of N120E, while fracture family 2 has an average strike of N20E (Figure 2d). Both fracture families form an average angle of 100 • to each other and an angle of 30 • and 70 • with respect to the main fault's strike. For both families, we consider a 10 • standard deviation in the strike and dip angle. Distributions in both strike and dip of the fractures follow a Fisher distribution (Fisher, 1995).

Fracture Shape
Detailed knowledge of 3D fracture geometry in nature is scarce due to incomplete outcrop data and limited resolution in 3D reflection-seismic data. The surface of small fractures is commonly assumed to be a penny-shaped disk of infinitesimal thickness (Priest & Hudson, 1981;Laslett, 1982;Piggott, 1997;Berkowitz & Adler, 1998). The main control on fracture aspect ratio may be the mechanical anisotropy due to rock layering (Nicol et al., 1996). The aspect ratio of fractures is typically assumed to vary from 0.5 -3.5, irrespective of slip direction. For simplicity, we assume elliptical fractures with aspect ratio 2. We assume pure shear fractures (mode II and mode III slip) and do not allow for fracture opening (mode I).

Material Properties
We assume uniform material properties to model fractures embedded in a 3D elastic half-space. The main fault and fractures are buried and do not intersect the free surface. This setup allows us to focus on interactions of fractures and the main fault without rupture and wave-propagation complexities introduced by material interfaces. We choose a Poisson solid with v p = 6000 m/s, v s = 3464 m/s, ρ = 2670 km/m 3 , and λ = µ = 32 GPa.

Initial Stresses, Friction Law, and Fault Strength
We assume a transitional strike-slip to normal faulting stress regime. Local stress conditions for each fracture and the main fault are modulated by their geometry and orientation within the regional A transitional stress regime may promote different faulting styles during a single compounded rupture (e.g., 1992 Landers and 2016 Norcia, 2016 Kaikōura earthquakes (Hauksson, 1994;Wollherr et al., 2018;Tinti et al., 2021;Ulrich et al., 2019), seen also for induced earthquakes (Schoenball et al., 2014;Palgunadi et al., 2020)).

Initial Stresses
We use a Cartesian initial stress tensor to load the fault and all fractures. We constrain the stress regime by the relative stress magnitude A φ (Simpson (1997) given as where n indicates the faulting style (n = 0 for normal faulting, n = 1 for strike-slip faulting, and n = 2 for reverse faulting), and φ denotes the stress shape ratio given by Here, σ 1 , σ 2 , and σ 3 are the maximum, intermediate, and minimum principal stress magnitudes, respectively.
We consider a stress shape ratio of φ = 0.9 and a predominantly normal stress regime (n = 0, (Simpson, 1997)), and hence enforce A φ = 0.9. This transitional stress regime implies that two faulting styles, normal and strike-slip, may be activated due to S v ≥ SH max SH min , where S v = σ 1 , SH max = σ 2 , and SH min = σ 3 are vertical overburden stress, maximum horizontal stress, and minimum horizontal stress, respectively.
Our model includes a depth-dependent effective normal stress, with pore pressure P f (z) = γρgz, where g is the gravitational force, z is depth, and γ is the fluid pressure ratio (Ulrich et al., 2019), which value is given by ρ water /ρ = 0.37 when pore fluid pressure is hydrostatic. The pore pressure counteracts the fault normal stress in the effective normal stress σ n = σ n − P f . We vary γ in the range γ = [0.37; 0.80] to explore varying fluid-overpressure scenarios.

Friction law and fault strength
The fault strength (τ p ) is defined by the relation between friction law and effective normal stress: We adopt the laboratory-based rate-and-state friction law with rapid velocity weakening (Lapusta et al., 2000;Rice, 2006;Noda et al., 2009;Dunham et al., 2011): where V is slip rate, V 0 is the reference slip rate, and a is the direct effect. The state variable θ evolves according to: -9-manuscript submitted to JGR: Solid Earth where L is the state-evolution slip distance for rate and state friction law. θ SS at steady-state is: where f SS (V ) is the friction coefficient at steady-state: where the low-velocity steady-state friction f LV (V ) is defined as: with the evolutional effect given by b.
Friction parameters are chosen to generate realistic stress drops and frictional resistance. Aside from the state-evolution slip distance (L, see section 2.4), all frictional parameters are constant (Table 1). We assume that the main fault and all fractures are frictionally unstable (b − a > 0). We use a characteristic weakening velocity V w = 0.1 m/s based on the experimentally observed onset of rapid decay of the effective friction coefficient, ranging between 0.01 -1 m/s (e.g., Rice (2006); Beeler et al. (2008); Di Toro et al.
We choose a steady-state friction coefficient of f 0 = 0.6 at reference slip velocity value V 0 = 10 −6 m/s. The steady-state weakened friction coefficient f w is set to 0.1, similar to Rice (2006).
Earthquake rupture dynamics are largely controlled by the relative pre-stress ratio R that describes the ratio of the maximum possible stress drop and frictional strength drop (Aochi & Madariaga, 2003) In Eq. 10, τ 0 is the initial shear traction, τ d = f w σ n is dynamic stress, and ∆τ d is dynamic stress drop. τ p = f p σ n is the peak dynamic stress, and τ p − τ d is the maximum dynamic strength reduction.
The peak value of the friction coefficient (f p ) depends on the rupture dynamics (Garagash, 2021) but is approximated in evaluating Eq. 10 by the reference value f 0 . The value of f p in simulations varies along the fault and fractures and may exceed f 0 but rarely falls below it. R 0 is the maximum possible value of R for a fracture at the most-optimal orientation. Hence, the fault-local R is always smaller than or equal to R 0 . We prescribe the maximum pre-stress ratio, R 0 , to be constant across the model to constrain the initial stress state (Section 2.3.1). R 0 depends on the proximity to failure of an optimally oriented fault (R 0 = 1 corresponds to the maximum degree of stress criticality, i.e., an optimally oriented fault is at "failure" for the given initial stress state). We consider R 0 = 0.8 in all simulations. In nature, faults and fractures are -10-manuscript submitted to JGR: Solid Earth typically not optimally oriented. In a dynamic rupture scenario, only a small part of a fault or fracture must reach failure to initiate a sustained rupture.

Fault-Size-Dependent State Evolution Slip Distance (L) and Fracture Energy (G)
Dynamic rupture modeling across different fault or stress-heterogeneity scales may use variable characteristic slip distances (Bizzarri & Cocco, 2003;Ando & Yamashita, 2007;Galvez et al., 2021;Ulrich et al., 2022). The critical nucleation size required to initiate spontaneous or runaway rupture (i.e., ruptures across the entire main fault (Galis et al., 2019)), scale with characteristic slip distance L under rate-and-state friction (Rubin & Ampuero, 2005), or with critical slip distance D c under linear slip-weakening friction (Andrews, 1976;Galis et al., 2015). To model dynamic rupture on fractures of different sizes, we adopt a fracture-scale dependent L. Our scaling of L  is linked to the scale-dependence of the fracture energy emerging from seismological observations (e.g., Abercrombie and Rice (2005)) and earthquake physics.
The elastic energy, stored in the rock volume and available to be released and dissipated as frictional heat, fracture energy G, expended to propagate the rupture, and radiated in seismic waves, scales with the fault size R (e.g., Madariaga (1976)). Thus, fracture energy has to scale with fault size for the fault to be able to host dynamic rupture . Such linear scaling of G with fault size R has been previously proposed in relation to fault growth (Scholz et al., 1993) and co-seismic rupture with off-fault plasticity (Andrews, 2005). Recently, Gabriel et al. (2023) inferred from a multi-weakening fracture energy decomposition applied to seismologic estimates of fracture energy (Abercrombie & Rice, 2005;Tinti et al., 2005;Mai et al., 2006;Causse et al., 2014;Viesca & Garagash, 2015) that: Here, we assume G ≈ G c (R) and apply an equivalent scaling of the state evolution slip distance L following an analytical approximation of fracture energy for a rate-and-state governed fault (Garagash, 2021): We use the reference value f 0 to approximate the peak friction f p at the rupture front, and σ n ∼ 40 MPa representative of lithostatic and hydrostatic gradients at the median fault depth (∼ 3 km) to define:

Results
In the following, we present eighteen 3D dynamic rupture simulations, each on the same multiscale fracture network with an embedded listric main fault. All simulations are carried out using the opensource software SeisSol (Section 6 and Appendix A). Ten simulations are scenarios under varying pre-stress conditions. In five simulations, we shift the rupture nucleation location within the fault damage zone and vary fluid overpressure levels. Finally, in three scenarios, we use higher-resolution simulations for generating reliable high-frequency seismic waveforms (Appendix C). Figure 3 shows the computational high-resolution mesh overlain by a snapshot of absolute slip and vertical particle velocity at time t = 6 s for a pure rupture cascade (explained below). To improve visibility and aid interpretation, 3D views of our results are presented in an "exploded" view, generated by displacing each fracture centroid from coordinates (x, y, z) to (3.5x, 3.5y, z).
We first analyze fault and fracture strength for varying maximum horizontal stress (SH max = σ 2 ) orientations without performing dynamic rupture simulations (Section 3.1). In Section 3.2, we describe ten dynamic rupture simulations for variable SH max orientation (Ψ), and analyze cascading dynamic rupture within the fracture network and on the main fault. In Section 3.3, we describe the earthquake kinematics of all fracture network-main fault rupture scenarios. Next, we compare two scenarios of cascading ruptures initiated at different locations and with different mechanisms in Section 3.4: (i) a dynamic rupture cascade initiating on a fracture due to overstress at ∼1 km away from the main fault, and (ii) a dynamic rupture cascade initiating at the same location as (i) assuming an elevated (in excess of hydrostatic) pore fluid pressure regime within the fracture network. Finally, in Section 3.5, we describe seismic waveform characteristics.

Fault and Fracture Strength for Variable Maximum Horizontal Stress Orientation
We examine different orientations of maximum horizontal stress, ranging from Ψ = 40 • − 120 • with respect to North (normal to the main fault strike; Figure 2e). Changes in Ψ gradually alter the criticality of the main fault and all smaller fractures, ranging from favorably to unfavorably oriented fractures with respect to ambient stress conditions. Recall that favorably oriented fractures or fault segments are characterized by higher values of R, whereby R depends on how the normal and shear tractions are resolved on fractures with varying orientations within the ambient stress field.
The static slip tendency analyses to examine how variations in Ψ result in favorably or unfavorably oriented dynamic rupture planes (fractures, main fault, or both) do not yet involve dynamic rupture simulations, but still provide a valuable preliminary assessment of possibly mechanically viable conditions for earthquake rupture (Palgunadi et al., 2020). They consist of quantifying the relative pre-stress ratio R on the main fault and on the fractures for a given R 0 . However, the static slip tendency analyses cannot fully anticipate the dynamic processes or potential interactions between multiscale fractures, particularly if fractures are unfavorably oriented. As we show later, in several cases, dynamic stress interactions between multisegmented fractures facilitate sustained rupture even on unfavorable oriented fractures.
For the assumed A φ = 0.9 and using Eq. 10 with f 0 = 0.6 and f w = 0.1, we characterize and display the failure propensity of fractures and the main fault for different Ψ ( Figure 4). Based on this analysis, we identify three prominent cases for subsequent dynamic rupture modeling: Case 1 has unfavorably oriented fractures and main fault plane for Ψ = 40 • ; Case 2 has favorably oriented fractures and unfavorably oriented main fault for Ψ = 65 • ; Case 3 has unfavorably oriented fractures but favorably oriented main fault plane for Ψ = 120 • . Considering the variability in strike and dip orientation directions of the small-scale fractures, some fractures will remain unfavorably oriented even under generally favorable orientation (e.g., Ψ = 65 • ).

Cascading dynamic rupture in multiscale fracture networks
To investigate the conditions under which cascading dynamic rupture in the fracture network may occur, we conduct ten 3D dynamic rupture simulations including the three cases (Case 1-3) defined above.
The seven additional simulations are defined as "rupture on unfavorably to favorably oriented fractures but unfavorably oriented main fault" for Ψ = 50 • (Case 1a) and Ψ = 60 • (Case 1b) and "rupture on unfavorably to well oriented main fault plane" for Ψ = 70 • , 80 • , 90 • , 100 • , 110 • (Case 2a). We prescribe a common hypocenter at x = −800 m, y = 600 m, and z = −2800 m. For all cases in this section, we initiate rupture inside a 400 m-radius sphere centered on the main fault that includes a subset of fractures. This radius is larger than the critical nucleation sizes on fractures (see Appendix D).

Case 1, Ψ = 40 • , unfavorably oriented fractures and main fault
We refer to this case as "failed rupture nucleation". For Ψ = 40 • , the static analysis ( Figure 4) shows unfavorable orientation for the two fracture families and the main fault. We observe that no self-sustained rupture is generated ( Figure S2), as the slip rate remains confined within the nucleation volume (Movies S1a, S1b). Increasing the nucleation radius (we tested radii of up to 600 m) and increasing nucleation overstress (up to R 0 =5) do not lead to self-sustained rupture initiation. Interestingly, during the nucleation phase, rupture branching occurs at the intersection of two fractures (red arrow in Figure S2). However, rupture is arrested subsequently, and slip remains small with an average of less than 0.02 m.
-14-manuscript submitted to JGR: Solid Earth  Figure S2; Movies S3a and S3b). We identify multiple backward-propagating slip episodes within the fracture network at the footwall side toward the nucleation area. The total rupture duration is t ≈ 3 s.
Late-stage rupture on fractures occurs mainly in the fracture network surrounding the eastern part of the main fault (green to yellowish colors in Figure S2). The dynamic rupture cascade terminates spontaneously on different fractures via three mechanisms: (1) abrupt rupture cessation at individual fracture boundaries, despite connected neighboring fractures (2) smooth rupture stopping after reaching unfavorably oriented connected fractures (R = 0.1 − 0.25), and (3) rupture termination on isolated fractures due to relatively large spacing between fractures that precludes dynamic rupture jumping (here > 80 m).
Maximum slip is 0.4 m, located on the deepest fractures at a depth of 5 km ( Figure S4), and average slip across all slipped fractures is 0.04 m. This scenario exhibits sustained cascading rupture or pure cascade behavior, which may be expected from the static analysis ( Figure 4) due to the favorable orientation of most fractures. We term the cascading rupture characteristics as "sub-optimal", since only ∼ 65%, that is, 536 out of 854 fractures dynamically slipped ( Figure S5). Slip is distributed on fractures within both the hanging wall and footwall fault zone of the main fault plane. Interestingly, fractures on the hanging wall slip first, followed by those on the footwall. However, the dynamic rupture cascade does not activate fractures located in the western region of the main fault footwall.
-16-manuscript submitted to JGR: Solid Earth We highlight two fracture-to-fracture cascading mechanisms observed in the simulation with Ψ = 60 • : (i) dynamic rupture transfer via rupture branching across connected fractures, or (ii) rupture jumping across unconnected fractures. We observe that rupture branching dominates the cascading process, activating approximately ∼ 94% of the slipped fractures. The remaining ∼ 6% of slipped fractures are in close proximity to the evolving cascade (with an average distance of < 80 m between fractures) and are activated by rupture jumping.
3.2.4 Case 2, Ψ = 65 • , dynamic rupture on favorably oriented fractures and an unfavorably oriented main fault We refer to this case as a "pure cascade". Most fractures are favorably oriented for Ψ = 65 • and sub-critically stressed (τ /σ n in the range 0.22 − 0.5), but the main fault is unfavorably oriented. Likewise, six fractures have a low relative pre-stress ratio (R < 0.3; Figure 4), which are fewer than in the case of Ψ = 60 • but enough to eventually arrest a dynamic rupture cascade. In this simulation, dynamic rupture propagates as a sustained cascade, generating dynamic slip across the entire fracture network ( Figure 6).
The main fault experiences only small, localized slip despite several fault-fracture intersections (see rupture time evolution in Figure 6 and Figure S3).
Dynamic rupture propagates in a zigzaging pattern, activating neighboring fractures of both fracture families. The cascade progressively evolves from one fracture to another, favoring those that are connected or  Figure 6a. The inability of fractures to activate the main fault is likely due to the fracture energy of the main fault being ∼ 10 times larger than that of the largest fractures.
-18-manuscript submitted to JGR: Solid Earth Abrupt rupture termination creates locally more widespread fracture-network slip at the main fault boundary.
Rupture jumping also occurs between parallel fractures adjacent to the main fault. Slip on fractures unconnected to the main fault is small and quickly self-arrests (Figure 7b). The percentage of slipped fractures is ∼ 58%, i.e., 477 fractures ( Figures S5 and S7). In this case, slipping fractures are located preferably at the dilatational side of the main fault, that is, at the hanging wall in the eastern and the footwall in the western main fault (Figures 7 and S7).

Rupture kinematics
In the following, we analyze the kinematic rupture properties of all ten dynamic rupture simulations with varying Ψ, including equivalent point source moment tensor estimates, seismic moment (magnitude), moment rate function, average rupture speed on individual fractures or the main fault, and stress drop (∆τ ).
Inspecting the kinematic parameters for cascading earthquakes may help to identify observable signatures of cascading rupture episodes in real fault zones.

Equivalent point source moment tensor
We examine the apparent far-field source mechanism by calculating an equivalent moment tensor for a

Moment rate function
We analyze the moment rate functions (MRFs) and their amplitude spectra for different simulations (Figures 8a and 8b). Cases with large Ψ = 100 • − 120 • that include runaway rupture on the main fault show a simple triangle-shaped moment rate function (Meier et al., 2017), despite the activated off-fault fractures (Ψ = 100 • − 120 • at t < 0.5s in Figure 8b). These MRFs include small amplitude variations shortly after rupture initiation, which can be observed more clearly using a logarithmic scale (Figure 8b).  Figure 8c shows our modeled seismic moment spectra that follow a classical spectral decay with ω −2 at high frequencies (e.g., Aki (1967); Vallée et al. (2011)).

Average rupture speeds and cascading speed
Analyzing rupture speed in complex 3D models is challenging (Bizzarri, 2013). We define average rupture speeds in three different ways: (1) the overall average rupture speed for the entire process or seismologically-inferrable average rupture speed (v R,tot ), (2) the surface-averaged rupture speeds across all fractures and the main fault (v R,ave ), and (3) the "cascading speed" (v C ). These three different rupture speeds are useful to distinguish rupture speeds for cascading and non-cascading ruptures. v R,tot (1) is calculated as the ratio of the distance of slipped fractures and the main fault relative to the hypocenter and the duration of the rupture. The cascading speed (v C ) quantifies the apparent rupture speed along the fracture network and is determined from the slope of a scatter plot showing the distance to the hypocenter versus rupture onset time at each slipped fracture ( Figure 9). We manually determine the cascading speed for each Ψ-case. We ignore Case 1, Ψ = 40 • , where slip occurs only within the nucleation volume. Table 2 shows that v R,tot ranges from 0.63c s to 0.80c s , where c s is the shear wave speed. For Ψ < 90 • , cascading rupture across the fracture network results in lower v R,tot , with v R,ave > v R,tot . We observe supershear rupture (c s ≤ v R,f < v p ) in parts of the fracture network, for instance for Ψ = 65 • (see Figure S9a), resulting in higher v R,ave .
Supershear rupture occurs on smaller fractures after rupture branches or jumps onto them from larger fractures or from the main fault during runaway rupture ( Figure S9b,c,d). Several larger fractures also host supershear rupture speeds. Although inferences of supershear rupture speeds are rare for natural earthquakes, recent studies suggest that supershear rupture may frequently occur at a local scale, i.e., over a small part -25-manuscript submitted to JGR: Solid Earth  Compared to the surface-averaged rupture speeds on all fractures and the main fault v R,ave , the value of the cascading speed v C is consistently lower for cascades propagating within the fracture network, reflecting the delayed rupture propagation due to rupture branching and jumping. When considering Ψ = 65 • , which involves the greatest number of slipped fractures, v C is notably smaller than v R,ave , at v C = 0.65c s ( Table   2). For Ψ = 70 • and Ψ = 80 • , v C increases to v C = 0.75c s . For fracture network ruptures, v C is almost equivalent to v R,tot . In the mixed case of Ψ = 90 • with limited cascading rupture, v C is comparable to v R,ave . If cascading rupture across the fracture network does not occur (as for Ψ = 100 • − 120 • ), v C is comparable to v R,ave and v R,tot .

Average stress drop
We define the spatial-averaged stress drop as the vector integral over all points that ruptured and experienced a dynamic stress drop (Noda et al., 2013). The stress drop is computed from the difference between the initial and final shear stress after the rupture terminates. Our simulations produce an average stress drop ∆τ comparable to observations of crustal earthquakes (Huang et al., 2017). ∆τ for (sub-shear, non-cascading) main fault rupture is lower than 6 MPa (e.g., Ψ = 120 • , ∆τ = 5.9 MPa) while cascading ruptures generate higher ∆τ (e.g., Ψ = 65, ∆τ = 9.6 MPa, Table 2). Higher ∆τ can be associated with overall favorably oriented fractures (Figure 4), and with the occurrence of supershear v R,f on several fractures, comparable to observations in laboratory experiments (Passelègue et al., 2013). Lower ∆τ when rupture on -27-manuscript submitted to JGR: Solid Earth In our model, a relatively high stress drop facilitates rupture transfer across the fracture network.
Although considerable uncertainties in calculating stress drop from seismic observations exist (Abercrombie, 2021), relatively high stress drop may be observed in earthquakes across geometrically complex fault systems, for example, for the 1992 M w 7.3 Landers earthquake, a high stress drop (averaging ∆τ > 10 MPa) is observed and modeled (Kanamori et al., 1992;Wollherr et al., 2019)  We consider rupture initiation at ∼1 km fault-normal distance to mimic an off-main fault disturbance at deeper hypocentral depth than in earlier cases. The hypocenter is changed to x = −45 m, y = 1900 m, and z = −3430 m. We initiate rupture on a single prescribed fracture which has multiple intersections with neighboring fractures. To initiate the rupture, we use a nucleation radius of r nuc = 150 m, smaller than in previous cases (Section 3.2). We examine five scenarios based on the same 1 km distant off-main fault nucleation: (1) a reference scenario with the same hydrostatic pore fluid pressure condition (pore fluid pressure ratio γ = 0.37), as before, and (2) and four scenarios testing higher pore fluid pressure ratio(γ > 0.37), leading to decreased effective normal stress.

Case 4: rupture initiation on a fracture away from the main fault
The initial stress and fault/fracture strength conditions follow Case 2 (section 3.2.4, Ψ = 65 • ). The only difference is the location of the hypocenter which is now placed on a hanging wall fracture at ∼ 1 km distance normal to the main fault (red circle and star in Figure 10a). As in Case 2, we observe a pure rupture cascade that does not lead to a runaway rupture on the main fault (Figure 10b). The rupture duration is t = 3 s, hence slightly shorter than Case 2. However, the number of slipped fractures is higher (584 fractures slip), with more slipped fractures are located at the main fault's footwall.
-29-manuscript submitted to JGR: Solid Earth In this scenario, we apply initial stress similar to Case 2 (Section 3.2.4) but increase the fluid pressure ratio γ to above hydrostatic levels, hence, effectively decreasing normal stress. We note that increasing fluid pressure does not increase fault criticality in this framework, as the maximum pre-stress ratio R is kept constant across simulations. We assume uniformly distributed γ within the fracture network, across all fractures and the main fault. We consider four cases of overpressurized pore fluids: (1) γ = 0.5, (2) γ = 0.6, (3) γ = 0.7, and (4) γ = 0.8. When assuming γ = 0.5, the spatial slip distribution is comparable to using γ = 0.37 in Case 4 (Figure 11a), however there a fewer slipped fractures (569 out of 854, i.e. ∼ 70%).
Rupture activation is delayed within the footwall in the western part of the main fault at t = 2.5 s ( Figure   11a, Movie S12), leading to overall longer rupture duration (t = 4 s) compared to Case 4. The average stress drop also decreases as γ increases.

High-resolution seismic waveforms and spectral characteristics
Cascading and non-cascading ruptures are expected to result in different seismic waveform characteristics, which if observable, may provide insight into cascading rupture processes. We apply time-domain and Fourier spectral analyses to assess if there are notable differences in seismic-radiation properties.
For high-resolution wavefield modeling, we reduce the mesh element edge length to at most 250 m within a refinement volume (26 × 26 × 10 km 3 ). The resulting 88 million element mesh increases the computational cost of each forward simulation to 18h on 512 nodes (295,000 CPUh). This mesh resolves frequencies up to ∼ 7 Hz within the refined volume and for constant c s = 3464 m/s. Due to the statically adaptive mesh, resolution reaches up to ∼ 13 Hz at station 81, closest to the fracture network ( Figure   12). We analyze synthetic seismograms and Fourier amplitude spectra of ground velocity for three different scenarios using the high-resolution computational models: (1) cascading rupture within the fracture network  high-frequency radiation than the other two cases (noticeable already in the waveforms), which in contrast have higher energy in the low-frequency band. This suggests that cascading rupture may generate stronger high-frequency content due to the dynamic complexities governing its rupture processes, including abrupt rupture termination, acceleration, and deceleration of rupture fronts when branching and jumping across multiscale fractures.

Discussion
Our 3D dynamic rupture simulations in a geometrically complex fault network generate a rich set of results that raise a number of questions and implications. Below, we discuss conditions that lead to cascading rupture, implications of the connectivity and distribution of fractures, and overall source characteristics.

Static and dynamic conditions leading to cascading ruptures
Our study identifies at least three conditions that promote cascading dynamic rupture within the fracture network. First, the state evolution slip distance must scale with fracture and fault size to furnish the minimum fracture energy needed for rupture growth . Second, at least one fracture family must have a favorable relative pre-stress ratio (R ≥ 0.6), while the other family of fractures must have at least a conditionally favorable pre-stress ratio (0.3 ≤ R < 0.6). Higher R implies a more favorable orientation towards the ambient stress and initial shear stresses closer to critical. Third, fractures must be connected or densely packed to allow for sufficient stress transfer. Finally, when considering only one fracture family, the dynamic rupture cascade is suppressed ( Figure S11).
Cascading dynamic rupture can occur independent of the prescribed hypocenter location if all of the above conditions are fulfilled. The hypocenter may be located near the main fault or at the periphery of the fracture network and generate equally sustained rupture cascades. We show two hypocenter locations producing volumetric earthquakes of moment magnitude M w ≈ 5.5. The cascading earthquake magnitude is restricted by the dimension and distribution of the fracture network. Hence, rupture size may grow for a larger fault and its associated fracture network.

Multiple rupture fronts and re-nucleation
We observe that during fracture network rupture cascades, fractures may experience repeated nucleation while being ruptured by multiple rupture fronts, particularly if a fracture intersects with more than one other -33-manuscript submitted to JGR: Solid Earth -34-manuscript submitted to JGR: Solid Earth fracture (see Figure S12, Movies S3b, S4b, S5b, S6b). We identify three distinct dynamic mechanisms (1) near-simultaneous nucleation at two or more locations on a single fracture ( Figure S12a), (2) sequential nucleation at two or more locations on a single fracture (Figure S12b), and (3) repeated nucleation at the same point on a fracture ( Figure S12c). In the first mechanism, rupture fronts from two neighboring fractures simultaneously reach the edge of the same intersecting fracture. In the second mechanism, rupture initiates at different times in response to the timing of neighboring intersecting fractures slipping. The last mechanism occurs when two or more rupture fronts sequentially pass through the same fracture-fracture intersection.
Mechanisms 2 and 3 often occur together at the same fracture. In rare cases, for ∼ 1% of the ruptured fractures during a cascade, we observe more than two rupture initiations on a single fracture.

Fracture connectivity
Our simulations demonstrate that connected fractures facilitate cascading rupture. Compared to rupture jumping, rupture branching is significantly more effective: > 90% of all fractures that slip during a rupture cascade connect by direct branching. As a result, dynamic rupture cascades arrest when encountering unconnected fractures, even though self-sustained rupture can continue toward connected fractures ( Figure S12d). In our model, dynamic rupture jumping does not occur beyond an inter-fracture spacing of 80 m. This distance is observed as the largest fracture-fracture distance where wave-transmitted dynamic stresses from a neighboring unconnected fracture are sufficient to nucleate self-sustained rupture, i.e., sufficiently overstress an area matching the critical nucleation size of a distant unconnected fracture. While this distance depends on our frictional and geometrical model parameterization, we infer that densely spaced fractures favor dynamic rupture cascades.

Distribution of slip within a fracture network rupture cascade
The two conjugate fracture families produce left-and right lateral slip, with slip-direction (rake angle) 0 • and 180 • for fracture families 1 and 2, respectively. For main fault rupture-induced slip within the fracture network, slip accumulates on the dilatational side of the rupture propagation direction ( Figure   S13). This slip distribution resembles deformation patterns occurring in dynamic rupture simulations due to off-fault damage or off-fault plastic yielding (e.g., Dalguer et al. (2003); Ando and Yamashita (2007) Slip across the fracture network is driven by the intricate interaction of static and dynamic factors, which include dynamic stresses due to seismic waves, static Coulomb stresses due to evolving slip, and clamping and unclamping due to variations in normal stress. Hence, different realizations of fracture network dynamic rupture models regarding distribution, size, spacing, and connectivity of fractures may result in different slip patterns. However, we expect the major slip pattern characteristics of the cascading rupture to remain the same for a different but statistically similar fracture network model, such as in Gabriel et al. (2023), which includes a planar, vertical main fault under strike-slip loading.

Source characteristics of cascading ruptures
The equivalent point-source moment tensor representation of cascading and non-cascading ruptures changes from strike-slip to thrust faulting. Earthquakes across multiple large faults can promote non-DC moment tensor solutions by superposing distinct DC components (Julian et al., 1998;Palgunadi et al., 2020).
Combining several focal mechanisms may still produce a DC moment tensor solution (Julian et al., 1998).
For our set of cascading ruptures within the fracture network, we find that the equivalent moment tensors are characterized by insignificant non-DC components (∼ −1.8% to −5%). We note that both fracture families are conjugated, and therefore have similar average focal mechanisms. However, if the main fault is activated during cascading rupture a significant non-DC component emerges.
Multiple peaks in the MRF have also been attributed to other seismic source complexities either due to heterogeneous pre-stress on a planar fault (Ripperger et al., 2007), fractally rough fault surfaces (Shi & Day, 2013;Zielke et al., 2017;Tal et al., 2018;Danré et al., 2019), or non-uniform frictional parameters (e.g., variable characteristic slip distance (D c ) on a single fault plane (Renou et al., 2022)). Our study reveals that a small-scale fracture network can also generate multi-peak MRFs corresponding to multiple sub-events on fractures.
One of the main features of cascading rupture is the slow cascading speed v C , despite localized occurrences of supershear rupture speed within the fracture network. We hypothesize that the cascading speed v C may appear as the "true" single-fault rupture speed v R,tot if observed from a distance, i.e., observations of low rupture speeds may be at least partially explained by cascading rupture on a complex fracture network.
The 2019 M w 7.1 Ridgecrest earthquake is a good example. Its inferred low rupture speed (1.8 -2.0 km/s, (Chen et al., 2020), which is still faster than our pure cascading speed) could be affected by the fact that rupture propagated through a geometrically complex fault system and activated several off-fault fractures (Ross et al., 2019;X. Xu et al., 2020;Taufiqurrahman et al., 2023). We expect that v C may differ for other fracture-network geometries, distributions in fracture size and density. However, the analysis of alternative fracture-network configurations is beyond the scope of this study. Intuitively, we expect cascading rupture involving complicated fault geometries to generate lower v C than the surface-averaged rupture speed of the slipping fractures and main fault (v R,ave ), as rupture accelerates and decelerates during branching and jumping across many fractures.

Limitations
Due to the high computational demands of each simulation and the challenges of generating the 3D computational mesh accounting for multiscale intersecting fractures and faults, the results in this study are limited to one specific fracture network configuration, including fixed fracture orientation, distribution and size. Determining an exact fracture network parameterization for a particular rock volume or geo-reservoir is challenging. Often, statistical methods are used to constrain a fracture network because there is no reliable direct method for measuring the field-scale 3D fracture distribution. The geometry and distribution of fractures can vary significantly based on geological, tectonic, and mechanical factors. These factors likely affect the dynamics of a cascading earthquake. For example, removing one fracture family causes a cascading rupture to cease prematurely ( Figure S11).
For simplicity, we consider a homogeneous elastic-isotropic material. We here want to ensure that the subsurface structure does not affect the rupture dynamics and seismic radiation and does not mask the first-order physical signatures we intend to examine. However, natural fault zones comprise not only multi-scale fractures but also a damage zone around the fault core with lower elastic moduli and seismic wave speeds than the surrounding primary host rock, all of which may affect dynamic rupture and seismicwave radiation (Harris & Day, 1997;Huang & Ampuero, 2011;Huang et al., 2014). We do not account for co-seismically induced distributed damage or off-fault plastic deformation that may interact with discrete fractures (Andrews, 2005;Gabriel et al., 2013;S. Xu et al., 2015;Gabriel et al., 2021).
We also ignore aseismic and poroelastic processes that may additionally affect static and dynamic stress conditions in a fault zone, particularly for overpressurized fluid conditions (e.g., Segall and Lu (2015) (2022)), in our simulation framework will be readily possible in future work.

Conclusions
We present eighteen 3D dynamic rupture simulations within a complex fracture network of more than 800 intersected multiscale fracture planes surrounding a listric main fault. We vary prestress conditions, hypocenters, and fluid overpressure and analyze general aspects of seismic wave radiation. Our dynamic models reveal three mechanisms that promote cascading rupture: (1) the state evolution slip distance scales with fracture size , (2)

Acknowledgments
The authors thank the members of the Computational Earthquake Seismology (CES) group at KAUST for many fruitful discussions and suggestions. We thank SeisSol's core developers (see www.seissol.org). DIG acknowledges support by the Natural Sciences and Engineering Research Council (Discovery Grant 05743). Part of the analysis was implemented using Obspy (Beyreuther et al., 2010). Figures were prepared using Paraview (Ahrens et al., 2005) and Matplotlib (Hunter, 2007).

Data and Resources
The version of SeisSol used in this study is described in https://seissol.readthedocs.io/en/ latest/fault-tagging.html#using-more-than-189-dynamic-rupture-tags with commit version 917250fd.
Another alternative can be retrieved from SeisSol for hundreds of fault tagging in branch SeisSol64FractureNetwork
We use high-order basis functions of polynomial degree p = 4 achieving O5-accuracy in wave propagation in space and time for all simulations. Achieving high spatial and temporal resolution is crucial for resolving the detailed spatiotemporal evolution of the rupture processes governed by the variable process zone size in our frictional parameterization and the geometrically complicated fault planes, including the listric main fault geometry and numerous fault-fracture intersections. Since the process zone size varies considerably across our fractures, we measure and ensure to resolve the minimum process zone size following Wollherr et al. (2018). Our highest on-fault resolution is 4 m (for the smallest fractures of size 100 m), resolving the minimum cohesive zone of 8 m with on average 2 fifth-order accurate elements, that is 12 Gaussian integration points, ensuring highly accurate results. One high-resolution dynamic rupture simulation with the 45 million cell mesh requires 12.5h on 512 nodes on Shaheen II which is equivalent to 204,800 CPUh.

Appendix B Equivalent Moment Tensor Calculation
We determine an equivalent moment tensor representation for each dynamic rupture simulation. We calculate an element-local moment tensor following Lay and Wallace (1995). The equivalent moment tensor is then defined as the summation of all fault-local moment tensors of each fault element i.

Appendix C Mesh Generation
For this study, we construct the geometry of the fracture network using the third-party commercial software FRACMAN (Dershowitz et al., 2019). The boundary of the numerical domain is defined by employing the open-source mesh generator gmsh (Geuzaine & Remacle, 2009) in a Cartesian coordinate system. A large numerical domain of 40 × 40 × 20 km 3 is used to mitigate the effect of expected reflected waves from the absorbing boundary. We select the shortest fracture length to limit small-scale fracture intersections, whose size restricts the smallest time step width. We discretize the unstructured tetrahedral mesh using Simmodeler (Simmetric Inc., 2020). Mesh-element edge lengths vary depending on the size of the process zone to ensure convergence and to numerically resolve the dynamic fault strength drop behind the rupture front. Following Wollherr et al. (2018), the on-fault mesh-element edge lengths vary from 4 m for the smallest (R=100 m) to 45 m for the largest fracture size (R=500 m). We gradually increase the mesh-element edge size of the tetrahedral mesh by a factor of 6% away from each fracture and fault plane to save computational cost while avoiding reflection from the domain boundary. Consequently, the edge length on the main listric fault reaches a maximum size of 100 m due to its connected and intersected surface with small fractures.
In total, the computational domain is composed of 49 million volume elements for fifteen dynamic rupture simulations (Cases 1, 1a, 1b, 2, 2a, 3, 4, and 5 ). For Cases 6, 7, and 8, to improve the temporal resolution of seismic wavefield, we reduce the mesh element edge length to a constant value 250 m within refinement volume 26 × 26 × 10 km 3 , resulting in 88 million volume elements.

Appendix D Rupture Initiation
Dynamic earthquake rupture simulations are commonly initiated by assigning a small area on the fault as time-dependent overstressed or reduced in strength; this area is the predefined nucleation zone (or hypocenter). We apply a time-dependent over-stress centered at the hypocenter location selected for each scenario ( Figure Appendix D.1). We choose the nucleation radius based on the numerical solution provided by Galis et al. (2015). Given their mathematical expressions and our initial pre-stress loading conditions, the estimated nucleation radius (r nuc ) is 400 m. The time-dependent stress increase within nucleation area R nuc is calculated by increasing relative pre-stress ratio R 0 as where Ω(r ) is a Gaussian step function, r is the radius from the hypocenter, and S(t) is a smoothed step function. The Gaussian step function is given by Ω(r ) = ξ exp r 2 r 2 − r 2 nuc , ∀ r < r nuc Ω(r ) = 0, otherwise (D2) ξ is the initial pre-stress ratio inside the nucleation patch. We set ξ = 3. The smoothed step function is formulated as , for 0 < t < T T indicates the nucleation time when the overstress is applied, chosen here as T = 0.1s. We apply a similar nucleation procedure for Cases 4 and 5, but with a smaller nucleation size of 150 m (Section 3.4).

Appendix E Rupture Kinematics for Cases 4 and 5
Case 4 differs from Case 2 only in the nucleation procedure (on a single fracture vs. volumetric) and nucleation location (distant vs near main fault). In Case 4, the cascading rupture produces kinematic source parameters similar to those of Case 2, such as v R,ave and total moment magnitudes. Only slight variations are observed in the average stress drop (∆τ ) and the equivalent point source moment tensor (Table E1) with Case 2. Thus, even though the cascading rupture occurs at two different locations with identical Ψ, the point-source parameters are similar.
In the five Case 5 simulations, we increase the fluid pressure ratio γ, leading to decreasing v R,ave and ∆τ (Table E1)

Contents of this file
1. Figure S1: Fracture size distribution.
3. Figure S3: Slip only on the main fault.
5. Figure S5: Stereonet plot of slipped fractures overlain by relative prestress ratio (R) for different Ψ shown in lower hemisphere projection.
6. Figure  8. Figure S8: Snapshot focuses on the rupture front of the main fault without showing fractures for scenario Ψ = 120 • .
10. Figure  11. Figure S11: Slip distribution if only considering one fracture family.