Elastic Impact Consequences for High‐Frequency Earthquake Ground Motion

A fundamental question of earthquake science is what produces damaging high‐frequency ground motion, with the classic Brune‐Haskell model postulating that abrupt fault slip causes it. However, even when amended with heterogeneous rupture, the model fails to explain observations of different sized repeating earthquakes and has challenges explaining high‐frequency radiation patterns. We propose an additional cause for high‐frequency earthquake spectra from elastic collisions of structures within a rupturing fault zone. The collision spectrum is set by an impact contact time proportional to the size of colliding structures so that spectra depend on fundamentally different physical parameters compared with slip models. When added to standard models, collisions can reconcile the discrepant observations since the size, shape, and orientation of structures vary between different fault zones but remain constant within a fault segment. High‐frequency earthquake ground motions and damage may therefore be an outgrowth of fault‐zone structure rather than sudden initiation of slip.


Introduction
It is commonly assumed that the spectra of earthquakes are determined by their spatiotemporal slip history (Kanamori & Anderson, 1975). The standard Brune-Haskell model (Brune, 1970;Haskell, 1964) assumes that earthquakes rupture a rectangular fault patch and that the decay of energy at high frequencies is related to the average rise time of slip on each subpatch (T r ) during a total duration of slip (T), with the spectrum falling off as f −2 at high frequencies above a corner frequency f 0 ≈ 1/T r (Figure 1a; Haskell, 1964). Within the context of this model, including factors related to dynamic rupture (Madariaga, 1976), the stress drop on the fault is related to the corner frequency and seismic moment (Allmann & Shearer, 2009;Kanamori & Anderson, 1975), enabling measurements of stress drop from seismological observations. In this framework, f 0 is interpreted to be proportional to rupture velocity and to the cube root of stress drop (Allmann & Shearer, 2009), and the spectral behavior is due to the sudden initiation (and ending) of fault slip (Brune, 1970;Haskell, 1964;Madariaga, 1976). It is widely acknowledged that the Brune-Haskell model is overly simplistic, and newer models have since added complexity to the Brune-Haskell picture by accounting for smaller-scale slip heterogeneity that is either stochastic (e.g., Aki, 1967;Andrews, 1980;Frankel, 1991;Haskell, 1966) or deterministic (e.g., Dunham et al., 2011); however, these newer models still fundamentally rely on the same physics as the Brune-Haskell model for producing high-frequency ground motion. For example, in Mai and Beroza (2002), Shi and Day (2013), and Graves and Pitarka (2016), it is still the abrupt changes in presumed frictionally governed slip on a large number of fault patches that determine both the corner frequency and the spectral falloff above the corner frequency. In other works, ground motions are matched empirically with assumed spectral shapes (e.g., Boore, 2003;Hanks & McGuire, 1981); these models are not based on source physics and thus have limitations for learning about the physics of earthquakes or making predictions about ground motions that go beyond the assumed empiricism.
While this conceptual picture has dominated seismological interpretation for the last five decades, certain repeating earthquake observations are difficult to reconcile within the fault slip framework. For example, although stress drop is interpreted to be independent of magnitude over seven orders of magnitude of slip globally, regional studies display a significant correlation between stress drop and magnitude (Abercrombie, 1995;Abercrombie & Rice, 2005). Some of this correlation is ascribed to observational bias (Ide & Beroza, 2001); however, observations of repeating earthquakes that span a wide range of magnitudes, but have nearly constant corner frequencies (Harrington & Brodsky, 2009;Lengline et al., 2014;Lin et al., 2012;Lin et al., 2016), demonstrate that stress drop determined through a Brune-Haskell analysis cannot be constant for certain groups of events (Figure 2c). This is challenging to explain since dynamical models  Brune-Haskell model (left) and the proposed elastic impact model (right). The top row shows schematics of the two models, with the spectrum determined by fault slip (moment tensor source) in the Brune-Haskell model (a) and by elastic impact of fault-zone structures (Faulkner et al., 2003;Swanson, 2006) (idealized as particles) in the impact model (b). The middle row shows the time history of slip on a patch for the Brune-Haskell model (rise time T r ) and the time history of impact force for the impact model (contact time T c ). The bottom row shows the idealized predicted far-field ground motion displacement spectra for the two models (rupture duration T, without attenuation). predict only modest variability in stress drop (Chen & Lapusta, 2009), and other explanations like extreme differences in geometry or rupture velocity (Lin & Lapusta, 2018) do not have compelling observational support (Harrington & Brodsky, 2009).
As discussed more in later sections, interpretations of earthquake spectra using fault slip models also fail to explain observed high-frequency radiation patterns (Graves & Pitarka, 2016;Somerville et al., 1997), stress drop variations between mature (interplate) and immature (intraplate) fault zones (Kanamori & Anderson, 1975), and shallow earthquakes with high stress drops (Abercrombie, 2014). Furthermore, fault zones are observed to be complex (Chester & Chester, 1998;Faulkner et al., 2003;Swanson, 2006), with a range of structures ( Figure 1b) that may collide with the rough sides of the fault zone during slip, preventing rupture from being governed purely by frictionally mediated fault slip (Figure 1a).
We propose that elastic collisions of structures within a rupturing fault zone could contribute to high-frequency seismic spectra. The physics of these elastic collisions is distinct from all previously considered physical models and may provide a physical basis for certain empirical stochastic ground motion models (Boore, 2003;Hanks & McGuire, 1981). We discuss how this new model, in conjunction with traditional slip models, can help explain the aforementioned discrepant observations. We discuss the implications for seismological interpretation of corner frequencies as well as for ground motions more generally.

Model
We propose that high-frequency seismic spectra may have significant contributions from the elastic impacts between rough surfaces of structures within a rupturing fault zone ( Figure 1b) (Faulkner et al., 2003;Swanson, 2006). To evaluate this contribution, we approximate the structures within the fault zone as individual particles that collide during slip, in analogy to a debris flow. Although fault slip occurs under much higher confining stress than debris flows, processes such as thermal pressurization (Rice, 2006), acoustic fluidization (Melosh, 1979), and brecciation/hydrofracture (Sibson, 1986;Sibson, 1989) may allow for structures to lose frictional contact and, subsequently, freely impact each other (despite limited void space) in a similar fashion to how shearing during debris flows causes bulk dilation and particle impacts rather than purely frictional slip (Iverson, 1997). Such impacts may be necessary to accommodate large-scale fault motion in the presence of small-scale geometric incompatibilities caused by rough fault-zone structure like sidewall ripouts (Swanson, 2006). While this evidence is only suggestive, and should be assessed further in future work, we believe it is nonetheless useful to describe the implications of such an impact model without a complete understanding of the precise way in which the collisions might be accommodated.
Hertzian elastic contact theory (Johnson, 2004), developed originally to describe the deformation of glass lenses (Hertz, 1881), has been applied to particle impacts in debris flows (Farin et al., 2019;Lai et al., 2018) and predicts that the timescale of impact, T c , is given by where m is the mass of the particle, R is the effective radius of curvature, E * = E/[2(1 − ν 2 )], E is Young's modulus, ν is Poisson's ratio, V z is impact speed, and ρ is particle density, and the second approximation assumes an ellipsoidal impactor with aspect ratio ξ > 1 impacting a flat surface (see Supporting Information). Debris flow models (Farin et al., 2019;Lai et al., 2018) indicate that impact speed scales with slip speed, V s , so that V z ≈ V s ≈ 1 m/s, implying that T c is very weakly dependent on slip speed and proportional to particle size. Elastic contact theory also predicts the force time history to be nearly a truncated sinusoid (Johnson, 2004;equation 11.25), that is, with force time history proportional to i t ð Þ ¼ cos πt Tc ·rect Tc t ð Þ (where rect is a unit boxcar function) (Figure 1b), which has a spectrum given by (Bracewell, 1978) I For the sum of independent impacts, the far-field unattenuated ground displacement spectrum is proportional to the impact force spectrum (see Supporting Information). The predicted displacement spectrum therefore has a corner frequency f 0 = 1/T c and a high-frequency decay of f −2 (Figure 3a), just as in the Brune-Haskell model, and similar overall spectral shape but with f 0 determined by particle size and shape rather than stress drop or the spatiotemporal slip distribution or total earthquake duration (see Supporting Information).
In the Brune-Haskell model and its heterogeneous slip variants, all high-frequency energy comes from abrupt initiation and ending of local fault slip (Brune, 1970;Mai & Beroza, 2002). If changes in slip rate are relatively smooth, then the high-frequency spectral decay can be dominated by elastic impacts. For example, for a smooth Gaussian source time function (STF) or moment rate function s(t), with s(t) = s Gauss (t; T r ) * s Gauss (t; T), where * denotes the convolution of two Gaussians, then the source spectrum is so that the STF spectrum decays exponentially above 1/T r rather than as f −2 as in the unmodified Brune-Haskell or omega square model (Figure 3a). Heterogeneous slip results in stronger, more continuous high-frequency radiation but also predicts an STF that decays faster than f −2 if the slip spectrum decays in wave number faster than is typically posited (Andrews, 1980;Mai & Beroza, 2002) or if asperity STFs are smoother than the assumed abrupt triangular or boxcar slip functions (Mai & Beroza, 2002). Indeed, recent dynamical simulations suggest that the amount of fault roughness required for dynamical slip models to match the f −2 spectrum (fault parallel amplitude-to-wavelength ratios of 10 −2 ; Dunham et al., 2011;Shi & Day, 2013) is significantly higher than observed (fault parallel amplitude-to-wavelength ratios of~10 −3 ; Power & Tullis, 1991). This suggests that it is worthwhile to consider the elastic impact model as a potential alternative cause for high-frequency ground motion. The impact model may therefore also provide a physical basis and interpretation for empirical stochastic ground motion models that posit an omega square source spectrum (e.g., Boore, 2003).

Results and Implications
For a reasonable range of parameters (see Supporting Information; Lay & Wallace, 1995), we predict that the zero-frequency ground motion amplitude is dominated by the finite fault slip of a smooth STF, while the high-frequency corner and subsequent decay can be dominated by the elastic impacts (Figure 3a), depending on how smooth the slip is assumed to be. For example, as described above, dynamical simulations with fault roughness in the observed range predict smoother slip than is required to explain the highest frequency ground motions (Figure 3a).
With the elastic impact model, the corner-frequency-inferred stress drop observations (Figure 2c) can be reconciled by noting that individual fault-zone segments have a set average particle size, while-globally -average fault-zone particle size increases with fault length (Figure 2). Earthquakes of different sizes can therefore share the same high-frequency spectral shape if they share the same fault-zone properties, as repeating earthquakes should, with Figure 3b showing an example of predicted spectra for two earthquakes of different magnitudes and durations but the same particle size R and otherwise smooth Gaussian slip. Thus, the elastic impact model can explain the constant corner frequencies of repeating earthquakes with distinctly different contact physics than explanations relying on repetitive slip on asperities of vastly different amplitudes (e.g., Harrington & Brodsky, 2009), which (at face value) does not occur in state-of-the-art simulations (Chen & Lapusta, 2009). Importantly, the impact model's cause of high-frequency radiation is independent of the true stress drop (see equation (1)); the apparent inferred stress drops from this model shown in Figure 2c are an artificial construct created from the predicted corner frequencies (see equation (2)) and applying the traditional stress drop-corner-frequency scaling.
Observations suggest that fault roughness is either self-similar (Hurst exponent H = 1) or that relative roughness decreases modestly with fault length (H ≈ 0.6 − 1.0) Candela et al., 2011;Power et al., 1987;Renard et al., 2006). Since roughness determines the size of fault structures (Figure 1b), these structural observations indicate that elastic impact-generated corner frequencies should be determined indirectly by the roughness within a given fault-zone segment (Figures 2a and 2c). Conversely, since High-frequency falloff of best fit to the Brune-Haskell plus elastic impact model has a slope of −2.1, very nearly equal to the −2 slope of the omega square model. The parameters used are T c = 0.1 s, f 0 = 10 Hz, and T r = 1 s. (b) Predicted spectral amplitudes for two ruptures with Gaussian STFs of different durations but equal particle size R (and thus equal contact time T c ). T = T r = 2 s (blue), and T = T r = 1 s (red). Amplitudes are scaled by the cube of duration (Duputel et al., 2013). Two-corner-frequency fits to the two spectra are shown, with higher corner frequencies being approximately equal (and close to 1/T c ). Area left of dashed line is irresolvable without broadband instrumentation.
observed structures generally increase with fault length (Figure 2b), Brune-Haskell-inferred stress drops are predicted to be nearly magnitude independent globally for the elastic impact model, with H = 0.95 fitting the aggregate (nonrepeating earthquake) observations reasonably well (Figure 2c). The best-fitting H being somewhat larger than observed field values suggests that either estimates of H are underestimated (e.g., the largest length scales at a given outcrop may be systematically smoother than interfault-segment roughness) or earthquake catalogs are incomplete even in the midmagnitude range. For a fault zone with a distribution of particle sizes, particles much larger than the median dominate the predicted signal (see Supporting Information) with the relatively small Chi-Chi repeaters (f 0 ≈ 20 Hz) requiring 73rd percentile particle sizes of 3.6 m and the larger Northridge aftershocks (f 0 ≈ 2 Hz) requiring 73rd percentile sizes of 36 m to explain the observations (see Supporting Information; Engelder, 1974).
There are several other implications of the elastic impact model. Perhaps most importantly, it implies that corner frequency measurements may reflect the size of fault-zone structures rather than earthquake stress drop or rupture velocity. This has wide-ranging implications for both interpretations of corner frequency and the cause of high-frequency earthquake ground motion damage. In the impact model, high-frequency ground motion results from the physics and time scale of elastic impacts rather than the initiation and completion of local fault slip. Thus, high-frequency energy can occur during the entire rupture process (e.g., Ishii et al., 2005) even without abrupt frictional changes and would not only be governed by the physics of friction-mediated local slip (Andrews, 1980;Dunham et al., 2011;Zielke et al., 2017).
Another implication is that when the number of impacts is large, high-frequency radiation patterns would be expected to be more isotropic (see Supporting Information), consistent with observed radiation patterns that are difficult to explain with standard interpretations (Somerville et al., 1997). Even when fault roughness (e.g., Dunham et al., 2011;Shi & Day, 2013) or heterogeneous structure (Graves & Pitarka, 2016) is accounted for, heterogeneous slip models still tend to either underpredict the degree to which the magnitudes of high-frequency fault-normal and fault-parallel ground motions are similar (e.g., Somerville et al., 1997) or overpredict the degree of isotropy for low-frequency data (Takemura et al., 2009). Thus, an additional physical mechanism to further increase high-frequency isotropy may be warranted. If elastic impacts dominate high-frequency seismic spectra, stress drops could still be constrained by spectra but such information would need to be extracted from an intermediate frequency band, making stress drop inferences more challenging (Figure 3b). For example, rise times T r and durations T would generally be longer than predicted by high-frequency corner frequency measurements, which are supported observationally for some events that are time resolved (Asano & Iwata, 2012;Duputel et al., 2013) and may be irresolvable without broadband instrumentation (Figure 3b).
The impact model also provides a structural context to explain variability in corner-frequency-based estimates of stress drop between mature (interplate and high slip rate) and immature (intraplate and low slip rate) fault zones (Allmann & Shearer, 2009;Goebel et al., 2017;Kanamori & Anderson, 1975). Mature fault zones, with more foliated, smoother, high aspect ratio structures would be expected to have larger T c for fault zones with similar widths (see Supporting Information; Brodsky et al., 2011). Similarly, the model explains the surprising observation that shallow earthquakes have some of the highest inferred stress drops (Abercrombie, 2014) despite the expectation that overall stresses are lower (Allmann & Shearer, 2009); within the present model, this can be a consequence of faults being rougher at shallower depths (Cochran et al., 2009;Finzi et al., 2009;Sylvester, 1988;Unsworth et al., 1997). Earthquake hazards would accordingly be greater near faults that are less foliated and rougher. In summary, the impact model provides an alternative explanation for the high-frequency corner frequency, which suggests that standard estimates of stress drop from corner frequency could be systematically inaccurate.
Given the encouraging observational support for the elastic impact model and the challenges that traditional frictional slip models have in explaining the full range of observations, we suggest it is worthwhile to make further observations (e.g., of particle sizes or rupture initiation) that would help confirm or reject the elastic impact model and the assumptions inherent in more traditional models.