Stick-slip in a stack: how slip dissonance reveals aging

We perform physical and numerical experiments to study the stick-slip response of a stack of slabs in contact through dry frictional interfaces driven in quasistatic shear. The ratio between the drive's stiffness and the slab's shear stiffness controls the presence or absence of slip synchronization. A sufficiently high stiffness ratio leads to synchronization, comprising periodic slip events in which all interfaces slip simultaneously. A lower stiffness ratio leads to asynchronous slips and, experimentally, to the stick-slip amplitude becoming broadly distributed as the number of layers in the stack increases. We interpret this broadening in light of the combined effect of complex loading paths due to the asynchronous slips and creep. Consequently, the aging rate of the interfaces can be readily extracted from the stick-slip cycles, and it is found to be of the same order of magnitude as existing experimental results on a similar material. Finally, we discuss the emergence of slow slips and an increase in aging-rate variations when more slabs are added to the stack.


I. Introduction
Multiple frictional interfaces coupled by elasticity are ubiquitous in everyday objects including books [1,2], textiles [3][4][5], and multilayer composites [6,7].In geology, systems comprising multiple frictional interfaces are the norm rather than an exception.For example, layered rocks such as shale can show multiple sliding interfaces under shear [8,9].At the larger scales relevant for terrestrial faults, slips producing earthquakes are usually not isolated but embedded into complex fault networks [10].The mechanical response of such assemblies of frictional interfaces involves the coupling between the elastic deformation of the layers and the barriers to sliding of the interfaces.
Predicting the onset of slipping is a long-standing problem, even for a single frictional interface [11].Physical insight and understanding of this class of problems have been driven primarily by high-precision experiments of sliding PMMA blocks whose optical transparency enabled the spatiotemporal tracking of the local contact area [12].These pioneering experiments have elucidated that the onset of slip involves a rupture front that 'unzips' the interface.A correlation with strain measurements close to the interface showed that the stress field and dynamics of the front are well described by fracture mechanics with the fracture energy as the sole fitting parameter [13].However, the mechanism underlying the nucleation of the rupture front remains elusive, primarily due to experimental limitations, for which novel protocols are being proposed [14].
From a theoretical perspective, the most common models for the onset of frictional slippage [15][16][17][18] capture the phenomenology that sliding starts, and after a transient, continues in a steady state when the shear force F balances the friction forces µN , where N is the normal force and µ is the "friction coefficient".In these "rate-andstate" models, µ depends nonlinearly on the slip rate v and history θ: µ = µ(v, θ).At intermediate values of v, the friction coefficient is usually assumed to display slip weakening (µ is a decreasing function of v) such that the interface is unstable.During slip nucleation, the elasticity and inertia of the bulk have a stabilizing effect [19][20][21][22], such that there exists an effective flow curve whose steady-state displays a minimum µ c at v = v c [19].Consequently, any perturbation decays and vanishes if the applied load is F/N < µ c [23].At higher applied loads, the interface destabilizes if a slip patch reaches a critical size [21] beyond which its dynamics are well described by a sharp rupture front [19,21,24] that can be modeled by linear elastic fracture mechanics [13,25,26].A significantly debated question is the nature of the instability and consequently, how the critical size diverges as a function of F/N − µ c , e.g.[12,21,22,27].
A direct consequence of the phenomenology described above is that the interface can display stick-slip behavior when driven quasistatically (at a rate V ≪ v c ).The link between the stick-slip amplitude and the parameters of the rate-and-state models is debated [14,21,22,[28][29][30][31].One of the authors [22,29] recently proposed an encompassing theory linking the cited rate dependence and the nucleation of an instability by avalanches of microscopic failures.Beyond this athermal view, it is a wellknown experimental fact that the initial onset to sliding is history-dependent and increases with the time that the interfaces were at rest [17,[32][33][34].This aging behavior is associated with creep [35] and described by rate-andstate models where the variable θ introduced above is regarded as time.Creeping of the interfaces must affect the stick-slip amplitude, but disentangling its contribution is challenging because slip events occur at a narrowly distributed interval (on which statistical fluctuations of the stick-slip amplitude overwhelm the effect of aging [12]).
Beyond a single frictional interface, when multiple frictional interfaces are present, the elasticity of the bulk may potentially lead to a non-trivial coupling.For example, elastic interactions between faults may strongly affect their slip dynamics [36,37].In addition, acoustic waves transmitted through the elastic bulk may lead to remote triggering of earthquakes [38,39], though the large temporal separation suggests a complex coupling.Predicting the mechanical response of an assembly of elastic frictional interfaces is then a formidable but important challenge.In particular, identifying the key parameters coupling the layers together and elucidating the role of the number of interfaces are open questions.
Here, we report results from a combined experimental and numerical investigation on the quasistatic stick-slip response of a stack of elastic slabs in contact through frictional interfaces.Based on the ratio between the stiffness of the drive and the shear stiffness of the slab, we distinguish two regimes: "stiff" and "compliant" driving.In the stiff-driving regime, our numerical results exhibit periodic slips involving all the layers leading to narrowly distributed force drops.By contrast, in the compliantdriving regime, we observe both numerically and experimentally a decoupling of the slip events along the different layers, with interfaces sliding one by one.In the experiments, we find that this loss of periodicity is accompanied by a broadening of the distribution of the stickslip amplitudes with the number of layers.We highlight the role of interface creep in this broadening, exposed by the complex loading paths of the interfaces induced by mechanical coupling between the increasing number of layers.Overall, the stick-slip response of a stack to shear is controlled by the stiffness ratio between the drive and the layers whereby in the stiff-driving case the stack acts as one layer with periodic slips, while in the compliantdriving case a rich coupling between the layers makes slips much more unpredictable.

II. Definition of the problem
We assess the shear response of a model system comprising a stack of n ∈ [1,5] identical slabs of thickness h resting on a surface whose position is fixed.In Fig. 1, we present a schematic diagram of the system and a photograph of our experimental setup.Each slab, and its lowermost frictional interface, are numbered as i = 1, 2, . . ., n from below.We impose homogeneous shear between all the slabs by connecting each slab through identical springs of stiffness K to a lever that is driven to rotate around a fixed axis at a set rate (Fig. 1a).The spring connecting to the i-th layer is attached to the lever at a distance ih from the rotation axis.
For the drive, we impose the lever's top horizontal displacement U (t) = V t, where t represents time, and V is the imposed velocity, taken to be small enough for the interfaces to display stick-slip (such that V < v c [22,33] has to hold).Thus, our drive imposes a shear rate γ ≡ V /H, where H is the height of the lever, driving each spring i at a velocity v i = ih γ.During the periods in which the interfaces are 'stuck', this drive causes a monotonically increasing shear stress at each of the interfaces.
Our study seeks to address the following questions: (1) What are the relevant parameters controlling the slip synchronization of the interfaces?(2) How does the shear response evolve with an increasing number of layers n?

III. Stiff vs. Compliant Driving
For a system with a single frictional interface, the stickslip instability occurs only if the driving stiffness K < K c , where K c depends on the flow properties of the interface and the applied driving rate [15,16,33,40]; a summary of the calculation for a rate-and-state model is provided in Ref [33].The experimental and numerical systems are taken in the stick-slip regime (such that K < K c has to hold).We will argue below that with multiple interfaces the rigidity of the drive also controls the degree of slip synchronization.
To gain insight into the effect of the drive on the shear response of a stack, we regard the driving 'springs' as a parabolic potential energy imposing the mean position of each slab, such that the slab is free to build up shear.We now discuss what happens in the limit of stiff and compliant driving, defined next.
Rigidity ratio Φ.We define the rigidity ratio Φ ≡ K/K s , where K is the rigidity of the driving springs, and K s = AG/h is the shear rigidity of the slabs, with G the shear modulus, A the surface area of the frictional interface, and h the slab's height.This ratio Φ then quantifies the relative deformation of the driving springs in comparison to the shear deformation of the slabs.Below, we investigate and discuss two limit regimes: stiff (high-Φ) and compliant (low-Φ) driving.In the first, significant deformations occur within the layers (and the springs are stiff), whereas, in the second, the deformations occur in the springs (and the layers are stiff).The loading and slip behaviors of the stacks in these two limit regimes are schematically represented in Fig. 2. Stiff driving (high-Φ).Let us suppose that the interface a = 1 starts slipping.As the mean position of slab a is fixed by the stiff spring, the slab can only react with a negative shear deformation.Consequently, the shear stress on the interface above, i = a + 1, is increased.This can trigger a slip on that interface, which in turn causes a stress increase and slip on the interface above for the same reason, and so on for the rest of the stack.The cascade results in a multi-slip event that erases all memory of the system.In this case, the multi-layered system thus acts as a system with a single interface with effective properties [41], showing periodic stick-slip cycles.
Compliant driving (low-Φ).With soft springs, the system can respond to a slip at interface a by advancing the mean position of slabs i > a and compressing the driving springs.Therefore, the stress on the interfaces i > a relaxes, making a macroscopic multi-slip event unlikely.A sequence of single-slip events is thus to be expected.With an increasing number of layers, the slip sequence of the multiple interfaces may lose its periodicity.

IV. Numerical simulations A. Numerical model
We implement the model system shown schematically in (Fig. 1a) into numerical simulations.The numerical model consists of n + 1 identical elastic layers separated by frictional interfaces.Following [29], we idealize the frictional contact problem in order to focus on the dis-order in the shear response along the frictional interface.We assume that the interface is disordered but perfectly flat.Furthermore, we assume that contacts detach when a critical shear stress is reached.As such, the system is macroscopically at zero temperature (contacts do not detach spontaneously due to random fluctuations caused by an external temperature).Moreover, fluctuations of normal stress -e.g.due to an inhomogeneous load [12], inhomogeneous pore pressures [42], or partial slip [43] do not affect the contact strength.In particular, we consider a mesoscopic scale on which an effective 'block' of a finite width resists elastically to shear up to a threshold, after which it yields.The local slip then propagates until a new 'contact' is formed (i.e., it is again elastic but with a new threshold).In this framework, each block represents a frictional contact (or a patch of contacts that are so strongly coupled by elasticity that they act as an effective contact) that, upon yielding, forms a new contact with a new yielding threshold.
Geometrically we do not seek to precisely model Fig. 1a as its numerical treatment, together with the disorder, requires an intractably large number of blocks.Therefore, we consider periodic boundary conditions in the horizontal direction and control the mean position of each slab through a parabolic potential energy.
The details of the numerical model are as follows: each frictional interface consists of n x equal-sized square blocks of linear size l 0 that are completely linear elastic under volumetric deformation but yield under shear (deviatoric deformation) when a set yield stress is reached.Assuming that the yield threshold is isotropic in principal deviatoric strain space, this model now corresponds to a deviatoric potential energy that consists of a sequence of parabolic potentials in equivalent deviatoric strain space.The disorder arises from independently randomly drawing the yield strain sequence of each block.We assume that the blocks and the bulk have the same elastic moduli.
Differently from [29], we add a parabolic potential (with curvature K) to the mean horizontal position of each of the elastic slabs i > 0, adding a homogeneous force density per layer.The bottom layer is not driven through its mean position; instead, the position of the bottom edge is fixed.We set the mean horizontal position of a slab i equal to γih, such that γ represents the lever rotation (see Fig. 1a).
A key feature of the model is that shear can be applied according to the quasistatic protocol.Moreover, the linear elastic response permits an event-driven protocol.In alternation, we increase the shear by a finite amount ∆γ(t) that is maximized with the constraint that no microscopic yielding takes place (preserving mechanical equilibrium), and then add an infinitesimal shear δγ to trigger a microscopic event (after which we minimize energy).This allows us to perfectly separate events.
We choose n x = 2 × 3 6 , which is still tractable to simulate, but of the minimal order not to be dominated by finite size effects, as we checked for a single frictional in-terface [22,29].Furthermore, we take h/ℓ 0 ≈ n x /4 based on balancing h/(ℓ 0 n x ) small enough to have acoustic interactions while avoiding driving the blocks in a fixed displacement such that collective effects are suppressed if h ≪ n x ℓ 0 (e.g.[44]).
The above model predicts stick-slip behavior [29] when full inertial dynamics are considered (using overdamped dynamics, this model predicts the abundantly studied depinning transition [45]).We consider such inertial dynamics by applying the finite element method to discretize space.Along the frictional interface(s), elements coincide with the mesoscopic blocks.In the elastic slabs away from the frictional interface, the elements are coarsened to gain numerical efficiency (such that the height h is only approximated as we fix the aspect ratio of elements to one, see Appendix E).We use the velocity Verlet algorithm to integrate discrete time (with a time step significantly smaller than the time of a single oscillation in a well in one block).We remark that assuming periodicity requires us to add a small damping term to the inertial dynamics such that waves with a wavelength equal to the horizontal size of the system (n x l 0 ) are critically damped.Consequently, we must take h/ℓ 0 < n x to have acoustic coupling between the interfaces.
We note that following the underdamped dynamics in such a disordered system is very costly -it requires on the order of one billion time steps per realization.We are able to perform such simulations by using a dedicated optimized code, but mostly because we make two important assumptions.First, we coarsen elements in the elastic bulk as we discuss above.This is known to lead to spurious wave reflections [46], but we presume that these are not important because disorder already leads to a broad spectrum of wavelengths.Second, we assume small strains and rotations, such that the integration volume is assumed constant -reducing the computational cost of numerical integration.To still solve significant slips, the yield strains of each block are scaled by a small factor (all presented results have been renormalized to eliminate this factor).However, this still inherently limits the description to moderate slips.As such, we cannot choose an arbitrarily small driving stiffness K.

B. Numerical results
Our numerical model allows us to first illustrate the simple argumentation on the role of driving that we made above.We consider a driving rigidity such that Φ ≃ 10 −3 (stiff driving) and Φ ≃ 10 −6 (compliant driving).In the two-dimensional model, A = n x ℓ 0 , such that K s = 4G for our geometry; we use K = 10 −3 and K = 10 −6 and G = 1/2.In Fig. 3a and Fig. 3b, for stiff and compliant driving, respectively, we plot a typical macroscopic stress Σ (volume-averaged stress) as a function of applied shear γ.Note that the stress is shown in units of the typical yield stress of one block, and the rotation in units of the rotation needed to yield a typical block at i = 1.
Macroscopic slip events are defined when all blocks along one or more layers yield at least once.Below, we will refer to sliding interfaces and associated quantities by an index a, while i will be kept as the running index for the layers.Slip events correspond to macroscopic stress drops in Figs.3a and 3b, and we distinguish between 'single-slip' events (all blocks on a single layer yield at least once) and 'multi-slip' events (all blocks on more than one layer yield at least once).Stress drops produced by single-slip events are labeled following the color code introduced in Fig. 1, while multi-slip events are kept black.These slip events are separated by 'stick' intervals during which only microscopic events are observed, where one or several blocks yield at least once, as indicated with markers (black dots).
The results confirm that stiff driving causes a periodic stick-slip sequence with many slip events corresponding to multi-slip events (Fig. 3a) while compliant driving results in a seemingly less periodic sequence of single-slip events in Fig. 3b.This finding is supported by plotting the fraction of slip events involving s = 1, . . ., n interfaces in Fig. 3c for different n (see legend in Fig. 3d).On the one hand, stiff driving results in single-and multislip events for a comparable fraction of loading history (we discuss in Appendix F that sequences of single and multi-slip events alternate).On the other hand, compliant driving shows single slip in the large majority of slip events.
In the compliant regime, a direct measurement of the stress drop along the slipping interface a, ∆µ a , displays no n dependence in Fig. 3d.The quantity µ a is defined as the volume average stress on the blocks corresponding to weak layer a, also shown in units of the typical yield stress of one block.Given that, by construction, normal stress plays no role in our model, here µ a is akin to a friction coefficient.The finite width of the distribution is attributed to the inherent disorder of the interfaces.
As we have shown that, for high-Φ (stiff driving), multilayer stick-slip is apparently similar to that of a single interface, next, we concentrate on the low-Φ regime to explore a potential influence of the number of plates n.

V. Experiments
We proceed by proposing an experimental realization of the sheared-multilayer model system of Fig. 1a, adapted to measure the effect of the number of sliding layers n on the slip synchronicity and amplitude; see Fig. 1b.Similarly to the numerical model, the position of each slab is driven by connecting it to the driving lever through linear springs (see schematic in Fig. 1a).Naturally, connecting the spring to the edges of the slabs might introduce boundary effects, but this is mitigated by the fact that our experimental system is effectively much larger than our numerical model (given that it presumably has much more local contact patches).

A. Experimental apparatus
The experimental setup shown in Fig. 1b comprises a stack of frictional plates (color-coded from purple to orange), an actuating lever (green), and driving springs We indicate all (microscopic) yielding events with a black dot marker.Slip events on a single layer are indicated in color (see legend), while slip events in black involve more than one interface.(c) The fraction ρ(s) of macroscopic slip events involving s = 1, . . ., n layers, for stiff (dashed) and compliant (solid) driving; see legend in (d) for color-map and markers.(d) Distribution of stress drops at the slipping interface for different n in the compliant regime (for slip events on a single layer, for which s = 1 in (c)).See the main text for definitions and units.
(pink).Each component of the setup is detailed in Fig. 4 and below.
The stack is made of a set of rectangular PMMA slabs (Snow WH10 DC by Röhm), each of dimensions h = 10 mm, L = 150 mm and an out-of-plane width of 80 mm.A normal force N is applied on the topmost slab by a dead weight of 5 kg (N = 49 N).To ensure a spatially homogeneous contacting surface at this relatively low normal force (compared to other PMMA-PMMA friction experiments [13,14,47]), we use acrylic plates whose surface was pre-roughened with asperities of size ∼ 25 µm that are larger than potential natural height variations of PMMA (Fig. 4(d)).We assume that the normal force is uniformly distributed and that it is the same for each layer (the weight of each slab is less than 3% of that of the dead weight).
The stack is sheared by imposing the displacement at the top of the lever (H = 100 mm) at a constant speed V = 10 µm/s (i.e.γ = V /H = 10 −4 s −1 ), using a DC lin-ear actuator (L-220.70DG,Physiks Instruments) that is attached via a steel junction assumed rigid (Fig. 4(a.1)).The PMMA lever, made of two 6 mm thick and 100 mm wide parallel slabs, is sufficiently wide not to bend while pulling the slabs and rotates smoothly on ball bearings around its rotation axis.The total horizontal force F needed to rotate the lever (Fig. 1) is measured using a uniaxial force sensor (LRM200 25 lb, Futek) placed between the steel junction and the actuator (Fig. 4(a.2)).
The springs connecting the slabs to the lever are curved beams laser-cut from PMMA (colored in pink in Fig. 1b), with an equivalent stiffness of K = 55 N/mm when pulled or compressed along the horizontal axis.The beams have a rectangular 5 × 5 mm cross-section and are precurved with a radius of curvature of 100 mm over a cord of 150 mm in the middle and extended on both sides by two 25 mm straight portions.
The ends of the springs are attached to both the slabs and lever via ball-bearing links to ensure a free rotation and, thus, horizontal driving forces.The links are inserted between the two slabs of the lever and are made of a central steel pin press-fitted at the end of the spring (Fig. 4(b.1), orange).This pin is embedded into small ball bearings on each side (green), allowing a free relative rotation between the springs and the lever.The bearings are then press-fitted into acrylic components (pink) with two threaded holes.These components are screwed directly on both sides of the lever.The rotation axis of the pins is aligned precisely along a line passing through the center of rotation of the lever.Figs.4(b.1) and 4(b.2) show the side and top-side views of the link; and Fig. 4(b.3)shows the link prior to attachment to the lever.The other sides of the springs are linked to the frictional slabs.The spring-slab links have a similar design, but here the acrylic components are press-fit into the slabs.The height of the acrylic components is slightly smaller than the height of the slabs (8 mm vs. 10 mm) to avoid any contact with the interfaces.Fig. 4(c.1) and Fig. 4(c.2) show top-side views of the links, respectively, attached and out of a frictional slab.These links ensure a rigid connection between all the components, so that only the springs deform during actuation.Moreover, the ball bearings ensure that the force transmission remains horizontal as the lever rotates.The steel junction between the lever and force sensor is also attached via analogous links to remain horizontal at all times.
The experimental set-up, with its springs and PMMA slabs, corresponds to the compliant driving limit with Φ ≃ 6 × 10 −5 , of the same order as for the compliant regime in the numerics.Indeed, Φ ≡ K/K s , with K s = AG/h the shear stiffness of the slabs and G ≡ E/(2(1 + ν)), with, for PMMA, Young's modulus E = 2 GPa and Poisson's ratio ν = 0.3.
In addition to the global force measurement F , we also measure the absolute average horizontal position of each slab x i , by tracking a red marker placed on their side (Fig. 4(e.1)), from photographs taken at a rate of 5 fps using a digital camera (Flea3 FL3-U3-20E4C, Flir, lin-ear pixel size: 70 µm).Using a color threshold, the red markers are extracted, as shown in Fig. 4(e.2), and their geometric center is located (red crosses).Owing to the large size of the markers (154 mm 2 , corresponding to approximately 30000 pixels), x i can be determined with a sub-pixel resolution of 5 µm.The relative displacement between slabs is R i ≡ x i −x i−1 (see Fig. 1a), which serves as a proxy for the total slip at the interface i (neglecting the shear deformation of the slabs).To vary the number of sliding interfaces n, we keep the same number of slabs (5) but remove 5 − n springs, starting from the top (see Fig. 1b where n = 4).This procedure ensures robust image detection and reduces external contamination of the interfaces by keeping them in contact.Each time the slabs are disassembled to vary n, the interfaces are cleaned with isopropanol and quickly dried using compressed air.For each value of n, we perform 10 runs during which we drive over a range ∆γ = 0.6 rad, starting at γ = −0.30rad, each time excluding γ between −0.3 rad and −0.27 rad (300 s) to ensure measuring in a steady state.After each run, the lever is reset back to γ = −0.30rad.On average, each connected layer is forced to move by a total relative distance of R tot = h∆γ = 6 mm during a run that lasts 6000 s in total.

B. Measurements and lever kinematic
For a stack with n = 4, we present in Fig. 5a a typical time series extract of the force F (t) required to actuate the lever (top left plot), together with the corresponding relative position of the slabs R i (t) (bottom-left plot).The experiments exhibit stick-slip, with stick periods when the slabs are immobile (R i ≈ constant) and F increases monotonically, punctuated by macroscopic slip events.These slip events are identified by a sudden position jump, ∆R a > 0 (with a denoting the sliding interface), accompanied by an abrupt force drop ∆F > 0, cf.Fig. 5a.On all occasions, we find that only one layer slips at a time, recovering similar dynamics as in the numerical model in the compliant driving regime with a similar value for Φ (Fig. 3b).However, we note that during the stick periods, we observe what seems to be 'slow slip' where an interface moves gradually, leading to a non-linear force response.These are out of our primary focus but are discussed at the end of the section.
For each value of n, we acquire an ensemble of at least 100 slip events per layer, such that the slip quantities associated can be represented as probability distributions.For example, in Fig. 5b, we show the probability distribution of force drops, P (∆F ), occurring on the interface a = 1, for all cases of n considered.Starting from the peaked distribution for n = 1, as n increases, the distributions broaden and take higher average values.
In contrast with more classic stick-slip experiments with a single interface [32], the global measure ∆F is not a direct quantification of the frictional properties of the interface but couples with the specific kinematic of the lever.Still, the fact that only one interface slips allows us to extract a jump in a friction-like quantity ∆µ a (or stick-slip amplitude) from ∆F .We define the friction coefficient µ i of an interface i as the horizontal force acting on this interface divided by the normal force.Considering the horizontal force balance on an interface i, the interface has to resist the combined forces of the pulling springs of the slabs j ≥ i, such that where f j is the force due to the driving spring on slab j (see Fig. 1a for a visual representation of µ i and f j ).
When an interface a slides by ∆R a , the relative positions of the other interfaces remain unchanged (we only observe single-slip events).Consequently, the absolute horizontal position x i of the layers above a is increased (a) Top plot: extract of a time series of the macroscopic force F (t) for a system of n = 4 frictional interfaces.The color of the force drops ∆F follows the color code in Fig. 1 and indicates the index a of the slipping interface.Bottom plot: corresponding relative displacement (total slip) Ri(t) − Ri(t0) (with t0 = 5400 s an arbitrary value) of each interface i.Each slip event is characterized by ∆Ra.We denote Ta, the time between subsequent slip events on the same interface.Note that we show i = 5 only for completeness, by definition, R5 = 0 if n = 4.(b) Probability distribution function P (∆F ) for slip at interface a = 1 for an increasing number of layers n.(c) Comparison between a direct measurement of ∆Ra, and the computed ∆Ra(∆F ), obtained through Eq. ( 5), for each detected slip event.
by ∆R a : This sliding induces a drop in the spring forces: ∆f i = K∆x i = K∆R a for i ≥ a, and ∆f i = 0 for i < a.Note that for consistency, we define ∆f i to be positive.From Eq. ( 1), we can then express ∆µ a as a function of the slip distance ∆R a : We proceed by linking ∆R a to the global force drop ∆F , using the fact that only one interface slips at a time.Through moment balance on the lever, we obtain: Combining Eqs. ( 2) and ( 4), we obtain a relation between the global quantity ∆F and the local ∆R a : (5) (Note that i is the only varying term in the sum and n i=a i = (n(n + 1) − a(a + 1))/2 = (n + a)(n − a + 1)/2.)This result is verified in Fig. 5c.Indeed, a direct measurement of ∆R a is very close to the inversion of Eq. ( 5), in which ∆R a follows from the measured ∆F without any fitting parameter.
Finally, we combine Eqs. ( 3) and ( 5) to obtain the sought relation between ∆µ a and ∆F : We have thereby disentangled the friction properties of the interface from the kinematics of the lever.Using Eq. ( 6), we can now obtain a measure of the stick-slip amplitude of the interfaces ∆µ a , extracted directly from the global force ∆F .The position measurements R i are used only to identify the slipping interface a.
The lever kinematics introduces a strong coupling between the interfaces.Via Eqs. ( 1) and ( 2), a slip on an interface a will induce drops of friction coefficient µ i on all the other interfaces, even if no slip occurs on them.
Central experimental result.Next, we assess the effect of having multiple sheared interfaces on their frictional properties.Fig. 6 shows the probability distributions P (∆µ a ) associated with the different sliding interfaces a (different panels) and the increasing number of total active interfaces n (different colors).Each interface is compared to its response when sliding individually (n = 1 in black, see Appendix A for experimental protocol).For all the interfaces, the stacks exhibit significantly enriched statistics when compared to a single sliding layer (n = 1), in contrast to the numerical predictions reported above (Fig. 3c) where ∆µ a was independent of n.With increasing n, the location of the major peak of P (∆µ a ) shifts to lower values of ∆µ a and the respective distributions become broader as secondary peaks emerge.

C. Interpretation
We seek to interpret the above experimental findings evidencing a variation of the frictional properties as more layers are added to the stack (Fig. 6), whereas they are independent of n in the numerics (Fig. 3d).Our methodology is to identify the possible mechanisms contributing to variations of ∆µ a and quantify their potential signature in the physical and numerical experiments.First, we will attribute the finite width of the peaks in the P (∆µ a ) distributions, even for n = 1, to the interface disorder also present in the numerics.Then, we will argue that the major changes with n observed in P (∆µ a ) result from the combined effect of increasingly complex loading paths and aging of the interfaces.Finally, we speculate how creep of the interfaces, at the origin of ag- ing and slow slip, could be influenced by the increase of an effective temperature with n.
Interface disorder.Even when sliding individually (n = 1), the frictional properties of the interfaces are distributed: P (∆µ a ) has a finite width, see black curves with circles in Fig. 6.These underlying statistical fluctuations, also present in the numerical model (Fig. 3d), are considered to be related to the disorder of the contacting interfaces.The rough interface induces a broad distribution of barriers, leading to collective events with non-trivially distributed sizes.These collective events nucleate the macroscopic slip [22,29], such that the stress at which slip is nucleated is distributed.
Let us now verify experimentally that, in the case of individually sliding layers (n = 1), the measured fluctuations of ∆µ a correspond to distinct frictional interface strength.In the individual configuration, the spring drives the layer at a constant rate ḟa = Kh γ (in practice, the shear rate imposed by the lever is adapted for each a and set to γ/a to account for the difference in height, see Appendix A).The shear applied to the interface then grows at a rate μa = ḟa /N = Kh γ/N .As such, we expect the stick-slip amplitude to be proportional to the time between consecutive slips T a , following: This is consistent with our data in Fig. 7a, thus confirming that the finite width of the primary and secondary peaks in P (∆µ a ) results from statistical fluctuations of frictional interface strength.However, these fluctuations do not account for the shifts of the peaks and the appearance of secondary peaks.
In the following, we argue that coupling via the lever is responsible for increasingly complex loading paths with n, leading to broadly distributed waiting times between slips T a (I).Then, aging of the interface, present in the experiments but not in the numerics, becomes significant in broadening ∆µ a (II).
(I) Complex loading path.The coupling of the interfaces via the driving lever has two effects.
First, the loading rate at a given interface increases with n.Using Eq. ( 1) and ḟi = Kih γ while no interfaces are sliding, an interface a in a stack of n layers undergoes a loading rate of μa = K γh(n+a)(n−a+1)/(2N ), which is an increasing function of n.With this increased loading rate with n, we thus expect the time between slips T a to typically decrease with n for all interfaces (note that this effect is not captured by our quasistatic numerics in which γ = 0).
Second, the lever couples the load on an interface to slips on other interfaces via Eqs.( 1) and ( 2).Slip on one interface unloads the interface itself, but also all the other interfaces in the stack, even if no slip occurs on them.In particular, if no slip occurs anywhere, µ a is a linear function of time and reaches its frictional strength with a time T a , which is a decreasing function with n.But as n increases, slips occurring on other interfaces decrease µ a instantaneously (with a K-and n-dependent amplitude [48]).This delays the time when the interface reaches its frictional strength, leading to various T a .
In Fig. 7b, we plot the probability distribution function of T a for a = 1 and increasing n.We indeed recover that the peaked distribution for n = 1 shifts to lower values of T a with increasing n (when loading the interface without interruption by other slips).Moreover, secondary peaks in P (T a ) start to appear, which we interpret to be due to slip events on the other layers.These observations are robust for the other interfaces (see Appendix B).The change in loading rate with n, together with the complex loading path, allows our experimental system to probe a broad distribution of T a on all its interfaces.
(II) Aging.It is a known experimental fact that the macroscopic stress required for the onset of sliding, characterized by µ s (the 'static friction coefficient' in Amontons-Coulomb's terminology [11,33]), depends on the duration T that the interface was static: µ s = B ln(T /T 0 ) [17,[32][33][34], where the aging rate of the interface B remains a constitutive parameter and T 0 a microscopic time-scale.Let us now consider ∆µ a as a proxy for µ s , assuming that a slip event unloads the interface to a well-defined and constant quantity (µ d the 'dynamic friction coefficient'), as is supported by [32].We expect to find, for each sliding interface a and over a wide range of T a , that the stick-slip amplitude follows: Thereby T 0 is an unknown microscopic time-scale that we take equal to one second following common experimental literature (summarized e.g. in [33]).In Fig. 7c, we assess this expectation experimentally by plotting ∆µ a versus T a in a semi-logarithmic scale.To capture the general trend, we bin logarithmically T a , corresponding to the black markers with corresponding error bars in Fig. 7c.These averaged values of ∆µ a are indeed consistent with a straight line in the semi-log plot, and we extract the slope B = 0.053 ± 0.005, which is of the same order of magnitude as measured in classical stop-and-go experiments that are in the 10 −2 order [33], and of a direct surface observation on PMMA at room temperature that reports B = 0.009 ± 0.001 [47].Aging of the interface's contacts then translates the large peak shifts and the emergence of new ones in the T a distributions (Fig. 7b) into qualitatively similar changes in the ∆µ a distributions as observed in Fig. 6.In Fig. 7d, we schematically represent the coupled role of (I), (II), and disorder, following the evolution of the interfacial stress of an interface µ a with time, starting from the last slip event.The layer slips when it reaches µ a = µ s (T a ), whereby µ s is distributed in some way for fixed T a because of disorder (illustrated as a red-shaded area, where for simplicity, we lump all fluctuations in the threshold to sliding) and increases logarithmically with time because of aging.For n = 1 (black line), µ a increases linearly at the same rate for all the events, thus exploring a narrow region of µ s (the shaded red region due to disorder) and T a , following Eq.( 7).In the case of multiple active interfaces (n > 1, green lines), µ a increases faster given that μa is an increasing function of n, and several scenarios arise.If no other slip events occur in the stack, µ a linearly reaches µ s , resulting in a lower value of T a and ∆µ a .However, if sliding events occur elsewhere during loading, µ a will drop before linearly increasing again, delaying slip and thus increasing T a , and consequently ∆µ a , because of aging.
Creep effective temperature.During stick intervals, microscopic events occur on the interfaces, propagating elastic waves across the system [29].As we increase the number of interfaces in the system, we can expect that the overall mechanical noise created by the microscopic events also increases.If we speculatively interpret this mechanical noise as an effective temperature, we would expect a change of aging rate B with n.Let us define the aging rate for a single event as B a ≡ ∆µ a / ln(T a /T 0 ) and associated probability distributions for all associated events (see Fig. 8a for a = 1, and Appendix B for the other interfaces).Although the mean of P (B a ) does not change with n, we do find that the width of the distribution P (B a ) is an increasing function of n mainly for the lowermost interfaces (i ≤ 2), as shown in Fig. 8b.
For the same interfaces (i ≤ 2), we also observe distinctly different slip dynamics when n > 1.In particular, as n increases, we find that interfaces i = 1 and i = 2 are increasingly subject to slow slip, defined as sliding significantly slower than the slip events (see Fig. 8c for an example and Appendix D for a quantitative characterization).Slow slip is not accompanied by a sudden macroscopic stress drop but rather by a slow decrease of Ḟ .In Fig. 8d, we estimate the fraction of slow slip R slow compared to the total sliding distance R tot .It is computed by comparing R tot to the total accumulated slip during individual slip events (sum of all slip events :slip elsewhere  7).For multiple sliding interfaces (n ≥ 1): (b) Probability distribution of the waiting time Ta between two consecutive slip events at interface a = 1, for increasing n.(c) For each detected slip event, correlation between ∆µa, and its corresponding waiting time Ta (semi-logarithmic scale).The black markers correspond to the mean values for a logarithmic binning of Ta (error bars indicate the standard deviation for that bin), and the dotted line a fit of Eq. ( 8) with B = 0.053 ± 0.005.Alternative representations of this data-set are displayed in Appendix C. (d) Schematic of the proposed mechanism leading to multimodal and wider distributions of Ta (and thus ∆µa) as n increases.
∆R a ), such that R slow ≡ 1 − ∆R a /R tot .Once a slow slip event starts, it appears to be stopped only by slip events occurring either on the same interface or on any other interface.An increase in the effective temperature of the interface with n could also act as a potential destabilization factor of the contacts at the interfaces, increasing the occurrence of slow slips with n.

A. Summary
We have explored the stick-slip response of a system with multiple interfaces by proposing a model system comprising n vertically stacked slabs, each connected to a lever whose rotation is imposed.The interfaces were driven in quasistatic (homogeneous) shear.We proposed a dimensionless quantity Φ as the ratio between the driving stiffness and the elastic shear stiffness of the slabs.We have argued and demonstrated numerically that the system displays synchronization if Φ is sufficiently large (Φ ≳ 10 −3 ).In that case, the system acts close to a single frictional interface with effective properties.If Φ is small (Φ ∼ 10 −6 ), interfaces slip one by one, as also confirmed experimentally.We expect non-trivial collective effects with increasing n only in the low-Φ limit, which we addressed through experiments.In the numerics, the stick-slip amplitude of the interfaces ∆µ a displays a distribution with finite width because of statistical fluctuations of the interfaces, but no measurable changes with n.By contrast, we measured experimentally that the probability distribution of stick-slip amplitude ∆µ a shows a general broadening with n, with peaks shifting to lower values and secondary peaks appearing.The interfaces are coupled via the lever, exposing them to a complex loading path, and leading to a broad distribution of the waiting times T a between two slip events on an interface.We find that T a is now spanning two decades, such that the aging of the interfaces plays a crucial role in the broadening of ∆µ a .The complex distributions of ∆µ a can then be interpreted as the combined effect of interface disorder, also observed numerically, and aging.For narrow waiting times T a , multiple slips explore the statistical fluctuations of the contacting interfaces, giving the finite width of the peaks in the distributions.In addition, creep-induced aging gives a long-time general trend over widely distributed T a .
Furthermore, we observe that increasing n has a significant impact on creep-induced phenomena like aging rate and slow slip.We suggest that these additional consequences of adding more layers to the stack might be evidence of an increase in an effective interface tempera-ture due to mechanical noise of microscopic events.
In conclusion, the relative rigidity of the drive against the layers dictates whether a stack of interfaces responds or not.When layers slide one by one, increasing their number leads to a complex coupling, making the prediction of the next slip more challenging.

B. Limitations and outlook
It is pertinent to discuss some limitations of our model system and provide suggestions for future work.
Stiffness ratio.We have defined the stiff and compliant driving regimes, as characterized by the relative order-of-magnitude estimation of the respective stiffness ratio Φ. Identifying an equivalent of Φ in systems with more intricate geometries, such as fault networks, might contribute to clarifying their dynamics and help slip predictions.Hence, a more systematic exploration of the response of stacks with varying Φ would be of great interest.
However, we are currently restricted to a limited range of Φ.For our numerics, low-Φ are challenging due to a combination of the assumption of small rotations and finite machine precision.To be able to continue sliding indefinitely, our model should be extended with the possibility to reset the local deformation along the frictional interface to a zero average while keeping the identical stress state.In contrast, high-Φ are challenging experimentally with the current setup.If the driving stiffness is increased, the motor/lever system can no longer be considered rigid, invalidating the relation between global (∆F ) and local (∆µ a ) force drops (Eq.( 5)).Instead, if the shear stiffness of the slab materials is significantly lowered, a different class of frictional properties is expected as soft materials tend to be adhesive [49].
Creep.Our proposed model system allowed us to measure the aging rate B of the interfaces thanks to complex stick-slip sequences without the need for stop-andgo experiments.However, while our measured value of B is compatible with previous experiments on PMMA [33,47], it differs by a factor of five.Possible sources of differences are roughness, inter-realization variations, and stress inhomogeneities.First, the PMMA plates used here have a much higher surface roughness than in [33,47].If our model is extended with thermal fluctuations it likely displays creep such that the relationship between the distribution of barriers, linked to surface roughness, and creep could be investigated.Second, we measure B on an ensemble of n = 5 interfaces.Between interfaces it is estimated that B differs by a factor of about two.Third, recent experimental observations find a relationship between the aging rate B and the applied shear load [50].Our setup naturally imposes a broad range of shear loads on the interface.However, to measure the empirical law by [50], our setup would need to be augmented to also provide stress measurements per interface by measuring individually the internal forces f i of the driving springs.To study this effect numerically, on top of adding temperature to capture creep, the model would likely have to be made sensitive to pressure inhomogeneities that may arise from the normal load or partial slip events, commonly thought to be a significant perturbation [43].
Slow slip.The origin of the experimentally observed slow slip is unknown.A tempting hypothesis is that smooth sliding is the result of activated yielding events due to increasing mechanical noise on other interfaces.However, it is not clear why this interpretation would lead to slow slip occurring predominantly on the lowermost interfaces (which is qualitatively robust to changing the order of the slabs, see Appendix D).Furthermore, slow slip is not observed numerically.This could, however, in part be due to our small homogeneous background damping term (currently chosen to avoid nonphysical periodic wave propagation).An alternative hypothesis is that, by increasing n, the loading rate becomes sufficiently high to drive the interfaces away from the stick-slip regime (if v ≥ v c ).This second interpretation is consistent with slow slip being predominantly observed on the lowermost interfaces.Note that, differ-ent from the experiments, our numerical model drives infinitely slowly.
Aftershocks.Aftershocks appear if creep is added to the drive of a simple spring-block model [51,52].A creeping drive is often associated with the high temperatures in Earth's core [51].Our experimental system displays slow sliding of 'deep' layers already at room temperature.A key question is if aftershocks appear in the top layers of our system as well.Answering this question experimentally would require exposing microscopic events, which would likely involve studying acoustic emissions (for which PMMA may not be the optimal choice).All reported stresses have been normalized by the typical yield stress.For the numerical results, γ is normalized by γ 0 , which is defined as the shear corresponding to a yield event in a system in which all blocks of the first interface (i = 1) have a yield strain equal to the typical yield strain.

F. Slip sequences -numerics
We argue that the waiting time between slip on an interface, T a , is a crucial quantity.Since the numerics run a quasistatic protocol T a is not defined.The closest equivalent quantity is the increment is the shear ∆γ a between two slip events on a layer a.We plot this quantity in Fig. 14 for a = 1 for both stiff and compliant driving.This plot should be qualitatively compared with Fig. 7b of the main text.However, only a coarse qualitative comparison is possible due to the following differences between the numerics and the experiments.Firstly, the numerics follow a quasistatic loading protocol such that γ ≡ 0 and is independent of n.Secondly, numerically we use a force density to control the mean position of each plate.Consequently, the applied force does not increase with n.Lastly, the simulations are athermal and hence without aging.What remains is that the higher fraction of single-slip events at lower K leads to a slightly differently distribution P (∆γ) for n > 1 in Fig. 14b.A better comparison with Fig. 7b would require a lower K, currently outside our reach.In addition the relationship between T a and ∆γ a would have to be calibrated for the quasistatic protocol for each n.
We comment on the fraction of events that are part of a sequence of single-or multi-slip events using Fig. 15.For stiff driving, in Fig. 15a, sequences of single-and multislip events alternate.As n increases, such sequences alternate more often as reflected by the decreasing fraction of events that are part of any sequence (blue line).Of these sequences, there is an equal partitioning in the number of single-and multi-slip sequences (the two highlighted regions have an approximately equal height for every n).For compliant driving, the situation is sim-ple: there are almost exclusively single slip events, see Fig. 15b.

G. Front dynamics -numerics
We briefly discuss the dynamics of a multi-slip eventin which multiple layers slip during a macroscopic event.
For simplicity we consider a system with n = 2. Figs.16a and 16b show event maps of two typical system-spanning events (in which all blocks of each interface fail at least once).As observed, the applied infinitesimal increase of applied shear triggers an avalanche on one of the layers.This avalanche propagates as a fractal object in spacetime until it transitions to a ballistic event that propagates through the entire system [22,29] with a front velocity c f of about twice the shear wave speed c s (see fit in Fig. 16a).With a delay of the order of the time a shear wave takes to propagate through the layer above or below (c s /h, with h the height of each layer), the other layer starts to fail.In contrast to the first slipping layer,   Fraction of slip events that are part of a sequence of single-or multi-slip events for (a) stiff and (b) compliant driving.The plot is constructed such that the blue line (square markers) corresponds to the total fraction of events that is part of a sequence of consecutive single-or multi-slip events.The red line (circular markers) is the fraction of events that are part of sequences of single-slip events.The area between the red and blue curves, highlighted in blue, thus corresponds to the fraction of events that are part of a sequence of multi-slip events.
by chance avalanches can nucleate at multiple locations (e.g.Fig. 16a).These avalanches start to propagate independently before linking-up.We note that nucleating a single avalanche on the secondary slipping layer is also frequently observed (e.g.Fig. 16b).We quantify the delay of activity on the secondary slipping layer.We measure the duration τ between the first event on the first layer with activity and the first event on the other layer.
The distribution of τ is shown in Fig. 17, confirming that the delay is of the order of the time it takes a shear wave to propagate through a layer.

Figure 1 .
Figure 1.(a) Schematic of our model system, shown here with n = 4 active (driven) layers.The color code of the layers is used throughout the figures.(b) Photograph of the corresponding experimental apparatus.

Figure 2 .
Figure 2. Schematic of the loading and consecutive slip of the interfaces in the compliant (a) and stiff (b) regimes.

3 Figure 3 .
Figure 3.Numerical results: (a, b) Typical steady-state global stress Σ response as a function of applied shear γ for n = 3 for (a) stiff driving and (b) compliant driving.We indicate all (microscopic) yielding events with a black dot marker.Slip events on a single layer are indicated in color (see legend), while slip events in black involve more than one interface.(c) The fraction ρ(s) of macroscopic slip events involving s = 1, . . ., n layers, for stiff (dashed) and compliant (solid) driving; see legend in (d) for color-map and markers.(d) Distribution of stress drops at the slipping interface for different n in the compliant regime (for slip events on a single layer, for which s = 1 in (c)).See the main text for definitions and units.

Figure 4 .
Figure 4. Technical details on the main components of the experimental apparatus.Top: overview of the set-up with the location of the different parts: (a) force measurement, (b) lever-spring link, (c) spring-slab link, (d) frictional slabs surface and (e) position measurement.Relevant details are shown in the bottom panels (numbered (a.1), etc.), see text for details.
Figure 5.(a) Top plot: extract of a time series of the macroscopic force F (t) for a system of n = 4 frictional interfaces.The color of the force drops ∆F follows the color code in Fig.1and indicates the index a of the slipping interface.Bottom plot: corresponding relative displacement (total slip) Ri(t) − Ri(t0) (with t0 = 5400 s an arbitrary value) of each interface i.Each slip event is characterized by ∆Ra.We denote Ta, the time between subsequent slip events on the same interface.Note that we show i = 5 only for completeness, by definition, R5 = 0 if n = 4.(b) Probability distribution function P (∆F ) for slip at interface a = 1 for an increasing number of layers n.(c) Comparison between a direct measurement of ∆Ra, and the computed ∆Ra(∆F ), obtained through Eq. (5), for each detected slip event.

Figure 7 .
Figure 7. (a) For an interface sliding individually (n = 1), stick-slip amplitude ∆µa as a function of the waiting time since the last slip event Ta.The black line corresponds to the prediction of Eq.(7).For multiple sliding interfaces (n ≥ 1): (b) Probability distribution of the waiting time Ta between two consecutive slip events at interface a = 1, for increasing n.(c) For each detected slip event, correlation between ∆µa, and its corresponding waiting time Ta (semi-logarithmic scale).The black markers correspond to the mean values for a logarithmic binning of Ta (error bars indicate the standard deviation for that bin), and the dotted line a fit of Eq. (8) with B = 0.053 ± 0.005.Alternative representations of this data-set are displayed in Appendix C. (d) Schematic of the proposed mechanism leading to multimodal and wider distributions of Ta (and thus ∆µa) as n increases.

Figure 8 .
Figure 8.(a) Distribution of the aging rate Ba = ∆µa/ ln(Ta/T0) for all events on the interface a = 1 and all n.(b) Standard deviation of Ba for all layers a, normalized by that quantity for n = 1, as a function of n.(c) For n = 5, a typical time evolution displaying slow slip of the relative position Ra(t) − Ra(t0), of the layer a = 2.The sliding distance Rtot is categorized as slip events of amplitude ∆Ra, and slow slip.The fraction of slow slip is defined by R slow = 1 − ∆Ra/Rtot.(d) Fraction of slow slip R slow for all layers a as a function of the number of layers n.Error bars represent the minimal and maximal values of R slow obtained over 10 experimental runs.

Figure 9 .
Figure 9. (a) Verification of ∆Ra(∆F ) (Eq. (A1)) for individual interfaces.(b) The probability distributions of the drop of force ∆F for the different interfaces a are shown using a different color and marker.(c) Probability distributions of the corresponding drop in frictional resistance ∆µa as estimated from ∆F using Eq.(A2).

Figure 11 .Figure 12 .
Figure 11.For each slip event, its stick-slip amplitude ∆µa as a function of waiting time Ta, corresponding to the duration between the current slip event and the last slip event on the same interface.The dotted line in all the panels is the aging law ∆µa = B ln(Ta/T0) (Eq.(8), dotted line) with B = 0.053.(a) Reproduction of Fig.7c, where each event is colored by the interface a it occurred, and its symbol reflects the number of interfaces n.The black markers are the average ∆µa per bin of Ta, used to fit the aging law.(b) The same data as previously are presented, with this time colors indicating n, and different symbols the sliding interface.(c) On top of the previous data, the event measured on the interface sliding individually are included using black markers.Over this narrow range of Ta, disorder at the interface is sampled following the proportionality given by Eq. (7) (solid line).(d) The data of (a) is biased to only events where there was no detected slow slip after the last slip event on the same interface.

Figure 13 .
Figure 13.Finite element mesh of the numerical model for n = 4.The colors indicate the different layers (color coding is the same as in Fig.1).The positions of the bottommost nodes are fixed as indicated.The shear γ is imposed by a parabolic potential of curvature K associated to the mean position of each layer (the representation of the lever and springs is thus schematic).

Figure 14 .
Figure 14.Shear increment ∆γa between slip events on the first interface (a = 1) as a function of n for stiff driving (a) and compliant driving (b).Note that we consider any slip event involving a = 1, i.e. also those in which multiple interfaces slip.

Figure 15 .
Figure15.Fraction of slip events that are part of a sequence of single-or multi-slip events for (a) stiff and (b) compliant driving.The plot is constructed such that the blue line (square markers) corresponds to the total fraction of events that is part of a sequence of consecutive single-or multi-slip events.The red line (circular markers) is the fraction of events that are part of sequences of single-slip events.The area between the red and blue curves, highlighted in blue, thus corresponds to the fraction of events that are part of a sequence of multi-slip events.

Figure 16 .
Figure 16.(a, b) Event-map of two typical simulations for n = 2 .Each panel contains two plots above each other: one for each interface i.In each plot, the horizontal axis is the horizontal position r normalized by the size of one block along each of the weak layers ℓ0.The vertical axis is the time t since triggering the first failure anywhere in the system, normalized by the time it takes a shear wave to travel the height h of each layer (cs/h, with cs the shear wave speed).A point is drawn every time there is a failure (failure in the direction of applied shear in black, and in the opposite direction in blue).

4 PFigure 17 .
Figure 17.Probability distribution of the duration τ between activity on the primary and secondary slipping layer for n = 2, normalized by the time it takes a shear wave to travel one layer.