The influence of faulting style on the size-distribution of global earthquakes

Earth and Planetary Science Letters 527 (2019) 115791

Contents lists available at ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

The influence of faulting style on the size-distribution of global earthquakes A. Petruccelli a,b,∗ , D. Schorlemmer c , T. Tormann b,1 , A.P. Rinaldi b , S. Wiemer b , P. Gasperini a , G. Vannucci d a

Dipartimento di Fisica e Astronomia, University of Bologna, Italy Swiss Seismological Service, ETH Zurich, Switzerland c Helmholtz Centre Potsdam, GFZ Potsdam, Germany d Istituto Nazionale di Geofisica e Vulcanologia, Bologna, Italy b

a r t i c l e

i n f o

Article history: Received 4 November 2018 Received in revised form 22 August 2019 Accepted 24 August 2019 Available online 9 September 2019 Editor: M. Ishii Keywords: statistical seismology earthquake size-distribution faulting styles Anderson’s theory of faulting Mohr-Coulomb failure criterion

a b s t r a c t We derive a unifying formulation, reliable at all scales, linking Anderson’s faulting theory with the earthquake size-distribution, whose exponent is known as the b-value. Anderson’s theory, introduced in 1905, related fault orientation to stress conditions. Independently, laboratory measurements on acoustic emissions have established that the applied differential stress controls their b-value. Our global survey revealed that observed spatial variations of b are controlled by different stress regimes, generally being lower in compressional (subduction trenches and continental collisional systems) and higher in extensional regimes (oceanic ridges). This confirmed previous observations that the b-value depends on the rake angle of focal mechanisms. Using a new plunge/dip-angles-based b-value analysis, we also identified further systematic influences of faulting geometry: steep normal faults (also typical of the oldest subduction zones) experience the highest proportion of smaller events, while low-angle thrust faults (typical of youngest subduction zones) undergo proportionally larger, more hazardous, events, differently from what would be expected by only allowing for rake-angle dependency. To date, however, no physical model has ever been proposed to explain how earthquakes size-distribution, differential stress and faulting styles relate to each other. Here, we propose and analytically derive a unifying formulation for describing how fault orientation and differential stresses determine b-value. Our formulation confirms that b-values decay linearly with increasing differential stress, but it also predicts a different dipdependent modulation according to the tectonic environment, opening up new ways of assessing a region’s seismic hazard. © 2019 Elsevier B.V. All rights reserved.

1. Introduction One of the ongoing debates in seismology concerns the understanding of the physical implications of earthquake size-distribution, its potential variations on different scales and its use in earthquake hazard assessment and forecasting. Empirically, the behaviour of the number of detected earthquakes N ( M ) at a magni-

*

Corresponding author at: Swiss Seismological Service SED, ETH Zurich, Switzerland. E-mail addresses: [email protected] (A. Petruccelli), [email protected] (D. Schorlemmer), [email protected] (T. Tormann), [email protected] (A.P. Rinaldi), [email protected] (S. Wiemer), [email protected] (P. Gasperini), [email protected] (G. Vannucci). 1 Now at PartnerRe, Zurich, Switzerland. https://doi.org/10.1016/j.epsl.2019.115791 0012-821X/© 2019 Elsevier B.V. All rights reserved.

tude M is generally well expressed by the Gutenberg-Richter (GR2 ) (Gutenberg and Richter, 1944) relation:

log N ( M ) = a − bM

(1)

In this GR relation, the b-value (commonly a value of around one) quantifies the relative proportion of larger to smaller earthquakes: the lower the b-value, the higher the relative frequency of large magnitude events and, in general, the higher the seismic hazard. The large number of publications on GR b-values reflects the importance of a proper understanding of this parameter.

2 Abbreviations used throughout the paper: b: b-value; CMT: centroid moment tensor; FM: focal mechanism; FS: faulting style; GR: Gutenberg and Richter law; NR: normal faulting mechanisms; S2005: Schorlemmer et al. (2005); SHA: seismic hazard assessment; SS: strike-slip mechanisms; TH: thrust faulting mechanisms.

2

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

Fig. 1. Dip-slip faulting schemes, principal axes orientations and ternary diagram. If principal stress (σ1 or σ3 ) is always assumed to occur along the vertical z and intermediate stress (σ2 ) lies in one of the horizontal directions y, the dip-slip faulting mechanisms (NR – green – and TH – blue) are easily inferred and the tectonic stress σ H V acting on the fault is the differential stress σ = σ1 − σ3 . Vectors: σn = normal stress, τ = shear stress, n = fault plane normal, d = slip vector. Angles: β = fault dip angle, ϕ = fault strike angle (set to 0◦ in our model), λ = rake angle, θ = complements β , δ = plunge angles, and P, T, B = the main deformation axes. (For interpretation of the colours in the figure(s), the reader is referred to the web version of this article.)

Stress distribution in the Earth’ crust is a critical parameter for understanding earthquake nucleation, but the sparsity of insitu measurements means that stressing conditions on faults are probably the key unknown for advancing earthquake forecasting. Laboratory measurements repeatedly conducted since the 1960s (Amitrano, 2003; Goebel et al., 2013; Scholz, 1968) have established that the differential stress (i.e. the difference between the maximum, σ1 , and minimum, σ3 , principal stress in the stress tensor) applied to a rock sample determines the b-value of acoustic emissions: the higher the applied differential stress, the lower the b-value. On the other hand, other factors such as confining and pore-fluid pressure (which may entail effective stress) (Sammonds et al., 1992), and material heterogeneity (Mogi, 1962), which itself determines fracture toughness, might play a role in causing b-values to vary. A wide range of observations is consistent with the laboratory results: asperities (locked patches of faults) are, by definition, areas of high differential stress, and many observational studies have shown them to exhibit low b-values (Ghosh et al., 2008; Schorlemmer and Wiemer, 2005; Tormann et al., 2015, 2014, 2012). By contrast, weak zones such as creeping sections of faults, mid-ocean ridges and volcanic regions typically have high b-values (Kagan, 1997; Okal and Romanowicz, 1994; Roberts et al., 2015; Schorlemmer et al., 2004; Tormann et al., 2014). In parallel, another interpretation, based on the theory of critical phase transition, explains the apparent decrease in b-value as being attributable to the truncation of the power law by an exponential, whose characteristic size increases as the rupture approaches (Amitrano, 2012). Anderson’s theory of faulting (Anderson, 1905) defines three possible styles, depending on the orientation of the three principal stresses. For a given faulting style, the differential stress required for reactivation is strictly linked to frictional properties and fault orientation. In the case of reactivation, for optimally oriented faults, the expected dip varies with the friction coefficient μ. For example, at a coefficient of 0.6, the dip angle is approx. 30◦ , as opposed to around 60◦ for compressive and extensive faults respectively, whereas dipping transcurrent faults are characterised by roughly 90◦ and intermediate stress conditions. Indeed, a systematic dependence of b-values on faulting style –

via the rake angle λ of FM – for global and regional datasets has been proposed (Gulia and Wiemer, 2010; Petruccelli et al., 2018; Schorlemmer et al., 2005; Yang et al., 2012). Extensional regimes dominated by normal faulting (NR, λ = −90◦ ) are characterised by the highest b-values (approx. 1.1–1.2), whereas regimes dominated by thrust faulting (TH, λ = 90◦ ) are subject to the highest differential stress and have low b-values (around 0.7-0.9), whereas strike-slip faults (SS, λ = 0, ±180◦ ) come in the middle (roughly 1). For subduction zones, b-values reportedly depend on local tectonic properties, such as the age of the subducting plate or slab buoyancy (Nishikawa and Ide, 2014). So far, however, there has been no quantitative comparison and model that unifies Anderson’s faulting theory, differential-stress dependency and earthquake size-distribution. To develop such a framework, we conducted a global survey on earthquake sizedistribution and tested whether the spatial patterns of b-values reflect major global seismotectonic structures, as predicted by Anderson’s theory. We then developed a new analytic approach linking b-values to faulting geometry, using a plunge-based ternary diagram (Frohlich, 2001; Frohlich and Apperson, 1992). Finally, we combined those analyses of b-values with fault modelling to derive a relationship linking b-values to the differential stress applied when modelling faults. 2. An empirical model: theory and calculation We developed a model to explain the observed variations of b-values with plunge angle δ T . Ranging from 0◦ to 90◦ , δ angles (Fig. 1) correspond to the dip (with respect to the horizontal direction) of the P , B, and T axes (i.e. the moment tensor eigenvectors), with sizes (eigenvalues) ranging from lowest to highest respectively. These axes form an orthogonal system and, as a first approximation, permit a reasonable estimation of principal stress directions s1 , s2 , s3 (Célérier, 2010). In the ternary diagram (or triangle), each corner corresponds to a ‘pure’ tectonic style (δ = 90◦ ), and each point is unequivocally defined by the three δ angles (Frohlich, 2001; Frohlich and Apperson, 1992; Serpelloni et al., 2007). The original classification (Frohlich and Apperson, 1992) defines SS and NR mechanisms as those where the

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

values δ B and δ P are greater than 60◦ , and TH as having a value δT greater than 50◦ . Mechanisms that satisfy none of these criteria are classified as ‘odd’ and occupy the central part of the ternary diagram. Anderson’s faulting theory (Anderson, 1905) describes the stress conditions needed to activate a fault in a given orientation, given the frictional properties of the rock. For pure dip-slip faults, the difference between maximum and minimum tectonic stress (σhV = σh − σ V ), corresponds in absolute values to the differential stress (σ = σ1 − σ3 ), i.e. difference between maximum and minimum principal stress. Applying the Mohr-Coulomb failure criterion (Coulomb, 1776), given a friction coefficient μ, the differential stress required to reactivate a dip-slip fault with dip β can be expressed as a function of plunge angle δ T as (see Appendix A for details and Fig. 1):

σ (δT | μ) = =

2μσ  lith

± sin 2β(δT ) − μ[1 − cos 2β(δT )] 2μρ gzref [1 − ΛTH ] NR

± sin 2β(δT ) − μ[1 − cos 2β(δT )]

(2)

where the plus sign (+) applies to TH events (δ T → 90◦ ) and the  = σ (1 − minus sign (−) applies to NR events (δ T → 0◦ ). σlith lith   Λ) = ρ gzref [1 − Λ TH ] is the ‘effective’ lithostatic stress (corrected NR

by a faulting-style-dependent factor Λ (ranging from 0 to 1) that considers the reduction in effective stress due to pore pressure), ρ is the assumed crustal density of 2,700 kg/m3 , g is gravity and zref is a reference depth depending on tectonic style). We assumed a range of 0.6-0.8 for the friction coefficient and a cohesionless fault because this range is the most commonly encountered in the literature, though recent analyses infer that a fault zone may actually present much lower frictional coefficients (Collettini et al., 2019). Our model assumed two different lithostatic effective stress states, i.e. two different reference depths, one deeper TH (30 km) and one shallower NR (15 km). This choice implied two different values for factor Λ, i.e. with pore pressure varying according to the faulting style. For NR, we considered a hydrostatic stress gradient (Λ = 0.3), while for TH we set an ‘over-pressurised’ Λ equal to 0.4 (see Appendix). In fact, portions of upper crust are mostly hydrostatic (Townend and Zoback, 2000), while higher fluid overpressures are common in compressional regimes, allowing the reactivation of fault at low shear stress (Sibson, 2004, 2014; Copley, 2018). Although the stress levels in the Earth’s crust acting on faults are generally not known, the differential stresses inferred from equation (2) can be compared with observed b-values. As a general trend, b-value and inferred differential stress are inversely proportional, as already proposed (Scholz, 2015), the parameters being a reference value br (1.23) and a global stress gradient k (0.0012 MPa−1 ) computed from a linear regression procedure carried out on depth-dependent b-value data (Spada et al., 2013) from different parts of the world. The general, linear relation was obtained by assuming a global friction coefficient of 0.75 and different stress gradients according to the local dominant tectonic regime: for compressional regions dominated by TH a vertical gradient of differential stress of 45 MPa/km, for extensional regions dominated by NR a gradient of 11.25 MPa/km, and for SS regions a gradient of 20 MPa/km. Actually, bearing in mind the region under consideration (Scholz, 2002), both the stress gradient and the reference value br could be a constant kFS , depending on the faulting style (FS), whereby:

b(σ ) = br ,FS − kFS σ

(3)

3

This suggested that b-values of global earthquakes are indeed determined by the respective differential stress levels set by the faulting regime. We could then generalise the relationship proposed for dip-slip regimes (Scholz, 2015), combining equations (2) and (3) as:

b(δ T | bFS , kFS , μ) = bFS ∓ kFS

[2μρ gzref [1 − ΛTH ]] NR

sin 2β(δ T ) ∓ μ[1 − cos 2β(δ T )]

(4)

where the upper sign refers to TH and the lower sign to NR fault zones, while bFS and kFS are faulting-style-dependent constants of the system. The unifying formulation (4) combines the empirically-observed inverse proportionality between b with differential stress σ , which in turn is modulated by a frictional term that includes variations in dip angle β . Equation (4) is only valid in a dip-slip configuration (NR-TH, λ = ±90◦ ), where one of the principal axes is assumed to be vertical (σ1 or σ3 respectively) and the intermediate one (σ2 ) is always horizontal (see Fig. 1). 3. Global survey of b-values (and M c ) If differential stresses impact b-values not just in the laboratory, but also in nature, and do so in accordance with equation (4), we must expect b-values to vary systematically between different tectonic provinces around the globe. To test this hypothesis, we conducted a comprehensive global survey of b-values, based on the Global CMT catalogue (Dziewonski et al., 1981; Ekström et al., 2012) covering the period from 1 January 1980 to 30 September 2016 (Fig. 2). We limited our analysis to earthquakes with a hypocentral depth of 0–50 km (Schorlemmer et al., 2005) and used moment magnitudes MW binned to M = 0.1. The GCMT catalogue provides strike, dip, and rake values for both nodal planes (Ekström et al., 2012), referring to the first plane as plane 1 and to the second as plane 2. Before estimating local b-values, we had to find a way of computing the local magnitude of completeness, M c , from the closest seismicity (Fig. 2). M c is the value above which the GR relation holds, defining the complete part of the dataset. M c is a critical parameter for determining b-values, because if it is underestimated, b-values will be systematically underestimated, too (e.g., Woessner and Wiemer, 2005). Since catalogue completeness depends on the distribution of seismic stations, which changes over time, spatiotemporal variations in completeness are expected in every catalogue. Overall, completeness is expected to improve over time as network density increases. Separate local assessments of M c are thus required for different periods. For each earthquake, we estimated the local M c (Fig. 2A) from the closest seismicity and plotted this value at the epicentre. Spatial b (and Mc ) distributions are often computed by selecting events within sampling volumes on regular grids (Wiemer and Wyss, 2000). However, in view of the non-uniformity of global earthquake density, we adopted a slightly different approach, taking the epicentres of all earthquakes in the top 50 km as nodes. Around each node, we selected all earthquakes within cylinders of different radii (ranging from 200 km to 1,500 km, with a step width of 100 km) down to a depth of 50 km (Fig. 2B). To ensure reliable parameter estimates, we only considered cylinders containing at least 200 events for M c estimations and more than 100 events above the completeness magnitude for subsequently estimating the b-values. To arrive at the best possible spatial resolution, we used the smallest of the above cylinders exceeding the required event numbers and estimated M c using the ‘maximumcurvature’ method (adding 0.2 magnitude units to be conservative, Wiemer and Wyss, 2002; Woessner and Wiemer, 2005). When assessing standard deviation, computations were bootstrapped 100 times.

4

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

Fig. 2. A global map of M c (A) and of selection radius R (B) for the GCMT catalogue (1980–2016, depth = 0 to 50 km). Dotted black line: global plate boundary model PB2002 (Bird, 2003); solid black line: subduction zones according to the SLAB 1.0 model (Hayes et al., 2012). A, B, C and D are sample points (see Fig. 3).

Steady improvement of the global seismographic network has raised the GCMT catalogue’s completeness level from values of 5.4 and 5.5 (between 1976 and 2004) to 5.0 (after 2010) (Ekström et al., 2012). However, the M c map (Fig. 2A) for the GCMT catalogue over the entire period since 1980 shows that most M c values fall between 5.4 and 5.5, with only minor deviations in a few particular areas: along most oceanic ridges and in the Himalayan area the M c values are slightly lower (between 5.1 and 5.3), probably due to an arbitrary choice by the catalogue’s authors to select lower minimum cut-off magnitudes for such zones. The sampling radii and thus also spatial resolution vary, depending on earthquake density. Subduction zones typically emerge most clearly due to the higher seismic activity there, whereas continental and especially mid-ocean-ridge areas are least clearly resolved. Once M c had been estimated, we could compute b-values using the standard maximum-likelihood method (Aki, 1965). We found statistically significant and highly systematic spatial variations of b-values (Fig. 3A), consistent with the hypothesis and in line with previous works (Kagan, 1997; Nishikawa and Ide, 2014). The three major tectonic categories on Earth are distinctly different (Fig. 3A, B). We found the highest b-values along mid-ocean ridges, the

lowest b-values in subduction zones and intermediate b-values in continental zones typically characterised by a mixture of faulting styles. For example, when using Utsu’s test (Utsu, 1966), populations differed markedly at significance levels below 0.01 (see Table 1). Fig. 3B shows earthquake size-distributions in four selected locations. We also checked that the observed values were accurately described by a power-law behaviour (see Appendix A), performing a goodness-of-fit test for A, B, C, D distributions using a GR model (Aki, 1965). The fit-in-magnitude percentage obtained for all of them was greater than or equal to about 90%. To test whether the tectonic regime is indeed a statistically significant predictor of b-values, we first derived the distribution of b-values separately for each regime, following a previously published classification schema (Bird, 2003). The resulting distributions (coloured lines in Fig. 3C), including the global one (black), were not normally distributed, as confirmed by a KolmogorovSmirnov test (see Table 2). We then applied the non-parametric Wilcoxon signed rank test to identify whether the earthquake populations had been sampled from continuous distributions with equal medians (at a significance level of 5% − see Table 2).

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

5

Fig. 3. A) Map of b-values for the Global CMT (1980–2016) catalogue (for underlying R and M c maps see Fig. 2). B) Earthquake size-distributions corresponding to points A (R = 200 km, M c = 5.3, fit = 96.1776%), B (R = 500 km, M c = 5.4, fit = 89.1967%), C (R = 600 km, M c = 5.2, fit = 92.7667%), and D (R = 1200 km, M c = 5.2% fit = 93.0818%) of a) and Fig. 2. C) Cumulative density functions (CDFs) of b-values of a) sampled, moving along the PB2002, within a radius of 50 km. D) Mean b-values for SLAB 1.0 subduction zones, error bars: standard deviation, red: Mariana zone (high b-values), blue: Chilean zone (low b-values). Table 1 Utsu test P b , where P b is the probability of the b-values for A, B, C and D columns in Figs. 3A and B (second column) being equal. Figures in bold and italic indicate that the null hypothesis is 5% significant; bold alone indicates a 1% significant null hypothesis (i.e. ‘significantly’ and ‘highly significantly’ different b-values). The first column contains the number of complete events used to estimate b-values. P b of 1 denotes a consistency check while 0 indicates a value lower than 0.01. Number of events

b-Value

Points

A

B

C

D

309 152 126 102

0.83 ± 0.10 1.02 ± 0.07 1.12 ± 0.02 1.66 ± 0.06

A (subduction) B (continent) C (ocean) D (ocean)

1 0.192 0.02 0

0.192 1 0.586 0.004

0.02 0.586 1 0.015

0 0.004 0.015 1

We determined that, statistically, subduction earthquake sizedistribution (blue) is significantly different from oceanic distribution (green), which in turn differs from continental distribution (red). Accordingly, tectonic regime is indeed a predictor of b-values. However, subduction zones vary greatly in their struc-

Table 2 Wilcoxon and Kolmogorov-Smirnov tests (last column) for the PDFs of Fig. 3C (5% significance level). The Wilcoxon test computes the probability of being wrong in rejecting the null hypotheses that two distributions stem from continuous distributions with the same medians, while the Kolmogorov-Smirnov test assumes normally distributed data as null hypotheses. Probabilities of 1 denote a consistency check (100% probability of success if we discard the null hypothesis that two identical b-values are instead different), while 0 (0% probability of success) indicates a value lower than 0.01. Wilcoxon test PDF

Continents

Oceanic

Subduction

Kolmogorov-Smirnov test

Continents Oceanic Subduction

1 0 0.06

0 1 0

0.06 0 1

0 0 0

ture and characteristics. A first-order classification (Uyeda, 1982; Uyeda and Kanamori, 1979) distinguishes between Chilean-type subduction (fast-moving, young subducting lithosphere with strong compression along strongly coupled interfaces) and Mariana-type

6

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

Fig. 4. Detailed views of different seismotectonic structures, illustrating b-values in a map view (from Fig. 1A) and ternary diagrams. A) Mariana subduction zone: mostly NR and TH events, all with high b-values. B) Chilean subduction zone: most TH events display low b-values (high b-value SS and NR events occur along oceanic ridge faults in the southwest). C) Oceanic ridge: mostly NR, some SS events with high and very high b-values. D) Continental (Himalayan) collision zone: mixture of event types with intermediate b-values.

subduction (slow-moving, old subducting lithosphere accompanied by extensional roll-back mechanisms and back-arc spreading, but not associated with major earthquakes due to low coupling). We found statistically distinctive global b-values for these specific zones (Fig. 2D), with Chilean-type subduction zones tending to have low b-values, and Mariana-type subduction zones higher values, in accordance with the differential-stress b-value hypothesis.

4. Subduction zones, continental and mid-oceanic b-values Zooming in on some of the main seismotectonic structures depicted in Fig. 3 further backed our hypotheses. The Mariana subduction zone (Fig. 4A) is part of a convergent oceanic margin roughly 2,800 km in length. Here, the old seafloor is subducting deep into the crust, giving rise to widespread volcanic activity and hydrothermal emissions and a high level of

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

7

Fig. 5. A) Map view of Global CMT FMs (first nodal plane). Green: NR mechanisms, blue: TH mechanisms, red: SS mechanisms. B) Variations in b-values with rake angle λ for bin widths γ = 40◦ (left) and γ = 20◦ (right). Grey: results from S2005 (1980–2004, z = 0–50 km, M c ≥ 5.5, γ = 40◦ ). Vertical bars: Shi and Bolt uncertainties (Shi and Bolt, 1982). Horizontal grey lines: b-values for the entire S2005 datasets (dashed) and for the updated catalogue (solid).

very deep earthquake activity compared to other subduction zones (Stern et al., 2003). By using the ternary diagram (Fig. 4A) both compression and extension can coexist in such zones. Indeed, down-dip tension in the lower zone and down-dip compression in the upper zone have been previously described (Samowitz and Forsyth, 1981): thermal stresses and inelastic unbending of the upper part of the slab cause low-magnitude shallow seismicity, consistent with generally high b-values. The Chilean subduction zone (Fig. 4B) originates from the subduction of the Nazca Plate beneath the South America Plate and is associated with widespread and relatively high-magnitude TH seismicity, including the largest ever observed earthquake, with a magnitude of M = 9.6. In this regime, strong coupling between the overriding plates probably causes very high differential stresses, so low b-values are observed. The offshore spreading ridge and fracture zone to the south (characterised by NR and SS seismicity) exhibit distinctly higher b-values. However, since subduction zones can also slip aseismically, distinct subduction zones may exhibit differences in their seismiccoupling convergence rate relation (Stein and Okal, 2007). Seismicity along the spreading rifts and transform faults of midoceanic ridges (Fig. 4C) has previously been associated with higher b-values than the global average (Kagan, 1997). Given these systems’ low coupling coefficient, the bulk of deformation along them occurs aseismically. While spreading ridges are characterised by volcanic activity and associated NR seismicity, transform faults generate numerous slow earthquakes, their length and linearity triggering rather small, strike-slip earthquakes (Boettcher and Jor-

dan, 2004). High pore pressures and the extensional regime suggest low differential stresses in these regions, which are consistently reflected by high measured b-values. Continental collision zones are typified by strong tectonic heterogeneity on spatial scales that exceed the resolution capabilities of this global study. Consequently, no distinction can be made here between individual local regimes and their imprints on local b-values. Instead, we observe overall intermediate b-values for continental collision boundaries, partially with a tendency towards lower b-values e.g. in the Himalayan region (e.g., Fig. 4D), but we also note a tendency towards higher b-values in the European area. Continental rift systems like the East African Rift are characterised by high b-values, as expected in an extensional regime. 5. Global b dependence on tectonic styles: rake and plunge angles Focal mechanisms (FMs) by themselves can also be used to test the influence of differential stress on b-values (Fig. 5A). This approach was first proposed by Schorlemmer et al. (2005) (hereafter referred to as S2005), who demonstrated b-values’ systematic dependence on the FM’s rake angle, λ. To test whether faulting style and hence differential stress systematically influence b-values, we validated the b(λ) dependency observed by S2005, using independent data collected by global seismic networks over the past 10 yrs. We used the longer period of data now available, compared to S2005 (2004 update), to reassess the systematic dependency of b-values on earthquakes’ rake angle. We followed the identical procedure to S2005, barring one modification: instead of using over-

8

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

Fig. 6. a) The Global CMT FM ternary-based b-values diagram. b). A) b-value mapping inside the FM triangle, using the nearest neighbours approach (see MM): black dotted lines delineate areas of almost ‘pure’ NR, SS, and TH events. B) and D) Association between rake values (dotted vertical lines ±90◦ , ±45◦ , 0◦ ) with a sampling bin of 20◦ (as in Fig. 5B, right) and representation in the ternary plot for those selected rake bins. C) Earthquake size-distributions for subsets of B: the a-value for TH-SS is doubled for plotting purposes.

lapping rake bins with a width, γ , of 40◦ (Fig. 5B, left panel), the higher number of events (12,448 compared to the original 7,636) allowed us to reduce this value to 20◦ (Fig. 5B, right panel), revealing that within a given style there is a dependency on differential stress. For this analysis, we cut off the entire catalogue at the overall completeness magnitude of M c = 5.5 and show the analysis for plane 1 over the entire study period (1980–2016). Global seismotectonic structures are clearly reflected in the distribution of FM rake angles (Fig. 5A): besides abundant plateinterface TH seismicity in subduction zones (blue, λ ∼ 90◦ ), NR events (green, λ ∼ −90◦ ) are characteristic for ridges, shallow offtrenches and continental rifts, while SS events (red, λ ∼ 0◦ , approximately 180◦ ) occur partially along ridges and continental margins. Reassessing the rake-angle dependence of b-values, the results obtained by S2005 are shown in grey and are almost identical to the rake-angle dependence obtained when data from the last 10 yrs are added (γ = 40◦ , Fig. 3B, left). When γ = 20◦ (Fig. 5B, right) instead, we found that peak b-values are offset from ‘pure’ NR events (λ = −90◦ ) by approximately 45◦ . Using the Wilcoxon test, we formally tested the null hypothesis (no rake dependence for b-values) against the alternative scenario hypothesized according to S2005. We compared the residual distributions (summed on all λi ) of our data (γ = 40◦ ) with respect to S2005 and ‘constant b’ models, by bootstrapping 1,000 times. The two histograms (see Fig. S1) differ by a level of significance smaller than 0.01, meaning there is a 99% probability they are different. We have also verified the same b–λ dependence for both nodal planes and different periods individually (see Fig. S2). Similar to S2005, we used a depth layer of 0 to 50 km, while varying the period (1980–2004, 1980–2016 and 2005–2016 on different rows) and parameter γ (±20◦ , ±30◦ , ±40◦ in different columns). All frames display b-values associated with magnitude data for the nodal planes 1 and 2 (as indicated arbitrarily by the authors of the catalogue Ekström et al. (2012)). The analysed pe-

riods were the original one in S2005 (1980–2004, for which we set M c at 5.5, indicated with a grey frame), the updated one (1980–2016, where again M c = 5.5) and the new, independent period (2005–2016, where we set M c at 5.2). The same overall shape emerged for all periods. Being able to use smaller γ -values of ±20◦ revealed (Fig. 3b) that b-rake dependency is not simply sinusoidal (Petruccelli et al., 2018) or multimodal (Roberts et al., 2016), but peaked for normal events at −150◦ and −45◦ . We also analysed the faulting-style dependence of b in detail by introducing a more advanced b-value analysis. Instead of the rake angle, the three plunge angles δ P , δ T , and δ B can also be used to infer the tectonic style of a FM. We plotted each earthquake above the overall completeness level (5.5) on the ternary diagram and computed a corresponding b-value by sampling its 500 closest neighbouring events on the diagram (Fig. 6). Since events in a ternary diagram tend to cluster, we did not use samples from an equal tessellation, which might have caused major differences in the number of selected events, but instead adopted the closest neighbour selection strategy to preserve b-value pattern continuity. We sampled neighbouring earthquakes (Fig. 6A, the smoothing effects of our procedure showed is illustrated in Fig. S3) in a ternary diagram representation of plunge angles δ P T B , allowing us not only to see rake angle λ dependency (S2005), but also the full complexity of focal mechanisms (FM), including also the possible combination of two styles (mixed style). This analysis confirmed the first-order pattern seen in the global survey, with higher, intermediate, and lower b-values for NR (left corner), SS (top corner), and TH (right corner) events respectively. However, this representation reveals previously unknown systematic b-value variations along the bottom edge of the ternary diagram, when transitioning from NR and TH earthquakes. The very lowest b-values (b = 0.7) were observed for low-angle thrust faults (50◦  δ T  60◦ , 5◦  β  15◦ ) while the highest b-values (b = 1.4) occurred for steep normal faults (30◦  δ T  40◦ , 75◦  β 

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

9

Fig. 7. A) Analytical comparison of differential stress (Eq. (2)) scaled by lithostatic load for ‘pure’ dip-slip (λ = ±90◦ ). B) β(δ T ) relations: optimal planes (NR and TH corners) are expected to dip at β = 45◦ , while, moving towards the centre of the diagram, one orientation rotates towards the vertical (β = 90◦ ) and the other towards the horizontal (β = 0◦ ) (see text). C) NR (δ B ≤ 5◦ , δ T ≤ 40◦ , see triangle above) GCMT b-values of Equation (4) with different μ on plunge angles δ T or dip angles β1,2 (bottom, see E)). Parameters: bNR = 0.14, kNR = 0.005, friction values μ = 0.6 (solid), 0.7, 0.8 (dashed), ΛNR = 0.3, zref = 15 km. D) Same as C) but for TH (δ B ≤ 5◦ , δ T ≥ 55◦ ) regime (bTH = 1.20, kTH = 0.00002, ΛTH = 0.4, zref = 30 km). E) Dip-slip faulting schema for NR and TH regimes: arrows describe the reactivation processes of dip-slip faults with dip β : β1 refers to a low-dip plane, β2 to a high-dip plane.

85◦ ). Figs. 6B and 6D illustrate how the ternary diagram (Fig. 6A) and rake-based analyses of FM introduced in 2005 (Fig. 5B, right) are related: the FMs’ sampled rake windows roughly span triangular δ sections inside the ternary diagram, thus equally reflecting the differences in b-values for NR and TH events (Fig. 6C). Note that ‘pure’ dip-slip faults (λ = ±90◦ ) are located at the bottom of the triangle, whereas ‘pure’ strike-slips (λ = 0, ±180◦ ) vertically align in the diagram. 6. A unified framework to link b-values, faulting style and differential stress To validate our theoretical formulation against our observations, we first plotted the dependence of equation (2), illustrating that TH regimes require a larger differential stress than NR regimes for reactivation (Copley, 2018). The comparison is shown in Fig. 7A,

where we plot the ratio σ /σlith as function of δ T for a coefficient of friction (μ = 0.6). More specifically, this dependence could be separated into four different scenarios (Fig. 7B) because the β –δT relationship leads to this number of possibilities (see equation (A.17)-(A.18)): steeply dipping planes (grey curves, β > 45◦ ) and shallowly dipping planes (black curves, β < 45◦ ) for each of the two tectonic styles (NR, with δ T < 45◦ and TH, with δ T > 45◦ ). Corresponding to the NR–TH corners (δ T ∼ 0–90◦ ), both planes tend to dip 45◦ (Fig. 7B) but, as they move towards the centre of the diagram, one solution (grey) tends towards 90◦ , while the complementary one (black) tends towards 0◦ . In other words, our model describes the reactivation cycle of a dip-slip fault with β (or δ T ): a high-dip NR is likely to fail with respect to a low-dip NR, while physically a low-dip TH is more likely to fail with respect to a high-dip TH.

10

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

Figs. 7C and D compare the proposed analytical solutions (see equation (4), corresponding lines and regimes of Fig. 7A) for b-values as a function of δ T or β for three different friction values and for NR and TH regimes respectively. In this case, as described above, we use a zref and Λ specific for the given regime (30 km and 0.4 for TH and 15 km and 0.3 for NR, respectively). The variation of δ T or β is plotted in Fig. 7E, where grey and black refer to high and low angles, respectively (Fig. 7B). Observations are plotted as coloured dots and overall match the analytical solutions very closely. Equation (4) specifically predicts the highest b-values (Fig. 7C, grey curves) for steep normal faults (Fig. 7E, left side, grey faulting scheme), and the lowest ones (Fig. 7D, black curves) for shallow thrust faults (Fig. 7E, right side, black faulting scheme), also found in the data. The results also explain the observation along the bottom row of the ternary analysis (Fig. 6). In addition, our observations also match the Thrust-fault paradox: over-pressurised (ΛTH > ΛNR ) TH results in the lowest b-values, i.e. in the potential of nucleating megathrust earthquakes. The assumption of two different reference depths for NR and TH also suggests a further, overall analysis of the b-value dependence on depth to cross-check if the Earth’s strength profile is confirmed on global scale (see Fig. S4). Note that data from δ T = 40◦ to δ T = 55◦ (grey points in the ternary diagram of Fig. 7) should not be considered, because purely horizontal or vertical motion cases (δ T ∼ = 45◦ ) represent a singularity of the analytical problem. In fact, here the differential stresses required for fracturing would provide unrealistic solutions, as both mechanisms tend towards the vertical or horizontal (β = 90◦ , 0◦ , Fig. 7B) and the tectonic stress, σ , required for faulting, tends towards infinity (Fig. 7A). Strike-slip faults fill the upper part of the ternary diagram with intermediate b-values (0.9–1.0, Fig. 6), matching the intermediate differential stress levels predicted by equations (3) and (4). However, while for dip-slip fault zones we can rely on the fact that one of the principal stress is vertical (or almost vertical), for strike-slip faults no direct, unequivocal link can be established between the strike of the slip vector (which determines the amount of differential stress for a SS) and the plunges (Célérier, 2010). Therefore, we cannot derive a unique equation relating differential stress to faulting style for primarily strike-slip earthquakes (see Appendix A). 7. Conclusions We developed a unifying, empirical framework that links Anderson’s faulting theory and differential stress (Fig. 1) with the earthquake size-distribution. The results presented here are fully consistent with an inverse dependency of b-value on stress initially obtained in the controlled and reproducible conditions of rock laboratories (Scholz, 1968; Amitrano, 2003) for acoustic emissions of magnitudes −4 and ruptures of micro-metres. Thus, our study suggests that the stress dependency of the earthquake sizedistribution universally scales up to the largest earthquakes of magnitude 5.5-9 investigated here. Based on independent data, our results confirm the hypothesis made 14 yrs ago by S2005, that b-values depend on rake (Fig. 5), making it extremely likely that global b-values vary systematically with faulting style and in a statistically-significant way (Fig. S1). The analytical link derived between stress and b-values (see equations (2) and (4) and Fig. 7) is consistent with the pattern of globally observed b-values (Figs. 3 and 4) and explains first-order observations (Fig. 5) of b-values with faulting style (Petruccelli et al., 2018). Our new approach also specifically shows that, contrary to previous assumptions (S2005), the lowest b-values should not be observed for ‘pure’ TH events with a 45◦ dip, but should rather occur for very shallowly dipping TH events (Figs. 6A and 6D), such as Chilean subduction zones (Fig. 4B). Likewise, our model predicts

that the highest b-values should occur for NR to SS regimes, not for ‘pure’ NR events (Figs. 6A and 7C). In addition, our analytical solution accounts for a variable scenario with pressurised fluid, as highlighted by the so-called ‘thrust-fault paradox’ (Hubbert and Rubey, 1959). The analysis of b-values and faulting mechanisms is prone to many uncertainties in the data, regarding location and magnitude or systematic biases (Kagan, 1999, 2010). Moreover, faulting styles and tectonic regimes are sometimes poorly known and highly heterogeneous. However, we find it remarkable overall how systematically b-values vary with faulting style (Figs. 6 and 7) and tectonic region (Figs. 3, 4 and 5). Even more remarkable, in our view, is the fact that empirical observations can be explained by a simple geomechanical formulation based solely on Anderson’s faulting theory from 1905 and the Mohr-Coulomb failure criterion, which dates back to 1776. We established this framework based on global data, because their broadest range of tectonic provinces and FM diversity covered a homogeneous dataset. However, the underlying principles are universal and should apply to earthquakes on all scales. Accordingly, one of the first implications of our research is that we provide a quantitative framework for regional and local applications for linking earthquake size-distribution and tectonic regimes. The model we derived can and should now be tested using independent data, either using regional, high-quality earthquake catalogues or in lab-based studies linking stress and the faulting styles of acoustic emissions. We also suggest that the ternary diagram analysis introduced in Fig. 6 constitutes a powerful new way of studying b-values and their relationship to faulting on all scales. In conclusion, b-values are used in all studies of seismic hazard as a tool to determine the rate of recurrence of events of all sizes, in particular larger events that may dominate the hazard. We suggest that the overwhelming evidence of b-value dependency on differential stress (with appropriate hypotheses) and tectonic style necessitates a shift in how seismic hazard assessment (SHA) is conducted at the local, regional and global scales. Essentially, all SHA studies conducted to date, and most earthquake forecast models, use only empirical observations of b-values for the region, because so far no theory has been able to link it to geological properties. We believe that our results open up a pathway for a much deeper, geologically inspired interpretation of observed spatially varying b-values, enabling the construction of different kinds of seismogenic source model that combine knowledge of tectonics, observed faulting styles and past observed earthquakes in a single, consistent model. This will be especially important for regions with low-to-moderate seismicity, where empirical data provide poor constraints on the frequency of the largest events. Appendix A A.1. Dip-slip fault modelling According to Anderson’s theory, (as an absolute value) the applied shear stress τ for dip-slip faults is:

1

|τ | = σ sin 2θ

(A.1)

2

where σ is the tectonic stress and the angle θ is the complementary angle to the dip angle β = π /2 − θ . With respect to β , this results in:

1

1

1

2

2

2

|τ | = σ sin 2θ = σ sin[π − 2β] = σ sin 2β

(A.2)

We express Anderson’s faulting criterion for dip-slip faults in terms of δ T along the bottom of the ternary diagram (Figs. 1, 6, 7). The plunge angle δ T is linked to the vertical component of the T -axis vector t z (Gasperini and Vannucci, 2003):

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

δT = arcsin t z

(A.3)

that can be computed by calculating the sum of vertical components of the outwardly normal n z and of the slip vector d z :

tz =

nz + dz

√

(A.4)

2

Such vectors are a function of both the dip and rake angles:

11

So, ‘pure’ NR and TH dip at 45◦ , while the dip of one plane, moving towards the centre, increases to 90◦ while the other decreases to 0◦ . We can assume that the vertical stress σz is always lithostatic pressure, whereas horizontal stress σy varies with tectonic stress (assuming that stresses are positive for TH and negative for NR):

σz = ρ gz

(A.19) (A.20)

n z = − cos β

(A.5)

σ y = ρ gz + σ

d z = − sin β sin λ

(A.6)

Since, in this instance, vertical and horizontal stress also equal the maximum and minimum principal stress, tectonic and differential stress are the same (i.e. σ = σ1 − σ3 ). Assuming this configuration (Fig. 1), deriving the normal and shear stress along the fault zone is quite a straightforward calculation. We consider a dip-slip fault dipping at an angle β ranging from (0◦ to 90◦ ) (vertical directed downward). The components of the outward normal versor n are then:

Accordingly, the plunge angle δ T , as a function of both the dip and rake angles, is:



δT = arcsin

− cos β − sin β sin λ √



(A.7)

2

Along the bottom line connecting the NR vertex and the TH vertex, the rake angle is λ = ±90◦ (Célérier, 2010), so:



δT = arcsin





− cos β ∓ sin β cos β ± sin β = arcsin − √ √ 2

2



(A.8)

where the plus sign (+) applies to TH and the minus sign (−) applies to NR faulting. Equation (A.8) can be rewritten as follows:



√

2 sin δ T =

− cos β + sin β NR − cos β − sin β TH

(A.9)

The linear combination, or harmonic addition, of sine and cosine waves is equivalent to a single sine wave with a phase term ϕ and an amplitude factor A:

a1 sin α + a2 cos α = A sin(α + ϕ )

⎛

⎞ − sin ϕ sin β nˆ = ⎝ cos ϕ sin β ⎠ − cos β

(A.21)

where angle ϕ is the strike direction. Since differential stress for dip-slip configurations does not depend on strike direction (see main text), unlike in the strike-slip scenario we are allowed, for simplicity’s sake, to set ϕ = 0◦ . The traction components T i may then be written as:

T i = σi j n j =



σx σ y σ z



⎛

⎞

⎛

⎞

0 0 ⎝ sin β ⎠ = ⎝ σ y sin β ⎠ − cos β −σz cos β

(A.10)

(A.22)

where the original amplitudes a1 and a2 summed in quadrature provide A:

The scalar product of the traction and the outward normal vector provides the normal stress σn :

A=

σn = T · nˆ = σ y sin2 β + σz cos2 β



a21 + a22

and



ϕ = atan

a2

(A.11)

= σy

(A.12)

a1

Thus, TH provides two solutions for β as a function of δ T (remembering that sin α = sin(π − α ))

√

√

√

√

2 sin δ T =





2 sin β + 45◦ ⇒ β = δ T − 45◦



  2 sin δ T = 2 sin 180 − β + 45◦ √   = 2 sin 135◦ − β ⇒ β = 135◦ − δT

(A.13)

◦

(A.14)

=

1 − cos 2β 2

σz + σ y 2

+

+ σz

1 + cos 2β

σz − σ y 2

2 (A.23)

cos 2β

where we used duplication formulas for sine and cosine. We assume that the shear stress versor τ is oriented as the slip vector d (see Fig. 1):

⎛

⎞

⎛

Similarly, for the NR regime the phase terms are −45◦ and 135◦

√

√

√

√



 ◦

⇒ β = δT + 45◦ (A.15) √ 

 ◦  ◦ ◦ 2 sin δ T = 2 sin β + 135 = 2 sin 180 − 45 − β √   (A.16) = 2 sin 45◦ − β ⇒ β = 45◦ − δT 2 sin δ T =

2 sin β − 45



The four solutions (A.13)–(A.16) for β (two dip angles for two nodal planes and two tectonic styles), which describe how to switch from plunge δ T to β (and vice versa) through linear relations, are displayed in Fig. 5B and can be summarised as:





NR 0◦ ≤ δ T ≤ 45◦ β(δ T ) =



◦

 ◦



TH 45 ≤ δ T ≤ 90 β(δ T ) =

45◦ − δ T 45◦ + δ T



δT

− 45◦

135◦ − δ T

(A.17)

(A.24) Shear stress is then derived from the scalar product of the traction and shear-stress versor:

τ = τ · nˆ = −σ y sin λ sin β cos β + σz sin λ sin β cos β =

1 2

sin 2β sin λ(σz − σ y )

(A.25)

where we used a multiplication formula for the sine. For ‘pure’ dip-slip mechanisms (λ = ±90◦ , lower part of the triangle, upper sign for TH, lower sign for NR), (A.25) becomes:

1

τ = ± sin 2β(σz − σ y ) 2

(A.18)

⎞

sin λ cos β sin ψ + cos λ cos ψ cos λ τ = ⎝− sin λ cos β cos ψ + cos λ sin ψ ⎠ = ⎝− sin λ cos β ⎠ − sin β sin λ − sin β sin λ

By taking the expression of the stress in (A.19)-(A.20):

(A.26)

12

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

σn = ρ gz + τ =∓

σ 2

σ 2

(1 − cos 2β)

(A.27) (A.28)

sin 2β

(which we assume as being the lithostatic stress) is expressed by the constant:

Λ=

p

σzz

=

p

(A.38)

σlith

assuming positive stresses for TH (σ > 0) and negative stresses for NR (σ < 0). Using a Mohr-Coulomb criterion (Coulomb, 1776) defined as:

the lithostatic stress in equation (2) turns into ‘effective’ lithostatic stress:

|τ | = C + μ(σn − p )

 σlith = σlith ( zref )[1 − Λ] = ρ gzref [1 − Λ]

(A.29)

where C is cohesion, μ is the frictional coefficient (between 0.6 and 1 for most rocks Turcotte and Schubert, 2002) and p is pore pressure (considered for the time being as hydrostatic and used for completeness). Substituting the values for σn and τ as a function of σ gives us:

   σ  σ |τ | = ∓ sin 2β  = ± sin 2β 2

(A.30)

2

with the plus sign (+) for TH (σ > 0) and the minus sign (−) for NR (σ < 0). Then:

±

σ 2

sin 2β = C + μ(ρ gz − p ) + μ

σ 2

(1 − cos 2β)

(A.31)

Finally, calculating σ gives us the equation:

σ =

2μ(ρ gz − p ) ± sin 2β − μ(1 − cos 2β)

(A.32)

By setting σlith = (ρ gz − p ), we obtain the plot for σσ shown in lith Fig. 7A. Equation (A.32) admits real values if:

± sin 2β − μ(1 − cos 2β) = 0

(A.33)

This can be rewritten as:

± sin 2β + μ cos 2β = μ

(A.34)

According to harmonic addition – and remembering that the tangent has a period of 180◦ –, the left-hand term in (A.28) can be separated in:

 

1 + μ2 sin[±2β + atan μ] = μ

(A.35)

1 + μ2 sin[±2β + atan μ − π ] = μ

(A.36)

μ) = μ/√(1 + μ2)   sin[±2β + atan μ] = sin[atan μ] β = 0◦ ⇒ sin[±2β + atan μ − π ] = sin[atan μ] β = 90◦

(A.39)

where Λ ranges from 0 to 1 and lithostatic stress is determined at a specific reference depth zref . Hubbert and Rubey (1959) also show that the minimum dip angle β for gravitational sliding is related to the angle for optimal failure Φ :

tan β = (1 − Λ) tan Φ

(A.40)

Equation (A.40) describes the overthrust phenomenon: as Λ → 1 (p → σzz ) then β → 0, meaning that low-dip (thrust) fault, due to overpressure, is paradoxically locked, because it has already failed before completing the slip. When Λ → 0 (p → 0) instead, then β → 45◦ and tan β ∼ tan , meaning that fault zones are optimally oriented for failure, so they creep. So over-pressure is higher for low-dip TH faulting, dipping more than the ‘optimal’ 30◦ (if we assume a friction coefficient of about 0.6). Typical values for Λ, for low-dipping TH, range from 0.2 to 0.6 for different subduction zones. We set the Λ-factor for TH at 0.4, while for NR it is legitimate to assume hydrostatic pore pressure (for a rock density of 2,700 kg/m3 ), where Λ is approx. 0.3. Admitting over-pressurisation to the TH faulting regime (see (A.32)) means that equation (2) (main text) can be rewritten as:

σ =

 2μσlith

± sin 2β − μ(1 − cos 2β)

=

2μρ gzref [1 − ΛT H  ] NR

± sin 2β − μ(1 − cos 2β) (A.41)

 comes from (A.39) and factor Λ is determined by faultwhere σlith ing style. The reference depths (assumed here to be 15 km and 30 km), for TH and NR regimes, are estimated by averaging the two populations of FM in Fig. 7C and D.

A.3. Strike-slip fault modelling

But remember that sin (atan

(A.37)

for both tectonic styles. Using (A.17)-(A.18), these values correspond to δ T = 45◦ , representing a singularity of the problem. A.2. Over-pressurisation for thrust-faulting regime (the thrust-fault paradox) Some geological models portray the effect of pore-pressure p as hydrostatic (Anderson, 1905). In some cases, like in continental environments or along oceanic ridges, this assumption is justified. But in subduction zones, many authors tend to agree that TH faulting requires elevated pore pressures to explain frictional sliding, an assertion commonly referred to as the ‘thrust-fault paradox’ (Hubbert and Rubey, 1959). Hubbert and Rubey (1959) showed that if the ratio between pore pressure and vertical stress σzz

As already reported in the manuscript, dip-slip modelling is feasible because in this configuration at least one of principal axes is assumed to be vertical. Unfortunately, the same cannot apply to other zones in the ternary diagram, because there are no univocal relationships (see Célérier, 2010). The centre of the diagram represents a mixture of different styles, while the inclined edges rely on different hypotheses for the stress states (NR-SS and TH-SS), as explained above. For example, for strike-slip modelling (the top of the triangle), where both σ1 and σ3 lie on the horizontal plane, the differential stress depends on a different angle (strike), say ψ (Anderson, 1905), as opposed to on dip β (close to 90◦ ):

σ S S =

2μ(ρ gz − p )

± sin 2ψ + μ cos(2ψ)

(A.42)

Using the Anderson theory of faulting, faulting starts ( ddψ [σ S S ] = 0) as soon as the condition

cot ψ = ±μ

(A.43)

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

is satisfied (see Turcotte and Schubert, 2002, equation 8.36). Considering sign conventions (ψ positive in top quadrants and negative in bottom quadrants), just for ‘pure’ pure-strike-slip (λ = 0, ±180◦ ), we ended up with a total of 8 possible expressions for linking a plunge (δ B in this case) and ψ (four quadrants and two possible slip directions, positive and negative). Moreover, it is impossible to distinguish between left-lateral mechanisms and right lateral mechanisms in the ternary diagram, as they inevitably fall in the same region (the top centre) (see also Célérier, 2010). In short, rather than being related to how much the fault can dip (as is the case for a dip-slip fault), the differential stress needed to reactivate a strike-slip fault will depend on the (horizontal) directions in which the fault might slip, if a strike angle (with respect to σ1 , for instance, ranging from −π to π ) is defined. The reactivation process associated with a strike-slip fault does not establish a direct, univocal link between the strike of the slip vector (which determines the amount of differential stress for a SS) and plunges on the B axis (for which b-values are known, from Fig. 3) differently than in the dip-slip scenario, where relations for the fault dip and plunges on the T axis are possible. Therefore, the limit in the SS scenario is represented by the metric of the ternary space, since its usage implies that the position of each event in the diagram has to depend on plunge angles, which are in turn related to the vertical components of the moment tensor principal axes, as well as on the impossibility of knowing the real directions of stress principal axes. A.4. GR b-value goodness-of-fit-test According to the GR model (Aki, 1965), the cumulative number of earthquakes detected at each binned magnitude M is:

N ( M ) = 10a−bM

(A.44)

For each frequency-magnitude distribution, by setting the maximum magnitude of the model as equal to the maximum magnitude detected, percentages of fit can be estimated as:



p FIT = 1 −

| M i − mi |  Mi

(A.45)

(where M and m are empirical and theoretical magnitudes). Appendix B. Supplementary material Supplementary material related to this article can be found online at https://doi.org/10.1016/j.epsl.2019.115791. References Aki, K., 1965. Maximum likelihood estimate of b in the formula log N = a − bM and its confidence limits. Bull. Earthq. Res. Inst. Univ. Tokyo. Amitrano, D., 2003. Brittle-ductile transition and associated seismicity: experimental and numerical studies and relationship with the b value. J. Geophys. Res. 108, 1–15. https://doi.org/10.1029/2001JB000680. Amitrano, D., 2012. Variability in the power-law distributions of rupture events. Eur. Phys. J. Spec. Top. 205, 199–215. https://doi.org/10.1140/epjst/e2012-01571-9. Anderson, E.M., 1905. The dynamics of faulting. Trans. Edinb. Geol. Soc. 8, 387–402. https://doi.org/10.1144/transed.8.3.387. Bird, P., 2003. An updated digital model of plate boundaries. Geochem. Geophys. Geosyst. 4. https://doi.org/10.1029/2001GC000252. Boettcher, M.S., Jordan, T.H., 2004. Earthquake scaling relations for mid-ocean ridge transform faults. J. Geophys. Res., Solid Earth. https://doi.org/10.1029/ 2004JB003110. Célérier, B., 2010. Remarks on the relationship between the tectonic regime, the rake of the slip vectors, the dip of the nodal planes, and the plunges of the P, B, and T axes of earthquake focal mechanisms. Tectonophysics 482, 42–49. https:// doi.org/10.1016/j.tecto.2009.03.006. Collettini, C., Tesei, T., Scuderi, M.M., Carpenter, B.M., Viti, C., 2019. Beyond Byerlee friction, weak faults and implications for slip behavior. Earth Planet. Sci. Lett. 519, 245–263. https://doi.org/10.1016/J.EPSL.2019.05.011.

13

Copley, A., 2018. The strength of earthquake-generating faults. Q. J. Geol. Soc. Lond. https://doi.org/10.1144/jgs2017-037. Coulomb, C.A., 1776. Essai sur une application des regles des maximis et minimis a quelquels problemesde statique relatifs, a la architecture. Mem. Acad. Roy. Div. Sav. 7, 343–387. Dziewonski, A.M., Chou, T.-A., Woodhouse, J.H., 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. J. Geophys. Res. 86, 2825–2852. https://doi.org/10.1029/ JB086iB04p02825. ´ Ekström, G., Nettles, M., Dziewonski, A.M., 2012. The global CMT project 2004-2010: centroid-moment tensors for 13,017 earthquakes. Phys. Earth Planet. Inter. 200–201, 1–9. https://doi.org/10.1016/j.pepi.2012.04.002. Frohlich, C., 2001. Display and quantitative assessment of distributions of earthquake focal mechanisms. Geophys. J. Int. 144, 300–308. https://doi.org/10.1046/j.1365246X.2001.00341.x. Frohlich, C., Apperson, K.D., 1992. Earthquake focal mechanisms, moment tensors, and the consistency of seismic activity near plate boundaries. Tectonics 11, 279–296. https://doi.org/10.1029/91TC02888. Gasperini, P., Vannucci, G., 2003. FPSPACK: a package of FORTRAN subroutines to manage earthquake focal mechanism data. Comput. Geosci. 29, 893–901. https:// doi.org/10.1016/S0098-3004(03)00096-7. Ghosh, A., Newman, A.V., Thomas, A.M., Farmer, G.T., 2008. Interface locking along the subduction megathrust from b-value mapping near Nicoya Peninsula, Costa Rica. Geophys. Res. Lett. 35. https://doi.org/10.1029/2007GL031617. Goebel, T.H.W., Schorlemmer, D., Becker, T.W., Dresen, G., Sammis, C.G., 2013. Acoustic emissions document stress changes over many seismic cycles in stick-slip experiments. Geophys. Res. Lett. 40, 2049–2054. https://doi.org/10.1002/grl.50507. Gulia, L., Wiemer, S., 2010. The influence of tectonic regimes on the earthquake size distribution: a case study for Italy. Geophys. Res. Lett. 37. https://doi.org/10. 1029/2010GL043066. Gutenberg, B., Richter, C.F., 1944. Frequency of earthquakes in California. Bull. Seismol. Soc. Am. 34, 185–188. Hayes, G.P., Wald, D.J., Johnson, R.L., 2012. Slab1.0: a three-dimensional model of global subduction zone geometries. J. Geophys. Res., Solid Earth 117. https:// doi.org/10.1029/2011JB008524. Hubbert, B.M., Rubey, W.W., 1959. Role of fluid pressure in mechanics of overthrust faulting, 1: Mechanics of fluid-filled porous solids and its application to overthrust faulting. Kagan, Y.Y., 1997. Seismic moment-frequency relation for shallow earthquakes: regional comparison. J. Geophys. Res. 102, 2835. https://doi.org/10.1029/ 96JB03386. Kagan, Y.Y., 1999. Universality of the seismic moment-frequency relation. Pure Appl. Geophys. 155, 537–573. https://doi.org/10.1007/s000240050277. Kagan, Y.Y., 2010. Earthquake size distribution: power-law with exponent β ≡ 1/2? Tectonophysics 490, 103–114. https://doi.org/10.1016/j.tecto.2010.04.034. Mogi, K., 1962. Study of elastic shocks caused by the fracture of heterogeneous materials and its relation to earthquake phenomena. Bull. Earthq. Res. Inst. Univ. Tokyo 40, 125–173. Nishikawa, T., Ide, S., 2014. Earthquake size distribution in subduction zones linked to slab buoyancy. Nat. Geosci. 7, 904–908. https://doi.org/10.1038/ngeo2279. Okal, E.A., Romanowicz, B.A., 1994. On the variation of b-values with earthquake size. Phys. Earth Planet. Inter. 87, 55–76. https://doi.org/10.1016/0031-9201(94) 90021-3. Petruccelli, A., Vannucci, G., Lolli, B., Gasperini, P., 2018. Harmonic fluctuation of the slope of the frequency–magnitude distribution (b-value) as a function of the angle of rake. Bull. Seismol. Soc. Am. 108. https://doi.org/10.1785/0120170328. Roberts, N.S., Bell, A.F., Main, I.G., 2016. Mode switching in volcanic seismicity: El Hierro 2011-2013. Geophys. Res. Lett. 43, 4288–4296. https://doi.org/10.1002/ 2016GL068809. Roberts, N.S., Bell, A.F., Main, I.G., 2015. Are volcanic seismic b-values high, and if so when? J. Volcanol. Geotherm. Res. https://doi.org/10.1016/j.jvolgeores.2015. 10.021. Sammonds, P.R., Meredith, P.G., Main, I.G., 1992. Role of pore fluids in the generation of seismic precursors to shear fracture. Nature 359, 228–230. https://doi.org/10. 1038/359228a0. Samowitz, I.R., Forsyth, D.W., 1981. Double seismic zone beneath the Mariana Island Arc. J. Geophys. Res., Solid Earth 86, 7013–7021. https://doi.org/10.1029/ JB086iB08p07013. Scholz, C.H., 1968. The frequency-magnitude relation of microfracturing in rock and its relation to earthquakes. Bull. Seismol. Soc. Am. 58, 399–415. Scholz, C.H., 2002. The Mechanics of Earthquakes and Faulting, 2nd edition. Cambridge University Press. Scholz, C.H., 2015. On the stress dependence of the earthquake b value. Geophys. Res. Lett. 42, 1399–1402. https://doi.org/10.1002/2014GL062863. Schorlemmer, D., Wiemer, S., 2005. Earth science: microseismicity data forecast rupture area. Nature 434, 1086. https://doi.org/10.1038/4341086a. Schorlemmer, D., Wiemer, S., Wyss, M., 2005. Variations in earthquake-size distribution across different stress regimes. Nature 437, 539–542. https://doi.org/10. 1038/nature04094.

14

A. Petruccelli et al. / Earth and Planetary Science Letters 527 (2019) 115791

Schorlemmer, D., Wiemer, S., Wyss, M., 2004. Earthquake statistics at Parkfield, 1: stationarity of b values. J. Geophys. Res., Solid Earth 109, 1–17. https://doi.org/ 10.1029/2004JB003234. Serpelloni, E., Vannucci, G., Pondrelli, S., Argnani, A., Casula, G., Anzidei, M., Baldi, P., Gasperini, P., 2007. Kinematics of the Western Africa-Eurasia plate boundary from focal mechanisms and GPS data. Geophys. J. Int. 169, 1180–1200. https:// doi.org/10.1111/j.1365-246X.2007.03367.x. Shi, Y., Bolt, B., 1982. The standard error of the magnitude-frequency b-value. Bull. Seismol. Soc. Am. 72, 1677–1687. Sibson, R.H., 2004. Controls on maximum fluid overpressure defining conditions for mesozonal mineralisation. J. Struct. Geol. https://doi.org/10.1016/j.jsg.2003. 11.003. Sibson, R.H., 2014. Earthquake rupturing in fluid-overpressured crust: how common? Pure Appl. Geophys. 171, 2867–2885. https://doi.org/10.1007/s00024014-0838-3. Spada, M., Tormann, T., Wiemer, S., Enescu, B., 2013. Generic dependence of the frequency-size distribution of earthquakes on depth and its relation to the strength profile of the crust. Geophys. Res. Lett. 40, 709–714. https://doi.org/ 10.1029/2012GL054198. Stein, S., Okal, E.A., 2007. Ultralong period seismic study of the December 2004 Indian Ocean earthquake and implications for regional tectonics and the subduction process. Bull. Seismol. Soc. Am. 97. https://doi.org/10.1785/0120050617. Stern, R.J., Fouch, M.J., Klemperer, S.L., 2003. An overview of the Izu-Bonin-Mariana subduction factory. Insid. Subduction Fact. 138, 175–222. https://doi.org/10. 1029/138GM10. Tormann, T., Enescu, B., Woessner, J., Wiemer, S., 2015. Randomness of megathrust earthquakes implied by rapid stress recovery after the Japan earthquake. Nat. Geosci. 8, 152–158. https://doi.org/10.1038/ngeo2343. Tormann, T., Wiemer, S., Hardebeck, J.L., 2012. Earthquake recurrence models fail when earthquakes fail to reset the stress field. Geophys. Res. Lett. 39. https:// doi.org/10.1029/2012GL052913.

Tormann, T., Wiemer, S., Mignan, A., 2014. Systematic survey of high-resolution b value imaging along Californian faults: inference on asperities. J. Geophys. Res., Solid Earth. https://doi.org/10.1002/2013JB010867. Townend, J., Zoback, M.D., 2000. How faulting keeps the crust strong. Geology. https://doi.org/10.1130/0091-7613(2000)028<0399:HFKTCS>2.3.CO;2. Turcotte, D.L., Schubert, G., 2002. Geodynamics, 2nd edition. Cambridge University Press. Utsu, T., 1966. A statistical significance test of the difference in b-value between two earthquake groups. J. Phys. Earth 14, 37–40. https://doi.org/10.4294/jpe1952.14. 37. Uyeda, S., 1982. Subduction zones: an introduction to comparative subductology. Tectonophysics 81, 133–159. https://doi.org/10.1016/0040-1951(82)90126-3. Uyeda, S., Kanamori, H., 1979. Back-arc opening and the mode of subduction. J. Geophys. Res. 84, 1049. https://doi.org/10.1029/JB084iB03p01049. Wiemer, S., Wyss, M., 2000. Minimum magnitude of completeness in earthquake catalogs: examples from Alaska, the Western United States and Japan. Bull. Seismol. Soc. Am. 90, 859–869. https://doi.org/10.1785/0119990114. Wiemer, S., Wyss, M., 2002. Mapping spatial variability of the frequency-magnitude distribution of earthquakes. Adv. Geophys. 45. https://doi.org/10.1016/S00652687(02)80007-3. Woessner, J., Wiemer, S., 2005. Assessing the quality of earthquake catalogues: estimating the magnitude of completeness and its uncertainty. Bull. Seismol. Soc. Am. 95, 684–698. https://doi.org/10.1785/0120040007. Yang, W., Hauksson, E., Shearer, P.M., 2012. Computing a large refined catalog of focal mechanisms for southern California (1981-2010): temporal stability of the style of faulting. Bull. Seismol. Soc. Am. 102, 1179–1194. https://doi.org/10. 1785/0120110311.