# Inelastic deformation accrued over multiple seismic cycles: Insights from an elastic-plastic slider-and-springboard model

Pierre Dublanchet, Jean-Arthur Olive

► To cite this version:

Pierre Dublanchet, Jean-Arthur Olive. Inelastic deformation accrued over multiple seismic cycles: Insights from an elastic-plastic slider-and-springboard model. Seismica, 2024, 3 (2), 10.26443/seismica.v3i2.1345. hal-04740233

HAL Id: hal-04740233

https://minesparis-psl.hal.science/hal-04740233v1

Submitted on 3 Nov 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
# Inelastic deformation accrued over multiple seismic cycles: Insights from an elastic-plastic slider-and-springboard model

Pierre Dublanchet ID * 1, Jean-Arthur Olive ID 2

¹Mines Paris, PSL University, Centre for geosciences and geoengineering, Fontainebleau, France, ²École normale supérieure, PSL University, Laboratoire de Géologie CNRS, Paris, France

Author contributions: Conceptualization: Pierre Dublanchet, Jean-Arthur Olive. Software: Pierre Dublanchet. Formal Analysis: Pierre Dublanchet, Jean-Arthur Olive. Investigation: Pierre Dublanchet, Jean-Arthur Olive. Writing - Original draft: Pierre Dublanchet, Jean-Arthur Olive.

**Abstract** We study a toy model designed to build physical insight into the problem of slow accumulation of non-recoverable strain in fault blocks over multiple earthquake cycles. The model consists of a thin, horizontal elastic-plastic plate (springboard) in frictional contact with a vertical, rigid wall moving downward at a steady speed. Our model produces stick-slip cycles consisting of interseismic plate downwarping and coseismic plate upwarping as long as the moment of the frictional force at the contact does not exceed the maximum (purely plastic) bending moment the plate can sustain. We show that the duration of individual earthquake cycles and the spatial pattern of interseismic deflection are controlled by two stress ratios involving the peak yield stress of the plate, the frictional strength of the fault and the coseismic stress drop. We show that non-recoverable plate deflection accumulates over successive earthquake cycles if the plate's yield strength decreases through time, causing a progressive decrease of the aforementioned stress ratios. We derive scaling relations between the rate of accumulation of inelastic deformation, the relative tectonic plate velocity, and the rate of lithospheric weakening. Our results are consistent with observations of long-term permanent deformation of natural fault regions.

Production Editor:

Carmine Galasso

Handling Editor:

Mathilde Radiguet

Copy & Layout Editor:

Hannah F. Mark

Received:

March 29, 2024

Accepted:

September 23, 2024

Published:

October 15, 2024

## 1 Introduction

Faults that dissect the Earth's brittle upper crust can accommodate large displacements, yet spend most of the time in a locked state. This allows the relative motion of tectonic plates to cause slow, distributed deformation of fault blocks that store elastic energy, ultimately released when faults slip seismically or aseismically (Reid, 1910). In subduction zones, for instance, the dominant deformation signal measured by geodetic methods is related to mega-earthquake cycles along the plate interface (Savage, 1983; Wang et al., 2012), which consist of 3 phases (Prawirodirdjo et al., 2010; Watanabe et al., 2021; Sagiya and Meneses-Gutierrez, 2022; Rolandone, 2022). During the interseismic phase, elastic strain builds up around locked portions of the megathrust for tens to hundreds of years (Lay, 2015; Avouac, 2015). This manifests as offshore subsidence and onshore uplift at rates ranging from millimeters to centimeters per year (e.g., Burgette et al., 2009). A fraction of the accumulated energy is eventually released in a matter of seconds to minutes during the coseismic phase (i.e., the earthquake), which involves sudden offshore uplift – a potential source of tsunamis – and offshore subsidence. Additional energy may be released in the following days to months during the postseismic phase as portions of the megathrust experience afterslip, and stresses slowly relax in the asthenosphere (Wang, 2007; Wang et al., 2012). Associated ground displacement rates may then exceed several centimeters per year (Ozawa et al., 1999; Fletcher et al., 2001; Wang et al., 2003; Bürgmann et al., 2005; Métois et al., 2012; Stevens and Avouac, 2015; Klein et al., 2016; Loveless and Meade, 2016; Jolivet et al., 2020; Zhao et al., 2022; Periollat et al., 2022).

Overall, the highly transient deformation pattern that defines a megathrust cycle is well explained by elastic models in which slow interseismic deformation is entirely recoverable. Specifically, the widely used back-slip model of Savage (1983) postulates that interseismic and coseismic deformation perfectly balance each other, and result in no residual deformation other than an increment of slip on the plate interface.

A growing number of observations, however, challenge this basic assumption. Many point out similarities between the spatial pattern of interseismic deformation measured over tens of years and morphological features of the upper plate that have been shaped over hundreds of thousands of years. In Chile, for example, peninsulas seem to

<footnote>*Corresponding author: pierre.dublanchet@minesparis.psl.eu</footnote>
coincide with areas of reduced interseismic locking (Saillard et al., 2017). In the Himalayas, the elevation profile and incision rates of major rivers suggest that they have been subjected to vertical displacements persisting over hundreds of kyrs which mimic current interseismic displacements (Lavé and Avouac, 2001; Meade, 2010; Dal Zilio et al., 2021a).

Overall, this suggests that coseismic and interseismic deformation are not perfectly balanced and instead add up to non-recoverable strain that resembles the interseismic deformation field, and leaves a distinct signature in subduction landscapes (Lavé and Avouac, 2001; Meade, 2010; Cattin and Avouac, 2000; Avouac, 2003). Non-recoverable (i.e., inelastic) deformation is also evident in subduction upper plates either in the form of coseismic fractures (Baker et al., 2013), or diffuse seismicity above the locked megathrust during the interseismic phase of the megathrust cycle (Madella and Ehlers, 2021; Oryan et al., 2024). Whether inelastic deformation accumulates predominantly in a quasi-static manner during interseismic loading or through dynamic coseismic damage is nevertheless still debated (Simpson, 2023).

Measuring the associated deformation rates remains extremely challenging with standard geodetic methods, but can be achieved through geomorphological markers that track cumulative displacements of the upper plate over thousands of seismic cycles ($10^5 – 10^6$ years) (Lavé and Avouac, 2001; Mouslopoulou et al., 2016). Marine terraces, for example, indicate sustained uplift on the order of $0.1 \text{ mm.yr}^{-1}$ along the Chilean coast over the past million year, with faster rates found in areas that experience faster interseismic uplift today (Jolivet et al., 2020). There, the estimated “permanent” uplift rate amounts to 4–8% of the interseismic uplift rate as determined by geodesy. Overall, rates of permanent uplift may vary across settings between $0.1$ and almost $10 \text{ mm.yr}^{-1}$ (Lavé and Avouac, 2001; Mouslopoulou et al., 2016), but are typically at least an order of magnitude slower than interseismic uplift rates measured at the same location (Jolivet et al., 2020; Dal Zilio et al., 2021a). Due to the sparse coverage of geomorphological measurements, the spatial pattern of non-recoverable displacements across a subduction upper plate is poorly known. Some authors have speculated that it could be a small fraction of the interseismic uplift field (Mouslopoulou et al., 2016; Saillard et al., 2017; Meade, 2010), but so far a theoretical framework describing the relationship between permanent and elastic deformation is lacking. One avenue to move this debate forward is to consider the mechanical processes through which interseismic and coseismic deformation may cause increments of non-recoverable strain that accumulate over multiple cycles (Baden et al., 2022; Mallick et al., 2021; Oryan et al., 2024).

How permanent deformation builds up over geological timescales is a question that goes far beyond the specific case of subduction zones where it is well documented: it can be raised for any kind of fault system (King et al., 1988), and studied with models incorporating either purely plastic rheologies (Davis et al., 1983), or visco-elastic-plastic rheologies (Cattin and Avouac, 2000; Ruh and Vergés, 2018; Menant et al., 2020). In these approaches, the entire plate boundary is modeled as a continuous medium with possibly temperature-dependent properties and the effects of gravity, and accounting for surface processes such as erosion-deposition. However, because seismic cycles do not spontaneously emerge in such models, one cannot straightforwardly compare long-term deformation patterns to interseismic or coseismic deformation (Cattin and Avouac, 2000). On the other hand, seismic cycle models spontaneously reproduce the different stages of the earthquake cycle, and the associated deformation pattern. This class of models assumes at least a balance between elasticity of the bulk rock and dynamic friction along a pre-existing interface. In this framework, a plate boundary or a major fault zone can be modeled either as Burridge-Knopoff spring-and-slider systems (Burridge and Knopoff, 1967; Carlson et al., 1994), or as a frictional interface between elastically deformable solids (Rice, 1993; Lapusta et al., 2000). The crucial ingredient allowing the spontaneous occurrence of earthquake cycles is rate-and-state friction (Dieterich, 1979; Ruina, 1983). As long as pure linear elasticity is considered, permanent deformation never accumulates over successive cycles.

Several earthquake cycle models account for the inelastic behavior of the rock mass by assuming visco-elastic or visco-elasto-plastic rheology in addition to rate-and-state friction (Allison and Dunham, 2018, 2021; Dal Zilio et al., 2021b; Barbot, 2018; Lambert and Barbot, 2016). The study of such visco-elastic models has shown how the different features of the earthquake cycle (interseismic or postseismic deformation) can be influenced by the rheological behavior of the crust and upper mantle. Depending on the thermal balance within these layers, the depth, and the distance to the faults, the contribution of viscous strain can be significant in the total observed deformation. Moreover, such approaches reproduce complex sequences of earthquake and aseismic slip observed at the scale of a single fault.

Off-fault plasticity can also be incorporated in the same framework: Erickson et al. (2017) developed a coupling between quasi-dynamic elasticity, rate-and-state friction and a Drucker-Prager yield criterion at the scale of a 2D antiplane planar fault. Again, this approach allowed simulation of the evolution of inelastic deformation in the bulk rock surrounding the fault zone. Simpson (2023) simulated the spontaneous development and activation of faults in a 2D plain strain, fully dynamic, elasto-plastic, Mohr Coulomb domain undergoing slow tectonic loading. In this latter model, increments of plastic strain accumulate both during the interseismic and coseismic stages, over several
earthquake cycles. Mia et al. (2022) recently developed a similar model incorporating full elastodynamics, and was able to reproduce complex sequences of rapid and aseismic slip events. All these approaches however concentrated on the effects of viscoelasticity and off-fault plasticity on the earthquake cycle and the earthquake rupture properties, and generally do not discuss what controls the long-term accumulation of inelastic strain.

Recent attempts have also been made at designing numerical models that incorporate both seismic cycles and the build-up of long-term inelastic deformation (van Dinther et al., 2013; Mallick et al., 2021). Due to their high computational cost, these models do not yet lend themselves to detailed parameter explorations. Here we present a complementary, simpler approach that captures the essential physics of a general fault system producing seismic cycles interacting with off-fault inelastic strain. Our model couples a simplified elasto-plastic rheology with rate-and-state friction in a simple 1-D framework consisting of a thin elastic plate scraping a vertical wall. This simple analog of a fault block subjected to cycles of loading and unloading allows us to derive scaling relations for seismic cycle duration as well as the rates and spatial pattern of inelastic strain build-up, in relation to the plate yield stress and the characteristics of the plate interface (geometry, dynamic friction). Through this simple slider-and-springboard model, we specifically relate the accumulation of inelastic strain to progressive weakening of upper tectonic plates due to loading fatigue that accrues cycle after cycle.

# 2 Slider-and-springboard model

## 2.1 Model setup and plate strength

Our model is illustrated in Figure 1a. It consists of a thin horizontal plate of thickness $H$, length $L$, and density $\rho$ scraping a vertical rigid wall that moves downward at a steady speed $v_0$, representing relative plate motion. We also assume that a normal stress $\sigma = (1 - \lambda) \rho g H / 2$ acts on the plate-wall contact. The $(1 - \lambda)$ term accounts for the effect of fluid pressure, with $\lambda = 0$ and $\lambda \sim 1$ corresponding to an interface that is completely dry and one with near-lithostatic pore fluid pressure, respectively. Friction along the contact causes the plate to bend leading to a deflection profile $w(x)$. $w$ is by definition positive for downward bending. We also define the curvature $\omega(x)$ of the plate as:

$$ \omega = \frac{d^2 w}{dx^2}. \qquad (1) $$

As will be shown in the following, the plate's curvature $\omega$ uniformly increases during loading (down-warping), and decreases during unloading.

This geometry is an idealized representation of a plate boundary, where the thin plate represents the behavior of a fault block accumulating distributed permanent deformation as offset accumulates on the fault. The contact between the plate and the wall can be seen as the fault zone where most of the seismic cycle takes place, and relative plate motion is accommodated in part by successive earthquakes. This geometry could be seen as the most simple improvement to the classical spring-and-slider fault model (Burridge and Knopoff, 1967) that considers a fault normal length-scale through the plate length $L$.

We define the relative slip $\delta$ on the interface as follows:

$$ \delta(t) = v_0 t - w(L, t), \qquad (2) $$

where $t$ is time, and $w(L, t)$ the deflection at the right end of the plate. The frictional force at the interface is positive in the downward direction $-y$. Its magnitude $F$ per unit length along the $z$ direction is given by:

$$ F = f \sigma H, \qquad (3) $$

where $f$ is the dynamic friction coefficient, to be defined in Section 2.4.

In the limit of small deflections, namely $dw/dx \ll 1$, the dominant component of the stress tensor within the plate is the fiber stress $\sigma_{xx}$ (Turcotte and Schubert, 2002). The plate is assumed to deform in an elastic-plastic manner. For simplicity, we adopt the diamond-shaped yield stress envelope proposed by Buck (1988) and shown in Figure 1b. In this model, the plastic yield stress vanishes at the free surfaces of the plate, and its magnitude increases linearly towards the neutral fiber ($y = 0$ in Figure 1a). Introducing $\sigma_Y$, the maximum magnitude of the yield stress at $y = 0$, the plastic yield stress $\sigma_p(y)$ can be written as:

$$ \sigma_p(y) = \pm \sigma_Y k_Y(t) \left( 1 - \frac{2|y|}{H} \right). \qquad (4) $$

$k_Y(t)$ in equation (4) is a factor that accounts for the possibility of a temporal evolution of the yield stress. Here we will either assume that the yield stress remains constant ($k_Y = 1$), or decreases through time because of damage or
fatigue effects (Cerfontaine and Collin, 2017), so that $0 < k_Y \le 1$. When fiber stresses inside the plate fall below the yield stress envelope, the plate behaves elastically, so that changes in $\sigma_{xx}$ are linearly related to changes in curvature $\omega$. When $\sigma_{xx}$ reaches the envelope, it remains equal to $\sigma_p$ as long as $\dot{\omega}$ does not change its sign. Note that because the yield envelope is symmetric, the fiber stress can saturate both during loading or unloading, in compression ($\sigma_{xx} < 0$) or tension ($\sigma_{xx} > 0$).

## 2.2 Linking plate curvature and bending moment

The distribution of fiber stresses across the plate gives rise to a bending moment $M(x)$ given by:

$$M(x) = \int_{-H/2}^{H/2} \sigma_{xx}(x, y)ydy. \qquad (5)$$

Let us first consider the evolution of stress and bending moment within the plate during a loading phase. To this end, we assume that friction at the contact is infinite so that the right end of the plate follows the rigid wall (no relative slip $\delta = 0$). Because its left end is fixed, the plate bends downward and its curvature increases (Figure 1a). A typical profile of fiber stresses inside the plate is shown in Figure 1b, and can be written as:

$$\sigma_{xx} = \begin{cases} E'\omega y & \text{if } |y| < h/2 \\ \text{sign}(y)\sigma_Y k_Y (1 - 2|y|/H) & \text{if } |y| > h/2 \end{cases} \qquad (6)$$

where $h$ is the thickness of the elastic core of the plate (i.e. the plate domain where $\sigma_{xx}(y)$ does not saturate at $\sigma_p(y)$), and $E'$ is the modified Young's modulus of the plate defined from the Young's modulus $E$ and the Poisson ratio $\nu$ as:

$$E' = \frac{E}{1 - \nu^2}. \qquad (7)$$

At $y = h/2$, $\sigma_{xx}$ is continuous so that:

$$h = \frac{H}{1 + E' H \omega / 2 \sigma_Y k_Y} = \frac{H}{1 + R_0 \omega / k_Y}. \qquad (8)$$

In the latter expression, $R_0$ is a characteristic radius of curvature given by:

$$R_0 = \frac{E' H}{2 \sigma_Y}. \qquad (9)$$

From Equation (8), $h$ is close to $H$ for $\omega \ll 1/R_0$, corresponding to a quasi-elastic behavior. Then $h$ decreases with increasing $\omega$, and goes to 0 as $\omega \gg 1/R_0$, which corresponds to a purely plastic behavior. $R_0$ is thus the critical radius of curvature that marks a significant deviation from purely elastic behavior. Equation (6) can be integrated following (5), leading to the following expression of the bending moment:

$$M = D\omega \frac{(1 + R_0\omega/2k_Y)}{(1 + R_0\omega/k_Y)^2}, \qquad (10)$$

where $D$ is the elastic flexural rigidity defined as:

$$D = \frac{E' H^3}{12}. \qquad (11)$$

At small strain ($\omega \ll 1/R_0$), Equation (10) reduces to the purely elastic case: $M = D\omega$. For larger strain, namely when $\omega$ becomes significantly larger than $1/R_0$, the bending moment deviates from the purely elastic solution, it increases more slowly, and plateaus at a maximum value $M_{max}$ as $\omega$ goes to infinity. This maximum possible bending moment (corresponding to the completely saturated stress profile shown in blue in Figure 1b) is given by:

$$M_{max} = \frac{Dk_Y}{2R_0} = \frac{\sigma_Y k_Y H^2}{12}. \qquad (12)$$

Equation (10) is illustrated in Figure 1c for $k_Y = 1$ (black curve), and compared to the purely elastic and the purely plastic case.

As will be shown later, a loading phase cannot last indefinitely because reaching the finite, static frictional strength of the fault (contact at $x = L$) will eventually cause cycles of plate loading and unloading. This implies that Equations (6) and (10) will no longer be valid as soon as the first loading phase ends. The fiber stress profile will then depend on the residual stress profile ($\sigma_{xx}^0$) and the residual curvature ($\omega_0$) attained at the end of the previous loading (or
<figure><img src="image_2.png" /><figcaption>**Figure 1** (a) Model geometry. The elastic-plastic springboard is represented in gray. The black dotted line along *x* is the neutral fiber. The frictional contact (i.e. the plate interface) is represented as a thick black line. (b) Yield stress envelope at position *x* along the plate (gray), along with a typical purely elastic fiber stress profile (dashed red), a purely plastic fiber stress profile (dashed blue), and an elastic-plastic stress profile (black). The yield envelope is inspired by Buck (1988). (c) Bending moment vs. curvature during plate loading for a purely elastic (red), a purely plastic (blue) and the elastic-plastic plate considered here (black). The color code is the same as in panel (b).</figcaption></figure>

unloading) phase. It can be written succinctly as:

$$ \sigma_{xx}(x, y) = \text{sign} \left( [\sigma_{xx}^0 + E'(\omega(x) - \omega_0(x))y] \right) \min \left[ |\sigma_{xx}^0 + E'(\omega(x) - \omega_0(x))y|, |\sigma_p| \right]. \quad (13) $$

The corresponding bending moment can then be obtained by integrating this stress profile according to (5). Details of the calculation are provided in Appendix A. We end up with the following expression for the bending moment:

$$ \mathcal{M} = \mathcal{M}_0 + D(\omega - \omega_0) \frac{(1 + R_0|\omega - \omega_0|/4k_Y)}{(1 + R_0|\omega - \omega_0|/2k_Y)^2}, \quad (14) $$

where $\mathcal{M}_0$ is the residual bending moment, i.e., the bending moment reached at the end of the previous loading or unloading phase. During loading ($\omega > \omega_0$), $\mathcal{M}$ increases above $\mathcal{M}_0$, and it decreases below $\mathcal{M}_0$ during unloading ($\omega < \omega_0$). Note that the change in bending moment (second term on the right-hand side of Equation 14) is anti-symmetric with respect to changes in curvature $\omega - \omega_0$.

These examples of stress evolution during loading and unloading show that the model allows for quasi-static yielding, in particular during loading (or interseismic) phases. However, it does not capture any dynamic coseismic yielding accumulating close to the fault zone, even if unloading is not purely elastic.

Using Equation (14), it is possible to determine the evolution of the plate's bending moment to a given history of
changes in plate curvature. This requires keeping track of individual phases of loading (increasing $\omega$) and unloading (decreasing $\omega$), and the residual moment and curvature at the transitions between these phases. To illustrate this procedure, we consider a simple curvature evolution chosen to mimic the patterns that will be exhibited later on by our complete model. This consists of an initial linear increase in curvature, followed by periodic oscillations superimposed on another linear trend:

$$R_0\omega = \begin{cases} Ct/T & \text{for } t < T \\ A(t-T) + B \sin \frac{2\pi(t-T)}{T} + C & \text{for } t > T \end{cases} \quad (15)$$

In Equation (15), $t$ is time, $T$ the period of imposed curvature oscillations, and $A, B$ and $C$ are arbitrary constants. The resulting bending moment history calculated from (10) and (14) is illustrated in Figure 2 for $B = 0.9$, $C = 1$ and two values of $A = 0$ and $A = 0.1$. Since the plate conserves a memory of its past stress history, the bending moment follows a complex evolution (Figures 2b and 2c). As soon as the first loading phase ends ($t = T$), the same bending moment can be achieved with different values of absolute curvature $\omega$. However, in the case of a purely cyclic loading ($A = 0$), only two curvatures can occur for the same bending moment, and the system remains on a cyclic trajectory in the ($M, \omega$) space (Figure 2b). By contrast, a steady linear increase in curvature ($A > 0$) will lead to steady shifts in the admissible bending moments, which manifests as drifting cycles in Figure 2c.

Of course in practical applications of our model, the spatio-temporal evolution of plate curvature is not known a priori. It must instead be obtained by solving equilibrium equations in the plate, coupled with the dynamic evolution of the frictional resistance on the plate-wall interface.

## 2.3 Moment and force balance within the plate

In order to find the deflection history within the plate under steady motion of the wall, we invoke the force balance and the moment balance derived by Turcotte and Schubert (2002). The moment balance on a plate element between $x$ and $x + dx$ can be written as:

$$M(x + dx) - M(x) + Pdw + Vdx = 0, \quad (16)$$

where $P(x)$ is a horizontal compressive force acting within the plate, and $V(x)$ is the vertical shear force, both acting on a vertical section of the plate at position $x$ and $x + dx$. The force balance projected on the $x$-axis indicates that the horizontal force $P(x)$ is constant along the plate, and thus equal to its value at $x = L$, i.e., $P = \sigma H$. In these conditions, the moment balance (16) becomes (after dividing by $dx$ and taking the limit $dx \to 0$):

$$\frac{dM}{dx} + \sigma H \frac{dw}{dx} + V = 0. \quad (17)$$

The force balance projected on the $y$-axis can be written as:

$$\frac{dV}{dx} = -q(x) = 0, \quad (18)$$

where $q(x)$ describes vertical loads applied on the plate, which we set to zero here. This implies that $V$ does not depend on position $x$. $V$ can thus be eliminated from the moment balance by taking the derivative with respect to $x$, finally leading to:

$$\frac{d^2 M}{dx^2} + \sigma H \frac{d^2 w}{dx^2} = 0. \quad (19)$$

We next consider boundary conditions, beginning at $x = 0$. The plate being pinned on its left edge (Figure 1) ensures that

$$w(0, t) = \frac{dw}{dx}(0, t) = 0. \quad (20)$$

At $x = L$, we first assume that the curvature vanishes (no moment is applied) so that:

$$\omega(L, t) = \frac{d^2 w}{dx^2}(L, t) = 0, \quad (21)$$

and that the shear force $V$ acting on the vertical section of the plate at $x = L$ (i.e., the fault) satisfies a quasi-dynamic balance, so that:

$$V = F + \eta H v. \quad (22)$$
<figure><img src="image_3.png" /><figcaption>**Figure 2** Bending moment history at an arbitrary position *x* resulting from an imposed succession of loading and unloading phases of the plate, for illustrative purposes only. (a) Imposed curvature history defined in Equation (15) with *B* = 0.9, *C* = 1 and two different values of *A*. (b) Bending moment vs. curvature for *A* = 0 (purely cyclic loading). (c) Bending moment vs. curvature for *A* = 0.1 (oscillating and linearly increasing curvature).</figcaption></figure>

In Equation (22), *F* is the frictional force defined in (3) and *v* = $\dot{\delta}$ is the instantaneous slip rate on the fault. Note that, like *F*, *V* is positive in the downward direction. $\eta = \mu/2c_s$ is the damping parameter introduced by Rice (1993), where $\mu$ is the shear modulus of the plate, and $c_s$ the shear wave speed. Enforcing the moment balance Equation (17) at *x* = *L* yields our fourth and final boundary condition:

$$V = -\frac{d\mathcal{M}}{dx}(L,t) - \sigma H \frac{dw}{dx}(L,t) = F + \eta H v. \quad (23)$$
<table><caption>**Table 1** Parameters used for simulations shown in Figures 3 and 4</caption><thead><tr><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td>Young's modulus E = 2μ(1 + ν)</td><td>75 GPa</td></tr><tr><td>Poisson ratio ν</td><td>0.25</td></tr><tr><td>Shear modulus μ</td><td>30 GPa</td></tr><tr><td>Plate material density</td><td>2700 kg.m<sup>-3</sup></td></tr><tr><td>Shear wave speed c<sub>s</sub> = √μ/ρ</td><td>3.33 km.s<sup>-1</sup></td></tr><tr><td>Damping parameter η = μ/2c<sub>s</sub></td><td>4.5 MPa.s.m<sup>-1</sup></td></tr><tr><td>Plate length L</td><td>10 km</td></tr><tr><td>Plate thickness H</td><td>3 km</td></tr><tr><td>Gravitational acceleration g</td><td>9.81 m.s<sup>-2</sup></td></tr><tr><td>Normal stress σ = 0.5ρgH</td><td>39.7 MPa</td></tr><tr><td>Reference friction coefficient f<sub>0</sub></td><td>0.6</td></tr><tr><td>Plate rate v<sub>0</sub></td><td>10<sup>-9</sup> m.s<sup>-1</sup> ≈ 3.15 cm.yr<sup>-1</sup></td></tr><tr><td>Initial slip rate v<sub>i</sub></td><td>10<sup>-30</sup> m.s<sup>-1</sup></td></tr><tr><td>Direct effect parameter a</td><td>0.008</td></tr><tr><td>State evolution parameter b</td><td>0.012</td></tr><tr><td>Critical slip d<sub>c</sub></td><td>0.01 m</td></tr></tbody></table>

## 2.4 Dynamic friction and deflection history

The set of Equations (1) to (5), (13), and (19) to (23) are solved with a finite difference method to obtain the plate's deflection history $w(x, t)$ driven by steady motion of the wall. Equation (23) is analogous to the force balance in a classical slider-and-spring model (Helmstetter and Shaw, 2009), where the spring has been replaced by an elastic-plastic thin plate. The last missing ingredient is the evolution law for the dynamic friction coefficient $f$ contained in $F$ (Equation 3).

We assume that $f$ obeys rate-and-state friction (Dieterich, 1979; Marone, 1998a), so that:

$$f = f_0 + a \ln \frac{v}{v_0} + b \ln \frac{v_0 \theta}{d_c}, \qquad (24)$$

where $v = \dot{\delta}$ is the instantaneous slip rate at the interface, and $\theta$ the time-dependent state variable incorporating the dependence on the past slip history. We assume the state variable evolves according to the ageing law (Ruina, 1983):

$$\frac{d\theta}{dt} = 1 - \frac{v\theta}{d_c}. \qquad (25)$$

In Equations (24) and (25), $f_0$, $a$, $b$ and $d_c$ are respectively a reference constant friction coefficient, the direct effect parameter, the state evolution parameter, and the critical slip necessary for frictional evolution (Marone, 1998a).

Details about the algorithm used to solve our system of equations are presented in Appendix B.1. In doing so we assumed no initial deflection ($w(x, 0) = 0$), an initial slip rate $v_i \ll v_0$ so that the fault is initially locked, and a state variable $\theta_0$ such that the initial shear stress $V = 0$. From Equation (24) and the definition of $V$, the initial state variable $\theta_0$ is given by:

$$\theta_0 = \frac{d_c}{v_0} \left( \frac{v_i}{v_0} \right)^{a/b} \exp \left[ - \frac{\eta v_i}{b\sigma} - \frac{f_0}{b} \right] \qquad (26)$$

The algorithm is based on the non-dimensional form of the governing equations, which are provided in Appendix B.2. Our numerical results include both the slip history along the fault and the deflection profile across the plate, which will have an elastic and an inelastic (plastic) component.

# 3 Elasto-plastic seismic cycles

## 3.1 Constant plate strength

### 3.1.1 Numerical simulation results

As a first step, we assume that the plate's yield stress envelope does not evolve through time, which amounts to setting $k_Y = 1$ in Equation (4). We define a reference case with the parameters listed in Table 1. The deflection history of the plate in this reference case is shown in Figures 3 and 4. Figure 4 also shows the evolution of the fault's slip rate, state variable, and frictional shear stress, as well as the maximum plate deflection $w_m = w(L, t)$ for different values of the peak yield stress $\sigma_Y$. In Figure 4, the purely elastic case corresponds to $\sigma_Y \to \infty$.
<figure><img src="image_4.png" /><figcaption>**Figure 3** (a), (b) Plate deflection for the reference elastic-plastic simulation (parameters of Table 1) with $\sigma_Y = 1.6\sigma^* = 1.52$ GPa. The characteristic stress $\sigma^*$ is defined in Equation (27). The normalized deflection $w/L$ is indicated by the color scale. The red box indicates the time frame detailed in the right panels. Black dashed lines in left panels highlight earthquake occurrence (coseismic unloading). (c), (d): Evolution of elastic core thickness $h$ at three different positions $x$ along the plate. Red dots indicate the times corresponding to panels (e) to (j). (e) to (j): Plastic yielding within the plate at 10%, 90% of interseismic loading phases, and during coseismic unloading. The elastic core with stress below the yield envelope is represented in gray. The black regions indicate a stress at the yield envelope. Vertical red dashed lines indicate locations where elastic core thickness history is represented in panels (c) and (d).</figcaption></figure>

The slow tectonic loading imposed by the motion of the wall is accommodated by the deflection of the plate. At any time, the plate deflection monotonically increases from $x = 0$ to $x = L$ (Figure 3a and b). All the simulations start with a loading phase where deflection steadily increases at the plate rate $v_0$ (Figure 4d), indicating that virtually no relative slip occurs along the fault. The contact between the plate and the wall is stuck, the slip rate is negligible ($v \ll v_0$, Figure 4a), and the state variable increases (Figure 4b). The shear stress on the fault increases at the same time (Figure 4c). This increase is linear for a purely elastic plate, indicating a constant effective stiffness. In elastic-plastic plates ($\sigma_Y < +\infty$), the stressing rate decreases, indicating a progressive reduction of effective stiffness, accompanied by progressive yielding of the top and bottom of the plate, that manifests as a reduction of its elastic core thickness particularly near its left edge (Figure 3c, d, e and f). This initial phase stops when the deflection reaches a threshold, which we call $w_s$ (or $w_e$ in the purely elastic case). At that point, a stick-slip cycle starts, in which deflection slowly increases (the plate bends down) at the loading rate $v_0$ for approximately 1000 years before suddenly decreasing (the plate springs back up). The stick slip oscillation shown in Figures 3 and 4 is our earthquake cycle. The two successive steps correspond to interseismic loading (the fault is stuck, the slip rate at the contact is much slower than $v_0$, the shear stress increases), and coseismic unloading (rapid slip on the order of $m.s^{-1}$ occurs on the fault, associated with a drop in shear stress $F/H$ on the order of 5 MPa). In these simulations, the deflection oscillates around the mean level $w_s$ (or $w_e$). We also observe a periodic yielding of the top and bottom parts of the plate (Figures 3c, d, g, h, i and j): interseismic loading is associated with a reduction of the elastic core of the plate $h$, essentially in its left part. Then $h$ suddenly increases close to $H$ at the onset of the coseismic unloading, and slightly
decreases before the start of the next interseismic period. Coseismic unloading thus occurs almost purely elastically (h remains close to H), even if a very limited coseismic yielding occurs (Figures 3g, j).

For lower values of $\sigma_Y$, however, we do not observe stick slip oscillations and find that our algorithm diverges. In Appendix C, we show that this situation corresponds to a contact that remains stuck, because the frictional strength of the fault is greater than the bulk strength (the plastic limit) of the plate. In these conditions, the maximum bending moment the plate can sustain (12) is not able to balance the moment of the frictional force at the right end of the plate (i.e., on the fault), causing failure in the bulk. The moment of the frictional force being approximately $f_0\sigma HL$, we conclude that stick-slip is only possible when $M_{max}$ exceeds $f_0\sigma HL$, which can be recast as a minimum admissible value for the peak yield stress called $\sigma^*$:

$$ \sigma_Y > \sigma^* = \frac{12f_0\sigma L}{H}. \qquad (27) $$

This condition for stick-slip occurrence is in overall agreement with the numerical results shown in Figure 4.

<figure><img src="image_5.png" /><figcaption>**Figure 4** (a) Slip rate on the fault for a purely elastic plate (blue), and an elastic-plastic plate characterized by two different peak yield stresses $\sigma_Y$ (green and red). The characteristic stress $\sigma^*$ is defined in Equation (27). Other parameters are listed in Table 1. The reference simulation shown in Figure 3 is plotted in red. (b) State variable. (c) Shear stress. (d) Maximum deflection. The heavy black dashed line indicates the loading rate $v_0$. We call $T$ and $T_e$ the stick-slip period (earthquake cycle duration) for the elastic-plastic and purely elastic cases, respectively. $w_s$ and $w_e$ correspond to the mean level of deflection under elastic-plastic and purely elastic rheology respectively.</figcaption></figure>

In the earthquake cycle regime ($\sigma_Y > \sigma^*$), we also note a decrease of the mean level of deflection $w_s$ with increasing yield stress above $\sigma^*$ (Figure 4d). A smaller yield stress close to $\sigma^*$ means that the plate experiences plastic yielding more readily, and can develop larger strains in response to the same load. This manifests as larger mean deflection for weaker plates. This dependence is confirmed in Figure 6a for a wider set of simulations. The mean level of elastic-plastic deflection $w_s$ decreases rapidly towards the purely elastic deflection $w_e$ as the yield stress increases. Another feature shown by our simulations is the change in seismic cycle duration with peak yield stress $\sigma_Y$ (Figure 4d). This dependence is also shown in Figure 6b: elastic-plastic seismic cycles are always longer than purely elastic seismic cycles ($T \ge T_e$). As detailed later, elastic plate accumulates less elastic energy for the same increment of deflection, leading to a slower increase of stress on the fault during interseismic loading, and thus a longer time
is needed to recover the coseismic stress drop. We also find that $T$ decreases rapidly toward $T_e$ as $\sigma_Y$ increases.

Furthermore, when considering greater and greater values of $\sigma_Y$ above $\sigma^*$, we find that the spatial pattern of plate deflection remains qualitatively similar to that shown in Figure 3, albeit with slight differences. This is illustrated in Figure 5, where the interseismic deflection in the reference elastic-plastic case (Table 1 and $\sigma_Y = 1.6\sigma^*$) and in a purely elastic case (table 1 and $\sigma_Y \to +\infty$) are compared. Although the patterns are to first order similar (Figure 5a), the elastic-plastic pattern is characterized by a slightly more pronounced deflection in the middle of the plate, on the order of 0.2% of the total deflection in this case (Figure 5b).

<figure><img src="image_6.png" /><figcaption>**Figure 5** (a) Interseismic deflection profiles for the elastic-plastic reference simulation (solid lines), and the purely elastic case (dots). (b) Inelastic interseismic deflection. Colored profiles indicate the difference between the purely elastic and the elastic-plastic deflection profiles shown in Figure (a) throughout the interseismic phase. Dashed black lines correspond to the approximate solution obtained from Equations (28) and (31).</figcaption></figure>

### 3.1.2 Analytical approximations and scaling relations

Our simulations lead us to conclude that the plate's peak yield stress $\sigma_Y$ controls important features of the earthquake cycle: its duration $T$, the mean level of deflection $w_s$ and the spatial pattern of interseismic deflection. To build further intuition into these controls, we design analytical models and scaling relations that capture key ingredients of our simulations, beginning with the interseismic deflection profile.

As shown in Appendix C, as long as the horizontal loading term $\sigma H dw/dx$ in the moment balance Equation (17) is negligible compared to the shear force $V$, the moment balance can be integrated to yield an approximate closed form solution for the deflection profile $w$, provided the history of the fault's shear force $V(t)$ is known. In the case of interseismic loading, we approximate it by assuming the fault is perfectly locked and the right edge of the plate moves downward at velocity $v_0$. With details given in Appendix C.1, we obtain the following expression for changes in plate deflection ($\Delta w(x, t)$, measured relative to the deflection at the end of the previous coseismic phase that marks $t = 0$):

$$ \Delta w = v_0 t^* \left\{ - \left(\frac{x}{L}\right)^2 - \frac{4}{\alpha(t/t^*)} \frac{x}{L} \sqrt{1 - \alpha(t/t^*)} - \frac{8}{3\alpha(t/t^*)^2} \left[ (1 - \alpha(t/t^*))^{3/2} - (1 - \alpha(t/t^*) + \alpha(t/t^*) \frac{x}{L})^{3/2} \right] \right\}. \quad (28) $$

This equation introduces a characteristic time $t^*$:

$$ t^* = \frac{2\sigma_Y L^2}{E' H v_0} = \frac{L^2}{R_0 v_0}. \quad (29) $$

It also involves a function $\alpha(t/t^*)$ which is defined as:

$$ \alpha(s) = \frac{-8 + 4\sqrt{4 + 27s(s+1)^2}}{9(s+1)^2}. \quad (30) $$

$\alpha$ corresponds here to the moment exerted by the vertical force $V$ at the right end of the plate, normalized by $\sigma_Y H^2 / 4$
(Appendix C.1). As interseismic deflection accumulates, V increases and so does α. With a purely elastic rheology, V and thus α are expected to increase linearly with deflection, i.e., with time under constant loading speed $v_0$. In an elasto-plastic plate, however, a smaller increase in V or α is required to achieve the same amount of deflection. This is because plasticity causes a non-linear response that is all the more significant as s gets close to or greater than 1, that is: as t gets closer to $t^*$. In the limit $t \ll t^*$ (or $s \ll 1$), Equation (30) predicts a linear increase of the normalized moment with time. $t^*$ can thus be understood as the characteristic time needed to observe a significant difference between a purely elastic and an elasto-plastic plate, or, said differently, to generate a curvature of $1/R_0$ if the right end of the plate is steadily moved at speed $v_0$.

Similarly, in a purely elastic plate, the interseismic change in deflection can be approximated as:

$$ \Delta w_e = \frac{v_0 t}{2} \left(\frac{x}{L}\right)^2 \left(3 - \frac{x}{L}\right). \qquad (31) $$

Again, details leading to Equation (31) are provided in Appendix C.1. The difference between the closed form solutions (28) and (31) is in rough agreement with the numerical results as shown in Figure 5. It provides a reasonable approximation for the spatial pattern and evolution of the inelastic component of interseismic deflection in our simulations. We conclude that the accumulation of inelastic deflection during the interseismic phase is essentially controlled by the loading rate $v_0$ and the characteristic time $t^*$. From Equation (29), a small plate yield stress $\sigma_Y$ leads to a small value of $t^*$, and a rapid accumulation of inelastic deflection during the interseismic phase. As the yield stress increases, $t^*$ becomes larger, and the inelastic component builds up more slowly.

Equations (28) and (31) can also be used to derive scaling relations for the mean level of plate deflection $w_s$ and the duration of a seismic cycle $T$. The mean level of deflection can first be assessed by assuming the plate's bending moment when it is deflected at its right end by an amount $w_s$ roughly balances the moment of the shear force $V$ averaged over many earthquake cycles, which is approximately $f_0\sigma H L$. We show in Appendix C.2 that $w_s$ is then given by Equation (28) evaluated at $x = L$, with $\alpha = \sigma^*/\sigma_Y$. Similarly, the mean plate deflection in the elastic case $w_e$ is obtained from (31), bearing in mind that $\alpha$ is a linear function of time in this special case (details about the derivation are provided in Appendix C.2). The ratio between $w_s$ and $w_e$ can then be shown to depend solely on the ratio of $\sigma_Y$ to the characteristic stress $\sigma^*$:

$$ \frac{w_s}{w_e} = -3 \frac{\sigma_Y}{\sigma^*} - 12 \left(\frac{\sigma_Y}{\sigma^*}\right)^2 \sqrt{1 - \frac{\sigma^*}{\sigma_Y}} - 8 \left(\frac{\sigma_Y}{\sigma^*}\right)^3 \left[ \left(1 - \frac{\sigma^*}{\sigma_Y}\right)^{3/2} - 1 \right] = f\left(\frac{\sigma_Y}{\sigma^*}\right). \qquad (32) $$

As shown in Figure 6, our scaling relation (32) captures the decrease of $w_s$ with increasing yield stress $\sigma_Y$. As $\sigma_Y \to +\infty$, the functional $f$ reaches 1 so that $w_s$ converges to $w_e$, as expected.

<figure><img src="image_7.png" /><figcaption>**Figure 6** (a) Mean level of deflection $w_s$ vs. normalized yield stress. $w_e$ is the mean level of deflection expected under purely elastic rheology ($\sigma_Y \to +\infty$). (b) Earthquake cycle duration $T$ vs. normalized yield stress. $T_e$ is the earthquake cycle duration for a purely elastic plate. Colored symbols are simulation results. The black solid line in (a) and (b) is our analytical approximation, which corresponds to the functional $f$ defined in Equations (32) and (33). We explored here different values of $\sigma_Y$, of $f_0$ (0.4 and 0.6), of $H$ (1 and 3 km), and of $b$ (0.01 and 0.012). Other parameters are listed in Table 1.</figcaption></figure>

We finally turn to estimating the duration $T$ of elastic-plastic seismic cycles, which can be done with the help of Equation (28). At the end of an interseismic phase, $\Delta w$ evaluated at $x = L$ is simply $v_0T$, and $\alpha$ corresponds to the
normalized moment of V assuming V is on the order of the earthquake stress drop Δσ multiplied by H. This allows us to derive an expression for T, and for T<sub>e</sub> in the purely elastic case (using Equation 31). We detail this reasoning in Appendix C.3, and show that:

$$ \frac{T}{T_e} = f\left(\frac{\sigma_Y}{\Delta\sigma^*}\right), \qquad (33) $$

where f is the function introduced in Equation (32), and Δσ* is the characteristic stress drop given by:

$$ \Delta\sigma^* = 6\Delta\sigma\frac{L}{H} = 6(b-a)\sigma \ln\frac{v_{sis}}{v_0}\frac{L}{H}. \qquad (34) $$

Δσ in (34) is the coseismic stress drop associated with an earthquake on the fault. The stress drop is controlled by the frictional strength, and following Rice and Tse (1986), we write it as $(b-a)\sigma \ln v_{sis}/v_0$, $v_{sis}$ being an order of magnitude of the coseismic slip rate (typically between 1 cm.s⁻¹ and 1 m.s⁻¹). The scaling of Rice and Tse (1986) predicts a stress drop of 5 MPa, which is in agreement with the stress drop produced by our simulations (Figure 4c).

Again, scaling relation (33) is well supported by our simulation results (Figure 6). This demonstrates that the duration of earthquake cycles with an elasto-plastic plate is essentially controlled by the stress ratio $\sigma_Y/\Delta\sigma^*$.

In the simulations presented so far, before stick-slip cycles can occur, the plate must flex down to a point where the shear force on its right edge matches the static friction of the fault. In doing so, inelastic strain develops near the top and bottom of the plate, in areas with large curvatures (Figure 3). The plate, however, retains an elastic core, which can store and release elastic energy during the interseismic and coseismic phases. It essentially behaves as a damaged elastic spring, which is slightly less stiff that a perfectly elastic plate. Interestingly, while non-recoverable deformation has developed in the plate, there is no net accumulation of “new” non-recoverable deformation over successive cycles, as illustrated in Figures 3c to j and in Figure 4d.

As shown in Appendix C, to first order, the bending moment balances the moment exerted by the frictional force on the fault. Then, the maximum and minimum possible bending moments are essentially controlled by the static frictional strength and the residual frictional strength (after a coseismic strength drop), respectively. Both of these quantities do not change from one seismic cycle to the next. The minimum and maximum bending moments therefore do not change either. This is similar to the situation illustrated in the moment vs. curvature diagram of Figure 2b: because the maximum and minimum bending moments stay the same, the maximum and minimum curvature also stay the same, and the system's trajectory remains fixed in moment-curvature space.

In light of this result, a net accumulation of inelastic strain over successive cycles requires changes in either the frictional properties of the plate interface, or the yield strength/elastic properties of the plate through time. Next we focus on the effect of damage accumulation resulting from cyclic loading and unloading of the upper crust, i.e., cyclic fatigue. This mechanism could plausibly result in a steady weakening of the plate's yield strength envelope, as well as its elastic properties, while frictional strength remains constant.

## 3.2 Decreasing plate strength

### 3.2.1 Numerical simulations with cyclic fatigue

Cyclic fatigue manifests as a decrease in the strength of a material caused by periodic loading and unloading (Cerfontaine and Collin, 2017). One underlying mechanism can be the slow, subcritical growth of small-scale cracks that lengthen incrementally during the loading phase of each cycle (Scholz, 1972; Brantut et al., 2013; Cerfontaine and Collin, 2017; Scholz, 2002). Long-term changes to the internal state of rocks (increased damage) then result in softening of their mechanical properties (Bhat et al., 2012). In our model, the plate experiences cycles of interseismic loading and coseismic unloading, controlled by dynamic friction along the fault. The plate may therefore experience cyclic fatigue over successive earthquake cycles. In the following, we will assume that this effect only leads to a decrease of the peak yield stress $\sigma_Y$ with time. We account for this through the $k_Y$ parameter (see Equation 4 for a definition of $k_Y$), which is now allowed to decrease from one cycle to the next. We specifically adopt the functional form of Cerfontaine and Collin (2017), which after n cycles gives:

$$ k_Y = 1 - \phi \log n. \qquad (35) $$

The non-dimensional constant $\phi$ in the above equation will hereafter be called the fatigue parameter. From Cerfontaine and Collin (2017), $\phi$ is of the order of 0.02 to 0.03 for rock materials. Here we explore values ranging from 0.01 to 0.1. Our numerical methodology otherwise remains unchanged. The maximum deflection history obtained for the reference case (parameters in Table 1) with an intact yield stress $\sigma_Y = 2\sigma^*$ and for different values of $\phi$ is shown in Figure 7. Accounting for cyclic fatigue leads to a slow accumulation of permanent deflection on top of the
stick-slip (earthquake) cycles. In this case, coseismic unloading does not fully compensate interseismic deflection. This effect becomes more pronounced as the fatigue parameter $\phi$ increases. It should be noted that for a given $\phi$ the rate of permanent deflection accumulation is fastest during the very first cycles, and slows down with time.

<figure><img src="image_8.png" /><figcaption>**Figure 7** Maximum deflection $w_m = w(L, t)$ of an elasto-plastic plate undergoing cyclic fatigue, for different values of the fatigue parameter $\phi$. Solid lines are the numerical solutions, dashed lines are the approximate solutions from Equation (36). The initial loading phase (from zero deflection to the first earthquake) is truncated, and the different curves are shifted horizontally to the same origin ($t = 0$ set to the time of the first earthquake) and vertically so that $w_m = 0$ at $t = 0$. These simulations use $\sigma_Y = 2\sigma^*$ along with the parameters listed in Table 1.</figcaption></figure>

### 3.2.2 Scaling relations for the accumulation of non-recoverable deformation

The permanent deflection accrued here over successive earthquake cycles can be understood as a direct consequence of the relationship between plate deflection and yield stress presented in Section 3.1.2. First, as the yield stress progressively decreases due to fatigue, the ratio $\sigma_Y/\sigma^*$ decreases, and the plate is driven towards a new mean level of deflection $w_s$ that is larger than the preceding one (see Figure 6a). Then, since the ratio $\sigma_Y/\Delta\sigma^*$ also decreases, we expect that fatigue progressively leads to an increase of the seismic cycle duration (Figure 6b). The combination of these two effects is what drives the permanent deflection in our model. In Appendix C.5 and Figure 7, we show that the first effect (increase of the mean level of deflection) is enough to explain the rate of permanent deflection observed in the simulations. Namely, the envelope of maximum deflection ($w_p(t)$, measured at $x = L$) is well approximated by the following equation (directly derived from Equation 28):

$$w_p(t) = v_0 t^* \psi \left[ \frac{\sigma_Y}{\sigma^*}, \phi, n(t) \right]. \qquad (36)$$

In Equation (36), $t^*$ is the characteristic time defined in Equation (29), $n(t)$ is the total number of earthquake cycles at time $t$ ($t = 0$ here corresponds to the first earthquake on the fault, when permanent deflection starts to accumulate, as shown in Figure 7), and $\psi$ is the functional defined in Appendix C.5 (Equation 92). $n(t)$ can be approximated as:

$$n(t) = \frac{t}{T_e} = \frac{\sigma_Y}{t^* \Delta\sigma^*}. \qquad (37)$$

where $T_e$ is the purely elastic earthquake cycle duration. We show in Appendix C.5 that $T_e$ writes $t^*\Delta\sigma^*/\sigma_Y$, leading to Equation (37). A more accurate value for $n$ would have been obtained if we had accounted for the change in cycle duration $T$ with time, but as shown in the preceding section, this change only represents a small fraction of $T_e$ (Figure
6b) that is neglected here, so that $T \simeq T_e$. $w_p$ from Equation (36) is represented with dashed lines in Figure 7. It is in good agreement with the numerical solution. The permanent inelastic deflection is therefore controlled by the stress ratios $\sigma_Y/\sigma^*$ and $\sigma_Y/\Delta\sigma^*$, as well as the fatigue parameter $\phi$, the long-term plate rate $v_0$ and the characteristic time $t^*$.

To compare our models to natural systems, it is useful to define proxies for the inelastic strain rate, e.g., a rate of permanent deflection accumulation. Since this rate is not constant through time (Figure 7), we can only define an average rate $R_p$ between the $k^{th}$ and $n^{th}$ earthquake cycle ($n > k$). Figure 8a shows this rate as a function of the fatigue parameter for a broad set of simulations. We find $R_p$ ranging from less than 0.1 mm.yr$^{-1}$ to about 1.4 mm.yr$^{-1}$ for the range of parameters considered here. In each case, $R_p$ increases with the fatigue parameter $\phi$, and it is also larger as $\sigma_Y$ gets closer to the minimum yield stress required for stick-slip occurrence $\sigma^*$. Frictional parameters $f_0$ and $b-a$ also influence the rate of permanent deflection, essentially because they contribute to $\sigma^*$ and $\Delta\sigma^*$ (Equations (27) and (34)). Using Equation (36), we show in Appendix C.5 that $R_p$ between cycle k and cycle n can be approximated as:

$$R_p = \frac{3v_0}{(n-k)} \frac{\sigma_Y}{\Delta\sigma^*} \chi \left[ \frac{\sigma_Y}{\sigma^*}, \phi, n, k \right], \quad (38)$$

where $\chi$ is the functional defined in Equation (96). Equation (38) is in overall agreement with our numerical results (Figure 8b). Again, the rate of permanent deflection increases as the fatigue parameter $\phi$ increases, and as $\sigma_Y$ gets closer to $\sigma^*$. It also linearly increases with the stress ratio $\sigma_Y/\Delta\sigma^*$.

The fatigue effect is the cause of permanent strain accumulation over successive cycles, through progressive weakening of the bulk material. At some point, however, the weakened yield stress $k_Y\sigma_Y$ eventually reaches $\sigma^*$, the minimum value necessary for the occurrence of stick-slip on the fault. At that point, the earthquake cycle is replaced by bulk failure within the plate. The system therefore has a limited lifetime: the maximum number of earthquake cycles $n_{max}$ that can be accumulated before bulk failure is attained assuming $k_Y(n_{max})\sigma_Y = \sigma^*$, which yields:

$$n_{max} = 10^{(1-\sigma^*/\sigma_Y)/\phi}. \quad (39)$$

<figure><img src="image_9.png" /><figcaption>**Figure 8** (a) Rate of permanent (inelastic) deflection accumulation $R_p$ vs. fatigue parameter $\phi$. Symbols are numerical estimates of average $R_p$ between cycles $k = 10$ and $n = 15$. The shape of the symbols indicate different frictional parameters ($f_0$ and $a-b$), the colors refer to different values of the peak yield stress $\sigma_Y$. $\sigma^*$ is the characteristic stress defined in Equation (27). (b) Same as (a) for normalized $R_p$. $\Delta\sigma^*$ is the characteristic stress defined in equation (34). Colored symbols are the same numerical results as in panel (a). The black dashed lines correspond to the approximate solutions obtained from Equation (38), for the different $\sigma_Y/\sigma^*$ considered here.</figcaption></figure>

The spatial distribution of permanent deflection is represented in Figure 9a for a particular simulation with $f_0 = 0.4$, $b-a = 0.002$, $\sigma_Y = 2\sigma^*$ and $\phi = 0.05$ (the other parameters are listed in Table 1). The permanent deflection is
calculated here as the difference between the residual deflection at the end of earthquake cycle *n* and the residual deflection at the end of the very first earthquake cycle. The permanent deflection pattern is also compared to the interseismic deflection profile and the coseismic deflection profile. Interseismic and coseismic deflections clearly do not fully compensate each other here, leading to a small increase in deflection at each cycle, reaching about 14% of a typical interseismic deflection after 35 cycles (0.4% of the total interseismic deflection per earthquake cycle, Figure 9a). When normalizing the profiles by the maximum deflection (Figure 9b), we note that the shape of the permanent deflection is different from the interseismic deflection: the curvature of the permanent deflection is more pronounced near *x* = 0 than for the interseismic deflection. The coseismic deflection is on the other hand about the same magnitude as the interseismic one (with opposite sign).

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

**Figure 9** (a) Deflection profiles under cyclic fatigue. The permanent component is evaluated at the time indicated by the green circles in panel (c), assuming zero permanent deflection at the end of the first earthquake. The interseismic and coseismic components indicate the deflection accumulated during the blue and red intervals in panel (c). (b) Same profiles as in (a), but normalized by the maximum interseismic deflection (taken at *x* = *L*). (c) Maximum deflection history *w*<sub>m</sub> (black line). The green circles indicate when permanent deflection is evaluated in panel (a). The blue and red lines indicate the interseismic and coseismic phases used to determine the interseismic and coseismic deflection profiles in panels (a) and (b). The blue circle is the end of the interseismic phase, the red circle marks the end of the coseismic phase. The results presented in this figure were obtained with *f*<sub>0</sub> = 0.4, *b* − *a* = 0.002, *σ*<sub>Y</sub> = 2*σ*<sup>*</sup> and a fatigue parameter *φ* = 0.05. The other parameters are listed in Table 1. Note that in (a) and (b), we have reverted the sign of the deflection to represent the actual shape of the plate (e.g., downward bending during the interseismic phase).

# 4 Discussion

## 4.1 Impact of elasto-plasticity on the earthquake cycle

### 4.1.1 Occurrence of stick slip and analogy with the classical slider-and-spring system

In our analysis, when $\sigma_Y$ exceeds $\sigma^*$, the system enters a stick-slip regime. The classical rate-and-state spring-block model (Ruina, 1983; Rice and Tse, 1986; Helmstetter and Shaw, 2009) generally produces stick-slip oscillations under
velocity-weakening conditions ($a - b < 0$) when the stiffness of the spring $k$ is lower than the critical value $k_c = (b-a)\sigma/d_c$. The purely elastic plate is equivalent to a spring-block model in which the spring stiffness is given by $k_e = 1/H\partial V/\partial w(L) = 3E'H^2/8L^3$ (Appendix C.4). In the range of parameters considered here, this purely elastic stiffness is always much smaller than $k_c$. Since the effective stiffness of an elasto-plastic plate $k$ typically decreases from $k_e$ as curvature develops in the plate (Appendix C.4), we also have $k < k_c$. It is therefore not surprising to obtain stick-slip oscillations in our model under velocity-weakening conditions.

The absence of stick slip for $\sigma_Y < \sigma^*$ is a feature that cannot be reproduced by a purely elastic spring-and-slider system, but arises as soon as a finite yield stress is assumed for the spring (Mia et al., 2023). For the elasto-plastic spring-and-slider system, the bulk failure regime corresponds to the locked state described by Mia et al. (2023). The transition to the locked regime occurs either at constant stiffness if the yield stress is too small, or at constant yield stress with a reduction of the stiffness (Mia et al., 2023). This behavior is analog to our springboard model, bearing in mind that the effective stiffness of the plate decreases as $L/H$ and thus $\sigma^*$ increases. The transition therefore occurs either if $\sigma_Y$ is decreased at constant $\sigma^*$ (or effective stiffness), or if $\sigma^*$ increases (effective stiffness decreases) while keeping $\sigma_Y$ constant.

Building on this analogy with the spring-block model, we expect to generate stable aseismic slip on the fault under velocity-strengthening conditions ($a - b > 0$), or if the geometric ratio $H^2/L^3$ that appears in the above expression of $k_e$ increases. It should be noted that only increasing the thickness $H$ of the plate would result in a linear increase in the average normal stress acting on the plate interface ($\sigma \sim \rho gH$). This would therefore cause $k_c$ to increase. That said, because $k$ scales as $H^2$, we expect that $k$ would increase faster than $k_c$ as $H$ increases. If the plate can be thick enough that $k > k_c$, the plate would first bend towards an equilibrium configuration where it can slip stably at rate $v_0$. As is the case in the stick-slip regime, no permanent accumulation should accumulate past the initial bending phase, unless $R_\sigma = \sigma_Y/\sigma^*$ somehow decreases through time.

Although our model shares many analogies with the elasto-plastic spring-and-slider system, it is also a significant improvement of this classical approach owing to the possibility of computing a spatial pattern of deflection, in the fault normal direction.

### 4.1.2 Elasto-plasticity with constant yield stress

If the plate can retain its yield strength through time (constant $\sigma_Y$), the effects of its elasto-plastic rheology on earthquake cycles are fully determined by a characteristic time $t^*$ (Equation 29, the time past which inelastic deformation becomes significant) as well as two dimensionless stress ratios: $R_\sigma = \sigma_Y/\sigma^*$ and $R_{\Delta\sigma} = \sigma_Y/\Delta\sigma^*$. The first is the ratio of the peak yield stress $\sigma_Y$ to the minimum yield stress required for the occurrence of stick-slip cycles $\sigma^*$ (Equation 27). This ratio can be recast as a moment ratio: $M_{max}/M_f$, $M_{max}$ being the maximum bending moment the plate can sustain (Equation 12), and $M_f = f_0\sigma H L$ the moment exerted by the frictional force (of the order of $f_0\sigma H$) at the right edge of the plate (fault). Elasto-plastic effects become significant as the moment of the frictional force approaches the maximum bending moment of the plate. For example, values of $R_\sigma$ closer to (but above) 1 lead to plate oscillations around a greater mean deflection (Figures 4d and 6a), with wider portions of the plate being brought to plastic yielding (Figure 3c-j). Plasticity also has a very subtle effect on the spatial pattern of interseismic deflection. Figure 5 for example shows that interseismic deflection profiles have a slightly different shape in an elasto-plastic vs. a purely elastic plate. The difference is however very small (less than 0.2%), suggesting that treating off-fault deformation measured at the time scale of one seismic cycle or less as entirely elastic, as is common practice in seismo-geodesy, is valid to first order. Finally, if the moment of the frictional force exceeds the plate's maximum bending moment ($R_\sigma \le 1$), the bulk of the plate fails and seismic cycles cannot occur on the plate interface.

Earthquake cycle models accounting for elastic-plastic bulk in 2D also show that the fault slip regime (locked, slow slip or earthquakes) is to some extent controlled by the ratio between bulk yield stress and fault strength (Erickson et al., 2017; Mia et al., 2022, 2023), which could be one interpretation of the ratio $R_\sigma$. A straightforward extension of the present study would be to assess whether $R_\sigma$ ratio is generalizable in a 2D or 3D elasto-plastic continuum with a finite fault.

The second stress ratio $R_{\Delta\sigma}$ primarily influences the duration of individual seismic cycles, provided $R_\sigma > 1$. As before, $R_{\Delta\sigma}$ can be recast as $M_{max}/\Delta M_f$. $\Delta M_f$ denotes the coseismic change in the moment exerted by the frictional force, which is primarily driven by the reduction in fault friction from its static to dynamic value. Again, the lengthening of the seismic cycle due to elasto-plasticity effects become significant as $\Delta M_f$ approaches $M_{max}$, i.e., as $R_{\Delta\sigma}$ approaches 1 (Figures 4d and 6b). We however note that $R_{\Delta\sigma}$ should be significantly larger than unity, essentially because $\Delta M_f < M_f$ or, said differently, because the change in friction coefficient during the earthquake cycle is a small fraction of the absolute friction coefficient ($\Delta f < f_0$). In the limit of an infinitely strong plate (very large $\sigma_Y$), the maximum bending moment the plate can sustain is infinite, and the plate behaves elastically.
### 4.1.3 Elasto-plasticity with decreasing yield stress

If the plate's yield stress envelope and elastic moduli as well as the fault's frictional parameters remain constant through time, our model cannot produce a net accumulation of non-recoverable deflection from cycle to cycle. Instead, it produces elastic oscillations around a partially yielded state (Figure 3d). It could be argued, however, that the effective strength of the plate should be different during a loading and an unloading phase. Coseismic deformation is fast (seconds to minutes), and much shorter than the characteristic deformation time needed to activate viscous-like inelastic processes (e.g., the Maxwell time of the lithosphere). In this case, the only coseismic plastic deformation would then be related to dynamic damage creation, which is generally limited to the vicinity of the fault (Rodriguez Padilla et al., 2022), or to the earth surface (Baker et al., 2013). By contrast, an interseismic phase lasting > 100 yrs would be long enough to activate viscous-like inelastic deformation mechanisms in the fault blocks, characterized by an effective viscosity below ~ $10^{20}$ Pa.s (e.g., with a Maxwell time below 100 yrs). The difference between coseismic and interseismic behavior regarding viscous-like inelastic deformation could be captured by a highly asymmetric yield strength envelope instead of the diamond-shaped one we have assumed (Figure 1b). That said, even with an asymmetric envelope, the only way to return to the pre-earthquake frictional strength after a coseismic phase is to return to the deflection level at the end of the previous interseismic stage. Thus, we do not expect that an asymmetric envelope would enable any more build-up of permanent strain than a symmetric envelope.

In order to get permanent deformation that accrues over successive earthquake cycles in our model, it is necessary to alter the fault's frictional properties, or the plate's elastic parameters or yield stress envelope through time.

It is well known from field studies and laboratory experiments that the mechanical properties of rocks can weaken during earthquakes, as a consequence of coseismic damage (Thomas and Bhat, 2018). This effect, however, tends to occur in close proximity to the fault zone, where coseismic fracturing is most intense (Rodriguez Padilla et al., 2022). Wider extents of lithosphere may also undergo progressive damage at slow strain rates, as a result of distributed brittle creep. This process results from the sub-critical lengthening of small-scale cracks (e.g., Atkinson, 1984), made possible by slip on small frictional defects such as grain boundaries (Brantut et al., 2013). Macroscopically, it manifests as a change in elastic properties, and brings the rock closer to failure (Bhat et al., 2012). Under periodic loading, brittle creep translates into cyclic fatigue, i.e., a progressive decrease in yield strength (Cerfontaine and Collin, 2017). Here, we restricted our investigations to the case of cyclic fatigue. Other sources of damage would however need to be considered. We leave the study of their contribution for future work.

Considering fatigue or damage effects requires additional parameters to fully characterize inelastic deformation. In our case, we introduced the fatigue parameter $\phi$, which is defined empirically as the decrease of the macroscopic yield strength per loading cycle, which we assimilated here to the earthquake cycle. From experimental studies, $\phi$ is on the order of a few percent (Cerfontaine and Collin, 2017). To account for the fact that experiments are typically conducted on small samples at strain rates much larger than tectonic rates, we explored the effect of a wider range of $\phi$ (1 – 10%). We find that greater values of $\phi$ promote faster accumulation of non-recoverable deflection (Figure 7). Large values of $\phi$ however promote a shorter overall lifespan for the system, as the plate strength will decrease below $\sigma^*$ more quickly, resulting in bulk failure and the cessation of earthquake cycles after $n_{max}$ cycles (Equation 39). Interestingly, with the range of parameter values investigated here, our simulations produce rates of permanent plate deflection between 0.1 and 1.5 mm.yr$^{-1}$, which is a reasonable order of magnitude (Jolivet et al., 2020). We further detail in the next section the implications of our model concerning the rate of inelastic deformation within natural fault zones.

## 4.2 Mechanical control on the accumulation of permanent deformation

Field estimates of non-recoverable deformation rates typically rely on the uplift of datable geomorphic markers such as marine terraces or paleoshorelines (Melnick, 2016; Mouslopoulou et al., 2016). They may also rely on estimates of fluvial incision rates under the assumption that they perfectly compensate rock uplift rates (Lavé and Avouac, 2001; Meade, 2010). Depending on the geodynamic context and the erosion model considered, permanent deformation rates roughly vary between 4 and 8% of the rate of interseismic deformation, which is between 0.1 and 5 mm.yr$^{-1}$ (Melnick, 2016; Jolivet et al., 2020).

Because our model oversimplifies the geometry of a natural fault zone, the comparison with field estimates of non-recoverable deformation rates is not straightforward. We must first assume that the ratio of permanent to interseismic deflection measured at the plate-wall contact ($R_p/v_0$) is representative of the permanent to interseismic deflection ratio at any point in the plate. The second difficulty in comparing our simulations to natural systems is to choose the proper observation window, that has to be long enough to capture the slow accumulation of permanent deformation. ~ 100 kyrs is thus a representative time scale for the geomorphic processes that shape the landscape as a result of non-recoverable uplift. This amounts to hundreds of seismic cycles, each lasting hundreds of years
(Figure 9c). However, because cyclic fatigue first causes rapid weakening of an intact plate followed by slower and slower weakening, the rate of permanent deflection will not be the same if it is measured early on in the evolution of the upper plate (e.g., during its very first 100 seismic cycles: Scenario A for a “young plate”), or late, in an already weakened plate (e.g., between the 900th and 1000th cycle: Scenario B for an “old plate”). Our model predictions of $R_p/v_0$ for the two different scenarios are represented in Figure 10a (Scenario A) and 10b (Scenario B). We assumed here for simplicity a constant ratio $\sigma^*/\Delta\sigma^* = 2f_0/\Delta f \simeq 20$, as suggested by Equations (27), (34).

<figure><img src="image_11.png" /><figcaption>**Figure 10** Mechanical controls on the rate of inelastic deflection accumulation $R_p/v_0$. The color scale corresponds to the ratio between the inelastic deflection accumulation rate $R_p$ and the long-term loading (plate) rate $v_0$, for a range of plausible values for $\phi$ and intact stress ratio $\sigma_Y/\sigma^*$. (a) $R_p$ is averaged over the first hundred seismic cycles of each simulation ($k = 1, n = 100$), or (b) over a hundred cycles occurring later on in the plate's evolution: between cycle $k = 900$ and $n = 1000$. Solid black lines mark the limit that separate the earthquake cycle regime, in which the plate can experience stick-slip, from the regime in which the plate fails under a lower stress than needed to slip on the interface (bulk failure). Each black line assumes a different value of the plate lifespan $n_{max}$ (Equation 39). The dashed red contours outline a plausible range of $R_p/v_0$ for natural fault systems as estimated from field studies (Jolivet et al., 2020).</figcaption></figure>

In both scenarios, the ratio between permanent and interseismic rates ($R_p/v_0$) increases as the fatigue parameter increases and the intact yield stress decreases. In each panel of Figure 10, we indicate with solid black lines the limit between the earthquake cycle regime and the bulk plate failure regime, for two different hypotheses on the lifetime of the plate ($n_{max}$ as defined in Equation (39)) equal to 100 or 1000 cycles). As an example, consider Scenario A for a young plate (Figure 10a). Let us assume that the plate has a finite lifespan $n_{max} = 100$ cycles, meaning that material fatigue will have decreased the plate's strength from its intact value $\sigma_Y$ to $\sigma^*$ (the minimum stress required for stick-slip to occur) after 100 seismic cycles. This implies that earthquakes will cease to happen as bulk failure of the plate replaces stick-slip along the wall. In our example, if the plate had an intact stress ratio $\sigma_Y/\sigma^*$ of 1.5 and a fatigue parameter of 0.3, earthquakes would cease after fewer than 100 cycles. The white field in Figure 10a thus represents a region of parameter space that fails to produce at least 100 earthquakes on the plate interface. The colored field, by contrast, allows at least 100 earthquakes. In those conditions, values of $R_p/v_0$ compatible with the observations outline a narrow region of parameter space indicated with red dashed contours.

The situation analogous to bulk failure of the plate and cessation of seismicity at the main fault could be that damage accumulation localizes in a new fault zone within the plate, leaving the first abandoned. Our model does not allow for the creation of new faults, and thus the bulk failure regime has to be seen as a limit of our model. Also, for that reason, $n_{max}$ can be understood as a proxy for the time over which the current plate interface remains the dominant, seismically-active structure in the system. In natural contexts, this time is related to the very long-term structural evolution of the fault zone, and may be well above several Myrs (1000s of cycles). If we assume $n_{max} = 1000$ cycles, then plausible combinations of stress ratios and fatigue parameters must fall within the dashed red contours and above the $n_{max} = 1000$ curve, which is fairly restrictive. In that case, the fatigue parameter cannot exceed ~0.1, and $\sigma_Y/\sigma^*$ must be less than ~1.5. A trade-off however remains between a strong intact plate (high $\sigma_Y/\sigma^*$) experiencing rapid weakening (high $\phi$) and a weaker intact plate experiencing slow weakening.

In Scenario B, where the permanent deflection rate is measured later, between the 900th and 1000th cycles (Figure 10b), the weakening rate has already had time to slow down significantly. Thus, to produce the rates measured in nature, we need greater values of $\phi$ (faster weakening) than in Scenario A (Figure 10a). We find, however, that in this configuration, the overriding plate would be very close to bulk failure at the end of our measurement window ($n_{max} = 1000$). This illustrates another trade-off between the weakening rate (which decreases through time) and the time interval over which permanent deflection is estimated.
In spite of these trade-offs, Figure 10a shows that with reasonable parameter values, our toy model can produce the right order of magnitude for the ratio of permanent to interseismic deformation. Larger ratios correspond to the lowest yield stress, and/or the largest fatigue parameter. In other words, achieving a rapid build-up of inelastic strain requires a bulk strength (yield stress) that is (or rapidly becomes) close to the frictional strength of major faults.

## 4.3 Future directions

Our slider-and-springboard model oversimplifies fault zones by design, to retain only the essential physics of reversible and irreversible deformation and build intuition into what controls the behavior of these complex systems. Here we nonetheless identify key areas of possible model improvements.

The cyclic fatigue model we have used to describe the progressive decrease of the plate yield strength is entirely empirical, and does not capture any specific mechanism of damage accumulation. The values of the $\phi$ parameter, in particular, rely on experiments performed on small samples, at loading frequencies (< 1 day) and strain rates much larger than those associated with seismic cycles (100s of yrs). As noted by Cerfontaine and Collin (2017), the loading period can be a strong influence on the build-up of material fatigue. This feature might eventually explain why we must invoke values of $\phi$ that are typically an order of magnitude larger than those determined experimentally to obtain plausible rates of permanent deflection. Another explanation could be that we have not considered any kind of dynamic coseismic damage, that would enhance the plate weakening rate. This issue could be mitigated if instead of a quasi-static system with ad-hoc cyclic fatigue model, we considered a dynamic micro-mechanical model for the accumulation of brittle damage in the plate (Bhat et al., 2012). In this class of models, the loss of bulk strength would be indexed on a damage variable that describes the evolution of the rock's internal state, e.g., as pre-existing defects lengthen when loaded steadily above a certain threshold (but below the peak strength of the material: Atkinson, 1984). Such descriptions would also allow us to self-consistently account for changes in the plate's elastic moduli, which will surely affect its deflection pattern.

Since coseismic dynamic damage originating from the fault was neglected, we also assumed that the yield stress only depends on time and not on the distance from the fault. Coseismic stress redistribution could however generate non-uniform damage patterns, as shown by continuum based models (Okubo et al., 2019; Mia et al., 2023; Simpson, 2023). Elasto-plastic earthquake cycle models (Erickson et al., 2017; Mia et al., 2022, 2023) in particular predict a decrease of the plastic strain with distance to the fault. In our model, plastic strain preferentially accumulates at $x = 0$, where the plate is pinned leading to maximum curvature. This limitation will have to be resolved in future studies.

So far we have only considered the possibility of strength-loss of the plate (cyclic fatigue decrease of yield stress), neglecting the possibility of strength recovery. Strengthening of the plate could occur as a consequence of thermal, or chemical crack healing processes (Kirby, 1984; Marone, 1998b; Kanagawa et al., 2000; Tenthorey et al., 2003). The competition between damage and healing would then control the accumulation rate of non-recoverable deflection, which would be maximized if damage has a dominant effect or if the timescale for healing processes is much longer than typical earthquake cycle duration $T$. Here again, this would require to use a proper micromechanical approach to describe healing. Note however that frictional healing on the main fault (plate-wall contact) is to some extent accounted for by the rate-and-state friction law, where the state variable $\theta$ and thus the friction coefficient increases during the interseismic stage (figure 4b).

In comparing our numerical results with our analytical approximations, we have shown that the horizontal force term $\sigma H d^2w/dx^2$ can be neglected to first order in the moment balance Equation (19). This is because plate strains (or gradients of deflection $dw/dx$) typically remain much smaller than the friction coefficient $f$, such that the vertical shear force $V \sim f\sigma H$ dominates the compressive force. This condition originates from the assumption of a frictional contact, which explicitly relates the horizontal and vertical forces that act on the board. $\sigma Hdw/dx$ being small, negligible buckling occurs in our system. Buckling could nevertheless be significant in a different geometry, and should not be ruled out from the present study. Buckling induced by long-term shortening at a compressive plate boundary clearly participates in the accumulation of permanent deflection (e.g., folding in an accretionary wedge). Our toy model could prove useful in understanding how this potential source of permanent strain interacts with inelastic deformation specifically driven by the earthquake cycle.

Finally, in our balance of forces acting on the plate, we have neglected for simplicity any buoyancy effect such as the restoring force that arises because of density contrasts $\Delta\rho$ between the crust and the mantle (Turcotte and Schubert, 2002). Accounting for gravity would be necessary towards applying this model to a thrust or normal fault system accommodating vertical offset. Interestingly, Simpson (2015) showed that permanent deformation should be expected even in purely elastic crustal blocks around a dip-slip fault. This is because relative vertical motion, which is maximized at the fault, must vanish far away. The corresponding gradient in vertical displacement can be thought of as non-recoverable deformation that necessarily accompanies the build-up of fault offset. It is likely that
such warping of the fault blocks would increase flexural stresses to the point of plastic yielding and contribute to the accumulation of inelastic strain, on top of the deformation patterns related to the seismic cycle.

Considering buoyancy in the model would introduce a characteristic horizontal length scale to the problem, i.e., the flexural wavelength $\lambda = (4D/\Delta\rho g)^{1/4}$, where $\Delta\rho$ denotes the density contrast between the crust and the mantle. A point load situated on the fault would only lead to significant elastic deflection at distances smaller than a multiple of $\lambda$ away from the fault. With the parameters used here, and assuming $\Delta\rho = 600 \text{ kg.m}^{-3}$, $\lambda$ is about 15 to 35 km, for $H = 1$ and $H = 3$ km respectively. We would thus expect that for $L$ larger than $\lambda$, the stress ratios $R_\sigma$ and $R_{\Delta\sigma}$ controlling the behavior of the system will no longer depend on $L$ but on $\lambda$. Using the definition of the flexural parameter $\lambda$, we can eliminate $L$, and express the stress ratios in terms of the plate thickness $H$ only, leading to:

$$R_{\sigma} = \frac{\sigma_{Y}(3\Delta\rho gH)^{1/4}}{12f_{0}\sigma E'^{1/4}}, \quad R_{\Delta\sigma} = \frac{\sigma_{Y}(3\Delta\rho gH)^{1/4}}{6\Delta\sigma E'^{1/4}}. \qquad (40)$$

This suggests that the maximum yielding (minimum thickness of the plate's elastic core) would occur roughly at a distance $\lambda$ from the fault, and that as $H$ increases, the behavior of the plate gets closer to a purely elastic behavior. However, further numerical investigation accounting for buoyancy effects are required to validate this hypothesis.

# 5 Conclusion

Our elastic-plastic slider-and-springboard model provides new insights into the physics behind the long-term evolution of seismogenic plate boundaries. By coupling elasticity, dynamic friction, and plastic rheology in a simple 1-D framework, we can derive simple scaling relations that describe the partitioning of recoverable vs. non-recoverable deformation in the lithosphere adjacent to a major fault zone.

We show that the elasto-plastic rheology of the lithosphere influences the earthquake cycle and the deformation pattern around the fault zone. This influence is all the more important when the yield strength of the bulk rock is close to the frictional strength of the fault itself. In our model, permanent deflection only accumulates over successive earthquake cycles when the lithosphere is progressively damaged and weakened by cycles of loading and unloading. Here again, the rate of permanent inelastic deformation is a decreasing function of bulk yield strength, and an increasing function of the rate of damage creation. The theoretical scalings derived in this model constitute a first step towards a quantitative, mechanical interpretation of short- and long-term deformation rates measured at plate boundaries through geodesy and geomorphology.

# A Elasto-plastic bending moment under cyclic loading

Let us write $\sigma_{xx}^0$, $\omega_0$ and $M_0$ the residual stress, curvature and bending moment at the end of a loading (or unloading) phase. At later times, the stress profile within the plate obeys Equation (13). Note that the stress profile is always symmetric with respect to the neutral fiber $y = 0$. Figure 11 provides a schematic representation of the evolution of fiber stresses during loading or unloading. Again we write $h$ the thickness of the elastic core of the plate, i.e., the region where fiber stresses remain below $\sigma_p$. As illustrated in Figure 11, wherever $|y| < h/2$, the change in stress profile is purely elastic, and it is simply a rotation of the $\sigma_{xx}^0$ vs. $y$ line by an angle $\epsilon$ related to the change in curvature $\omega - \omega_0$. For $|y| > h/2$, the stress profile switches sign and goes from the yield envelope in compression to the yield envelope in tension (or vice-versa), such that the change in stress $\sigma_{xx}$ at position $y$ is twice the magnitude of $\sigma_p(y)$.

In Figure 11, distance $OD$ is $h/2$, $OA$ is $k_Y\sigma_Y$, and $OC$ is $H/2$. Along line $BD$ (at the upper limit of the elastic core), the total change in stress $\sigma_{xx} - \sigma_{xx}^0$ is given by $E'(\omega - \omega_0)h/2$. We conclude that $h$ relates to the curvature change $(\omega - \omega_0)$ as:

$$h = \frac{H}{1 + R_0 |\omega - \omega_0| / 2k_Y}, \qquad (41)$$

and that $\sigma_{xx}$ can therefore be rewritten as:

$$\sigma_{xx} = \sigma_{xx}^0 \begin{cases} E'(\omega - \omega_0)y & \text{if } |y| < h/2 \\ \operatorname{sign}(y)\operatorname{sign}(\omega - \omega_0)2\sigma_Y k_Y (1 - 2|y|/H) & \text{if } |y| > h/2 \end{cases} \qquad (42)$$

The stress profile (42) can then be integrated according to (5), and after making use of (41), we obtain:

$$M = M_0 + D(\omega - \omega_0) \frac{(1 + R_0 |\omega - \omega_0| / 4k_Y)}{(1 + R_0 |\omega - \omega_0| / 2k_Y)^2}. \qquad (43)$$
(a) unloading (coseismic)

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

(b) loading (interseismic)

<figure><img src="image_13.png" /><figcaption>**Figure 11** Schematic evolution of the stress profiles $\sigma_{xx}$ (red line) during coseismic unloading (a) and interseismic loading (b). The blue line indicates the residual stress profile $\sigma_{xx}^0$. The vertical axis is the y axis (only the upper half of the plate is shown), the horizontal axis is the stress. The oblique dotted lines indicate the yield envelope $\sigma_p$ defined in Equation (4). O, A, B, C and D are reference points used in the main text. $\epsilon$ is the angle of rotation corresponding to elastic stress changes.</figcaption></figure>

# B Numerical evaluation of the deflection history

## B.1 Numerical evaluation

Here we present the numerical scheme used to solve for the plate deflection history. As a first step, we define the non-dimensional variables $\Phi$ and $\Theta$ as:

$$ \Phi = \ln \frac{v}{v_0}, \qquad \Theta = \ln \frac{v_0 \theta}{d_c}. \qquad (44) $$

We also write $u(t) = w(L, t)$ the deflection at the right end of the plate (along the fault). Using these substitutions, we take the time derivative of Equation (2), leading to:

$$ \dot{u} = v_0(1 - e^{\Phi}). \qquad (45) $$

The state evolution Equation (25) becomes:

$$ \dot{\Theta} = \frac{v_0}{d_c} [e^{-\Theta} - e^{\Phi}]. \qquad (46) $$

We will show in the next paragraph that $\Phi$ can be seen as a function of $u$ and $\Theta$ that can not be written explicitly and requires a numerical evaluation. Therefore we write the two preceding Equations (45) and (46) as:

$$ \begin{aligned} \dot{u} &= v_0 (1 - e^{\Phi[u, \Theta]}) \\ \dot{\Theta} &= \frac{v_0}{d_c} (e^{-\Theta} - e^{\Phi[u, \Theta]}). \end{aligned} \qquad (47) $$

Here, we solve the system (47) numerically through an adaptive time step Runge-Kutta (RK) Fehlberg algorithm (Fehlberg, 1969). To do so, we need to evaluate the functional $\Phi$ at each step of the RK algorithm. The procedure is described in the next paragraph, where we also show that it requires finding the distribution of the deflection $w(x, t)$ across the plate.

Let us assume that at a particular time $t$, $u$ and $\Theta$ are known. Then the moment balance (19), could be solved for $w$ across the plate, assuming the following boundary conditions:

$$ \begin{aligned} w(0, t) &= \frac{dw}{dx}(0, t) = 0 \\ \frac{d^2 w}{dx^2}(L, t) &= 0 \\ w(L, t) &= u \end{aligned} \qquad (48) $$

Note that the first 3 conditions correspond to Equations (20) and (21). To solve these, we discretize the plate along the x axis into $N-1$ elements ($N$ being the number of nodes) of constant size $\Delta x$, and write $M_i$, $\omega_i$ and $w_i$ the bending
moment, curvature and deflection at nodes $i = 1, ..., N$. Here we use $N = 100$ nodes. The moment balance (19), the definition of curvature (1), the definition of the bending moment (5), the expression for fiber stress (13), and the boundary conditions (48) then become a set of $3N$ non-linear algebraic equations that can be solved with a standard Newton-Raphson algorithm. In doing so, we used second-order centered finite differences to approximate the first and second spatial derivatives of the variables.

Obtaining the distribution of $w$ and $M$ across the plate then allows us to evaluate the left-hand side of Equation (23), again by making use of finite difference approximations of the spatial derivatives. The right-hand side of Equation (23) then only depends on $v$ and $\theta$ (from the definition of $F$ in Equation 3 and the expression of the rate-and-state friction coefficient in Equation 24), or: on $\Phi$ and $\Theta$. Thus, if $\Theta$ is known at time $t$ in addition to $u$, this equation can be used to find $\Phi$, again using a Newton-Raphson algorithm.

It should be noted that at each time step, we need to keep track of the residual curvature $\omega_0$ and the residual stress profiles $\sigma_{xx}^0$, defined as the distributions of $w$ and $\sigma_{xx}$ when $\dot{\omega}$ changes sign, that is: when a transition between a loading and an unloading phase (or vice-versa) occurs.

## B.2 Non-dimensional equations

In this section, we derive the non-dimensional form of the equations governing the deflection of the plate, which are used in our numerical model. For this, we make the following substitutions:

$$ 
\begin{aligned}
t &\rightarrow tv_0/d_c \\
x &\rightarrow x/L \\
y &\rightarrow 2y/H \\
(\sigma_{xx}, \sigma_p) &\rightarrow (\sigma_{xx}, \sigma_p)/\sigma_Y \\
(w, u) &\rightarrow (w, u)/d_c \\
\omega &\rightarrow \omega L^2/d_c \\
M &\rightarrow 4M/(\sigma_Y H^2)
\end{aligned} 
\quad (49) $$

Equations (4) and (13) for the plate's yield envelope $\sigma_p$ and fiber stress $\sigma_{xx}$ become (in non-dimensional form):

$$ 
\begin{aligned}
\sigma_p &= \pm k_Y (1 - |y|) \\
\sigma_{xx} &= \operatorname{sign} ([\sigma_{xx}^0 + \bar{E}(\omega - \omega_0)y]) \min [|\sigma_{xx}^0| + \bar{E}(\omega - \omega_0)y|, |\sigma_p|],
\end{aligned} 
\quad (50) $$

where $k_Y$ is the ratio between the yield stress and the reference yield stress ($k_Y = 1$ in the absence of fatigue) and $\bar{E}$ is the non-dimensional modulus given by:

$$ \bar{E} = \frac{E' H d_c}{2 \sigma_Y L^2}. \quad (51) $$

$\bar{E}$ can also be seen as the non-dimensional radius of curvature $R_0$ defined in (9). The definition of the bending moment (5) becomes:

$$ M = \int_{-1}^{1} \sigma_{xx} y dy, \quad (52) $$

and the moment balance (19) writes:

$$ \frac{d^2 M}{dx^2} + \xi \frac{d^2 w}{dx^2} = 0, \quad (53) $$

where $\xi$ is a non-dimensional parameter defined as:

$$ \xi = \frac{4\sigma d_c}{\sigma_Y H}. \quad (54) $$

The boundary conditions (48) required to solve the mechanical problem write:

$$ 
\begin{aligned}
w(0, t) &= \frac{dw}{dx}(0, t) = 0 \\
\frac{d^2 w}{dx^2}(1, t) &= 0 \\
w(1, t) &= u.
\end{aligned} 
\quad (55) $$
The moment balance along the fault (23) becomes:

$$-\frac{d\mathcal{M}}{dx}(1,t) - \xi \frac{dw}{dx}(1,t) = \gamma (f_0 + a\Phi + b\Theta + \beta e^{\Phi}), \quad (56)$$

where $\Phi$ and $\Theta$ are the non-dimensional variables defined in equation (44), and $\gamma$ and $\beta$ are non-dimensional parameters given by:

$$\gamma = \frac{4\sigma L}{\sigma_Y H}, \quad \beta = \frac{\eta v_0}{\sigma}. \quad (57)$$

In writing (56) we made use of the definition of the friction coefficient (24). The system (47) then becomes:

$$\begin{cases} \dot{u} = 1 - e^{\Phi[u,\Theta]} \\ \dot{\Theta} = e^{-\Theta} - e^{\Phi[u,\Theta]}. \end{cases} \quad (58)$$

In summary, the behavior of the system is governed by a set of 7 non-dimensional parameters: $f_0$, $a$, $b$, $\bar{E}$, $\beta$, $\gamma$, and $\xi$.

The bending moment (10) corresponding to a single (interseismic) loading phase starting from zero deflection becomes:

$$\mathcal{M} = \frac{2}{3} \bar{E} \omega \frac{(1 + \bar{E} \omega / 2k_Y)}{(1 + \bar{E} \omega / k_Y)^2}. \quad (59)$$

After the first loading phase, for a cyclic loading history, the bending moment expression (14) becomes:

$$\mathcal{M} = \mathcal{M}_0 + \frac{2}{3} \bar{E}(\omega - \omega_0) \frac{(1 + \bar{E}|\omega - \omega_0|/4k_Y)}{(1 + \bar{E}|\omega - \omega_0|/2k_Y)^2}. \quad (60)$$

Note that Equation (60) is similar to (59) during the (interseismic) loading phase ($\omega > \omega_0$), provided that $\omega_0 = \mathcal{M}_0 = 0$ and $k_Y$ becomes $k_Y/2$. Thus, during the loading phase of a cyclic deflection, we can rewrite the bending moment as:

$$\mathcal{M} = \mathcal{M}_0 + \frac{2}{3} \bar{E}(\omega - \omega_0) \frac{(1 + \bar{E}(\omega - \omega_0)/4\kappa)}{(1 + \bar{E}(\omega - \omega_0)/2\kappa)^2}, \quad (61)$$

where $\kappa = k_Y/2$, $\omega_0 = \mathcal{M}_0 = 0$ for the first loading, and $\kappa = k_Y$ otherwise. In the following, we will make use of (61) to derive scaling relations that describe the plate deflection history.

Finally, from (61), the maximum bending moment the plate can sustain $\mathcal{M}_{max}$ (Equation 12) becomes in non-dimensional form:

$$\mathcal{M}_{max} = \mathcal{M}_0 + \frac{2\kappa}{3}, \quad (62)$$

which reduces to $1/3$ for a single loading phase with $\kappa = k_Y/2 = 1/2$.

The purely elastic case is obtained with $\kappa$ or $k_Y \to \infty$, that is:

$$\mathcal{M} = \mathcal{M}_0 + \frac{2}{3} \bar{E}(\omega - \omega_0). \quad (63)$$

## C Earthquake cycle scalings

In this section, we derive an approximate closed form solution for the interseismic deflection. We then derive approximate scaling relations from this solution for the mean deflection of the plate (the mean level of deflection averaged over several earthquake cycles), for the duration of the earthquake cycle, and for the rate of permanent deflection accumulation. To do so, we will assume that the deflection of the plate remains small, i.e., a fraction of the fault length, so that $dw/dx \ll 1$, and in particular $dw/dx \ll f$, $f$ being the friction coefficient, that only slightly deviates from $f_0$. In such conditions, the second term on the left-hand sides of Equations (19), (23), (53) and (56) becomes negligible. Another way of showing this property is to note that the ratio between this term and the shear force $V$ is (from Appendix B.2) approximately $\xi/\gamma \simeq f_0 d_c / L$, which is typically much smaller than unity. Integrating the moment balance (53) twice, one finally obtains (in non-dimensional form):

$$\mathcal{M} = \gamma(f + \beta e^{\Phi})(1 - x). \quad (64)$$

All the developments presented in Appendix C will be based on the non-dimensional variables defined in Appendix (B.2). The dimensional version of all the results are provided in the main text.
## C.1 Interseismic deflection profile

We first derive a closed-form approximate solution for the interseismic deflection profile. From Equation (64), the bending moment at time $t$ of the interseismic period is given by:

$$ \mathcal{M} - \mathcal{M}_0 = \frac{2}{3}\alpha(t)(1-x), \qquad (65) $$

where $\mathcal{M}_0$ is the residual bending moment at the start of an interseismic period and $\alpha(t)$ is the normalized moment of the frictional force defined as:

$$ \alpha(t) = \frac{3}{2}\gamma\Delta f(t). \qquad (66) $$

In (66), $\Delta f(t)$ is the difference between the friction coefficient at time $t$ and the friction coefficient at the start of the interseismic period. Substituting the expression of the interseismic bending moment change (61) into (65) yields a polynomial expression that can be solved to obtain the following expression for the change in curvature $\Delta\omega = \omega - \omega_0$:

$$ \frac{\bar{E}\Delta\omega}{2\kappa} = -1 + \frac{1}{\sqrt{1 - \alpha(1-x)/\kappa}}. \qquad (67) $$

Note that we have chosen the positive root for $\Delta\omega$ because the geometry of the model imposes a downward deflection. Then, the solution (67) requires that the term under the square root is always positive. This latter condition is satisfied as long as the moment of vertical shear force on the fault $2/3\alpha(1-x)$ does not exceed the maximum bending moment $\mathcal{M}_{max} = 2\kappa/3$, that is: as long as $\alpha/\kappa < 1$. Equation (67) can then be integrated twice, taking into account the boundary conditions $w(0, t) = dw/dx(0, t) = 0$. We end up with the following interseismic deflection profile $\Delta w = w - w_0$:

$$ \frac{\bar{E}\Delta w(x, \alpha)}{2\kappa} = -\frac{x^2}{2} - \frac{2\kappa x}{\alpha} \sqrt{1 - \frac{\alpha}{\kappa} - \frac{4\kappa^2}{3\alpha^2} \left[ \left(1 - \frac{\alpha}{\kappa}\right)^{3/2} - \left(1 - \frac{\alpha}{\kappa} + \frac{\alpha}{\kappa}x\right)^{3/2} \right]}. \qquad (68) $$

The variation of $\alpha(t)$ during the interseismic period can be obtained if we assume no slip on the fault. Writing $t_0$ the start time of the interseismic period, we have $\Delta w(1, \alpha) = t - t_0$, and expression (68) can be inverted to get:

$$ \alpha(t) = \frac{-8 + 4\sqrt{4 + 27(t - t_0)/t^* \left[ (t - t_0)/t^* + 1 \right]^2}}{9 \left[ (t - t_0)/t^* + 1 \right]^2}. \qquad (69) $$

where $t_0$ is the start time of the interseismic period, and $t^*$ is the (non-dimensional) characteristic time given by:

$$ t^* = \frac{1}{\bar{E}} \qquad (70) $$

Assuming $t_0 = 0$, we write Equations (68) and (69) in the more compact form:

$$ \left\{ \begin{aligned} \Delta w(x, t) &= \kappa t^* \left\{ -x^2 - \frac{4\kappa x}{\alpha(t/t^*)} \sqrt{1 - \frac{\alpha(t/t^*)}{\kappa}} - \frac{8\kappa^2}{3\alpha^2(t/t^*)} \left[ \left(1 - \frac{\alpha(t/t^*)}{\kappa}\right)^{3/2} - \left(1 - \frac{\alpha(t/t^*)}{\kappa}(1-x)\right)^{3/2} \right] \right\} \\ \alpha(s) &= \frac{-8+4\sqrt{4+27s(s+1)^2}}{9(s+1)^2}. \end{aligned} \right. \qquad (71) $$

The functional $\alpha$ is plotted against $t/t^*$ in Figure 12a, along with two different approximations obtained in the limit of small $t/t^*$. $\alpha$ can be thought of as the change of normalized bending moment during the interseismic phase. Initially it follows a linear trend, but as $t$ gets closer to $t^*$, it deviates from the linear elastic trend, and increases more slowly. $t^*$ is thus the characteristic time past which elasto-plastic effects become significant during the interseismic bending of the plate. In Figure 12b, we show the evolution of $\alpha$ for the typical duration of interseismic loading discussed in the main text, and for different values of yield stress (and thus $t^*$). The other parameters are taken from Table 1.

The same derivation can be done for a purely elastic case. Then, Equations (65) and (63) lead to:

$$ \bar{E}\Delta\omega_e = \alpha_e(t)(1-x). \qquad (72) $$

Note that we write here $\alpha$ with a $e$ subscript, since it is different from the change in moment in an elasto-plastic plate at the same time $t$. Upon integrating twice, we get the purely elastic interseismic deflection profile $\Delta w_e$ as:

$$ \bar{E}\Delta w_e(x, \alpha_e) = \alpha_e(t) \frac{x^2}{2} \left(1 - \frac{x}{3}\right). \qquad (73) $$
Again assuming that $\Delta w_e(1, t) = t - t_0$ during the interseismic stage, we can solve for $\alpha_e(t)$, and inject the result in Equation (73). We finally get:

$$ \alpha_e(t) = 3 \frac{(t - t_0)}{t^*}, \qquad (74) $$

and:

$$ \Delta w_e(x, t) = (t - t_0) \frac{x^2}{2} (3 - x). \qquad (75) $$

<figure><img src="image_14.png" /><figcaption>**Figure 12** Normalized change in bending moment during interseismic loading $\alpha$, vs normalized time $t/t^*$ (a) and vs. time (b), for different values of $\sigma_Y/\sigma^*$ (and thus $t^*$). Solid lines correspond to Equation (71), dashed and dotted lines indicate two approximations of (71), assuming $t/t^* \ll 1$. Other parameters are defined in Table 1.</figcaption></figure>
## C.2 Mean level of deflection

The closed form solutions (68) and (73) can be used to derive a scaling relation for the mean level of deflection $w_s$ around which the earthquake cycle operates (Figures 4 and 6). For that, we assume that the envelope of deflection history shown in Figure 4 results from a single elasto-plastic loading (interseismic) stage ending up in a steady equilibrium configuration with maximum deflection $w_s$. Meanwhile, the value of $\alpha$ increases from 0 to approximately $\alpha_w = 3\gamma f_0/2$.

$$w_s = \Delta w(1, \alpha_w) = \kappa t^* \left\{ -1 - \frac{4\kappa}{\alpha_w} \sqrt{1 - \frac{\alpha_w}{\kappa}} - \frac{8\kappa^2}{3\alpha_w^2} \left[ \left(1 - \frac{\alpha_w}{\kappa}\right)^{3/2} - 1 \right] \right\}, \quad (76)$$

where $\Delta w$ is given by (68). Similarly, in the purely elastic case, the mean level of deflection stabilizes at $w_e$, which from Equation (73) writes:

$$w_e = \Delta w_e(1, \alpha_w) = \frac{t^* \alpha_w}{3}. \quad (77)$$

Since we consider a single loading phase, $\kappa = k_Y/2$. We also have:

$$\frac{\alpha_w}{\kappa} = \frac{\sigma^*}{\sigma_Y k_Y}, \quad \sigma^* = 12 \frac{f_0 \sigma L}{H}. \quad (78)$$

Ultimately, we obtain:

$$\frac{w_s}{w_e} = -3 \frac{\sigma_Y k_Y}{\sigma^*} - 12 \left( \frac{\sigma_Y k_Y}{\sigma^*} \right)^2 \sqrt{1 - \frac{\sigma^*}{k_Y \sigma_Y}} - 8 \left( \frac{k_Y \sigma_Y}{\sigma^*} \right)^3 \left[ \left( 1 - \frac{\sigma^*}{k_Y \sigma_Y} \right)^{3/2} - 1 \right], \quad (79)$$

which is a decreasing function of the stress ratio $k_Y \sigma_Y / \sigma^*$. It is again noteworthy that the stick-slip regime is observed as long as the moment exerted by the frictional force on the right end of the plate does not exceed the maximum bending moment the plate can support ($\alpha_w / \kappa < 1$), which translates into:

$$k_Y \sigma_Y > \sigma^*. \quad (80)$$

## C.3 Earthquake cycle duration

Similarly, the duration of the earthquake cycle can be assessed from Equations (68) and (73). First, we assume that the interseismic phase is much longer than the co-and postseismic stages, so that the duration of the earthquake cycle $T$ (or $T_e$ in the purely elastic case) is given by the duration of the interseismic phase. Then, since no slip occurs on the fault during the interseismic loading, the (normalized) deflection accumulated interseismically is given by the (normalized) cycle duration, so that:

$$T = \Delta w(1, \alpha_T) = \kappa t^* \left\{ -1 - \frac{4\kappa}{\alpha_T} \sqrt{1 - \frac{\alpha_T}{\kappa}} - \frac{8\kappa^2}{3\alpha_T^2} \left[ \left(1 - \frac{\alpha_T}{\kappa}\right)^{3/2} - 1 \right] \right\} \\ T_e = \Delta w_e(1, \alpha_T) = \frac{t^* \alpha_T}{3}, \quad (81)$$

where $\Delta w$ and $\Delta w_e$ are given by Equations (68) and (73) respectively, and $\alpha_T = 3\gamma\Delta f/2$ is the interseismic change in the moment of the frictional force. $\Delta f$ is here the amplitude of the variation of the friction coefficient through an earthquake cycle. Since after the initial loading $\kappa = k_Y$, we can write (from the definition of $\gamma$):

$$\frac{\alpha_T}{\kappa} = \frac{\Delta \sigma^*}{k_Y \sigma_Y}, \quad \Delta \sigma^* = 6 \frac{\Delta f \sigma L}{H}. \quad (82)$$

Note that $\Delta f \sigma$ is the typical stress drop $\Delta \sigma$ associated with an earthquake on the fault. Therefore, we conclude that:

$$\frac{T}{T_e} = -3 \frac{\sigma_Y k_Y}{\Delta \sigma^*} - 12 \left( \frac{\sigma_Y k_Y}{\Delta \sigma^*} \right)^2 \sqrt{1 - \frac{\Delta \sigma^*}{k_Y \sigma_Y}} - 8 \left( \frac{k_Y \sigma_Y}{\Delta \sigma^*} \right)^3 \left[ \left( 1 - \frac{\Delta \sigma^*}{k_Y \sigma_Y} \right)^{3/2} - 1 \right]. \quad (83)$$

Note that scaling relations (79) and (83) have the same form: they invoke the same decreasing function of $k_Y \sigma_Y / \sigma^*$ for $w_s / w_e$, and of $k_Y \sigma_Y / \Delta \sigma^*$ for $T / T_e$. Note that the ratio between $\Delta \sigma^*$ and $\sigma^*$ is of the order of $\Delta f / 2f_0 \ll 1$, so that $\Delta \sigma^* \ll \sigma^*$. Thus, as $k_Y \sigma_Y / \sigma^*$ ranges from 1 to $+\infty$, $k_Y \sigma_Y / \Delta \sigma^*$ takes higher values, and is generally much
larger than 1. A first-order Taylor expansion of scaling (83) gives:

$$ \frac{T}{T_e} \approx 1 + \frac{9}{16} \frac{\Delta\sigma^*}{k_Y\sigma_Y}. \qquad (84) $$

As $\sigma_Y \to +\infty$, we recover the purely elastic solution $T = T_e$.

## C.4 Effective interseismic stiffness

In the preceding section, it has been shown that the deflection is controlled by the normalized moment of the shear force $\alpha$, and how it evolves in time (Equation (71)). As illustrated in Figure 12, the increase in interseismic deflection is associated with an increase of $\alpha$, which implies an increase in the vertical shear force acting on the right end of the plate (the fault), since the plate length is constant. We can thus define a normalized stiffness $k$ as:

$$ k = \frac{\partial\alpha}{\partial\Delta w(1, \alpha)}. \qquad (85) $$

In the purely elastic case, Equation (73) leads to:

$$ \bar{k}_e = 3\bar{E} = \frac{3}{t^*}. \qquad (86) $$

Similarly, in the elasto-plastic case, we get from Equation (71):

$$ \bar{k} = \frac{3}{t^*} \frac{(\alpha/\kappa)^3 \sqrt{1 - \alpha/\kappa}}{2 \left[ 8 - (\alpha/\kappa)^2 - 4(\alpha/\kappa) - 8\sqrt{1 - \alpha/\kappa} \right]}, \qquad (87) $$

which, from Equation (86) can be rewritten as:

$$ \bar{k} = \bar{k}_e \frac{(\alpha/\kappa)^3 \sqrt{1 - \alpha/\kappa}}{2 \left[ 8 - (\alpha/\kappa)^2 - 4(\alpha/\kappa) - 8\sqrt{1 - \alpha/\kappa} \right]}. \qquad (88) $$

The ratio $k/k_e$ is one for $\alpha = 0$ (no shear stress applied), and decreases with increasing $\alpha$. Note that we always have $\alpha < 1$ in the stick-slip regime.

## C.5 Permanent deflection with fatigue

In this section we derive a scaling relation for permanent deflection accumulation across many earthquake cycles, when the plate undergoes cyclic fatigue. Here again, we will neglect the second term describing the horizontal compressive force in the moment balance (19) or (53), for simplicity.

To derive such a relation, we build on the scaling (79). At each earthquake cycle, the yield stress $\sigma_Y$ decreases, so that, according to Equation (79), the mean level of deflection increases. Here we are interested in how much this level of deflection changes after $n$ earthquake cycles. Writing $w_s(n)$ the mean level of deflection after $n$ earthquake cycles, the permanent deflection $w_p(n)$ can thus be defined as:

$$ w_p(n) = w_s(n) - w_s(1) = w_s[\sigma_Y(n)] - w_s[\sigma_Y], \qquad (89) $$

where $\sigma_Y(n)$ is given by:

$$ \sigma_Y(n) = \sigma_Y k_Y(n) = \sigma_Y (1 - \phi \log n). \qquad (90) $$

From Equations (76), (78), (89) and (90), we get:

$$ w_p(n) = t^* \psi \left[ \frac{\sigma_Y}{\sigma^*}, \phi, n \right], \qquad (91) $$
where the functional $\psi$ is defined as:

$$ \psi[z, \phi, n] = \frac{1}{2}\phi \log n + 2z \left[ \sqrt{1 - \frac{1}{z}} - (1 - \phi \log n)^2 \sqrt{1 - \frac{1}{z(1 - \phi \log n)}} \right] + \frac{4}{3}z^2 \left\{ \left(1 - \frac{1}{z}\right)^{3/2} - 1 - (1 - \phi \log n)^3 \left[ \left(1 - \frac{1}{z(1 - \phi \log n)}\right)^{3/2} - 1 \right] \right\}. \quad (92) $$

The permanent deflection given in (91) can be expressed as a function of time, assuming that we approximately have:

$$ n \approx \frac{t - t_0}{T_e} = 3 \frac{(t - t_0)}{t^*} \frac{\sigma_Y}{\Delta\sigma^*}, \quad (93) $$

$t_0$ being the time of the first earthquake, $T_e$ the purely elastic cycle duration defined in (81), and $t^*$ the characteristic timescale defined in Equation (29). In doing so, we neglect the change in cycle duration caused by a decrease of the yield stress. This is motivated by the scaling (83) shown in Figure 6: for any value of $\sigma_Y$ larger than $\sigma^*$, the cycle duration only deviates from $T_e$ by a few percent.

The mean rate of permanent deflection between cycle $k$ and cycle $n$ is therefore approximately:

$$ R_p = \frac{w_p(n) - w_p(k)}{(n-k)T_e} = \frac{3\sigma_Y}{(n-k)\Delta\sigma^*} \left\{ \psi\left[\frac{\sigma_Y}{\sigma^*}, \phi, n\right] - \psi\left[\frac{\sigma_Y}{\sigma^*}, \phi, k\right] \right\}. \quad (94) $$

This can be recast as:

$$ R_p = \frac{3}{(n-k)} \frac{\sigma_Y}{\Delta\sigma^*} \chi \left[ \frac{\sigma_Y}{\sigma^*}, \phi, n, k \right], \quad (95) $$

where $\chi$ is given by:

$$ \chi[z, \phi, n] = \frac{1}{2}\phi \log \frac{n}{k} + 2z \left[ (1 - \phi \log k)^2 \sqrt{1 - \frac{1}{z(1 - \phi \log k)}} - (1 - \phi \log n)^2 \sqrt{1 - \frac{1}{z(1 - \phi \log n)}} \right] + \frac{4}{3}z^2 \left\{ (1 - \phi \log k)^3 \left[ \left(1 - \frac{1}{z(1 - \phi \log k)}\right)^{3/2} - 1 \right] - (1 - \phi \log n)^3 \left[ \left(1 - \frac{1}{z(1 - \phi \log n)}\right)^{3/2} - 1 \right] \right\}. \quad (96) $$

Here again, the change in cycle duration has been neglected. The rate of permanent deflection is controlled by the fatigue parameter $\phi$, and by the two non-dimensional stress ratios $\sigma_Y/\sigma^*$ and $\sigma_Y/\Delta\sigma^*$.

## Acknowledgements

The authors thank the editor Mathilde Radiguet and two anonymous reviewers for their insightful comments that improved the manuscript. J.-A.O. acknowledges funding from an Emergence(s) Ville de Paris grant, as well as ANR Tremplin grant *GeoSigMA*.

## Data and code availability

To ensure full reproducibility and ease-of-use of our framework, we provide the code used to perform the simulations at https://zenodo.org/doi/10.5281/zenodo.7858863.

## Competing interests

The authors have no competing interests.
References

Allison, K. L. and Dunham, E. M. Earthquake cycle simulations with rate-and-state friction and power-law viscoelasticity. *Tectonophysics*, 733:232-256, May 2018. doi: 10.1016/j.tecto.2017.10.021.

Allison, K. L. and Dunham, E. M. Influence of Shear Heating and Thermomechanical Coupling on Earthquake Sequences and the Brittle-Ductile Transition. *Journal of Geophysical Research: Solid Earth*, 126(6), June 2021. doi: 10.1029/2020jb021394.

Atkinson, B. K. Subcritical crack growth in geological materials. *Journal of Geophysical Research: Solid Earth*, 89(B6):4077-4114, June 1984. doi: 10.1029/jb089ib06p04077.

Avouac, J.-P. *MOUNTAIN BUILDING, EROSION, AND THE SEISMIC CYCLE IN THE NEPAL HIMALAYA*, page 1-80. Elsevier, 2003. doi: 10.1016/s0065-2687(03)46001-9.

Avouac, J.-P. From Geodetic Imaging of Seismic and Aseismic Fault Slip to Dynamic Modeling of the Seismic Cycle. *Annual Review of Earth and Planetary Sciences*, 43(1):233-271, May 2015. doi: 10.1146/annurev-earth-060614-105302.

Baden, C. W., Shuster, D. L., Aron, F., Fosdick, J. C., Bürgmann, R., and Hilley, G. E. Bridging earthquakes and mountain building in the Santa Cruz Mountains, CA. *Science Advances*, 8(8), Feb. 2022. doi: 10.1126/sciadv.abi6031.

Baker, A., Allmendinger, R. W., Owen, L. A., and Rech, J. A. Permanent deformation caused by subduction earthquakes in northern Chile. *Nature Geoscience*, 6(6):492-496, Apr. 2013. doi: 10.1038/ngeo1789.

Barbot, S. Asthenosphere flow modulated by megathrust earthquake cycles. Apr. 2018. doi: 10.31223/osf.io/g8ahe.

Bhat, H. S., Rosakis, A. J., and Sammis, C. G. A Micromechanics Based Constitutive Model for Brittle Failure at High Strain Rates. *Journal of Applied Mechanics*, 79(3), Apr. 2012. doi: 10.1115/1.4005897.

Brantut, N., Heap, M., Meredith, P., and Baud, P. Time-dependent cracking and brittle creep in crustal rocks: A review. *Journal of Structural Geology*, 52:17-43, July 2013. doi: 10.1016/j.jsg.2013.03.007.

Buck, W. R. flexural rotation of normal faults. *Tectonics*, 7(5):959-973, Oct. 1988. doi: 10.1029/tc007i005p00959.

Burgette, R. J., Weldon, R. J., and Schmidt, D. A. Interseismic uplift rates for western Oregon and along-strike variation in locking on the Cascadia subduction zone. *Journal of Geophysical Research: Solid Earth*, 114(B1), Jan. 2009. doi: 10.1029/2008jb005679.

Burridge, R. and Knopoff, L. Model and theoretical seismicity. *Bulletin of the Seismological Society of America*, 57(3):341-371, June 1967. doi: 10.1785/bssa0570030341.

Bürgmann, R., Kogan, M. G., Steblov, G. M., Hilley, G., Levin, V. E., and Apel, E. Interseismic coupling and asperity distribution along the Kamchatka subduction zone. *Journal of Geophysical Research: Solid Earth*, 110(B7), July 2005. doi: 10.1029/2005jb003648.

Carlson, J. M., Langer, J. S., and Shaw, B. E. Dynamics of earthquake faults. *Reviews of Modern Physics*, 66(2):657-670, Apr. 1994. doi: 10.1103/revmodphys.66.657.

Cattin, R. and Avouac, J. P. Modeling mountain building and the seismic cycle in the Himalaya of Nepal. *Journal of Geophysical Research: Solid Earth*, 105(B6):13389-13407, June 2000. doi: 10.1029/2000jb900032.

Cerfontaine, B. and Collin, F. Cyclic and Fatigue Behaviour of Rock Materials: Review, Interpretation and Research Perspectives. *Rock Mechanics and Rock Engineering*, 51(2):391-414, Oct. 2017. doi: 10.1007/s00603-017-1337-5.

Dal Zilio, L., Hetényi, G., Hubbard, J., and Bollinger, L. Building the Himalaya from tectonic to earthquake scales. *Nature Reviews Earth & Environment*, 2(4):251-268, Mar. 2021a. doi: 10.1038/s43017-021-00143-1.

Dal Zilio, L., Lapusta, N., Avouac, J.-P., and Gerya, T. Subduction earthquake sequences in a non-linear visco-elasto-plastic megathrust. *Geophysical Journal International*, 229(2):1098-1121, Dec. 2021b. doi: 10.1093/gji/ggab521.

Davis, D., Suppe, J., and Dahlen, F. A. Mechanics of fold-and-thrust belts and accretionary wedges. *Journal of Geophysical Research: Solid Earth*, 88(B2):1153-1172, Feb. 1983. doi: 10.1029/jb088ib02p01153.

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

Erickson, B. A., Dunham, E. M., and Khosravifar, A. A finite difference method for off-fault plasticity throughout the earthquake cycle. *Journal of the Mechanics and Physics of Solids*, 109:50-77, Dec. 2017. doi: 10.1016/j.jmps.2017.08.002.

Fehlberg, E. Low-order classical Runge-Kutta formulas with stepsize control and their application to some heat transfer problems. 1969.

Fletcher, H. J., Beavan, J., Freymueller, J., and Gilbert, L. High interseismic coupling of the Alaska Subduction Zone SW of Kodiak Island inferred from GPS data. *Geophysical Research Letters*, 28(3):443-446, Feb. 2001. doi: 10.1029/2000gl012258.

Helmstetter, A. and Shaw, B. E. Afterslip and aftershocks in the rate-and-state friction law. *Journal of Geophysical Research: Solid Earth*, 114 (B1), Jan. 2009. doi: 10.1029/2007jb005077.

Jolivet, R., Simons, M., Duputel, Z., Olive, J., Bhat, H. S., and Bletery, Q. Interseismic Loading of Subduction Megathrust Drives Long-Term Uplift in Northern Chile. *Geophysical Research Letters*, 47(8), Apr. 2020. doi: 10.1029/2019gl085377.

Kanagawa, K., Cox, S. F., and Zhang, S. Effects of dissolution-precipitation processes on the strength and mechanical behavior of quartz gouge at high-temperature hydrothermal conditions. *Journal of Geophysical Research: Solid Earth*, 105(B5):11115-11126, May 2000. doi: 10.1029/2000jb900038.

King, G. C. P., Stein, R. S., and Rundle, J. B. The Growth of Geological Structures by Repeated Earthquakes. 1. Conceptual Framework. Journal
of Geophysical Research: Solid Earth, 93(B11):13307-13318, Nov. 1988. doi: 10.1029/jb093ib11p13307.

Kirby, S. H. Introduction and digest to the Special Issue on Chemical Effects of Water on the Deformation and Strengths of Rocks. *Journal of Geophysical Research: Solid Earth*, 89(B6):3991-3995, June 1984. doi: 10.1029/jb089ib06p03991.

Klein, E., Fleitout, L., Vigny, C., and Garaud, J. Afterslip and viscoelastic relaxation model inferred from the large-scale post-seismic deformation following the 2010Mw8.8 Maule earthquake (Chile). *Geophysical Journal International*, 205(3):1455-1472, Mar. 2016. doi: 10.1093/gji/ggw086.

Lambert, V. and Barbot, S. Contribution of viscoelastic flow in earthquake cycles within the lithosphere-asthenosphere system. *Geophysical Research Letters*, 43(19), Oct. 2016. doi: 10.1002/2016gl070345.

Lapusta, N., Rice, J. R., Ben-Zion, Y., and Zheng, G. Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate- and state-dependent friction. *Journal of Geophysical Research: Solid Earth*, 105(B10):23765-23789, Oct. 2000. doi: 10.1029/2000jb900250.

Lavé, J. and Avouac, J. P. Fluvial incision and tectonic uplift across the Himalayas of central Nepal. *Journal of Geophysical Research: Solid Earth*, 106(B11):26561-26591, Nov. 2001. doi: 10.1029/2001jb000359.

Lay, T. The surge of great earthquakes from 2004 to 2014. *Earth and Planetary Science Letters*, 409:133-146, Jan. 2015. doi: 10.1016/j.epsl.2014.10.047.

Loveless, J. P. and Meade, B. J. Two decades of spatiotemporal variations in subduction zone coupling offshore Japan. *Earth and Planetary Science Letters*, 436:19-30, Feb. 2016. doi: 10.1016/j.epsl.2015.12.033.

Madella, A. and Ehlers, T. A. Contribution of background seismicity to forearc uplift. *Nature Geoscience*, 14(8):620-625, June 2021. doi: 10.1038/s41561-021-00779-0.

Mallick, R., Burgmann, R., Johnson, K. M., and Hubbard, J. A unified framework for earthquake sequences and the growth of geological structure in fold-thrust belts. Mar. 2021. doi: 10.1002/essoar.10506463.1.

Marone, C. LABORATORY-DERIVED FRICTION LAWS AND THEIR APPLICATION TO SEISMIC FAULTING. *Annual Review of Earth and Planetary Sciences*, 26(1):643-696, May 1998a. doi: 10.1146/annurev.earth.26.1.643.

Marone, C. The effect of loading rate on static friction and the rate of fault healing during the earthquake cycle. *Nature*, 391(6662):69-72, Jan. 1998b. doi: 10.1038/34157.

Meade, B. J. The signature of an unbalanced earthquake cycle in Himalayan topography? *Geology*, 38(11):987-990, Oct. 2010. doi: 10.1130/g31439.1.

Melnick, D. Rise of the central Andean coast by earthquakes straddling the Moho. *Nature Geoscience*, 9(5):401-407, Mar. 2016. doi: 10.1038/ngeo2683.

Menant, A., Angiboust, S., Gerya, T., Lacassin, R., Simoes, M., and Grandin, R. Transient stripping of subducting slabs controls periodic forearc uplift. *Nature Communications*, 11(1), Apr. 2020. doi: 10.1038/s41467-020-15580-7.

Mia, M. S., Abdelmeguid, M., and Elbanna, A. Spatio-temporal clustering of seismicity enabled by off-fault plasticity. Apr. 2022. doi: 10.31223/x50p8b.

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

Mouslopoulou, V., Oncken, O., Hainzl, S., and Nicol, A. Uplift rate transients at subduction margins due to earthquake clustering. *Tectonics*, 35(10):2370-2384, Oct. 2016. doi: 10.1002/2016tc004248.

Métois, M., Socquet, A., and Vigny, C. Interseismic coupling, segmentation and mechanical behavior of the central Chile subduction zone. *Journal of Geophysical Research: Solid Earth*, 117(B3), Mar. 2012. doi: 10.1029/2011jb008736.

Okubo, K., Bhat, H. S., Rougier, E., Marty, S., Schubnel, A., Lei, Z., Knight, E. E., and Klinger, Y. Dynamics, Radiation, and Overall Energy Budget of Earthquake Rupture With Coseismic Off-Fault Damage. *Journal of Geophysical Research: Solid Earth*, 124(11):11771-11801, Nov. 2019. doi: 10.1029/2019jb017304.

Oryan, B., Olive, J.-A., Jolivet, R., Malatesta, L. C., Gailleton, B., and Bruhat, L. Megathrust locking encoded in subduction landscapes. *Science Advances*, 10(17), Apr. 2024. doi: 10.1126/sciadv.adl4286.

Ozawa, T., Tabei, T., and Miyazaki, S. Interplate coupling along the Nankai Trough off southwest Japan derived from GPS measurements. *Geophysical Research Letters*, 26(7):927-930, Apr. 1999. doi: 10.1029/1999gl900145.

Periollat, A., Radiguet, M., Weiss, J., Twardzik, C., Amitrano, D., Cotte, N., Marill, L., and Socquet, A. Transient Brittle Creep Mechanism Explains Early Postseismic Phase of the 2011 Tohoku-Oki Megathrust Earthquake: Observations by High-Rate GPS Solutions. *Journal of Geophysical Research: Solid Earth*, 127(8), Aug. 2022. doi: 10.1029/2022jb024005.

Prawirodirdjo, L., McCaffrey, R., Chadwell, C. D., Bock, Y., and Subarya, C. Geodetic observations of an earthquake cycle at the Sumatra subduction zone: Role of interseismic strain segmentation. *Journal of Geophysical Research: Solid Earth*, 115(B3), Mar. 2010. doi: 10.1029/2008jb006139.

Reid, H. Elastic rebound theory. *Univ Calif Publ. Bull Dept Geol Sci*, 6:413-433, 1910.

Rice, J. R. Spatio-temporal complexity of slip on a fault. *Journal of Geophysical Research: Solid Earth*, 98(B6):9885-9907, June 1993. doi: 10.1029/93jb00191.
Rice, J. R. and Tse, S. T. Dynamic motion of a single degree of freedom system following a rate and state dependent friction law. *Journal of Geophysical Research: Solid Earth*, 91(B1):521-530, Jan. 1986. doi: 10.1029/jb091ib01p00521.

Rodriguez Padilla, A. M., Oskin, M. E., Milliner, C. W. D., and Plesch, A. Accrual of widespread rock damage from the 2019 Ridgecrest earthquakes. *Nature Geoscience*, 15(3):222-226, Feb. 2022. doi: 10.1038/s41561-021-00888-w.

Rolandone, F. *The Seismic Cycle: From Observation to Modeling*. John Wiley & Sons, 2022.

Ruh, J. B. and Vergés, J. Effects of reactivated extensional basement faults on structural evolution of fold-and-thrust belts: Insights from numerical modelling applied to the Kopet Dagh Mountains. *Tectonophysics*, 746:493-511, Oct. 2018. doi: 10.1016/j.tecto.2017.05.020.

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

Sagiya, T. and Meneses-Gutierrez, A. Geodetic and Geological Deformation of the Island Arc in Northeast Japan Revealed by the 2011 Tohoku Earthquake. *Annual Review of Earth and Planetary Sciences*, 50(1):345-368, May 2022. doi: 10.1146/annurev-earth-032320-074429.

Saillard, M., Audin, L., Rousset, B., Avouac, J.-P., Chlieh, M., Hall, S. R., Husson, L., and Farber, D. L. From the seismic cycle to long-term deformation: linking seismic coupling and Quaternary coastal geomorphology along the Andean megathrust: Interseismic Coupling/Coastal Morphology. *Tectonics*, 36(2):241-256, Feb. 2017. doi: 10.1002/2016tc004156.

Savage, J. C. A dislocation model of strain accumulation and release at a subduction zone. *Journal of Geophysical Research: Solid Earth*, 88 (B6):4984-4996, June 1983. doi: 10.1029/jb088ib06p04984.

Scholz, C. H. Static fatigue of quartz. *Journal of Geophysical Research*, 77(11):2104-2114, Apr. 1972. doi: 10.1029/jb077i011p02104.

Scholz, C. H. *The Mechanics of Earthquakes and Faulting*. Cambridge University Press, May 2002. doi: 10.1017/cbo9780511818516.

Simpson, G. Accumulation of permanent deformation during earthquake cycles on reverse faults. *Journal of Geophysical Research: Solid Earth*, 120(3):1958-1974, Mar. 2015. doi: 10.1002/2014jb011442.

Simpson, G. Emergence and growth of faults during earthquakes: Insights from a dynamic elasto-plastic continuum model. *Tectonophysics*, 868:230089, Dec. 2023. doi: 10.1016/j.tecto.2023.230089.

Stevens, V. L. and Avouac, J. P. Interseismic coupling on the main Himalayan thrust. *Geophysical Research Letters*, 42(14):5828-5837, July 2015. doi: 10.1002/2015gl064845.

Tenthorey, E., Cox, S. F., and Todd, H. F. Evolution of strength recovery and permeability during fluid-rock reaction in experimental fault zones. *Earth and Planetary Science Letters*, 206(1-2):161-172, Jan. 2003. doi: 10.1016/s0012-821x(02)01082-8.

Thomas, M. Y. and Bhat, H. S. Dynamic evolution of off-fault medium during an earthquake: a micromechanics based model. *Geophysical Journal International*, 214(2):1267-1280, May 2018. doi: 10.1093/gji/ggy129.

Turcotte, D. L. and Schubert, G. *Geodynamics*. Cambridge university press, 2002.

van Dinther, Y., Gerya, T. V., Dalguer, L. A., Mai, P. M., Morra, G., and Giardini, D. The seismic cycle at subduction thrusts: Insights from seismo-thermo-mechanical models. *Journal of Geophysical Research: Solid Earth*, 118(12):6183-6202, Dec. 2013. doi: 10.1002/2013jb010380.

Wang, K. 17. *Elastic and Viscoelastic Models of Crustal Deformation in Subduction Earthquake Cycles*, page 540-575. Columbia University Press, Dec. 2007. doi: 10.7312/dixo13866-017.

Wang, K., Wells, R., Mazzotti, S., Hyndman, R. D., and Sagiya, T. A revised dislocation model of interseismic deformation of the Cascadia subduction zone. *Journal of Geophysical Research: Solid Earth*, 108(B1), Jan. 2003. doi: 10.1029/2001jb001227.

Wang, K., Hu, Y., and He, J. Deformation cycles of subduction earthquakes in a viscoelastic Earth. *Nature*, 484(7394):327-332, Apr. 2012. doi: 10.1038/nature11032.

Watanabe, S.-i., Ishikawa, T., Nakamura, Y., and Yokota, Y. Co- and postseismic slip behaviors extracted from decadal seafloor geodesy after the 2011 Tohoku-oki earthquake. *Earth, Planets and Space*, 73(1), Aug. 2021. doi: 10.1186/s40623-021-01487-0.

Zhao, B., Bürgmann, R., Wang, D., Zhang, J., Yu, J., and Li, Q. Aseismic slip and recent ruptures of persistent asperities along the Alaska-Aleutian subduction zone. *Nature Communications*, 13(1), June 2022. doi: 10.1038/s41467-022-30883-7.

The article *Inelastic deformation accrued over multiple seismic cycles: Insights from an elastic-plastic slider-and-springboard model* © 2024 by Pierre Dublanchet is licensed under *CC BY 4.0*.