资源环境科技发展态势分析平台

Mechanics and Statistics of Postseismic Shaking

https://doi.org/10.1029/2025GL116673
2025-09-12
查看原文
Geophysical Research Letters . Volume 52 , issue 18
T. Clements, E. S. Cochran, S. E. Minson, N. J. van der Elst, C. E. Yoon, A. Baltay, M. T. Page

Abstract

Analysis of 2 weeks of continuous post-seismic shaking after the 2019 M7.1 Ridgecrest, CA earthquake sequence using 4 nearby borehole seismometers reveals that continuous ground motions decay as Omori's law in time and follow the Gutenberg-Richter distribution in logarithmic amplitude. The measured temporal decay in amplitudes agrees with predictions of the rate-and-state framework and indicates shaking amplitudes are proportional to the velocity of afterslip. Our ground motion-based statistical framework provides a basis to forecast shaking intensity in the minutes to hours after a large earthquake.

Plain Language Summary

If an earthquake ruptures a fault and no one feels it, does it produce shaking? To answer this question, we measured seismic ground motion in the 2 weeks following the 2019 Ridgecrest earthquake in California. During this time, shaking was intermittently perceptible by humans, whereas nearby seismometers recorded continuous seismic shaking above pre-event noise levels. The amplitude of post-earthquake shaking decreased at a rate proportional to the inverse of time since the peak ground motion. This suggests that post-earthquake shaking is related to the gradual sliding that occurs along faults after an earthquake. Our findings provide a new method to forecast shaking in the minutes and hours following a major earthquake, giving communities and first responders valuable information about ongoing seismic hazard.

Key Points

  • Continuous seismic ground motion obeys the Omori and Gutenberg-Richter laws following a large earthquake

  • Shaking statistics using continuous seismograms do not require earthquake detection, location, or source parameter estimation

  • A rate and state model predicts a power-law decay in postseismic shaking with time

1 Introduction

An earthquake is almost always followed by postseismic shaking and when it is strong and damaging, noticeable shaking may be felt intermittently for days to weeks. When postseismic shaking is not reported, it is likely because the seismic source was either limited in strength, very deep, or far from the observer (Ōmori, 1894). During an aftershock sequence, it is desirable to quantify shaking at any particular point of observation, independent of the original seismic energy (Richter, 1935; Wald et al., 1999). At a particular site, combining measurements of the frequency of occurrence of shaking as a function of amplitude and the rate of decay of shaking with time could form the basis of a statistical shaking forecast.

Obtaining reliable forecasts of postseismic shaking in the minutes to hours after strong shaking is mostly limited by the reliance of current earthquake forecasts on earthquake catalogs, which have several limitations. First, earthquakes in nature are not point processes (Michael, 2005); they have spatial and temporal dimensions that can only be estimated once the event is over. Second, further increasing the spatial or temporal resolution of an earthquake catalog is limited by incompleteness of earthquake catalogs following a large earthquake, when strong shaking from the mainshock and large aftershocks obscures weaker shaking from smaller aftershocks (Kagan, 2004). Third, assumptions about attenuation, 3D structure, radiation pattern, the constancy of stress drops, and rupture propagation must be made to estimate the absolute location, timing, and magnitude of any particular earthquake (Brune, 1970; Thatcher & Hanks, 1973). Fourth, earthquake source parameters are inferred quantities, as opposed to directly measurable quantities, such as ground displacement, velocity, or acceleration. To address these issues, we construct a model and statistics of postseismic shaking using directly observed continuous seismic ground motion.

2 Distribution of Seismic Amplitudes

Usually the b $b$ -value in Gutenberg and Richter (1944)'s (GR) formula
f ( M ) = k 1 0 − b M ( k , b ; constant ) $f(M)=k1{0}^{-bM}\,(k,b;\text{constant})$ (1)
expressing the magnitude-frequency relation is determined by the maximum likelihood estimate using a catalog of earthquake magnitudes, M $M$ (Aki, 1965; Utsu, 1965). Recent developments in ground motion-based earthquake early warning (Kodera et al., 2018) and earthquake forecasting (Clements et al., 2024; Enescu et al., 2007; Lippiello et al., 2016, 2019), however, suggest an alternative method to forecast future shaking from continuous seismograms. Contrary to earthquake catalogs, continuous seismograms capture the full physics of earthquake rupture, record the summation of ground motion from multiple overlapping earthquakes, do not require assumptions about attenuation, and are available in real time. Recently, Sawazaki (2021, 2025) developed a statistical framework to predict the exceedance probability of aftershock ground motion using continuous seismic data from a single seismometer, and Shearer and Shakibay Senobari (2025) showed that the distribution of ground motion amplitudes recorded after the M7.1 Ridgecrest, California earthquake follow the Reasenberg and Jones (1989) model for aftershocks. Here, we extract time-varying statistics of an earthquake sequence from continuous seismic amplitude time series recorded at a single location. We follow Ishimoto and Iida (1939) (II), who showed that the distribution of peak ground motion amplitudes at a single location during an aftershock sequence follow a power law of the form
f ( a ) = k a − m $f(a)=k{a}^{-m}$ (2)
where a $a$ is the peak trace displacement, velocity, or acceleration of seismic ground motion for a single earthquake, f ( a ) $f(a)$ represents the number of earthquakes with amplitude between a $a$ and a + d a $a+da$ , and k $k$ and m $m$ are constants. The GR and II-laws have similar maximum likelihood estimates (Aki, 1965; Clauset et al., 2009; Utsu, 1965)
b = log 10 e M ‾ − M c $b=\frac{{\log }_{10}e}{\overline{M}-{M}_{c}}$ (3)
m = 1 + log 10 e log 10 a ‾ − log 10 a c $m=1+\frac{{\log }_{10}e}{\overline{{\log }_{10}a}-{\log }_{10}{a}_{c}}$ (4)
respectively, where log 10 a ‾ $\overline{{\log }_{10}a}$ and M ‾ $\overline{M}$ denote the average log amplitude and magnitude, respectively, M c ${M}_{c}$ is the magnitude of completeness, and a c ${a}_{c}$ is the minimum amplitude. If we use a simple ground motion model of the form
A = g 0 ( r ) + g M $A={g}_{0}(r)+gM$ (5)
where A = log 10 a $A={\log }_{10}a$ is the logarithm of the seismic amplitude, g 0 ${g}_{0}$ is a constant that depends on the distance r $r$ from a seismic source to an observer and g $g$ is a constant that relates amplitude and magnitude. Suzuki (1959), who assumed g = 1 $g=1$ , and later Utsu (1972) showed that,
m = b g + 1 . $m=\frac{b}{g}+1.$ (6)

The major implication of Equation 6 is that the b $b$ -value can be measured solely from the power-law distribution of ground motion amplitudes recorded at a single location (Page, 1968). Here, we assume g = 1 $g=1$ , following Richter (1935). We provide a derivation for the probabilistic inference of b $b$ -value from peak ground motion amplitudes in Text S1 in Supporting Information S1. A derivation for the b $b$ -value of a continuous seismogram, including body, surface, and coda waves, is beyond the scope of this paper.

3 Temporal Decay of Seismic Amplitudes

Temporal decays in post-seismic phenomena are ubiquitous: the rate of aftershocks (Ōmori, 1894; Utsu et al., 1995) and the velocity of afterslip (Ingleby & Wright, 2017; Perfettini & Avouac, 2004) decay like a power-law following a large earthquake. Here, we develop a mechanical model for the temporal decay in postseismic shaking, which we call “aftershaking,” and define as the continuous seismic ground motion recorded in the seconds to days following a large earthquake. To do so, we suppose that aftershaking and the two other main after-seismic quantities, aftershocks and afterslip, all originate from the same physical process: stick slip instability in the rate and state friction framework (Dieterich, 1979, 1981; Rice, 1983; Ruina, 1983). We start with the Perfettini and Avouac (2004); Perfettini and Avouac (2007) model, hereafter PA, in which afterslip in velocity-strengthening regions drives unruptured, locked regions to failure, causing aftershocks. PA showed that the slip, δ ( t ) $\delta (t)$ , and velocity, V ( t ) $V(t)$ , on a fault with a velocity-strengthening rheology and rate-dependent friction law following a shear stress increase Δ τ ( t ) ${\Delta }\tau (t)$ , are given by
δ ( t ) = V L t r log 1 + V + V L exp t t r − 1 $\delta (t)={V}_{L}{t}_{r}\,\log \left[1+\frac{{V}_{+}}{{V}_{L}}\left(\mathrm{exp}\left(\frac{t}{{t}_{r}}\right)-1\right)\right]$ (7)
V ( t ) = V + exp t t r 1 + V + V L exp t t r − 1 $V(t)={V}_{+}\frac{\mathrm{exp}\left(\frac{t}{{t}_{r}}\right)}{1+\frac{{V}_{+}}{{V}_{L}}\left(\mathrm{exp}\left(\frac{t}{{t}_{r}}\right)-1\right)}$ (8)
where V + ${V}_{+}$ and V L ${V}_{L}$ are the sliding velocities before and after the shear stress step, respectively, V + / V L = exp Δ τ a r s σ ${V}_{+}/{V}_{L}=\mathrm{exp}\left(\frac{{\Delta }\tau }{{a}_{rs}\sigma }\right)$ , a r s ${a}_{rs}$ is the rate constant, σ n ${\sigma }_{n}$ is the effective normal stress, and t r = a r s σ n τ ̇ ${t}_{r}=\frac{{a}_{rs}{\sigma }_{n}}{\dot{\tau }}$ is a timescale representing the period of postseismic relaxation, where τ ̇ $\dot{\tau }$ is the background stressing rate. The PA model explains both the logarithmic increase in surface displacement with time due to afterslip (Marone et al., 1991) and the Omori-like decay in seismicity rates (Dieterich, 1994) following a large earthquake. For times t ≪ t r $t\ll {t}_{r}$ , Equation 8 behaves like a power law, similar in form to Ōmori (1894)'s law.
How then, does the sliding velocity on the fault relate to the seismic amplitude on the surface? Let us begin with a thought experiment with two seismic sensors. The first is buried on the fault at depth and the second is placed on the surface in the far-field, more than a few wavelengths from the fault. For the first sensor, the sliding velocity during rupture or afterslip is its seismic velocity. For the far field sensor on the surface, its seismic velocity is proportional to the slip velocity on the fault through the representation theorem (Aki & Richards, 2002). Theoretical and empirical evidence strongly support this proportionality: (a) dynamic rupture models have a linear relation between rupture velocity and slip velocity (Weng & Ampuero, 2022); (b) there is a linear relation between rupture velocity and peak ground velocity in near-fault regions (McGarr & Fletcher, 2007); and (c) laboratory friction experiments find that the slip velocity is proportional to acoustic emission amplitudes (Bolton et al., 2022). Building on these lines of evidence, we substitute the PA model for sliding velocity, V ( t ) $V(t)$ , with a seismic amplitude on the fault, a ( t ) $a(t)$ , that is proportional to the sliding velocity, a ( t ) ∝ V ( t ) $a(t)\propto V(t)$ , which gives,
a ( t ) = a + exp t t r 1 + a + a L exp t t r − 1 . $a(t)={a}_{+}\frac{\mathrm{exp}\left(\frac{t}{{t}_{r}}\right)}{1+\frac{{a}_{+}}{{a}_{L}}\left(\mathrm{exp}\left(\frac{t}{{t}_{r}}\right)-1\right)}.$ (9)
where a + ${a}_{+}$ and a L ${a}_{L}$ are the seismic amplitudes on the fault before and after the shear stress step. The seismic amplitude observed on the surface at a distance r $r$ from a slip patch at depth, a ( r , t ) $a(r,t)$ will be modified by site and path effects
a ( r , t ) = S ⋅ a ( t ) ⋅ e − r r , $a(r,t)=S\cdot a(t)\cdot \frac{{e}^{-r}}{r},$ (10)
where S $S$ accounts for site effects, a ( t ) $a(t)$ is the seismic amplitude on the fault at depth, and e − r r $\frac{{e}^{-r}}{r}$ accounts for geometric spreading and anelastic attenuation. Similar to studies of aftershock rate and afterslip, we take a ( t ) $a(t)$ as the time-average amplitude of seismic waves generated by slip velocity V ( t ) $V(t)$ over a reasonably short time, for example, minutes, rather than the instantaneous seismic amplitude. For times t ≪ t r $t\ll {t}_{r}$ , Equation 10 reduces to
a ( r , t ) = K ( r ) c + t $a(r,t)=\frac{K(r)}{c+t}$ (11)
where c = a L t r a + = a r s σ n τ ̇ exp − Δ τ / a r s σ n $c=\frac{{a}_{L}{t}_{r}}{{a}_{+}}=\left(\frac{{a}_{rs}{\sigma }_{n}}{\dot{\tau }}\right)\mathrm{exp}\left(-{\Delta }\tau /{a}_{rs}{\sigma }_{n}\right)$ and K = S ⋅ a L ⋅ t r e − r r = S ⋅ a r s σ n τ ̇ a L e − r r $K=S\cdot {a}_{L}\cdot {t}_{r}\frac{{e}^{-r}}{r}=S\cdot \left(\frac{{a}_{rs}{\sigma }_{n}}{\dot{\tau }}\right){a}_{L}\frac{{e}^{-r}}{r}$ (derivation in Text S2 in Supporting Information S1). Equation 11 has a similar form to the PA and Dieterich formulations for afterslip and seismicity rate, respectively, and the original (Ōmori, 1894) and modified (Utsu et al., 1995) Omori laws for the decay of aftershock rate with time with two notable exceptions: (a) Equation 11 describes ground motion everywhere on the surface, rather than on the fault at depth and (b) our formulation decays as t − 1 ${t}^{-1}$ , rather than t − p ${t}^{-p}$ . Regarding the exponent of the power-law fall-off, Sawazaki (2021) found that the maximum amplitude of a continuous seismogram after a large earthquake follows a power-law decay with time as t − p g b ${t}^{\tfrac{-pg}{b}}$ , where p $p$ , g $g$ , and b $b$ are all constants about equal to one.

4 Measuring Temporal Ground Motion Statistics

Seismic ground motion following a large earthquake is essentially a convolution of finite duration, band-limited, white Gaussian noise signals (Hanks & McGuire, 1981). This has led to difficulty in detecting individual seismic phases in the subsequent seconds to hours after the mainshock (Enescu et al., 2007; Kagan, 2004; Peng & Zhao, 2009). This so-called “short-term incompleteness” suggests a statistical approach to postseismic forecasting that fits the combined distribution of continuous seismic ground motion and noise, rather than the peak ground motion from distinct aftershocks. To do so, we repurpose Ogata and Katsura (1993); Ogata and Katsura (2006)'s full-magnitude range model, which combines the GR-law with a detection rate function, with a formulation based on logarithmic ground motion amplitudes. Our full-amplitude model has two parts: (a) a detection function which assumes that seismic noise amplitudes above 1 Hz over short time-periods follow a log-normal distribution (Figures 1d and 1g) (Ringdal & Bungum, 1977) and (b) a seismic amplitude distribution that follows the II-law (Figures 1c and 1f) in linear amplitude or the GR-law in log-amplitude with the assumption that g = 1 $g=1$ in Equation 5. The detection rate function, q ( A | μ , σ ) $q(A\vert \mu ,\sigma )$ , for seismic amplitude thus follows a cumulative log-normal distribution (Ringdal, 1975),
q ( A | μ , σ ) = ∫ − ∞ A 1 2 π σ e − ( x − μ ) 2 2 σ 2 d x , $q(A\vert \mu ,\sigma )=\int \nolimits_{-\infty }^{A}\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\tfrac{{(x-\mu )}^{2}}{2{\sigma }^{2}}}dx,$ (12)
where μ $\mu $ is the logarithmic seismic amplitude, A $A$ , that has a 50 % $\%$ detection rate and σ $\sigma $ controls the amplitude range over which full detection occurs. We define A c = μ + σ ${A}_{c}=\mu +\sigma $ as the “amplitude of completeness,” which in this case is the 8 4 t h $8{4}^{th}$ percentile of logarithmic ground motion within a particular time window, T $T$ . Similar to the time-varying magnitude of completeness (Hainzl, 2016), A C ( t ) ${A}_{C}(t)$ measures the minimum seismic noise level generated by overlapping earthquakes and decays through time following a large earthquake (Shearer & Shakibay Senobari, 2025). A C ${A}_{C}$ is equivalent to the seismicity rate in the Dieterich model, the afterslip velocity in the PA model, and the log of amplitude a ( t ) $a(t)$ in Equation 11: high seismicity rates generate measurably higher A C ${A}_{C}$ . Within a single time window the distribution of logarithmic seismic amplitudes above A C ${A}_{C}$ follows the GR distribution, that is, G R ( A ) = k 1 0 − b A − A C $GR(A)=k1{0}^{-b\left(A-{A}_{C}\right)}$ , whereas over many time windows A C ${A}_{C}$ decays like Omori's law. The product G R ( A ) q ( A | μ , σ ) $GR(A)q(A\vert \mu ,\sigma )$ matches histograms of continuous seismic amplitudes well, as shown in Figures 1e and 1h. We fit all ground motion time series amplitudes recorded by a single sensor to the probability density
f A | b A , μ , σ = 1 0 − b A q ( A | μ , σ ) ∫ − ∞ ∞ 1 0 − b A q ( A | μ , σ ) d A = β e − β ( A − μ ) − β 2 σ 2 / 2 q ( A | μ , σ ) \begin{align*}\hfill f\left(A\vert {b}_{A},\mu ,\sigma \right)& =\frac{1{0}^{-bA}q(A\vert \mu ,\sigma )}{\int \nolimits_{-\infty }^{\infty }1{0}^{-bA}q(A\vert \mu ,\sigma )dA}\hfill \\ \hfill & =\beta {e}^{-\beta (A-\mu )-{\beta }^{2}{\sigma }^{2}/2}q(A\vert \mu ,\sigma )\hfill \end{align*} (13)
where β = b A ln 10 $\beta ={b}_{A}\,\ln\,10$ and b A ${b}_{A}$ is the amplitude-based b $b$ -value, or the falloff in logarithmic ground motion amplitudes. We fit the parameters β , μ , a n d σ $\beta ,\mu ,\mathrm{a}\mathrm{n}\mathrm{d}\,\sigma $ to seismic amplitude data at a particular station by maximizing the log likelihood function
ln L b A , μ , σ = ∑ 0 < t i < T ln f A t i | b A , μ , σ $\ln\,L\left({b}_{A},\mu ,\sigma \right)=\sum\limits _{0< {t}_{i}< T}\ln\,f\left(A\left({t}_{i}\right)\vert {b}_{A},\mu ,\sigma \right)$ (14)
where t i ${t}_{i}$ are the sample times of all amplitude samples in the time interval ( 0 , T ) $(0,T)$ (Ogata & Katsura, 2006), using a constrained optimization algorithm.
Details are in the caption following the image
Figure 1
Open in figure viewerPowerPoint

A summary of the aftershaking statistical model. (a) Horizontal component velocity seismogram for station PB.B921 for 1–10 min after the 7/6/2019 M7.1 Ridgecrest earthquake. (b) Envelope of the 3-component velocity seismogram for station PB.B921 for same time period as in (a). (c) Power law and (f) straight line indicate the Ishimoto-Iida relation in linear and logarithmic scales, respectively. (d) Detection rate (solid black line) of seismic ground motion as a function of linear and (g) logarithmic seismic amplitude. Dashed gray lines indicates μ $\mu $ , the amplitude with 50 % $50\%$ detection rate. Dotted gray line indicates σ $\sigma $ , the amplitude range over which full detection occurs. (e) Full-amplitude model, f A | b A , μ , σ $f\left(A\vert {b}_{A},\mu ,\sigma \right)$ (black line), fit to amplitude time series in (b). Blue dots in (e) and (h) represent duration of amplitudes in (b) in logarithmic amplitude bins in linear and logarithmic time, respectively. Amplitudes are shown in bins for display purposes - model is fit to continuous amplitudes. Vertical solid gray line indicates the amplitude of completeness, A C ${A}_{C}$ . (c)–(h) modified from Ogata and Katsura (2006) Figure 2.

5 Application to Continuous Earthquake Ground Motion

Here, we apply our full-amplitude range method to measure shaking statistics during the 2019 Ridgecrest, CA earthquake sequence. The sequence began with a M6.4 foreshock on July 4, followed by a M7.1 mainshock 34 hr later on July 6. Over 100,000 recorded aftershocks (Ross et al., 2019), extensive field mapping of surface offsets (DuRoss et al., 2020), and multiple satellite geodetic measurements (Barnhart et al., 2019) illuminated a complex system of right-lateral northwest-striking and left-lateral southwest-striking faults. The Ridgecrest sequence is an interesting test case of our method for two reasons: (a) drops in GR b $b$ -value were retrospectively measured in the area and time preceding the M7.1 earthquake (Dascher-Cousineau et al., 2020; Gulia et al., 2020; Huang et al., 2020; Nanjo, 2020), and (b) the sequence was one of the best-recorded continental earthquakes to date. To test our full-amplitude method, we use continuous velocity waveforms recorded by four borehole seismometers (2 Hz geophones) in the Plate Boundary Observatory network (Mattioli et al., 2020) in the 3 days preceding and 14 days following the M6.4 earthquake on July 4 as shown in Figure 2a. The stations are located within 5–30 km of the mainshock surface rupture and at depths of 140–190 m below the surface. We apply minimal processing to day-long, raw seismic waveforms with the goal of preserving seismic amplitudes while removing ambient noise signals and instrumental irregularities: we detrend, apply a 20-s cosine taper, carefully remove time gaps and clipped sections (one station clipped during the M6.4 earthquake and two during the M7.1 earthquake), apply a bandpass filter between 2 and 10 Hz, downsample to 20 Hz, remove the instrument gain, and convert 3-component time series to a 1D seismic amplitude envelope as,
A ( t ) = log 10 E ( t ) 2 + N ( t ) 2 + Z ( t ) 2 , $A(t)={\log }_{10}\left(\sqrt{E{(t)}^{2}+N{(t)}^{2}+Z{(t)}^{2}}\right),$ (15)
where E ( t ) $E(t)$ , N ( t ) $N(t)$ , and Z ( t ) $Z(t)$ are the east, north, and vertical component time series, respectively.
Details are in the caption following the image
Figure 2
Open in figure viewerPowerPoint

(a) Map showing location of four borehole seismic stations (triangles) and location of Ridgecrest earthquake locations from Quake Template Matching (QTM) Catalog (Ross et al., 2019). Dashed circle indicates spatial extent of the QTM catalog. (b) Power-law decay in the amplitude of completeness, A C ${A}_{C}$ , following the 2019 M7.1 Ridgecrest earthquake. Time is measured relative to P-wave arrival at each station and the first window is 13 s long. Colored lines indicate A C ${A}_{C}$ for each station in panel (a) measured in expanding exponential time windows. Black and gray dashed lines are fit of Equation 11 and K t − p $K{t}^{-p}$ to A C ${A}_{C}$ averaged over stations, respectively.

Our model predicts a power-law t − 1 ${t}^{-1}$ decay in ground motion amplitudes following a large earthquake. To measure this decay, we calculate the minimal resolvable seismic ground motion, A C ${A}_{C}$ , from continuous seismograms immediately following the M7.1 earthquake. Starting 50 s after the time of peak ground motion, when aftershock ground motion becomes identifiable from coda waves, we sample A C ${A}_{C}$ in exponentially expanding time windows, from 13 s long initially to 3 days long 2 weeks later. We bootstrap sample each window 1,000 times and use 10,000 amplitude samples for windows longer than 500 s for computational expediency. As shown in Figure 2b, A C ${A}_{C}$ at each station follows a power law decay with time in accordance with Omori's law. Our fit for the falloff in A C ${A}_{C}$ using Equation 11 gives c = 0 s $c=0\hspace*{.5em}\mathrm{s}$ when constraining c > 0 $c > 0$ or c = − 24 s $c=-24\hspace*{.5em}\mathrm{s}$ when allowing for c < 0 $c< 0$ . We note, though, that the falloff in A C ${A}_{C}$ for this particular sequence is better fit for later times with functional form A ( t ) = K t p $A(t)=\frac{K}{{t}^{p}}$ , with p = 1.12 $p=1.12$ , which is close to the p = 1.1 $p=1.1$ of the most recent U.S. Geological Survey operational earthquake forecast for the 2019 Ridgecrest, CA earthquake sequence (U.S. Geological Survey, 2025). Directly measuring the c $c$ -value from continuous waveforms is likely impossible, though, given the influence of surface and coda waves in the tens of seconds following peak ground motion at a particular site (Enescu et al., 2009; Kagan, 2004; Kagan & Houston, 2005; Peng et al., 2006, 2007).

To assess temporal trends in shaking statistics, we calculate A C ${A}_{C}$ and b A ${b}_{A}$ values using moving 4-hr windows of ground motion throughout the entire sequence. For comparison, we calculate b $b$ -values and time-dependent magnitude of completeness, M c ( t ) ${M}_{c}(t)$ , using the high-resolution earthquake catalog of Ross et al. (2019) with the β + ${\beta }^{+}$ (van der Elst, 2021) and maximum curvature (Wiemer & Wyss, 2000) techniques, respectively. For β + ${\beta }^{+}$ estimates, we round magnitudes to the nearest 0.1 magnitude units, use a 900-event moving window (equivalent to 3.9 hr event windows, on average), and impose a positive magnitude difference cutoff of d M c = + 0.5 $d{M}_{c}=+0.5$ , while we define the magnitude of completeness as the mode of the distribution of rounded magnitudes over a 400-event window plus a 0.2 magnitude unit correction following van der Elst (2021). These parameter choices mainly affect absolute, rather than relative, estimates of the b $b$ -value (Figure S1 in Supporting Information S1). As shown in Figure 3a, M c ( t ) ${M}_{c}(t)$ and A C ( t ) ${A}_{C}(t)$ closely track each other through time with a distance-dependent offset. Among the four stations, A C ${A}_{C}$ falls off with distance from the fault, as expected by distance-dependent attenuation. Interestingly, σ $\sigma $ , the range over which amplitudes are fully detected, is similar across stations in the 1–2 days following the M7.1, then decays faster in time at further away stations (Figure S2 in Supporting Information S1). Prior to the M6.4 foreshock, fitting Equation 13 measures the mode of the background noise level of ≈ 5 × 1 0 − 6 c m / s ${\approx} 5\times 1{0}^{-6}\mathrm{c}\mathrm{m}/\mathrm{s}$ , rather than the minimum seismic ground motion generated by background seismicity. Establishing a pre-event A C ${A}_{C}$ is important for placing a lower bound on amplitude steps. Distinct steps in A C ${A}_{C}$ and M c ${M}_{c}$ occur following three earthquakes in the Ridgecrest sequence: M6.4 on July 4, M5.4 on July 5, and M7.1 on July 6. The M7.1 had a minimum of 4 orders of magnitude increase in A C ${A}_{C}$ above the background noise level, while the steps in A C ${A}_{C}$ and M c ${M}_{c}$ are one and two orders of magnitude lower following the M6.4 and 5.4, respectively, than after the M7.1. In the Dieterich and PA models, the increase in seismicity rate, or jump in A C ${A}_{C}$ in our case, is directly proportional to the stress step, Δ τ ${\Delta }\tau $ , whereas the aftershock duration is independent of earthquake magnitude or Δ τ ${\Delta }\tau $ . A C ${A}_{C}$ approaches the background level at the furthest station (B916) after 2 weeks. Aftershocks of the M7.1 did not noticeably raise A C ${A}_{C}$ above the expected t − 1 ${t}^{-1}$ decay. b $b$ -values computed with β + ${\beta }^{+}$ and b A ${b}_{A}$ show similar dynamics through the sequence (Figure 3b). b $b$ -values are lowest immediately following the M6.4 but then rise before the M7.1, as noted by Gulia et al. (2020). Low b A ${b}_{A}$ values in the 4 hr following the M6.4, M5.4, and M7.1 earthquakes are related to the difficulty in applying the full-amplitude model for large jumps in seismic amplitude (see Figure S3 and Text S3 in Supporting Information S1 for details). Even though β + ${\beta }^{+}$ accounts for changing completeness, b A ${b}_{A}$ time series are more correlated to β + ${\beta }^{+}$ in the 10–20 days following the M7.1 than in the first 10 days (Figure S4 in Supporting Information S1). This difference in b $b$ -values suggests the Quake Template Matching (QTM) catalog is missing earthquakes early on in the sequence when the seismicity rate is highest (van der Elst & Page, 2023) and achieves full detection as the earthquake rate decays (Figure S5 in Supporting Information S1). b A ${b}_{A}$ time series oscillate around a value of 0.9 for the 2 weeks following the M7.1, then rise above b = 1 $b=1$ for two furthest stations when A C ${A}_{C}$ approaches its pre-event level. The main difference between catalog-based b $b$ -values and b A ${b}_{A}$ is that b A ${b}_{A}$ is a local and surficial b $b$ -value. b A ${b}_{A}$ is sensitive to nearby seismic sources and easiest to constrain when the rate of seismic shaking is high.

Details are in the caption following the image
Figure 3
Open in figure viewerPowerPoint

Comparison of shaking statistics to earthquake catalog statistics for the 2019 Ridgecrest, CA earthquake sequence. (a) Fall-off in amplitude of completeness, A C ${A}_{C}$ (colored lines), and magnitude of completeness, M c ${M}_{c}$ (black line). Magnitude of completeness is calculated using the Quake Template Matching (QTM) catalog (Ross et al., 2019). Vertical gray bars indicate influence of M6.4, M5.4, and M7.1 earthquakes, respectively. Horizontal dashed line indicates pre-sequence median A C ${A}_{C}$ . (b) Comparison of amplitude-based b $b$ -value, b A ${b}_{A}$ (colored lines), at each station and b $b$ -value computed with β + ${\beta }^{+}$ method using QTM catalog.

6 Discussion

Gutenberg-Richter (Ishimoto-Iida) and Omori scaling laws can be extracted from the statistics of continuous ground motion. This suggests the possibility of monitoring the b $b$ -value and temporal decay of an on-going earthquake sequence in real time. In that case, amplitude forecasting with continuous data could bridge the fields of earthquake early warning and earthquake forecasting in the minutes to hours following a large earthquake. At a forecasting location, aftershaking will obey a power-law decay in time and amplitude distribution. The amplitude distribution will follow the Ishimoto-Iida relation in linear amplitude or the Gutenberg-Richter distribution in logarithmic amplitude, which through Equation 5 gives an equivalence between the GR and amplitude-based b $b$ -value, the magnitude and logarithmic seismic amplitude, and the magnitude of completeness and the amplitude of completeness. Using continuous seismograms, we infer that the c $c$ -value in the modified-Omori law is close to zero. This suggests that velocity-strengthening driven afterslip begins immediately after the end of rupture. By extension, aftershock times should be measured with respect to the end of rupture, rather than the origin time.

A shaking-based approach to earthquake statistics has several limitations. First, seismometers must be recording near the fault or in an area of interest before a seismic sequence. Second, processing continuous waveform data, rather than peak amplitudes, is computationally expensive. Third, anthropogenic and environmental seismic noise mask the minimum resolvable seismic ground motion, A C ${A}_{C}$ , earlier in an earthquake sequence than the peak seismic amplitudes that are used to create earthquake catalogs. These issues could be ameliorated with the judicious placement of high-quality surface or borehole stations near faults of interest and methodological improvements. Additionally, shaking statistics do not require earthquake detection, location, or source parameter estimation, which avoids the potential error propagation that comes from estimating intermediate quantities. Finally, since shaking-based statistics could be derived from streaming waveforms in near real time, they have the potential to be available well before high-resolution catalogs, allowing for sooner forecasts of ground motion.

Our results can be generalized to the following statements: (a) the minimum resolvable ground motion at the surface is proportional to the afterslip velocity at depth, (b) linear and logarithmic seismic ground motion amplitudes obey the Omori and Gutenberg-Richter laws, respectively, following a large earthquake, and (c) a seismogram is a continuous earthquake catalog and can be utilized in that context.

Acknowledgments

TC was partially supported by the U.S. Geological Survey Mendenhall fellowship program. The authors thank Jeanne Hardebeck and two anonymous reviewers for reviews that improved the manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

    Conflict of Interest

    The authors declare no conflicts of interest relevant to this study.

    Data Availability Statement

    Earthquake waveforms were downloaded from the Southern California Earthquake Data Center (SCEDC, 2013). The QTM earthquake catalog for the Ridgecrest sequence (Ross et al., 2019) was downloaded from the Southern California Earthquake Data Center (https://scedc.caltech.edu/data/qtm-ridgecrest.html). The figures were created using matplotlib (Hunter, 2007) and PyGMT (Uieda et al., 2021) software.

      相关资源