Ground‐Motion Characteristics of Cascading Earthquakes in a Multiscale Fracture Network Open Access

Geological faults develop complex multiscale features during their evolutionary stages, ranging from fine‐scale topographic variations (Power and Tullis, 1991; Candela et al., 2010) to large‐scale branching and segmentation (Chester et al., 1993; Ben‐Zion and Sammis, 2003; Mitchell and Faulkner, 2009, 2012), or to smaller‐scale fracture networks and damage zones. These geometrical irregularities induce variations in stresses around the faults due to tectonic deformation and both aseismic and seismic events (Faulkner et al., 2010; Griffith et al., 2012). They play a significant role in fault mechanics and earthquake rupture dynamics, including nucleation, propagation, and termination (Ide and Aochi, 2005; Bhat et al., 2007; Dunham et al., 2011; Kyriakopoulos et al., 2019; Okubo et al., 2019), resulting in distinct radiated seismic wavefields compared to simple fault geometries.

Although numerous studies have examined the role of geometrical complexity in dynamic rupture and its influence on the seismic wavefield, two fundamental questions remain unanswered: (1) How do cascading ruptures on complex fault geometries influence the generation of high‐frequency seismic radiation? (2) How do cascading ruptures across multiscale fractures affect seismologically inferred corner‐frequency‐based static stress drop and near‐source ground‐motion characteristics? Addressing these questions is essential for advancing our understanding of fracture‐network cascading ruptures and for potentially improving seismic hazard assessment.

Complex fault geometries, such as fault bends, branches, or rough surfaces, create stress concentrations along the fault and may lead to high‐frequency seismic wave radiation. When a rupture propagates through these areas, it undergoes rapid changes in stress conditions, resulting in abrupt acceleration or deceleration (Madariaga, 1977; Shi and Day, 2013). Fault bends or stepovers may cause the rupture to change direction or split, resulting in fluctuations in rupture velocity and direction (Kame and Uchida, 2008). Off‐fault damage acts as an additional source of geometrical complexity and high‐frequency radiation, occurring over small spatial and temporal scales (Andrews, 2005; Okubo et al., 2019; Wollherr et al., 2019; Gabriel et al., 2021; Zhao et al., 2024).

Despite remarkable past efforts, observing and modeling the effects of complex fault geometry in detail remain a challenge. Common physics‐based simulations often simplify the source as one or a few infinitesimally thin, planar surfaces. To approximate fault‐zone complexity, modelers have used proxies such as lower rigidity zones (Huang and Ampuero, 2011; Huang et al., 2014), fault roughness (Roten et al., 2011; Shi and Day, 2013; Mai et al., 2018), and heterogeneous source properties (e.g., spatially varying rupture speed, slip velocity, rise time, and critical slip distance, $Dc$⁠) (Ide and Aochi, 2005; Ripperger et al., 2007; Gallović and Valentová, 2023). However, these proxies may not adequately capture the complete multiscale geometrical complexity observed in fault zones. The main difficulty lies in computational challenges and the lack of detailed fault models that explicitly incorporate multiscale fractures.

Over the past decade, facilitated by increased computation power and improved numerical methods, more advanced 3D dynamic rupture simulations that account for geometrical irregularity, spatial variability in fault friction and strength, and heterogeneous media have demonstrated that these complexities can generate ground motions with high‐frequency content (Withers, Olsen, Shi, and Day, 2019; Taufiqurrahman et al., 2022; Gallović and Valentová, 2023). Physics‐based simulations also suggest that the distribution of ground motions varies depending on earthquake rupture characteristics, highlighting the importance of conducting more realistic earthquake simulations, especially in the near‐source region (Ide and Aochi, 2005; Vyas et al., 2016; Ando et al., 2017; Wollherr et al., 2019).

Current applications of these advancements include simulations that capture high‐frequency ground motions up to 5–10 Hz in the near‐source region by considering fault and material heterogeneity (Shi and Day, 2013; Withers, Olsen, Day, and Shi, 2019; Taufiqurrahman et al., 2022; Vyas et al., 2024). Recent studies by Palgunadi, Gabriel, Garagash, Ulrich, and Mai (2024) and Gabriel et al. (2024) show that using dynamic rupture simulations and explicitly accounting for a geometrically complex fault network, moderate earthquakes may result from a cascading rupture across a fracture network, leading to distinct ground‐motion characteristics. Potential advantages of these models include a more accurate representation of rupture dynamics and improved calculation of ground‐shaking intensity, which are critical for seismic hazard assessments.

Seismological estimates of corner frequency and stress drop are often based on the assumption of a circular single‐plane rupture model (Brune, 1970; Madariaga, 1977). However, measurements, especially for moderate‐sized earthquakes (⁠$Mw<6$⁠), are not well resolved and contain significant random and potentially systematic uncertainties (Abercrombie, 2021). Unfortunately, real earthquake ruptures are more complex than the single‐plane rupture model due to geometrical and source‐related complexities. As indicated by their point‐source representations of the source time function, earthquakes with complex source properties exhibit distinct displacement spectra, such as intermediate‐frequency fall‐off (Abercrombie, 2021). Therefore, seismologically inferred corner frequency and stress‐drop estimates may only be applicable to earthquakes with relatively simple source geometries and properties. Investigating highly ideal cases using numerical modeling of fracture‐network cascading earthquakes may provide insights into the consequence of the seismological approach to determine more complex earthquake source properties.

In this study, we leverage advanced 3D dynamic rupture simulations that explicitly incorporate geometrically complex, multiscale fractures to investigate how high‐frequency ground motion is generated during cascading earthquakes. Unlike previous work, we model the fault network in detail, capturing realistic geometric features and using kinematic ground‐motion simulations to isolate the contributions of individual slipping fractures to the radiated seismic wavefield. We also analyze seismically inferred stress drops, investigate spatial variations in the frequency content of the near‐source wavefield, and evaluate ground‐shaking intensities. Our results are then validated against earthquake observations in central Italy. By incorporating explicit modeling of multiscale fault complexity, this study offers new insights into the influence of fault geometry on ground‐motion characteristics.

Our model is not specifically designed to replicate the central Italy earthquakes, but these events serve as excellent validation cases for our study. The fault systems in central Italy exhibit listric normal‐fault geometry with complex fracture networks analogous to our model. Furthermore, these earthquakes were recorded by a dense seismic network, providing high‐quality observational data that allows us to compare our theoretical findings with real‐world events.

In this section, we summarize our approach to generating ground motions from dynamic ruptures occurring within a fracture network. For complete descriptions of the fracture network construction and the rupture parameters, we refer the reader to Palgunadi, Gabriel, Garagash, Ulrich, and Mai (2024). To examine the individual contributions of each fracture to the overall seismogram, we employ a kinematic simulation approach that utilizes rupture parameters derived from the outputs of dynamic rupture simulations.

This study employs a fault model consisting of a listric fault with depth‐dependent dip angles ranging from 30° to 80° and a strike direction of N270E. The main fault measures 8 km along strike and 4 km along dip, with a burial depth of 1 km, ensuring the fault does not intersect the free surface (Fig. S1, available in the supplemental material to this article). We choose a buried fault to focus on interactions within the fracture network itself and to avoid additional physical mechanisms due to free‐surface interactions that further complicate the rupture dynamics. This type of listric fault geometry, common in transitional regimes between normal and strike‐slip faulting, is often associated with subsurface reservoirs of oil, gas, or geothermal energy (Onajite, 2013; Ward et al., 2016; Dixon et al., 2019). The fractures are confined to a volumetric damage zone of 10 km × 1 km × 6 km in the x (along strike, east–west), y (normal to the main fault, north–south), and z (depth, positive upward) directions (inset in Fig. 1a).

The density and size distribution of the fractures are statistically constrained by field observations using power‐law relationships (Mitchell and Faulkner, 2009; Savage and Brodsky, 2011). Fracture density is quantified as the number of fractures per unit length (⁠$P10$⁠), assumed to follow a power‐law decay characterized by a decay constant m = 0.8 and a coefficient C = 2.5 as a function of fracture distance (d) (⁠$P10=2.5d−0.8$⁠, Palgunadi, Gabriel, Garagash, Ulrich, and Mai, 2024). Fractures are confined to the designated network region, with densities decreasing from the vicinity of the main fault outward. The fracture sizes (R) follow a power‐law distribution with radii ranging from 100 to 500 m, decaying as $∝R−2$ (Lavoine et al., 2019). The orientations of fractures are constrained by the observed orientations of large‐scale fractures and conjugate faults within fault damage zones, typically ranging from 25° to 30° relative to the main fault (Mitchell and Faulkner, 2009). The network includes two dominant fracture families with strike orientations of N120E and N20E, consisting of 423 and 431 fractures, respectively. Among all fractures, 78% are connected to more than one fracture, 3% are connected to only one fracture, and 19% are unconnected (Fig. S1 or Fig. 1b).

To conduct 3D dynamic rupture simulations, we use the open‐source software SeisSol (see Data and Resources), which solves the dynamic rupture problem on a specified surface of potentially complex 3D geometry embedded in a 3D Earth‐structural model and subjected to an initial regional stress field. These simulations allow modeling of the complete space–time dynamic earthquake rupture process in a physically self‐consistent manner.

SeisSol employs fully adaptive, unstructured tetrahedral meshes, enabling the representation of geometrically complex geological structures (Pelties et al., 2014; Uphoff et al., 2017). Recent computational optimizations version used in this study have significantly enhanced the accuracy and scalability of SeisSol (Uphoff and Bader, 2020; Dorozhinskii and Bader, 2021), including through an efficient local timestepping algorithm (Breuer et al., 2016; Uphoff et al., 2017).

The computational domain for this study is a $40×40×20km3$ cube, discretized into an unstructured tetrahedral mesh with varying element edge lengths. We employ high‐order basis functions of polynomial degree p = 4, achieving fifth‐order accuracy (⁠$O5$⁠) in both space and time for seismic wave propagation. To ensure the resolution of the minimum breakdown zone size, the smallest fracture (100 m) is resolved with a maximum on‐fault resolution of 4 m, resulting in a minimum cohesive zone width of 8 m (following Wollherr et al., 2018). This is achieved using an average of two fifth‐order accurate elements, each with 12 Gaussian integration points. The mesh is refined within a smaller cube measuring $26×26×10km3$⁠, with the largest mesh element edge length set at 250 m, resulting in over 88 million tetrahedral cells. We choose a homogeneous medium to isolate the effects from rupture dynamics across a fracture network from any wave‐propagation effects. However, in reality, fault damage zones exhibit spatially varying material properties that differ from those of the surrounding host rock to mimic multiscale fracture networks (from hundreds to micrometer scale). These variations can generate distinct seismic wavefields, including reverberating waves within the damage zone and nonlinear interactions with rupture dynamics (Harris and Day, 1997; Duan, 2008; Huang and Ampuero, 2011; Huang et al., 2014). We model the medium as a homogeneous 3D half‐space with P‐wave velocity $VP=6.0km/s$⁠, S‐wave velocity $VS=3.464km/s$⁠, rock density $ρ=2670kg/m3$⁠, and rigidity $μ=λ=32GPa$⁠. This mesh configuration resolves frequencies up to 7 Hz everywhere within the refined volume and ∼13 Hz at all stations close to the fracture network.

We adopt a laboratory‐based rate‐and‐state friction law (Rice, 2006) with rapid velocity weakening (Noda et al., 2009; Dunham et al., 2011). The frictional parameters are held constant, except for the state evolution slip distance L, which scales linearly with fracture size (Gabriel et al., 2024). The model includes the direct effect parameter (a = 0.01) and the evolution parameter (b = 0.014), ensuring the main fault and fractures remain frictionally unstable. The characteristic weakening velocity $Vw=0.1m/s$ is based on the experimentally observed set of rapid decay of the effective friction coefficient, ranging between 0.01 and 1 m/s (Rice, 2006; Beeler et al., 2008; Di Toro et al., 2011). The steady‐state friction coefficient (⁠$f0$⁠) is fixed at 0.6 at a reference slip velocity $V0=10−6m/s$⁠, with a weakened steady‐state friction coefficient (⁠$fw$⁠) of 0.1 (Rice, 2006).

In this study, we focus on three specific cases described in Palgunadi, Gabriel, Garagash, Ulrich, and Mai (2024): case 1, a pure cascading rupture within the fracture network, is obtained by assuming a maximum horizontal stress orientation (⁠$Ψ$⁠) of 65°. Case 2, a rupture with off‐fault fracture slip and case 3, a main‐fault rupture without the fracture network in the model, are obtained with $Ψ=120°$ (see Palgunadi, Gabriel, Garagash, Ulrich, and Mai, 2024). Fault and fracture strength are characterized using the relative prestress ratio $R$ (Aochi and Madariaga, 2003; Ulrich et al., 2019). $R$ is defined as $(τ0−fwσn′)/((fp−fw)σn′)$⁠, in which $τ0$ is the initial shear stress, $σn′$ is the effective normal stress, and $fp$ denotes the peak friction coefficient during dynamic rupture, which here is approximated by its reference value $f0$⁠. A constant $R0=0.8$ is set for fractures at the maximum possible value of $R$⁠, with shear stress to effective normal stress loading ratios (⁠$τ0/σn′$⁠) ranging from 0.2 to 0.5, indicating subcritical stress conditions.

Case 1 features a favorable relative prestress ratio (⁠$R≥0.6$⁠) for the fractures and unfavorable relative prestress ratio (⁠$R<0.3$⁠) for the main fault (⁠$Ψ=65°$⁠); case 2 favorable $R$ for the main fault but unfavorable $R$ for fractures (⁠$Ψ=120°$⁠); case 3 has the same stress conditions as case 2 but does not incorporate the fracture network, hence, the rupture is restricted to the single main‐fault surface. Here, we focus our attention on the radiated seismic wavefield and the implications of the intricate rupture dynamics on the ground‐motion properties.

Throughout this study, we use the synthetic waveforms from the dynamic rupture models and analyze characteristic frequencies and other ground‐motion parameters. Furthermore, we perform kinematic ground‐motion calculations to explore the contribution of individual fractures to the overall radiated seismic wavefield.

To understand which slipped fracture contributed to which parts of the radiated seismic wavefield, it is necessary to decompose the complete wavefield. However, in dynamic rupture simulations, decomposition of the seismic wavefield into the individual contributions of the slipped fractures is impossible due to the inherent nonlinearity. Therefore, we employ kinematic rupture simulations as a means to efficiently decompose the dynamic rupture waveforms into individual contributions (e.g., case 1 in Fig. 2).

For the kinematic ground‐motion computation, we use the open‐source software AXITRA (Cotton and Coutant, 1997), which models wave propagation from the source to the receiver generated by a prescribed earthquake rupture model. AXITRA computes the stress field radiated by the six components of the moment tensor, employing the reflectivity method alongside discrete wavenumber decomposition Green’s functions for an axis‐symmetric medium (Bouchon, 1981). In our study, we apply the kinematic source parameters as obtained from the dynamic rupture simulations, including final slip (see Fig. 1), moment‐rate function, seismic moment, rake angle, and rise time (refer to Fig. 3, Fig. S2).

The Green’s functions assume a homogeneous medium, consistent with the dynamic rupture model, and target a frequency of 13 Hz. Each fault element (triangle in a tetrahedron), with a minimum slip of 0.001 m, is treated as a point source with its individual kinematic parameters and source time functions. The seismogram generated from the kinematic model integrates waveform contributions from all slipped fractures and the main fault. For the fracture‐network cascading rupture scenario, the kinematic source model incorporates ∼8 million point sources, corresponding to each slipped cell on the fractures and the main fault. For other scenarios, the kinematic rupture is sampled with 5 million and 1.5 million point sources, respectively, aligning with cases involving rupture with (case 2) and without (case 3) off‐fault fracture slip.

Following Schliwa and Gabriel (2024), we use Brune’s model to examine the spectral decay characteristic of the near‐field waveforms, recognizing that equation (1) has been developed under the far‐field assumption. However, as proposed in Schliwa and Gabriel (2024), the spatial variations of near‐field amplitude spectra may carry important information on the details of the rupture. We calculate $fc$ for each virtual receiver at the Earth’s surface to assess potential spatial variations in $fc$⁠. Analyzing the source‐spectral behavior from our simulations is subject to the complication that the duration of rupture typically exceeds the arrival‐time difference between P and S waves. Consequently, clearly distinguishing between P‐ and S‐wave spectral properties is often impossible, which presents a deviation from some of Brune (1970) underlying theoretical assumptions.

It is important to note that in the idealized setup of this study, neither intrinsic nor scattering attenuation are modeled, nor are site effects considered. Consequently, the spectral shape of the simulated ground motions reflects only source effects and geometrical spreading in a homogeneous half‐space. We explored both the original Brune model (⁠$γ=1$⁠) and the sharper Boatwright (1980) model (⁠$γ=2$⁠), finding that the former provides a better fit with our synthetic waveforms. We focus on station distances less than 15 km and determine $fc$ using the grid‐search method of Schliwa and Gabriel (2024), the best fit’s 95% confidence interval considered as the uncertainty in $fc$⁠.

The cascading earthquake generates pulses (i.e., alternating periods of slip) on the main fault and fractures in the macroscale view, distinguished by their short rise time compared to the rupture time (see Figs. S2a,b) (Heaton, 1990) and slipped fractures across the fracture network (see Fig. 1). The cascading rupture generates localized high slip patches at several fault–fracture or fracture–fracture intersections, as detailed in Palgunadi, Gabriel, Garagash, Ulrich, and Mai (2024). In contrast, the non‐cascading earthquake produces slip concentrated on the main fault and its directly connected fractures.

To analyze differences in seismograms or ground motions, we compare three scenarios: (1) a cascading rupture within the fracture network (case 1), (2) a rupture with off‐fault fracture slip (case 2), and (3) a main‐fault rupture without the fracture network being incorporated into the model (case 3), resulting in magnitudes of $Mw$ 5.5, 6.0, and 6.0, respectively. Figure 4 displays three‐component velocity seismograms for all three cases at five stations, respectively. Cases 1 and 2 are discussed in detail in Palgunadi, Gabriel, Garagash, Ulrich, and Mai (2024).

To compare the three cases, we use the outputs generated by dynamic rupture models. Figure 4a highlights significant differences in ground‐motion properties between fracture‐network cascading (case 1) and non‐cascading earthquakes (cases 2 and 3), as described in Palgunadi, Gabriel, Garagash, Ulrich, and Mai (2024). Notably, case 1 exhibits higher frequency content compared to cases 2 and 3, as illustrated in Figure 4b. The waveforms of case 1 are characterized by short‐wavelength amplitude variations across all components and stations (highlighted in red in Fig. 4a). Specifically, the east–west (fault‐parallel) component shows higher amplitudes than the north–south (fault‐normal) component. In addition, the waveforms in case 1 include late‐arriving waves that resemble “coda” waves typically generated by seismic scattering, here attributed solely to rupture (i.e., source) complexity and associated increased rupture duration compared to cases 2 and 3 (Fig. 2).

The seismograms for cases 2 and 3 display nearly identical low‐frequency signatures across all components. However, case 3 has lower frequency content than case 2, attributed to the dynamic rupture that occurs only on the main fault, without involving the fracture network in the model (depicted in blue in Fig. 4b). In both cases, higher amplitudes are observed on the vertical component (UD) at stations 2 and 81, located on the hanging wall of the listric main fault (Fig. 4a).

To investigate the origin of the source coda‐wave‐like signatures observed in case 1, we reconstruct the overall seismograms by combining contributions from each fracture. Dynamic rupture simulations are nonlinear, which does not permit the direct isolation of contributions from individual seismic sources (i.e., ruptured fractures) within a fracture network. Although isochrone analysis is possible for simpler fault geometries (Spudich and Cranswick, 1984; Schliwa and Gabriel, 2024), we here conduct kinematic simulations to isolate individual fracture contributions from the seismograms. This method is commonly used to generate the seismic wavefield from dynamic rupture simulations, for example, from dynamic rupture models based on the boundary integral element method (Aochi and Fukuyama, 2002; Goto et al., 2010) or for generating teleseismic synthetics (van Driel et al., 2015; Taufiqurrahman et al., 2023) and was applied by Ripperger et al. (2007), Mai et al. (2018), and Gallović and Valentová (2023). All elements on the slipped fractures and the main fault are treated as point sources with their respective source parameters.

To validate the kinematic approach, we present a comparison of seismograms generated directly in a dynamic rupture model (shown in red solid line) and through kinematic rupture (shown in blue dashed line) for case 1, specifically for near‐field stations across the north–south direction (red circle in the inset of Fig. 3). The kinematic approach well reproduces the wave characteristics observed in the dynamic rupture seismograms in the frequency range 0.001–6 Hz. Figure 2 illustrates a sketch of how different parts of the rupture radiate seismic waves, and how these then combine to form the complete seismograms at a site of interest (Fig. 3).

The first wave arrivals correspond to the P wave, marking the onset of the fracture‐network cascading rupture. The subsequent S‐wave arrivals are not attributed to specific fractures. The distinctive source coda‐wave‐like signatures arise from sequences of fractures that slip later in the shaking sequence, leading to a gradual decrease in amplitude as the event progresses.

The equivalent near‐field characteristic frequency (⁠$fc$⁠, Schliwa and Gabriel, 2024) is calculated using Brune’s model, as described earlier in the Equivalent near‐field characteristic frequency section. This frequency is crucial for characterizing earthquake sources by analyzing seismic waveform spectra (Kaneko and Shearer, 2014; Trugman and Shearer, 2017). We assess and compare the waveform spectra and $fc$ across the cascading and non‐cascading ruptures, providing statistics and spatial distributions of $fc$⁠.

Figure 5 shows $fc$‐values for the fracture‐network cascading earthquake (case 1) for the vertical component; EW and NS components are depicted in Figures S3 and S4 (stations within 4.8–5.2 km), respectively. On average, case 1 generates a broad $fc$ distribution with a mean of 1.1 Hz and a standard deviation of 0.5 Hz (Fig. S5a). Comparable mean and standard deviation values are observed across all components. In contrast, cases 2 and 3 exhibit a narrower distribution of $fc$‐values, predominantly below 0.5 Hz (0.21 and 0.17 Hz). Case 2 shows a higher average $fc$ than case 3, with standard deviations of 0.09 Hz (Fig. S5b) and 0.07 Hz (Fig. S5c), respectively. Considering the rupture durations from dynamic rupture of cases 1, 2, and 3, which are 3.5, 2.6, and 2.5, respectively (Palgunadi, Gabriel, Garagash, Ulrich, and Mai, 2024), the corresponding far‐field corner frequencies can be estimated using the relation $fc∼1/T$⁠, in which T is the rupture duration. This yields approximation values of 0.28, 0.34, and 0.4 Hz for the three cases. These $fc$ estimates are relatively comparable, except for case 1.

Despite these variations, the mean $fc$‐values from all three cases align with estimates based on Allmann and Shearer (2009) (Fig. 5b). The seismologically inferred stress‐drop estimates for cases 2 and 3 closely match the mean values from the simulations. For the cascading rupture case 1, the mean $fc$ is higher than the typical value for a moment magnitude of $Mw$ 5.5. It corresponds to an estimated stress drop (equation 3) of ∼145 MPa, which is almost an order of magnitude higher than the average stress drop directly calculated on the fractures and fault, which is 9.6 MPa (Palgunadi, Gabriel, Garagash, Ulrich, and Mai, 2024). This result suggests that in the near‐field region, spectral characteristics are sensitive to multifracture rupture behavior. Hence, the average stress‐drop calculation estimated by Brune’s model across all stations is biased by slip on the fractures, and associated high $fc$ observed at stations with azimuths aligned with the fracture trends. Although the average stress drop calculated on the fractures and main fault is 9.6 MPa, on a local scale, some patches on a fracture experience a maximum stress drop of up to 12 MPa, especially at the intersections between two or more fractures (Palgunadi, Gabriel, Garagash, Ulrich, and Mai, 2024).

For case 1, optimizing the high‐frequency fall‐off parameter (n) in equation (1) yields lower values of $fc$ compared to using a fixed value of n = 2 (Fig. S6). The average $fc$ across all virtual stations is 0.44 ± 0.15 Hz (blue dot in Fig. 5b), corresponding to an average seismologically inferred stress drop of $Δσ=9.28MPa$⁠.

We extract $fc$‐values within a 4.8 and 5.2 km distance from the epicenter for all azimuths to examine their azimuthal dependence. The highest $fc$‐values, reaching up to 5 Hz, are observed in the vertical component near the epicenter. $fc$ demonstrates an azimuthal dependence (Kaneko and Shearer, 2015), with higher frequencies in the direction of the slipped fracture network (Fig. 5c). Some $fc$ estimations, particularly those exceeding 4 Hz, show high uncertainties (⁠$Δfc∼2Hz$⁠). On average, azimuths of around 20°, 120°, 220° , and 300° correspond to high $fc$‐values (⁠$fc∼4Hz$⁠). These azimuths approximately align with the fault planes of the equivalent point‐source moment tensor solution, characterized by pure strike‐slip faulting mechanisms (strike = 114°, dip = 88°, rake = 0°) in case 1 (table 2 in Palgunadi, Gabriel, Garagash, Ulrich, and Mai, 2024), with $Ψ=65°$⁠).

To put our corner‐frequency results into perspective with detailed observations, we analyze earthquake recordings from the three largest earthquakes of the well‐documented recent central Apennines earthquake sequence, in Italy: the 2009 $Mw$ 6.2 L’Aquila earthquake, the 2016 $Mw$ 6.2 central Italy, Amatrice, earthquake, and the 2016 $Mw$ 6.6 Norcia earthquake. These events were selected due to their normal‐fault geometry with complex fracture networks, the dense distribution of strong motion seismometers, and the availability of freely accessible data (see Data and Resources). To calculate the corner frequency, each seismogram is analyzed using Brune’s model. All three components are treated separately and then averaged. We identify an azimuthal dependence of the equivalent near‐field characteristic frequency for all three earthquakes: the Norcia earthquake (see Fig. 6), the Amatrice earthquake (Fig. S7a), and the L’Aquila earthquake (Fig. S7b). Despite variations in the epicentral distances of the stations, higher corner‐frequency values tend to be oriented toward the fault strike. For these earthquakes, the corner‐frequency ranges from 0.1 to 0.85 Hz, with an average of ∼0.3–0.4 Hz (Maercklin et al., 2011, Calderoni et al., 2017; Supino et al., 2019).

Based on the synthetic data from case 1 and observed strong motion data, the radiation pattern and forward directivity may lead to an increase in the corner frequency or $fc$ in the fracture direction, at least for stations within an epicentral distance of less than 100 km. This observation is particularly emphasized for stations close to the earthquake source and in agreement with large‐scale dynamic rupture simulations results in Schliwa and Gabriel (2024).

We analyze the relationship between the equivalent near‐field characteristic frequency (⁠$fc$⁠) and source complexity. Figure 7 illustrates that $fc$‐values vary across different azimuths, with higher frequencies observed in the dominant fracture strike directions. The higher frequencies inferred from the vertical component correlate with the SH‐wave radiation pattern (Fig. 7), representing the nodal lobes in which seismic energy consists of incoherent radiation scattered by the randomized fracture orientations. The relatively high $fc$ persists up to a resolved distance of 12.5 km across all azimuths. Conversely, the distribution of $fc$‐values in the east–west and north–south components appear less pronounced.

The azimuthal variation of $fc$ remains consistent even when n is allowed to vary during the inversion (Fig. S6). The primary difference is a systematic reduction in corner frequency across all azimuths, with an average decrease of ∼0.65 Hz. Interestingly, n also exhibits azimuthal dependence, tending to be higher in the direction of the fracture strike. On the vertical component, n shows greater variability, ranging from 1 to 3, whereas on the horizontal component, n is generally greater than 1.8.

High‐resolution dynamic rupture simulations for the 1992 $Mw$ 7.3 Landers earthquake and the 2019 $Mw$ 7.1 Ridgecrest earthquake also reveals spatial variations in corner frequency or $fc$ (Schliwa and Gabriel, 2024). These two earthquakes were hosted by geometrically complex fault segments in the primary large‐scale fault network. $fc$ variations are caused by the radiation pattern, rupture directivity, and local dynamic source effects, for example, facilitating the identification of rupture gaps and rake rotations due to free‐surface interaction (e.g., Kearse et al., 2019). Our findings further emphasize the consistency across various seismic events for using the near‐field $fc$ variability for source characterization.

The presence of two pairs of peaks corresponding to nearly orthogonal fracture families indicates the characteristics of the cascading rupture across these two dominant fracture families. This $fc$ pattern is also reflected in the moment tensor solution of case 1, which is misaligned with the strike of the main fault (Gabriel et al., 2024; Palgunadi, Gabriel, Garagash, Ulrich, and Mai, 2024).

We apply a site‐amplification correction to the ground‐motion synthetics before comparing our simulations with empirical ground‐motion models (GMMs). This correction accounts for the fact that our simulations are conducted in a homogeneous high shear‐wave velocity (⁠$VS$ of 3.464 km/s), in contrast to data‐based GMMs that include site conditions that most often are characterized by much lower seismic wave speeds (⁠$VS30<1.5km/s$⁠). We apply an amplification factor to all seismograms to scale amplitudes to represent an assumed $VS30$‐value of 0.76 km/s, following the methodology outlined by Borcherdt (1994, 2002). This approach has been applied in previous studies (e.g., Mai et al., 2010; Moczo et al., 2018). Figure S8 illustrates the filter and its effect on ground‐motion waveforms. The filter amplifies ground motion in the frequency range of 0.1–3 Hz, while slightly decreasing amplitudes above 3 Hz (Figs. S8a,b). The waveform differences before and after correction are minimal (Fig. S8c).

We then compare the peak ground velocity (PGV) and peak ground acceleration (PGA) of 81 stations against the predictions of selected GMMs (Fig. 8) by Abrahamson et al. (2014), Boore et al. (2014), Campbell and Bozorgnia (2014), and Chiou and Youngs (2014). The three cases exhibit distinct peak amplitude characteristics. In case 1, the PGA values are generally higher than the median GMM predictions but remain within the GMMs’ standard deviation. Some values, particularly at stations very close to the main fault (Joyner–Boore distance, $RJB<5km$⁠), tend to overestimate the GMMs. In case 2, the PGA values are slightly below the GMMs’ median value for an $Mw$ 6.0 (represented by dashed lines in Fig. 8a). In case 3, the PGA values generally fall below the median GMM predictions, indicating an overall underestimation of the synthetics.

For PGV, case 1 aligns well with the median GMM value for $RJB>5km$ and underestimates GMMs for $RJB≤5km$⁠. Cases 2 and 3 yield comparable PGV values compared to GMMs. Both cases show two clusters of PGV values, corresponding to ground shaking on the hanging wall and footwall, respectively. The lower PGV values may be attributed to the fact that our simulations consider a homogeneous velocity model and small‐scale source effects, whereas GMMs incorporate data recorded under a wide range of geological conditions. Figure 8 confirms that the ground motions from our 3D dynamic rupture simulations are in general agreement with empirical predictions. Therefore, we proceed to discuss directional patterns of shaking and further ground‐motion properties.

We assess PGA and PGV values in different azimuthal directions and compare these with GMMs to examine azimuthal dependencies of shaking levels (Fig. 9). The azimuth angles range from 0° to 360°, increasing clockwise, with 0° corresponding to north. In case 1, a clear distinction is observed between the violet and green dots, representing azimuthal angles of ∼300° and 120°, respectively. The PGA and PGV values in the rupture direction (120°, green squares) tend to overestimate the GMMs for $RJB>5km$⁠, whereas the region corresponding to the backward rupture direction aligns well with the GMMs (violet squares) and is underestimated for $RJB<5km$⁠. In cases 2 and 3, the PGA values generally underestimate the GMMs, except in specific azimuthal directions (110°–130°) in which simulated ground motions exceed GMM predictions at distances beyond 10 km. The PGV values show good agreement with the GMMs in the hanging‐wall region.

Ground motions exhibit complexity due to variations in the spatiotemporal evolution of ruptures, as demonstrated in various studies (e.g., Mai, 2009; Dunham et al., 2011; Shi and Day, 2013). The cascading rupture, involving complicated dynamic interactions within the fracture network (Palgunadi, Gabriel, Garagash, Ulrich, and Mai, 2024), results in ground motions that significantly differ from those of noncascading ruptures. Figure 10 presents maps comparing PGV using GMRotD50, a metric independent of the sensor orientations (Boore et al., 2006) for three scenarios: cascading rupture (case 1, Fig. 10a,d), rupture with off‐fault fractures slip (case 2, Fig. 10b,e), and rupture without the fracture network in the model (case 3, Fig. 10c,f).

Case 1 exhibits greater variability in PGA and PGV amplitudes across all azimuths compared to cases 2 and 3. The distributions for cases 2 and 3 indicate pronounced shaking on the hanging‐wall side of the listric fault, corroborating findings by Abrahamson and Somerville (1993) and Somerville et al. (1996). This observation is consistent with kinematic ground‐motion simulations along listric faults (Ofoegbu and Ferrill, 1998; Passone and Mai, 2017; Rodgers et al., 2019; Moratto et al., 2023).

The variability in PGA for case 2 is largely attributable to off‐fault fracture slip, evident from the relatively higher PGA values on the footwall side compared to case 3 (Fig. 10b,c). Although both cases 1 and 2 involve multiple fractures, the amplitude variability of PGA between the two cases is significant. Case 1, with 561 slipped fractures representing 70% of the total fractures, exhibits higher PGA values. In contrast, case 2, with 477 slipped fractures accounting for 58% of the total, shows less variability. This less than 20% increase in the number of slipped multiscale fractures results in a substantial rise in PGA values. This dramatic difference in PGA is mainly due to a significant increase in rupture speed (⁠$Vr$⁠) within individual fractures (Fig. S9). In case 1, the rupture speed is predominantly $Vr>2900m/s$ (Fig. S9a), whereas case 2 features slower $Vr$‐values (Fig. S9b). According to the kinematic and dynamic approach, $Vr$ plays a major role in producing high peak amplitudes, especially in the near‐field region (Robinson et al., 2006; Zollo et al., 2006; Bizzarri and Spudich, 2008).

Our comparative study of cases 1, 2, and 3 reveals significant differences in the azimuthal dependence of ground‐motion intensities for the near‐field region, highlighting ground‐motion patterns between cascading and non‐cascading ruptures. To further analyze the distribution of PGA and PGV values across all azimuths, we calculate the mean and standard deviation for every 10° azimuthal bin for an epicentral distance of <12.5 km. In Figure 11a,b, the red color represents case 1, which shows no apparent azimuthal dependence of PGA and PGV values. Instead, they remain relatively constant across all azimuths. Case 1 suggests that ground motions are uniformly distributed across all azimuths. In contrast, cases 2 and 3 exhibit clear azimuthal dependence and demonstrate similar patterns. The green and blue colors in Figure 11a,b indicate high PGA and PGV values on the hanging wall (270° < azimuth ≤ 360° and 0° ≤ azimuth <90°) and low values on the footwall (90° ≤ azimuth ≤270°). Interestingly, the standard deviation of case 1 is significantly higher on the footwall side compared to cases 2 and 3 for both PGA and PGV.

Pseudospectral acceleration (PSA) is a widely used measure for assessing the severity of ground motion on structural systems in earthquake engineering. Based on single‐degree‐of‐freedom spring‐mass‐damping systems, it significantly aids in estimating the peak response of structures to seismic events (Baker et al., 2021). We calculate PSA values with a damping factor of $ζ=5%$ over a period range of T = 0.1–10 s. We then compare these PSA values from the three scenarios against several GMMs (Abrahamson and Silva, 1997; Boore et al., 1997, 2014; Campbell and Bozorgnia, 2014).

Figure 12a compares our obtained PSA values against the GMMs. Case 1 PSA values overall align with the GMMs’ median values, but exhibits higher values at shorter periods and lower values at longer periods, compared to the GMMs. Specifically, case 1 shows elevated PSA values for shorter periods (T = 0.1–0.2 s), indicative of intense high‐frequency ground‐motion radiation. In contrast, cases 2 and 3 demonstrate higher PSA values at longer periods and lower PSA values at shorter periods relative to the median GMM values. This pattern is highlighted in Figure 12a, with case 2 (green dots) showing higher PSA values at short periods compared to case 3 (blue squares). However, at longer periods (T > 1 s), PSA values for cases 2 and 3 converge, suggesting that seismic wave radiation from the main fault predominates at these periods, whereas contributions from slipping off‐fault fractures are more significant at shorter periods.

Figure 12b–d illustrates the distance‐dependent characteristics of PSA for each case over various periods. PSA values for the three cases align well with respective GMMs (⁠$Mw$ 5.5 for the fracture‐network cascading earthquake and $Mw$ 6.0 for cases 2 and 3). At longer periods, the PSA values for all cases are comparable and fall within the median estimates provided by the GMMs.

Figure 13 explores the distribution of PSA values at T = 0.3 s, shedding light on the extent of ground shaking across different azimuths. In case 1, significant ground shaking is evident across all azimuths, much like the values for PGA and PGV in Figure 10. PSA values at T = 0.3 s exhibit relatively constant levels across all directions, as depicted by the red dots in Figure 11c, with consistent standard deviations.

For cases 2 and 3, intense ground shaking primarily occurs on the hanging‐wall side of the main fault. Notably, case 2 also displays elevated PSA values to the southeast of the epicenter, reflecting the activation of off‐fault fractures within the damage zone. This spatial distribution underscores the impact of off‐fault fracture activation in generating high ground‐motion amplitudes across various azimuths. Both cases 2 and 3 show higher PSA values on the hanging‐wall side and reduced values toward the footwall. Intriguingly, case 2 exhibits higher PSA values on the footwall side than case 3, presumably owing to the additional off‐fault slipping fractures (represented by green triangles in Fig. 11c).

Our high‐resolution 3D dynamic rupture cascade simulations within a geometrically complex fracture network provide valuable insights into the characteristics of high‐frequency seismic wave generation and radiation. In this section, we discuss the properties of seismic waves resulting from fracture‐network cascading ruptures, exploring their implications, significance, and inherent limitations.

Cascading earthquakes within a fracture network generate high‐frequency seismograms across all components and distances, as demonstrated in our near‐field numerical simulations. However, for real earthquakes, this high‐frequency content is further influenced by complex Earth structures such as medium heterogeneity, ground‐surface topography, and fault roughness, echoing findings from previous studies (Bouckovalas and Papadimitriou, 2005; Sato et al., 2012; Imperatori and Mai, 2015; Vyas et al., 2021; Park et al., 2022; Taufiqurrahman et al., 2022; Gallović and Valentová, 2023). In practical observations, discerning whether an earthquake has occurred as a rupture cascade remains challenging due to Earth’s heterogeneities that obscure the origins of high‐frequency radiated energy.

The cascading rupture is marked by an indistinct onset of S waves across all components, attributed to the concurrent arrival of multiple P and S waves generated by slip on the various fractures, as demonstrated by our kinematic approach. Because these waves propagate over longer distances, the distinction between wave types becomes clearer owing to increased travel times (station 81 versus 45 in Fig. 4).

Our numerical simulations (both dynamic and kinematic approaches) indicate that fracture‐network cascading earthquakes generate seismic waves with higher equivalent near‐field characteristic frequencies (⁠$fc$⁠). The distribution of $fc$‐values across different azimuths serves as an indicator for identifying the characteristic behavior of fracture‐network cascading earthquakes. Our findings show that $fc$‐values demonstrate an azimuthal dependence, aligning with the radiation pattern of the fracture direction (Trugman et al., 2021; Schliwa and Gabriel, 2024). Although possible as in this study, confirming this azimuthal dependence in actual observations which have limited recording networks poses practical challenges due to 3D scattering in heterogeneous media, topographic effects, site‐specific conditions, and intrinsic attenuation.

To further validate our findings on the azimuthal variations of $fc$ as a proxy for corner frequency, employing near‐field measurements from a station array of diverse azimuths would be beneficial, such as those available from dense seismic networks like Italy’s strong motion network (see Fig. S7), Türkiye’s AFAD network that recorded the 2023 Kahramanmaraş earthquakes, Japanese (Natural Research Institute for Earth Science and Disaster Prevention [NIED], 2019) and U.S. National (Central Geological Survey [CGS] and U.S. Geological Survey [USGS], 2005) strong motion networks. Temporary dense strong motion arrays could also provide valuable data. Moreover, because seismic instrumentation becomes increasingly denser and more widely available, observations of spatially varying corner frequencies could provide critical insights into rupture propagation and directivity effects, thereby enhancing our understanding of earthquake source mechanics.

The fracture‐network cascading rupture is also characterized by a higher average $fc$⁠. The far‐field corner frequency estimated using seismologically inferred methods (Brune’s model spectra) is considerably higher than expected based on moment magnitude scaling relations (Fig. 5c). Consequently, this higher average $fc$ corresponds to a high stress drop, which is estimated to differ by almost an order of magnitude compared to the average stress drop calculated directly on the fault from the dynamic rupture simulation. In real earthquake events, the direct on‐fault estimation of stress drop is not possible. Nevertheless, our simulations suggest that using the conventional single‐plane rupture model is not appropriate for complex earthquake sources like rupture cascades. In the conventional model (i.e., seismologically inferred $fc$⁠), the stress drop is proportional to $M0/R3$ (Eshelby, 1957), in which the source radius $R∼VS/fc$⁠. In this relation, the estimation of $M0$ is presumably valid, whereas the inferred source radius is more representative of a single fracture or fault plane, not the collective rupture of a fracture‐network cascade involving hundreds of fractures. Nevertheless, this high average $fc$ and stress drop using Brune’s model estimation could be potentially observable features of fracture‐network cascading ruptures.

Considering high‐frequency fall‐off (n) as an optimized variable during the inversion in equation (1) yields a considerably lower average $fc$ than fixing n = 2. Including n as a free variable in the inversion process produces stress‐drop estimations that are closer to the on‐fault estimates. This relationship holds at least for the model we considered in this study, in which the material is homogeneous and complexity is solely related to fault geometry.

Treating the high‐frequency fall‐off parameter (n) as an additional free variable in the inversion yields a lower average $fc$ compared to fixing n = 2 (Fig. S6). Incorporating n as a free parameter also results in stress‐drop estimates that are more consistent with on‐fault dynamic stress‐drop values. This relationship, where treating n as a free parameter leads to lower $fc$ estimates, holds under the assumption of a homogeneous medium, with source complexity primarily governed by fault geometry and rupture dynamics. To further validate this relationship, it is important to account for material heterogeneity within the fracture network and to investigate its impact on ground‐motion characteristics.

In the context of a point‐source approximation, $fc$ serves as an indicator of source duration and typically scales as 1/T. The rupture durations in our cases 1 to 3 are 3.5, 2.6, and 2.5 s, respectively (Palgunadi, Gabriel, Garagash, Ulrich, and Mai, 2024), corresponding to expected corner frequencies of ∼0.286, 0.385, and 0.4 Hz. However, despite the longer rupture duration and lower expected corner frequency in case 1, the cascading rupture generates higher‐frequency wavefields. These results suggest that Brune’s model, based on far‐field assumptions and a simplified circular crack geometry, may not be appropriate for near‐field estimates of cascading rupture events due to: (1) the lack of clear separation between P and S waves, and (2) the presence of a geometrically complex source. The observation of a higher average $fc$ combined with longer rupture duration—this may indicate the influence of cascading rupture along a complex fault geometry. In contrast, the corner‐frequency estimates for cases 2 and 3 are more consistent with the expected values, suggesting better alignment with the assumptions of the Brune model. A more detailed evaluation of stress‐drop estimates under near‐field conditions remains beyond the scope of this study.

The ground shaking generated by fracture‐network cascading and non‐cascading ruptures exhibits distinct spatial amplitude patterns. Cascading ruptures result in heterogeneously distributed strong ground shaking across all azimuths, whereas non‐cascading ruptures typically generate intense shaking concentrated on the hanging‐wall side of the listric fault. Because the orientation of maximum ambient horizontal stress (⁠$Ψ$⁠) plays a critical role in dictating the rupture process, it also determines the actual shaking values and their distributions given a specific fault and fault system under consideration. $Ψ$ is the regional background stress state under natural conditions, without a specific direction with respect to a particular fracture orientation. In essence, for improved ground‐motion estimation, the regional stress field must be considered. However, accurately determining local stress orientations remains challenging due to significant uncertainties in current estimates (Heidbach et al., 2018). As indicated by Palgunadi, Gabriel, Garagash, Ulrich, and Mai (2024), a slight 20° difference in $Ψ$ can shift the rupture style from cascading to non‐cascading within the fracture network.

Fracture‐network cascading ruptures radiate more high‐frequency seismic energy compared to ruptures with off‐fault fracture slip, suggesting that current GMMs might underestimate PGA values, especially near faults surrounded by a multiscale fracture network. Our simulations indicate that the average PGA in cascading cases is nearly twice as large as that of non‐cascading cases, highlighting the limited effectiveness of non‐cascading ruptures in radiating high‐frequency waves. In addition, the activation of fractures limits the extent and affected area of high PGV values. Therefore, incorporating measurable geometrical fault complexity may improve the accuracy of ground‐shaking estimates.

The presence of strong high‐frequency shaking in the near‐field region raises concerns about the structural safety of buildings and facilities located near fracture networks, particularly in highly fractured fault zones close to population. Such fault networks may be activated by local stress perturbation, thus increasing the hazard. For instance, Türkiye’s $Mw$ 7.7 Kahramanmaraş earthquake, where the rupture started on an off‐fault segment, located ∼15 km east of the east Anatolian fault, and generated a devastating earthquake (Jia et al., 2023; Mai et al., 2023).

With the availability of increasing computational resources, advanced numerical methods, and efficient modern codes, it has become feasible and affordable to integrate dynamic earthquake simulations into physics‐based shaking assessments for seismic hazard analysis. However, this practice is not yet widespread. Adopting a rigorous physics‐based approach that accounts for complex fault geometry, topography, and medium heterogeneity will help significantly improve ground‐shaking estimates.

All results obtained from our simulations are contingent upon the model assumptions detailed in Palgunadi, Gabriel, Garagash, Ulrich, and Mai (2024). We acknowledge that specific findings, such as the spatial distribution of shaking or equivalent near‐field characteristic frequency estimates, may vary with different configurations of fracture networks and fault models. Nonetheless, the fundamental physics underlying the dynamic interactions among fractures will remain consistent.

Although our model is specifically designed for moderate‐size earthquakes, our findings may also offer insights relevant to larger earthquake studies. In particular, they show how fault‐zone geometrical complexity can either facilitate or hinder rupture growth by changing the orientation of the ambient horizontal stress (Gabriel et al., 2024; Palgunadi, Gabriel, Garagash, Ulrich, and Mai, 2024). In larger events (e.g., $Mw>6$⁠), the dynamic rupture processes are typically more complex, involving additional physical mechanisms such as thermal pressurization, as discussed in Gabriel et al. (2024). Because this is the first study to explicitly incorporate individual dynamic slip on off‐fault fractures, our objective is to isolate the effects of geometric complexity within the fracture network from those of other physical processes.

Our study primarily focuses on the geometrical complexity of faults while deliberately excluding several other factors. These include medium heterogeneities, low‐velocity zones, milli‐to‐microscale fractures, topography, off‐fault plastic deformation, fault roughness, site effects, and rupture dynamic interactions between faults and the free surface.

Seismic tomography studies of fault‐zone structures reveal significant variations in material properties relative to the surrounding host rock. Fault zones typically exhibit lower rigidity because of multiscale fracture networks, especially at millimeter‐to‐micrometer scale, which cannot be modeled explicitly. This can lead to high‐frequency oscillations in slip velocity near the rupture front (Harris and Day, 1997; Schliwa et al., 2025) and influence rupture speed (Huang et al., 2014). These low‐velocity zones also give rise to reverberating waves that may result in multiple slip pulses and shorter rise times (Huang and Ampuero, 2011). Incorporating a low‐velocity zone or material with reduced rigidity into our models could produce even more complex dynamic rupture behaviors, complicating interpretation. Such complexities might include additional rupture front episodes due to head‐wave reflections onto fracture surfaces, increased instances of localized supershear rupture along fractures, altered timing of fracture activation or arrest, and possibly enhanced high‐frequency content in the resulting seismic wavefield.

Our models do not consider fault and fracture networks that intersect the free surface. This negligence can be significant because free‐surface effects may substantially alter rupture dynamics and the resulting seismic wavefield. For instance, interactions between the free surface and cascading rupture may enhance local dynamic effects on individual fractures and the main fault, such as rake rotations (Kearse et al., 2019), potentially trigger a transition to supershear rupture (Kaneko and Lapusta, 2010), or even generate higher ground motions in the hanging wall (Oglesby et al., 2000).

Another shortcoming of our study is the lack of observational data for spontaneous cascading fracture‐network ruptures, which impedes more direct comparisons of our simulation outcomes with real earthquakes. A potential case study that may provide insights into such phenomena is the 2019 $Mw$ 7.1 Ridgecrest earthquake. However, the sparse distribution of seismic stations in the near‐field region limits our capacity to investigate the signatures characteristic of fracture‐network cascading earthquakes thoroughly. The search for the fracture‐network cascading earthquake may be possible in a well‐recorded moderate earthquake, such as induced or triggered earthquakes in a geo‐reservoir.

In this study, we investigate the near‐source ground‐motion properties of three high‐resolution 3D dynamic rupture simulations within a geometrically complex fault zone featuring a main listric fault and over 800 multiscale fractures. The three analyzed cases include a cascading rupture within the fracture network, a main‐fault non‐cascading rupture with secondary slip in the fracture network, and a main‐fault non‐cascading rupture without a fracture network. Our results demonstrate that cascading ruptures within fracture networks generate high‐frequency seismic radiation due to simultaneous slip on multiple fractures. The incoherence of these dynamic rupture processes produces source coda‐wave‐like signatures. The inferred high‐resolution near‐field characteristic frequency (⁠$fc$⁠) and stress drop are significantly elevated in the fracture‐network cascading rupture cases, approximately an order of magnitude higher than the estimation directly on the fault of the dynamic rupture simulation. In contrast, the non‐cascading main‐fault rupture scenarios show comparable $fc$ and stress‐drop values. The fracture‐network cascading ruptures show distinct azimuthal dependence of $fc$ for epicentral distances <15 km, with locally increased $fc$ of vertical components reflecting the complex radiation pattern and rupture directivity. These results are consistent with observations from large‐scale geometrically complex dynamic rupture simulations of real earthquakes and the well‐recorded 2009–2016 central Apennines, Italy, earthquake sequence. We identify several characteristics highlighting the unique ground‐shaking patterns associated with cascading ruptures in fracture networks. Simulated ground motions from cascading ruptures show consistently higher levels of PGA, PGV, and PSA across azimuths, with PGA nearly double that of non‐cascading scenarios. Our study emphasizes how fault‐zone complexities profoundly affect rupture dynamics and seismic wave radiation, and hence ground‐motion amplitude and distribution. High‐performance computing enables detailed modeling of these effects, emphasizing the importance of incorporating geometrical complexity into physics‐based seismic hazard assessment.

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| (last accessed February 2025). Patched meshing software PUMGen can be cloned from github branch PUMGen FaceIdentification64bit (https://github.com/palgunadi1993/PUMGen/tree/PUMGenFaceIdentification64bit, last accessed February 2025). The strong motion data from three central Italy’s earthquakes can be found in https://esm-db.eu/ (last accessed January 2023) (Luzi et al., 2020). The AXITRA software can be downloaded in https://github.com/coutanto/axitra/tree/master (last accessed February 2025). Various ground‐motion models are extracted from https://github.com/bakerjw/GMMs (last accessed February 2025). The supplemental material contains eight additional figures mentioned in the article. All input and mesh files are available in the Zenodo repository at Palgunadi, Gabriel, Garagash, Ulrich, Schliwa, and Mai (2024).

The authors acknowledge that there are no conflicts of interest recorded.

The authors thank the computational earthquake seismology (CES) group at KAUST for fruitful discussions and suggestions. The authors thank SeisSol’s core developers (see www.seissol.org, last accessed February 2025). The authors thank Deputy Editor‐in‐Chief Junghyun Park, Associate Editor Adrien Oth, and two reviewers for their thorough comments and suggestions. Computing resources were provided by KAUST, Thuwal, Saudi Arabia (KAUST, Project k1587, k1488, and k1343 on Shaheen II; and Project k10043 and k10044 on Shaheen‐3). The work presented in this article was supported by KAUST Grants (FRacture Activation in Geo‐reservoir‐physics of induced Earthquakes in complex fault Networks [FRAGEN], URF/1/3389‐01‐01, and BAS/1339‐01‐01). A.‐A. G. acknowledges support by Horizon Europe (ChEESE‐2P Grant Number 101093038, DT‐GEO Grant Number 101058129, and Geo‐INQUIRE Grant Number 101058518), the National Science Foundation (Grants Number EAR‐2121568, EAR‐2121568, OAC‐2139536, and OAC‐2311208), the National Aeronautics and Space Administration (80NSSC20K0495), and the Southern California Earthquake Center (SCEC Awards 22135, 23121). T. U. and A.‐A. G. acknowledge support from the Bavarian State Ministry for Science and Art in the framework of the project Geothermal‐Alliance Bavaria. D. I. G. 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).