Title : Viscous fault creep controls the stress-dependence of modelled earthquake statistics

The ability to estimate the likelihood of particular earthquake magnitudes occurring in a given region is critical for seismic hazard assessment. Earthquake size and recurrence statistics have been empirically linked to stress state, however there is ongoing debate as to which fault-zone processes are responsible for this link. We numerically model combined viscous creep and frictional sliding of a fault-zone, where applied shear stress controls the interplay between these mechanisms. This model reproduces the stress-dependent earthquake magnitude distribution observed in nature. At low stress, many fault segments creep and impede ruptures, limiting earthquake sizes. At high stress, more segments are close to frictional failure and large earthquakes are more frequent. Contrasts in earthquake statistics between regions, with depth and through time, may be explained by stress variation, which could be used in the future to further constrain probabilistic models of regional seismicity.


41
The cornerstone of modern seismic hazard assessment is the ability to model the statistical heterogeneous combinations of localized frictional sliding and distributed viscous creep ((27-29), 98 Fig. 1a). At the 100 km scale, some fault segments appear to accommodate permanent 99 deformation (at present day) without accumulation of elastic strain in geodetic inversions (26, 30, 100 31), and are interpreted as weak, creeping zones. Subduction zone megathrusts with a large 101 proportion of creeping material commonly occur where the hanging wall is experiencing 102 extension (32), implying low tectonic differential stress. Weak phyllosilicates have been found in 103 a creeping section of the San Andreas fault (33). Creeping regions may also act as 'barriers' to 104 earthquakes, promoting their arrest, which is also consistent with relaxed stress (30). 'Locked' 105 regions that do not host inter-seismic creep instead may act as 'asperities', which fully rupture in 106 potentially large earthquakes (34). As loading stress controls whether creep or frictional failure 107 occur, the distribution of creeping and locked regions may be modified by changing tectonic 108 stress conditions (9, 35).
the seismogenic zone inferred in nature (28, 37) and from microphysical models (38). Each fault element is composed of frictional and viscous mechanisms in series, such that the weakest one 150 dominates deformation. At low stress, viscous deformation is dominant (Fig. 1b). A higher stress 151 is required for relatively slow frictional sliding to dominate, as is required for earthquake 152 nucleation to occur. At steady-state seismic slip velocities, frictional sliding occurs at reduced 153 stress due to dynamic weakening, however the fault must strengthen again for subsequent 154 earthquakes to nucleate.

156
Shear-zone thickness is the only parameter that is varied, which acts to control strain-rate ̇ 157 (for a given slip rate , ̇= / ) and thus viscous strength ( =). A large lowers the 158 strain-rate, making creep more efficient and lowering the stress that the fault-zone can deform at 159 (Fig. 1b,c). While we change to control fault background stress, this may also represent the  At relatively low stress (reproduced using large W), such as for Fig. 3a, earthquakes are generally 172 limited to isolated regions of high viscosity and are consequently restricted to low magnitudes 173 ( ! < 5.7 and on average ! ∼ 4). Defining an asperity as a friction-dominated area, where 174 earthquakes nucleate, 'effective asperities' are predicted as areas where the stress required for 175 creep to accommodate the loading slip rate is greater than the steady-state frictional strength, over 176 lengths larger than the earthquake nucleation length-scale (the black stripes at the bottom of each 177 panel in Fig. 3). At low stress, earthquakes are predominately limited to each of these effective 178 asperities, rarely propagating into adjacent creeping regions or spanning multiple asperities.

180
With decreasing (and correspondingly increasing stress), the effective asperity sizes increase 181 (Fig. 3b-c). The larger effective asperities correspond to fault regions hosting larger earthquakes, 182 that also occasionally span multiple asperities. Small earthquakes also persist, both hosted on 183 small asperities and occurring as partial ruptures of larger asperities, nucleating on patches with 184 particularly high viscosity. When "#$ = 47 MPa ( = 50 m), earthquakes with ! < 6.9 185 occur, due to ruptures that propagate over large effective asperities as well as small intervening 186 regions (~ 1 km wide) of low viscosity that are otherwise dominated by inter-seismic creep (e.g.

187
at 15-20 km in Fig. 3c). Large events occur less frequently than small events and with greater displacement, as occurs in natural scaling relationships. The largest events in the reference model-189 set are ! < 7.4 and are limited by the fault length. We ran the reference model-set three times 190 with different randomized viscosity distributions. Both maximum stress "#$ and maximum 191 magnitude ! can vary for a given , though we will show that the characterization of 192 seismogenic behavior in terms of "#$ is robust. The model statistics are described in terms of "#$ (a model output), the range of which is 207 produced by varying , as the seismogenic behavior is later shown to be generalized in terms of 208 fault stress for all models (which also holds for time-averaged stress, but not , Fig. S2) and 209 long-term fault stress is the quantity linked to tectonics. We firstly describe the maximum ! and 210 seismic coupling , which is calculated as the ratio of the total accumulated seismic slip to the 211 total loading displacement, over the analyzed model period. Both maximum ! and show a 212 clear trend of increase with stress ( Fig. 4a-b). This increase is gradual over the range "#$ = 20 -213 40 MPa and more rapid at greater stress. We define an aseismic model as being characterized by 214 < 0.3 and maximum ! < 6. While there is no agreed definition of aseismic behavior, as 215 seismicity is pervasive in the Earth's crust, these values are typical of less seismogenic subduction 216 zones (2, 31). Models with "#$ < 40 MPa are then relatively aseismic, while becoming 217 increasingly seismogenic at higher stress.

219
We calculate b-values for the reference model-set by combining the randomized model 220 realizations into single, more statistically complete, catalogues for each modelled W. The 221 recurrence times of seismic events are reasonably approximated by the G-R relationship (Fig. 5).

222
The b-value increases with increasing and correspondingly decreasing "#$ , reflecting the 223 decreased likelihood of large events with decreasing stress. The b-value becomes relatively 224 constant at ~1.5 for "#$ ≤ 34 MPa ( ≥ 300 m), within uncertainty (Fig. 4c), which we 225 characterize as aseismic, in agreement with the characterization based on maximum ! and . At

234
The models can be characterized more generally in terms of the non-dimensional pre-stress ratio  We tested the sensitivity of our results to the imposed viscosity probability distribution, using   We also quantify how visco-frictional interplay influences seismogenic behavior by plotting the 267 statistical data as a function of the proportion of the fault-zone that is locked ( Fig. 4e-g), , The b-value variation with stress can be related to the growth of effective asperities with 285 increasing stress. While the imposed patch sizes follow a power-law, different combinations of 286 these patches combine into the effective asperities for varying "#$ . The size distribution of the 287 effective asperities still follows a power-law for each , with the power-law exponent decreasing 288 with increasing stress (Fig. S4). The effective asperity size distributions can be converted to 289 predicted event catalogues by assuming that each earthquake is confined to an isolated effective 290 asperity, as is generally the case at low "#$ (Fig. 3), that all events have a constant stress drop 291 and that the stressing rate of an asperity is inversely proportional to its size (such that small 292 asperities host events more regularly). The b-values for these predicted events agree with the 293 modelled events, indicating that the effective asperity size distributions play a role in the stress-294 dependence of the b-value (Fig. S4).

296
We use this method of predicting the b-value from the effective asperity distribution to test the  To further examine this rupture behavior we use an analytic energy balance calculation, which 325 predicts that rupture arrest occurs when there is insufficient available stress drop (i.e. ̅ % ) at the 326 rupture front to drive further rupture propagation (15). This calculation broadly reproduces the 327 dependence of rupture penetration on both matrix viscosity and patch width, for relatively narrow 328 patch widths, and dependence primarily on viscosity for the widest patches. However, ruptures arrest more rapidly in the numerical models than estimated by the energy balance calculation.

330
This is because ̅ % tends to be lower than expected for steady-state creep, as stress shadows form 331 adjacent to locked patches. respectively. In the absence of depth-dependent viscosity, the lower and higher stress models 362 have "#$ = 43 MPa and 47 MPa respectively, near the lower and higher ends of the pre-stress 363 range required to reproduce a variety of seismogenic behavior ( ̅ % = 0.26 and 0.43; Fig. S3).

365
The upper half (above 27.5 km) of the lower stress model (Fig. 7a) hosts earthquakes ranging 366 from ! ~ 3 to 7, with a b-value of 1.1 and = 0.4 (Fig. S5). The viscosity is lowered with

378
The higher stress subduction thrust model is shown in Fig. 7b. The fault shear stress is more  The modeling presented here has analyzed the underlying mechanism responsible for variation of rupture collectively at high stress (9), corresponding to our modelled effective asperities. In this case, the distribution of earthquake sizes and nucleation sites depends on a combination of 438 inherited properties at small scales and tectonic stress at larger scales. There is subsequently 439 uncertainty in using the extent of asperities to constrain the maximum ! , as they may change 440 effective size or link together with changing stress, particularly following changes to the regional 441 fault stress state throughout the earthquake cycle.

443
These models indicate deterministic relationships between b-value, , maximum ! and stress, 444 where the greatest variation in seismogenic behavior occurs over a range of 40 to 50 MPa shear 445 stress, or ̅ % = 0.1 to 0.5. This stress range is smaller than previous ~100 MPa order estimates 446 based on regional tectonic stress (8, 46), reflecting instead the lower shear stress that plate 447 interfaces locally operate at (47). Geodynamic estimates indicate a ~20% variation in plate  MPa), which is a sufficiently small stress contrast for seismogenic behavior to vary spatially 485 between tectonic settings. These model applications are preliminary, but highlight the potential to apply earthquake cycle models that incorporate inter-seismic creep in understanding regional Viscous deformation is assumed to follow Newtonian viscous creep, for shear zone thickness W 507 and viscosity : Frictional deformation is assumed to follow the regularized form of rate and state friction (58), 509 which allows the frictional slip rate to vanish (when creep dominates) and can be written as: Where is the effective normal stress, is the direct effect parameter, is the evolution effect 513 parameter, : is a characteristic slip distance which governs the evolution of the state parameter 514 , and % is a reference velocity at the corresponding shear stress % . Finally, we take the aging 515 law to describe the state evolution: 517 518 QDYN adopts the quasi-dynamic approximation by Rice (1993), with the associated force In which ;< is a stiffness matrix that relates the change in shear stress on the -th fault element 523 ( ; ) to slip on the -th fault element ( < ) (the Einstein summation convention is adopted here). The 524 last term on the right-hand side represents the stress change due to seismic wave radiation normal to the fault plane, and comprises the shear modulus of the medium and the shear wave speed ( , 526 which are both uniform and constant in the simulations.

528
In order to connect the force balance with the point-wise fault rheology, we expand the equation 529 for slip rate into its partial derivatives: 532 533 Substitution into the force balance and rewriting then gives: The partial derivatives in this relation are readily obtained from the rheological model, and read: Finally, we note that for parallel operation of frictional slip and ductile creep, it is necessary to 548 solve for and and compute = ( , ), rather than solving for and and computing = 549 ( , ), which is a common choice in the modelling community. The reason for this is that it is 550 non-trivial to satisfy the hard constraint = -= 1 while solving for in the latter approach. Each patch has a uniform , sampled from a log-uniform distribution, for ";A = 10 6C Pa s, 561 "#$ = 10 )% Pa s and a uniformly random variable n (varying between 0 and 1): We assume the frictional parameters = 9e-3 and = 2e-2, characteristic slip distance : = 1e-564 2 m, reference friction coefficient % = 0.6 and reference velocity % = 1e-9 m/s. These The maximum average shear stress "#$ is measured over the final 1000 years of the model run, 586 in order to avoid the viscous relaxation that is still occurring at the start of some aseismic models 587 (Fig. 2), though there is no discernible change in seismogenic behavior during this relaxation 588 period. "#$ is chosen as the primary metric to characterize models as it corresponds to the pre- Four additional model-sets (Table 1) were run in order to test the sensitivity of our modelled 626 earthquake statistics to the probability distribution used to generate the viscosity distribution.

627
Model-sets B and C were identical to the reference model-set, however assuming a lower and  The trends of decreasing b-value and increasing and maximum ! with increasing stress 644 generally hold for these alternative viscosity distributions (Fig. 4). A lower viscosity contrast 645 (model-set B) tends to result in the possibility of some relatively aseismic models at high stress 646 (Fig. 4b), while a higher viscosity contrast (model-set C) makes little difference from the  Fig. 4c).

653 654
Single asperity rupture test models

656
A suite of simplified models were run to constrain how far a rupture can propagate through a 657 creeping region, to provide further context for the primary model-set. A patch of high viscosity 658 (10 20 Pa s), with width # , was surrounded by homogenous material with viscosity " , assuming 659 =10 m (example model shown in Fig. S7). Both " and # were varied over 35 models, in 660 order to explore their influence on rupture propagation distance into the creeping region (Fig. 6).

662
For most models, the rupture propagation is a small fraction of # , except for a subset of models

666
(mostly still at the km scale) that decreases with increasing ̅ % . This critical ̅ % corresponds to 667 when stress is relaxed to 10 MPa below ( .

669
A simple energy balance model was calculated to further understand why rupture length is so 670 limited at low ̅ % , even for large patches. A rupture has a stress intensity factor that evolves 671 with rupture length and available stress drop (Fig. S8). Rupture arrest will occur if < : . The The matrix is assumed to be creeping prior to earthquake nucleation, such that the maximum The analytical solution is shown in Fig. 6 and predicts a plateau in rupture length with increasing