An important seismological question is whether detectable preparatory phases exist and can be observed before major earthquakes. In case of a positive response, the next question would be whether the same preparatory phase is common to all mainshocks, and many recent studies went deep into these issues. For example, Bletery and Nocquet (2023) analyzed the very final moments of the seismic cycle on a set of 90 earthquakes (M7); they found that a ∼2-hour-long exponential acceleration characterizes the slip before the ruptures, suggesting that large earthquakes start with a precursory phase of slip. Other recent laboratory observations focused on the nucleation process: for example, Bolton et al. (2023) pointed to key differences in the nucleation process of slow and fast labquakes, showing that the evolution in space and time of the foreshock activity is linked to fault slip velocity, similar to the increase and coalescence of foreshocks/seismic activity that occur prior to large crustal earthquakes (Dodge et al., 1996; Bouchon et al., 2011; Kato et al., 2016; 2012; Chen and Shearer, 2013; Brodsky and Lay, 2014; Sugan et al., 2014; Ellsworth and Bulut, 2018; Yoon et al., 2019; Ben-Zion and Zaliapin, 2020; Cattania and Segall, 2021).

Another important study on laboratory earthquakes is the one by Vasseur et al. (2017), who documented systematic evidence for a detectable preparatory phase to failure in rock samples during laboratory experiments, showing that an accurate prediction of the time of occurrence of the catastrophic failure of the sample was possible in most cases, and that its accuracy depended on the typical inter-flaw distance of the material that was being probed. They showed that there are materials with few defects (e.g., glass) for which the time of occurrence of the catastrophic failure is inherently unpredictable. This is because in a sample with a dense distribution of defects, cracks can connect adjacent defects when the material is under stress, and the process of catastrophic failure is at the end of a somewhat smooth coalescence of a myriad of individual cracks. A very homogeneous material, on the contrary, is characterized by a small number of defects, and for a crack to connect two distant flaws the sample needs to undergo a large stress that may easily surpass its very strength.

But how can we go from the centimetric laboratory scale to the tens or the hundreds of kilometers of a large crustal mainshock? Assuming that such a scaling operation is feasible, during the latter portion of the seismic cycle we expect that crustal rocks are subjected to a dramatic process of coalescence of suitably oriented cracks within the volume that immediately surrounds the fault plane. As a consequence, we expect substantial changes in crack density, with spatial and temporal fluctuations of the population of cracks of different dimensions, which would significantly impact the permeability of rocks, the anelastic attenuation of seismic waves, and possibly any deep non-volcanic tremors in an area (NVT), when they exist (Nadeau and Dolenc, 2005; Nadeau and Guilhem, 2009).

Although we think that detecting the exponential slip acceleration observed by Bletery and Nocquet (2023) would be very useful in practice, monitoring the evolution of other physical parameters during the entire earthquake cycle, in absence of slip, could provide a deeper understanding of the evolution of the stress on the fault, and, more generally, on the physics of earthquakes and faulting. One of such parameters of interest is the non-geometric attenuation: Q S − 1 t , f .

If we want to use anelastic attenuation to monitor stress changes and crack coalescence within crustal rocks, we need to understand the role of crustal fluids in seismic attenuation. In fact, pore fluids in crustal rocks define the anelastic attenuation of seismic waves (O’Connell and Budiansky, 1977), and more importantly they may dominate the evolution of the strength of a fault during the latest portion of its seismic cycle through diffusion (Sibson, 2009; Lucente et al., 2010; Scholtz, 2019; Malagnini et al., 2022). The variability of anelastic attenuation over time is dominated by the fluctuations of rock’s bulk permeability (Malagnini and Parsons, 2020), and can be used to study interesting aspects of single- and multi-mainshock seismic sequences (Malagnini et al., 2022; Lucente et al., 2010; D’Amico et al., 2010).

Previous work by Malagnini et al. (2019) and Malagnini and Parsons (2020) demonstrated that fluctuations of non-geometric crustal attenuation (not including site-related effects) are intimately related to crustal rocks’ bulk permeability, and thus to the state of stress they undergo. The cited authors explored the following causal relationships: (i) between perturbations of stress (tectonic, seasonal, tidal) and perturbations of Q S − 1 , that is: δ Q S − 1 t , f , and (ii) between shaking-induced changes in rock bulk permeability and perturbations δ Q S − 1 t , f . Such damage and healing processe were noted by Kelly et al. (2013), and processes of removal/redeposition of colloidal deposits in rock pores and cracks may also be influential (Roeloffs, 1998; Liu and Manga, 2009).

About relationships of the kind (i) listed above, the modeling effort undertaken by Johnson et al. (2017) about hydrospheric, thermal, and tidal loading cycles on the faults of California allowed Malagnini and co-workers to interpret the specific spectral peaks of the attenuation time histories, and, more importantly, to indirectly estimate the extreme sensitivity of seismic attenuation to stress changes. About relationships of the second kind (ii), their interpretation was based on the studies by Brodsky, Manga, and co-workers (e.g., Brodsky et al., 2003; Manga and Brodsky, 2006; Liu and Manga, 2009; Manga et al., 2012), and were also backed by results of numerical experiments (Barbosa et al., 2019).

Our hypothesis is that the changes in rock bulk permeability that take place during the pre-critical stage of the seismic cycle (dilatancy) must have detectable signatures on the anelastic contribution to the total seismic attenuation. It is conceivable that the increased permeability along the very fault plane under subcritical conditions (Townend and Zoback, 2000) allows pore pressure to diminish immediately below the fault, at the upper limit of NVT depths (in Central California, tremor activity characterizes a planar patch at depths between 10 and 40 km, see Guilhem and Nadeau, 2012). Pore pressure diminution at NVT depths would correspond to a decrease in NVT, and to an upward pore fluid migration that would take place along the fault plane. In turn, the latter could weaken the fault segment studied here, and facilitate the occurrence of a mainshock. In what follows we will try to substantiate this hypothesis by looking for correlations between Q S − 1 and NVT activity.

When the fault belongs to a system like the one responsible for the 2016 sequence of Amatrice-Visso-Norcia (Central Italy), a major issue about recognizing the signature of criticality on seismic attenuation is represented by its geometric complexity. In fact, multiple individual structures of such a system could go critical one after the other, and fail in a cascade of multiple mainshocks, each one with its own preparatory phase, shaking-induced damage, and healing phase, with all these phenomena partially overlapping in time.

Our understanding of the phenomena described in this study is based on the hypothesis that the fluctuations δ Q S − 1 t , f , and those of deep NVT activity, provide information about the state of stress of the crustal rocks that are in the immediate vicinity of the fault surface and for the NVT deep below the seismogenic zone, especially about how they evolve over time. The correlation between 1/Q and the stress acting on the fault was investigated in detail by Malagnini and Parsons (2020), and by Malagnini et al. (2022), who related the variability of Q S − 1 to that of the bulk permeability of crustal rocks.

Within the limits of our hypothesis, a preparatory phase for the Parkfield mainshocks could be detected based on the observation of changing NVT activity and of some monotonic variations of the attenuation parameter. The idea is that the crack population within crustal rocks has a specific (yet unknown) statistical distribution of dimensions, and that its characteristics vary as a function of time following the fluctuations of the stress acting upon the investigated volume.

In an initial stage in which the crustal rocks gradually undergo an increasing deviatoric stress from an unstressed initial state, the crack density is expected to grow (e.g., Scholz, 2019), the rock bulk permeability is expected to increase, and the same is true for the attenuation parameter. In a later stage of the seismic cycle, when the fault reaches a subcritical stress level, crack coalescence becomes the dominant phenomenon: some cracks would grow and coalesce, forming the structure where the catastrophic fracture will eventually take place. The phenomenon of fracture coalescence was studied by Bolton et al. (2023), who recorded acoustic emissions (AE) in lab samples under stress and concluded that a transition from pervasive to localized deformation must occur within the sample in the lead-up to failure.

In general, the strain localization produced by crack coalescence would relax the deviatoric stress away from the volume in which coalescence takes place, so that the volume surrounding the growing fracture would go into compaction, with a substantial reduction of crack density (Scholz, 2019). It is conceivable that the unknown crack distribution would change dramatically during the whole process of crack coalescence and growth of the macroscopic fracture, and such changes would be drastic for the portion of the crack population characterized by some preferred lengths and orientations, with the longer cracks contributing the most to the process. It is conceivable that these processes are linked to the activity rate changes observed by the NVT, at least partially.

The clear bifurcations shown by δ Q S − 1 t , f in two frequency bands in Figures 5A, B right before the Parkfield manshock tell us that the crack population within different frequency bands reacts in two opposite ways to the stress reaching a subcritical level. In fact, after a 7 months-long time window, started on 22 December 2003, of smooth monotonic increase of δ Q S − 1 t , f , the 30–50 Hz stack shows a decreasing attenuation (starting in the month of July 2004), whereas the 2–5 Hz stack shows an opposite trend of increasing attenuation. In both cases, the July-September 2004 fluctuations that lead to the bifurcations of Figures 5A, B are much larger than the standard errors of the individual data points of δ Q S − 1 t , f . We point out that since our dataset can only yield information about the crustal volume in the immediate vicinity to the fault plane, our estimates of δ Q S − 1 t , f may primarily be used as a proxy to fault strength.

Below 5 Hz the wavefield sampled in this study is dominated by surface waves (see Supplementary Figure S3). This makes the bifurcation even more interesting, because long wavelengths from surface waves (0.2–1.5 km) do not sample a deeper portion of the crust than the 30–50 Hz short wavelengths (20–60 m) from direct S-waves. In fact, it is just the opposite, and the increased attenuation in the 30–50 Hz band is partially influenced by the deeper portion of the SAF, where the nucleation of the 2004 earthquake took place, and just above the wide area where the NVT activity takes place. In fact, the total attenuation is the integral of the local attenuation experienced along the entire travel path. If the increased attenuation for low-frequency seismic waves is the result of crack opening and interconnection related to strain localization through the crack coalescence phenomenon, the decreasing attenuation in the high-frequency band is likely the result of the closing cracks in the 20–60 m length range, as the stress localizes on the growing (will-be) catastrophic fracture and relaxes elsewhere.

Likely, the bifurcations seen in Figures 5A, B between the low- and the high-frequency estimates of the stacked attenuation on both sides of the SAF mark the transition of the Parkfield’s asperity to criticality, when the coalescence of cracks gradually evolves into a catastrophic fracture (Scholz, 2019). The 2–5 Hz band describes the attenuation at wavelengths between 0.2 and 1.5 km, which likely represent the minimum crack lengths participating in the coalescence that takes place along the SAF.

From Figures 5A, B we see that the splitting between the low- and the high-frequencies is common to both sides of the SAF in the 6–8 weeks before the Parkfield mainshock. Moreover, the estimates of the attenuation parameter on the Pacific side of the SAF in the two frequency bands (2–5 and 30–50 Hz) fluctuate synchronously, whereas on the North American side of the SAF the two frequency bands show a more irregular behavior, with asynchronous fluctuations. These differences are the results of the different fluctuations of two subsets of cracks with specific average lengths in the general crack population and are likely due to the structural differences existing between the crust on the two sides of the SAF. Finally, a synchronous behavior characterizes both sides of the SAF in the two intermediate frequency bands of 5–10 and 12–26 Hz, where the bifurcations are not observed (Supplementary Figure S1), indicating that all the action (coalescence/compaction due to stress localization/relaxation) is confined within the two extreme bands of the available spectrum.

Another interesting feature is that the San Simeon earthquake produces a change in all the time histories of the stacked δ Q S − 1 t , f , regardless of the presence of bifurcations: they all suddenly become smooth, with a tendency for a monotonic trend of increasing attenuation parameter. Our hypothesis is that the time window that starts with the occurrence of the San Simeon earthquake and ends at the onset of bifurcation (July 2004) contains the first part of the crack coalescence process, where the drastic reduction of δ Q S − 1 t , f observed on both sides of the SAF is the signature of stress localization in a very narrow volume along the fault surface. We argue that the pre-critical state is marked by the onset of the bifurcation: it starts when the crack coalescence (stress localization) is strong enough to allow compaction away from the growing fracture, where the stress that acts within the narrow volume surrounding the macroscopic fracture can relax.

About the differences between the time histories δ Q S − 1 t , f on the two sides of the SAF, we point out that the SAF at Parkfield puts in contact two very different geological structures, with granite on the west and the metamorphosed Franciscan formation on the east. However, at the end of a period of smoother variability of the attenuation parameters, the same bifurcation appears on both sides of the SAF, reinforcing our hypothesis of stress localization and volume compaction.

This study aims at the recognition of the pre-critical state of the transitional segment of the SAF at Parkfield, in the hypothesis that the fluctuations of the attenuation parameter represent a proxy for the stress conditions in the crustal volume crossed by an active fault. We argue that, along the fault segment under scrutiny, the signature of the pre-critical state of the 2004 mainshock is represented by the variations in opposite directions of the low- and the high-frequency attenuation parameters as a function of time (low-frequency attenuation increases sharply, high-frequency attenuation decreases). The bifurcations are observed on both sides of the fault, and last 6–8 weeks. The pre-critical state of the Parkfield asperity emerges after a period of smooth variability of the attenuation parameter that lasts for about 1 year.

The bifurcation of the stacked seismic attenuation is the result of the different trends followed by the populations of cracks in the two ranges of frequency at both ends of the available spectrum: 2–5 and 30–50 Hz. To interact with the seismic energy, cracks’ lengths must form permeability structures with dimensions comparable to the wavelengths carried by the seismic radiation: 400–1,500 m for the 2–5 Hz band, and 60–100 m for the 30–50 Hz band. Cracks in the permeability structures whose dimensions are comparable to the short-wavelength radiation tend to close because crack coalescence causes strain localization into the macroscopic fracture, and compaction and stress relaxation in the surrounding volume. The long-wavelength radiation, on the contrary, gets more attenuated because the longer cracks that dominate the coalescence process increase the bulk permeability in the volume immediately around the fault plane by enhancing their aperture and interconnection. Finally, the 2–5 Hz band is dominated by surface waves that sample a shallower crustal structure than the direct S-waves of the 30–50 Hz band.

Another important property of the attenuation parameter in 2004 is the sudden reduction of its variance σ 2 δ Q S − 1 t , f on both sides of the SAF, starting right after the occurrence of the San Simeon earthquake. A similar reduction of 2 [(1/Q)] can be observed starting at the beginning of 2021 on the sole Pacific side of the fault. If we currently are at the end of the seismic cycle, the fact that the attenuation signature is different from what observed in 2004 may be due to one, two, or all the three following reasons:

1. The sudden drop in variance in 2004, which occurred on both sides of the SAF, was due to a dramatic acceleration in the process of crack coalescence from the stress push from the San Simeon earthquake. Once the coalescence is mature, further changes in stress would not change the permeability structure, resulting in smoother time histories. Nothing like that occurred recently, and it is possible that the variance will go down gradually on both sides of the SAF in the next few months, following a gradual increase in the tectonic load.

About the intertwined behavior of δ Q S − 1 t , f and the NVT activity, it is possible that the increased permeability along the very fault plane under subcritical conditions allows pore pressure to diminish at NVT depth (in our case, this would correspond to the observed decrease in NVT activity that starts in 2020), and to an upward pore fluid migration along the fault plane. In turn, the latter may weaken the fault and facilitate the occurrence of the mainshock. Within the limits of this hypothesis, monitoring the fluctuations of NVT activity along the SAF would be useful in recognizing the pre-critical state of the fault.

The quiet period preceding the 2004 mainshock is consistent with our idea of upward fluid migration, and with the twofold role of the latter: (i) a lower pore pressure at NVT depth that inhibits the tremor activity; (ii) the increased pore-fluid pressure upward that lowers the fault strength and facilitates the nucleation of the mainshock. Likely, the pre-mainshock variability of the attenuation parameter is small due to the stabilizing action of strain localization upon δ Q S − 1 t , f . What causes a foretremor (i.e., a sharp increase in NVT with limited duration, occurring after a quiet period) is not clear, and the time lag between a foretremor and the mainshock may correspond to the diffusion time needed to reach a sufficient pore-fluid pressure along the fault. Consistently with our hypothesis, larger events (that primarily occur above the tremor zone) do not occur during the actual times of most tremor episodes. Rather, there is a definite tendency for the larger events to occur after the tremor episodes by a few days throughout the longer-term period (Figure 5 of Guilhem and Nadeau, 2012; Supplementary Figure S4 of the Supplementary Material of the same paper).

In the next few months (or years) the Parkfield segment of the SAF will provide most of the answers to our questions. In the meantime, it may be worthwhile applying our technique on different tectonic settings characterized by more complex fault geometries and ongoing NVT activity. Aside from more advanced theoretical developments and laboratory experiments on seismic attenuation at increasing differential stresses, the most important ingredient for a successful use of attenuation and NVT information would be the use of data from dense seismic networks equipped with borehole sensors, recording at high sampling rate (250–500 sps for borehole stations; max 200 sps for surface instruments, see Malagnini et al., 2021).