# Validation of a Multilayered Analog Model Integrating Crust-Mantle Visco-Elastic Coupling to Investigate Subduction Megathrust Earthquake Cycle

Yannick Caniven, Stéphane Domínguez

► To cite this version:

Yannick Caniven, Stéphane Domínguez. Validation of a Multilayered Analog Model Integrating Crust-Mantle Visco-Elastic Coupling to Investigate Subduction Megathrust Earthquake Cycle. Journal of Geophysical Research: Solid Earth, 2021, 126 (2), pp.e2020JB020342. 10.1029/2020JB020342. hal-03185699

HAL Id: hal-03185699

https://hal.science/hal-03185699

Submitted on 30 Mar 2021

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.
# Validation of a Multilayered Analog Model Integrating Crust-Mantle Visco-Elastic Coupling to Investigate Subduction Megathrust Earthquake Cycle

Yannick Caniven<fnref n="1" /> and Stéphane Domínguez<fnref n="2" />

## Key Points:

* Analog model rheology allows for reproducing pre-, inter- and coseismic phases, followed by after-slip and postseismic relaxation

* A broad variability of slip events, from aseismic slow slips to dynamic lab quakes, is observed depending on the boundary conditions

* The visco-elastic coupling between the crust and the mantle wedge partly controls the frictional behavior and the earthquake size

## Abstract

We have developed a scaled analog model of a subduction zone simulating seismic cycle deformation phases. Its rheology is based on multilayered visco-elasto-plastic materials to account for the mechanical behavior of a continental lithospheric plate overriding a subducting oceanic plate. The seismogenic zone displays unstable slip behavior, extending at depth into a weak interface with stable slip properties. The model succeeds in reproducing interseismic phases interrupted by coseismic ruptures and followed by after-slip. The experimental data catalog shows a broad variability of slip events from aseismic slow slips to fast dynamic lab quakes. Results also show the occurrence of both isolated and precursory slow-slip events arising before the mainshocks. Given the absence of fluids in the model, the broad variability in slip event velocity can be attributed to fault roughness complexity. The model rheology induces also a key visco-elastic coupling between the elastic overriding plate and the mantle wedge allowing, for the first time, to reproduce experimentally a realistic postseismic visco-elastic relaxation phase. Preliminary results reveal that the tectonic loading rate modulates this visco-elastic coupling. A low loading rate weakens it, which increase the amount of storable interseismic elastic deformation, and favors the occurrence of large megathrust events. A high loading rate strengthens it, which minimize the accumulation of interseismic elastic deformation, the slip-event sizes, and promote aseismic creep. This new scaled-analog subduction model is a complementary tool to investigate earthquake mechanics and improve the interpretation of geodetic and seismological records.

<footnote n="1">¹Department of Earth, Environmental, and Planetary Sciences, Rice University, Houston, TX, USA</footnote>

<footnote n="2">²Géosciences Montpellier, Université de Montpellier, CNRS, Montpellier, France</footnote>
# 1. Introduction

The recent Japanese and Sumatra major subduction earthquakes have caused dramatic economic and human losses (Information [NCEI], n.d.). It is, then, of high scientific and societal interest to improve our understanding of how these catastrophic natural events nucleate and generate crustal deformations. To achieve this challenge, the difficulties to investigate earthquake deformations undersea and to assess the specific temporal scales of seismic cycle processes must be overcome (Subarya et al., 2006).

Scientific investigations conducted on recent megathrust earthquakes illustrate the impact of these limitations (e.g., Goldfinger et al., 2013; Ozawa et al., 2011). During the Tohoku-Oki earthquake (2011, $M_w$ = 9.1), experimental offshore geodetic measurements, recorded for the first time motion near the tip of the overriding plate, and revealed that coseismic slip amplified seaward to reach a maximum (>50 m) close to the trench, as confirmed by tsunami data inversion (e.g., linuma et al., 2012; Koketsu et al., 2011). Until then, inversions of inland geophysical data predicted that coseismic slip distributions steeply tapered upward to the trench (e.g., Ito et al., 2011). This discovery questions our degree of knowledge and the robustness of the mechanical models used to constrain the main properties of the seismogenic zone (SZ) and the tsunamogenic potential of subduction faults.

Open scientific questions also remain regarding deep deformation processes. Among those, how friction heterogeneity controls the interplate mechanical coupling appears as a crucial issue. In many subduction zones, like Cascadia, geodetic data reveal complex spatio-temporal patterns with slow sleep events and quiet seismic activity phases. These transient events are likely related to variations of deep coupling controlled by upper plate (UP) stress state, interplate friction properties (e.g., Pacheco et al., 1993), and pore fluid pressurization processes (Bizzarri, 2009; Faulkner & Rutter, 2001; Warren-Smith et al., 2019). The brittle/ductile transition, which marks the downdip limit of the SZ, is also a key feature notably controlling the seismicity and surface deformation (e.g., Hyndman et al., 1997). Several studies especially show that transient behavior may also be related to the interaction of unstable brittle and stable ductile behavior in mélanges in the transitions zones of the megathrust fault zone (Fagereng & Sibson, 2010; Hayman & Lavier, 2014; Saffer & Wallace, 2015; Tong & Lavier, 2018).

Despite these useful studies and promising investigations of exhumed subduction zones (e.g., Ujiie et al., 2018), the mechanical transition zone remains difficult to constrain because of its deep location (>30 km), which prevents collecting relevant information (Figure 1). Similarly, geophysical data partly cover the interseismic phase duration but are frequently extrapolated over longer periods assuming tectonic steady-rate conditions. These observations highlight the time-scale extent issue in the earthquake cycle, which is far longer than the span of modern and even historical scientific observations.
<figure><img src="image_1.png" /><figcaption>**Figure 1.** Schema of an idealized Japan type subduction zone configuration.</figcaption></figure>

To overcome part of these time and spatial scale constraints, new experimental approaches emerged to investigate the physics of earthquakes (Rosenau et al., 2017, and references therein). Among those, scaled analog models dedicated to the simulation of subduction zone seismic cycle and associated deformation phases (e.g., Corbi et al., 2013, 2017, 2019, 2020, 2011; Rosenau et al., 2019, 2009; van Dinther et al., 2013). These models use elasto-plastic analog material wedges (gelatin in the case of Corbi's experiments, described by Di Giuseppe et al., 2009 and rubber-sugar mixture for Rosenau and Van-Dinther's experiments, described by Rosenau et al., 2009) slipping on a low angle tilted plane. The SZ consists of a frictional interface limited toward the surface and at depth by two low friction zones (Corbi et al., 2013; Rosenau et al., 2009). These models succeed in reproducing stick-slip frictional behavior and seismic cycle phases.

Following these previous works, and inspired by pioneer studies (e.g., Anooshehpoor & Brune, 2003; Brune, 1973; Brune & Anooshehpoor, 1998; Day et al., 2008; McBean et al., 2015), we developed a new type of scaled analog models designed to simulate strike-slip fault earthquakes, and seismic cycle phases, including for the first time realistic postseismic deformations (Caniven et al., 2015). To achieve this goal, we used a multilayered visco-elasto-plastic model reproducing the rheology and associated mechanical behavior of the upper and lower crust. This choice of rheology appeared as a critical building-block to simulate visco-elastic coupling between the layers, postseismic deformation, and far-field stress transfers. Numerous experiments revealed interesting results, highlighting the role of key parameters on the fault slip patterns and seismic cycle dynamics such as the tectonic loading rate (Caniven, 2014) or the along-fault normal stress distribution (Caniven et al., 2017).

Following the same experimental methodology and multilayered rheology, we developed a scaled analog model of a subduction zone to simulate megathrust events and realistic earthquake cycle phases. In such a geodynamic context, the viscoelastic coupling between the crust and the mantle wedge is a crucial parameter. Without it, the postseismic
viscoelastic relaxation phase won't be observable as well as other associated deformation processes.

The main objectives of the present study are to report and validate a new analog modeling method by comparing experimental results to natural cases, discussing model capabilities, and limitations. With these aims, we present also preliminary results concerning the impact of loading rate on seismic cycle properties because the deformation processes and couplings that emerge here are directly related to the multilayered rheology of the analog model.

By studying model kinematics and mechanical behavior, we expect to better constrain the effects of main physical couplings and boundary conditions on the seismic cycle and earthquake dynamics, leading to a better interpretation of geophysical data including satellite measurements and seismology.

## 2. Experimental Model

### 2.1. Experimental Set-Up, Boundary Conditions, and Model Rheology

The experimental set-up is constituted of a 2 m long rigid aluminum frame supporting a 2 cm thick and 30 cm wide constrained and bent PVC plate simulating the bending of the subducting oceanic lithosphere (Figure 2). The modeled subducting slab is pasted on a 1 mm thick and 6 cm wide aluminum plate slipping on the PVC plate, via a viscous coupling to ensure a constant and stable velocity. The subducting plate is pulled beneath a rigid, 3 cm thick and 30 cm wide, vertical aluminum backstop using a computerized steeping-motor and mechanical reductor. Using this device, the subduction velocity can be set from 0.2 to 10 µm/s (~0.1–3.6 cm/h). The rear part of the 6 cm wide modeled UP lithosphere is pasted on a 5 mm thick vertical aluminum plate, equipped with a strain gauge fixed to the rigid backstop. The UP has a wedge shape and rests freely on the subducting slab trenchward and on a ductile material reservoir landward (Figure 2). Note that the device is not equipped with lateral glass walls, and consequently, it is not affected by lateral boundary bias.
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Experimental setup and monitoring system.</figcaption></figure>

Based on feedback from the strike-slip fault model (Caniven et al., 2015, 2017), we use a similar multilayered visco-elasto-plastic rheology to reproduce the mechanical behavior of a subduction zone at the lithospheric scale (Figure 2). The layer simulating the ductile mechanical behavior of the UP lithospheric mantle is constituted by a silicone compound (PDMS-SGM 36, Dow Corning) characterized by a viscosity of 40 kPa.s at room temperature (Rudolf et al., 2016). The layers simulating the elasto-plastic seismogenic UP crust and subducting oceanic crust are composed of high resilience polyurethane (PU) foams with a Young modulus, a shear modulus, and a Poisson ratio of 95 ± 10 kPa, 50 ± 5 kPa, and 0.15 ± 0.05, respectively. Note that the Poisson ratio is calculated in the range of PU foam deformation undergone during the experiment (<5%). The base of the elastic UP soaks into the first millimeters of the silicone wedge, ensuring a tectonic loading rate-dependent mechanical coupling linking the viscous mantle wedge and the elastic crust (see Figure S1 for a detailed view of the contact foam-silicone). This is referred to here as the visco-elastic coupling. The rheology of this part provides a transitional behavior with a shear modulus $G_{Foam-Si} \ge G_{foam}$ (50 kPa) and a dynamic viscosity $\eta_{Foam-Si} >> \eta_{Si}$ ($10^4$ Pa.s$^{-1}$). A thin upper layer (<3–4 mm) of granular material (glass or lead microbeads) is deposited on the top of the UP to adjust the imposed vertical load (60–400 Pa).
The subduction interface is made up of two parts. The upper part simulates the SZ along the contact between the upper and subducting foam plates. It dips 10°–11° and is 0.5 m long (± 150 km). The two fault surfaces (foam-foam) are covered by a specific epoxy resin coating (Caniven et al., 2015) to adjust the static and dynamic friction coefficients to realistic values (about 0.6 ± 0.1 and 0.3 ± 0.1, respectively). Low or high friction asperities can be simulated by sprinkling silica microspheres or by pasting angular quartz particles. The deeper part of the subduction interface extends 20 cm at depth (about 60–70 km) under the ductile silicone wedge layer that is separated from the subduction interface by a very thin and smooth Teflon sheet. This both ensures the silicone reservoir tightness and imposes a stable slip behavior, as commonly observed at such depth (>30 km). The transition between these two distinct parts of the interface is referred to as the Stable-Unstable Transition Zone (SUTZ).

## 2.2. Model Scaling

Analog model Scaling is achieved using well-known theoretical concepts developed by Hubbert (1937) and Cobbold and Jackson (1992). To scale the subduction model, we used the same key physical and mechanical parameters as for our strike-slip fault model (Caniven et al., 2015, 2017). At the seismic time scale, the crust elasticity is one of the main mechanical parameters controlling its deformation. In numerical/analytical modeling, the whole crust is considered as an infinite elastic half-space (e.g., Okada, 1985). The viscosity of the lower crust and the lithospheric mantle, as well as the frictional properties along the seismogenic plate interface, are also important parameters strongly affecting model dynamic behavior.

Spatial scaling was determined considering the maximum available model set-up size (~1.4 m). We impose, then, a geometric factor ($L^*$) of $3 \times 10^{-6}$ to model a ~400 km long subduction zone (Table 1). To scale model elasticity, the polyurethane foam has been selected with a shear modulus ($G$ ~ 50 kPa) to scale order value for nature ($G$~10–15 GPa) using $L^*$ too. For scaling viscosity, we consider a dynamic viscosity ($\eta$) for mantle wedge in the order of $10^{19}$ Pa s (Manea & Gurnis, 2007, and references therein). With $\eta = 4 \times 10^4$ Pa.s for the silicone (Rudolf et al., 2016), it gives a viscosity scaling factor $\eta^* = 4 \times 10^{-15}$. Viscosity scaling is then considered valid if $\eta^*$ matches with the ratio of stress factor ($\sigma^*$) to strain rate factor ($\varepsilon^*$), with $\sigma^*$ the ratio between $L^*$ and the density factor $\rho^*$ in a natural gravity field. Considering a mantle density of 3.3 (Duesterhoeft et al., 2014), an $\varepsilon^*$ of $10^{10}$ (Table 1) in our experiments, and a silicone density of 0.97, we find $\eta^* = 10^{-15}$, hence consistently in the order of the expected value. To scale model frictional properties, we considered the main rate-and-state constitutive parameters (Scholz, 1998). We examined the critical slip distance $D_c$, the critical nucleation size $L_c$, and the friction rate parameter (a–b). The related theory, definition, and determination of these parameters are detailed in Caniven et al., (2015) and Caniven et al., (2017). The model scaling is correct if $D_c$ and $L_c$ scale to nature using $L^*$ and if (a–b) remains identical between model and
<table><caption>**Table 1**
Model Scaling Parameters</caption><thead><tr><th>Parameters</th><th>Symbol</th><th>Dim. (MLT)</th><th>Unit</th><th>Model (M)</th><th>Nature (N)</th><th>Scaling factor (M/N)</th></tr></thead><tbody><tr><td><strong>Physic</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Crust density (mean)</td><td>ρ<sub>c</sub></td><td>M/L<sup>3</sup></td><td>Kg.m<sup>-3</sup></td><td>39</td><td>3,000</td><td>1.3 × 10<sup>-2</sup></td></tr><tr><td><strong>Geometric</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Length (Fault length)</td><td>L<sub>f</sub></td><td>L</td><td>m</td><td>1</td><td>3 × 10<sup>5</sup></td><td>3.3 × 10<sup>-6</sup></td></tr><tr><td>Coseismic fault slip (mean)</td><td>D<sub>mean</sub></td><td>L</td><td>m</td><td>6 × 10<sup>-5</sup></td><td>~5</td><td>1.2 × 10<sup>-5</sup></td></tr><tr><td><strong>Mechanic</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Crust Shear modulus</td><td>G<sub>c</sub></td><td>M/LT<sup>2</sup></td><td>Pa</td><td>5 × 10<sup>4</sup></td><td>~10<sup>10</sup></td><td>5 × 10<sup>-6</sup></td></tr><tr><td>Crust Poisson ratio</td><td>ν</td><td>—</td><td>—</td><td>~0.15</td><td>0.25</td><td>0.6</td></tr><tr><td>Crust Static friction coefficient</td><td>μ</td><td>—</td><td>—</td><td>0.6</td><td>0.5</td><td>1.2</td></tr><tr><td>Friction Rate</td><td>a-b</td><td>—</td><td>—</td><td>~ -0.02</td><td>-0.015</td><td>~1</td></tr><tr><td>Crust stress</td><td>σ</td><td>M/LT<sup>2</sup></td><td>Pa</td><td>500</td><td>1.5 × 10<sup>8</sup></td><td>3 × 10<sup>-6</sup></td></tr><tr><td>Coseismic slip gradient (mean)</td><td>γ<sub>co</sub></td><td>—</td><td>—</td><td>5 × 10<sup>-4</sup></td><td>5 × 10<sup>-5</sup></td><td>10</td></tr><tr><td>Lithospheric mantle viscosity</td><td>η<sub>lc</sub></td><td>M/LT</td><td>Pa.s</td><td>4 × 10<sup>4</sup></td><td>10<sup>19</sup></td><td>4 × 10<sup>-15</sup></td></tr><tr><td><strong>Kinematic</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Gravitational acceleration</td><td>g</td><td>L/T<sup>2</sup></td><td>m.s<sup>-2</sup></td><td>9.8</td><td>9.8</td><td>1</td></tr><tr><td>Strain rate (mean)</td><td>ε<sup>o</sup></td><td>T<sup>-1</sup></td><td>s<sup>-1</sup></td><td>10<sup>-4</sup></td><td>10<sup>-14</sup></td><td>1 × 10<sup>10</sup></td></tr><tr><td>Mean velocity (interseismic)</td><td>V<sub>i</sub></td><td>L/T</td><td>m.s<sup>-1</sup></td><td>3 × 10<sup>-6</sup></td><td>2 × 10<sup>-9</sup></td><td>1.5 × 10<sup>3</sup></td></tr><tr><td>Mean velocity (coseismic)</td><td>V<sub>c</sub></td><td>L/T</td><td>m.s<sup>-1</sup></td><td>3 × 10<sup>-3</sup></td><td>1.5</td><td>2 × 10<sup>-3</sup></td></tr><tr><td>Rupture velocity (mean)</td><td>V<sub>r</sub></td><td>L/T</td><td>m.s<sup>-1</sup></td><td>>10</td><td>3 × 10<sup>3</sup></td><td>3 × 10<sup>-3</sup></td></tr><tr><td>Coseismic D<sub>max</sub>/D<sub>mean</sub></td><td>γD</td><td>—</td><td>—</td><td>~1.7</td><td>~1.8</td><td>1</td></tr><tr><td>Time (interseismic duration)</td><td>T<sub>i</sub></td><td>T</td><td>s</td><td>100</td><td>4 × 10<sup>10</sup></td><td>2.5 × 10<sup>-9</sup></td></tr><tr><td>Time (coseismic duration)</td><td>T<sub>c</sub></td><td>T</td><td>s</td><td>0.05</td><td>25</td><td>2 × 10<sup>-3</sup></td></tr><tr><td><strong>Energetic</strong></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Seismic moment</td><td>Mo</td><td>ML<sup>2</sup>/T<sup>2</sup></td><td>N.m</td><td>2.4 × 10<sup>-1</sup></td><td>1.4 × 10<sup>20</sup></td><td>2 × 10<sup>-21</sup></td></tr><tr><td>Coseismic stress drop</td><td>Δτ</td><td>M/LT<sup>2</sup></td><td>Pa</td><td>~20</td><td>3 × 10<sup>6</sup></td><td>6.6 × 10<sup>-6</sup></td></tr></tbody></table>

nature. In the model, (a-b) is estimated to -0.017 (Caniven et al, 2015), comparable to usual estimates for rock properties providing velocity-weakening conditions favoring stick-slip frictional behavior (e.g., Scholz, 1998). We estimate $D_c$ to $3 \times 10^{-5}$ m, upscaling to ~10 m, thus in the common determined range for earthquakes (0.1-10 m) (Fukuyama & Suzuki, 2016). Finally, we determined $L_c$ to $10^{-3}$ m, upscaling to 300 m, therefore in the usual estimates of earthquake nucleation size (100-1000 m) (Kaneko et al., 2016; Ohnaka, 2013).

Note that the density of the used foam plates is low ($\rho$ = 40 kg.m$^{-3}$) compared to the expected density of seismogenic UP and subducting oceanic crust (3,000 kg.m$^{-3}$). This may induce model scaling inconsistencies (stress ratio) by lowering the lithostatic and normal stress along the subduction interface. To overcome this effect, we adjusted the
average lithostatic stress σᵢ along the interface (~300 Pa) to fit common natural values (~100 MPa) using L*.

Kinematic scaling of the model follows the double time-scale definition proposed by Rosenau et al. (2009) and Corbi et al. (2011). This implies considering separately coseismic and interseismic temporal scaling (Caniven et al., 2015). During the interseismic phase, as viscous forces dominate inertial forces in the silicone wedge, the interseismic temporal factor (Ti*) is determined using the ratio between η* and σ*. We find Ti* = 2.5 × 10⁻⁹ (1 s upscales to ~10 years). Assuming that the Froude number must be similar between model and nature, the coseismic temporal factor (Tc*) corresponds to the root square of L*. We obtained Tc* = 2 × 10⁻³ (1 s upscales to ~500 s). Interseismic and coseismic velocity scaling factors are then determined as the ratios between L* and their respective temporal scaling factor. As a result, during the interseismic period, 5 mm/h upscales to 3 cm/yr and, during the coseismic stage, 3 cm/s upscales to 1.5 m/s.

## 2.3. Model Metrology and Data Processing

Model kinematic evolution is quantified using a high-performance monitoring system, based on the sub-pixel spectral correlation of high-resolution digital pictures (Canon 5D MarkII camera, 25 megapixels resolution) and a FullHD SONY video camera. This device enables very accurate measurements of model deformations in side view with a spatial resolution in the range of 1–6 mm and a precision in the order of a few microns. Using L*, this measurement density can be compared to a geodetic network consisting of ~2 km spaced GPS stations. High-resolution images and measurements are performed each 5–6 s using the CANON5D camera and each 1/25s with the SONY video camera. Note that, because video data storage exceeds hundreds of Giga-bytes per minute, only some portions of each experiment are monitored using the video camera. The kinematic evolution of the model is tracked by recording surface deformation, determining mechanical, and seismic coupling, analyzing the spatiotemporal organization of slip events, and postseismic relaxation stages. In addition, to provide live-detection of coseismic slip events during the experiments, a strain gauge (Figure 2) records the force applied to the backstop, the lab quakes being often associated with small (<0.4 N) but detectable force drops (see Figure S2 for an example of force data record). This live monitoring is also used to detect when the system has reached the seismic cycle steady-state after a first phase of generalized model elastic loading. Because the strain gauge records the integrated effect of various deformation processes acting over the whole model, a deeper quantitative analysis of the force data is yet technically challenging and, then, beyond the scope of this paper.

Using this 2D device, 12 experiments, representing a huge (>3 To) database, have been carried out in the form of several short parametric studies. We investigated more particularly the role of four main parameters on seismic cycle properties: the plate
subduction velocity, the loading (normal stress), the elastic thickness of the slab, and the friction heterogeneity (strong or weak asperities, see Figure S3).

<figure><img src="image_3.png" /><figcaption>**Figure 3.** Kinematics of a typical analog earthquake cycle sequence.</figcaption></figure>

# 3. Results

## 3.1. General Kinematics and Mechanical Model Behavior

Figure 3 illustrates the typical model kinematics observed in our experiments through a short sequence (~45 s) extracted from a 2 h duration experiment monitored every 5 s. Experiments are characterized by a first stage of generalized elastic loading before the fault begins to experience a stick-slip behavior. Typically, hundreds of slip events and related seismic cycles are generated during a single experiment. Along the SZ, the evolution of interplate coupling controls the occurrence of fast coseismic incremental slip events and slow aseismic creep (slip slower or equal to the convergence rate).

After a 1 mn long interseismic period (~600 years), characterized by elastic loading (steps 1 and 2 in Figure 3 and previous steps), the plate interface starts to unlock, and creep triggers along the deeper part of the SZ (step 3). The creep rate accelerates and evolves into a low amplitude slip event (~10 microns) propagating trenchward, up to the central part of the SZ (step 4). Then, the whole subduction interface suddenly ruptures into a megathrust lab quake with a maximum coseismic slip of several hundred microns (step 5). The following postseismic stage shows residual after-slip extending all along the plate interface and viscous relaxation in the rear part of the model (steps 6 and 7). Finally (steps 8 and 9), the plate interface relocks and a new interseismic period begins. The seismic cycle phases are described in more detail in the following sections.
Using fast sampling rate measurements (video camera, 25 Hz), the analysis of model kinematics reveals that lab quakes exhibit a broad variability both in amplitude and slip rate (Figure 4). Slip event types range from aseismic slow slips (duration greater than several tens of seconds) to fast dynamic lab quakes (duration lower than a few seconds associated with an acoustically detectable rupture phase lasting a few tens of milliseconds). Although the events exhibit a continuum, in this study, we distinguish fast events as those that show a distinct slip pulse. The duration of such fast events generally does not exceed 8 s. Longer duration slip events correspond to aseismic and slow slip events.

This analysis reveals that such variability is at least partly controlled by the boundary conditions (tectonic loading rate and imposed initial regional normal stress). More details about the role of boundary conditions on model slip behavior are presented and discussed later.

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

**Figure 4.** Analyses of slip properties and definition of events. (a) Example of lab quake duration versus horizontal displacement (total and coseismic), as measured near the model surface, just above the SZ. Experiments with different boundary conditions (loading rate LR and average applied normal stress NS) are plotted. "Patch" indicates the presence of a frictional heterogeneity (asperity) along the plate interface (green). The rugosity of the fault plane is locally modified by pasting granular material (see Figure S3). The yellow shaded background indicates a continuum from fast to slow-slip events. In this study, event durations shorter than 8 s with an impulsive phase are defining as fast events. The red and blue curves provide examples of measurements discriminating slow (no-impulsive phase) and fast (impulsive phase) events. (b) Details of the measurement procedure.

## 3.2. Comparisons of the Model to Natural Examples

In this section, the analog model results are compared to geophysical data sets coming from different regions and selected for their quality and type to make the comparisons as relevant as possible. Our main selection criteria of each data set require (1) a sufficiently
high density to cover the entire subduction zone and (2) a relevant distinction from the other seismic cycle phases. For example, although we compare the analog coseismic stage to the Tohoku-Oki earthquake data, we believe that the region is not suitable to compare the interseismic stage because (1) seafloor measurements over the overriding plate are lacking and (2) according to several studies, the nucleation of the Tohoku-Oki event likely occurred on a decadal timescale (Ito et al., 2013; Kato et al., 2012; Meng et al., 2011) which may have altered the interseismic stage.

### 3.2.1. Interseismic Phase

During the interseismic period (Figure 5), the friction along the SZ increases the mechanical coupling between the two plates, inducing a partial or a total locking of the subduction interface. Consequently, slab continuous landward motion forces the UP to move in the same direction. UP horizontal displacements reach a maximum, equal to the slab subduction velocity, near the trench, then linearly decrease toward zero landward, indicating that elastic deformation is accumulating. Interestingly, despite the relative simplicity of the analog model, the interseismic phases exhibit a notable kinematic variability. Episodic slip events are commonly observed along the SZ, except when the subduction interface is fully locked. These local low-velocity trenchward slip events can be interpreted as transitory aseismic creep or, even, as slow-slip lab quakes (step 4 in Figure 3).

Figure 5 details a typical case of a high interseismic interplate coupling. Model kinematics is analyzed from displacement field measurements and displacement profiles. The experimental data are compared to interseismic geodetic measurements (Figures 5e and 5f) acquired across the Vanuatu and Nankai subduction zones (Aoki & Scholz, 2009; Bergeot et al., 2009). We selected the Vanuatu subduction zone since the numerous islands, extending from the trench along a 200 km transversal transect, have been monitored for 15 years by GPS stations (Bergeot et al., 2009). Such a data set is of great interest because islands are generally sparsely distributed along subduction zones, strongly limiting the development of dense geodetic networks. One can notice that in both model and Vanuatu's natural case, the horizontal and vertical surface motions of the overriding plate exhibit similar kinematics. The horizontal surface displacement decreases landward as it is theoretically expected for interseismic elastic loading. The vertical component shows subsidence above the upper part of the SZ and evidence of uplift above its deeper portion, close to the SUTZ. These patterns appear to be controlled by the evolution of the interplate mechanical coupling along the SZ. The interplate coupling (C) can be defined as the ratio between the overriding landward plate velocity (Rover) and the subducting plate velocity (Rslab) along the SZ such as $C = \text{Rover}/\text{Rslab}$; with $C = 0$ for a fully unlocked plate interface and $C = 1$ for a fully locked one.

In the experiments, the variation of C along the SZ can be directly determined since model kinematics is monitored in side-view/cross-section (Figure 5g). In Nature, the interplate coupling is quantified by inverting geodetic data (GPS and InSAR) generally situated at far
distances from the SZ and concentrated inland rather than being uniformly distributed. The poor geometry of the geodetic network leads to significant uncertainty on the quantification of C along the SZ (e.g., Scholz and Campos, 2012).

INTERSEISMIC PHASE : A Strong Interplate Coupling Case

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

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

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

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

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

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

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

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

**Figure 5.** The analog interseismic phase compared to geodetic measurements. (a) Displacement vector field; (b) displacement amplitude field; (c and d) displacement components; (e and f) surface profiles of the model and in Natural cases, respectively; (g and h) slip distribution along the interface and mechanical coupling compare with Nankai case in (i). The uncertainty of experimental data is 2 microns. SUTZ, Stable-Unstable Transition Zone.
For the considered interseismic snapshot example, C is equal to 0.3 close to the (SUTZ) and it increases up to 0.8 at half-distance along the SZ, reaching 0.9 near the trench. As a comparison with a natural example, Figure 5h shows the variations of C along the SZ for the Nankai subduction zone, deduced from geodetic data inversions (Aoki & Scholz, 2009). A similar pattern between the analog model and nature can be observed with a steep increase of C starting from the SUTZ region up to the middle part of the SZ. C continues to increase but at a slower rate to reach 0.9 close to the trench.

<figure><img src="image_13.png" /><figcaption>**Figure 6.** The analog coseismic phase and geophysical comparisons. (a) Displacement vector field; (b) displacement amplitude field; (c and d) displacement components; (e and f) surface profiles of the model and in Natural cases, respectively; (g and h) slip distribution along the interface compare with the Tohoku-Okı case in (i). The uncertainty of experimental data is 2 microns. SUTZ, Stable-Unstable Transition Zone.</figcaption></figure>
### 3.2.2. Coseismic Phase

Figure 6 shows an example of one of the largest coseismic slip events that typically occur during the experiments. The incremental displacement field components ($D_{hor}$ and $D_{ver}$) (Figure 6c) and associated profiles (Figure 6e) exhibit an asymmetrical kinematic pattern. UP Horizontal displacements increase rapidly toward the trench, reaching a maximum value, $D_{hor}$ = 600 µm, located about 75 mm landward of the trench. The vertical surface displacement component shows a subsidence of ($D_{ver}$ = -25 µm, upscales to the range -2; -8 m), located trenchward of the SUTZ and above the deeper part of the SZ, and a two-step uplift pattern above the upper part of the SZ. The uplift reaches a first step above the central part of the SZ ($D_{ver}$ = 100 µm) and it increases again significantly toward the trench ($D_{ver}$ = 200 µm). This pattern arises mainly from the interaction of several processes; (1) the deformation induced by the elastic rebound that compacts vertically the UP (subsidence) while it is stretched horizontally, (2) UP trenchward displacements along the 10°-11° dipping plate interface that induce positive vertical motion (uplift), (3) the surface free edge effect which amplifies the uplift displacement near the trench.

These experimental data can be interestingly compared to the 2011 $M_w$ = 9 Tohoku-Oki earthquake (Figure 6f) considered as the best instrumentally monitored large event. Despite the lack of geodetic measurements of seafloor displacements less than 50 km away from the trench, the selected data set here presents

the advantage to use differential bathymetry before and after the earthquake, providing relevant estimates of the near-trench coseismic motion (Fujimoto, 2014; linuma et al., 2012; Sun et al., 2017). Both horizontal and vertical geodetic patterns show significant similarities with the analog model but also relevant quantitative matches when using the scaling factor $L^*$. One interesting feature concerns the shift in the location of the horizontal and vertical maximum surface motions. In the analog model, this shift is about 45 mm. Using $L^*$, it corresponds to about 10–15 km in nature, close to the measured value of about 15 km in the Tohoku-Oki earthquake case. Again, the free edge effect at the trench can explain this pattern, involving a maximum vertical displacement where the UP tip locally straightens. This contributes to a local reduction of the horizontal displacement at the trench and the shift of its maximum value slightly landward.

The coseismic slip distribution along the SZ displays strong similarities too (Figures 6g and 6h). First, for the Tohoku-Oki earthquake, the coseismic slip distribution exhibits an arctangent shape with the steepest part extending from 150 km back to the trench. Using $L^*$, it corresponds to a trench-distance of about 600 mm in the analog model which exactly corresponds to the region where the coseismic slip strongly increases toward the trench in the model. Second, the coseismic slip gradients are very similar, ranging between $10^{-4}$ and $10^{-3}$. Finally, the coseismic slips reach the maximum value over a taper part near the trench, over a distance of 40 km to the Japanese Trench (Figure 6h) that correctly upscales to the observed 150 mm in the analog model (Figure 6g).
### 3.2.3. Postseismic Phase

Following moderate to large coseismic slip events, postseismic deformation can be observed and analyzed. After megathrust lab quakes, it remains detectable up to ~20 s after the coseismic phase (~several hundreds of years in nature) (Figures 4 and 7 right). The increase of mechanical interplate coupling and the relocking of the SZ occur gradually, via a specific displacement pattern similar to the one following natural subduction earthquakes: During the early stage of the postseismic phase (postseismic 1), interplate coupling remains null and the model kinematics is dominated by low amplitude afterslip motion along the SZ (~tens of microns). Indeed, as for the coseismic stage, surface horizontal displacement kinematics still indicates a general trenchward motion of the whole UP and maximum slip values (~35 µm) located close to the trench (Figures 7a and 8a). Vertical displacements reveal also two main zones of uplift, one situated at the tip of the UP and a second one above the SUTZ (Figures 7 and 8a). This pattern is comparable to the postseismic geodetic cumulative measurements acquired after the 2010 Mw 8.8 Maule earthquake (Peña et al., 2019) (Figure 8b). We selected this natural example because, unlike the Tohoku event case, the ruptured subduction zone is very close to the South American coast where available geodetic stations are located (Klein et al., 2016; Peña et al., 2019). This configuration provides better constraints to support the relevance of postseismic models. In the analog model, the general displacement pattern is explained by the delayed mechanical response of the coupled silicone layer to the coseismic elastic traction. Indeed, during the early postseismic stage, the UP can still expand horizontally between the trench and the SUTZ because afterslip occurs along the SZ. As a consequence, the maximum displacement values remain at the trench, although with one order of lower amplitude compare to the main coseismic slip. At the same time, between the backstop and the SUTZ, ductile shearing occurs at the top of the silicone wedge and now allows releasing tractions in the overlying part of the UP. During this postseismic stage 1, the vector field above the silicone wedge evidences a horizontal displacement localized at the SUTZ zone (Figure 7). This delayed strain release transfer through the UP shifts the coseismic maximum vertical displacement value at the model surface from the trench to the SUTZ point, causing the second displacement peak.

This general pattern evolves through the following stage (Postseismic 2), decreasing by one order in amplitude. The direction of some displacement vectors of the UP reverses along the deepest part of the interface, meaning that the SZ starts to relock in this region. This feature is illustrated in Figure 8d by a time-space map showing the evolution of the interplate mechanical coupling C along the plate interface determined by the ratio between the upper and the lower plate displacements (C = Rover/Rslab).

We find that the fault relocking starts during the postseismic stage 2 at the downdip part of the SZ with a coupling value that reaches about 0.5. Figure 8c illustrates the cumulative slip along the interface through this postseismic sequence and confirms the lowest amount of slip along the fault to the postseismic stage 2 close to the SUTZ.
<figure><img src="image_14.png" /><figcaption>**Figure 7.** Kinematic of the postseismic phase following a large coseismic event in the model. The labels A, B, C, D (stage 3) delimit the quadrilateral moving area mentioned in the text.</figcaption></figure>

During postseismic stage 3, the tip region of the UP starts subsiding and some of the displacement vectors into this area show a landward motion (Figure 7). This is consistent with an ongoing relocking of the plate interface evidenced by an interplate coupling greater than 0.5. Note that while the upper part of the SZ relocks, the previously downdip relocked region unlocks again through a slow-slip event. Indeed, from the backstop to the middle part of the SZ, the vector field still indicates a global slow motion of the UP toward the trench. This is expressed more clearly along the foam base in contact with the silicone layer (Figure 7— Postseismic 3). The displacement pattern in this region exhibits a quadrilateral moving area embedded between the labeled points A, B, C, and D. This pattern suggests that ductile shearing drives the displacements observed along the base of the UP in contact with the silicone wedge. The general displacement pattern reflects a contraction in the central part of the UP resulting from the combination of (1) the ongoing SZ relocking and (2) the delayed visco-elastic relaxation of the rear part of the model (Govers et al., 2018; Shi et al., 2020; Sobolev & Muldashev, 2017; Wang et al., 2012). This
kinematic pattern highlights the key role of the ductile layer in the kinematic of an active tectonic system during the postseismic stages.

The displacement pattern of the following stage (sub-interseismic stage 4) indicates the SZ tends to relock along its entire length with an interplate coupling up to ~0.75 (Figure 8). The lowest coupling is measured close to the SUTZ with a value of ~0.25. At the surface of the UP, the signal becomes almost similar to that of the interseismic phase (Figure 7a). From this stage, the following phases show a general increase of interplate coupling with values reaching close to 1 in the central region. The cumulative slip along the SZ (Figure 8c) only shows small displacements occurring close to the SUTZ while the rest of the plate interface remains mostly locked. We consider this stage as the end of the postseismic phase and the beginning of a new interseismic phase.

POSTSEISMIC PHASE : A typical sequence

<figure><img src="image_15.png" /><figcaption>**Figure 8.** Kinematics analyses of the postseismic sequence of Figure 7 following a large coseismic event in the model and comparison with geodetic data. (a and b) Surface displacement profiles; (c) cumulative slip along the interface in the model and (d) evolution of the associated interplate mechanical coupling. The uncertainty of experimental data is 2 microns. SUTZ, Stable-Unstable Transition Zone.</figcaption></figure>
### 3.3. Influence of Boundary Conditions

In our previous analog models investigating the strike-slip fault seismic cycle, the tectonic loading rate appeared as a critical boundary condition controlling model kinematics and the slip properties of the SZ (Caniven, 2014; Caniven et al., 2013, 2015). We performed a similar parametric study using the subduction model configuration. Figure 9 illustrates the comparison of cumulative slip along the SZ between two experiments performed at two different loading velocities (0.6 and 1.26 cm/h). Additional statistical analyses of these two experiments (maximum slip distribution of events) are provided in supporting information (Figure S4).

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

**Figure 9.** Cumulative slip along the subduction interface compared between (a) slow (0.6 cm/h = 1.7 microns/s) and (b) faster (1.26 cm/h = 3.5 microns/s) loading rate cases. Extracts from 3h long experiments. All black curves are separated by 5 s. Red color highlights largest fault slip events. (c and d) Examples of cumulative strain $\epsilon_{xx}$ (horizontal X axis component) during the interseismic period. For the two experiments, the extracted period corresponds to a cumulative displacement of 2 mm imposed by the general loading. For the slow experiment (SUB2D4), a large event occurs immediately after this period.
For the slow experiment (Figure 9a), the SZ exhibits very large characteristic megathrust like events with maximum coseismic slip reaching 600 µm (see an example in Figure 6). We also observed that the amount of cumulative slip in the lower part of the SZ was slightly increasing before the major slip events meaning the plate interface was gradually unlocking. This slow and progressive unlocking initiates during the interseismic period close to the SUTZ then, propagates up to the middle part of the seismogenic region before the triggering of a large coseismic slip instability. Between two coseismic cycles, most of the plate interface slip occurs between the SUTZ zone and the middle part of the SZ. The only few slip events recorded upward corresponds to postseismic slip (after-slip phase). The upper part of the SZ is thus strongly coupled ($C > 0.9$) during the interseismic periods leading to the storage of a large amount of elastic strain energy in the surrounding media (Figure 9c).

For a higher loading rate (Figure 9b), the fault slip behavior is mainly dominated by aseismic creep and low interplate frictional coupling. The kinematics of the interseismic stages reveals a gradual increase of the cumulative slip along the entire SZ and the occurrence of small slip-events, meaning that the fault never locks. Therefore, the amount of stored elastic deformation is considerably minimized during the interseismic phase (Figure 9d).

As evidenced in the strike-slip analog model of Caniven et al. (2015), this suggests that the fault behavior mainly depends on the visco-elastic coupling between the upper foam plate and the silicone wedge. Because of the strain-rate-dependent rheology of the viscous material, at low loading velocity, the viscosity of the silicone wedge decreases as well as its mechanical coupling with the base of the upper foam plate where a viscous shearing develops. This weak coupling favors the accumulation of elastic strain in the freely deformed UP during longer interseismic periods released then, by larger coseismic events (Figure 9c). Reversely, fast loading velocity increases the viscous strength of the silicone wedge as well as its mechanical coupling with the base of the upper foam plate. This strong coupling minimizes the accumulation of elastic strain in the UP and forces the frictional interface to slip at a velocity close to the tectonic loading rate (Figure 9d).

# 4. Discussion

## 4.1. Model Limitations and Potentialities

As for any model, the extrapolation of our experimental results to nature faces some limitations related to the complexity of geological processes.

The analog model rheology was designed to simulate the main physical and mechanical properties of an idealized subduction zone. Thermal and pressure gradients, controlling phase changes, and metamorphism processes are not reproduced. The fault interface geometry remains fixed during the experiments and no off-fault brittle deformation is
allowed (Sylvester, 1988). Fluid pressure and fault damage processes, known to influence fault slip stability (Bizzarri, 2009; Faulkner & Rutter, 2001), are not simulated. Interestingly, despite the lack of fluids or phase change processes, aseismic creep, and slow-slip events occur, questioning the specific role of such processes on fault slip dynamics. Model scaling is perfectible. The density contrast between elastic foam plates and the analog ductile mantle wedge is out of scale although we compensated this limitation by adjusting a scaled lithostatic load.

Finally, due to technical considerations, high-resolution kinematic measurements are acquired at a 5–6 s time interval. Consequently, the resulting strain fields correspond to the sum of all the slip events and elastic deformations that occurred during this period. This feature only affects the smallest lab quakes that can be masked by the interseismic reloading included in the same measurement time interval.

Consequently, the interpretation of analog model results should be handled with care by focusing mainly on first order tendencies.

Nevertheless, the model reveals significant similitudes with natural cases in regards to the kinematics and the mechanics of the system.

By using visco-elasto-plastic stratified rheology, the model is capable to simulate all the seismic cycle phases: the interseismic, the preseismic, the coseismic, and the postseismic stages, including unstable and stable fault-slip behaviors. The postseismic stage is correctly simulated by including realistic after-slip and visco-elastic kinematics and mechanical behavior, representing a real advance in the analog modeling of the subduction zone seismic cycle (Rosenau et al., 2017).

The visco-elastic mechanical coupling between the elastic overriding foam plate and the ductile silicone wedge highlights interesting deformation processes that could not be evidenced with uniform elastic model rheology. The change in this visco-elastic coupling depends on the loading rate and conduct to drastic changes in the seismic cycle dynamics and the earthquake energy that are discussed in the following section.

## 4.2. Role of the Visco-Elastic Coupling Between the Layers

In our strike-slip and subduction analog models, lab quake size distribution, interplate coupling, and more generally the mechanical behavior of the SZ are strongly modulated by the imposed loading rate. The lower is the loading rate, the lower is the amount of aseismic fault creep, and the larger are the lab quakes. One should note that the amplitude of slip events also depends on the general applied normal load on the plate interface (Figure 4). One possible explanation for this relationship, linking loading rate and lab quake coseismic slip, could be that the static friction along the SZ is time-dependent as already evidenced for rocks and metals (Dieterich, 1972). Frictional sliding experiments, using two rock samples, demonstrated that the static friction coefficient evolves with the
logarithm of the time duration while the two blocks are held in stationary contact before sliding (Dieterich, 1972; Li et al., 2011; Marone, 1998).

<figure><img src="image_17.png" /><figcaption>**Figure 10.** Schematic illustration of the crust-mantle visco-elastic coupling effects on the earthquake size and slip behavior. (a) Weak visco-elastic coupling case. In our model, this corresponds to slow convergent experiments. (b) Strong visco-elastic coupling case. In our model, this corresponds to faster convergent experiments. The size of the yellow stars illustrates the size of earthquakes that can be generated. Gray arrows represent the elastic strain and the red color schematically highlights the volume of elastically deformed crust during the interseismic loading.</figcaption></figure>

To investigate a possible time-dependent effect of the experimental SZ frictional properties, we performed tribological measurements using a specific linear shear testing device (principle similar to the one used by Hencher & Richards, 2015). Results revealed that, as already observed with the strike-slip analog model (Caniven, 2014; Caniven et al., 2015), this process alone is not enough to explain the observed lab quake amplitude variability. Furthermore, loading rate tests conducted without the visco-elastic coupling (i.e., without the silicone layer), confirmed this conclusion.

Based on numerical simulation integrating rate-and-state friction laws, Cao and Aki (1986) suggested it would be necessary to divide the elastic loading rate by 1000 to get a coseismic stress-drop increase of just 30%. The time-dependence of the frictional properties appears, then, not sufficient to explain the observations both in nature and our
model. Several studies proposed supplemental healing processes, involving pore-fluid pressure, to explain observed stress-drop variability (Dieterich, 1972; Marone, 1998; Muhuri et al., 2003). However, the absence of fluids in our model rules out this explanation.

Recently, several studies revealed a possible relationship between earthquake magnitude and strain rates by showing that the strongest earthquakes frequently locate in low strain rate regions (Doglioni et al., 2015; Riguzzi et al., 2012). It is expected that high strain rate regions host strongly coupled upper and lower crusts, favoring quasi-rigid continuous displacements along faults rather than stick-slip behavior. In such a context, the storage of elastic deformation is low, and only small to moderate earthquakes can trigger. Doglioni et al. (2015) discussed the so-called earthquake energy accumulated during the interseismic period above the SUTZ zone as a function of the volume of the upper crust involved in the seismic cycle.

In our subduction analog model, we propose that the elastically deformed volume of the upper-crust depends on the visco-elastic mechanical coupling that is modulated by the tectonic loading rate. Figure 10 illustrates schematically the proposed relationships between crust-mantle visco-elastic coupling and earthquake behavior. In our experiments, the volume of the elastic energy reservoir increases as the loading rate and the visco-elastic coupling decrease, lengthening the interseismic periods and increasing seismic event magnitudes (Figure 10a). Reversely, increasing the loading rate and so, the visco-elastic coupling tends to minimize the volume of stored elastic strain energy in the UP favoring, then, the occurrence of aseismic fault creep and small to moderate lab quakes (Figure 10b). Note that in our experiments, the visco-elastic coupling is essentially modulated by the strain rate. In nature, it depends

on additional parameters such as pressure-temperature conditions, and fluid content controlling the rheological properties of rocks (e.g., Evans & Kohlstedt, 1995; Hyndman & Peacock, 2003; Karato & Wu, 1993; Keefner et al., 2011). Heimpel (2006) investigated the control of the visco-elastic coupling parameter on the earthquake size with numerical simulations. Interestingly, the conclusion of this work states that the earthquake magnitude can be increased when the seismogenic crust is decoupled from the deeper layers, which is consistent with our inferences. Our results suggest that, in addition to the intrinsic frictional fault processes, the whole rheology of a tectonic system and the involved rock's volume around the fault plane must be considered in the analysis of the slip behavior.

### 4.3. Precursory Stages

In the last decades, the development of geophysical investigations favored the discovery of slow-slip events preceding great coseismic ruptures (Kato et al., 2012; Ruiz et al., 2014). Precursory slow-slip earthquakes have been stated before the outbreak of the 2011 $M_w = 9.1$ Tohoku-Oki (Ito et al., 2013; Kato et al., 2012; Meng et al., 2011) and the 2014 $M_w$
= 8.1 Iquique (Kato & Nakagawa, 2014; Ruiz et al., 2014) earthquakes. For the two cases, data revealed that large slow-slip events started nucleating on a decadal timescale before the occurrence of the mainshocks.

Using the whole experimental data set, we evidenced that most of the large experimental lab quakes are also preceded by an early gradual unlocking of the deeper part of the SZ on the form of accelerated creep and precursory slow-slips events (Figures 9 and 10). Interestingly, these precursory stages can also be detected at the surface of the model with a more discreet signal, depending mainly on the location of the measurement site. This observation could potentially have a significant impact on the implementation and interpretation of geodetic data.

Figure 11 shows the comparison between model surface horizontal displacements measured at three distinct sites along the UP surface (Figure 11a) and the evolution of the instantaneous slip along the SZ (Figure 11b). The three sites are situated close to the backstop (pt1, red), a little bit trenchward above the SUTZ zone (pt2, green), near the trench (pt3, blue).

In many cases, during the interseismic period, the site 2 (green) recorded an inversion of motion several tens of seconds before the occurrence of a large coseismic event. These displacements recorded at the second site are related to precursory slow-slip events occurring along the deeper portion of the SZ (Figure 11b). Note that some isolated slow-slip events can also be evidenced during the interseismic period between large coseismic events.

Precursory slip behavior has been also investigated experimentally (e.g., Ohnaka, 1992; Rubinstein et al., 2007). The onset of the instability during this predynamic phase has been described by rate-and-state dependent friction (Dieterich, 1992; Marone, 1998; Yamashita & Ohnaka, 1991). To provide a more physical explanation about the mechanical processes controlling slip and its stability, new laboratory and numerical models have been developed. Experiments of frictional sliding with rock samples (Harbord et al., 2017) and analog seismic cycle models (Caniven et al., 2017; Corbi et al., 2017; Rosenau et al., 2019) suggested that fault roughness and related spatial distribution of normal stress at the local scale significantly act on the frictional behavior. These explanations are also supported by numerical models investigating the role of geometrical fault complexity on offset distribution and nucleation processes (Blank & Morgan, 2019; Fournier & Morgan, 2012; Romanet et al., 2018; Tal et al., 2018; Zielke et al., 2017).

The spatial organization of contacts related to variations in fault surface morphology is proposed to be the controlling factor of the fault slip stability. Especially, Blank and Morgan (2019) provided an interesting explanation of precursory slow-slip using Discrete Element Method numerical modeling to explore relationships between slow-slip events and dynamic earthquakes. This study revealed gradual and localized stress drops, occurring before the fast slip stage, attributed to the fault dilation and the associated progressive unlatching of geometrical asperities.
<figure><img src="image_18.png" /><figcaption>**Figure 11.** (a) Surface deformation at several sites. The green site reveals interesting geodetic precursors. The blue site indicates progressive interface unlocking and relocking, respectively before and after large events. SUTZ, Stable-Unstable Transition Zone. (b) Spatio-temporal evolution of the instantaneous slip along the interface during the sequence contained within the dashed frame in (a). The blue curve on the right side is the synchronized corresponding cumulative horizontal displacement at the blue site.</figcaption></figure>

In our analog model, the fault surface morphology is characterized by heterogeneous micro-topographic variations with wavelengths of 1–3 mm and amplitudes in the range of 0.1–1 mm (see details of the roughness analysis in Caniven et al., 2015). Similarly to Blank and Morgan (2019), considering such a complexity of the analog fault surface in our
model and the lack of pore fluid interaction, we believe that the roughness can be considered as the key dominant factor in the control of fast versus slow slip and the broad variability of slip-events. The fact that geodetic precursors recorded at the surface are not always systematic (both in the model and in nature) is probably due to the complexity of strain patterns at depth, masked, and averaged through the media. This also suggests that during this precursory stage, the strain and stress change patterns along the fault could

bring out significant constraints on earthquake early forecast approaches (Caniven et al., 2019).

**5. Conclusion**

We have developed a new analog model to investigate experimentally subduction zone seismic cycle and megathrust earthquake dynamics. Our model integrates realistic boundary conditions and stratified visco-elasto-plastic rheology to simulate the overriding plate and subducting slab mechanical behaviors. The model succeeds in replicating the strain field of all seismic cycle phases. The coseismic phase is characterized by dynamic fault ruptures presenting wide ranges of magnitudes, slip velocities, and locations along the SZ. Episodic aseismic creep and slow-slip events occur during the main interseismic phase depending on the amount of frictional coupling along the SZ.

As in Nature, quasi-instantaneous coseismic slips (duration < to 200 milliseconds) up to very slow-slip events (duration > tens of seconds) are observed in our experiments. The simulated postseismic deformation is also consistent with natural observations. It starts with after-slip propagating along the SZ immediately followed by visco-elastic relaxation, affecting the portion of the UP resting on the ductile mantle wedge, induced by slow slip horizontal shearing at the foam-silicone transition. The visco-elastic coupling between crust and mantle induces a delayed release of elastic stress in the UP providing a realistic postseismic surface strain pattern consistent with available geodetic measurements.

Based on these results, we believe that model scaling is acceptable to extrapolate analog results to natural cases. First results exhibit interesting relationships between imposed kinematical boundary conditions and fault slip behavior. Especially, the tectonic loading rate appears as a critical parameter controlling seismic cycle properties and slip stability in the model. Our results also highlight the fundamental role of the visco-elastic mechanical coupling between the different rheologic layers. We propose the visco-elastic coupling controls the volume of the elastic upper crust involved in the seismic cycle, and thus, the elastic deformation storable during the interseismic phase.

A fast loading rate strengthens the visco-elastic coupling, induces a lower amount of accumulated elastic deformation mainly released, then, by aseismic creep, slow slip-events, and small to moderate lab quakes. A slow loading rate weakens the visco-elastic coupling, induces a larger amount of accumulated elastic deformation, released mostly by
large megathrust coseismic events. Another interesting point is the evidence of precursory slow-slip events before the mainshocks in our experiments. We argue that the broad variety of events and slip behavior might be related to the fault roughness complexity rather than to fluid processes, these latter being absent in the model. These results support the scientific interest of the analog modeling approach we developed that could afford new insights into the mechanics of the subduction seismic cycle and associated earthquake hazard.

## Acknowledgments

This project has been funded by CNRS-INSU (TelluS-Aleas). We thank the editor Rachel Abercrombie and the associate editor, Michel Cooke and other anonymous reviewers for helping to improve the manuscript, C. Romano for highly effective technical assistance, Julia Morgan, and David Blank for fruitful discussions.

## Data Availability Statement

All data sets and materials underlying this study are published open access in Caniven and Domínguez (2020) (DOI:10.5281/zenodo.3820381) and additional details are provided in Supporting Information.

## References

* Anooshehpoor, R., & Brune, J. (2003). Rupture directivity in a foam rubber physical model. American Geophysical Union, Fall Meeting 2003.

* Aoki, Y., & Scholz, C. H. (2009). Imaging interseismic locking at the Nankai subduction zone, southwest Japan. In S. Lallemand, & F. Funiciello (Eds.), *Subduction zone geodynamics, frontiers in Earth sciences* (pp. 159–171). Berlin, Heidelberg: Springer. https://doi.org/10.1007/978-3-540-87974-9_9

* Bergeot, N., Bouin, M. N., Diament, M., Pelletier, B., Régnier, M., Calmant, S., & Ballu, V. (2009). Horizontal and vertical interseismic velocity fields in the Vanuatu subduction zone from GPS measurements: Evidence for a central Vanuatu locked zone. *Journal of Geophysical Research, 114*, B06405. https://doi.org/10.1029/2007JB005249

* Bizzarri, A. (2009). What does control earthquake ruptures and dynamic faulting? a review of different competing mechanisms. In S. Vinciguerra, & Y. Bernabé (Eds.), *Rock physics and natural hazards, pageoph topical volumes* (pp. 741–776). Basel, Switzerland: Birkhäuser.

* Blank, D. G., & Morgan, J. K. (2019). Precursory Stress changes and fault dilation lead to fault rupture: Insights from discrete element simulations. *Geophysical Research Letters, 46*, 3180–3188. https://doi.org/10.1029/2018GL081007

* Brune, J. N. (1973). Earthquake modeling by stick-slip along precut surfaces in stressed foam rubber. *Bulletin of the Seismological Society of America, 63*(6–1), 2105–2119.

* Brune, J. N., & Anooshehpoor, A. (1998). A physical model of the effect of a shallow weak layer on strong ground motion for strike-slip ruptures. *Bulletin of the Seismological Society of America, 88*(4), 1070–1078.
Caniven, Y. (2014). *Kinematics and mechanics of strike-slip faults at the seismic cycle time-scale: Insights from an experimental approach* (Doctoral dissertation): University of Montpellier. Retrieved from http://www.theses.fr/2014MON20057

Caniven, Y., & Dominguez, S. (2020). *Modeling subduction megathrust earthquake cycle: Insights from a visco-elasto-plastic analog model*. Zenodo. http://doi.org/10.5281/zenodo.3820381

Caniven, Y., Dominguez, S., Soliva, R., Cattin, R., Peyret, M., Chéry, J., & Romano, C. (2013). *Dynamics of a strike-slip fault analog model: Effects of the tectonic loading rate*. San Francisco, CA: AGU. https://ui.adsabs.harvard.edu/abs/2013AGUFM.T51C2478C.

Caniven, Y., Dominguez, S., Soliva, R., Cattin, R., Peyret, M., Marchandon, M., et al. (2015). A new multilayered visco-elasto-plastic experimental model to study strike-slip fault seismic cycle. *Tectonics, 34*, 232–264. https://doi.org/10.1002/2014TC003701

Caniven, Y., Dominguez, S., Soliva, R., Peyret, M., Cattin, R., & Maerten, F. (2017). Relationships between along-fault heterogeneous normal stress and fault slip patterns during the seismic cycle: Insights from a strike-slip fault laboratory model. *Earth and Planetary Science Letters, 480*, 147–157. https://doi.org/10.1016/j.epsl.2017.10.009

Caniven, Y., Morgan, J., Blank, D. G., & Cattin, R. (2019). *DEM simulations of earthquake cycles along strike-slip faults: Controls on the slow-slip nucleation of earthquakes*. Presented at the AGU Fall Meeting 2019, AGU. Retrieved from https://ui.adsabs.harvard.edu/abs/2019AGUFM.S21F0572C.

Cao, T., & Aki, K. (1986). Effect of slip rate on stress drop. *Pure and Applied Geophysics, 124*, 515–529. https://doi.org/10.1007/BF00877214

Cobbold, P. R., & Jackson, M. P. A. (1992). Gum rosin (colophony): A suitable material for thermomechanical modelling of the lithosphere. *Tectonophysics, 210*, 255–271. https://doi.org/10.1016/0040-1951(92)90325-Z

Corbi, F., Bedford, J., Sandri, L., Funiciello, F., Gualandi, A., & Rosenau, M. (2020). Predicting imminence of analog megathrust earthquakes with machine learning: Implications for monitoring subduction zones. *Geophysical Research Letters, 47*(7). e2019GL086615. https://doi.org/10.1029/2019GL086615

Corbi, F., Funiciello, F., Brizzi, S., Lallemand, S., & Rosenau, M. (2017). Control of asperities size and spacing on seismic behavior of subduction megathrusts. *Geophysical Research Letters, 44*(16), 8227–8235. https://doi.org/10.1002/2017GL074182

Corbi, F., Funiciello, F., Faccenna, C., Ranalli, G., & Heuret, A. (2011). Seismic variability of subduction thrust faults: Insights from laboratory models. *Journal of Geophysical Research, 116*, B06304. https://doi.org/10.1029/2010JB007993

Corbi, F., Funiciello, F., Moroni, M., van Dinther, Y., Mai, P. M., Dalguer, L. A., & Faccenna, C. (2013). The seismic cycle at subduction thrusts: 1. Insights from laboratory models. *Journal of Geophysical Research-Solid Earth, 118*(4), 1483–1501. https://doi.org/10.1029/2012JB009481

Corbi, F., Sandri, L., Bedford, J., Funiciello, F., Brizzi, S., Rosenau, M., & Lallemand, S. (2019). Machine learning can predict the timing and size of analog earthquakes. *Geophysical Research Letters, 46*(3), 1303–1311. https://doi.org/10.1029/2018GL081251

Day, S. M., Gonzalez, S. H., Anooshehpoor, R., & Brune, J. N. (2008). Scale-model and numerical simulations of near-fault seismic directivity. *Bulletin of the Seismological Society of America, 98*(3), 1186–1206. https://doi.org/10.1785/0120070190

Di Giuseppe, E., Funiciello, F., Corbi, F., Ranalli, G., & Mojoli, G. (2009). Gelatins as rock analogs: A systematic study of their rheological and physical properties. *Tectonophysics, 473*, 391–403. https://doi.org/10.1016/j.tecto.2009.03.012

Dieterich, J. H. (1972). Time-dependent friction in rocks. *Journal of Geophysical Research, 77*, 3690–3697. https://doi.org/10.1029/JB077i020p03690

Dieterich, J. H. (1992). Earthquake nucleation on faults with rate-and state-dependent strength. *Tectonophysics, 211*, 115–134. https://doi.org/10.1016/0040-1951(92)90055-B
Doglioni, C., Barba, S., Eugenio, C., & Riguzzi, F. (2015). Fault on-off versus strain rate and earthquakes energy. *Geoscience Frontiers, 6*, 265–276. https://doi.org/10.1016/j.gsf.2013.12.007

Duesterhoefft, E., Quinteros, J., Oberhänsli, R., Bousquet, R., & de Capitani, C. (2014). Relative impact of mantle densification and eclogitization of slabs on subduction dynamics: A numerical thermodynamic/thermokinematic investigation of metamorphic density evolution. *Tectonophysics*, 637, 20–29. https://doi.org/10.1016/j.tecto.2014.09.009

Evans, B., & Kohlstedt, D. L. (1995). Rheology of rocks. In T. J. Ahrens (Ed.), *Rock physics and phase relations: A handbook of physical constants* (pp. 149–165): AGU.

Fagereng, A., & Sibson, R. (2010). Melange rheology and seismic style. *Geology*, 38, 751–754.

Faulkner, D. R., & Rutter, E. H. (2001). Can the maintenance of overpressured fluids in large strike-slip fault zones explain their apparent weakness? *Geology*, 29, 503–506. https://doi.org/10.1130/0091-7613(2001)029<0503:CTMOOF>2.0.CO;2

Fournier, T., & Morgan, J. (2012). Insights to slip behavior on rough faults using discrete element modeling. *Geophysical Research Letters*, 39, L12304. https://doi.org/10.1029/2012GL051899

Fujimoto, H. (2014). Seafloor geodetic approaches to subduction thrust earthquakes. *Monographs on Environment, Earth and Planets*, 2(2), 23–63. https://doi.org/10.5047/meep.2014.00202.0023

Fukuyama, E., & Suzuki, W. (2016). Near-fault deformation and Dc" during the 2016 Mw7.1 Kumamoto earthquake. *Earth, Planets and Space*, 68, 194. https://doi.org/10.1186/s40623-016-0570-6

Goldfinger, C., Ikeda, Y., Yeats, R. S., & Ren, J. (2013). Superquakes and supercycles. *Seismological Research Letters*, 84, 24–32. https://doi.org/10.1785/0220110135

Govers, R., Furlong, K. P., van de Wiel, L., Herman, M. W., & Broerse, T. (2018). The geodetic signature of the earthquake cycle at subduction zones: Model constraints on the deep processes. *Reviews of Geophysics*, 56(1), 6–49. https://doi.org/10.1002/2017RG000586

Harbord, C. W. A., Nielsen, S. B., Paola, N. D., & Holdsworth, R. E. (2017). Earthquake nucleation on rough faults. *Geology*, 45, 931–934. https://doi.org/10.1130/G39181.1

Hayman, N. W., & Lavier, L. L. (2014). The geologic record of deep episodic tremor and slip. *Geology*, 42(3), 195–198. https://doi.org/10.1130/G34990.1

Heimpel, M. (2006). Earthquake scaling: The effect of a viscoelastic asthenosphere. *Geophysical Journal International*, 166(1), 170–178. https://doi.org/10.1111/j.1365-246X.2006.02815.x

Hencher, S. R., & Richards, L. R. (2015). Assessing the shear strength of rock discontinuities at laboratory and field scales. *Rock Mechanics and Rock Engineering*, 48, 883–905. https://doi.org/10.1007/s00603-014-0633-6

Hubbert, M. K. (1937). Theory of scale models as applied to the study of geologic structures. *Bulletin of the Geological Society of America*, 48, 1459–1520.

Hyndman, R. D., & Peacock, S. M. (2003). Serpentinization of the forearc mantle. *Earth and Planetary Science Letters*, 212(3–4), 417–432.

Hyndman, R. D., Yamano, M., & Oleskevich, D. A. (1997). The seismogenic zone of subduction thrust faults. *Island Arc*, 6, 244–260. https://doi.org/10.1111/j.1440-1738.1997.tb00175.x

Iinuma, T., Hino, R., Kido, M., Inazu, D., Osada, Y., Ito, Y., et al. (2012). Coseismic slip distribution of the 2011 off the Pacific Coast of Tohoku Earthquake (M9.0) refined by means of seafloor geodetic data. *Journal of Geophysical Research*, 117, B07409. https://doi.org/10.1029/2012JB009186

Ito, Y., Hino, R., Kido, M., Fujimoto, H., Osada, Y., Inazu, D., et al. (2013). Episodic slow slip events in the Japan subduction zone before the 2011 Tohoku-Oki earthquake. *Tectonophysics*, 600, 14–26. https://doi.org/10.1016/j.tecto.2012.08.022
Ito, Y., Tsuji, T., Osada, Y., Kido, M., Inazu, D., Hayashi, Y., et al. (2011). Frontal wedge deformation near the source region of the 2011 Tohoku-Oki earthquake. *Geophysical Research Letters*, 38, L00G05. http://dx.doi.org/10.1029/2011GL048355

Kaneko, Y., Nielsen, S. B., & Carpenter, B. M. (2016). The onset of laboratory earthquakes explained by nucleating rupture on a rate-and-state fault. *Journal of Geophysical Research: Solid Earth*, 121, 6071–6091. https://doi.org/10.1002/2016JB013143

Karato, S., & Wu, P. (1993). Rheology of the upper mantle: A synthesis. *Science*, 260(5109), 771–778. https://doi.org/10.1126/science.260.5109.771

Kato, A., & Nakagawa, S. (2014). Multiple slow-slip events during a foreshock sequence of the 2014 Iquique, Chile Mw 8.1 earthquake. *Geophysical Research Letters*, 41, 5420–5427. https://doi.org/10.1002/2014GL061138

Kato, A., Obara, K., Igarashi, T., Tsuruoka, H., Nakagawa, S., & Hirata, N. (2012). Propagation of slow slip leading up to the 2011 Mw 9.0 Tohoku-Oki earthquake. *Science*, 335, 705–708. https://doi.org/10.1126/science.1215141

Keefner, J. W., Mackwell, S. J., Kohlstedt, D. L., & Heidelbach, F. (2011). Dependence of dislocation creep of dunite on oxygen fugacity: implications for viscosity variations in Earth's mantle. *Journal of Geophysical Research*, 116, B05201. https://doi.org/10.1029/2010JB007748

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

Koketsu, K., Yokota, Y., Nishimura, N., Yagi, Y., Miyazaki, S., Satake, K., et al. (2011). A unified source model for the 2011 Tohoku earthquake. *Earth and Planetary Science Letters*, 310, 480–487. https://doi.org/10.1016/j.epsl.2011.09.009

Li, Q., Tullis, T. E., Goldsby, D., & Carpick, R. W. (2011). Frictional ageing from interfacial bonding and the origins of rate and state friction. *Nature*, 480, 233–236. https://doi.org/10.1038/nature10589

Manea, V., & Gurnis, M. (2007). Subduction zone evolution and low viscosity wedges and channels. *Earth and Planetary Science Letters*, 264(1), 22–45. https://doi.org/10.1016/j.epsl.2007.08.030

Marone, C. (1998). Laboratory-derived friction laws and their application to seismic faulting. *Annual Review of Earth and Planetary Sciences*, 26, 643–696. https://doi.org/10.1146/annurev.earth.26.1.643

McBean, K. M., Anderson, J. G., Brune, J. N., & Anooshehpoor, R. (2015). Statistics of ground motions in a foam rubber model of a strike-slip faultstatistics of ground motions in a foam rubber model of a strike-slip fault. *Bulletin of the Seismological Society of America*, 105(3), 1456–1467. https://doi.org/10.1785/0120140276

Meng, L., Inbal, A., & Ampuero, J.-P. (2011). A window into the complexity of the dynamic rupture of the 2011 Mw 9 Tohoku-Oki earthquake. *Geophysical Research Letters*, 38, L00G07. https://doi.org/10.1029/2011GL048118

Muhuri, S. K., Dewers, T. A., Scott, T. E., & Reches, Z. (2003). Interseismic fault strengthening and earthquake-slip instability: Friction or cohesion? *Geology*, 31, 881–884. https://doi.org/10.1130/G19601.1

Ohnaka, M. (1992). Earthquake source nucleation: A physical model for short-term precursors. *Tectonophysics*, 211, 149–178. https://doi.org/10.1016/0040-1951(92)90057-D

Ohnaka, M. (2013). *The physics of rock failure and earthquakes*: Cambridge University Press. https://doi.org/10.1017/CBO9781139342865. Okada, Y. (1985). Surface deformation due to shear and tensile faults in a half-space. *Bulletin of the Seismological Society of America*, 75(4), 1135–1154.

Ozawa, S., Nishimura, T., Suito, H., Kobayashi, T., Tobita, M., & Imakiire, T. (2011). Coseismic and postseismic slip of the 2011 magnitude-9 Tohoku-Oki earthquake. *Nature*, 475, 373–376. https://doi.org/10.1038/nature10227
Pacheco, J. F., Sykes, L. R., & Scholz, C. H. (1993). Nature of seismic coupling along simple plate boundaries of the subduction type. *Journal of Geophysical Research*, 98, 14133–14159. https://doi.org/10.1029/93JB00349

Peña, C., Heidbach, O., Moreno, M., Bedford, J., Ziegler, M., Tassara, A., & Oncken, O. (2019). Role of lower crust in the postseismic deformation of the 2010 Maule earthquake: Insights from a model with power-law rheology. *Pure and Applied Geophysics*, 176, 3913–3928. https://doi.org/10.1007/s00024-018-02090-3

Riguzzi, F., Crespi, M., Devoti, R., Doglioni, C., Pietrantonio, G., & Pisani, A. R. (2012). Geodetic strain rate and earthquake size: New clues for seismic hazard studies. *Physics of the Earth and Planetary Interiors*, 206–207, 67–75. https://doi.org/10.1016/j.pepi.2012.07.005

Romanet, P., Bhat, H. S., Jolivet, R., & Madariaga, R. (2018). Fast and Slow Slip Events Emerge Due to Fault Geometrical Complexity. *Geophysical Research Letters*, 45, 4809–4819. https://doi.org/10.1029/2018GL077579

Rosenau, M., Corbi, F., & Dominguez, S. (2017). Analogue earthquakes and seismic cycles: Experimental modelling across timescales. *Solid Earth*, 8, 597–635. https://doi.org/10.5194/se-8-597-2017

Rosenau, M., Horenko, I., Corbi, F., Rudolf, M., Kornhuber, R., & Oncken, O. (2019). Synchronization of great subduction megathrust earthquakes: Insights from scale model analysis. *Journal of Geophysical Research: Solid Earth*, 124(4), 3646–3661. https://doi.org/10.1029/2018JB016597

Rosenau, M., Lohrmann, J., & Oncken, O. (2009). Shocks in a box: An analogue model of subduction earthquake cycles with application to seismotectonic forearc evolution. *Journal of Geophysical Research*, 114, B01409. https://doi.org/10.1029/2008JB005665

Rubinstein, S. M., Cohen, G., & Fineberg, J. (2007). Dynamics of precursors to frictional sliding. *Physical Review Letters*, 98, 226103. https://doi.org/10.1103/PhysRevLett.98.226103

Rudolf, M., Boutelier, D., Rosenau, M., Schreurs, G., & Oncken, O. (2016). Rheological benchmark of silicone oils used for analog modeling of short- and long-term lithospheric deformation. *Tectonophysics*, 684, 12–22. https://doi.org/10.1016/j.tecto.2015.11.028

Ruiz, S., Metois, M., Fuenzalida, A., Ruiz, J., Leyton, F., Grandin, R., et al. (2014). Intense foreshocks and a slow slip event preceded the 2014 Iquique Mw 8.1 earthquake. *Science*, 345, 1165–1169. https://doi.org/10.1126/science.1256074

Saffer, D. M., & Wallace, L. M. (2015). The frictional, hydrologic, metamorphic and thermal habitat of shallow slow earthquakes. *Nature Geoscience*, 8(8), 594–600. https://doi.org/10.1038/ngeo2490

Scholz, C. H. (1998). Earthquakes and friction laws. *Nature*, 391, 37–42. https://doi.org/10.1038/34097

Scholz, C. H., & Campos, J. (2012). The seismic coupling of subduction zones revisited. *Journal of Geophysical Research*, 117, B05310. https://doi.org/10.1029/2011JB009003

Shi, Q., Barbot, S., Wei, S., Tapponnier, P., Matsuzawa, T., & Shibazaki, B. (2020). Structural control and system-level behavior of the seismic cycle at the Nankai Trough. *Earth, Planets and Space*, 72(1), 27. https://doi.org/10.1186/s40623-020-1145-0

Sobolev, S. V., & Muldashev, I. A. (2017). Modeling seismic cycles of great megathrust earthquakes across the scales with focus at postseismic phase. *Geochemistry, Geophysics, Geosystems*, 18(12), 4387–4408. https://doi.org/10.1002/2017GC007230

Subarya, C., Chlieh, M., Prawirodirdjo, L., Avouac, J.-P., Bock, Y., Sieh, K., et al. (2006). Plate-boundary deformation associated with the great Sumatra–Andaman earthquake. *Nature*, 440, 46–51. https://doi.org/10.1038/nature04522

Sun, T., Wang, K., Fujiwara, T., Kodaira, S., & He, J. (2017). Large fault slip peaking at trench in the 2011 Tohoku-Oki earthquake. *Nature Communications*, 8(1), 14044. https://doi.org/10.1038/ncomms14044

Sylvester, G. (1988). Strike-slip faults. *Geological Society of America Bulletin*, 100, 1666–1703. https://doi.org/10.1007/s00024-08-02097-6
process. *Journal of Geophysical Research: Solid Earth*, 123, 437–456. https://doi.org/
10.1002/2017JB014746

Tong, X., & Lavier, L. L. (2018). Simulation of slip transients and earthquakes in finite thickness shear zones with a plastic formulation. *Nature Communications*, 9(1), 3893. https://doi.org/10.1038/s41467-018-06390-z

Ujiie, K., Saishu, H., Fagereng, Å., Nishiyama, N., Otsubo, M., Masuyama, H., & Kagi, H. (2018). An explanation of episodic tremor and slow slip constrained by crack-seal veins and viscous shear in subduction mélange. *Geophysical Research Letters*, 45, 5371–5379. https://doi.org/10.1029/2018GL078374

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

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

Warren-Smith, E., Fry, B., Wallace, L., Chon, E., Henrys, S., Sheehan, A., et al. (2019). Episodic stress and fluid pressure cycling in subducting oceanic crust during slow slip. *Nature Geoscience*, 12, 475–481. https://doi.org/10.1038/s41561-019-0367-x

Yamashita, T., & Ohnaka, M. (1991). Nucleation process of unstable rupture in the brittle regime: A theoretical approach based on experimentally inferred relations. *Journal of Geophysical Research*, 96, 8351–8367. https://doi.org/10.1029/91JB00106

Zielke, O., Galis, M., & Mai, P. M. (2017). Fault roughness and strength heterogeneity control earthquake size and stress drop. *Geophysical Research Letters*, 44, 777–783. https://doi.org/10.1002/2016GL071700