The Moment Duration Scaling Relation for Slow Rupture Arises From Transient Rupture Speeds

The relation between seismic moment and earthquake duration for slow rupture follows a different power law exponent than subshear rupture. The origin of this difference in exponents remains unclear. Here, we introduce a minimal one‐dimensional Burridge‐Knopoff model which contains slow, subshear, and supershear rupture and demonstrate that different power law exponents occur because the rupture speed of slow events contains long‐lived transients. Our findings suggest that there exists a continuum of slip modes between the slow and fast slip end‐members but that the natural selection of stress on faults can cause less frequent events in the intermediate range. We find that slow events on one‐dimensional faults follow M¯0,slow,1D∝T¯0.63 with transition to M¯0,slow,1D∝T¯32 for longer systems or larger prestress, while the subshear events follow M¯0,sub‐shear,1D∝T¯2 . The model also predicts a supershear scaling relation M¯0,super‐shear,1D∝T¯3 . Under the assumption of radial symmetry, the generalization to two‐dimensional fault planes compares well with observations.


Introduction
Over the last decades, an increasing number of slow slip events on faults have been reported (Bürgmann, 2018). A measure that is viewed as a key to unraveling the mechanism of slow and fast rupture is the relation between seismic moment M 0 and slip event duration T. Regular fast earthquakes have long been known to follow a moment duration scaling relation of M 0 ∝T 3 . Ide et al. (2007) suggested that slow events follow a unified scaling relation M 0 ∝T. They suggested that the linear relation between moment and duration for slow events can be explained in two ways: (1) the average slip is proportional to the fault length as for fast propagation, and the stress drop is constant for all events, which gives the relation M 0 ∝T. (2) The slip amount is constant for all events, and the fault area increases linearly with time L 2 ∝T, which results in M 0 ∝T. Peng and Gomberg (2010) elaborated on the ideas of Ide et al. (2007) and reached a different conclusion that rupture should span a continuum between fast and slow velocity end-members. However, almost 10 years after the suggestion of Peng and Gomberg (2010), observations of events between the fast and slow end-members are still sparse. Later studies have reported on a variety of scalings between moment and duration ranging from M 0 ∝T to M 0 ∝T 2 (Aguiar et al., 2009;Frank et al., 2018;Gao et al., 2012;Ide et al., 2008), although the definition of events can vary somewhat between different studies.
The shape of slip patches can also influence the observed scaling. Ben-Zion (2012) argued that fractal slip patches can result in a scaling relation M 0 ∝T 2 =logðTÞ because the average displacement is approximately ©2019. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. constant rather than proportional to the rupture dimension. Bounded propagation can also play an important role (Ben-Zion, 2012;Gomberg et al., 2016). Gomberg et al. (2016) suggested that the scaling relation between moment and duration is the same for slow and fast events but that a transition occurs between a two-dimensional scaling and a one-dimensional scaling when the rupture propagation switches from unbounded to bounded in one direction. Assuming that the fault displacement can be approximated using dislocation theory, this results in a transition from T 3 to T. They suggest that there should be a bimodal but continuous distribution of slip modes and that a difference in scaling relations alone does not imply a fundamental difference between fast and slow slip. The above-mentioned theoretical considerations implicitly assume constant rupture velocity. However, this contradicts observations by Gao et al. (2012) that show that the average rupture speed for slow events decreases with increasing seismic moment, which is a strong indication of transient rupture speeds.

10.1029/2019GL084436
Slow slip events emerge in numerical models with rate-and-state friction. Colella et al. (2011) simulated a Cascadia-like subduction zone using rate-and-state friction. They analyzed a large number of slip events and found that the seismic moment M w scales as M w ∝T 1.5 for M w ≤ 5.6 with a transition to M w ∝T 2 for M w >5.6. Shibazaki et al. (2012) modeled the subsuction zone of southwest Japan with rate-and-state friction. For slow events, they found a scaling M 0 ∝T 1.3 . Liu (2014) used rate-and-state friction on a 3-D subduction fault model and found a scaling M 0 ∝T 1.85 . Romanet et al. (2018) highlighted the role of interactions between faults. They argue that the scaling relationships of slow slip events and earthquakes emerge from geometrical complexities due to interactions between fault segments. The moment duration scalings have not only been addressed using rate-and-state friction. Ide (2008) introduced a Brownian walk model for slow rupture, where the assumption is that there is a random expansion and contraction of the fault area, so that its radius can be described as a Brownian walk with a damping term. This model predicts M 0 ∝T for large T.
Here, our goal is to answer the following two questions: (1) Is there a separation of two distinct classes (Ide et al., 2007), or is there a continuum of possible scaling relations between the fast and slow end-members (Peng & Gomberg, 2010)? and (2) Can a difference in M 0 -T scaling relations alone be attributed to different physical mechanisms behind slow and fast rupture? We address both (1) and (2) through a Burridge-Knopoff-type model with Amontons-Coulomb friction with a velocity strengthening friction term that has previously been shown to contain a large variety of rupture phenomena, including subshear, supershear, and slow propagation . Velocity strengthening friction has been shown to be a generic feature of dry friction (Bar-Sinai et al., 2014) and has been reported in Halite shear zones at low slip speeds or large confining pressures (Shimamoto, 1986). The friction law can also be interpreted as a transition from a dry contact to a lubricated sliding regime with increasing velocity (a Stribeck curve) under the additional assumption that the transition from dry contact to lubricated sliding occurs at a small sliding speed (Gelinck & Schipper, 2000;Olsson et al., 1998).
For homogeneously stressed faults, the model can be reduced to only two dimensionless parameters τ and α representing the prestress and a velocity strengthening friction term, respectively. The advantage of such approach is that the simplicity of the model allows us to calculate moment duration scaling relations both through numerical simulations and through analytical calculations. Through numerical simulations, we demonstrate that there exists a continuum of rupture modes between the slow and fast end-members but that the most likely selection of τ in nature produces two distinct classes separating subshear and slow rupture velocities. Through analytical calculations, we show that the scaling relation for slow fronts arises due to long-lived transients in the rupture velocity. Such transient rupture velocity has been observed in nature through a dependence on the average rupture speed on the seismic moment for slow fronts (Gao et al., 2012). In addition, the model predicts a separate scaling for supershear rupture not previously reported in the literature.

A One-Dimensional Burridge-Knopoff Containing Slow and Fast Rupture
We solve the one-dimensional Burridge-Knopoff model (Burridge & Knopoff, 1967) with Amontons-Coulomb friction with a linear velocity strengthening term. The dimensionless equation of motion for a chain of N blocks can be written as (a detailed derivation can be found in the supporting information) which is integrated using the Euler-Cromer method with dt ¼ 10 −3 . u is the dimensionless displacement, and is the dimensionless prestress where σ N is the normal stress, τ is the initial shear stress, and μ s and μ k are the static and dynamic coefficients of friction, respectively. The symbol ± denotes the sign of the block velocity.
For positive velocities, we only need to consider τ þ , but negative velocities can occur in a small subset of our simulations. In such situations, we need to prescribe the relation between μ s and μ k , which we set to μ s =2μ k , so that τ − ¼ τ þ þ 2. In the rest of the paper we will use τ as a reference to τ þ . The second dimensionless parameter is a viscous term, where ρ is the density, E is the elastic modulus, and α is a velocity strengthening term with units (Pa s/m). α can range from 0 to infinity, where α ¼ 0 recovers the ordinary Amontons-Coulomb friction without viscosity. τ has an upper limit of 1, where the prestress equals the static friction threshold. For τ <0, steady-state propagation does not occur (Amundsen et al., 2015). This corresponds to a prestress below the dynamic friction level. We thus simulate τ ∈½10 −7 ; 1 and α∈½10 −3 ; 10. The boundary conditions assuming triggering through soft tangential loading (small driving velocity V and driving spring stiffness K) are given by The rightmost block is fixed so that u N ¼ 0. Blocks start to move once the static friction threshold is reached, which in dimensionless units can be written as Moving blocks restick if the velocity changes sign. The system is sketched in Figure 1a. This model has previously been used to determine the steady-state rupture velocity which includes subshear, supershear, and slow rupture, as well as an arresting region at low τ and intermediate α . The steady-state front speed v c;∞ can be found exactly when α ¼ 0 (Amundsen et al., 2015). For α>0 we can use the empirical expression )

Model Results
For each simulation we predefine τ and α and each simulation consists of a single event. From each event we extract the displacement u and the fault length L, which is found from the position of the rightmost block that has ruptured. We run the simulations until all blocks are immobile or until the average velocity reaches 0.1% of the steady-state slip speed τ =α . In dimensionless units the zeroth-order moment for rupture propagation along a line is where ⟨u⟩ is the average displacement on a fault of length L. The seismic moment and the duration are measured when 99% of the total displacement has been reached. Figure 1 shows how M 0;1D and the event duration T are measured in the simulations. For each simulation, we measure the duration, as well as the fault length L and the average block displacement ⟨u⟩. This results in a single point in the ðM 0;1D ; T Þ diagram for each event. . We also measure the block displacement uðl; tÞ (f), from which we extract the average displacement ⟨u⟩ (g). Using equation (6), we obtain the seismic moment and the duration of the event marked with a star in (h). of Peng and Gomberg (2010). The model also predicts a second scaling relation for supershear rupture, which is found at large τ , that has not previously been reported (Figure 2b).

Origin of Scaling Relations-Analytical Calculations
The simplicity of the model allows an analytical treatment of several aspects which helps explain the various scaling relations between seismic moment and event duration. We summarize the analytical predictions for slip, front speed, and event duration and explain why the different scaling relations appear. A detailed derivation is given in the supporting information.
First, we can determine the average slip on a fault. If the stress is at the dynamic level after rupture (the stress drop equals τ ), we can calculate ⟨u⟩ exactly Equation (7) is derived for soft tangential loading, and we stress that a different boundary condition could lead to different dependencies between L, ⟨u⟩, and τ. Combining equation (7) with equation (6), we find that the seismic moment can be written as which only depends on the prestress τ and the length of the fault L. To obtain the moment duration scaling relation, we need to determine LðT Þ and thus have to combine equation (8) with information about the rupture propagation and the afterslip (i.e., the amount of slip after the propagation has stopped).
A key observation on the rupture propagation is shown in Figure 3. While fast fronts exhibit short transients and quickly reach the steady-state propagation speed given by equation (5)

Fast Regime
The short transients in the fast regime indicate that we can approximate the fault length by where t rupture is the time it takes for a rupture front to reach the end of the fault. The time it takes to arrest completely depends upon the existence of a backward propagating arresting front. If we assume that this backward propagation occurs at roughly the same speed as the forward propagation, we obtain so that This relation implies that there is a separate scaling for subshear (τ →0) and supershear rupture (τ →1) in our simulations Note that we also predict that M 0;sub -shear;1D will transition to a T 3 scaling for large L if τ is small but nonzero. The moment duration, in the fast regime using equation (11), is shown in the two bottom lines of Figure 4. This figure also shows numerically obtained values for the moment duration. The agreement between the numerical simulations and equation (11) is good.

Slow Regime
For slow fronts, v c ðtÞ is transient, and we observe that v c ðtÞ is well described by a power law with exponent β for large α and small τ . To obtain an approximation for v c ðt; α; τ Þ, we plot v c ðtÞ for a selection of τ and α in Figure 3. All curves collapse when we subtract the steady-state front velocity v c;∞ and scale with the initial rupture velocity v c;0 given in equation (S32). We can then write down v c;slow ≈ðv c;0 −v c;∞ Þðv c;0 tÞ −β þ v c;∞ (14) Figure 3 shows that this relation fits well with simulations for small τ and large α, and we measure empirically the exponent β=0.685±0.002 (using τ ¼ 10 −3 and α ¼ 10). To obtain a parametric equation for M 0 and T , we need to find T ðLÞ. T has two main components: the time it takes to rupture a fault of length L , t rupture , and the time it takes for all motion to stop. Unlike for fast fronts, the arresting phase in the slow regime is not governed by a backward propagating arresting front but rather a slow exponential decay in velocity. We denote this time t afterslip and define t rupture can be found from equation (14) (17) and (18)), while the bottom two use the fast scaling relation (equation (11)). The colored regions highlight the different scaling exponents discussed in the text.
The afterslip time can also be found analytically, and the detailed calculation is given in the supporting information. The result is where logð100Þ indicates that we take the time when 99% of the slip has been accumulated (which is necessary because the afterslip is exponentially decaying). The prediction of seismic moment versus duration can then be found using equation (8) for the seismic moment, equation (17) for t rupture (this has to be solved numerically for nonzero τ), and equation (18) for t afterslip , with T ¼ t rupture þ t afterslip . The excellent agreement between the analytical approach and the numerical simulations is demonstrated in Figure 4.
We can determine the bound on the slow front scaling relation by noting that for infinitesimal τ, v c;∞ ≈0 and . This yields where the first term will dominate over the second term (negligible afterslip) for large L because 1 1−β >2 for the measured β=0.685±0.002. We can then solve for with 2(1−β)≈0.63. We also observe a transition to a different scaling at large M 0;1D when τ is nonzero. To obtain the exponent in this regime, we note that in this limit the steady-state rupture velocity is reached, so that For large L and nonzero τ, the afterslip will dominate, so that L∝T 1 2 . The cubic term in equation (8) This means that the moment duration scaling relation in the slow regime is expected to follow a power law with exponent 2(1−β) with a possible transition to 3 2 at large M 0;1D .

Discussion
We have demonstrated that a simple Burridge-Knopoff model with Amontons-Coulomb friction is capable of reproducing the range of power law scaling relations between seismic moment and duration observed in nature. The simplicity of the model means that we can calculate the scaling relations analytically, and we find the one-dimensional exponents 2(1−β) with a transition to 3/2 for large seismic moments for slow rupture, 2 for subshear rupture, and 3 for supershear rupture, where β≈0.685 is the power law exponent of the transient slow rupture velocity.
In this letter, we aimed to address two questions. First, whether there is a separation of two distinct classes or is there a continuum of possible scaling relations between the fast and slow end-members. We argue that the most likely value for τ is close to 0, which corresponds to shear stress at the dynamic level or to ruptures where the stress drop is small compared to the background stress like in faults (Shearer et al., 2006). If this is indeed the case, the moment duration scaling should contain a continuum of slip modes between the slow and fast end-members. However, because large τ would in this case be unlikely, it would result in a distinction of fast and slow scalings simply because this is more likely. This would indicate that both the interpretations by Ide et al. (2007) and by Peng and Gomberg (2010) hold in the sense that there is a continuum of slip modes, but the natural variation of τ could result in more frequent events along the end-member scalings. In our simulations, the separation into the slow and subshear scaling relations occurs spontaneously under the assumption that τ ≈0.
The second question we aimed to address was whether a difference in M 0 -T scaling relations alone could be attributed to different physical mechanisms behind slow and fast rupture. Our model contains only two dimensionless parameters, which highlights that the observed scaling relations do not necessitate complex underlying mechanisms. The same friction law with different values for the coefficients and a varying prestress can explain the entire range of scaling relations, and the slow scaling regime arises simply because slow rupture speeds are transient. We have previously shown that fast rupture is governed by inertia, while slow rupture is noninertial , which has consequences for whether the slow and fast scaling relations can be attributed to different underlying physical mechanisms. While the derivation of the scaling relations presented in this letter does not require specification of the underlying physical mechanism causing transient rupture speeds, transient rupture speeds are only observed in the noninertial regime. This suggests that the different scalings observed in the model originate because fast rupture is inertial while slow rupture is not.
To compare our results to observations on faults, it is instructive to discuss relations that would be obtained for rupture on a 2-D plane. If we can assume radial symmetry, we can use the same expression for ⟨u⟩ as in 1-D, but M 0;2D ¼ ⟨u⟩L 2 , which changes the scaling by a term L. This changes the scaling relations to where 3(1−β)≈0.945 is the exponent that is dominant for τ ¼ 0 at large L. This is remarkably close to the hypothesized exponent of 1 from observations (Ide et al., 2007). However, it should be noted that the prestress is not expected to be radially symmetric, which puts limitations on this extension. We stress that future studies should incorporate two-dimensional simulations to address these scaling relations without such limitation.
The transition from 3(1−β) to 2 also indicates that a simple linear scaling relation between seismic moment and duration for slow events is not appropriate, because it is only valid at τ ¼ 0. We find it likely that a scaling in the approximate range M 0 ∝T to T 2 should be observed for slow events, depending also on the decaying exponent β. For a constant α, this variation in the power law exponent occurs due to changes in the stress state of the interface. This is in line with observations, where different studies have reported on scaling exponents ranging from approximately 1 to 2 (Aguiar et al., 2009;Frank et al., 2018;Gao et al., 2012;Ide et al., 2007Ide et al., , 2008. From our results in Figure 3 we are in a position to explain the observed relation between average rupture speed and seismic moment (Gao et al., 2012). , which indicates that β≈0.6±0.025. Using equation (26) yields a moment duration scaling relation for slow rupture following M 0 ∝T f1:1;1:3g , which is fully consistent with their observed linear relationship between seismic moment and duration.
Here, we have assumed that propagation is not bounded. Gomberg et al. (2016) demonstrated that there will be a change from a two-dimensional scaling to a one-dimensional scaling when the rupture propagation goes from unbounded to bounded in one of the directions. While we have demonstrated that different scalings can originate without such effect, a bounded system would add a number of possible transitions in moment duration and would in principle allow for scaling relations following both the two-dimensional and the one-dimensional exponents.

Conclusion
Linear elasticity and Amontons-Coulomb friction with a viscous term is sufficient to produce a large variety in scaling exponents between seismic moment and duration. This suggests that different scaling relations for fast and slow slip events do not require complex underlying physical mechanisms. However, our findings do suggest that whether rupture is dominated by inertia or not plays an important role because fast inertial rupture fronts propagate at fairly constant speeds, while slow noninertial rupture fronts contain long-lived transients. Our findings also suggest that there exists a continuum of slip modes between the slow and fast slip end-members but that the natural selection of stress on faults can cause less frequent events in the intermediate range. We find that the subshear scaling follows M 0 ∝T 2 (which corresponds to T 3 in 2-D), while the slow scaling follows T 2(1−β) (which corresponds to T 3(1−β) in 2-D) with a transition to T 3 2 (T 2 in 2-D) for larger seismic moments depending on the prestress. β≈0.685 corresponds to the power law decay in the slow rupture velocity with time. The model also predicts a separate scaling for supershear rupture with M 0 ∝T 3 (T 4 in 2-D).