Synchronization of Great Subduction Megathrust Earthquakes: Insights From Scale Model Analysis

The magnitude of great subduction megathrust earthquakes is controlled mainly by the number of adjacent asperities failing synchronously and the resulting rupture length. Here we investigate experimentally the long‐term recurrence behavior of a pair of asperities coupled by static stress transfer over hundreds of seismic cycles. We statistically analyze long (c. 500 ka) time series of M8‐9 analogue earthquakes simulated using a seismotectonic scale model approach with two aims: First, to constrain probabilistic measures (frequency‐size distribution, variability) useful for hazard assessment and, second, to relate them with geometric observables (coseismic slip pattern, locking pattern). We find that the number of synchronized asperity failures relative to the number of individual asperity failures as well as the coefficients of variation of recurrence intervals and seismic moment scale with the logarithm of stress coupling between the asperities. Accordingly, tighter packed asperities tend to recur more periodically and with a more characteristic magnitude while more distant asperities show clustering of more variable sized events. The probability of synchronized failures seems to be controlled to first order by geometrical relations (i.e., spacing and offset of asperities). The effects of rheological properties are evident but it remains to be explored to which extent they vary in nature and how sensitive the system is to those.


Introduction
Giant magnitude 9 earthquakes unzip up to 1,000-km-long segments of active plate margins. Such long ruptures include failure of several asperities. Prerequisites to fail synchronously (or sequentially in short succession, i.e., within seconds) are a homogeneous high stress level along the margin (i.e., a late interseismic stage in different segments along the megathrust) and a trigger for nucleation which might be very small depending on the state of synchronization of neighboring fault segments (e.g., Scholz, 2010). Ruff (1996) introduced the idea of synchronization of the seismic cycle "clocks" in subduction zones by static stress transfer (stress coupling) leading to giant earthquakes. He developed and analyzed a simple mechanical model consisting of two frictional spring-sliders coupled by a spring as an analog of a segmented subduction zone with segments interacting by means of stress coupling (Figure 1a). He hypothesized that while individual recurrence intervals may initially be different (controlled by the individual frictional strength and spring stiffness) stress coupling may introduce variability and cause synchronization over multiple seismic cycles as mutual coseismic loading causes "clock advance" of the neighboring segments (Figure 1c). In a modern view Ruff's (1996) idea is based on static Coulomb failure stress (CFS) transfer between asperities on a fault (megathrust) bounded by an elastic medium (upper plate wedge, Figure 1b). In this model, clock advances shorten the seismic cycle of neighboring asperities gradually and trigger, if sufficiently late in the seismic cycle, synchronized events ( Figure 1c). Larger clock advances (higher steps in such a figure) should potentially trigger more synchronous events as the sensitive fraction of the seismic cycle is proportionally longer. The model thus predicts that synchronization scales with stress coupling.
Here we test this model and the effects of asperity distribution on stress coupling, earthquake recurrence behavior, and synchronization. The first to model such a system realistically were Kaneko et al. (2010). They came up with a fully dynamic simulation of a pair of velocity-weakening asperities separated by a velocity-strengthening barrier. This simulation demonstrated the role of the asperity distribution and frictional properties of the barrier in controlling rupture propagation across it. Because of the computational costs of such fully dynamic numerical simulations, the lengths of the simulated earthquakes where rather limited to few tens of cycles and most of the model runs were performed in 2-D.
Here we realize similar models by means of seismotectonic scale modeling  which allows a realistic simulation of comparatively long analogue earthquake sequences with up to 500 events at a rather low experiment and time cost compared to numerical simulation. We simulate a subduction zone megathrust system in an archetypical setup with two square seismogenic patches (asperities) characterized by velocity-weakening and therefore unstable stick-slip frictional behavior. The asperities are surrounded by velocity-strengthening material displaying stable creep and acting as a barrier to seismic slip. Stress coupling by means of static Coulomb failure stress (CFS) transfer is realistically implemented by the elastic wedge overlying the megathrust and quantified using elastic dislocation modeling. While frictional and elastic properties are kept constant we change the asperity distribution by varying the relative position of the two seismogenic patches along-strike and across-strike allowing us to explore the effects of variable stress coupling and strength contrasts between the two asperities.
Our study complements and extends recent analogue models by Corbi et al. (2017) who tested the geometric aspects of Kaneko et al. (2010)'s simulation using a seismotectonic scale model similar to the one we use. They were able to verify experimentally the major role of the geometric relation between the asperities in synchronization. While they were able to reproduce both the numerical results by Kaneko et al. (2010) as well as natural observations from Japan, the significance of frictional properties remained unexplored by Corbi et al. (2017).
Here we complement these studies first by providing an analogue scale model with a different set of frictional properties compared to Corbi et al. (2017) to allow testing their significance more specifically. Figure 1. The concept of stress coupling and synchronization in subduction zones: (a) A system of coupled spring sliders as depicted by Ruff (1996); (b) the corresponding system of asperities coupled by static Coulomb failure stress transfer in an elastic wedge; (c) time series of stress (or strain) history (cyclic loading) of the two sliders/asperities. Note that the model predicts shortening of cycles due to clock advance early in the cycle and triggering of synchronized events by clock advance sufficiently late in the seismic cycle.
Second, we introduce a strength contrast between the two asperities, a factor which has not been tested experimentally or numerically so far. Third, we generated about 10 times longer analogue sequences (up to 0.5 Million years long including several hundreds of M8-M9 events) allowing a more rigorous statistical analysis and more reliable tests for statistical significance. All data underlying this study are published open access in Rosenau et al. (2019). 2. Modeling and Analysis Methods 2.1. Seismotectonic Scale Modeling of a Subduction Megathrust Setting 2.1.1. Experimental Setup and Scaling Seismotectonic scale modeling is a cost-effective method to simulate long earthquake sequences in a fully three-dimensional, dynamic, and spatiotemporally quasi-continuous framework (e.g., Caniven et al., 2015Caniven et al., , 2017Corbi et al., 2013Corbi et al., , 2017Rosenau et al., , 2017. Here we recall the basics of the approach and report modifications specific to the present study. The experimental setup used in this study is a development from an earlier quasi-two-dimensional setup used for seismotectonic scale modeling by Rosenau et al. ( , 2010 where the method has been explained in detail. The setup used in the current study is six times wider with respect to previous experiments and therefore truly 3-D and allows simulating along-strike analogue earthquake ruptures (e.g., Pipping et al., 2016). The experimental device consists of a glass-sided box (100 cm across-strike, 60 cm along-strike, and 50 cm deep) with a 15°dipping basal conveyer plate equipped with blades to ensure the transfer of basal stresses into the overlying material. On top of the conveyer plate an elastoplastic wedge (subduction forearc model) is set up at appropriate scale and compressed against a rigid and fixed backwall (Figure 2a).
Dynamic similarity of the laboratory scale model with the natural prototype requires the ratios of forces, which are expressed as dimensionless numbers, to be the same as in nature (Hubbert, 1937). We use the following set of dimensionless numbers to ensure similarity with respect to strength τ, gravity G, and inertia I: 1. The ratio Γ between gravitation and strength (either elastic, frictional, or viscous) is (1) where ρ is the material density, l is a characteristic length, g is the gravitational acceleration, and τ is the elastic, frictional, or viscous strength.

The Froude Number Fr relates gravitation and inertia and is
where v′ is a characteristic coseismic velocity (e.g., rupture velocity).

The Cauchy Number Ca relates inertia and elasticity and is
where E is Young's modulus.
By keeping these dimensionless numbers the same in an experiment executed in the Earth's gravity field as in nature, the following scaling relationships are derived from equations (1) to (3): where "*" marks the model numbers and values. The ratios between model and natural prototype values are known as the scaling factors (Hubbert, 1937).
These scaling relationships dictate the experimental conditions and material properties (Table 1) for a given length scale and material density. The model materials (i.e., sugar, rubber, rice, silicone) used here are about three times less dense than their rock equivalent and designed at a length scale (l*/l) = 3.3·10 −6 such that 1 cm in the scale model corresponds to 3 km in nature. According to equations (4)-(6) it follows that the scale model has to be weaker than the natural prototype by a factor (τ*/τ) = 1.1·10 −6 and should deform 500 times slower during analogue earthquakes in order to properly scale the body forces. The corresponding coseismic time scale is (t*/t) = 1.8·10 −3 (i.e., 0.1 s in the lab corresponds to about a minute in nature). Because this dynamic time scale would result in unsuitable long recurrence intervals of analogue earthquakes in the laboratory and because inertial forces can be neglected during the quasistatic interevent period we scale the interseismic periods based on the Ramberg Number (Ramberg, 1967) which is derived from equation (1) for the viscous regime where the strength is strain rate controlled: where η is viscosity and v a characteristic interseismic velocity. According to the material properties (i.e., viscosity, density) and convergence velocity used in our experiments an interseismic time scaling factor of 1.3·10 −10 (1 s in the lab scales to~250 years) is derived (Table 1). of megathrust setup with calculated CFS changes ΔCFS on receiver patch (normalized to stress drop Δτ on trigger patch). Note the logarithmic fall-off with distance from the trigger asperity (see also Figure S2). Da and Db refer to the parameters used by Corbi et al. (2017). (c) Parameter space of this study: Asperity spacing (dx) and offset (dy) and corresponding stress coupling log[ΔCFS/Δτ] in color and isolines. For this plot slip is imposed on asperity #1 (trigger patch) and mean stress change calculated on asperity #2 (receiver patch). Gray shaded area corresponds to the subspace realized experimentally. Size of the asperities has not been changed in this study. CFS = Coulomb failure stress.
Note that scale models embody strong simplifications with respect to the natural prototype. Therefore, their application may be limited to those aspects realistically implemented (i.e., frictional behavior, elastic loading, etc.). See Rosenau et al. (2017) for a review of the seismotectonic scale modeling approach.

Scale Model Configuration and Material Properties
The generalized subduction zone model presented here is analogous to a 300-km-wide and 180-km-long forearc section from the trench to the volcanic arc ( Figure 2a). The scale model is made up of a granular wedge of an elastic-frictional plastic (elastoplastic) mixture of EPDM (Ethylene-propylene-diene-monomer) rubber pellets with refined sugar and flavored rice representing the brittle forearc lithosphere. The wedge overlies PDMS (Polydimethyl-siloxane) silicone oil representing the viscoelastic asthenosphere. We generalize the natural subduction geometry by considering a planar, 15°-dipping megathrust between an upper plate made up of~60-km-thick lithosphere and~20-km-thick asthenosphere below the arc and an oceanic plate. The latter is represented by a conveyer plate pulled constantly via a spring-loaded thrust pad at 50 μm/s simulating plate convergence at a long-term rate of about 60 mm/a in nature.
The model megathrust is defined by a few millimeters wide granular shear zone which forms at the base of the wedge (analog to a "subduction channel," Shreve & Cloos, 1986). It is characterized by rate-and statedependent frictional behavior similar to nature (Scholz, 1998). In particular, it includes two square, velocity-weakening patches of basal area A (20 × 20 cm~60 × 60 km) displaying stick-slip deformation and mimicking a pair of seismogenic asperities separated by a velocity-strengthening, aseismic barrier. The seismogenic depth range is limited between~5 and 40 km as constrained by geothermal models (e.g., Oleskevich et al., 1999) and geodetic inversions (e.g., Li et al., 2015) of natural prototypes. The two stick-slip patches (asperities) are actually thin cuboid volumes of flavored rice embedded into the wedge material during the manual sieving procedure. During an experiment the cuboids are sheared into trapezoids without significant geometric change of the actual shear plane. The friction rate-parameter a-b within the asperities is approximately −0.02 ). The wedge material including the barrier separating the two asperities as well as updip and downdip regions of the asperities are characterized by "aseismic" slip or stable sliding (creep) controlled by the velocity-strengthening behavior (a-b~+0.02, ) of frictional slip in the sugar-EPDM mixture. A viscous wedge made of PDMS (Polydimethylsiloxane) silicone oil (Korasilon G30M) corresponds to the mantle wedge in nature. Material properties of this seismotectonic  Rosenau et al. ( , 2010 and Rudolf et al. (2016) and are reported in Table 1.
The two asperities have an along subduction zone strike center-to-center distance (hereafter called spacing) dx (25, 30, 35, 40 cm~75, 90, 105, 120 km in nature) and a relative shift across subduction zone strike (hereafter called offset) dy (0, 10, 20 cm~0, 30, 60 km in nature, Figure 2b). This allows exploring the effects of asperity distribution on stress coupling (as defined below in section 2.2.2) as well as strength contrast. We define the latter as the shear strength of the weaker (shallower) patch (asperity #2) relative to the stronger (deeper) patch (asperity #1): Strength contrast therefore ranges theoretically from close to 0 to 1. Note the somewhat counterintuitive effect that low strength contrasts are reflected by high τ 2 /τ 1 values. In total 12 configurations of asperity distributions have been realized in which we vary the strength contrast from 0.6 to 1.0 and the stress coupling from a few ppm to few percent ( Figure 2c). The experimental runs took place under normal gravity conditions and in a dry room climate (22-23°C, 30-40% humidity).

Experimental Monitoring and Strain Analysis
For strain analysis of the evolving model wedges we use an optical image acquisition and correlation system (particle image velocimetry, PIV StrainMaster by LaVision, Germany, see Adam et al., 2005;Rosenau et al., , 2010Rosenau et al., , 2017, for applications in analogue tectonic and earthquake modeling).
During an experiment, the locations of particles on the model surface (i.e., within the x-y plane of the model, Figure 2) are recorded by sequential 11 MPx digital images of a 14-bit monochrome chargecoupled device camera acquired at a frequency of 10 Hz. The x-y displacement vector field between successive images is then determined by cross-correlation of textural differences (i.e., gray values) formed by groups of particles using a Fast Fourier Transform algorithm. The spatial resolution of the final displacement vector grid is~3 mm or about 1 km in nature. For each grid-cell, an average horizontal displacement vector is determined at micrometer precision (~decimeter scale in nature). This allows for observing episodic surface deformation events corresponding to earthquakes of moment magnitude M w > 8. Analogue earthquakes are characterized by episodic, usually more than one order-of-magnitude increased strain rates and a change in polarity of the wedge deformation from "landward" motion (in negative y-direction) and shortening of the wedge during the interseismic stage to "trenchward" motion and extension of the wedge during the coseismic stage ( Figure 3). Earthquakes typically occur within a 0.1-s time interval, that is, are captured by a single image.

Elastic Dislocation Modeling
We use elastic dislocation modeling following Okada (1985Okada ( , 1992 for coseismic slip inversion and static Coulomb failure stress (CFS) transfer calculation. We employ the Matlab-based software package "Coulomb" by Toda et al. (2011). The model setup for elastic modeling uses the up-scaled geometry of the analogue model and according prototype material properties (see Figure 2b).

Slip Inversion
Surface deformation during analogue earthquakes as captured by PIV is converted into coseismic slip along the megathrust using inversion factors derived by forward elastic dislocation modeling. Accordingly, we find the factors relating horizontal surface deformation in the trenchward (y-) direction above the dislocation center to dislocation slip D at depth to range between 0.2 and 0.5 depending nonlinearly on the trench distance (or depth) of the dislocation ( Figure S1). Shallow dislocations show larger factors, that is, are less attenuated. We do not aim at a formal inversion or distributed slip modeling. Instead, we consider here mean coseismic surface displacement over the projected surface area of the seismogenic patch to be a valuable proxy for mean coseismic slip over the asperity at depth.

Stress Coupling
For quantifying the interaction by means of stress coupling between the asperities we follow the principles of static CFS transfer modeling as established by King et al. (1994), Toda and Stein (2002), and Lin and Stein (2004).
The model setup for CFS modeling is such that we impose uniform thrust slip on a trigger patch and average the predicted CFS change (ΔCFS) for thrust faulting on the receiver patch (e.g., Figures 2b, 2c, and S2). To

10.1029/2018JB016597
Journal of Geophysical Research: Solid Earth . Note the evolution from an initial "aseismic stage" toward a "seismic" steady state with a transient stage in between. We used only data from the seismic steady state for further analysis. Note also the asymmetry in displacements above shallow and deep asperity which is related to the free-surface effect. (d and e) Zoom into incremental surface displacement time series of the two event clusters A and B, respectively.

10.1029/2018JB016597
Journal of Geophysical Research: Solid Earth derive at a dimensionless parameter describing the asperities' interaction ("stress coupling") we normalize ΔCFS on the receiver patch to static stress drop on the trigger patch.
Stress coupling ¼ ΔCFS=Δτ: We calculate the static stress drop Δτ as proportional to the product of strain (i.e., mean slip D over patch dimension A 1/2 ) and Young's modulus E (e.g. Lay & Wallace, 1995): Because of the uncertainties in natural material properties and actual rupture geometry, we use a Young's modulus E of 100 GPa and a shape factor S, which is close to one in nature, of 1 for calculation. In the presence of a strength contrast we do this calculation for both cases where either the deep or shallow asperity act as a trigger patch and average the resulting stress changes on the receiver patches. Note that stress coupling as defined here is a parameter independent of the actual slip on the trigger patch.
In the present setup stress coupling is in the order of less than a ppm up to 1% similar to nature (Figures 2b,  2c, and S2). Stress coupling falls off exponentially with distance and varies nonlinearly across-strike of the megathrust as a function of asperity spacing (dx) and offset (dy).

Surface Deformation Time Series Analysis: Event Detection and Classification
Experimental time series of surface deformation consist of typically a sequence of 30,000 images and corresponding incremental displacement vector fields for each experiment. To detect analogue earthquakes from such a big data set we usually rely on computational algorithms sensitive to accelerations validated by visual inspection. However, because of experimental noise such a kinematic approach based on thresholding velocity usually has a high detection limit. Instead of thresholding velocities to detect earthquakes we here employ a numerical time series analysis technique developed in computational statistics. This allows us to automatically detect and classify events which can be below the detection threshold of classical kinematic approaches.
As input we use the surface deformation time series at the orginal frequency of 10 Hz. In particular we use across-strike (y-) displacements averaged over the surface projected area of the two individual asperities ( Figure 3). Those data typically show an initial aseismic phase without much activity reflecting stress buildup and reorganization within the scale model ( Figure 3c) followed by a transient phase of increasing activity leading over into a "seismic" steady state. We use observations from this steady state for further analysis.
To analyze the obtained experimental time series, we deploy a nonparametric time series analysis methodology called Finite-Element-Method with Bounded Variation of model parameters (FEM-BV; Horenko, 2009Horenko, , 2010Metzner et al., 2012). Although it is computationally more expensive than the common methods, FEM-BV has several important conceptual advantages that were recently illustrated for various time series analysis applications in geosciences Kaiser et al., 2015;O'Kane et al., 2016;Risbey et al., 2015;Vercauteren & Klein, 2015). This nonparametric method is automatized, does not rely on any tunable user-defined parameters (like thresholds values for the event identification) and allows to go beyond strong parametric assumptions (like linearity, Gauss or Poisson distribution assumptions for observed densities, stationarity, or Markovianity)-assumptions that are a constitutive part of the more common statistical time series analysis approaches like multilinear regression, Hidden Markov Models, or clustering methods (e.g., Shearer & Stark, 2012). Going beyond these assumptions is especially important since analyzed data exhibits a strong regime-transition behavior, is nonstationary, non-Markovian, and non-Gaussian in the regimes. Moreover, defining ad hoc threshold values for the events could potentially introduce a user-defined bias. We refer to Metzner et al. (2012) for mathematical/statistical details of the FEM-BV methodology-as well as for its computational comparison with more common time series analysis methodologies.
The result of this analysis is a catalogue of analogue earthquake events classified into "solo" events which occur either in asperity #1 (deep/left with respect to Figure 2b) or in asperity #2 (right/shallow) versus "synchronized" (double) events where both asperities fail synchronously within one time frame (0.1 s; Figure 4a and data publication Rosenau et al., 2019). A further classification of solo events is based on interevent time. We refer to two events which occur separately in the two asperities within successive frames (i.e., T rec = 0.1 s)

10.1029/2018JB016597
Journal of Geophysical Research: Solid Earth as "clustered" events, with the first being the "main event" and the second event being an "aftershock." All events can also be divided based on their faulting mechanism (i.e., thrust versus normal faulting), however, it appears that normal faulting is only observed in aftershocks. A flow-chart describing the event classification is shown in Figure S3.

Statistical Analysis of Analogue Earthquake Sequences
Here we are ultimately interested in the control of extrinsic parameters (stress coupling and strength contrast) on the variability of earthquake sequences in general, and the probability of synchronized events in special. A simple measure of probability of any class of events used by earlier studies (e.g., Corbi et al., 2017;Kaneko et al., 2010) is the number of events of that class relative to the total number of events. Based on the long sequences of analogue earthquakes in our synthetic catalogs we can go, however, beyond this simple approach and explore the intrinsic (stochastic) and extrinsic (parameter-controlled) variability of the frequency-size distributions by means of univariate and bivariate statistics.
In particular, we use probability density functions (PDFs) and their coefficients of variation (CV) to describe the intrinsic variability of individual earthquake sequences. PDFs of moment magnitudes M w (Figures 4b  and S4) mimic the high-magnitude (M8+) end of natural frequency-size relationships while PDFs of the seismic moment M 0 and the recurrence interval T rec (Figures 4c and S4) and their CV are used to differentiate between periodic and aperiodic (e.g., random, clustered) occurrence of events.
The CV is defined as the standard deviation σ of a data population normalized to their mean value M: and serve as a first-order proxy for the intrinsic variability of the frequency-size distributions in individual earthquake sequences: A CV of 1 characterizes a random distribution of recurrence intervals or  Figure S4 for PDFs of all experiments. PDF = probability density function.

10.1029/2018JB016597
Journal of Geophysical Research: Solid Earth magnitudes. A CV <1 suggests a (quasi-) periodic recurrence or a characteristic magnitude of events. A CV >1 is characteristic of clustering in time or a dominance in small magnitudes (e.g., Kuehn et al., 2008;).
To test the control of extrinsic parameters (stress coupling and strength contrast) on recurrence behavior we finally analyzed the correlation of mean values and of CVs of M 0 and T rec with (the logarithm of) stress coupling and with strength contrast by means of linear regression. We use probability (p) values and coefficients of determination (R 2 ) to describe the goodness-of-fit of, and the extrinsic variability explained by, (log)linear regression models, respectively. A p value threshold of 0.05 is used to reject or accept the (log)linear model. P values between 0.05 and 0.1 are considered marginal.

Seismic Performance of the Scale Model
All experiments followed a consistent characteristic evolution from aseismic to seismic behavior ( Figure 3c): An early initial phase (about minutes 0-15) of aseismic deformation is dominated by distributed wedge shortening and minor thickening and followed by a transient phase of increasing seismicity (about minutes [15][16][17][18][19][20]. After about 20 min the wedge is in a steady state characterized by stationary time series of cyclic surface accelerations reflecting analogue earthquakes. The seismic evolution reflects the progression of the model deformation from dominantly plastic to dominantly elastic . A typical earthquake catalogue simulated by our scale model in the steady state consist of up to 500 events of moment magnitude M w 8-9 which occur over a time-period of about 500 ka in nature (Figure 4a). M8 events usually involve only one asperity while a synchronized failure of both asperities typically results in M9 events. Analogue earthquakes are always followed by afterslip lasting for not more than one frame (0.1 s or~25 years in nature) surrounding the asperities (Figures 3a and 3b). Generally, the shallow asperity generates more surface displacement than the deep one: this is related to static effects as predicted by elastic dislocation modeling ( Figure S1). The picture inverts when the correction for depth of dislocation is applied. In that case, deeper asperities show larger slip. This is consistent with higher loads causing higher frictional strength at greater depth as predicted by Mohr-Coulomb theory. As a consequence, the deeper asperities are mechanically stronger, fail at higher shear stresses and are therefore able to accumulate more slip deficit in the correspondingly longer interseismic periods compared to the shallow asperities.
A minority of aftershocks are relatively small normal faulting events. We interpret those as a result of dynamic overshoot during the preceding thrust event. Normal events occur almost exclusively in the shallow asperity. We include those rare normal events in our analysis since they represent an integral part of the long-term slip budget. Accordingly, they show up with a negative seismic moment in Figure 4a.

Event Probability
When analyzing the probability of occurrence of the different classes of events as a function of stress coupling ΔCFS/Δτ and strength contrast τ2/τ1 a clear picture emerges (Figure 5 and Table S1). Accordingly, synchronized events increase in number from 20% to 80% of total events as stress coupling increases by two orders of magnitude (from less than a ppm up to a percent, Figure 5a). Consistently, solo events decrease (Figure 5b). This testifies a higher degree of synchronization in strongly coupled systems. Regarding clustered events, only pairs with thrust aftershocks show a clear correlation with stress coupling (Figure 5c). Those clustered pairs with normal faulting aftershocks show no clear correlation with stress coupling but a negative correlation with τ2/τ1 (Figure 5d). This is consistent with overshoots occurring preferentially in shallow regions of the wedge. All other classes of events show, however, no significant correlation with strength contrast. An apparent increase of the range of proportions of individual events with τ2/τ1 reflects the systematically wider range in stress coupling realized for models with high τ2/τ1 (e.g., Figure 6).

Frequency-Size Distributions
Frequency-size distributions of simulated earthquakes share similar shapes. The probability density functions (PDFs) of moment magnitude are generally skewed negatively (toward the left) and very peaked as exemplified in Figure 4b. The PDFs of recurrence intervals are generally bimodal, characterized by a peak at short periods (0.1 s or 25 years) and a quasi-normally distributed peak around the mean recurrence interval of large events (at~1000-1500 years) as exemplified in Figure 4c.

Intrinsic and Extrinsic Variability
Plotting mean recurrence intervals and mean seismic moments of all events in individual models and their variability in terms of CV into the parameter space ( Figure 6) shows the following: Mean recurrence interval and seismic moment both increase with an increase in stress coupling (Figures 6a and 6b). At the same time their CVs decrease (Figures 6c  and 6d).
We interpret this correlation of M 0 and T rec with stress coupling as reflecting a dynamic interaction causing higher slip in case of more strongly coupled asperities. Larger slip per event consistently lengthens the interseismic period resulting in longer recurrence intervals in order to keep the long-term seismic slip rate constant. The increase in magnitude seems also to regulate the seismic cycle, decreasing the CV of M 0 and T rec by about 20%.
A weak positive correlation seems to exist between T rec and strength contrast but the p value is marginal (p = 0.1, Table S1) and only about one fourth of the variability is explained by extrinsic control (R 2 = 0.25, Table S1). Accordingly, earthquake frequency increases as the weak asperity becomes weaker. We interpret this as being the behavior predicted by Ruff (1996) where the weaker asperity, which has intrinsically shorter recurrence times, causes clock advance of the stronger asperity, which has intrinsically longer recurrence intervals. No correlation exists between M 0 as well as the associated CVs with strength contrast.
The correlations of M 0 and T rec with ΔCFS/Δτ are replotted in Figure 7 with a differentiation between all events and synchronized events only to explore the effects of stress coupling on intrinsic and extrinsic variability in more detail. Consistently, consideration of only synchronized events, which are generally larger than solo events, increases mean seismic moment and mean recurrence interval and decreases the associated CVs.
More interestingly, however, is the observation that the trends differ for the two groups of events: For example, the positive correlation of T rec with stress coupling observed for all events is inverted to a negative correlation if only synchronized events are considered (Figure 7a). This is the result of synchronized events being systematically rarer in more weakly coupled systems as has been predicted by Ruff (1996). At the same time, recurrence intervals of synchronized events are more sensitive to stress coupling having a higher R 2 = 0.44 than the recurrence intervals of all events (R 2 = 0.29). In fact, synchronized events recur almost Figure 6. Visualization of the variation of mean recurrence interval and its CV (a and c) and mean seismic moment and its CV (b and d) in the parameters space, that is, with stress coupling (ΔCFS/Δτ) and strength contrast (τ2/τ1). See Table S2 for linear regression analysis results. CFS = Coulomb failure stress; CV = coefficient of variation. randomly for weakly coupled systems (CVT rec~0 .9) and periodically for strongly coupled systems (CVT rec~0 .3). On the other side, the magnitudes' variability is minimized for synchronized events as compared to all events (CV 0.2 vs 0.4-0.6) and independent of stress coupling (p = 0.5).
In conclusion, our analysis indicates that stress coupling has the effect of causing more periodic, characteristic events. Compared to a reference experiment with a single asperity by Rosenau et al. (2017), however, any system of two asperities realized here shows a more complex recurrence behavior indicated by generally more frequent and smaller events (lower T rec and M 0 ) and significantly higher CVs (Figure 7).

Relation Between Asperity Distribution and Recurrence Behavior: A Characteristic Length Scale in Nature?
Based on experimentally simulated long subduction earthquake records we are able to shed light on the control of asperity distribution on the variability of subduction earthquakes in terms of magnitude and recurrence intervals. Rosenau et al. (2017) showed that the transition from a single asperity to two weakly coupled asperities involves a principle change from periodic (recurrence interval's CV~0.2) toward more randomly occurring earthquakes (CV~1). This is consistent with spring-slider models suggesting a single Figure 7. Bivariate plots of the most significant trends from Figure 6 with two classes of events differentiated (all events = black dots, synchronized events = red dots): Variation of mean recurrence interval T rec (a) and its CV (c) as well as mean seismic moment M 0 (b) and its CV (d) with stress coupling (ΔCFS/Δτ). Linear regressions with p < 0.05 are considered significant correlations, those with p = 0.05-0.1 marginal. See Table S2 for all p and R 2 values and text for interpretation of the correlations. Horizontal stippled line indicates the corresponding value in an experiment with a single asperity  for reference. CFS = Coulomb failure stress; CV = coefficient of variation.
isolated spring-slider system to be periodic while a coupled pair of springsliders shows a complex behavior (e.g., Ruff, 1996). The system simulated here shows a strong correlation between the stress coupling (controlled by asperity distribution) and recurrence variability increasing from 0.3 to 1 as stress coupling decreases (Figure 7c). At the same time earthquakes become more variable in magnitude as stress coupling decreases suggesting a causal link between characteristic earthquake behavior and asperity distribution.
The range in CVs in recurrence intervals spans a considerable larger range than what is seen in natural examples which is usually characterized by a CV <0.4. For example, the Holocene tsunami record offshore western North America suggest that great M w 9 Cascadia subduction zone earthquakes have occurred about every 500 to 600 years during the past 10 kyr (Goldfinger et al., 2003) with a CV of 0.36-0.39 (Sykes & Menke, 2006). For the Nankai trough, Sykes and Menke (2006) report a CV of 0.26-0.27. In the Northern Chile-Southern Peru seismic gap which last broke in 1877 (M w 8.8) the reported historical recurrence interval for the past 500 years has been estimated at 111 ± 33 years (Comte & Pardo, 1991) resulting in a CV of 0.3. Similarly, in southern Chile, in the area of the great 1960 and 2010 earthquakes, leveling and dating of Holocene strandlines by Bookhagen et al. (2006) suggests that great earthquakes have occurred every 180 ± 65 years over the last 3 to 4 kyr, from which a CV = 0.36 can be calculated. Clearly the regularity of great megathrust earthquakes implied by their low CV breaks down below about M w 8. However, b-values significantly below 1 for M w 5-8 earthquakes in many subduction zone (e.g., Bilek & Lay, 2018) suggest a nonrandom behavior also at smaller magnitude.
Although the natural data base is limited, this rather narrow range of low CVs and b-values in nature in combination with the here suggested causal link between regularity and asperity distribution let us speculate that there might be a characteristic length scale in the asperity distribution in nature. In our models a CV <0.4 in recurrence intervals is reached only by the narrow configurations where barriers between asperities are significantly smaller than the asperities themselves. Such a narrow asperity configuration can be found for example in the region of the 1960 and 2010 Chile earthquakes (Moreno et al., 2009(Moreno et al., , 2010(Moreno et al., , 2011. More examples can be found, for example, in Hayes (2017) finite fault model data base, however, a rigorous review of natural examples with respect to this aspect is beyond the scope of this paper.

Predicting Asperity Interaction: Toward Proxies for Barrier Efficiency
We simulated long time series of analog subduction megathrust earthquakes in order to constrain the recurrence pattern of a simple system with two asperities coupled by static stress transfer. Similar experiments  and numerical simulations (Kaneko et al., 2010) have been carried out to find the critical parameters controlling the probability of a rupture bridging the barrier and causing a synchronized failure of the asperities. We here add experimental data representing a different set of material parameters and geometries which allows testing the existing concepts and to identify the minimum set of parameters needed. Kaneko et al. (2010) suggested a set of parameters combined in a proxy for barrier efficiency called B. B is the ratio of the stress increase required to bridge the barrier to the coseismic stress drop. B included parameters which are directly and indirectly (involving assumptions) observable in nature (geometric, kinematic, dynamic, and friction parameters). Given the complexity of B and the uncertainty in the choice of some of the parameters included (e.g., frictional parameters), Corbi et al. (2017) aimed at a more simple proxy based solely on first-order geometric relationships easy to observe in nature, that is, the barrier-to-asperity length ratio Db/Da. With respect to these two proxies, we consider the stress coupling as defined here as a proxy for barrier efficiency of intermediate complexity. Similar to Db/Da, stress coupling can be inferred primarily from geometric observations (size and location of asperities).

10.1029/2018JB016597
Journal of Geophysical Research: Solid Earth In Figure 8 we compare the three proxies based on the setup presented in this study. Obviously, there is a good correlation between stress coupling, B, and Db/Da. Db/Da seems slightly more sensitive to stress coupling than B as suggested by its steeper slope in this plot. In any case, a correlation coefficient (R 2 ) of 0.6 to 0.8 suggests general interoperability of the three proxies. Figure 9 shows the collapse of existing experimental and numerical as well as real world data from the Nankai subduction zone  in a plot of percentage of synchronized events versus B. In contrast, plotting those data against Db/Da separates the data into roughly parallel trends. Because the data used represent a wide spectrum of geometrical and frictional parameters, the collapse indicates the universal nature of the proxy B for anticipating the probability (i.e. percentage) of synchronized events.
On the other hand, the systematic offset trends suggest that while Db/Da seems to have a strong control on synchronization, material properties cannot be neglected. For instance, it appears that the setup used in the present study generates synchronized events more easily. While for the experiments by Corbi et al. (2017) and the natural example a threshold for synchronized events at Db/Da of 0.5 emerges, in the experiments presented here this threshold is significantly higher (>1). This suggests that the barrier in the Corbi et al. (2017) experiments as well as in the Nankai area are mechanically more effective than in our setup.
We conclude that for the moment, the full complexity of the proxy B by Kaneko et al. (2010) is needed to account for the variability of mechanical parameters present in the experiments. To which extent these parameters vary in nature and therefore control the threshold value of Db/Da remains to be explored.

Conclusions
Based on experiments with seismotectonic scale models generating long time series of analog subduction megathrust earthquakes we explored the process of interaction and synchronization of two velocityweakening asperities separated by a velocity-strengthening barrier. We found the following: 1. Synchronization of asperities is controlled by the static Coulomb failure stress (CFS) transfer, here quantified by the stress coupling ΔCFS/Δτ. Accordingly, the percentage of synchronized events scales with the logarithm of (normalized) CFS change. 2. A strength contrast between the two asperities has no significant effect on synchronization but decreases the recurrence interval of synchronized events because the weaker asperity dictates the system's recurrence intervals. 3. Analogue earthquakes in strongly coupled systems (narrower asperity distribution) recur more periodically and with a more characteristic magnitude than in weakly coupled systems. 4. A narrow asperity distribution might be typical for natural subduction zones characterized by quasiperiodic recurrence.
Three proxies for the barrier efficiency, B (Kaneko et al., 2010), Db/Da , and the newly defined stress coupling ΔCFS/Δτ have been cross-validated and tested for applicability: 1. Db/Da is the most simple and easiest to apply proxy and incorporates the most sensitive parameters to work first-order. It relies on geometries which-if they are stationary over multiple seismic cycles-we are able to constrain using interseismic locking and paleoseismological observations. 2. B is the most universal proxy and it captures the physics-but several parameters are needed and not well constrained or uncertain in nature.

ΔCFS/Δτ is of intermediate complexity and interoperable with Db/Da and B.
In order to arrive at a minimum set of parameters necessary to describe seismic hazard in subduction zones we suggest to further explore the variability of those parameters in B, which are not well known in nature, to define the sensitivity of simpler proxies and to aim at constraining their upper and lower bounds.