Aftershocks are defined as events of generally lesser magnitude following a larger main shock that are distributed around the fracture area of the main shock (e.g., Uidas, 1999). In 1894, Fusakichi Omori observed that the rate of aftershocks decays as the inverse of time following a main shock (Omori and Coll, 1894). An improved fit to observed aftershock rates can be made by applying Tokuji Utsu’s (1961) formulation as k/(c + t) ^p , where k, c, and p, are constants that can vary with different sequences, with the exponent p usually ranging from ∼0.7 to 1.5. Additional modifications to Utsu’s formula were made by Shcherbakov et al. (2004) to incorporate earthquake magnitude scaling laws, as did Hainzl and Marsan (2008). It was noted by Utsu et al. (1995) that the Epidemic Type Aftershock Sequence (ETAS) models by Yosihiko Ogata (e.g., Ogata, 1988) achieved an improved representation of aftershock characteristics because these models account for secondary triggering (cascades) of earthquakes by aftershocks of the main shock (e.g., Ogata, 1998; Ziv, 2003).

In this paper we explore fits of a blend of physics-based aftershock rate models to observations through regression methods. We assess the relative roles of acknowledged physical causes of earthquake triggering to better understand their role in Omori’s law under the possibility that more than one cause is operating in any given aftershock sequence. We consider three primary modes of earthquake triggering within aftershock sequences and their relationships to Omori’s law: (1) static stress and/or dynamic triggering caused by slip of the mainshock, (2) fluid diffusion triggering, and (3) secondary triggering through aftershock cascades by static and dynamic earthquake triggering. Each of these triggering modes have physics-based explanations for Omori’s law. We explicitly do not include models that identify parameters based on empirical observations, and the point of this study is to see if any one, or a blend of physical models can fit observations. We do not intend to solve aftershock occurrence globally. This paper is instead intended to show a method on how we might assess the relative roles of different physical models/concepts in aftershock sequences. We intentionally fit a single aftershock sequence that is associated with a very high-resolution catalog as a proof of concept.

We combine relations describing the three triggering modes into one equation and we minimize the number of free parameters by combining constants together to broadly represent the physics of aftershock triggering. We use that relation and its parameters to fit an observed aftershock sequence from a high-resolution earthquake catalog of the 2016 Central Apennines, Italy earthquake sequence that is complete to M1.04 (Tan et al., 2021). We focus on the 64-day series following the initial mainshock, the M6.0 Amatrice earthquake, because this well-studied, high-quality catalog occurred in a region where there is evidence that multiple physical processes influencing aftershock occurrences have occurred. We use two different nonlinear regression methods for fitting to assess method dependence.

Aftershocks occurring on or near the rupture plane have been explained as a result of incomplete mainshock rupture and/or heterogeneous slip (e.g., Bullen and Bolt, 1985). Elastic strain caused by the mainshock has also been invoked to explain aftershocks away from the primary rupture plane as early as Clarence Dutton’s (1904, p. 62) statement, ‘‘. . . portions of the raised or sunken mass are subject to stresses which from time to time give rise to further movements. . . .’‘. Dutton’s (1904) and Bullen and Bolt’s (1985) explanations are consistent with modern concepts of stress transfer caused by fault slip that mostly correlate with the spatial pattern of off-fault aftershocks (Figure 1) (e.g., Harris, 1998 and references therein; Stein, 1999; Freed, 2005). Mancini et al. (2019) showed that static stress-based forecasting during the 2016 Central Apennines sequence was as effective as statistical ETAS simulations, particularly when secondary triggering by M ≥ 3 earthquakes was included.

It is also well established that earthquakes can be triggered by changes in pore fluid pressure that can alter the frictional state and contact area on a locked fault by forcing it open (Biot, 1941; Hubbert and Rubey, 1959). Human induced pore fluid pressure changes are observed to trigger earthquakes (Evans, 1966; Ellsworth, 2013 and references contained therein), and observations have also been made of fluid pressure changes caused by an earthquake that in turn triggers other earthquakes (e.g., Nur and Booker, 1972; Malagnini et al., 2012; Tung and Masterlark, 2018; Albano et al., 2019; Miller, 2020; Kato, 2024). Fluid diffusion triggering is observed within the Central Apennines mainshock sequence (e.g., Albano et al., 2019; Chiarabba, et al., 2020; Malagnini et al., 2022) (Figures 2, 3). In Figure 2, earthquake clusters are identified that have distinct progressions through time (t) and space (d) as d t ∝ t in a similar mode as observed in the Apennines by Malagnini et al. (2012).

Additionally, elastic dislocation modeling by Malagnini et al. (2022) of dilatation induced by mainshocks within the 2016 Central Apennines sequence and modeled fluid flow changes imply that higher fluid pressures might be expected at sequential mainshock locations (Figure 3). While not proving a pore fluid pressure earthquake triggering mode, these models and observations are consistent with the process. That and the abundant observations of human induced earthquakes (e.g., Evans, 1966; Ellsworth, 2013 and references contained therein) support the inclusion of a fluid diffusion model in our exploration.

Earthquakes are also observed to be triggered dynamically by passing seismic waves (e.g., Hill et al., 1993). The predominant indicator of dynamic triggering is that its aftershock spatial decay is more gradual than that from a static stress change, which diminishes as a function of distance d, as d ^-3 . Dynamic triggering close to the mainshock is most likely caused by shorter period direct P- and S-wave phases with amplitudes that decay with distance more gradually than static stresses (e.g., Parsons and Velasco, 2009). Interpretations of observed spatial aftershock decay rates vary (e.g., Felzer and Brodsky, 2006; Richards-Dinger et al., 2010), making it difficult to determine how important dynamic triggering is in the near source region. Earthquakes can be triggered at global distances by longer period surface waves (e.g., Velasco et al., 2008; Parsons et al., 2014 and references contained therein). The proportion of earthquakes triggered dynamically was estimated to be 34% vs. 66% static triggering after 3 M ≥ 7.0 California mainshocks by Hardebeck and Harris (2022), based on the occurrence of aftershocks in static stress-reduced areas (also visible in Figure 1). Similarly, Parsons (2002) found a ratio of 31% dynamic vs. 69% static triggering from M ≥ 5 aftershocks of global M ≥ 7.0 mainshocks.

Post mainshock viscoelastic relaxation in the deep crust and upper mantle can increase stresses on upper crustal faults that trigger earthquakes, and the duration of these processes can be much longer than the early phases of aftershock production (e.g., Freed and Lin, 2001; Pollitz and Sacks, 2002; Morikami and Mitsui, 2020). Earthquakes also appear to be triggered by adjacent slow slip zones (e.g., Sammis et al., 2016; Segou and Parsons, 2018; 2020), this process can also take years to occur, though Moutote et al. (2023) observed a slow slip process that encompassed the 2017 Valparaiso M6.9 foreshock-mainshock-aftershock sequence. Aftershocks have been attributed to viscoelastic stress recovery on and adjacent to the mainshock fault (e.g., Yamashita, 1979; Zhang and Shcherbakov, 2016). Mikumo. (1979) assumed a distribution of viscoelastic relaxation times, which then numerically gave a power law decay in the rate of aftershocks.

Key questions addressed here are, how do the aftershock triggering modes discussed above translate into Omori’s law, and which are most influential? Dieterich (1994) showed that Rate/State friction laws (Ruina, 1983) can explain Omori’s law by the spatial distribution of times to instability of nucleation processes. The concepts are also consistent with changes in stressing rate after a mainshock caused by afterslip, and/or viscoelastic processes (e.g., Pranger et al., 2022). The aftershock rate R(t) can be expressed as R t = r exp − Δ τ a σ − 1 exp − t t a + 1 (1) where r is a reference seismicity rate, Δ τ is the static stress change induced by the mainshock, a is a dimensionless rate/state parameter, σ is normal stress, and t _a is the aftershock duration which, in this paper, is the time from the 24 August 2016 Amatrice mainshock up to the 26 October 2016 Visso mainshock (Figure 4). Dieterich (1994) describes the process as “Seismicity is modeled as a sequence of earthquake nucleation events in which the distribution of initial conditions over the population of nucleation sources and stressing history control the timing of earthquakes. The model is implemented using solutions for nucleation of unstable fault slip on faults with experimentally derived rate- and state dependent fault properties.” Rate/State friction can also be applied to fit dynamically triggered earthquakes to an Omori law decay if two assumptions are made: (1) the experimental observation of a power-law distribution of frictional contact areas (Dieterich and Kilgore, 1996) also occurs in natural fault zones, and (2) there is a constant change (damage) to the critical slip distance D _c of frictional contacts caused by passing seismic waves (Parsons, 2005). Additionally, Rate/State friction can be applied to the concept of time-varying viscoelastic stress recovery and relaxation times Yamashita, 1979, Mikumo. (1979).

Fluid diffusion triggering is observed to be time dependent, and Nur and Booker (1972) developed an expression based on dislocation models in a porous medium as d N d T = 1 Ω ∫ V ∂ P ∂ t d V (2)

Where P is pore fluid pressure, V is volume and Ω is a constant. Malagnini et al. (2012) derived an expression for 1-D fluid diffusion under a steady state fluid pressure change based on the Nur and Booker (1972) relation as d N d T = 1 α P 0 − P 1 2 π D t 3 ∫ 0 ∞ x e x p − x 2 4 D t d x = P 0 − P 1 D Ω π 1 t (3) which implies that d N d T ∝ 1 t (4)

Thus, the rate of fluid diffusion triggered aftershocks differs from the Omori-law time dependence of d N d T ∝ 1 t (5) meaning that fluid diffusion may need to be accompanied by other aftershock-generating mechanisms to achieve Omori time dependence.

Cascades of secondary aftershock triggering can have important effects on the Omori decay after a mainshock. They create multiple series of peaks and decays that overprint the primary decay curve and scale with the magnitude of each aftershock (e.g., Ogata, 1998; Helmstetter and Sornette, 2002; Ziv, 2003; Ouillon and Sornette, 2005; Hainzl and Marsan, 2008). We use a randomized scaling term to fit secondary aftershock triggering because, while the rate of triggered events by the mainshock is predictable, the magnitudes of these events are not. This is an issue because larger aftershocks can trigger enough secondary aftershocks to significantly perturb the aftershock rate curve at unpredictable magnitudes within a sequence.

Putting the three terms together (Equations 3–5) and combining constant terms results in a nonlinear mixed effects problem as d N d T = r r e f Δ − 1 exp − t f + 1 + D t + c R 0 , 1 t (6) where t is time. The first term is the static and/or dynamic triggering rate, with r _ref being the reference (background) steady state seismicity rate. We treat r _ref as an unknown constant because the high-resolution catalog contains only ∼10 days of pre-mainshock observations, which is not long enough to constrain the background rate. We simplify the static stress change term exp − Δ τ a σ from Equation 1 as an unknown constant Δ . The value f is the aftershock duration, which is constant. The second term is fluid diffusion triggering, and we apply an unconstrained constant D that represents diffusivity (Malagnini et al., 2012). The third term is the secondary triggering term with an unconstrained scaling constant c multiplied by a random number generator ranging between 0 and 1 (R ^0,1). Our philosophy is to allow unconstrained ranges of values for free parameters to minimize bias, and to see if there are consistent trends in the results that can teach us about the underlying physics of Omori’s law. Of particular interest to us is the relative contributions between static stress triggering with Rate/State decay, fluid diffusion, and secondary triggering, allowing for the possibility that one or more terms might have ∼ zero contributions to solutions that match the observed aftershock decay curve with minimum RMS (root mean square) misfits.

We use two regression approaches for fitting Equation 6 to an observed aftershock sequence, which in this case is the 64-day sequence after the 24 August 2016 M6.0 Amatrice earthquake up to and just prior to the 26 October 2016 M5.9 Visso and 30 October 2016 M6.5 Norcia mainshock earthquakes (Figure 4). We are fitting the high-resolution earthquake catalog of the 2016 Central Apennines, Italy earthquake sequence by Tan et al. (2021) that is complete to M_L>0.3 (M > 1.04). It is important to have as complete a catalog as possible to best represent Omori’s law because many smaller events can be undetected, especially during the earliest parts of the sequence when so many low-magnitude aftershocks are occurring that it can be difficult to resolve all of them.

Our hypothesis is that there is a mix of at least three physical causes of aftershocks that operate simultaneously. We thus have four unconstrained parameters to solve for in the equation for combined physical models of Omori’s law (Equation 6). The first method we use is a simulated annealing regression (Kirkpatrick et al., 1983; Černy, 1985) to find minimum RMS misfits to the observed sequence. Minimizing RMS error is a nonlinear and multivariable optimization procedure because the observed and theoretical decay curves (Equation 6) are nonlinear functions in time. Simulation of the observed curve is begun by using a stochastic sampling method for four free parameters. It starts from a random initial state and continues until a maximum of one million steps have been taken to minimize the RMS misfit between the model and the observed curve (Figure 4). All parameters are completely unconstrained in the simulated annealing regression and can be negative or positive, which can mean over or under predicting aftershock rates.

The underdetermined nature of the problem requires multiple iterations, so we calculate groups of 100 solutions minimized over one million attempts, which are then sorted for the lowest RMS misfit. Parameter values form the lowest RMS solution are then used as a starting point to explore neighboring points for better fits by perturbation within ±20% of their values. We allow those that increase the objective by returning higher RMS misfits to avoid trapping in local minima. Finally, the solutions are again refined by perturbing parameters within ±10% of their values and are limited to solution parameters that have lower RMS misfits than the input values. The compute cost is low enough that no cooling schedule is needed, which is sometimes used in simulated annealing to begin restricting the solution space to reduce the computational load. We run 1,000 independent regressions following the above-described steps to reasonably explore the solution space, and to identify trade-offs between parameters. We show the best fit solution, and the disaggregation of the three components from Equation 6 in Figure 5. This solution shows the contributions of physical models to be 114% static stress/Rate-State triggering, 9% fluid diffusion and −23% secondary triggering. We achieve low RMS misfits because the random secondary triggering term enables fitting to rate fluctuations associated with higher magnitude aftershocks; these fluctuations are not noise, but are an important part of the overall signal (e.g., Ogata, 1998; Helmstetter and Sornette, 2002; Ziv, 2003; Ouillon and Sornette, 2005; Hainzl and Marsan, 2008).

We employ a second regression to test if there is method dependence on the results. We use a local interior-point optimization method for each of 1,000 random realizations of the secondary triggering term to minimize misfit to Equation (6). Interior-point regression involves a gradient search method (Forsgren et al., 2002; Boyd and Vandenberghe, 2004) that achieves optimization by going through the middle of a solid (e.g., convex hull) defined by the constraints rather than around its surface. The constraints together with the objective function (minimizing misfit) are combined using a logarithmic barrier function with the necessary Karush-Kuhn-Tucker (KKT) conditions. The nonlinear system of equations is then solved using Newton’s method. Interior-point implementations rely heavily on very efficient code (e.g., Cholesky decomposition) for factoring sparse symmetric matrices. Convergence for the interior-point method is determined by an augmented Lagrangian merit function (Birgin and Martínez, 2014). The interior-point method differs slightly from our application of the simulated annealing regression because the Rate/State reference rate (r), parameter (Δ), and diffusivity constant (D) are subject to positivity constraints.

We found that regularization was required for the interior-point approach to get low misfit solutions. This is the process of adding an additional constraint to the model to reduce its complexity by forcing certain predictor variables to have a smaller impact on the outcome or, no effect at all. We use LASSO (Least Absolute Shrinkage and Selection Operator) (Tibshirani, 1996) to help with feature selection and reduce the complexity of the model. It does this by regularizing the coefficients of each predictor variable, meaning it applies an L1 penalty for large coefficient values to bring them down to a size that is more manageable. This helps overall model accuracy by reducing the number of variables used while improving predictive power. LASSO can set the coefficient of some predictors toward zero, effectively eliminating them from the model. This helps to reduce complexity further and improve interpretability by reducing the number of variables included in a model.

As an alternative to using independent random variables for the secondary triggering term as given in Equation 6, the interior point algorithm also applied 1,000 realizations of an Ornstein-Uhlenbeck (O-U) random walk. The O-U random walk has random variations from a central expectation (zero for the secondary triggering term) that grow with time, and a restoring force that pulls values back toward that expectation (reversion parameter) that controls the dominant period of random fluctuations. We experimented with different reversion values to minimize cumulative RMS misfits. We show the best fit solution in Figure 6, which shows the relative contributions of physical models to be ∼120% static stress/Rate-State triggering, and −20% secondary triggering, with fluid diffusion triggering being effectively zero. The best-fit solutions from the interior point methods come from applying independent random variables for the secondary triggering term, whereas applying the O-U random walk to solutions causes consistently larger RMS misfits (see Supplementary Figure S2).

We find multiple regression solutions with varying RMS misfits. We calculate 1,000 regressions with each method to understand the relative contributions to the overall model from each term representing a physical model. We find the contributions by summing up the number of predicted earthquakes made by each term and dividing that number by the total prediction of all three terms. In Figure 7 we show distributions of the relative contributions of each of the three terms in Equation 6. The range of contributions from both methods are remarkably consistent (Figures 7A,B) with distribution means and standard deviations that are nearly the same, with the Rate/State contribution mean of 1.170 for the simulated annealing regression vs. 1.169 for the interior point, with standard deviations of 0.064 and 0.063. For the secondary triggering term, the means are the same at −1.19 with standard deviations of 0.059 and 0.063. The fluid diffusion distributions are non-Gaussian. In general, the Rate/State model alone overpredicts the aftershock rate by ∼0–30%, while secondary aftershocks cause a negative contribution ranging between ∼0 and 30%, with a fluid diffusion contribution ranging from ∼0 to 20%. These results come from very different regression methods, implying that this may be a necessary outcome based on the how the mixed effects problem was posed in Equation 6, as well as our choice of the observed sequence. Additionally, Rate/State parameters can be distance dependent (e.g., Page et al., 2024), meaning that applying a fixed set of parameters for the entire sequence could cause misfit. Most aftershocks of the M6.0 Amatrice are contained within a 25 km radius of the mainshock (Figure 3A).

While the two regression methods show very similar distributions of relative physical-model contributions, we find that the distributions of RMS misfit differ (Figure 8) between methods. There is overlap but the interior point regressions tend to have higher misfit values than do the simulated annealing regressions.

We explore further to see if there are characteristic sets of parameters that cluster in a lower RMS state, which might lead us to further conclusions about the relative influence of different physical models for Omori’s law. We plot the relationships between the Rate/State static stress change term, secondary triggering term, RMS misfit, and number of solutions for both regression methods (Figure 9). These plots reveal that the most frequent solutions are not associated with the lowest RMS misfits.

We plot each parameter against RMS misfit for both regressions (Figures 10, 11) in isolation to determine if there are any weights or parameter values that have relatively stronger influences on fitting the observed aftershock sequence. We do not find what appear to be functional relationships between the free parameters and RMS misfits when they are isolated. We show linear regression fitting to the broadly scattered data for reference, but the data are not related through linear functions. We note that the simulated annealing regression results show parameters that are more scattered as function of RMS misfit, whereas the interior point values are more tightly clustered.

We find that all solutions from both regressions for the c parameter in the random secondary triggering term of Equation 6 yield negative values. It is the only parameter to have negative values amongst those that have no positivity constraints (positivity is limited to weighting the Rate/State reference rate (r), Rate/State parameter (Δ), and the diffusivity constant (D) in only the interior point method). Thus, because of the consistently negative c parameter in our solutions, the third term in Equation 6 is always subtracted from the first two. We provide a physical explanation of this in a later section of the paper.

We plot parameters as histograms (Figures 12, 13) that enable us to identify the most frequent values. If there is more than one peak on the histograms, then we have a suggestion of a tradeoff between parameters. We find that the same tradeoffs that were evident in Figure 7 are also seen in the histogram plots, between weighting of the Rate/State and secondary triggering components.

In this section we discuss two initially surprising results that are consistent across both regression methods. These are: (1) the uniformly negative c parameter values that scale the random secondary triggering factors of Equations 6, and (2) The relatively low weighting of fluid diffusion processes in fitting the observed aftershock sequence including near-zero weighting from the interior point algorithm, despite independent observations of these processes in the Apennines (Malagnini et al., 2012; Albano et al., 2019; Chiarabba et al., 2020; Malagnini et al., 2022). These results make sense to us upon reflection as described below.

We apply two regression methods towards fitting an aftershock sequence to an equation that represents three physical models/concepts for earthquake triggering to understand the relative influences of: (1) static/dynamic stress changes from the mainshock, (2) fluid diffusion triggering, and (3) the random magnitude distribution vs. time amongst secondary triggering by aftershocks. There are four unconstrained parameters in Equation 6, which combines the three physical model terms. These are a reference seismicity rate, a static stress change term, a fluid diffusivity term, and a scaling constant for random secondary triggering. These four constants are solved by minimizing misfit to the first 64 days of aftershocks following the 24 August 2016 M6.0 Amatrice mainshock (Figure 4). We calculate 1,000 solutions from each regression method and find some consistent patterns amongst the solutions.

In all solutions the Rate/State model alone overpredicts the aftershock rate by 100%–130%, which is balanced by the random secondary triggering model which underpredicts by 0% to −30%. The fluid diffusion model predicts a 0%–20% contribution. Both regressions methods we used returned very similar and consistent results (Figure 7), including negative scaling factors for secondary triggering that appears to be required to fit the observed sequence that were not imposed. We conclude this is a necessary subtraction because while each earthquake that follows the mainshock triggers subsequent aftershocks, the net effect of each of them is to remove stress from the crust, leading to a self-extinguishing process (Figure 15).

Fluid diffusion triggering is consistently weighted low in both regression methods despite evidence that it is an important process in the Apennines (e.g., Bonini, 2007; Malagnini et al., 2012; Albano et al., 2019; Chiarabba, et al., 2020; Malagnini et al., 2022) (Figures 2, 3). We conclude that this is because diffusion processes are limited in time and space and thus do not have a strong signal over a long aftershock sequence.

Unaddressed by our regressions are the physical causes of secondary triggered earthquakes that each have their own Omori Law sequences that contribute to the overall signal. We have modeled this as rate effects from a uniform aftershock magnitude distribution (Figure 15b) that is subtracted from the expected sequence of the primary mainshock. We note here that there may be more opportunities for dynamic triggering within the sequence than from a single mainshock, to which roughly a third of the aftershocks have been attributed if fluid diffusion causes are neglected (e.g., Parsons, 2002; Hardebeck and Harris, 2022). These relatively small rate fluctuations that might be construed as noise are influential in our solutions such that they have as much as a 35% contribution to the overall aftershock time series distribution. This likely explains why statistical (e.g., Ogata, 1998; Ziv, 2003) or static stress change aftershock forecasts (Mancini et al., 2022) that update and account for secondary triggering are more accurate, which is a consideration for operational aftershock forecasting.

Results from regressions have enabled us to quantify ranges of the relative influences of physical processes on an aftershock rate distribution provided that the three models expressed in Equation 6 are correct and comprehensive. The primary conclusion that we reach is that any one physical model cannot alone fit the observed sequence as well as the combination of all three that we investigated (Figure 17). The concept of inverting for parameter values can be adapted to additional and/or alternative physical models in the future.

Aftershocks are most numerous immediately after a mainshock occurs, and we have known since Fusakichi Omori’s work in 1894 that they decay roughly as a function of 1/time. There are multiple observations/models for why this happens, and we are curious whether a mixture of these concepts is responsible, and if so, what is the relative importance of each. We distill these ideas into an equation of three terms: (1) a model of time-dependent friction on faults that causes a delayed response to stress changes instigated by the mainshock, (2) increased fluid pressure caused by the mainshock that pushes open fault walls enabling aftershocks to occur, and (3) aftershocks that can trigger additional aftershocks (secondary triggering). We use two different computational methods to find the best combinations of models that fit an observed aftershock sequence. We find that time dependent frictional responses to stress changes is the dominant cause of aftershocks as modified by the effects of secondary triggering. We find the impact of fluids in the crust is smaller for a long aftershock sequence but can be more important over shorter time periods, and that secondary aftershocks reduce stress in the crust and suppress continued seismicity.