Dynamic rupture in the Earth is inherently complex. The material is heterogeneous, the physics is nonlinear, and the system response is intermittent and spatially and temporally variable. It is no wonder that the size of the largest event at a given time, or the precise time of catastrophic failure, can be hard to predict—even in laboratory ‘scale-model’ experiments (Vasseur et al., 2015).

Nevertheless, amongst all the potential chaos, there is an emergent order. For example, the population of events obey several well-defined scaling laws, including the relationship between the frequency and magnitude of seismic events, and the frequency and spacing between rupture locations. The former is known as the Gutenberg-Richter law, whose scaling exponent—the seismic ‘b-value’—is the slope on a plot of the logarithm of frequency and magnitude, and the latter the slope D2 on a plot of log frequency versus log of the distance between event locations. Magnitude is already a logarithmic measure of source size, so both are fundamentally power-laws, containing no characteristic length scale. This is consistent with the scale-free or self-similar nature of a plethora of geological structures, including mapped faults and fractures (Bonnet et al., 2001).

The relationship between seismic b-value and the underlying physical size of the source depends on the scaling of slip and rupture area and the characteristics of the sensor used to measure the maximum amplitude of the radiated wave. In the typical case of a sensor operating as a velocity transducer, and assuming a scale-invariant ratio of slip to source length, the b-value is also the exponent of the frequency-source area A relation: F(A) ∼ A–b, or F(l) ∼ l–D, where l is a characteristic rupture length, and the length exponent D = 2b (Kanamori and Anderson, 1975).

The exponent D2, known as the correlation dimension, measures the degree of localization: a uniform random distribution of event locations in three dimensions would result in a correlation dimension of 3, reducing to 2 for events randomly distributed on a plane, and 1 on a line. However, D2 can in principle take on any non-integer value between 0 and 3 (Hirata et al., 1987). The question is, what controls the exponents, and hence the degree of localization of deformation and the potential for large seismic ruptures?

The first clue came from laboratory experiments on initially intact materials with different degrees of heterogeneity (Mogi, 1962), where higher b-values were associated with more heterogeneous materials. This is intuitively appealing: heterogeneous materials have many more potential nucleation sites, but they also have many more potential barriers to rupture propagation, thereby favoring a greater proportion of smaller events (Segall and Pollard, 1980; Sammonds and Ohnaka, 1998). The second was the observation that for a given rock type, the b-value itself evolved during deformation, with a clear negative correlation between b and the differential stress on the sample boundary (Scholz, 1968). This is also intuitively appealing—a greater driving stress would be more likely to overcome local barriers during rupture, hence increasing the proportion of larger events. It is also consistent with the observation of systematic variations in b-value with focal mechanisms in field data from different tectonic stress regimes (Schorlemmer et al., 2005). To second order, Meredith and Atkinson (1983) showed that b-value was negatively correlated to the stress intensity, a measure of the degree of local stress concentration at the tip of the largest crack, normalized to its critical value at system-sized failure (the fracture toughness), where results from different materials collapse on the same linear trend. The slope of this trend depends on the chemical activity of the pore fluid, with higher partial pressures leading to higher b-values.

All of these experiments were carried out on initially-intact rock samples. However, Earth’s brittle crust already contains many large faults and fractures as sources of preexisting macroscopic structural heterogeneity. How do they control the relevant scaling exponents, and hence the potential for large events?

Goebel et al. (2017, p. 815 in this issue of Geology) address this issue in an ingenious set of controlled laboratory experiments. Using the same starting material, they first generate an ideally smooth through-going fault by sawing through the cylindrical sample at an optimal angle of ∼30° to the vertical axis, and then by polishing the two surfaces. In an intermediate case, they artificially roughen the two surfaces. For an ideally heterogeneous fault, they use a pre-fractured sample, which has both a rougher fault surface and a greater degree of off-fault damage. The motivation is that the rougher fault is by definition ‘young’, whereas the smoother saw-cuts may be more representative of more mature faults with greater degrees of wear (Stirling et al., 1996). Goebel et al. also introduce an important innovation. They use the variability of focal mechanisms, specifically the P-axis of the moment tensor, to reveal the degree of heterogeneity in the local stress orientation as a function of the starting structural heterogeneity. In order to isolate the effect of preexisting structural heterogeneity, they analyze results at boundary stresses near stick-slip failure.

Goebel et al. have taken a lot of care in experimental design and analysis of the results, including important details such as consideration of the threshold of completeness for the catalogues, analysis of the maximum principal component of strain inferred from the focal mechanisms, and examining the convergence of the model parameters.

The results are very clear. The rougher faults promote more spatially distributed deformation with higher correlation dimensions, higher b-values, and more variable focal mechanisms—reflecting a greater degree of local stress orientation heterogeneity in the pre-fractured samples. They also show that D2 ≈ 2b or D2D. This is a remarkable result. Other things being equal, this means the scaling exponents for source rupture length and the distance between ruptures are similar to each other. For example, if we distribute a set of non-overlapping source rupture areas of different sizes such as they each collectively occupy a plane, so that D2 ≈ 2, then we might also expect b = 1 or D = 2 (Kanamori and Anderson, 1975), as observed here in one of the intermediate cases. In the case of the polished fault surfaces, values of D2 < 2 imply a set of larger precursory fractures nucleating in clusters on the fault plane, consistent with the observed pattern of rupture locations imaged on the fault plane.

The outstanding question is how these results scale to the field case. There is now a clear hypothesis that, other things being equal, mature faults should have lower values of b and D2, and hence a greater potential for generating large ruptures. This is consistent with the results of Stirling et al. (1996) who present field evidence that the ratio of the recurrence rate of small to large earthquakes along a fault zone may decrease as slip accumulates and the fault becomes smoother, at least for strike-slip faults. We might also expect b to scale positively with D2. However, it may be difficult to isolate the effect of starting material heterogeneity in the field case. We cannot image it directly, and changes in the scaling exponent could reflect changes in the remote or local stress, or the pore pressure (Sammonds et al., 1992). As a consequence field results may not be so clear cut; for example, the correlation between b and D2 can be negative and/or quite weak in field examples (Henderson et al., 1992, 1994), possibly due to local stress concentrators or major asperities (Main, 1992). Structural heterogeneity may also include fault jogs or offsets not examined here, which often control rupture arrest (Wesnousky, 2006).

In the laboratory, the variations in b-value can be quite large and clearly statistically significant. However, the large literature on field studies of b-value is beset by questions of statistical significance. It is not always clear that the inferred b -value is representative of the long-term underlying value, i.e., that the statistics have converged (Frohlich and Davis, 1993). The same applies to identifying a real change in b-value, even in earthquake sequences where we would expect the stress field to have changed (Shcherbakov et al., 2012). This is exacerbated by the fact that the standard formula for the uncertainty in b-value uncertainty can significantly underestimate the total error, after propagating the contribution from the estimation of the magnitude threshold for complete reporting (Roberts et al., 2015).

Ultimately, the lack of control in the field case may make the problem of uniquely inferring the cause of changes in b-value difficult, but it is clear from the results of Goebel et al. that independent estimation of the correlation dimension and the variability in focal mechanisms can provide important constraints. Clearly, we need to take (account of) the rough with the smooth.