**Fault-zone damage promotes pulse-like rupture and**rapid-tremor-reversals B. Idini∗1and J.-P. Ampuero1,2 1Seismological Laboratory, California Institute of Technology 2Universit´ e Cˆ ote d’Azur, IRD, CNRS, Observatoire de la Cˆ ote d’Azur, G´ eoazur This manuscript is a non-peer reviewed preprint submitted to Nature Geoscience. Corresponding author: Benjamin Idini, bidiniza@gps.caltech.edu –1–Abstract Damage zones are ubiquitous components of faults that may affect the nucleation, propagation and arrest of earthquake ruptures. Dynamic rupture simulations show that waves trapped in a fault damage zone can induce pulse-like rupture, a mode of earthquake rupture propagation that is a primary candidate for the origin of fault slip complexity and for the apparent weakness of major faults. However, the efficiency of the mechanism revealed in previous studies depends strongly on arbitrarily pre- scribed initial stresses. Here, we investigate the promotion of pulse-like rupture by damaged fault zones through numerical simulations of multiple earthquake cycles in which the distribution of initial stress before each rupture is a self-consistent result of the earthquake cycle. We consider a fault bisecting a homogeneous low-rigidity layer embedded in an intact medium. Using scaling arguments we show that pulse-like ruptures are expected to appear in a highly compliant fault zone after the rupture has grown larger than the fault zone thickness. We confirm this result by conduct- ing quasi-dynamic earthquake cycle simulations on fault zones with varying degrees of damage and thickness. Over a wide range of fault zone properties, fault-zone ef- fects produce pulse-like ruptures with shorter rise-times and flatter slip profiles than ruptures in an intact homogeneous medium. We also find complex rupture patterns involving back-propagating secondary fronts that emerge from the primary rupture front and propagate in the opposite direction. These complex slip patterns robustly persist over multiple earthquake cycles. While such patterns are challenging to resolve in current seismological observations of large earthquakes, slow-slip numerical models show similar slip complexity, suggesting a connection between a basic structural fea- ture of faults and rapid-tremor-reversals observed during episodic tremor and slip in Cascadia and Japan. 1 Introduction Pulse-like rupture is a common mode of earthquake propagation in which the duration of slip at each point of the fault, known as the rise-time, is short compared to the total rupture duration (Heaton, 1990). Rupture pulses play a prominent role in the theory of earthquake mechanics because they can radically affect the earthquake energy balance (Nielsen & Madariaga, 2003), reduce the apparent strength of faults (Noda et al., 2009), enhance the spatial heterogeneity of earthquake slip and stress (Aagaard & Heaton, 2008), and promote complexity of seismicity manifested by a broad range of event magnitudes (Cochard & Madariaga, 1996). Yet their origin is not completely established. Several mechanisms of pulse generation have been proposed, involving healing fronts emerging from features in the friction law (Cochard & Madariaga, 1996; G. Perrin et al., 1995), from earlier arrest of one dimension of rupture (Day, 1982; Johnson, 1990) or from fault heterogeneities (Beroza & Mikumo, 1996). Here, we show that the mere presence of damage zones around faults provides a simple and persistent mechanism of pulse-like rupture. Faults are usually embedded in a damaged zone (Fig. 1a) characterized in field observations by distributed fractures and micro-cracks (Mitchell & Faulkner, 2009; Savage & Brodsky, 2011) and in seismological and geodetic observations by reduced wave speeds or elastic modulus relative to the host rock (Ben-Zion et al., 2003; Li et al., 2007; Yang & Zhu, 2010). Seismic imaging methods resolve fault zones of strike- slip faults as flower-structures with depth-varying thickness and damage (Ben-Zion et al., 2003). Hereafter, we refer to these structures as low-velocity fault zones (LVFZ). Recent dynamic rupture simulations show that the presence of a LVFZ can induce complex rupture patterns: pulse-like rupture promoted by healing fronts mediated by reflected waves, oscillations of slip-rate and rupture speed, and supershear rupture at low background stress (Huang & Ampuero, 2011; Huang et al., 2014, 2016). One –2–

shortcoming of these and other numerical studies (Harris & Day, 1997) is that they are based on single-rupture simulations, which make arbitrary assumptions on the initial stresses on the fault. The assumed initial stresses have a strong effect on the resulting rupture patterns, including on the selection between pulse-like and crack-like rupture (Huang & Ampuero, 2011). Thus, the robustness of the generation of pulses by a LVFZ is still an open question. Here we address this question through multi-cycle earthquake modeling, in which the fault stress at the onset of each rupture is not prescribed arbitrarily but is a consistent result of the long-term dynamics (Erickson et al., 2019). Previous multi-cycle simulations in LVFZ models, conducted under moderate damage levels and moderate fault lengths relative to the nucleation size (Kaneko et al., 2011), showed that a LVFZ reduces the nucleation length and increases the peak slip rate, recurrence period and seismic potency, but did not report on rupture complexities such as pulse-like rupture. Here we investigate if the most dramatic effects of a LVFZ identified in single rupture simulations persist after multiple earthquake cycles when high levels of damage and large faults are considered. We follow two complementary approaches. In Section 2, we use a scaling argument to develop insight on how a LVFZ changes the kinematics of rupture propagation relative to an undamaged fault zone. In Section 3, we conduct earthquake cycle simulations including the effects of a LVFZ on the transfer of static stress. 2 A simple fault-zone model We consider a minimalistic, tabular LVFZ model defined by a finite fault of length L bisecting a homogeneous low-rigidity layer at y = 0. The low-rigidity layer represents the damage zone and is embedded in an intact medium (Fig. 1). The LVFZ is specified by its half-thickness h and damage level: ∆ = 1–µd µ, (1) where µdand µ are the shear moduli of the LVFZ and intact medium, respectively. The only non-zero displacement component points in the z direction, outward from the x − y plane, therefore slip is antiplane: ∆u(x,t) = uz(x,y = 0+,t) − uz(x,y = 0−,t) (2) The model converges to two different homogeneous end-member models, depending on the fault zone thickness. When h/L is very small, the model approaches a homogeneous intact medium with shear modulus µ. When h/L is very large, the model tends to a homogeneous damaged medium with shear modulus (1 − ∆)µ. Key effects of a LVFZ on rupture propagation are highlighted by analyzing the limiting case when ∆ → 1, in which the medium outside the damage zone is rigid. We consider a rupture growing quasi-statically with prescribed uniform stress drop ∆τ and increasing rupture length L(t). The fault-zone thickness h is fixed. The resulting slip profiles, computed by solving numerically a static problem as described in Appendix A, are shown in Fig. 1c. The shape of the slip profile is indicative of the style of rupture: crack-like ruptures show an elliptical slip profile while pulse-like ruptures have a flat slip profile. While the rupture is small (L ? 2h), it only interacts with the damaged zone and therefore has a crack-like slip profile, as in a uniformly damaged infinite medium. Its slip grows proportional to rupture length as ∆u(t) ∼ the rupture grows large (L ? 2h), it interacts with a thin elastic slab of thickness h and develops a pulse-like slip profile. Its slip reaches a value independent of rupture length, ∆u ∼ constant stress drop in a highly-damaged LVFZ will initiate as a crack-like rupture ∆τ 2µ(1−∆)L(t). As ∆τ µ(1−∆)h. Connecting these two stages together, a growing rupture with –3–

and later transition to pulse-like rupture. The transition is characterized by saturation of slip caused by the LVFZ once the rupture grows larger than 2h. The above picture of crack-to-pulse transition provides insight into what controls rise-time in a damaged fault zone in the absence of wave propagation effects, which will help us interpret the results of the quasi-dynamic earthquake cycle models presented later. The rise-time at the hypocenter is the time required for the appearance of a healing front. This time corresponds kinematically to the emergence of pulse-like rupture, which is approximately the time required for the size of the initial crack to grow up to L = 2h. Assuming a constant rupture speed vr, the size of the rupture is L ∼ vrt, hence the rise time at the hypocenter follows: t ∼2h (3) vr This estimation of rise-time is valid at other locations beyond the hypocenter assuming that the propagation speed of the healing front is close to the rupture speed. Because rise-time can be shorter away from the hypocenter (Huang & Ampuero, 2011), equation (3) should be taken as an upper bound. The resulting upper bound for the pulse width, defined as the distance between the position of the rupture front and the healing front, is: l ∼ vrt ∼ 2h The previous simplified analysis predicts the emergence of pulse-like rupture from static effects only, independently of the presence of trapped and reflected waves in the LVFZ. (4) 3 Multi-cycle numerical simulations 3.1 Methods We characterize the effect of a LVFZ on rise-time and slip profile in earthquake cycle simulations covering a wide range of values of LVFZ thickness and damage. We adopt a spectral boundary integral equation method (SBIEM) (Luo et al., 2017) and a quasi-dynamic approximation where wave-related effects are crudely represented by a radiation-damping term (Rice, 1993). The modeling approach captures the static LVFZ effects described earlier, and its computational efficiency enables a comprehen- sive exploration of the parameter space. The fault strength is prescribed by the Dieterich-Ruina rate-and-state friction law coupled with the “ageing law” of state evolution (Dieterich, 1981; Ruina, 1980, 1983): τ/σ = µ∗+ aln ?V ? ?V∗θ ? + bln (5) V∗ Dc ˙θ = 1 −V θ (6) Dc where τ and σ are the shear and normal fault stresses, respectively, V is slip rate and θ a fault state variable. The parameter a quantifies the direct effect, b the evolution effect, and Dcis the characteristic slip related to the state evolution. Under steady- state, b > a leads to stick-slip behavior (velocity-weakening, VW) whereas b < a leads to stable sliding behavior (velocity-strengthening, VS). We represent a seismogenic zone driven by surrounding creep by prescribing a VW patch of length Lvw= L/4 in the middle of the fault, surrounded by two VS segments and, at further distance, by steady uniform creep at slip rate Vpl. The parameter values of our numerical model are given in Supplementary Table C1. Ruptures that start as a crack and later turn into a pulse require a minimum rupture distance to develop the transition, therefore a sufficiently large ratio between –4–

Lvwand the nucleation length Lnuc(Rubin & Ampuero, 2005): Lnuc=π µDcb σ(b − a)2 (7) 2 Previous earthquake cycle simulations including a LVFZ model did not show significant reductions in rise-time for Lvw∼ 1.5Lnucand ∆ = 0.36 (Kaneko et al., 2011), so we extended the seismogenic length to Lvw ∼ 10Lnuc. Moreover, we explored values of ∆ ranging from moderate damage to the upper bound of current seismological observations (∆ = 0.5 − 0.9). We prescribed a minimum of five elements within an effective cohesive zone, µ∗Dc σb Lb=9π , (8) 32 where µ∗is an effective shear modulus that accounts for the LVFZ (Appendix C). We only consider results after several cycles, to avoid dependence on the arbitrarily prescribed initial conditions. 3.2 Results Complex slip patterns appear when damage is high (∆ > 0.7) and the fault zone is thin compared to the length of the seismogenic zone (2h < Lvw) (Fig. 3). Two signatures characterize the slip complexity: the promotion of pulse-like rupture and the re-rupture of previously healed fault segments during the same event. Pulses are defined by a drastic reduction of slip rate at a short distance behind the rupture front, leading to a short rise-time (Fig. 2b). Under friction laws different than rate-and-state, pulses can have a strict healing front where V = 0 m/s (Huang & Ampuero, 2011; Huang et al., 2014). We observe a systematic reduction of the average rise-time over a wide range of LVFZ thickness and damage values (Fig. 3). Short rise-times occur roughly within the range of LVFZ parameters that produce flat slip profiles in the static rupture models computed in Section 2 (Fig. 3). The re-rupture of previously healed fault segments is characterized by the emer- gence of secondary fronts propagating in the opposite direction to the main rupture front (Fig. 2b). These back-propagating fronts have a short rise-time and can re- rupture multiple times the same fault segment. Models with large seismogenic zones (Lvw >> Lnuc) can promote back-propagating fronts without requiring a LVFZ, but their rise-time is longer and their number of re-ruptures is small (Fig. 2d with Lvw∼ 100Lnuc). Events comprising a wide range of sizes develop in thick and highly damaged fault zones (Fig. 3c), where small events (L < Lvw) partially break the seismogenic zone from the edges. Small, non-characteristic events are known to emerge in rate-and-state friction models in homogeneous media and large seismogenic zones (Lvw >> Lnuc) (Cattania, 2019; Barbot, 2019). The LVFZ thickness and damage values promoting variable event magnitudes in our models are well understood by the increase in the Lvw/Lnuc ratio due to the reduction in Lnuc induced by the LVFZ (Fig. 3c). The smallest nucleation length is achieved in models with ∆ = 0.9 and 2h > Lvw, which have Lvw∼ 100Lnuc. 4 Discussion 4.1 Short-range stress transfer and the origin of pulses in a LVFZ Models with nearest-neighbour stress transfer, such as the Burridge-Knopoff model (Burridge & Knopoff, 1967), have been often used as a mechanical analog to –5–

earthquake rupture that generates complexity (Carlson et al., 1994). Under our cur- rent model parameters (Supplementary Table C1), ruptures propagate as pulses both in a nearest-neighbour model (Fig 2c) and in a fault-zone model with large damage (∆ = 0.9) (Fig. 2b). Here we show that the emergence of pulses in a LVFZ can be related to stress interactions approaching the nearest-neighbour regime across a wide range of slip wavelengths. The static stress transfer in a fault-zone model due to spatially-harmonic slip with wavelength k and unit amplitude is (Appendix A): K(k) =1 2µ(1 − ∆)|k|coth(h|k| + atanh(1 − ∆)) (9) Asymptotic analysis (Appendix B, Fig. 4) shows that at low k the stress transfer in a LVFZ tends to that of an intact homogeneous medium, whereas at high k it tends to that of a damaged homogeneous medium. In an intermediate regime, at wavelengths between 2πh and 2πh/(1−∆), the stress transfer is approximately nearest- neighbour. As ∆ increases, the relative bandwidth of the nearest-neighbour regime, which is ∝ 1/(1−∆), broadens, and the short rise-time observed in the nearest-neighbor model (Fig. 2c) appears in the LVFZ model as well. In other words, increasing the LVFZ damage level extends the bandwidth where pulse-like rupture can exist. The limiting case where ∆ → 1 analyzed in Section 2 represents an elastic layer of thickness 2h bounded by an infinitely rigid medium (Horowitz & Ruina, 1989). Stress interactions in that case are nearest-neighbour at wavelengths larger than ∼ 2πh. Such model is completely nearest-neighbour if the process zone size, the smallest characteristic length scale of slip, is larger than 2πh. 4.2 Origin of secondary rupture fronts The static solutions introduced in Section 2 provide insight on the origin of sec- ondary pulses. We showed there is a transition from crack-like to pulse-like rupture when the rupture size exceeds 2h, in an idealized situation where the only deformable medium is within the LVFZ. In a more realistic situation the medium outside the LVFZ is deformable as well. As the rupture continues growing to sizes much larger than 2h it transfers stress to the outer medium. We then expect a reverse transition from pulse-like behavior to the crack-like behavior of an intact homogeneous medium. Beyond this transition, slip increases in regions that were previously healed. There- fore, slip reactivation is expected there, leading to secondary rupture fronts. Because these secondary ruptures start small, they would also go through a pulse-like phase. Essentially, in the presence of a LVFZ, secondary pulses are necessary to complete the slip budget of a very large rupture. 4.3 Implications for rupture imaging Despite the ongoing increase in observations revealing the inherent complexity of large earthquakes (Meng et al., 2012; Ross et al., 2019), back-propagating fronts have been rarely reported (Beroza & Spudich, 1988; Meng et al., 2011). The slow slip rates and short length scales of the secondary fronts in our models suggest it would be challenging to resolve back-propagating fronts using current source-imaging techniques. Our results further motivate the quest for higher temporal and spatial resolution of source processes in finite-fault slip inversions. They can also guide future laboratory experiments to study connections between complex rupture and heterogeneity of the medium surrounding the fault, targeting the parameter space where our LVFZ model approaches the nearest-neighbour regime (Fig. 4). –6–

4.4 Rapid-Tremor-Reversals Slow slip and tremor phenomena offer a unique and systematic opportunity to observe complex slip patterns in slow motion. Our results suggest a highly compli- ant LVFZ model provides an explanation for tremor migration patterns observed in Cascadia and Japan. The back-propagating pulses identified in Fig. 2 qualitatively mimic observations of Rapid Tremor Reversals (RTRs) made in Cascadia and Japan during slow-slip events (Houston et al., 2011). Alternatively to previous models of RTR that require frictional heterogeneities (Luo & Ampuero, 2018), here the RTR patterns emerge from a quasi- static elastic effect in a fault with uniform frictional properties. There is seismological evidence that subduction megathrusts are surrounded by low velocity zones (Nedimovi´ c et al., 2003; Audet & Schaeffer, 2018). It is also observed that the Low Frequency Earthquakes (LFEs) that constitute tremors have a characteristic size (Watanabe et al., 2007; Chestler & Creager, 2017; J. C. Hawthorne et al., 2018). Characteristic scales also emerge in our LVFZ models, suggesting the testable hypothesis that the scales of LFEs are related to the properties of a LVFZ. Additional simulations show that back-propagating secondary pulses also occur in slow slip models (Fig. 2e). Introducing strengthening at high slip rate is a known approach to model slow-slip events (J. Hawthorne & Rubin, 2013). We added a linear velocity-strengthening term into the friction law (i.e., the fault strengthens ∝ V ), which is stronger than the logarithmic strengthening term of the conventional rate-and-state friction in equation (5). We chose a velocity-strengthening coefficient 106times larger than the radiation damping coefficient. We find that back-propagating fronts emerge during slow-slip events in a LVFZ model with the modified friction, although they are less vigorous than those observed in our fast-rupture results (Fig. 2). Further work is required to examine how damage zones quantitatively affect tremor migration patterns in more realistic slow-slip models. The most compliant fault-zone structures observed in strike-slip faults reach ∆ ∼ 0.85 (Li et al., 2007; Yang & Zhu, 2010), which is close to the minimum damage required by our model to show significant slip complexity (Fig. 3). For ∆ ∼ 0.85 and a reasonable fault-zone thickness 2h ∼ 100 m, the rupture length required to develop pulses and back-propagating secondary fronts must be larger than ∼ 2 km (Fig. 3). It is likely then that the quasi-static LVFZ effects described here do not operate during very small events. The properties of fault zones where RTRs are observed are harder to be resolved compared to crustal faults due to the larger depths involved. Dimensions of fault zones in subduction environments have been inferred from observations in exhumed subduction zones (Rowe et al., 2013) but their elastic properties remain poorly constrained. Receiver functions suggest that the vp/vsratio may increase over ∼ 75% due to over-pressurization of fluids within the low-velocity zone that surrounds regions where tremors are generated (Audet & Schaeffer, 2018). 4.5 Model limitations Further research is warranted to investigate whether the effects observed in our idealized fault zone model remain after releasing some of the simplifying assumptions, in particular the quasi-dynamic approximation and the tabular LVFZ geometry. In the absence of a LVFZ, quasi-dynamic simulations qualitatively agree with fully-dynamic simulations under a conventional Dieterich-Ruina friction law (Thomas et al., 2014). In contrast to quasi-dynamic models, dynamic simulations promote trapped waves in the presence of a LVFZ, which can perturb the dynamic stress on the fault leading to the promotion of healing fronts and pulse-like rupture (Huang & Ampuero, 2011; Huang et al., 2014). The dynamic effects of trapped waves may allow –7–

the slip complexity revealed here to operate over a broader range of LVFZ property values, including the lower, commonly observed levels of fault-zone damage. The structure of damage zones observed in the field is more sophisticated than a simple tabular region, usually following so-called flower structures (Savage & Brodsky, 2011; Mitchell & Faulkner, 2009). Moreover, LVFZ properties are not expected to be uniform along strike as the fault-zone thickness can vary with along-strike changes in fault geometry and the total amount of slip locally accumulated over time (Mitchell & Faulkner, 2009; C. Perrin et al., 2016; Ampuero & Mao, 2017). 5 Conclusions Our analytical arguments and simulation results show that the mechanism of gen- eration of rupture pulses by a damaged fault zone is robust, persistent across multiple earthquake cycles. The mechanism is not exclusively enabled by the dynamic stresses carried by trapped waves, but also emerges under quasi-static deformation. A formal analogy between a fault zone model and a nearest-neighbor (Burridge-Knopoff) model explains the emergent complexity: the latter is known to produce pulse-like rupture behaviour and, within a certain range of length scales, the stress transfer in a damaged fault zone is approximately nearest-neighbor. If this mechanism of pulse generation is dominant in natural faults, earthquake rise times should be proportional to fault zone thickness. Fault-zone effects can produce complex slip patterns, including back-propagating secondary rupture pulses that re-rupture previously healed fault segments. These back- propagating fronts have not been observed yet in large earthquakes, but they resemble rapid tremor reversals observed in Cascadia and Japan. They are also present in our slow-slip models in strongly compliant fault zones. Earthquake cycle modeling that incorporates the effects of fault-zone trapped-waves can help evaluate if pulse-like rupture and back-propagating pulses occur over a wider range of fault zone properties. Overall, fault-zone effects provide a simple mechanism to promote and sustain earthquake complexity, and a mechanical link between structural fault properties and seismicity. Figures –8–

(a) (b) (c) damaged medium intact medium fault zone thickness bilateral growth rupture length Figure 1. (a) Schematic representation of a fault zone. (b) Conceptualization of a fault zone as a simple tabular Low Velocity Fault Zone (LVFZ) model. The damaged and intact media have constant shear modulus, (1 − ∆)µ and µ, respectively. (c) Quasi-static rupture growth with uni- form stress drop in a LVFZ, showing a transition from crack-like (elliptical) to pulse-like (flat) slip profiles when the rupture length exceeds the LVFZ thickness. –9–

Figure 2. Spatiotemporal evolution of slip rate in the characteristic event of earthquake cycle models assuming (a) an intact homogeneous medium, (b) a LVFZ with ∆ = 0.9 and h = 32 m, (c) a nearest-neighbor model with ∆ = 0.99 and h = 50 m, (d) an intact homogeneous medium with ten times smaller nucleation length than (a), and (e) a slow-slip model with the same LVFZ as in (b). –10–

Figure 3. Properties of ruptures and seismicity in fault-zone models. (a) The average rise- time normalized by the total rupture duration, (b) average number of re-ruptures during an event, and (c) number of characteristic events over the total number of events are plotted as a function of damage level ∆ and fault-zone thickness h normalized by the size of the velocity- weakening fault segment Lvw. The rise-time is defined as the duration of slip rate exceeding 1 cm/s. Black contour lines in (a) are a semi-analytic prediction of the flatness of the slip profile, defined as the proportion of the fault length where slip exceeds 80% of its maximum. Charac- teristic events are defined as those events that completely break Lvw. The white contours in (c) show the estimated reduction of the nucleation length due to the LVFZ (contours of Lnuc in LVFZ normalized by its value in a homogeneous intact medium). –11–

Figure 4. Comparison between the static stress transfer kernels of a LVFZ with ∆ = 0.99 (black), homogeneous intact medium (blue dashed) and homogeneous damaged medium (orange dashed) in Fourier domain, as a function of the wavenumber k of slip. The black circles indicate the wavenumbers 10(1 − ∆)/h and 1/h that define a regime where the LVFZ kernel is similar to the nearest-neighbor kernel. –12–

A Derivation of the static stress transfer kernel in a LVFZ model A.1 Problem statement Let us consider a 2D elastic medium where the fault is located on the line y = 0. The medium is heterogeneous with the shear modulus µ depending only on y. Slip is anti-plane in the direction out of the x − y plane. A static slip distribution D(x) produces a shear stress on the fault T (x). Because the elasticity problem is linear, slip and stress are related by a linear relation: Z∞ where K is the static stress transfer kernel. The minus sign is introduced to give K a meaning analogous to stiffness in a spring-block system. As the problem is invariant by translation along x (K(x,x0) = K(x − x0)), equation (A.1) is a convolution: Z∞ After a Fourier transform, the convolution simplifies into a product: K(x,x0)D(x0)dx0, T (x) = − (A.1) −∞ K(x − x0)D(x0)dx0 T (x) = − (A.2) −∞ T (k) = −K(k)D(k), (A.3) where k is the wavenumber along the fault direction, x. In order to connect stress and slip in the fault, the goal is to derive an expression for the static kernel in spectral domain, the so-called spectral stiffness K(k). The equation of anti-plane elasticity governing the displacement field u(x,y) parallel to z is: µ(y)u,xx+ (µ(y)u,y),y= 0 (A.4) A first boundary condition is the symmetric distribution of applied slip on each side of the fault: u(x,0±) = ±1 A second boundary condition requires that displacement u must be finite at y = ±∞. As the boundary conditions are symmetric, the resulting displacements are symmetric as well and follow u(x,−y) = −u(x,y). Therefore, we solve for the half-plane y ≥ 0 only. After applying a Fourier transform over equation (A.4), our task is reduced to find u(k,y) such that: 2D(x) (A.5) −k2µ(y)u + (µ(y)u,y),y= 0 u(k,0+) =1 2D(k) (A.6) u(k,∞) < 0 and then evaluate the fault shear stress T (k) = µ(0)u(k,0),y. The problem has analytical solutions only for certain shear-modulus distributions µ(y). In the following we address two cases: a homogeneous medium and a two-layer medium. An analytical solution for an exponential distribution of µ(y) is possible as well but not exposed here (Ampuero et al., 2002). A.2 Homogeneous medium In a homogeneous medium, equation (A.4) reduces to: −k2u + u,yy= 0, (A.7) –13–

Its well-known general solution is u(k,y) = Aexp(−|k|y) + B exp(|k|y). The finite displacement boundary condition imposes B = 0, and the fault boundary condition at y = 0 implies A =1 After evaluating the shear stress on the fault, T (k) = −1 is: K(k) = −T (k)/D(k) =1 2D. The resulting displacement is u(k,y) =1 2D(k)exp(−|k|y). 2µ|k|D(k), the spectral stiffness 2µ|k| (A.8) A.3 Two-layer medium Consider a fault in a homogeneous medium surrounded by two layers of uniform half-thickness h and homogeneous but reduced shear modulus (Fig. 1b). Within the layers the shear modulus is µ(1 − ∆) and outside the layers it is µ. The derivation of the kernel in a layered medium follows the same steps as the previously addressed case. The differential equation is equation (A.7). Its general solution is a combination of exponential functions for each layer, together with a total of four new constants analogous to A and B. Two new boundary conditions arise: continuity of displace- ment and stress across the interface between the layers, u(k,h+) = u(k,h−) and (1 − ∆)u(k,h−),y= u(k,h+),y. It is possible to obtain the displacement D(k,y) after some algebraic work, then derive the shear stress at the fault and finally the spectral kernel: K(k) =1 2µ|k|(1 − ∆)coth(h|k| + atanh(1 − ∆)) (A.9) A first application of equation (A.9) is to numerically compute the slip profiles of a rupture with prescribed constant stress drop propagating in a fault with a LVFZ (Fig. 1c). By applying an inverse fast Fourier transform to equation (A.9) over a very long fault, we obtained a static stress transfer kernel in space domain, K(x). Then assuming a uniform stress drop, we solved numerically the discretized version of equation (A.2) to obtain the slip profiles. With a more significant effort, our numerical implementation of a LVFZ on multi- cycle earthquake simulations consists of combining the time-domain kernel of a fault with finite length (Cochard & Rice, 1997) with equation (A.9) in the frequency do- main. The numerical models shown in Fig. 2 are based on this implementation. We verified that the values of the obtained kernel are similar to those obtained by the more expensive approach of applying equation (A.9) to a periodic homogeneous fault 32 times longer (Fig. A.1). B The stress transfer in a LVFZ and a Burridge-Knopoff model B.1 Burridge-Knopoff (BK) model In a BK model (Burridge & Knopoff, 1967), the quasi-static slip and stress at the base of the i-th block of area dx2relate to each other as: Tidx2= −kL(Di− DL) +¯K(Di−1− 2Di+ Di+1) We furthermore introduce a loading stiffness per unit area of block surface defined as: (B.1) ¯KL= kL/dx2 (B.2) Taking the continuum limit (dx → 0) in equation (B.1): T (x) = −¯KL(D(x) − DL) +¯KD,xx (B.3) The second term of the right-hand side in the equation above is derived by writing: Di−1− 2Di+ Di+1= D(x − dx) − 2D(x) + D(x + dx) (B.4) –14–

10-3 3.5 Finite fault Trimmed periodic fault 3 2.5 2 1.5 1 0.5 0 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 Figure A.1. Static stress kernel in a LVFZ versus wavenumber based on two numerical im- plementations. The blue circles represent the combination of the kernel of a finite fault with equation (A.9) in the wavenumber domain. The continuous orange line is an approximated kernel using equation (A.9) over a periodic fault 32 times longer. –15–

and expanding terms in a Taylor series up to second order at small dx. Taking the Fourier transform for non-zero wavenumbers (|k| > 0): T (k) = −¯KLD(k) −¯Kk2D(k) = −(¯KL+¯Kk2)D(k) The loading displacement DL is spatially uniform, hence it only contributes when k = 0. We get the following static kernel in spectral domain: (B.5) K(k) =¯KL+¯Kk2 (B.6) B.2 Comparison between LVFZ and BK kernels For a LVFZ model, we rewrite the kernel given in equation (A.9) as: 2µd|k|1 + (1 − ∆)tanh(1 − ∆) (1 − ∆) + tanh(h|k|) K(k) =1 (B.7) In a highly damaged fault zone where ∆ → 1, the kernel reduces to: K(k) ≈1 1 2µd|k| (B.8) (1 − ∆) + tanh(h|k|) The high-frequency regime is defined by tanh(h|k|) ? 1 − ∆ and leads to: K(k) ≈1 2µd|k|coth(h|k|) (B.9) Moreover, if h|k| ? 1, by Taylor expansion we obtain: K(k) ≈µd 2h+1 6µdhk2 (B.10) This shows that, under certain conditions for h|k|, the stress transfer of the LVFZ model is equivalent to that of the BK model, with the following formal analogies: ¯KL=µd (B.11) 2h ¯K =µdh (B.12) 6 In the low-frequency regime defined by tanh(h|k|) ? 1 − ∆, the LVFZ is too narrow to have an effect in the kernel, leading to the kernel of an homogeneous medium: K(k) ≈1 2µ|k| (B.13) C Estimation of the process zone size in a LVFZ The components of slip and stress drop at a wavenumber k∗are related by an effective shear modulus, µ∗, which represents an effective value of the shear modulus profile µ(y) up to an off-fault distance comparable to λ∗= 2π/k∗. As is the case of a homogeneous medium, the static stress transfer at this effective length scale can be written as: K(k∗) =1 Combining equations (A.9) and (C.1), we obtain: 2µ∗|k∗| (C.1) µ∗= µ(1 − ∆)coth(h|k∗| + atanh(1 − ∆)) (C.2) An estimate of the process zone size in a LVFZ, L∗ in an intact medium, Lb, by: b, is related to the process zone size b= Lbµ∗ L∗ (C.3) µ –16–

After replacing equation (C.2) into equation (C.3) and setting λ∗≈ L∗ zone in a LVFZ satisfies: b, the process ?2πh ? L∗ b= Lb(1 − ∆)coth + atanh(1 − ∆) (C.4) L∗ b This equation is solved numerically to obtain L∗ b/Lbas a function of ∆ and h/Lb. Acknowledgments This work was supported by the Southern California Earthquake Center (SCEC) and by the French government through the FAULTS R GEMS project (ANR-17-CE31- 0008) and the UCAJEDI Investments in the Future project (ANR-15-IDEX-01) man- aged by the National Research Agency (ANR). SCEC is funded by NSF Cooperative Agreement EAR-0529922 and USGS Cooperative Agreement 07HQAG0008. Author contributions Both authors contributed to the writing of the manuscript and the interpretation of the numerical results. B.I. developed the scaling argument, performed the numer- ical simulations, and prepared the figures. J.-P. A. designed the study, developed the expressions for the spectral kernel, the static crack numerical solutions and the asymptotic analysis connecting the BK and LVFZ models. References Aagaard, B. T., & Heaton, T. (2008). Constraining fault constitutive behavior with slip and stress heterogeneity. Journal of Geophysical Research: Solid Earth, 113(B4). Ampuero, J. P., & Mao, X. (2017). Upper limit on damage zone thickness controlled by seismogenic depth. Fault Zone Dynamic Processes: Evolution of Fault Prop- erties During Seismic Rupture, 227, 243. Ampuero, J.-P., Vilotte, J.-P., & Sanchez-Sesma, F. under slip dependent friction law: simple models of fault zone. Journal of Geo- physical Research: Solid Earth, 107(B12). Audet, P., & Schaeffer, A. J. (2018). Fluid pressure and shear zone development over the locked to slow slip region in cascadia. eaar2982. Barbot, S. (2019). Slow-slip, slow earthquakes, period-two cycles, full and partial ruptures, and deterministic chaos in a single asperity fault. 768, 228171. Ben-Zion, Y., Peng, Z., Okaya, D., Seeber, L., Armbruster, J. G., Ozer, N., ... Ak- tar, M. (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. Beroza, G. C., & Mikumo, T. (1996). Short slip duration in dynamic rupture in the presence of heterogeneous fault properties. Solid Earth, 101(B10), 22449–22460. Beroza, G. C., & Spudich, P. (1988). havior: Application to the 1984 morgan hill, california, earthquake. Journal of Geophysical Research: Solid Earth, 93(B6), 6275–6296. Burridge, R., & Knopoff, L. (1967). Model and theoretical seismicity. Bulletin of the seismological society of america, 57(3), 341–371. Carlson, J. M., Langer, J. S., & Shaw, B. E. (1994). Dynamics of earthquake faults. Reviews of Modern Physics, 66(2), 657. Cattania, C. (2019). Complex earthquake sequences on simple faults. Geophysical Research Letters. (2002). Nucleation of rupture Science advances, 4(3), Tectonophysics, Journal of Geophysical Research: Linearized inversion for fault rupture be- –17–

Chestler, S., & Creager, K. (2017). Evidence for a scale-limited low-frequency earth- quake source process. Journal of Geophysical Research: Solid Earth, 122(4), 3099–3114. Cochard, A., & Madariaga, R. (1996). Complexity of seismicity due to highly rate- dependent friction. Journal of Geophysical Research: Solid Earth, 101(B11), 25321–25336. Cochard, A., & Rice, J. R. (1997). A spectral method for numerical elastodynamic fracture analysis without spatial replication of the rupture event. the Mechanics and Physics of Solids, 45(8), 1393–1418. Day, S. M. (1982). Three-dimensional finite difference simulation of fault dynamics: rectangular faults with fixed rupture velocity. Bulletin of the Seismological So- ciety of America, 72(3), 705–727. Dieterich, J. H. (1981). Constitutive properties of faults with simulated gouge. Me- chanical behavior of crustal rocks: the Handin volume, 24, 103–120. Erickson, B., Jiang, J., Barall, M., Lapusta, N., Dunham, E., Harris, R., ... others (2019). The community code verification exercise for simulating sequences of earthquakes and aseismic slip (seas). Harris, R. A., & Day, S. M. (1997). Effects of a low-velocity zone on a dynamic rup- ture. Bulletin of the Seismological Society of America, 87(5), 1267–1280. Hawthorne, J., & Rubin, A. M. (2013). in a rate and state friction model with a velocity-weakening to velocity- strengthening transition. Journal of Geophysical Research: Solid Earth, 118(7), 3785–3808. Hawthorne, J. C., Thomas, A. M., & Ampuero, J.-P. of low frequency earthquakes near parkfield, ca. tional, 216(1), 621–639. Heaton, T. H. (1990). Evidence for and implications of self-healing pulses of slip in earthquake rupture. Physics of the Earth and Planetary Interiors, 64(1), 1– 20. Horowitz, F. G., & Ruina, A. (1989). Slip patterns in a spatially homogeneous fault model. Journal of Geophysical Research: Solid Earth, 94(B8), 10279–10298. Houston, H., Delbridge, B. G., Wech, A. G., & Creager, K. C. (2011). Rapid tremor reversals in Cascadia generated by a weakened plate interface. science, 4(6), 404. Huang, Y., & Ampuero, J.-P. (2011). Pulse-like ruptures induced by low-velocity fault zones. Journal of Geophysical Research: Solid Earth, 116(B12). Huang, Y., Ampuero, J.-P., & Helmberger, D. V. modulated by waves in damaged fault zones. Journal of Geophysical Research: Solid Earth, 119(4), 3133–3154. Huang, Y., Ampuero, J.-P., & Helmberger, D. V. shear earthquakes in damaged fault zones–theory and observations. Planetary Science Letters, 433, 109–115. Johnson, E. (1990). On the initiation of unidirectional slip. Geophysical Journal In- ternational, 101(1), 125–132. 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: Solid Earth, 116(B10). Li, H., Zhu, L., & Yang, H. (2007). High-resolution structures of the Landers fault zone inferred from aftershock waveform data. tional, 171(3), 1295–1307. Luo, Y., & Ampuero, J.-P. (2018). Stability of faults with heterogeneous friction properties and effective normal stress. Tectonophysics, 733, 257–272. Luo, Y., Ampuero, J.-P., Galvez, P., van den Ende, M., & Idini, B. QDYN: a Quasi-DYNamic earthquake simulator (v1.1). (doi:10.5281/zenodo.322459) Journal of Laterally propagating slow slip events (2018). Geophysical Journal Interna- The rupture extent Nature Geo- (2014). Earthquake ruptures (2016). The potential for super- Earth and Geophysical Journal Interna- (2017). Zenodo. –18–

Meng, L., Ampuero, J., Page, M., & Hudnut, K. (2011). Seismological evidence and dynamic model of reverse rupture propagation during the 2010 m7. 2 el mayor cucapah earthquake. In Agu fall meeting abstracts. Meng, L., Ampuero, J.-P., Stock, J., Duputel, Z., Luo, Y., & Tsai, V. (2012). Earth- quake in a maze: Compressional rupture branching during the 2012 mw 8.6 sumatra earthquake. Science, 337(6095), 724–726. Mitchell, T., & Faulkner, D. (2009). surrounding strike-slip fault zones with a wide range of displacements: A field study from the Atacama fault system, northern Chile. Geology, 31(8), 802–816. Nedimovi´ c, M. R., Hyndman, R. D., Ramachandran, K., & Spence, G. D. Reflection signature of seismic and aseismic slip on the northern cascadia sub- duction interface. Nature, 424(6947), 416. Nielsen, S., & Madariaga, R. (2003). On the self-healing fracture mode. Bulletin of the Seismological Society of America, 93(6), 2375–2388. Noda, H., Dunham, E. M., & Rice, J. R. (2009). Earthquake ruptures with thermal weakening and the operation of major faults at low overall stress levels. nal of Geophysical Research: Solid Earth, 114(B7). Perrin, C., Manighetti, I., Ampuero, J.-P., Cappa, F., & Gaudemer, Y. (2016). Lo- cation of largest earthquake slip and fast rupture controlled by along-strike change in fault structural maturity due to fault growth. Journal of Geophysical Research: Solid Earth, 121(5), 3666–3685. Perrin, G., Rice, J. R., & Zheng, G. (1995). surface. Journal of the Mechanics and Physics of Solids, 43(9), 1461–1495. Rice, J. R. (1993). Spatio-temporal complexity of slip on a fault. physical Research: Solid Earth, 98(B6), 9885–9907. Ross, Z. E., Idini, B., Jia, Z., Stephenson, O. L., Zhong, M., Wang, X., ... others (2019). Hierarchical interlocked orthogonal faulting in the 2019 ridgecrest earthquake sequence. Science, 366(6463), 346–351. Rowe, C. D., Moore, J. C., Remitti, F., & Scientists, I. E. T. (2013). The thickness of subduction plate boundary faults from the seafloor into the seismogenic zone. Geology, 41(9), 991–994. Rubin, A., & Ampuero, J.-P. (2005). Earthquake nucleation on (aging) rate and state faults. J. Geophys. Res.-Sol. Ea., 110(B11). Ruina, A. (1980). Friction laws and instabilities: a quasi-static analysis of some dry friction behaviour. Ph. D. thesis, Division of Engineering, Brown University. Ruina, A. (1983). Slip instability and state variable friction laws. J. Geophys. Res.- Sol. Ea., 88(B12), 10359–10370. Savage, H. M., & Brodsky, E. E. (2011). placement of fracture distribution and secondary fault strands in fault damage zones. Journal of Geophysical Research: Solid Earth, 116(B3). Thomas, M. Y., Lapusta, N., Noda, H., & Avouac, J.-P. versus fully dynamic simulations of earthquakes and aseismic slip with and without enhanced coseismic weakening. Journal of Geophysical Research: Solid Earth, 119(3), 1986–2004. Watanabe, T., Hiramatsu, Y., & Obara, K. (2007). Scaling relationship between the duration and the amplitude of non-volcanic deep low-frequency tremors. physical research letters, 34(7). Yang, H., & Zhu, L. (2010). Shallow low-velocity zone of the San Jacinto fault from local earthquake waveform modelling. Geophysical Journal International, 183(1), 421–432. The nature and origin of off-fault damage Journal of Structural (2003). Jour- Self-healing slip pulse on a frictional Journal of Geo- Collateral damage: Evolution with dis- (2014). Quasi-dynamic Geo- –19–