Scaling of Repeating Earthquakes at the Transition From Aseismic to Seismic Slip

Some observations of repeating earthquakes show an unusual, non‐self‐similar scaling between seismic moment and corner frequency, a source property related to rupture size. These observations have been mostly reported in regions at the transition from stable to unstable slip, in geothermal reservoirs and subduction zones. What controls the non self‐similarity of these ruptures and how this is linked to the frictional stability of the interface are still open questions. Here we develop seismic cycle simulations of a single unstable slipping patch to investigate the mechanisms underlying this behavior. We show that temporal changes of normal stress on a fault can produce ruptures that exhibit the observed anomalous scaling. Our results highlight the role of fault zone fluid pressure in modulating the effective normal stress and contributing to the sliding stability of the fault.

stress drop of an earthquake, one notable one would be a change in loading rate. However changing the loading rate produces mostly a change of l and a modest variation of M 0 which is not consistent with the reported observations (Chaves et al., 2020;Vidale et al., 1994) (See Figure S1 in Supporting Information S1). Changing Δσ by changing fluid pressure is another possible scenario, which we consider here. We further use a numerical modeling approach to reproduce quantitative aspects of the anomalous source scaling observations. In particular, we investigate how a pore pressure perturbation on a fault can modify the nucleation of micro-earthquakes and the proportion of aseismic slip. We analyze the variations of the source parameters of the simulated seismic events as a function of the effective normal stress perturbation. We find that the abnormal scaling relation can emerge from distinct pore pressure conditions and is quantitatively explained by the model.

A Rate-And-State Friction Model
In earthquake cycle models, the equations governing slip on a fault are often modeled assuming rate-and-state state friction (Dieterich, 1992). We follow that framework and describe the modeling equations in Text S1 in Supporting Information S1. Under this model assumption, the stability of slip is governed by the balance between the shear stress imposed on the fault and the frictional strength. In a 1D spring slider model, unstable slip occurs when the rigidity of the medium, k, is higher than a critical rigidity (Scholz, 1998) where a and b are rate-and-state friction parameters that control the response of the friction coefficient to a change of slip velocity and of fault state variable, respectively, D c is a characteristic slip distance for the evolution of the state variable, and ef f = − is the effective normal stress on the fault. No spontaneous instability occurs if (b − a) < 0, that is, if the fault is velocity-strengthening at a steady state. If ef f decreases then k* decreases and the fault is less prone to unstable sliding. The rigidity of a locked circular patch of radius L on a plane driven by a remote stress is where G is the shear modulus of the elastic medium (Ampuero & Rubin, 2008). Combining Equations (1) and (2) gives the minimum size of the circular patch for nucleation of an instability A decrease in the effective normal stress causes an increase of L c . Thus, increasing the pore-pressure on a fault may stabilize its seismogenic patches and turn them into aseismic patches. It may also lead to an intermediate behavior in which seismogenic patches become subdued seismic patches, producing earthquakes with smaller moments compared to their unperturbed pore-pressure state. Here we investigate these intermediate cases with the help of a numerical model.

Modeling the Fault
We consider a quasi-dynamic numerical model of a fault governed by rate-and-state friction (Luo & Ampuero, 2018;Luo et al., 2017), as described in Text S1 in Supporting Information S1. We investigate a simple model of 1D straight fault embedded in a 2D elastic medium, driven by a far field loading and slipping in the anti-plane direction. While there might be quantitative differences between quasi-dynamic and fully-dynamic models (Thomas et al., 2014), we consider that the dynamic stress changes carried by seismic waves will not modify the results obtained here because we focus on isolated asperities with simple geometry. Moreover, we are mostly interested in the variation of observed parameters rather than in their absolute values, such that our approximate representation of co-seismic stress transfer should not impact the observed relative trends. We consider an asperity described as a potentially unstable patch with a < b. The rate-and-state friction properties a, b, and D c are assumed constant and independent of pore-pressure. While such a dependency is possible, as reported by Scuderi and Collettini (2016), we hypothesize that its effect on the fault shear strength is small compared to the linear dependence of shear strength on pore pressure via the effective normal stress.
We adopt similar properties of the fault and of the nucleation patch as Chen and Lapusta (2009) except that we consider a 1D fault. The elastic medium has a shear modulus G = 30 GPa and the far field loading velocity is V l = 10 −9 m/s. The fault has a steady state friction coefficient μ ss = 0.6 at the sliding velocity of V ss = V l . We also run simulations with a higher shear modulus (G = 45 GPa), keeping all other parameters constant, and obtain a similar conclusion (See Text S2 and Figure S3 in Supporting Information S1). A velocity-weakening patch of length L lies in the middle of the fault with a = 0.015 and b = 0.019. Outside this patch, we consider a velocity strengthening zone with a = 0.019 and b = 0.015. The total modeled domain has a length of L x = 1,000 m and outside the domain, we impose stable sliding at the plate velocity. The critical slip distance, D c , is set to 160 μm. We do the simulations with the boundary element software QDYN (Luo et al., 2017). The fault is decomposed into N = 1024 elements. The element size, dx = L x /N, is between 3 and 16 times smaller than the process zone size L b = GD c /bσ eff for the range of values of σ eff we considered, ensuring that the numerical model has a proper spatial resolution. To avoid the influence of the initial conditions, we only considered the results after several earthquake cycles have occurred. Each simulation is performed under temporally constant and spatially uniform normal stress. We systematically study the influence of the normal stress on the sliding stability of seismic fault patches of different sizes L.

Results
We first consider a velocity weakening region of half-size R = L/2 = 100 m. We run several simulations by imposing various values of the effective normal stress (typically between 20 and 200 MPa) which are representative of the conditions of an earthquake in the shallow crust. We show in Supporting Information S1 (Text S3, Figure  S4 and Tables S1 and S2 in Supporting Information S1) that the alternative approach of applying a step in σ n at some time within a repeating earthquake cycle and looking at the properties of the next rupture that immediately follows this step change, leads to similar results as those presented here. Under constant normal stress and constant remote loading rate, the velocity-weakening patch may experience phases of slip acceleration (instabilities). We track the evolution of the maximum slip velocity in Figure 1. It remains close to the loading velocity most of the time, punctuated by transient increases in slip rate. These transients only occur for sufficiently large normal stress. Below this normal stress threshold, the simulated fault is stable and no instability developed. We define v th = 1 cm/s as the velocity threshold separating aseismic from seismic slip, as in previous studies (e.g., Chen & Lapusta, 2009). When the normal stress is too low the slip rate remains lower than v th ; we then consider that all the slip on the asperity, even during the slip event, is aseismic. We observe that the maximum slip speed increases with increasing normal stress and exceeds v th when σ n becomes larger than about 30 MPa. The results of Rubin and Ampuero (2005) suggest that, for an antiplane rupture, the minimum normal stress required to cause seismic slip on a velocity-weakening patch of half-size R is Given our parameters values, we have σ n = 36 MPa, close to the observed value. For higher values of normal stress, the maximum slip velocity remains around 1 m/s. Therefore, there exists a minimum normal stress value above which seismically detectable events exist (Figure 1). Such prediction is a well known behavior of the rate-and-state friction model, in which the slip behavior of an isolated asperity is controlled by the ratio between the nucleation length defined by Rubin and Ampuero (2005) and the patch length (Barbot, 2019;Rubin, 2008).
We now document how the seismic moment evolves when the normal stress changes during seismic cycles. We determine the starting time, t s , and ending time, t e , of each seismic event as the times when the maximum slip rate becomes larger than and lower than v th , respectively. For each event, we compute the distribution of co-seismic slip, D, as a function of the distance to the center of the asperity, x i , as where u(x i , t) and v(x i , t) are, respectively, the slip and the slip rate computed along the fault at position x i and at time t (see Figure 2). The position x i varies from −L x /2 to L x /2 in N steps dx. We convert this distribution of slip from an 1D fault to slip on a (2D) fault plane following the approach of (Lapusta & Rice, 2003;Rubin & Ampuero, 2005) which assumes a radial symmetry of the slip profile, centered at the middle of the asperity, such that the seismic potency, P 0 , is computed as Finally, the seismic moment is obtained from M 0 = GP 0 .
We report in Figure 3 the evolution of the seismic moment, M 0 as a function of the effective normal stress, σ n . From the aseismic/seismic transition up to the maximum achieved seismic moment, M 0 increases by two orders of magnitude. Thus the same asperity can produce earthquakes with various seismic moments depending on their normal stress. We also report on the same figure the rupture size of seismic events, r. We simply consider that = max( )| ( ) > ∈ [ ; ] . We observe that r grows monotonically with the effective normal stress. The rupture size is either smaller than the seismic patch (velocity weakening zone) or larger by penetrating into the velocity-strengthening zone. The values of r are distributed around the values of the half-size of the unstable patch, R, and typically vary over a range of ±20%. Both values of r and M 0 grow rapidly after the aseismic/seismic transition and then flatten as the normal stress is increased.
To document how this observation translates at the moment versus radius scaling of earthquakes, we run simulations with different asperity radii (R = L/2 = 30, 70, 100, and 300 m) and, for each simulation, we extract the seismic moment and rupture size of each event when varying the normal stress (from 20 to 110 MPa). For simulations performed with R = 300 m, we increased the total modeled domain to a length of 2,000 m and the number of elements to 2048. We show in Figure 4 the scaling of the moment versus rupture size. We observe that, globally, all results fall within the typical scaling 0 = 7 16 Δ 3 where Δσ is the stress drop of a circular crack model (Eshelby, 1957). The range of stress drops identified in Figure 4 varies over two orders of magnitude, between 1 and 100 MPa. The lowest stress drops are observed for the lowest values of normal stress and, as the normal stress increases, so does the stress drop, as previously observed by Kato (2012). As expected from Equation 4, for a given asperity size, R, a seismic rupture is only observed for a sufficiently high normal stress. Below this value, the rupture is entirely aseismic (based on our definition). We note that the smaller size earthquakes need only a slight decrease of normal stress to undergo a profound variation  of stress drop, while for larger ruptures the change of normal stress required to produce the same effect is more important. This suggests that the effect of normal stress is mainly visible in the smaller ruptures.
For a given seismic patch or asperity size, the scaling between moment and rupture size appears to depart from Equation 7, showing a sharp increase of moment with rupture size. We document this scaling of the moment with rupture size within a seismic patch (i.e., a repeating earthquake sequence corresponding to a set of events with a similar asperity size but various prescribed values of normal stress). To reveal this connection, we compute for each of the four seismic patch sizes, R, the normalized moment, M 0 /〈M 0 〉, and the normalized rupture size, r/〈r〉, normalizing by their mean values. Despite some scatter, the two orders of magnitude variation of the normalized moment are retrieved while the variation of the normalized rupture size is small ( Figure 5). These results clearly indicate that for repeating ruptures on a seismic patch with varying normal stress, a moment -size scaling emerges with an exponent much higher than the value of three expected from self-similarity. A linear least squares regression between the logarithm of the normalized moment and the logarithm of the normalized rupture size reveals that the two quantities are well related by a linear relation with a slope of 12 ( Figure 5). We acknowledge however that a power law model is of limited quality due to the very small range of normalized rupture size and one could also consider other forms of laws to fit these data. However, we assume this power-law expression as it is the one reported in studies on natural cases (Bostock et al., 2015;Cauchie et al., 2020;Farge et al., 2020;Harrington & Brodsky, 2009).

Discussion
Our results can be compared to observations from active faults. For example, Bostock et al. (2015) resolve a similar scaling as in Figure 5 linking the corner frequency to the moment of low frequency earthquakes in the Cascadia subduction zone. Assuming the corner frequency is inversely proportional to the rupture length, as in classical earthquake source models (Savage, 1972), the moment -size exponent for the low frequency earthquakes in Cascadia is around 10. Similarly, Farge et al. (2020) resolve such an anomalous scaling for low frequency earthquakes in the Mexican subduction zone with an exponent between 8 and 19 depending on the method used. In the Soulz-sous-Forêts geothermal area, Cauchie et al. (2020) show that the scaling relation within repeating earthquake families exhibits a similar scaling with an exponent close to 20 (see Figure 4 for an example). For earthquakes on the San-Andreas fault at Parkfield, Harrington and Brodsky (2009) get an exponent of 17.
All these observations show that for a single family the seismic moment varies at most over two orders of magnitude. Our simulation results are in good agreement with these reported observations and thus constitute a possible explanation for this observed anomalous scaling. We note that our observations are not necessarily in contradiction with the self-similar scaling of LFE as reported for example, by Supino et al. (2020). Indeed, in their study, Supino et al. (2020) investigate a complete LFE population and not specifically repeating families. As we observed in Figure 4, an LFE population comprising seismic patches of variable size is bounded by the self-similar scaling and the unusual scaling is only observed within repeating earthquake families.
The variation of stress drop observed in natural earthquakes and the associated non self-similar scaling could also arise from other considerations. Indeed, as noted by Kaneko and Shearer (2015), rupture directivity or effects of complex geometry compared to the simplistic circular rupture model could also give rise to a variation of the earthquake moment with an almost constant corner frequency. It is also possible, as demonstrated by numerical simulations, that a fault with heterogeneous strength can lead to seismic ruptures on the same fault patch displaying variable moments but nearly constant apparent rupture size (Lin & Lapusta, 2018). Considering the diversity of these effects and their randomness, it seems quite unlikely that they all favorably contribute to producing the anomalous scaling observed across various tectonic settings. We would rather expect that such complexities produce scattering around an average value without any systematic trend. The anomalous scaling inferred from earthquake observations has been also explained by invoking a totally different mechanism involving elastic collisions between fault gouge particles (Tsai & Hirth, 2020).
Our model is limited in several aspects. First, we simulated a 1D fault in a 2D medium; some deviations can arise compared to a rupture on a 2D fault. Our approach also requires some assumptions on the slip distribution in order to compute the potency of each rupture. However, Li et al. (2022) recently showed that numerous outcomes of 2D and 3D numerical earthquake cycle models, such as stress drop, are comparable. This supports the validity of the presented results for higher dimensions. Second, we did not incorporate the most recent advances in friction models and fault weakening mechanisms, for instance, thermal pressurization or flash weakening processes (Acosta et al., 2018;Lambert et al., 2021). We thus acknowledge that the additional physics contained in these models can give rise to results different than the ones reported here. Furthermore, we did not perform a systematic parametric study, varying the a, b and D c parameters to test their influence on the resolved scaling. The set of parameters tested in this study has been previously considered in other simulations whose results were found stable while perturbing these values (Chen & Lapusta, 2009). Furthermore, here we only considered an isolated asperity, not interacting with any other asperities. However, faults generally contain several seismic patches that can interact and trigger each other. It remains to be investigated how these interactions can influence the properties derived in this study and impact the observed scaling, Finally, we stress that the normal stress on an asperity is not-necessarily uniform (Schmittbuhl et al., 2006). This can lead to some important effects as the change of fluid pressure on the fault can lead to a change of the contact area of the asperity and thus redistribute stress locally and modify the asperity. Here we preferred to keep a rather simple model which, despite its limitations, offers a straightforward mechanism for interpreting the variation of the moment despite similar rupture sizes observed for numerous repeating earthquake sequences worldwide. Our model proposes that these characteristics can be well understood within the framework of a frictional fault with varying average normal stress. This model requires normal stress fluctuations at the location of the asperity. The most direct explanation for such fluctuations involves the presence of fluid pressure and its variation. The existence of fluid at the location of the seismogenic patch is well understood for geothermal reservoirs but is still debated as a necessary component for the generation of low frequency earthquakes in subduction zones (Saffer & Tobin, 2011). The change of fluid pressure on the fault can then arise from a variation of the fluid pressure directly from the source region, notably in geothermal systems, or because slip on a nearby portion of the fault modifies the fluid flow and locally enhances fluid pressure (e.g., Shapiro et al., 2018).
In conclusion, our study highlights that changes in the effective normal stress can cause a significant variation of the seismic moment on a repeating earthquake sequence. This variation of the moment leads to a variation in stress drop of at most two orders of magnitude. This important fluctuation of the stress drop is observed at the transition between the seismic and aseismic slip, such that the repeating earthquake sequence exhibits a peculiar scaling behavior that can be used as an indicator of proximity to the frictional regime change.

Data Availability Statement
No data were used in this study. Version 2.3 of the software QDYN used for modeling the seismic cyle is preserved at https://doi.org/10.5281/zenodo.322459, available via GNU General Public License v3.0 only and developed openly at https://github.com/ydluo/qdyn.