<figure><img src="image_1.png" /></figure>

HAL
open science

# Mechanisms and Seismological Signatures of Rupture Complexity Induced by Fault Damage Zones in Fully-Dynamic Earthquake Cycle Models

Joseph Flores Cuba, Elif Oral, Benjamin Idini, Chao Liang, Jean-Paul Ampuero

<figure><img src="image_2.png" /></figure>

## To cite this version:

Joseph Flores Cuba, Elif Oral, Benjamin Idini, Chao Liang, Jean-Paul Ampuero. Mechanisms and Seismological Signatures of Rupture Complexity Induced by Fault Damage Zones in Fully-Dynamic Earthquake Cycle Models. Geophysical Research Letters, 2024, 51 (11), pp.e2024GL108792. 10.1029/2024GL108792 . insu-04726373

HAL Id: insu-04726373

https://insu.hal.science/insu-04726373v1

Submitted on 8 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geophysical Research Letters®

## RESEARCH LETTER

10.1029/2024GL108792

**Key Points:**

* Reduction of nucleation size and pulse-crack transitions are two distinct damage zone effects that induce back-propagating rupture fronts

* Damage effects can enhance high-frequency radiation and complexity of source time functions, potentially observable in the far field

* Back-propagating fronts have potential signatures in near-field seismograms and can affect peak ground motions

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

E. Oral and C. Liang,
elifo@caltech.edu;
chao.liang@scu.edu.cn

**Citation:**

Flores-Cuba, J., Oral, E., Idini, B., Liang, C., & Ampuero, J. P. (2024). Mechanisms and seismological signatures of rupture complexity induced by fault damage zones in fully-dynamic earthquake cycle models. *Geophysical Research Letters*, *51*, e2024GL108792. https://doi.org/10.1029/2024GL108792

Received 10 FEB 2024

Accepted 14 MAY 2024

© 2024. The Author(s).
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_4.png" /></figure>

# Mechanisms and Seismological Signatures of Rupture Complexity Induced by Fault Damage Zones in Fully-Dynamic Earthquake Cycle Models

J. Flores-Cuba¹,², E. Oral¹,³ ID, B. Idini³ ID, C. Liang¹,⁴ ID, and J. P. Ampuero¹ ID

¹Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Valbonne, France, ²Institut des Sciences de la Terre de Paris, Sorbonne Université, Paris, France, ³California Institute of Technology, Pasadena, CA, USA, ⁴Institute for Disaster Management and Reconstruction, Sichuan University, Chengdu, China

**Abstract** Damage zones are common around faults, but their effects on earthquake mechanics are still incompletely understood. Here, we investigate how damage affects rupture patterns, source time functions (STF) and ground motions in 2D fully-dynamic cycle models. We find that back-propagating rupture fronts emerge in large faults and can be triggered by residual stresses left by previous ruptures or by damage-induced pulse-to-crack transitions. Damage-induced back-propagating fronts are modulated by slip rate oscillations, amplify high-frequency radiation, and sharpen the multiple peaks in STF even in the absence of frictional heterogeneity or fault segmentation. Near-field ground motion is predominantly controlled by stress heterogeneity left by prior seismicity, and further amplified within the damage zone by trapped waves and outside it by secondary rupture fronts. This study refines our knowledge on damage zone effects on earthquake rupture and identifies their potentially observable signatures in the near and far field.

**Plain Language Summary** Faults are surrounded by layers of fractured rocks, known as damage zones, which can affect earthquakes and related hazards, but in ways that are still not well understood. Here, by running computer simulations, we investigate how damage zones influence earthquake ruptures and consequent ground motions. Our models fully account for seismic wave effects, produce multiple earthquake cycles, and span a large range of fault lengths and damage zone properties that are representative of natural faults. We identify characteristic patterns of earthquake rupture produced by damage zones: back-propagating fronts that re-rupture the fault, and oscillatory fault motions that affect ground shaking amplitude and frequency content. We identify which of these effects might be observable in seismograms recorded near and far from the fault. Overall, our computational study highlights significant effects of damage zones on earthquakes and on the shaking they produce. These results can guide us to better interpret earthquake source and ground motion observations, and to predict the potential characteristics of future events.

## 1. Introduction

Faults are usually surrounded by damage zones which, as increasingly demonstrated in numerical and observational studies, can substantially affect earthquake rupture processes. Fault damage zones are characterized in geological surveys by distributed fractures and micro-cracks (e.g., Mitchell & Faulkner, 2009; Savage & Brodsky, 2011), and in geophysical studies by compliant or low-velocity fault zones (e.g., Huang & Ampuero, 2011; Yang, 2015). Geological, geodetic and seismological data have been used to image damage zones (e.g., Cochran et al., 2009; Fialko, 2004; Huang et al., 2014; Lewis & Ben-Zion, 2010, and references cited therein). For example, damaged, or low-velocity, fault zones in California have thicknesses that range from several hundred meters (e.g., 200 m in the San Andreas fault) to more than a kilometer (e.g., ~1.5 km in the Calico fault). Analyzing optical satellite images of the 1992 M 7.3 Landers earthquake, Milliner et al. (2015) captured the generation of such off-fault damage, and pointed to a correlation between the width of damage zones and slip variability. A correlation between along-strike variation of fault maturity and slip is also supported by analysis of a worldwide database (Perrin et al., 2016). Very recently, co-seismic and pre-seismic damage zones with notable along-strike and off-fault variability were identified on the fault system that hosted the 2019 M 7.1 Ridgecrest earthquake (e.g., Qiu et al., 2021).

Previous modeling studies show that in the presence of damage zones, fault zone reflected waves, head waves and trapped waves can interact with the rupture and promote a number of source phenomena: pulse-like
rupture, premature rupture arrest, periodic modulation of slip rate, periodic patterns of off-fault damage, transition to supershear rupture at relatively low background stress, and rupture speeds that are theoretically unexpected for steady ruptures in homogeneous media (Harris & Day, 1997; Huang & Ampuero, 2011; Huang et al., 2014, 2016; Pelties et al., 2014). Some of these predicted effects of damage on earthquake rupture have been increasingly supported by seismological and geological observations. For example, evidence for unexpectedly fast rupture was found in earthquakes occurring within damage zones in Big Bear, Southern California (Huang et al., 2016). A faster rupture in the direction of increasing fault maturity (Perrin et al., 2016) and the sustained “slow supershear” of the 2018 Indonesia earthquake at a speed between S-wave and Eshelby's speed (Bao et al., 2019; Oral et al., 2020) have been also attributed to damage effects. Modeling studies also identify damage-induced rupture features that persist across multiple earthquake cycles, in particular back-propagating rupture fronts (Abdelmeguid et al., 2019; Idini & Ampuero, 2020; Nie & Barbot, 2022; Thakur et al., 2020) that resemble rupture patterns observed in real earthquakes (e.g., Beroza & Spudich, 1988; Hicks et al., 2020; Vallée et al., 2023).

However, our understanding of damage zone effects on earthquakes is still incomplete, partly due to limitations of previous modeling studies. Studies based on single-rupture simulations (e.g., Harris & Day, 1997; Huang et al., 2014, 2016; Oral et al., 2020) strongly depend on initial stresses that are prescribed arbitrarily. This limitation is addressed by seismic cycle modeling, in which the initial fault stresses for each earthquake result from the previous seismic and aseismic slip on the fault. To keep the computational cost affordable, the most systematic earthquake cycle studies on damaged faults (Idini & Ampuero, 2020; Nie & Barbot, 2022) adopt the quasi-dynamic approximation, in which seismic wave effects are only crudely modeled. Such dynamic effects are known to be important (Thomas et al., 2014), especially in presence of damage zones, as highlighted in recent fully-dynamic cycle models (Abdelmeguid et al., 2019; Thakur et al., 2020). On the other hand, due to their high computational cost, fully-dynamic cycle studies have explored a limited range of model parameters (e.g., Kaneko et al., 2011). In particular, the ratio of fault length to nucleation size has not yet been pushed to the high values required in continuum fault models with homogeneous friction properties to produce realistic statistics of seismicity (Barbot, 2019; Cattania, 2019) and to promote damage-induced rupture complexity (Idini & Ampuero, 2020).

Here we investigate the effects of damage zones on rupture patterns in 2D fully-dynamic earthquake cycle simulations that span a broad range of parameter values, representative of natural fault zone properties and fault lengths. To efficiently explore the fully-dynamic models, we select model parameters based on the insights from previous quasi-dynamic modeling (Idini & Ampuero, 2020). In the following, we first present our model assumptions and simulation settings. Next, we analyze the emergence of back-propagating fronts in large faults with and without damage, and the potential signatures of damage effects in near- and far-field ground motions.

## 2. Model

We consider a fault bisecting a damage zone embedded in a homogeneous elastic medium. We focus on a 2-D anti-plane problem, which corresponds to a vertical section across a strike-slip fault (Figure 1). The damage zone has a thickness $2h$ and a damage level $\Delta$, defined as the relative contrast of shear modulus between damaged ($\mu_d$) and intact rocks ($\mu_h$): $\Delta = 1 - \mu_d/\mu_h$. In terms of S-wave speeds of damaged ($V_d$) and host rocks ($V_h$), the damage level is $\Delta = 1 - (V_d/V_h)^2$.

The fault shear strength is governed by the conventional rate-and-state friction law with state evolution following the aging law (Dieterich, 1979; Ruina, 1983):

$$f(V, \theta) = f_0 + a \ln\left(\frac{V}{V_0}\right) + b \ln\left(\frac{V_0 \theta}{D_c}\right) \quad (1)$$

$$\dot{\theta} = 1 - \frac{V\theta}{D_c} \qquad (2)$$

where $V$ is slip velocity, $\theta$ the state variable, $D_c$ the characteristic slip distance of state evolution, $f_0$ the steady-state friction coefficient at the reference velocity $V_0$, and $a$ and $b$ the coefficients quantifying the direct and evolution effects, respectively.
<figure><img src="image_5.png" /><figcaption>**Figure 1.** Illustration of fault and damage zone model. (a) Conceptual, 3D representation of a vertical strike-slip fault with a damage zone (low velocity fault zone, LVFZ). (b) 2D model built based on the vertical cross-section in (a). The fault comprises a central velocity-weakening (VW) segment that hosts earthquakes, surrounded by two velocity-strengthening (VS) segments that host transient aseismic slip, in turn surrounded by outer segments that slip aseismically at steady plate velocity $V_{pt}$.</figcaption></figure>

The model comprises spatially variable frictional parameters that represent a seismogenic zone surrounded by creeping segments (Figure 1b). The central segment of length $L_{vw}$, referred to as “the fault” hereafter, is seismogenic: its friction is velocity-weakening at steady state ($a - b < 0$). It is surrounded by two segments of length $L_{vs} = L_{vw}/2$ that are velocity-strengthening ($a - b > 0$) and host transient aseismic slip. The outermost segments slip aseismically and steadily at the prescribed plate velocity $V_{pt}$, which provides the tectonic loading.

A characteristic length scale of the problem is the nucleation size, $L_{nuc}$, which is the size of the area of aseismic slip that precedes dynamic rupture. For the aging law and $a/b > 0.5$, a range of $a/b$ values typically observed in laboratory experiments, in a homogeneous medium with shear modulus $\mu$ (Rubin & Ampuero, 2005):

$$L_{nuc} = \frac{2\mu D_c b}{\pi \sigma (b - a)^2} \quad (3)$$

here $\sigma$ denotes the effective normal stress, which is uniform and constant in our models. A theoretical estimate of the nucleation size in a damage zone was derived and validated numerically by Kaneko et al. (2011). It depends on damage zone thickness $h$ and damage level $\Delta$, and ranges between the values given by Equation 3 with $\mu = \mu_h$ for small $h$ and with $\mu = \mu_d$ for large $h$. Here, we evaluate $L_{nuc}$ as Kaneko et al. (2011), and normalize distances by $L_{nuc}$ and time by $L_{nuc}/V_s$ ($V_s$ standing for S-wave speed).

The problem primarily depends on four non-dimensional parameters: damage thickness to fault length ratio ($2h/L_{vw}$), damage level ($\Delta$), fault length to nucleation size ratio ($L_{vw}/L_{nuc}$), and the ratio between direct and evolution effects of rate-and-state friction ($a/b$). We consider values of damage thickness and level within ranges that led to distinctive rupture patterns in previous work by Idini and Ampuero (2020). We vary $\Delta$ between 30% and 90%, which corresponds to a velocity reduction between 17% and 68%, similar to the range observed in nature (Huang et al., 2014). We set values of damage thickness down to $2h/L_{vw} = 1/40$. Large values of $L_{vw}/L_{nuc}$ are found necessary to produce seismicity with a realistic distribution of magnitudes (Barbot, 2019; Cattania, 2019), as mentioned earlier. We thus consider fault lengths as large as possible, while computationally affordable, up to $L_{vw}/L_{nuc} = 15$ for damage models and $L_{vw}/L_{nuc} = 40$ for homogeneous cases. The ratio $a/b$ is also known to affect the complexity of slip. For instance, when it approaches 1, the range of conditions allowing for slow slip under the aging law is wider (Rubin, 2008) and slip patterns during nucleation are more complex (Viesca, 2016). A main effect of $a/b$ is via its effect on the nucleation size (see Equation 3). Here, to focus on fault zone effects, we fixed $a/b = 0.74$, within the range of typical values in laboratory experiments. In Table S1 of Supporting Information S1, we provide all parameter values.
We use the spectral element method for 2D fully-dynamic earthquake cycle simulations. The method of Kaneko et al. (2011) was implemented by Liang et al. (2022) in the software SEM2DPACK with further optimizations and parallelism (see Data Availability Statement). It handles alternating time periods of quasi-static and dynamic fault slip by adaptive time stepping. In dynamic periods, the bottom and side boundaries function as absorbing boundaries. In quasi-static periods, we prescribe on these boundaries displacements that are consistent with the plate velocity, using a back slip approach. The switch criterion is the same as Kaneko et al. (2011) for which we verified numerical stability: the threshold velocity to switch from quasi-static to dynamic regime is 5, and 2 mm/s from dynamic to quasi-static regime. To ensure resolution of the process zone, we also set the maximum slip rate to 10 m/s after Andrews (2005). The simulations reproduce the fundamental phases of earthquake cycles: interseismic, pre-seismic, co-seismic and post-seismic slip. Here we focus on the co-seismic phases.

## 3. Results

### 3.1. Back-Propagating Fronts in Large Faults With and Without Damage

Back-propagating fronts are one notable form of rupture complexity associated with damage zone effects (Idini & Ampuero, 2020). They are secondary rupture fronts that propagate in the opposite direction to the main rupture front. Their possible existence in real faults was first suggested in a finite source inversion study of the 1984 Morgan Hills earthquake (Beroza & Spudich, 1988). Since then, to mitigate the non-uniqueness or ill-posedness of the inverse problem, most finite source inversions have adopted source parameterizations restricted to a single rupture front, which limits the possible discovery of more back-propagating fronts. More recently, with the advent of teleseismic back-projection studies and more flexible source inversion approaches, back-propagating fronts have been robustly imaged on different events, including the 2010 El Mayor-Cucapah (Meng et al., 2011), the 2016 Romanche oceanic transform fault (Hicks et al., 2020), and the 2019 intermediate-depth northern Peru earthquakes (Vallée et al., 2020). Numerical studies point to damage effects (Idini & Ampuero, 2020) and fault size effects (Barbot, 2019) as possible origins of back-propagating fronts. In the following, we distinguish these two effects.

Regardless of the presence of damage, we find that stress concentrations near the edges of creeping sections or of previous partial ruptures in a large fault can generate back-propagating fronts. Faults that are much longer than their nucleation length generate seismicity with a wide range of rupture sizes (Cattania, 2019), resulting in a heterogeneous stress state prior to any large rupture. In our simulations, such stress heterogeneity emerges when $L_{vw}/L_{nuc} \ge 10$. In smaller faults, as those considered by Kaneko et al. (2011), all events break the entire fault and leave a relatively smooth state of stress. Figures S1 and S2 in Supporting Information S1 show the fault stresses before and after a full rupture, and the spatiotemporal distribution of slip rate in models without damage zone, for $L_{vw}/L_{nuc} = 10$ and $L_{vw}/L_{nuc} = 40$. In both cases, rupture nucleates near the bottom edge and propagates bilaterally at average speeds of 50%–80% of the S-wave speed $V_s$. Near the stress concentrations in either side of the fault, upon the arrival of main rupture fronts, new fronts emerge and propagate in the opposite direction at speeds near $V_s$. The emergence of such secondary fronts in the absence of damage supports the previous findings of Figure 9a in Barbot (2019), Figure 1e in Cattania (2019), and Figure 2d in Idini and Ampuero (2020).

The initiation of back-propagating fronts at residual stresses occurs also in the presence of damage. Indeed, damage favors this mechanism by reducing the nucleation size. The example in Figure 2a shows two such back-propagating fronts nucleating near the peaks of initial stress. They are modulated by damage zone effects: they interact with fault zone trapped waves and break up into multiple pulses, as further discussed in Section 3.2.

The presence of damage produces a separate driving mechanism for back-propagating fronts, related to transitions between pulse-like and crack-like rupture behavior. This mechanism, first identified by Idini and Ampuero (2020) in quasi-dynamic models, can be summarized as follows. In homogeneous media, our models produce crack-like ruptures. Also in a damage zone, rupture is crack-like initially: since its size is much shorter than the damage zone thickness $2h$, being far from the host rock, it behaves as in a —uniformly damaged— homogeneous medium. When its length exceeds $2h$ in a highly-damaged zone, however, the rupture becomes pulse-like as it would in an elastic slab of thickness $2h$ with rigid boundaries (Field & Baker, 1962). As rupture grows much larger than $2h$, it starts losing its sensitivity to the damage zone and behaves as in a homogeneous intact medium. Therefore, the pulse front becomes crack-like again. For this new crack propagates bilaterally, two crack fronts emerge from the
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Back-propagating fronts in damaged faults. (a) Simulation results in the damage model with $L_{vw}/L_{nuc} = 5$, $\Delta = 60\%$ and $L_{vvi}/2h = 40$. Left: initial (black) and final (red) stresses along the fault, $\tau$, normalized as $\tau/\sigma - f_0$. Right: spatiotemporal evolution of slip rate. Panel (b) same as panel (a) but for the damage model with $L_{vvi}/L_{nuc} = 15$, $\Delta = 90\%$, and $L_{vvi}/2h = 40$. (c) Slip rate at the position $2.4 L_{nuc}$ in the damage model shown in (a). (d) Spectrum of a time window of slip rate, shown in (c) by a horizontal gray line, containing fault-zone-induced oscillations, normalized by the spectral amplitude at lowest frequency.</figcaption></figure>

pulse front (Figure 2b, detailed below). One of these secondary crack fronts propagates in the opposite direction to the pulse: a back-propagating front. As they keep growing, the new cracks undergo crack-to-pulse and pulse-to-crack transitions; and the process repeating successively leads to the formation of multiple secondary fronts. A slip budget argument further explains the necessity of multiple fronts: the slip produced by cracks and damage-induced pulses is largely different, because the former scales with rupture length whereas the latter scales with $h$. Even though ruptures much larger than $h$ eventually become cracks, they pass through a stage of pulse-like rupture which leaves a slip deficit. A back-propagating front makes up for this deficit, but only partially, because it also eventually turns into a pulse. To completely fill the slip gap thus requires multiple secondary fronts.
<figure><img src="image_7.png" /><figcaption>**Figure 3.** Far-field view of damage-induced source complexity. Source time functions (STF) of models with damage zone ($L_{vw}/L_{nuc} = 15$, $\Delta = 90\%$, and $L_{vw}/2h = 40$) and without damage zone ($L_{vw}/L_{nuc} = 40$) (left), and their spectra (right). To facilitate the comparison, STF are normalized by their peak value, time by rupture duration, spectral amplitudes by seismic moment, and frequency by the inverse of rupture duration. The damage model is the same as in Figure 2b.</figcaption></figure>

Our fully-dynamic simulations confirm the existence of the damage-induced mechanism of back-propagating fronts in large faults. Because a sufficiently large $L_{vw}/2h$ and high damage are required for the crack-to-pulse transition to manifest, we set $L_{vw}/2h = 40$ and $\Delta = 90\%$. An example is shown in Figure 2b. Initially, rupture propagates bilaterally at speeds in the range of 50%-100% $V_d$. Secondary fronts nucleate at various locations along the fault, including multiple fronts that nucleate well after the passage of the main front and propagate bilaterally at speeds close to $V_d$. Their nucleation points do not coincide with the peaks of initial stress, but rather with stress heterogeneities forming during the previous stages of the rupture. These rich rupture patterns do not occur in models without damage, even with large $L_{vw}/L_{nuc}$ (see the homogeneous case in Figure S2 of Supporting Information S1), which thus counters the suggestion of Nie and Barbot (2022) that rupture style in a damage zone is simply controlled by the ratio of fault size to nucleation size. Moreover, by comparing quasi-dynamic and fully-dynamic models, we find that dynamic effects tend to increase the occurrence of secondary fronts and amplify their peak slip rates (Figures S3 and S4 in Supporting Information S1).

## 3.2. Source Modulations Caused by Fault Damage Zones

We find that interactions between rupture fronts and trapped waves in a damage zone cause slip rate oscillations, at frequencies that are characteristic of the damage zone. Figure 2c shows an example of such oscillations. Their spectrum prominently peaks at a frequency near $V_d/4(2h)$, the fundamental frequency of wave reverberations across the damage zone which constructively interfere to form trapped waves. Our analyses of cases with different damage levels confirm this interpretation (Figure S5 in Supporting Information S1).

Damage-induced rupture effects can sharpen the complexity of source time functions (STF). Figure 3a compares STFs of models with and without damage zone. While the STF in the homogeneous case has a single peak, the damage model produces multiple sharp peaks resulting from both back-propagating fronts and slip rate oscillations (as also found in other cases, Figure S6 in Supporting Information S1). Real STFs often exhibit multiple peaks that are usually interpreted as sub-events originating from different rupture segments, often associated with structural segmentation by frictional or geometrical barriers along the fault (e.g., Ross et al., 2019; Vallée, 2013). Our finding alternatively suggests that multiple peaks in a STF can originate from damage zone effects, even on faults with uniform frictional properties and simple geometry.
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Stress heterogeneity and damage effects on near field ground motion. (a) Initial stress and spatiotemporal change of slip rate in the homogeneous model of $L_{vw}/L_{nuc} = 10$, panel (b) same as panel (a) but for the damage model with $L_{vw}/L_{nuc} = 1.2$, $\Delta = 90\%$ and $L_{vw}/2h = 40$. (c) Comparison of peak ground velocities (PGV) as a function of fault distance between homogeneous and damage models. Stress drop and PGV are normalized by $\sigma$ and $V_s\sigma_{f0}/\mu$, respectively.</figcaption></figure>

Damage-induced rupture complexity also amplifies high-frequency radiation. Comparing STF source spectra of models with and without damage zone (Figure 3b) reveals a systematic amplification above the corner frequency in damage models relative to homogeneous models, up to a factor of ~10. This highlights the potential significance of damage effects on far-field observations. The damage-caused excessive high-frequency radiation occurs in the broad band above the corner frequency—not at a specific frequency that can be associated with damage zone properties. While this challenges the inference of damage-induced rupture effects from far-field data, in the next section, we investigate the potential signatures of damage-induced slip rate oscillations and back-propagating fronts in near-field observations.

### 3.3. Damage Effects on Near-Field Ground Motions

The effects of both initial stresses and damage zones on earthquake rupture affect ground motion and its spatial variability. A higher initial stress can result in stronger ground motion by increasing stress drop (e.g., Cotton et al., 2013) and rupture speed (e.g., Aagaard & Heaton, 2004), and initial stress heterogeneity can enhance high-frequency strong motion (Kame & Uchida, 2008; Madariaga, 1983). Damage can amplify near-source ground motion by trapped wave modulation of the source (Section 3.2, and e.g., Ben-Zion et al., 2003). Our models produce ground motion amplification by both factors, and here we assess their respective effects on peak ground velocities (PGV).

We find that near-field ground motion is governed by initial stress heterogeneity, and further affected by damage effects. Figures 4a and 4b shows the initial stresses and spatial distribution of PGVs for two homogeneous and damage models. In the homogeneous model, rupture nucleates near the bottom edge where the initial stress is the largest. Initial stress peaks are also present near the upper edge. In the damage model, the largest stresses are concentrated near both edges, and residual stress peaks are also present in the central portion. The largest PGVs (above 2 and 4 m/s in homogeneous and damage models, respectively) are concentrated near these high-stress areas. This spatial correlation is expected from the radiation of strong motion phases due to abrupt changes in rupture speed when rupture encounters residual stress concentrations (Kame & Uchida, 2008; Madariaga, 1983). Such ground motion amplification due to initial stress heterogeneity manifests as along-fault ground motion variability. Comparing the two cases in Figure 4c, PGV decreases with distance to the fault as expected. The PGVs in the damage model are smaller than in the homogeneous model outside the damage zone, but larger inside the damage zone, by a factor of ~2. This damage-induced amplification results in a sharp contrast between the regions inside and outside the damage zone all along the fault. In Figure S7 in Sup-

porting Information S1 we show similar findings for a case with smaller nucleation length. Overall, we find that initial stresses are the main control of the spatial variability of peak ground motion along the fault, while damage-induced amplification strongly affects the ground motion variability across the fault.

Back-propagating fronts are visible in the near-field seismograms and can locally affect the peak ground motion. In both homogeneous and damage cases, waves radiated by secondary fronts are present in the seismograms as later arriving pulses at various distances (Figure S8 in Supporting Information S1). At some distances (Figure S9b in Supporting Information S1) the largest peaks are in the first arriving pulses, which are radiated by the primary rupture front. At other distances (Figures S9c and S9d in Supporting Information S1), the later pulses generated by secondary fronts have the largest amplitudes. Considering the recent advances in near-fault monitoring techniques (e.g., Li et al., 2023; Qiu et al., 2021), our finding underpins the potential for the discovery of back-propagating fronts in dense arrays close to faults.
4. Conclusions and Discussion

We studied the effects of fault damage zones on rupture dynamics and ground motions by 2D fully-dynamic earthquake cycle modeling. Our simulations span a relevant range of fault sizes (relative to nucleation size, *L*<sub>nuc</sub>) and damage zone properties, and expand the insights from previous quasi-dynamic modeling studies.

We confirm that both damage zone properties and relative fault size control rupture complexity, and we identify their respective effects. In particular, we distinguish the mechanisms of secondary rupture front generation due to each. On large faults, regardless of the presence of damage, the emergence of heterogeneous stress states featuring residual stress concentrations induces back-propagating fronts. In the presence of damage, an additional mechanism owing to a pulse-to-crack transition (Idini & Ampuero, 2020) operates on faults that have sufficiently high damage levels and thicknesses, and are relatively large (15*L*<sub>nuc</sub> here).

Damage-induced rupture complexity potentially imprints seismological signatures both in the near and far field. Rupture fronts interact with damage zone trapped waves, leading to oscillations in slip rate at resonance frequencies that are characteristic of the damage zone. Damage-induced oscillations and secondary fronts amplify high-frequency radiation and enhance the complexity of STF, manifested by multiple moment rate peaks, which is potentially observable in the far field. Regarding near-field ground motions, residual stress concentrations predominantly shape the spatial variability of PGV along strike, while damage affects the variability across the fault by introducing a contrast between ground motions inside and outside the damage zone. Additionally, damage-induced secondary fronts can locally amplify peak ground motions far from the damage zone, and increase the hazard therein.

Exploring the (possibly competing) effects of other relevant processes on rupture dynamics constitutes a major prospect of our study. To facilitate our parametric investigation and to better isolate fault zone effects, we considered a basic empirical friction law without including realistic complexities of rock behavior as a function of depth, temperature and slip rates (Barbot, 2023, and references therein). Moreover, we modeled damage zones as compliant zones of constant thickness without accounting for geometrical irregularities and roughness (e.g., Cappa et al., 2014; Sun & Cattania, 2024), off-fault plasticity (e.g., Mia et al., 2022, 2023), and time-dependent damage evolution (e.g., Thakur et al., 2020). In particular, the properties of natural damage zones vary systematically along fault dip and strike, which can result in systematic variations of seismicity and source properties (Cappa et al., 2014; Ito & Kaneko, 2023; Perrin et al., 2016). While our idealization of uniform damage geometry and properties allowed for a systematic identification of damage effects, the role of damage segmentation across multiple cycles can be readily studied with our current computational tools and deserves future research. We also note the need for additional experimental, observational, and numerical studies on time-dependent evolution of damage and friction properties to advance our models.

Data Availability Statement

All data needed to reproduce this work is available online: 2D fully dynamic cycle modeling tools and quasi-dynamic comparison models can be found at Ampuero (2012), Luo et al. (2017) and Oral (2024), respectively.

References

Aagaard, B. T., & Heaton, T. H. (2004). Near-source ground motions from simulations of sustained intersonic and supersonic fault ruptures. *Bulletin of the Seismological Society of America*, 94(6), 2064–2078. https://doi.org/10.1785/0120030249

Abdelmeguid, M., Ma, X., & Elbanna, A. (2019). A novel hybrid finite element-spectral boundary integral scheme for modeling earthquake cycles: Application to rate and state faults with low-velocity zones. *Journal of Geophysical Research*, 124(12), 12854–12881. https://doi.org/10.1029/2019jb018036

Ampuero, J.-P. (2012). SEM2DPACK, a spectral element software for 2D seismic wave propagation and earthquake source dynamics, v2.3.8. [Software]. Zenodo. https://doi.org/10.5281/zenodo.230363

Andrews, D. (2005). Rupture dynamics with energy loss outside the slip zone. *Journal of Geophysical Research*, 110(B1), B01307. https://doi.org/10.1029/2004JB003191

Bao, H., Ampuero, J.-P., Meng, L., Fielding, E. J., Liang, C., Milliner, C. W. D., et al. (2019). Early and persistent supershear rupture of the 2018 magnitude 7.5 Palu Earthquake. *Nature Geoscience*, 12(3), 200–205. https://doi.org/10.1038/s41561-018-0297-z

Barbot, S. (2019). Slow-slip, slow earthquakes, period-two cycles, full and partial ruptures, and deterministic chaos in a single asperity fault. *Tectonophysics*, 768, 228171. https://doi.org/10.1016/j.tecto.2019.228171

Barbot, S. (2023). Constitutive behavior of rocks during the seismic cycle. *AGU Advances*, 4(5), e2023AV000972. https://doi.org/10.1029/2023av000972

Acknowledgments

This work has been supported by the French National Research Agency (ANR) through project FAULTS_R_GEMS (Grant ANR-17-CE31-0008) and Investments-in-the-Future project UCAJEDI (Grant ANR-15-IDEX-01), and the UCAJEDI Académies 2 and 3 through the project PERFAULT-3M "Physics of earthquake rupture and fault growth: multi-scale modeling of material failure". It was also supported by the National Natural Science Foundation of China (NSFC) through the early career research grant (Grant 42204059) and the Fundamental Research Funds for Central Universities disseminated by IDMR at Sichuan University. We thank Prithvi Thakur for providing an initial version of the sparse-matrix optimization, Caroline Ramel for her support in the use of the cluster in Geoazur, and Hojjat Kaveh for his assistance on using QDYN with Python. We also thank Isabelle Manighetti, and acknowledge the Resnick High Performance Computing Center of Caltech when finalizing the simulations. We are also grateful to the Associate Editor, Ahmed Elbanna, and an anonymous reviewer for their constructive revision that improved the quality of this manuscript.
Ben-Zion, Y., Peng, Z., Okaya, D., Seeber, L., Armbruster, J. G., Ozer, N., et al. (2003). A shallow fault-zone structure illuminated by trapped waves in the Karadere-Duzce branch of the North Anatolian fault, western Turkey. *Geophysical Journal International*, **152**(3), 699–717.
https://doi.org/10.1046/j.1365-246x.2003.01870.x

Beroza, G. C., & Spudich, P. (1988). Linearized inversion for fault rupture behavior: Application to the 1984 Morgan Hill, California, earthquake. *Journal of Geophysical Research*, **93**(B6), 6275–6296. https://doi.org/10.1029/jb093ib06p06275

Cappa, F., Perrin, C., Manighetti, I., & Delor, E. (2014). Off-fault long-term damage: A condition to account for generic, triangular earthquake slip profiles. *Geochemistry, Geophysics, Geosystems*, **15**(4), 1476–1493. https://doi.org/10.1002/2013gc005182

Cattania, C. (2019). Complex earthquake sequences on simple faults. *Geophysical Research Letters*, **46**(17–18), 10384–10393. https://doi.org/10.1029/2019gl083628

Cochran, E. S., Li, Y.-G., Shearer, P. M., Barbot, S., Fialko, Y., & Vidale, J. E. (2009). Seismic and geodetic evidence for extensive, long-lived fault damage zones. *Geology*, **37**(4), 315–318. https://doi.org/10.1130/G25306A.1

Cotton, F., Archuleta, R., & Causse, M. (2013). What is sigma of the stress drop? *Seismological Research Letters*, **84**(1), 42–48. https://doi.org/10.1785/0220120087

Dieterich, J. H. (1979). Modeling of rock friction: 1. Experimental results and constitutive equations. *Journal of Geophysical Research*, **84**(B5), 2161–2168. https://doi.org/10.1029/jb084ib05p02161

Fialko, Y. (2004). Temperature fields generated by the Elastodynamic propagation of shear cracks in the Earth. *Journal of Geophysical Research*, **109**(B1), B01303. https://doi.org/10.1029/2003JB002497

Field, F., & Baker, B. (1962). Crack propagation under shear displacements.

Harris, R. A., & Day, S. M. (1997). Effects of a low-velocity zone on a dynamic rupture. *Bulletin of the Seismological Society of America*, **87**(5), 1267–1280. https://doi.org/10.1785/bssa0870051267

Hicks, S. P., Okuwaki, R., Steinberg, A., Rychert, C. A., Harmon, N., Abercrombie, R. E., et al. (2020). Back-propagating supershear rupture in the 2016 m w 7.1 Romanche transform fault Earthquake. *Nature Geoscience*, **13**(9), 647–653. https://doi.org/10.1038/s41561-020-0619-9

Huang, Y., & Ampuero, J.-P. (2011). Pulse-like ruptures induced by low-velocity fault zones. *Journal of Geophysical Research*, **116**(B12), B12307. https://doi.org/10.1029/2011jb008684

Huang, Y., Ampuero, J.-P., & Helmberger, D. V. (2014). Earthquake ruptures modulated by waves in damaged fault zones. *Journal of Geophysical Research*, **119**(4), 3133–3154. https://doi.org/10.1002/2013jb010724

Huang, Y., Ampuero, J.-P., & Helmberger, D. V. (2016). The potential for supershear Earthquakes in damaged fault zones – Theory and observations. *Earth and Planetary Science Letters*, **433**, 109–115. https://doi.org/10.1016/j.epsl.2015.10.046

Idini, B., & Ampuero, J.-P. (2020). Fault-zone damage promotes pulse-like rupture and back-propagating fronts via Quasi-Static effects. *Geophysical Research Letters*, **47**(23), e2020GL090736. https://doi.org/10.1029/2020gl090736

Ito, R., & Kaneko, Y. (2023). Physical mechanism for a temporal decrease of the Gutenberg-Richter b-value prior to a large Earthquake. *Journal of Geophysical Research*, **128**(12), e2023JB027413. https://doi.org/10.1029/2023jb027413

Kame, N., & Uchida, K. (2008). Seismic radiation from dynamic coalescence, and the reconstruction of dynamic source parameters on a planar fault. *Geophysical Journal International*, **174**(2), 696–706. https://doi.org/10.1111/j.1365-246x.2008.03849.x

Kaneko, Y., Ampuero, J.-P., & Lapusta, N. (2011). Spectral-element simulations of long-term fault slip: Effect of low-rigidity layers on earthquake-cycle dynamics. *Journal of Geophysical Research*, **116**(B10), B10313. https://doi.org/10.1029/2011jb008395

Lewis, M. A., & Ben-Zion, Y. (2010). Diversity of fault zone damage and trapping structures in the Parkfield section of the San Andreas fault from comprehensive analysis of near fault seismograms. *Geophysical Journal International*, **183**(3), 1579–1595. https://doi.org/10.1111/j.1365-246x.2010.04816.x

Li, J., Kim, T., Lapusta, N., Biondi, E., & Zhan, Z. (2023). The break of earthquake asperities imaged by distributed acoustic sensing. *Nature*, **620**(7975), 800–806. https://doi.org/10.1038/s41586-023-06227-w

Liang, C., Ampuero, J.-P., & Pino Muñoz, D. (2022). The paucity of supershear earthquakes on large faults governed by rate and state friction. *Geophysical Research Letters*, **49**(22), e2022GL099749. https://doi.org/10.1029/2022gl099749

Luo, Y., Ampuero, J.-P., Galvez, P., van den Ende, M., & Idini, B. (2017). QDYN: A Quasi-DYNamic earthquake simulator (v1.1) (QDYN1.1) [Software]. Zenodo. https://doi.org/10.5281/zenodo.322459

Madariaga, R. (1983). High frequency radiation from dynamic earthquake. *Annals of Geophysics*, **1**, 17.

Meng, L., Ampuero, J. P., Page, M. T., & Hudnut, K. W. (2011). Seismological evidence and dynamic model of reverse rupture propagation during the 2010 M7.2 El Mayor Cucapah earthquake. In *Agu fall meeting abstracts* (Vol. 2011).

Mia, M. S., Abdelmeguid, M., & Elbanna, A. E. (2022). Spatio-temporal clustering of seismicity enabled by off-fault plasticity. *Geophysical Research Letters*, **49**(8), e2021GL097601. https://doi.org/10.1029/2021gl097601

Mia, M. S., Abdelmeguid, M., & Elbanna, A. E. (2023). The spectrum of fault slip in elastoplastic fault zones. *Earth and Planetary Science Letters*, **619**, 118310. https://doi.org/10.1016/j.epsl.2023.118310

Milliner, C. W., Dolan, J. F., Hollingsworth, J., Leprince, S., Ayoub, F., & Sammis, C. G. (2015). Quantifying near-field and off-fault deformation patterns of the 1992 mw 7.3 l Anders earthquake. *Geochemistry, Geophysics, Geosystems*, **16**(5), 1577–1598. https://doi.org/10.1002/2014gc005693

Mitchell, T., & Faulkner, D. (2009). The nature and origin of off-fault damage surrounding strike-slip fault zones with a wide range of displacements: A field study from the Atacama fault system, Northern Chile. *Journal of Structural Geology*, **31**(8), 802–816. https://doi.org/10.1016/j.jsg.2009.05.002

Nie, S., & Barbot, S. (2022). Rupture styles linked to recurrence patterns in seismic cycles with a compliant fault zone. *Earth and Planetary Science Letters*, **591**, 117593. https://doi.org/10.1016/j.epsl.2022.117593

Oral, E. (2024). Elifo/Qdyn: Dataset for Flores-Cuba et al. [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.11130354

Oral, E., Weng, H., & Ampuero, J. P. (2020). Does a damaged-fault zone mitigate the near-field impact of supershear Earthquakes? Application to the 2018 7.5 Palu, Indonesia, earthquake. *Geophysical Research Letters*, **47**(1), e2019GL085649. https://doi.org/10.1029/2019gl085649

Pelties, C., Gabriel, A.-A., & Ampuero, J.-P. (2014). Verification of an ADER-DG method for complex dynamic rupture problems. *Geoscientific Model Development*, **7**(3), 847–866. https://doi.org/10.5194/gmd-7-847-2014

Perrin, C., Manighetti, I., Ampuero, J.-P., Cappa, F., & Gaudemer, Y. (2016). Location of largest earthquake slip and fast rupture controlled by along-strike change in fault structural maturity due to fault growth. *Journal of Geophysical Research*, **121**(5), 3666–3685. https://doi.org/10.1002/2015jb012671

Qiu, H., Ben-Zion, Y., Catchings, R., Goldman, M. R., Allam, A. A., & Steidl, J. (2021). Seismic imaging of the mw 7.1 Ridgecrest earthquake rupture zone from data recorded by dense linear arrays. *Journal of Geophysical Research: Solid Earth*, **126**(7), e2021JB022043. https://doi.org/10.1029/2021jb022043
Ross, Z. E., Idini, B., Jia, Z., Stephenson, O. L., Zhong, M., Wang, X., et al. (2019). Hierarchical interlocked orthogonal faulting in the 2019 Ridgecrest earthquake sequence. *Science*, **366**(6463), 346–351. https://doi.org/10.1126/science.aaz0109

Rubin, A. M. (2008). Episodic slow slip events and rate-and-state friction. *Journal of Geophysical Research*, **113**(B11), B11414. https://doi.org/10.1029/2008JB005642

Rubin, A. M., & Ampuero, J.-P. (2005). Earthquake nucleation on (aging) rate and state faults. *Journal of Geophysical Research*, **110**(B11), B11312. https://doi.org/10.1029/2005jb003686

Ruina, A. (1983). Slip instability and state variable friction laws. *Journal of Geophysical Research*, **88**(B12), 10359–10370. https://doi.org/10.1029/jb088ib12p10359

Savage, H. M., & Brodsky, E. E. (2011). Collateral damage: Evolution with displacement of fracture distribution and secondary fault strands in fault damage zones. *Journal of Geophysical Research*, **116**(B3), B03405. https://doi.org/10.1029/2010jb007665

Sun, Y., & Cattania, C. (2024). Propagation of slow slip events on rough faults: Clustering, back propagation, and re-rupturing.

Thakur, P., Huang, Y., & Kaneko, Y. (2020). Effects of low-velocity fault damage zones on long-term earthquake behaviors on mature strike-slip faults. *Journal of Geophysical Research: Solid Earth*, **125**(8), e2020JB019587. https://doi.org/10.1029/2020jb019587

Thomas, M. Y., Lapusta, N., Noda, H., & Avouac, J.-P. (2014). Quasi-dynamic versus fully dynamic simulations of earthquakes and aseismic slip with and without enhanced Coseismic weakening. *Journal of Geophysical Research*, **119**(3), 1986–2004. https://doi.org/10.1002/2013jb010615

Vallée, M. (2013). Source time function properties indicate a strain drop independent of earthquake depth and magnitude. *Nature Communications*, **4**(1), 2606. https://doi.org/10.1038/ncomms3606

Vallée, M., Grandin, R., Nocquet, J.-M., Villegas, J.-C., Vaca, S., Xie, Y., et al. (2020). Rupture characteristics of the 2019 North Peru Intraslab earthquake (Mw8.0). In *Egu general assembly conference abstracts*. 10429. https://doi.org/10.5194/egusphere-egu2020-10429

Vallée, M., Xie, Y., Grandin, R., Villegas-Lanza, J. C., Nocquet, J.-M., Vaca, S., et al. (2023). Self-reactivated rupture during the 2019 Mw= 8 Northern Peru Intraslab earthquake. *Earth and Planetary Science Letters*, **601**, 117886. https://doi.org/10.1016/j.epsl.2022.117886

Viesca, R. C. (2016). Self-similar slip instability on interfaces with rate-and state-dependent friction. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, **472**(2192), 20160254. https://doi.org/10.1098/rspa.2016.0254

Yang, H. (2015). Recent advances in imaging crustal fault zones: A review. *Earthquake Science*, **28**(2), 151–162. https://doi.org/10.1007/s11589-015-0114-3