# Diverse slip behaviour of velocity-weakening fault barriers

Diego Molina-Ormazabal, Jean-Paul Ampuero, Andrés Tassara

► To cite this version:

Diego Molina-Ormazabal, Jean-Paul Ampuero, Andrés Tassara. Diverse slip behaviour of velocity-weakening fault barriers. Nature Geoscience, 2023, 16 (12), pp.1200-1207. 10.1038/s41561-023-01312-1. hal-04885070

HAL Id: hal-04885070

https://hal.science/hal-04885070v1

Submitted on 3 Oct 2025

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Diverse slip behaviour of velocity-weakening fault barriers

Received: 22 March 2022

Accepted: 29 September 2023

Published online: 9 November 2023

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

Check for updates

Diego Molina-Ormazabal*<fnref n="1" />,²,<fnref n="4" />✉, Jean-Paul Ampuero*³ & Andrés Tassara*⁵

Seismic barriers are fault portions that promote earthquake rupture arrest and fault segmentation. Despite their fundamental role in controlling the maximum magnitude of earthquakes, the nature of seismic barriers is still uncertain. A common interpretation of barriers as having velocity-strengthening friction—steady-state friction that increases with increasing slip velocity—is only partially consistent with the thermal control of friction observed in laboratory experiments, which implies that most relevant materials in subduction channels are velocity-weakening at seismogenic depths. Here we examine the possibility of velocity-weakening barriers by conducting earthquake cycle simulations along a velocity-weakening megathrust segmented by lateral variations of frictional properties and normal stress. We show that velocity-weakening fault segments display a wide range of behaviours, including permanent barrier behaviour. They can be locked during long periods and release their slip deficit either seismically or aseismically. We quantify the efficiency of velocity-weakening barriers in arresting ruptures using a non-dimensional parameter based on fracture mechanics theory that can be constrained by observations on natural faults. Our results provide a theoretical framework that could improve physics-based seismic hazard assessment.

Understanding the mechanisms controlling earthquake arrest is crucial to elucidate the physics governing seismicity patterns and maximum earthquake size. Faults are often segmented into 'asperities' experiencing large slip during earthquakes and 'barriers' that prevent earthquake propagation<fnref n="1" />,². The origin of barriers in subduction megathrusts has been attributed to stress redistribution by past ruptures³,<fnref n="4" />, ductile behaviour at high temperatures⁵, geometrical fault bends⁶, subducted bathymetric anomalies⁷,⁸ and upper-plate geology⁹,¹⁰. Barriers frequently overlap with geological structures of the overriding plate¹¹,¹² and positive gravity anomalies¹³, suggesting a long-term, persistent barrier behaviour. However, some barriers lie in zones of negative or null gravity anomalies, potentially revealing a role of heterogeneous friction¹⁴,¹⁵. Moreover, barriers can display non-stationary behaviour. Most regions acting as barriers during historical earthquakes present diverse slip behaviours in

modern instrumental records, ranging from aseismic to seismic slip. Geodetic observations indicate that some past barriers host aseismic slip, either interseismically or postseismically⁵,¹⁰,¹⁶, but others have then participated in earthquake ruptures¹⁷,¹⁸. Seismic barriers have been categorized into persistent, frequent and ephemeral, according to the frequency of earthquakes traversing them³. Models with uniform fault friction properties can exhibit spontaneous segmentation and heterogeneous stresses¹⁹,²⁰ but cannot explain barriers that persist across several seismic cycles. Persistent barriers arise more naturally from heterogeneities of material properties controlling fault strength.

A diversity of fault slip patterns can be numerically simulated adopting the rate-and-state friction (RSF) law²¹,²². In this framework, fault stability is controlled by the dependence of steady-state friction on slip velocity. In general, velocity-strengthening (VS) regions, where

<footnote n="1">¹Departamento Ciencias de la Tierra, Facultad de Ciencias Químicas, Universidad de Concepción, Concepción, Chile. ²Millennium Nucleus, The Seismic Cycle Along Subduction Zones CYCLO, Valdivia, Chile. ³Université Côte d'Azur, IRD, CNRS, Observatoire de la Côte d'Azur, Géoazur, Valbonne, France.</footnote>

<footnote n="4">⁴Present address: Université Grenoble Alpes, CNRS, ISTerre, Grenoble, France. ✉e-mail: diego.molina@univ-grenoble-alpes.fr</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1 | Conceptual model and numerical simulation set-up.** a, Conceptual sketch (not to scale) of a VW (a – b < 0) subduction megathrust where seismic asperities are dominated by unstable sliding. b, A 2.5-dimensional (2.5D) model set-up of a fault with two seismogenic asperities (red bands), a central VW</figcaption></figure>

patch (yellow) whose normal stress $σ_n$ and RSF critical slip distance $D_c$ are varied systematically, and two VS (a – b > 0) regions (white) at the extremes of the fault. The fault is infinitely long, with spatially periodic slip along strike, and is loaded by steady slip (constant slip rate $V_{pl}$) at depths beyond the seismogenic width.

steady-state friction increases with increasing slip velocity, favour aseismic slip and earthquake arrest, while velocity-weakening (VW) zones promote earthquake nucleation and rupture propagation. RSF models commonly assume a patchy distribution of friction properties comprising VW asperities segmented by VS barriers²³⁻²⁵. A mechanical understanding of VS barriers is well established²⁶ and is being integrated in hazard assessment²⁷.

The VS-barrier model implies that large regions segmenting megathrusts must remain VS across the entire seismogenic depth, at temperatures between 150 and 350 °C (refs. 28,29). However, frictional laboratory experiments show that few of the materials that are relevant in the subduction environment (mostly clay minerals) remain VS at these temperatures³⁰,³¹, whereas most others (unaltered and altered rocks of the oceanic and continental plate, subducted sediments and so on) are VW³²⁻³⁴. In addition, some natural barriers exhibit intermediate to high geodetically estimated interseismic coupling (ISC)³⁵,³⁶, but that is expected for VS materials only in narrow stress shadow areas around locked asperities³⁷. Thus, it seems unlikely that all barriers along subduction zones are VS. The alternative is to consider that the entire seismogenic zone is VW, but asperities and barriers emerge from lateral variations of frictional strength (Fig. 1a).

Previous earthquake cycle simulations show that VW segments might promote earthquake arrest³⁸,³⁹, but the spectrum of slip behaviour displayed by VW barriers has not been characterized or understood. In this Article, we investigate slip patterns emerging on barriers along VW faults in earthquake cycle simulations as a function of parameters controlling barrier strength. We further recur to fracture mechanics theory to identify factors governing the efficiency of VW barriers in arresting ruptures. Finally, we discuss important implications of our results on natural seismicity patterns and their possible link with structural features of upper and lower plates.

## Simulating earthquake cycles around VW barriers

We study how along-strike variations of friction properties and normal stress on a VW fault can promote rupture arrest by performing quasi-dynamic simulations spanning many earthquake cycles over 5 kyr (Methods). Based on the conceptual model of Fig. 1a, our numerical model (Fig. 1b) comprises a VW fault segment surrounded by VS regions and loaded by uniform steady slip at depths beyond the seismogenic depth. The fault contains two asperities of length $L_A$. They are separated by a central VW patch of length $L_B$ with different effective normal stress $σ_n$ and RSF critical distance $D_c$ than the asperities. Values of frictional
<table><caption>**Table 1 | Frictional parameters employed on seismic cycle numerical simulations**</caption><thead><tr><th>Parameter</th><th>Nomenclature</th><th>Value</th></tr></thead><tbody><tr><td>Asperity size</td><td>L<sub>A</sub></td><td>70 km</td></tr><tr><td>Barrier size</td><td>L<sub>B</sub></td><td>20–60 km</td></tr><tr><td>VS section length</td><td>L<sub>vs</sub></td><td>40 km</td></tr><tr><td>Critical distance asperity</td><td>D<sub>c</sub></td><td>0.02</td></tr><tr><td>Normal stress asperity</td><td>σ<sub>n</sub></td><td>40 mPa</td></tr><tr><td>Critical distance barrier</td><td>D<sub>c</sub></td><td>0.02–0.40 m</td></tr><tr><td>Normal stress barrier</td><td>σ<sub>n</sub></td><td>40–150 mPa</td></tr><tr><td>a value asperity</td><td>a</td><td>0.017</td></tr><tr><td>a value barrier</td><td>a</td><td>0.015</td></tr><tr><td>a value VS section</td><td>a</td><td>0.032</td></tr><tr><td>b value along fault</td><td>b</td><td>0.019</td></tr><tr><td>Seismogenic width QDYN</td><td>W<sub>2.5d</sub></td><td>1,000 km</td></tr><tr><td>Velocity plate</td><td>V<sub>pl</sub></td><td>1.5×10<sup>-9</sup> m s<sup>-1</sup></td></tr><tr><td>Shear modulus</td><td>μ</td><td>30 GPa</td></tr></tbody></table>

properties and other model parameters are defined in Methods and Table 1.

Our simulations produce earthquake ruptures that agree with empirical scaling relations between rupture length, maximum slip and magnitude (Extended Data Fig. 1) and a wide spectrum of slip patterns spanning aseismic to seismic behaviours on the central patch, like those observed in natural faults (Fig. 2 and Extended Data Fig. 2). For example, the central patch can remain locked for several decades, manifested by high integrated long-term ISC (Methods) (Fig. 2b and Extended Data Fig. 2b), then release its accumulated slip deficit seismically or aseismically (Fig. 2c and Extended Data Fig. 2c). The mechanisms of slip on the central patch differ between simulations with varying $D_c$ and $σ_n$, as summarized in the following. Under certain conditions, the VW central patch can behave as a barrier. For each simulation, we quantify the barrier efficiency by the probability $P$ that an earthquake nucleating on an asperity breaks the whole central patch. Adopting the nomenclature proposed in a recent work³, we define ephemeral barriers by $P \ge 75\%$, frequent barriers by $25\% < P < 75\%$ and permanent barriers by $P \le 25\%$.

Central patches with large enough $D_c$ contrast relative to surrounding asperities behave as barriers. The barrier efficiency $P$ decreases with increasing $D_c$, and permanent barrier behaviour emerges at $D_c$ values larger than ~0.1 m (Fig. 3a). These permanent barriers accommodate the slip deficit mostly aseismically (Fig. 3c). Their predominant mechanism of aseismic slip is afterslip of earthquakes arrested on the central patch (Fig. 2a,b). Afterslip propagates along the central patch, usually penetrating the opposite asperity and inducing new earthquakes there. It propagates faster when $D_c$ on the central patch is larger (Extended Data Fig. 3c). At high $D_c$ contrasts, slow slip events can also contribute to release slip deficit, often triggered by coalescence of two afterslip fronts coming from each side of the barrier (Extended Data Fig. 4a,b). Further, on large $D_c$ central patches, small earthquakes emerge near the barrier edges (Extended Data Fig. 4b). This predominant transient aseismic slip behaviour on the barriers leads to highly variable interseismic locking from one cycle to the next, ranging from intermediate to high integrated ISC values (Fig. 2d and Extended Data Fig. 4c,f). Central patches having small $D_c$ contrast can behave as frequent or ephemeral barriers (Extended Data Fig. 4d–f) or simply as another seismic asperity (Extended Data Fig. 4g–i).

Central patches defined by a larger $σ_n$ release their slip deficit mostly seismically (Extended Data Fig. 2a–c), during earthquakes that

break them completely. The barrier efficiency $P$ decreases with increasing $σ_n$, but permanent barrier behaviour remains rare (Fig. 3b). As $σ_n$ is increased, barrier behaviour becomes more frequent, but the dominant slip behaviour remains seismic (Fig. 3c). Central patches with large $σ_n$ contrast remain locked during all the cycles (Extended Data Fig. 2b,d) and produce large seismic slip, even earthquake nucleation. They can also experience a small amount of aseismic slip in the form of afterslip confined mainly at their edges; afterslip penetrates farther into the barrier as $σ_n$ is increased but generally without spanning the whole barrier length (Extended Data Fig. 2a,b and Extended Data Fig. 5c,d). Central patches with low $σ_n$ contrast are easily crossed by earthquake ruptures (Extended Data Fig. 5e,f), thus behaving as another asperity instead of as a barrier.

## Fracture mechanics of simulated VW barriers

To better understand the physics governing VW barriers and to quantify their efficiency, we recur to linear elastic fracture mechanics (LEFM) (Methods). An earthquake can propagate if its static energy release rate $G_o$ exceeds the fracture energy $G_c$ dissipated by slip; otherwise, it stops. Our results indicate that most earthquakes stopping in the barrier satisfy the arrest condition $G_o < G_c$ (Fig. 4a). Noting that $G_c$ is proportional to $D_c$ and $σ_n$ (ref. 40) helps to explain why barrier efficiency $P$ increases with increasing $D_c$ or $σ_n$. Thus, the persistence of natural barriers could be due to spatial variations of frictional properties or normal stress that remain unchanged throughout several seismic cycles.

However, our finding that very efficient permanent barriers ($P < 20\%$) cannot be achieved by increasing indefinitely their $σ_n$ requires a different explanation. A possible alternative could be the effect of normal stress on recurrence time. The intrinsic recurrence time of an isolated fault segment is $Ti = (b - a)σ_n \log(V_{dyn}/V_{pl})W/(πμV_{pl})$ where a and b correspond to coefficients controlling the velocity and state effects, respectively, while $V_{pl}$ and $V_{dyn}$ are the plate velocity and dynamic velocity during seismic ruptures, respectively. Noting that $Ti$ is proportional to $σ_n$, we can infer that during some cycles, a segment with high $σ_n$ (longer recurrence) would stop ruptures nucleated on a neighbour fault section with lower $σ_n$ (shorter recurrence), which occurs partially in some simulations (Extended Data Fig. 5c,d). Nevertheless, we also find counterexamples: a simulation with $σ_n = 350$ Pa in the barrier has an intrinsic ratio of $Ti_A/Ti_B \approx 0.21$ but a substantially higher real ratio of ~0.5 (Extended Data Fig. 5a,b); the strong interaction between asperity and barrier makes them break together more often than if they were isolated. In addition, noting that the nucleation length is inversely proportional to $σ_n$, a higher $σ_n$ makes the central patch more unstable, promoting seismic slip. Notably, several natural barriers are in overlapping rupture ends, often experiencing coseismic slip¹⁸,³⁶. Hence, barrier behaviour with $σ_n$ contrast is controlled mostly by a competition between high fracture energy and small nucleation length rather than by differences in recurrence time.

Properties of adjacent asperities also affect a barrier's efficiency. Energetic ruptures can penetrate farther into a barrier, which we quantify by a higher asperity energy $G_{oi}$ defined as the static energy release rate measured when the rupture reaches the end of the initial asperity (Fig. 4b). The barrier/asperity stress drop ratio $ΔTb/ΔTa$ also influences rupture propagation inside a barrier. For instance, if $G_c/G_{oi}$ is fixed at 1.5, a higher $ΔTb/ΔTa$ promotes further rupture penetration (Fig. 4b). Assuming uniform stress drop inside the asperity and inside the barrier, LEFM predicts that the barrier efficiency depends on three non-dimensional parameters: $L_A/L_B$, $G_c/G_{oi}$ and $ΔTb/ΔTa$ (Methods). On the basis of our ensemble of simulations (Methods), we found that the following 'barrier efficiency index' $λ$ combining these three non-dimensional parameters is an adequate predictor of the probability of earthquakes breaking through the barrier (Fig. 4c):

$$ \lambda = \frac{L_A}{L_B} \frac{0.025}{G_c} \frac{-0.25}{G_{oi}} \frac{\Delta\tau_b}{\Delta\tau_a} $$
<figure><img src="image_4.png" /><figcaption>**Fig. 2 | Results of seismic cycle simulations producing a permanent barrier.** The permanent barrier is defined by $D_c = 0.15$ m. The barrier size is $L_B = 45$ km. **a,** Accumulated slip (y axis) along the fault (x axis) during interseismic periods (blue, every two years) and during earthquakes (red, every five seconds); black lines depict the start and end of earthquakes. Vertical lines display the edges of the barrier. **b,** Temporal evolution (y axis) of the instantaneous ISC (Methods) along the simulated fault (x axis). Blue solid lines show the extent of earthquake</figcaption></figure>

ruptures and green stars their epicentres. Dashed white lines indicate the edges of the central patch. **c,** Maximum slip rate across the barrier (black solid lines, log scale) as a function of time. The horizontal red line indicates our chosen threshold to define seismic slip rates, 1 cm s⁻¹. Aseismic transients are those that do not reach this threshold. **d,** Instantaneous ISC averaged along the barrier (red solid lines). ISC in the central patch varies throughout the seismic cycle, although the values are most often near 1 (fully locked).

Barriers behave mostly as permanent if $\lambda < 0.6$, as frequent if $\lambda \approx 0.75$ and as ephemeral if $\lambda > 0.85$.

To evaluate the efficiency of fault barriers in nature, we propose an estimator of $\lambda$ that can be constrained by geodetically estimated ISC distributions (Supplementary Information). Our estimator $\lambda'$ relies on scaling relations between energy quantities and coseismic slip $S$. We follow a fracture energy scaling¹¹ $G_c \approx BS^l$, where $B \approx 3$ and $l \approx 0.7$ are constants constrained by seismological observations. We adopt a theoretical estimate of the asperity energy (Methods) as $G_{oi} = \frac{\pi L_A \Delta \tau_A^2}{4\mu}$ where $\Delta T_A \approx \mu S_A / L_A$ is the stress drop on the asperity and $S_A$ its average coseismic slip. Accordingly, we define a barrier efficiency proxy (Methods):

$$ \lambda' = \left[ \frac{L_A}{L_B} \right]^{0.025} \left[ \frac{4L_A BS_B^n}{\pi \mu S_A^2} \right]^{-0.25} \left[ \frac{S_B L_A}{S_A L_B} \right]^{-0.38} $$

where $S_B$ is the coseismic final slip averaged on the barrier.

We illustrate the approach by evaluating the efficiency of the Arauco barrier in the Chilean subduction zone, around the southern end of the 2010 moment magnitude 8.8 Maule earthquake rupture (Supplementary Fig. 1).

The chosen parameters based on a formal segmentation model¹⁵ result in $\lambda' = 0.61$ for the Arauco barrier, which according to our model corresponds to permanent barriers, supporting the observation that this area has persistently stopped ruptures in historical and prehistoric records¹².

## Interseismic coupling pattern of VW barriers

Our results show that areas of high and low ISC do not necessarily correspond to asperities and barriers, respectively, challenging a common interpretation of geodetic observations⁵. The geodetically observed ISC spans time windows much shorter than an entire seismic cycle and is not necessarily representative of the long-term behaviour. Our
<figure><img src="image_5.png" /><figcaption>**Fig. 3** | The diversity of slip behaviours in VW barriers. a,b, Probability (P in percentage) of earthquake ruptures crossing the barrier as a function of D_c in central patches defined by a D_c contrast (a) or in barriers defined by a σ_n contrast (b). c, Probability (P) of earthquake ruptures crossing the barrier as a function of the ratio of aseismic to seismic slip averaged inside the central patch, for barriers</figcaption></figure>

defined by high D_c and by high σ_n. Horizontal grey lines separate permanent (P ≤ 25%), frequent (25% < P < 75%) and ephemeral (P ≥ 75%) barriers. Permanent barriers are obtained by increasing D_c but not by increasing σ_n. Symbol sizes are scaled by L_B = length of the central patch. Simulations with D_c contrast are displayed as circles, and simulations with σ_n variations are shown as squares.

simulations spanning multiple earthquake cycles support the idea that ISC can be time dependent despite fault mechanical properties being constant, owing to slow slip transients and extended afterslip<sup>42,43</sup>.

Remarkably, the integrated long-term ISC in most of our simulations is high, in both permanent and frequent barriers (Fig. 5). In general, permanent VW barriers (*P* < 15%) feature ISC > 0.9 (Fig. 5), in contrast with the lower values (ISC < 0.85) produced by permanent VS barriers<sup>26</sup>. When the low ISC of a barrier is a transient feature, subsequent locking can lead to increased seismic hazard. Conversely, areas with high or intermediate ISC are not necessarily asperities as elevated ISC can also happen on barriers. Three broad persistent barriers globally identified<sup>3</sup> (two in Chile and one in Sumatra) have intermediate to high geodetic ISC<sup>36,44</sup>. Therefore, VW barriers can host complex slip behaviour.

Numerical studies<sup>20,45,46</sup> have shown that the ratio between the asperity length (*L*) and the characteristic nucleation length *L*<sub>c</sub> = μ*D*<sub>c</sub>/*(b − a)*σ<sub>n</sub> is a good predictor of slip complexity in isolated VW faults, but this ratio seems insufficient to predict barrier behaviour in our models, especially for simulations with σ<sub>n</sub> contrast. In isolated VW faults, low values of *L*/*L*<sub>c</sub> promote aseismic slip, often related to barriers. However, our models indicate that central patches in segmented faults can behave as frequent barriers even if *L*<sub>B</sub>/*L*<sub>c</sub> > 10 (Extended Data Fig. 6).

Our simulations with high D_c demonstrate that the central VW patch can experience high long-term ISC, and earthquakes and creep can be collocated. Part of this complexity in slip behaviour has been previously reported<sup>46</sup>, showing that for 1.5 < L<sub>B</sub>/L<sub>C</sub> < 14 the fault may host complex slip behaviour. However, barriers with 0.3 < L<sub>B</sub>/L<sub>C</sub> < 1.5 might also exhibit slip complexity such as high ISC accompanied with rupture arrest. These findings suggest that D_c or σ_n heterogeneities

along the seismogenic zone can generate seismic barriers, like those observed in nature, for a more ample range of values of L/L<sub>c</sub> than the ones previously reported.

Some natural barriers also slip aseismically (low ISC). In our simulations, wide barriers with very high D_c can present very low ISC during long periods (Extended Data Fig. 4b,c). However, D_c is unlikely to reach such very high values in nature; thus, long-term creeping barriers are more likely related to VS areas. However, if the aseismic behaviour observed geodetically on a natural barrier does not correspond to its long-term behaviour, it could also be interpreted as a VW permanent barrier with increased D_c currently undergoing a long period of aseismic slip. Hence, geodetic observations spanning longer times are needed to assess whether natural barriers are VS or VW, which ultimately affects their impact on seismic hazard.

Barriers and asperities have been widely associated with VS and VW fault areas, respectively<sup>23-26</sup>. Further research efforts should aim at assessing whether the segmentation of subduction zones is instead driven primarily by heterogeneities of VW materials and therefore re-evaluating the common assumption that barriers are VS areas.

## Implications for structural control of earthquake behaviour

Several observations suggest that geological structures of upper and lower plates may control the distribution of asperities and barriers along subduction zones<sup>9-13</sup>. One possible explanation is that tectonic features affect the distribution of strength of the megathrust, which in turn affects seismicity patterns. Properties controlling fault strength include σ_n and D_c, which, as shown in our analysis, modulate barrier efficiency via their effect on the ratios G_c/G_oi and ΔTb/ΔT_a.
<figure><img src="image_6.png" /></figure>

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

<figure><img src="image_8.png" /><figcaption>**Fig. 4 | Rupture arrest conditions and barrier efficiency based on fracture mechanics.** a, Fracture energy ($G_c$: Methods) and energy release rate ($G_o$: Methods) numerically evaluated at the final rupture tip of earthquakes that stop at barriers defined by high $D_c$ and by high $\sigma_n$. For reference, the solid grey line shows a 1/1 relation. In most cases, the rupture arrest condition $G_c \ge G_o$ is verified. b, Ratio between the fracture energy of the barrier ($G_c$) and the static energy release rate of the asperity rupture ($G_{oi}$) as a function of total rupture length ($L$) normalized by the asperity length ($L_\lambda$) for all the events that stopped within the barrier. Points are simulation results, with $G_c$ measured (using equation (3); Methods) at the end of the rupture and $G_{oi}$ evaluated (using equation (4); Methods) at the end of the asperity. Curves are theoretical estimates of $G_c/G_{oi}$ (equation (12) in Methods),</figcaption></figure>

assuming different values of the ratio of barrier stress drop to asperity stress drop (see legend). c, Barrier efficiency: probability $P$ of earthquakes crossing the barrier as a function of a dimensionless number $\lambda$ that encapsulates the effects of energy ratio, stress drop ratio and aspect ratios for barriers defined by a contrast of $D_c$ (circles) and $\sigma_n$ (squares). Thick horizontal light blue lines separate ephemeral ($P \ge 75\%$), frequent ($25\% < P < 75\%$) and permanent ($P \le 25\%$) barrier types. The vertical grey lines show the transition from permanent to frequent and frequent to ephemeral barriers as a function of $\lambda$. The green circle highlights permanent barriers for the Arauco region. The black curve represents the fitted polynomial function used to construct $\lambda$ (Methods). Simulations with $D_c$ contrast are displayed as circles, and simulations with $\sigma_n$ variations are shown as squares.

<figure><img src="image_9.png" /><figcaption>**Fig. 5 | Long-term interseismic coupling in the central patch.** For each simulation, the long-term ISC is computed at each point inside the central patch (Methods) and then averaged spatially across it. Most values are higher than 0.95, implying that in the long term, some barriers can display high ISC. The shaded area represents the range of results for VS barriers²⁶. Simulations with $D_c$ contrast are displayed as circles, and simulations with $\sigma_n$ variations are shown as squares.</figcaption></figure>

Linking these findings with geological features can give insights on barrier behaviour in nature, for example, recognizing that areas with large $\sigma_n$ can promote earthquake arrest. Density heterogeneities in the upper plate are an important source of lateral variations of $\sigma_n$ (refs. 14,15,47). On the basis of analysing gravity anomalies¹⁵, it has been proposed that the Arauco barrier in Chile has higher $\sigma_n$ than surrounding areas. Prehistoric and historic earthquake ruptures coming from north and south overlap underneath Arauco Peninsula, driving the zone to participate in coseismic slip¹⁸,³⁶. This is consistent with both the description of permanent barriers found in our work and our simulations showing seismic ruptures on central patches with increased $\sigma_n$. Normal stress may also increase due to dilatancy strengthening, which has been proposed to explain VW barriers along transform faults⁴⁸ and could be extended to subduction megathrusts.

There is yet no clear path to infer $D_c$ at natural scales from geological features. In nature, $D_c$ can depend on fault roughness, gouge thickness, $\sigma_n$ (ref. 22), thermally activated dynamic weakening⁴⁹ and fluid content⁵⁰. Some of these factors might be linked to subducted features that contribute to fault roughness and supply wet sediments to the subduction channel¹⁷, thus inducing along-strike variations in $D_c$. Indeed, barriers usually coincide with subducting seamount chains, fracture zones or aseismic ridges⁷,⁸,¹⁷. Seismicity observed at the flanks of aseismic ridges¹⁷, where a pronounced megathrust roughness could induce larger $D_c$, is consistent with small earthquakes nucleating often at the edge of barriers with high $D_c$ in our simulations.

We have emphasized the fundamental role that temperature could have in the frictional behaviour along megathrusts, which is not fully incorporated in our model. Indeed, barrier efficiency can be modified
by the competition between thermally activated friction⁵¹ and healing mechanisms³² and by elastodynamics and fault weakening processes relevant at seismic slip rates⁵³. Therefore, although our physics-based models capture important features of seismicity, additional ingredients could be incorporated to further advance our understanding of fault segmentation.

Finally, our study supports that geological characteristics of subducting and overriding plates can impact fault behaviour. It also advances physics-based seismic hazard assessment by providing a theoretical framework to use future understandings of how geological structures modulate the spatial distribution of properties controlling fault strength, such as *D*<sub>c</sub> and *σ*<sub>n</sub>, and to anticipate the efficiency of natural fault barriers.

## Online content

Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-023-01312-1.

## References

1. Aki, K. Characterization of barriers on an earthquake fault. *J. Geophys. Res.* **84**, 6140 (1979).

2. Lay, T. & Kanamori, H. in *Earthquake Prediction: An International Review* Vol. 4 (eds Simpson, D. W. & Richards, P. G.) 579–592 (AGU, 1981).

3. Philibosian, B. & Meltzner, A. J. Segmentation and supercycles: a catalog of earthquake rupture patterns from the Sumatran Sunda megathrust and other well-studied faults worldwide. *Quat. Sci. Rev.* **241**, 106390 (2020).

4. Konca, A. et al. Partial rupture of a locked patch of the Sumatra megathrust during the 2007 earthquake sequence. *Nature* **456**, 631–635 (2008).

5. Wang, K. & Bilek, S. L. Invited review paper: fault creep caused by subduction of rough seafloor relief. *Tectonophysics* **610**, 1–24 (2014).

6. Wesnousky, S. Predicting the endpoints of earthquake ruptures. *Nature* **444**, 358–360 (2006).

7. Contreras-Reyes, E. & Carrizo, D. Control of high oceanic features and subduction channel on earthquake ruptures along the Chile–Peru subduction zone. *Phys. Earth Planet. Inter.* **186**, 49–58 (2011).

8. Sparkes, R., Tilmann, F., Hovius, N. & Hillier, J. Subducted seafloor relief stops rupture in South American great earthquakes: implications for rupture behavior in the 2010 Maule, Chile earthquake. *Earth Planet. Sci. Lett.* **298**, 89–94 (2010).

9. Jara-Muñoz, J., Melnick, D., Brill, D. & Strecker, M. R. Segmentation of the 2010 Maule Chile earthquake rupture from a joint analysis of uplifted marine terraces and seismic-cycle deformation patterns. *Quat. Sci. Rev.* **113**, 171–192 (2015).

10. Saillard, M. et al. From the seismic cycle to long-term deformation: linking seismic coupling and Quaternary coastal geomorphology along the Andean megathrust. *Tectonics* **36**, 241–256 (2017).

11. Collot, J.-Y. et al. Are rupture zone limits of great subduction earthquakes controlled by upper plate structures? Evidence from multichannel seismic reflection data acquired across the northern Ecuador–southwest Colombia margin. *J. Geophys. Res.* **109**, B11103 (2004).

12. Melnick, D., Bookhagen, B., Strecker, M. R. & Echtler, H. P. Segmentation of megathrust rupture zones from fore-arc deformation patterns over hundreds to millions of years, Arauco Peninsula, Chile. *J. Geophys. Res.* **114**, B01407 (2009).

13. Bassett, D. & Watts, A. B. Gravity anomalies, crustal structure, and seismicity at subduction zones: 2. Interrelationships between forearc structure and seismogenic behavior. *Geochem. Geophys. Geosyst.* **16**, 1541–1576 (2015).

14. Tassara, A. Control of forearc density structure on megathrust shear strength along the Chilean subduction zone. *Tectonophysics* **495**, 34–47 (2010).

15. Molina, D., Tassara, A., Abarca, R., Melnick, D. & Madella, A. Frictional segmentation of the Chilean megathrust from a multivariate analysis of geophysical, geological, and geodetic data. *J. Geophys. Res. Solid Earth* **126**, e2020JB020647 (2010). (2021).

16. Plata-Martinez, R. et al. Shallow slow earthquakes to decipher future catastrophic earthquakes in the Guerrero seismic gap. *Nat. Commun.* **12**, 3976 (2021).

17. Bassett, D. & Watts, A. B. Gravity anomalies, crustal structure, and seismicity at subduction zones: 1. Seafloor roughness and subducting relief. *Geochem. Geophys. Geosyst.* **16**, 1508–1540 (2015).

18. Ojeda, J., Ruiz, S., del Campo, F. & Carvajal, M. The 21 May 1960 $M_w$ 8.1 Concepción earthquake: a deep megathrust foreshock that started the 1960 Central-South Chilean seismic sequence. *Seismol. Res. Lett.* **91**, 1617–1627 (2020).

19. Horowitz, Franklin G. & Ruina, Andy Slip patterns in a spatially homogeneous fault model. *J. Geophys. Res. Solid Earth* **94**, 10279–10298 (1989).

20. Cattania, C. Complex earthquake sequences on simple faults. *Geophys. Res. Lett.* **46**, 10384–10393 (2019).

21. Dieterich, J. H. Time-dependent friction and the mechanics of stick-slip. *J. Geophys. Res.* **116**, 790–806 (1978).

22. Marone, C. Laboratory-derived friction laws and their application to seismic faulting. *Annu. Rev. Earth Planet. Sci.* **26**, 643–646 (1998).

23. Lui, S. & Lapusta, N. Repeating microearthquake sequences interact predominantly through postseismic slip. *Nat. Commun.* **7**, 13020 (2016).

24. Wei, M. & Pengcheng, S. Synchronization of earthquake cycles of adjacent segments on oceanic transform faults revealed by numerical simulation in the framework of rate-and-state friction. *J. Geophys. Res. Solid Earth* **126**, e2020JB020231 (2021).

25. Lambert, V. & Lapusta, N. Resolving simulated sequences of earthquakes and fault interactions: implications for physics-based seismic hazard assessment. *J. Geophys. Res. Solid Earth* **126**, e2021JB022193 (2021).

26. Kaneko, Y., Avouac, J. P. & Lapusta, N. Towards inferring earthquake patterns from geodetic observations of interseismic coupling. *Nat. Geosci.* **3**, 363–369 (2010).

27. Michel, S., Jolivet, R., Rollins, C., Jara, J. & Dal Zilio, L. Seismogenic potential of the main Himalayan thrust constrained by coupling segmentation and earthquake scaling. *Geophys. Res. Lett.* **48**, e2021GL093106 (2021).

28. Oleskevich, D. A., Hyndman, R. D. & Wang, K. The updip and downdip limits to great subduction earthquakes: thermal and structural models of Cascadia, south Alaska, SW Japan, and Chile. *J. Geophys. Res.* **104**, 14965–14991 (1999).

29. Hippchen, S. & Hyndman, R. D. Thermal and structural models of the Sumatra subduction zone: implications for the megathrust seismogenic zone. *J. Geophys. Res.* **113**, B12103 (2008).

30. Saffer, D. M. & Marone, C. Comparison of smectite- and illite-rich gouge frictional properties: application to the updip limit of the seismogenic zone along subduction megathrusts. *Earth Planet. Sci. Lett.* **215**, 219–235 (2003).

31. Okamoto, A. S. et al. Frictional properties of simulated chlorite gouge at hydrothermal conditions: implications for subduction megathrusts. *J. Geophys. Res. Solid Earth* **124**, 4545–4565 (2019).

32. Liu, Y. & Rice, J. R. Slow slip predictions based on granite and gabbro friction data compared to GPS measurements in northern Cascadia. *J. Geophys. Res.* **114**, B09407 (2009).
33. Phillips, N. J., Belzer, B., French, M. E., Rowe, C. D. & Ujiie, K. Frictional strengths of subduction thrust rocks in the region of shallow slow earthquakes. *J. Geophys. Res. Solid Earth* **125**, e2019JB018888 (2020).

34. den Hartog, S. A. M., Niemeijer, A. R. & Spiers, C. J. New constraints on megathrust slip stability under subduction zone P–T conditions. *Earth Planet. Sci. Lett.* **353–354**, 240–252 (2012).

35. Métois, M., Vigny, C. & Socquet, A. Interseismic coupling, megathrust earthquakes and seismic swarms along the Chilean subduction zone (38°–18°S). *Pure Appl. Geophys.* **173**, 1431–1449 (2016).

36. Moreno, M., Rosenau, M. & Oncken, O. 2010 Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone. *Nature* **467**, 198–202 (2010).

37. Hetland, E. A. & Simons, M. Post-seismic and interseismic fault creep II: transient creep and interseismic stress shadows on megathrusts. *Geophys. J. Int.* **181**, 99–112 (2010).

38. Hori, T. & Miyazaki, S. A possible mechanism of M9 earthquake generation cycles in the area of repeating M7~8 earthquakes surrounded by aseismic sliding. *Earth Planet Space* **63**, 773–777 (2011).

39. Boatwright, J. & Cocco, M. Frictional constraints on crustal faulting. *J. Geophys. Res.* **101**, 13895–13909 (1996).

40. Rubin, A. M. & Ampuero, J.-P. Earthquake nucleation on (aging) rate and state faults. *J. Geophys. Res.* **110**, B11312 (2005).

41. Weng, H. & Ampuero, J. P. Continuum of earthquake rupture speeds enabled by oblique slip. *Nat. Geosci.* **13**, 817–821 (2020).

42. Meltzner, A. J. et al. Time-varying interseismic strain rates and similar seismic ruptures on the Nias–Simeulue patch of the Sunda megathrust. *Quat. Sci. Rev.* **122**, 258–281 (2015).

43. Tsang, L. L. H., Meltzner, A. J., Hill, E. M., Freymueller, J. T. & Sieh, K. A paleogeodetic record of variable interseismic rates and megathrust coupling at Simeulue Island, Sumatra. *Geophys. Res. Lett.* **42**, 10585–10594 (2015).

44. Prawirodirdjo, L., McCaffrey, R., Chadwell, C. D., Bock, Y. & Subarya, C. Geodetic observations of an earthquake cycle at the Sumatra subduction zone: role of interseismic strain segmentation. *J. Geophys. Res.* **115**, B03414 (2010).

45. Rubin, A. M. Episodic slow slip events and rate-and-state friction. *J. Geophys. Res.* **113**, B11414 (2008).

46. Barbot, S. Slow-slip, slow earthquakes, period-two cycles, full and partial ruptures, and deterministic chaos in a single asperity fault. *Tectonophysics* **768**, 228171 (2019).

47. Li, D. & Liu, Y. Modeling slow-slip segmentation in Cascadia subduction zone constrained by tremor locations and gravity anomalies. *J. Geophys. Res. Solid Earth* **122**, 3138–3157 (2017).

48. Liu, Y., McGuire, J. J. & Behn, M. D. Aseismic transient slip on the Gofar transform fault, East Pacific Rise. *Proc. Natl Acad. Sci. USA* **117**, 10188–10194 (2020).

49. Brantut, N. Dilatancy-induced fluid pressure drop during dynamic rupture: direct experimental evidence and consequences for earthquake dynamics. *Earth Planet. Sci. Lett.* **538**, 116179 (2020).

50. Blanpied, M. L., Lockner, D. A. & Byerlee, J. D. Frictional slip of granite at hydrothermal conditions. *J. Geophys. Res.* **100**, 13045–13064 (1995).

51. Wang, B. & Barbot, S. Pulse-like ruptures, seismic swarms, and tremorgenic slow-slip events with thermally activated friction. *Earth Planet. Sci. Lett.* **603**, 117983 (2023). 2023.

52. Barbot, S. A rate-, state-, and temperature-dependent friction law with competing healing mechanisms. *J. Geophys. Res. Solid Earth* **127**, e2022JB025106 (2022).

53. Thomas, M. Y., Lapusta, N., Noda, H. & Avouac, J.-P. Quasi-dynamic versus fully dynamic simulations of earthquakes and aseismic slip with and without enhanced coseismic weakening. *J. Geophys. Res. Solid Earth* **119**, 1986–2004 (2014).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor (e.g. a society or other partner) holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

© The Author(s), under exclusive licence to Springer Nature Limited 2023
Methods

Earthquake cycle model set-up

We consider a fault governed by the RSF law with the aging law for the state evolution:

$$ \frac{\tau}{\sigma_n} = \mu^* + a \ln \left( \frac{V}{V^*} \right) + b \ln \left( \frac{V^* \theta}{D_c} \right) \quad (1) $$

$$ \frac{d\theta}{dt} = 1 - \frac{V\theta}{D_c} \quad (2) $$

Here $\tau$ and $\sigma_n$ are the frictional strength and normal stress on the fault, respectively, $V$ is the sliding velocity and $\theta$ is a state variable, $\mu^*$ and $V^*$ are reference values for friction and velocity, respectively, $D_c$ is a critical slip distance for state evolution, and $a$ and $b$ are positive coefficients controlling the velocity and state effects, respectively. The difference ($a-b$) defines the frictional regime on the fault. If $a-b>0$, the steady-state friction is VS: slip acceleration strengthens the fault, which inhibits earthquake nucleation and propagation, therefore promoting stable slip. If $a-b<0$, steady-state friction is VW: friction drops if slip accelerates, promoting unstable slip leading to earthquake nucleation and rupture propagation. Simulations considering the RSF capture spontaneously all the stages of the seismic cycle, including interseismic coupling, nucleation of coseismic ruptures, postseismic afterslip and slow slip transients.

We employ the QDYN software, which simulates earthquake cycles under the quasi-dynamic approximation and using adaptive time-stepping⁵⁴. The deformation is purely dip slip (antiplane), and frictional properties and normal stress are spatially variable along strike. Beyond the fault portion shown in Fig. 1b, our computational method assumes spatial periodicity along strike. We set VS buffer zones near the ends of the VW segment that act as permanent barriers to prevent ruptures from interacting strongly with their periodic replicates. The fault is embedded in an elastic, linear and isotropic medium, and it is loaded by uniform steady slip beyond a width $W_{2,5D}$. We capture the essence of the three-dimensional (3D) fault problem while keeping the computations affordable by doing 2.5D simulations: we model a 1D fault in a 2D elastic medium with modifications of the static stress transfer kernel that incorporate approximately the effect of a finite seismogenic width. We set model parameters to typical values coming from geophysical observations. The imposed frictional width is $W_{2,5D} = 1,000$ km, implying a real seismogenic width of $W = 200$ km. We set the elastic shear modulus $\mu = 30$ GPa, plate velocity $V_{pl} = 1.5 \times 10^{-9}$ m s⁻¹ (also set as the reference velocity $V^*$).

We considered an elementary subduction zone model with two VW asperities separated by a central VW section with different $\sigma_n$ and $D_c$. In the asperities, normal stress is $\sigma_n = 40$ MPa and critical slip distance $D_c = 2$cm. The two asperities have a fixed length $L_A = 70$ km. On the central section, $\sigma_n$ and $D_c$ are systematically varied, and its length $L_B$ is varied from 20 to 60 km (Fig. 1). Along the whole fault, $b = 0.019$. Most of the fault is VW ($a-b<0$), with $a = 0.015$ on the asperities and $a = 0.017$ on the central patch, except for two VS zones ($a-b>0$) at the outer ends of the asperities of length $L_{as} = 40$ km and $a = 0.032$. We define seismic behaviour on the basis of a slip rate threshold: we say that a point on the fault experiences seismic slip if its slip rate exceeds $1 \text{ cm s}^{-1}$ (Extended Data Fig. 7). We define transient aseismic slip by slip rates $> 10^{-9} \text{ m s}^{-1}$. We say the fault is 'strongly locked' when slip velocities are less than $10^{-10} \text{ m s}^{-1}$, well below the assumed plate velocity $1.5 \times 10^{-9} \text{ m s}^{-1}$.

We first ran reference simulations spanning 5 kyr, assuming uniform $\sigma_n$ and $D_c$ and varying only $L_B$. No barrier for earthquake propagation was observed, and the ruptures broke the whole fault (Extended Data Fig. 8). We then performed one set of simulations varying $\sigma_n$ between 40 and 150 MPa and another set varying $D_c$ between 0.02 and 0.4 m. Although in laboratory experiments $D_c$ is on the order of $10^{-5}$ to $10^{-3}$ m (refs. 21,55) or from a few hundred micrometres to a few

centimetres (ref. 56), seismologically inferred values of $D_c$ are of order 1 to 10 m (ref. 57). This demonstrates the ample range that $D_c$ could have. Due to the high computational cost of simulations with small $D_c$ values, we focus on the $D_c$ range between $10^{-3}$ and $10^{-1}$ m, thereby maintaining sufficient numerical resolution and representative values inferred from laboratory experiments and seismological observations.

The smallest size of earthquakes that can be produced in rate-and-state faults (with the aging law) is controlled by the nucleation length $L_\infty$ (ref. 40), which is proportional to $D_c$ and inversely proportional to $\sigma_n$. Since we are interested in the seismic cycle of large earthquakes interacting with barriers, rather than in reproducing small earthquake patterns (swarms, foreshocks, aftershocks and so on), we chose parameter values such that the nucleation length is sufficiently small to generate whole-asperity earthquakes (Table 1). In our modelled asperities, $L_\infty = 5,644$ m, and the ratio between segment length and nucleation length is $L_A/L_\infty \approx 12.4$, inhibiting the occurrence of very small earthquakes and favouring whole-asperity ruptures. This allows the study of the co-evolution between big earthquakes and a central patch barrier and keeps the computational cost affordable. In the central patch, the ratio $L_B/L_\infty$ varies according to $L_B$, $\sigma_n$ and $D_c$. For the shortest central patch considered here, $L_B = 20$ km, it ranges from 0.04 to 3.5.

The number of elements on the fault is 2,048 in simulations with $\sigma_n < 110$ MPa and $D_c > 0.02$ m, and 4,096 in other cases. The grid-cell size $\Delta x$ ranges from 60 to 136 m. The ratio between the process zone size $L_B$ (ref. 40) and grid-cell size is $L_B/\Delta x \approx 14$ in the asperities, and in the barrier, it varies from ~4 to ~158, providing sufficient numerical resolution⁵⁸.

The range of parameter values used in this work ($a, b, D_c, \sigma_n$) have been considered in recent earthquake cycle simulation studies and correspond to values reported in laboratory experiment²³,⁵⁹⁻⁶¹.

Energy balance for earthquakes

LEFM provides a theoretical framework to study rupture propagation and arrest. Whether an earthquake can break through a barrier depends on the balance between the mechanical energy available in the incoming rupture (its static energy release rate $G_o$) and the energy dissipated by slip in the barrier (the fracture energy $G_c$). Earthquakes propagate if $G_o \ge G_c$; otherwise they stop. To prove that this energy condition is satisfied in our simulations, for those earthquakes stopping on the central VW patch, we compute $G_c$ and $G_o$ as:

$$ G_c = \int_0^{\delta u} (\tau - \tau_o) du, \quad (3) $$

where $\tau$ is the fault shear stress as a function of slip and $\tau_o$ as a function of the residual stress at the final slip, and

$$ G_o = \frac{K_o^2}{2\mu} \quad (4) $$

where $K_o$ is the stress intensity factor. For a rupture located between $x=0$ and $x=1$, $K_o$ is computed by a weighted integral of the spatial distribution of stress drop $\Delta\tau(x)$ as⁶²

$$ K_o = \sqrt{\frac{a}{\pi}} \int_0^{2a} \Delta\tau \left[ \frac{1}{\sqrt{a^2 - (x-a)^2}} + \frac{x-a}{a\sqrt{a^2 - (x-a)^2}} \right] dx \quad (5) $$

where $a=1/2$. The weight function gives a larger contribution to the stress drop in areas located close to the rupture front.

To obtain more insights on the parameters controlling earthquake arrest in the central VW barrier, we approximate Go by assuming uniform stress drop in the asperity, $\Delta\tau_a$, and in the barrier, $\Delta\tau_b$. Then, denoting $L_A$ the asperity length, equation (5) becomes
$$ K_o = \sqrt{\frac{a}{\pi}} \int_0^l \Delta \tau_a \left[ \frac{1}{\sqrt{1-(\frac{x-a}{a})^2}} + \frac{x-a}{a\sqrt{1-(\frac{x-a}{a})^2}} \right] dx/a \\ + \sqrt{\frac{a}{\pi}} \int_{L_\Lambda}^l \Delta \tau_b \left[ \frac{1}{\sqrt{1-(\frac{x-a}{a})^2}} + \frac{x-a}{a\sqrt{1-(\frac{x-a}{a})^2}} \right] dx/a \quad (6) $$

$$ \text{Using } \int \frac{dn}{\sqrt{1-n^2}} = \arcsin(n) + c \text{ and } \int \frac{n dn}{\sqrt{1-n^2}} = -\cos[\arcsin(n)] + c, $$ $$ K_o = \sqrt{\frac{a}{\pi}} (\Delta \tau_a [\arcsin(\frac{x-a}{a})]_{L_\Lambda}^l - \cos(\arcsin(\frac{x-a}{a}))]_{L_\Lambda}^l) \\ + \Delta \tau_b [\arcsin(\frac{x-a}{a})]_{L_\Lambda}^l - \cos(\arcsin(\frac{x-a}{a}))]_{L_\Lambda}^l] \quad (7) $$

$$ K_o = \sqrt{\frac{a}{\pi}} \left[ \begin{aligned} & (\arcsin(\frac{l_\Lambda - a}{a}) - \arcsin(-1)) - \cos(\arcsin(\frac{l_\Lambda - a}{a})) \\ & + \cos(\arcsin(-1))) \Delta \tau_a + (\arcsin(1) - \arcsin(\frac{l_\Lambda - a}{a})) \\ & - \cos(\arcsin(1)) + \cos(\arcsin(\frac{l_\Lambda - a}{a})) \end{aligned} \right] \Delta \tau_b \quad (8) $$

$$ K_o = \sqrt{\frac{a}{\pi}} \left[ \begin{aligned} & (\arcsin(\frac{l_\Lambda - a}{a}) - \cos(\arcsin(\frac{l_\Lambda - a}{a}))) (\Delta\tau_a - \Delta\tau_b) \\ & + \frac{\pi}{2} (\Delta\tau_a + \Delta\tau_b) \end{aligned} \right] \quad (9) $$

Finally, the static energy release rate is

$$ G_o = \frac{\frac{a}{\pi} \left[ \left( \arcsin\left(\frac{l_\Lambda - a}{a}\right) - \cos\left(\arcsin\left(\frac{l_\Lambda - a}{a}\right)\right) \right) (\Delta\tau_a - \Delta\tau_b) + \frac{\pi}{2} (\Delta\tau_a + \Delta\tau_b) \right]^2}{2\mu} \quad (10) $$

For reference, we denote as $G_{oi}$ the static energy release rate of a rupture that breaks only the asperity ($l = L_\Lambda$):

$$ G_{oi} = \frac{\pi l_{\Lambda} \Delta \tau_{\Lambda}^2}{4\mu} \quad (11) $$

## Barrier efficiency $\lambda$

We construct a dimensionless parameter to quantify the resistance of the VW patch to earthquake propagation. Using equations (10) and (11) and defining the barrier-to-asperity stress drop ratio $T = \frac{\Delta\tau_b}{\Delta\tau_a}$,

$$ G_o/G_{oi} = \frac{\frac{l}{\pi^2 l_o} \left[ \arcsin\left(\frac{2l_\Lambda}{l} - 1\right) - \cos\left(\arcsin\left(\frac{2l_\Lambda}{l} - 1\right)\right) (1-T) + \frac{\pi}{2} (1+T) \right]^2}{\frac{\pi^2 l_\Lambda}{l_o} \left[ \arcsin\left(\frac{2l_\Lambda}{l} - 1\right) - \cos\left(\arcsin\left(\frac{2l_\Lambda}{l} - 1\right)\right) (1-T) + \frac{\pi}{2} (1+T) \right]^2} \quad (12) $$

We define an energy ratio

$$ \Gamma = G_c/G_0 \\ = G_c/G_{oi} \frac{\pi^2 l_\Lambda / l}{\left[ \arcsin\left(\frac{2l_\Lambda}{l} - 1\right) - \cos\left(\arcsin\left(\frac{2l_\Lambda}{l} - 1\right)\right) (1-T) + \frac{\pi}{2} (1+T) \right]^2} \quad (13) $$

The equilibrium condition $\Gamma = 1$ gives a non-dimensional equation of the form

$$ f\left(\frac{l}{L_{\Lambda}}, T, \frac{G_c}{G_{oi}}\right) = 0, \quad (14) $$

Thus, the normalized rupture arrest length $\frac{l}{L_\Lambda}$ is controlled by the energy ratio $G_c/G_{oi}$ and the stress drop ratio $T$:

$$ \frac{l}{L_{\Lambda}} = g\left(\frac{G_c}{G_{oi}}, T\right) \quad (15) $$

Intuitively, we expect that the probability of the barrier being crossed by an earthquake ($P$) is a function of the ratio $\frac{l}{L_B}$ between the

rupture arrest length and the barrier size: the larger is $\frac{l}{L_B}$, the more likely is the whole barrier to break often, and the higher is $P$. Dividing equation (15) by $L_B/L_\Lambda$, we see that $\frac{l}{L_B}$ is a function of $G_c/G_{oi}$, $T$ and $L_B/L_\Lambda$.

This particular example suggests that the resistance of a barrier to rupture penetration in general depends on the three non-dimensional parameters $G_c/G_{oi}$, $T$ and $L_B/L_\Lambda$. However, a closed-form expression of equation (15) cannot be obtained for an arbitrary distribution of stress drop. To find a conveniently compact mathematical expression that adequately explains our results, we assumed that the barrier resistance can be approximated in the following form:

$$ \lambda_i = \frac{L_\Lambda}{L_B} \frac{G_c}{G_{oi}} \frac{m}{T^n} \quad (16) $$

For each combination of exponents $k$, $m$ and $n$, we calculated $\lambda_i$ for all simulations, then fitted the data of $P$ with a fourth-order polynomial $F(\lambda) = a\lambda + b\lambda^2 + c\lambda^3 + d\lambda^4$, but only considering the data with $P < 99\%$. We then put back the removed data ($P > 99\%$) and computed the sum of squared residuals between the whole data and the constrained function $P = \min(100, F)$. The smallest sum of squared residuals is obtained for the values $k = 0.025$, $m = -0.25$ and $n = -0.38$. Thus, we define the 'barrier efficiency index' as

$$ \lambda = \frac{L_{\Lambda}}{L_{B}} \frac{0.025}{G_{c}} \frac{-0.25}{G_{oi}} \frac{\Delta \tau_{b}}{\Delta \tau_{a}} \frac{-0.38}{\Delta \tau_{a}} \quad (17) $$

## Computation of the barrier efficiency index $\lambda$

The barrier efficiency index $\lambda$ depends on the static energy release rate of the asperity, $G_{oi}$, which we measured in our simulations for all the ruptures arresting in the barrier or breaking the whole central patch by evaluating equations (4) and (5) when the rupture reaches the end of the asperity where it nucleated, at a rupture length of $L_\Lambda$. For events that stop immediately when the rupture encounters the barrier or at some distance within it, we measure $G_c$ at the rupture end. For events that break the whole barrier and propagate into the neighbouring asperity, we compute $G_c$ when the rupture reaches the end of the barrier. In all the cases, we measure $G_{oi}$ at the end of the first asperity (at the asperity-barrier limit). We also measure in each rupture the mean stress drop in the asperity ($\Delta\tau_a$) and in the barrier ($\Delta\tau_b$). For each multi-cycle simulation, we compute averaged values of these quantities among all events and then compute $\lambda$ using equation (17).

## Calibration of effective rupture width in the 2.5D model

We have performed 2.5D models on QDYN, which approximately includes the seismogenic width effect while maintaining low computational costs. The 2.5D model approximates the depth profile of slip as a half-sinusoid across the seismogenic width $W$. The slip in 2.5D models corresponds to the peak slip at the middle of the seismogenic depth of a 3D model. However, since the sinusoid approximation is crude, the effective width $W_{2.5D}$ used as input in QDYN differs from the real $W$. In this section, we calibrate the relation between $W_{2.5D}$ and $W$ by determining the ratio $W_{2.5D}/W$ such that the maximum slip $D_{max}$ of a very long rupture with uniform stress drop $\Delta\tau$ is the same in both 2.5D and 3D models. In the 2.5D model, the stress transfer kernel $K(k) = \frac{\Delta\tau(k)}{D(k)}$ as a function of the spatial wavenumber $k$ is given in equation A.24 of ref. 63. For very long wavelengths ($k = 0$) it is $K(0) = \pi \frac{\mu}{W_{2.5D}}$, thus $D_{max} = \frac{1}{\pi} \frac{\Delta\tau W_{2.5D}}{\mu}$. We verified this relation numerically. In the 3D model, the peak slip of a very long dip-slip rupture is that of a mode II crack in 2D (vertical cross-section plane). From ref. 64, $D_{max} = \frac{C \Delta\tau W}{\mu}$, where $C = 1.5$ if the Poisson's ratio is $\nu = 1/4$. For arbitrary values of $\nu$, we infer that $C = 2(1 - \nu)$. Equating the 2.5D and 3D expressions of $D_{max}$, we obtain
$$W_{2,5D} = 2\pi(1 - \nu)W \qquad (18)$$

For a Poisson's ratio of $\nu = 1/4$, we get

$$W_{2,5D}/W \approx 4.7 \qquad (19)$$

Therefore, the $W_{2,5D}$ input of QDYN is larger than the real seismic width $W$ by a factor of about 5.

## ISC

The instantaneous value of ISC is defined here as the fraction of slip deficit rate accommodated by stable creep²⁶:

$$ISC(x, t) = 1 - \frac{V(x, t)}{V_{\text{pl}}} \qquad (20)$$

where $V$ is the slip rate and $V_{\text{pl}}$ the long-term slip velocity imposed by plate convergence. The instantaneous ISC($x, t$) varies throughout a seismic cycle. We can compute it at all times and at all fault points in our simulations, but we show it only in areas with positive values ($V < V_{\text{pl}}$).

The long-term ISC is calculated at any point on the fault, following ref. 26, as

$$ISC(x) = 1 - \frac{\delta_{\text{int}}^{\text{cum}}(x)}{V_{\text{pl}} T_{\text{int}}^{\text{cum}}(x)} \qquad (21)$$

where $T_{\text{int}}^{\text{cum}}$ is the sum of the duration of all interseismic periods, defined as periods when $V < V_{\text{pl}}$, and $\delta_{\text{int}}^{\text{cum}}$ is the slip accumulated over all the interseismic periods.

## Data availability

The raw data coming from simulation results of figures in the main text (Figs. 3–5) are available through the Source data provided on the repository at https://doi.org/10.5281/zenodo.8329497.

## Code availability

The rate-and-state earthquake simulator (QDYN) used in this work is free, and it is freely available at https://github.com/ydluo/qdyn. Plots of figures were performed in Python using the package Matplotlib⁶⁵.

## References

54. Luo, Y. & Ampuero, J.-P. Stability of faults with heterogeneous friction properties and effective normal stress. Tectonophysics 733, 257–272 (2017).

55. Ohnaka, M. A constitutive scaling law and a unified comprehension for frictional slip failure, shear fracture of intact rock, and earthquake rupture. J. Geophys. Res. 108, 2080 (2003).

56. Brantut, Nicholas et al. High-velocity frictional properties of a clay-bearing fault gouge and implications for earthquake mechanics. J. Geophys. Res. Solid Earth 113, B10401 (2008).

57. Fukuyama, E., Mikumo, T. & Olsen, B. Estimation of the critical slip-weakening distance: theoretical background. Bull. Seismol. Soc. Am. 93, 1835–1840 (2003).

58. Perfettini, H. & Ampuero, J.-P. Dynamics of a velocity strengthening fault region: implications for slow earthquakes and postseismic slip. J. Geophys. Res. Solid Earth 113, B09411 (2008).

59. Barbot, S., Lapusta, N. & Avouac, J.-P. Under the hood of the earthquake machine: toward predictive modeling of the seismic cycle. Science 336, 707–710 (2012).

60. Noda, H. & Lapusta, N. Stable creeping fault segments can become destructive as a result of dynamic weakening. Nature 493, 518–521 (2013).

61. Beall, A., van den Ende, M., Ampuero, J.-P., Capitanio, F. A. & Fagereng, Å. Linking earthquake magnitude-frequency statistics and stress in visco-frictional fault zone models. Geophys. Res. Lett. 49, 20 (2022).

62. Rice, J. R. The mechanics of earthquake rupture. In Physics of the Earth’s Interior, Proc. Int. School of Physics ‘Enrico Fermi’ (eds Dziewonski, A. M. & Boschi, E.) 555–649. (Italian Physical Society and North-Holland, 1980).

63. Luo, Y. & Ampuero, J.-P. Stability of faults with heterogeneous friction properties and effective normal stress. Tectonophysics 733, 257–272 (2018).

64. Dieterich, J. H. Earthquake nucleation on faults with rate- and state-dependent strength. Tectonophysics 211, 115–134 (1992).

65. Hunter, J. D. Matplotlib: a 2D graphics environment. Comput. Sci. Eng. 9, 90–95 (2007).

66. Allen, T. I. & Hayes, G. P. Alternative rupture-scaling relationships for subduction interface and other offshore environments. Bull. Seismol. Soc. America 3, 1240–1253 (2017).

## Acknowledgements

This work was partially funded by the National Agency for Research and Development (ANID)/Scholarship Program/DOCTORADO BECAS CHILE/2017—21171169 to D.M.-O. We also acknowledge partial support by the Millennium Scientific Initiative (ICM) of the Chilean government through grant NCN19_167 ‘Millennium Nucleus CYCLO The Seismic Cycle Along Subduction Zones’ to D.M.-O. and A.T. We are also thankful to the Geoazur Seismes team, which also partially funded the work. J.-P.A. has been supported by the French government through the UCAJEDI Investments in the Future project (ANR-15-IDEX-01) managed by the National Research Agency (ANR).

## Author contributions

D.M.-O. designed and carried out the numerical simulations. J.-P.A. analysed the numerical experiments. J.-P.A. and D.M.-O. developed the theoretical model and wrote the main text. A.T. provided the technical motivation and contributed to writing the text.

## Competing interests

The authors declare no competing interests.

## Additional information

Extended data is available for this paper at https://doi.org/10.1038/s41561-023-01312-1.

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41561-023-01312-1.

**Correspondence and requests for materials** should be addressed to Diego Molina-Ormazabal.

**Peer review information** Nature Geoscience thanks the anonymous reviewers for their contribution to the peer review of this work. Primary Handling Editor: Tom Richardson, in collaboration with the Nature Geoscience team.

**Reprints and permissions information** is available at www.nature.com/reprints.
a

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

b

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

**Extended Data Fig. 1** | Comparison of source scalings between numerical slip models and natural observations. Rupture maximum slip (b) and length (a) versus earthquake magnitude. Stars depict results of numerical models. Solid lines show the empirical scaling relations for subduction earthquakes⁶⁶, dashed

lines their uncertainty (standard deviation). The magnitude has been calculated considering $M_o = \mu A \Delta S$, where $A = L \times W$ (Rupture length and constant seismogenic depth of 200 km respectively). $\Delta S$ corresponds to the rupture averaged slip.
<figure><img src="image_12.png" /></figure>

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

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

**Extended Data Fig. 2 |** Results of seismic cycle simulations producing a frequent barrier. Frequent barrier characterized by $σ_n=100 MPa$. The barrier size is $L_B=45$ km. a) Accumulated slip (y-axis) along the fault (x-axis) during interseismic periods (blue, every 2 years) and during earthquakes (red, every 5 seconds); Black lines depict the start and end of earthquakes. Vertical lines display the edges of the barrier. b) Temporal evolution (y-axis) of the instantaneous interseismic coupling (ISC, calculated as explained in the Methods) along the simulated fault (x-axis). Blue solid lines show the extent of

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

earthquake ruptures and green stars their epicenters. Dashed white lines indicate the edges of the central patch. c) Maximum slip rate across the barrier (black solid lines, log scale) as a function of time. The horizontal red line indicates our chosen threshold to define seismic slip rates, 1 cm/s. Aseismic transients are those that do not reach this threshold. d) Instantaneous ISC averaged along the barrier (red dashed lines). ISC in the central patch varies throughout the seismic cycle, although the values are most often near 1 (fully locked).
<figure><img src="image_16.png" /></figure>

**Extended Data Fig. 3** | The diversity of slip behavior in VW barrier types. P as a function of the percentage of rapid creep slip (stars and squares) and slow creep slip (diamonds and circles) relative to the total accumulated slip (Methods). A large transient slip relative to the total accumulated slip means an increase in the amount of aseismic slip hosted by the central patch (x-axis). As D_c is increased, the aseismic slip increases too. The same is observed in simulation with normal stress contrast (with values up to twice the asperity σ_n), but when σ_n exceeds three times the asperities' value, the transient slip decreases, implying an

increment in the seismic slip within the barrier. Further, a smaller ratio between transient and total accumulated slip means major seismic slip within the barrier. Vertical gray solid lines display barrier type depending on probability P. For permanent barrier with D_c constrast, there is a bifurcation between rapid and slow creep. This is reflected as a rapid afterslip entering in barrier with high D_c, which occur more slowly in barriers with D_c contrast no that high (frequent barriers).
<figure><img src="image_17.png" /><figcaption>**Extended Data Fig. 4** | Seismic cycle simulations of barriers with $D_c$ contrast. **a,b,c,** The barrier has a length of 45 km and is defined by very high $D_c$ = 0.3m. **d,e,f,** The barrier has a length of 45 km and is defined by high $D_c$ = 0.07m. **g,h,i,** The barrier has a length of 45 km and is defined by small $D_c$ = 0.04m. **a,d,g:** Accumulated slip profiles along the fault during interseismic periods (blue, every 2 years) and during earthquakes (red, every 5 seconds); Black lines depict the start and ending of earthquakes. Solid black lines display the edges of the CVWP.</figcaption></figure>

**b,e,f):** space-time distribution of instantaneous ISC (Methods). Blue solid lines and green stars show the rupture length and epicenter of earthquakes, respectively. Dashed white lines indicate the edges of the barrier. **c,f,i):** Instantaneous ISC averaged across the central patch (red dashed lines) and slip rate averaged across the whole barrier (black solid lines) as a function of time. ISC varies throughout the seismic cycle, although the values are most often close to 1 (fully locked).
a

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

c

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

e

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

b

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

d

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

f

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

Extended Data Fig. 5 | See next page for caption.
**Extended Data Fig. 5** | Seismic cycle simulations of barriers with $\sigma_n$ contrast.

**a,b,** The barrier is defined by high $\sigma_n$ = 350MPa. **c,d,** The barrier is defined by high $\sigma_n$ = 140MPa. **e,f,** The barrier is defined by small $\sigma_n$ = 70MPa. All cases depict barrier length of 30 km. **a,c,e:** Accumulated slip profiles along the fault during interseismic periods (blue, every 2 years) and during earthquakes (red,

every 5 seconds); Black lines depict the start and ending of earthquakes. Solid black lines display the edges of the central patch. **b,d,f:** space-time distribution of instantaneous ISC (Methods). Blue solid lines and green stars show the rupture length and epicenter of earthquakes, respectively. Dashed white lines indicate the edges of the barrier.
<figure><img src="image_24.png" /><figcaption>Extended Data Fig. 6 | Probability of ruptures crossing the barrier as a function of ratio $L_B/L_{CB}$. Squares and dots represent simulations with variable $D_c$ and $\sigma_n$ respectively.</figcaption></figure>
a

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

b

<figure><img src="image_26.png" /><figcaption>**Extended Data Fig. 7** | The typical slip rate evolution across seismic cycle. a)</figcaption></figure>

The barrier has a length of 40 km and is defined by high $D_c = 0.15m$. b) The
barrier has a length of 40 km and is defined by $\sigma_n = 100MPa$. Maximum slip rate

across the whole barrier (black solid line) and seismic asperities (right and left asperities, blue and red line respectively). The horizontal light green line indicates our chosen threshold for seismic slip rates, 1 cm/s.
<figure><img src="image_27.png" /></figure>

**Extended Data Fig. 8 | Reference seismic cycle simulation.** Cumulative slip profiles along the fault during interseismic periods (blue, every 2 years) and during earthquakes (red, every 5 seconds); Black lines depict the start and ending of earthquakes. Reference model with uniform *D*<sub>c</sub>=0.02 m and *o*<sub>n</sub>=40 MPa (there is a small increase of *a* in a 45 km-long central patch).