Crack to pulse transition and magnitude statistics during earthquake cycles on a self-similar rough fault

Faults in nature demonstrate fluctuations from planarity at most length scales that are relevant for earthquake dynamics. These fluctuations may influence all stages of the seismic cycle; earthquake nucleation, propagation, arrest, and inter-seismic behavior. Here I show quasi-dynamic plane-strain simulations of earthquake cycles on a self-similar and finite 10 km long rough fault with amplitude-to-wavelength ratio α = 0.01. The minimum roughness wavelength, λmin, and nucleation length scales are well resolved and much smaller than the fault length. Stress relaxation and fault loading is implemented using a variation of the backslip approach, which allows for efficient simulations of multiple cycles without stresses becoming unrealistically large. I explore varying λmin for the same stochastically generated realization of a rough fractal fault. Decreasing λmin causes the minimum and maximum earthquakes sizes to decrease. Thus the fault seismicity is characterized by smaller and more numerous earthquakes, on the other hand, increasing the λmin results in fewer and larger events. However, in all cases, the inferred b-value is constant and the same as for a reference no-roughness simulation (α = 0). I identify a new mechanism for generating pulse-like ruptures. Seismic events are initially crack-like, but at a critical length scale, they continue to propagate as pulses, locking in an approximately fixed amount of slip. I investigate this transition using simple arguments and derive a characteristic pulse length, Lc = λmin/(4π α) and slip distance, δc based on roughness drag. I hypothesize that the ratio λmin/α 2 can be roughly estimated from kinematic rupture models. Furthermore, I suggest that when the fault size is much larger than Lc, then most space-time characteristics of slip differ between a rough fault and a corresponding planar fault. 1

where α is the amplitude-to-wavelength ratio. Faults that obey such self-similarity have a power spectral density (PSD) (Power and Tullis, 1991): where k = 2π/λ is the wavenumber (λ is the wavelength). Fault roughness is often charac- on a self-similar fault, the average resistance to sliding due to geometric complexity is given 52 by the roughness drag: where δ is slip magnitude and λ min is the minimum wavelength that is present in the fault 54 profile (other symbols are defined in Table 1). The spatial extent of the slip patch must be 55 much larger than λ min for this to be valid.
where the meaning of δ versus δ is discussed later. The matrices of influence coefficients 84 are compressed using the H-matrix approach of Bradley and Segall (2011). The frictional 85 interface is governed by rate-and-state friction and aging law, respectively: where V and θ represent the slip speed and state at the center of each dislocation respec-87 tively.

88
An infinite planar fault with the same frictional properties will oscillate around V 0 as 89 long as the long term average of the elastic stress transfer is τ = 0. This is reasonable; 90 otherwise, the long term average velocity of the fault would be changing, which can only 91 occur if the loading is changed. The problem is more complicated for a non-planar and/or 92    Table 1) are set such that the crack half-length, which 109 marks the transition from nucleation to a dynamic instability, is constant L ∞ ≈ 30 m and 110 is therefore also well resolved. The fault profile was generated with λ min ranging from L ∞ /3 111 to 10 · L ∞ , but in all cases with the same random seed such that the Fourier decomposition 112 at larger wavelengths is identical in both magnitude and phase. The n + 1 time-step of the simulations starts by computing shear and normal stress using 115 the slip of the previous time step: τ n+1 = G τ δ n and σ n+1 = G σ δ n , where τ n+1 = τ 0 + τ n+1 116 and σ n+1 = σ 0 + σ n+1 . Furthermore, a prediction of the n + 1 value of the state variable is 117 made: Now Eq. 5 can be rearranged to provide an approximation for the slip speed at time step 119 n + 1 given that the relevant fields are known at time step n.
At high slip speeds (∼ 1 cm/s) where inertial effects become important, Eq. 8 is not accurate 121 and generates numerical dispersion. At dislocation centers that satisfy (V n > 1 cm/s) I solve 122 for V n+1 by using Eq. 9.
Remarkably, Eq. 9 only needs to be solved approximately to suppress dispersion and attain

128
Next step updates the state-variable and slip using the following equations: The final step determines the new time-step where is adjusted such that stability and convergence is found, for this study it was set to 131 1/64. The problem is initialized such that τ = τ 0 , σ = σ 0 and θ = d c /V 0 (1 + N (0, 0.01)) at 132 all dislocation centers (See Table 1). The fault is thus approximately at steady state V = V 0 133 initially apart form small amplitude Gaussian white noise added to the initial state. The 134 noise is generated using the same random seed and is thus the same for all simulations. In 135 this study I explore the statistics of event sizes for simulations that spans multiple cycles.

136
Each simulation is run for 6.4 million time steps, but one simulation (α = 0) was extended to exceed 100 seismic events. The only difference between each simulation is λ min . I thus do  In this study, I investigate a limit where the fault is much larger than the nucleation 161 length, and chaotic behavior is expected (Barbot, 2019   From the slip profiles above, we observe that the initial rupture always propagates the 171 whole length of the fault. However, later events tend to be partial ruptures except when 172 λ min is large (Figure 4). Initially, the shear and normal stresses are selected to be spa-  We may now use details of the rate-and-state friction law to estimate the maximum slip  with the same dependence on λ min /α 2 as Eq. 14. However, their formulation included an 227 unknown fitting coefficient, whereas no fitting is done here.   (Figure 7), which well characterizes the fall-off with increased moment. It generally 269 appears λ min does not control the fall-off, but as has been previously noted, the truncation 270 of the distribution is changed by λ min . It is notable that even for the no-roughness limit, the 271 events follow the same power-law distribution. This is consistent with recent work (Cattania, 272 2019), which showed in simulations and theory that a planar fault that is sufficiently large 273 could manifest a power-law distribution of events (see further discussion in Section 4.1).

274
Some interesting differences are found in Figure 7, when comparing the cases of λ min L ∞ 275 to λ min = 10L ∞ and the no-roughness case. We notice that at low moment bins, the em-     damping is included in a quasi-dynamic simulation. For example, where τ el (t, δ) is the stress from elastic interaction or wave-mediated stress transfer on a 334 planar fault, and τ l (t) represents externally imposed loading of the fault (see Table 1 for 335 other definitions).

337
Roughness has an important influence on both individual ruptures and frequency and 338 magnitude characteristics of events. Events start as crack-like ruptures, but due to roughness 339 drag, they transition to pulse-like ruptures at a characteristic length-scale, L c , determined 340 by fault roughness alone (Eq. 13). I suggest that slip on faults much larger than L c cannot 341 be approximated using a planar fault without, at least, including roughness drag (Eq. 16).

342
Pulses lock in approximately spatially fixed slip distance. The maximum slip, δ c , during 343 pulse-like rupture is set by roughness drag but also depends on the assumed friction law and 344 material properties. I conclude that fault roughness offers a plausible and general mechanism 345 for earthquakes to transition from cracks to pulses as they grow. I find that decreasing λ min ,