Objective Estimation of Spatially Variable Parameters of Epidemic Type Aftershock Sequence Model: Application to California

The ETAS model is widely employed to model the spatio-temporal distribution of earthquakes, generally using spatially invariant parameters. We propose an efficient method for the estimation of spatially varying parameters, using the Expectation-Maximization (EM) algorithm and spatial Voronoi tessellation ensembles. We use the Bayesian Information Criterion (BIC) to rank inverted models given their likelihood and complexity and select the best models to finally compute an ensemble model at any location. Using a synthetic catalog, we also check that the proposed method correctly inverts the known parameters. We apply the proposed method to earthquakes included in the ANSS catalog that occurred within the time period 1981-2015 in a spatial polygon around California. The results indicate a significant spatial variation of the ETAS parameters. We find that the efficiency of earthquakes to trigger future ones (quantified by the branching ratio) positively correlates with surface heat flow. In contrast, the rate of earthquakes triggered by far-field tectonic loading or background seismicity rate shows no such correlation, suggesting the relevance of triggering possibly through fluid-induced activation. Furthermore, the branching ratio and background seismicity rate are found to be uncorrelated with hypocentral depths, indicating that the seismic coupling remains invariant of hypocentral depths in the study region. Additionally, triggering seems to be mostly dominated by small earthquakes. Consequently, the static stress change studies should not only focus on the Coulomb stress changes caused by specific moderate to large earthquakes but also account for the secondary static stress changes caused by smaller earthquakes.

💡 Research Summary

The paper addresses a long‑standing limitation of the Epidemic‑Type Aftershock Sequence (ETAS) model: the assumption that its governing parameters are spatially invariant. Recognizing that geological, thermal, and fluid‑related conditions vary across a region, the authors develop a statistically rigorous framework that allows ETAS parameters to change continuously in space. The core of the methodology combines the Expectation‑Maximization (EM) algorithm with ensembles of spatial Voronoi tessellations. In the EM step, each earthquake is probabilistically assigned to either the background seismicity component or to an aftershock cascade generated by preceding events. This yields posterior probabilities that serve as weights for updating the ETAS parameters. To capture spatial heterogeneity, the study constructs multiple ETAS models on Voronoi partitions of differing resolutions (i.e., varying numbers of cells). For each partition, the EM algorithm estimates a full set of parameters (background rate μ, productivity K, magnitude scaling α, temporal decay p, spatial decay d, and branching ratio n). Model selection is performed using the Bayesian Information Criterion (BIC), which penalizes model complexity while rewarding goodness‑of‑fit. The BIC‑optimal partitions are retained, and their parameter estimates are combined through a location‑specific weighted average, producing an ensemble map of each ETAS parameter across the study area.

The authors first validate the approach on a synthetic catalog where the true spatial variation of all parameters is known. The inversion accurately recovers the prescribed fields, and the BIC correctly identifies the optimal Voronoi resolution, demonstrating that the method can distinguish genuine spatial signal from over‑fitting.

The validated technique is then applied to real seismicity recorded in the ANSS catalog from 1981 to 2015 within a polygon encompassing California. The resulting parameter maps reveal pronounced spatial variability. The background rate μ shows relatively modest variation and no systematic correlation with surface heat flow, suggesting that far‑field tectonic loading contributes uniformly across the region. In contrast, the branching ratio n (the expected number of aftershocks per event) exhibits a strong positive correlation with heat‑flow measurements. Areas of high geothermal flux—often associated with volcanic arcs and active fluid pathways—display elevated n values, implying that fluid‑induced weakening enhances the efficiency of earthquake triggering. The productivity K and magnitude‑scaling α parameters also vary, but the dominant pattern is that small earthquakes (magnitudes below ~3) generate the majority of aftershocks; the estimated contribution of events below the magnitude of completeness exceeds 70 % of total triggering. This finding challenges the conventional focus on moderate‑to‑large events in static stress‑change studies and argues for incorporating the cumulative effect of numerous small ruptures.

Depth analyses reveal that neither n nor μ correlates with hypocentral depth, indicating that seismic coupling and triggering efficiency remain essentially constant from shallow to deeper portions of the crust within the surveyed area.

Overall, the study demonstrates that (1) spatially variable ETAS parameters can be robustly estimated using an EM‑Voronoi‑BIC framework; (2) the branching ratio is a sensitive proxy for fluid‑related processes, as evidenced by its heat‑flow correlation; and (3) small earthquakes dominate the cascade dynamics, implying that hazard models and stress‑change calculations should account for their collective impact. The authors conclude that incorporating spatial heterogeneity into ETAS modeling markedly improves the realism of seismicity simulations and offers a pathway to integrate geophysical observables (e.g., heat flow, fluid pressure) into probabilistic earthquake forecasting.