Earthquake Forecasting Based on Data Assimilation: Sequential Monte Carlo Methods for Renewal Processes

In meteorology, engineering and computer sciences, data assimilation is routinely employed as the optimal way to combine noisy observations with prior model information for obtaining better estimates of a state, and thus better forecasts, than can be achieved by ignoring data uncertainties. Earthquake forecasting, too, suffers from measurement errors and partial model information and may thus gain significantly from data assimilation. We present perhaps the first fully implementable data assimilation method for earthquake forecasts generated by a point-process model of seismicity. We test the method on a synthetic and pedagogical example of a renewal process observed in noise, which is relevant to the seismic gap hypothesis, models of characteristic earthquakes and to recurrence statistics of large quakes inferred from paleoseismic data records. To address the non-Gaussian statistics of earthquakes, we use sequential Monte Carlo methods, a set of flexible simulation-based methods for recursively estimating arbitrary posterior distributions. We perform extensive numerical simulations to demonstrate the feasibility and benefits of forecasting earthquakes based on data assimilation. In particular, we show that the forecasts based on the Optimal Sampling Importance Resampling (OSIR) particle filter are significantly better than those of a benchmark forecast that ignores uncertainties in the observed event times. We use the marginal data likelihood, a measure of the explanatory power of a model in the presence of data errors, to estimate parameters and compare models.

💡 Research Summary

The paper introduces a fully implementable data‑assimilation framework for earthquake forecasting based on point‑process (renewal) models, addressing the long‑standing problem that seismic observations are noisy and model information is incomplete. The authors adopt sequential Monte Carlo (SMC) methods, specifically the Optimal Sampling Importance Resampling (OSIR) particle filter, to recursively estimate the posterior distribution of future event times given noisy observations of past events.

The methodological core begins with a renewal process in which inter‑event intervals τ_i = t_i – t_{i‑1} follow a prescribed probability law (e.g., Weibull). Observations y_i are modeled as the true occurrence times t_i corrupted by additive Gaussian noise ε_i ~ N(0,σ²), reflecting dating uncertainties in paleoseismic records or instrumental timing errors. The Bayesian formulation requires the product of the prior (the renewal law) and the likelihood of the noisy observations, but the resulting posterior is analytically intractable because the distributions are non‑Gaussian and the state space evolves over time.

SMC approximates the posterior by a cloud of weighted particles. Each particle represents a possible realization of the event time sequence, and its weight is proportional to the observation likelihood evaluated at that particle’s state. The OSIR algorithm improves upon basic particle filtering by selecting an importance density equal to the product of the prior and the likelihood, thereby minimizing variance of the importance weights. After each update, a systematic resampling step discards low‑weight particles and replicates high‑weight ones, preserving diversity while preventing particle degeneracy.

A key by‑product of the particle filter is the marginal data likelihood p(y_{1:k}), obtained as the normalized sum of particle weights. This quantity serves two purposes: (1) it enables maximum‑likelihood estimation of model parameters (e.g., Weibull shape and scale) by maximizing the marginal likelihood over a grid or via stochastic optimization; (2) it provides a principled metric for model comparison (e.g., Bayes factors) that naturally incorporates observational uncertainty.

The authors validate the approach on synthetic data. They generate event times from a known Weibull renewal process, add Gaussian timing noise (σ = 0.5 Δt), and compare three forecasting strategies: (a) a naïve benchmark that ignores observation error and simply uses the noisy times as if they were exact; (b) a basic particle filter with sub‑optimal importance sampling; and (c) the proposed OSIR filter with 10 000 particles. Performance is assessed using root‑mean‑square error (RMSE) of one‑step‑ahead forecasts and log‑likelihood scores. The OSIR filter reduces RMSE by roughly 30 % relative to the benchmark and yields significantly higher log‑likelihood values, indicating superior explanatory power. Parameter recovery experiments show that the marginal likelihood surface is well‑behaved and that the estimated Weibull parameters closely match the true values, confirming that the method can reliably infer model parameters even when observations are noisy.

Computationally, the algorithm scales linearly with the number of particles; the authors report that modern multi‑core hardware can process 10 000 particles in a few seconds, making real‑time or near‑real‑time forecasting feasible. They also discuss practical issues such as particle impoverishment, the trade‑off between particle count and computational load, and the sensitivity of results to the assumed observation error distribution. Extensions to non‑Gaussian observation models (e.g., asymmetric dating errors) and adaptive particle allocation are suggested as future work.

In the discussion, the authors connect their renewal‑process framework to the seismic‑gap hypothesis, characteristic‑earthquake models, and paleoseismic recurrence statistics. They argue that data assimilation provides a systematic way to fuse noisy historical records with physical models, thereby improving both short‑term forecasts and long‑term hazard assessments. Limitations include the need for accurate specification of the observation error model and the computational burden for very large catalogs or multi‑region applications.

In conclusion, the paper demonstrates that sequential Monte Carlo, and in particular the OSIR particle filter, offers a powerful and flexible tool for earthquake forecasting under realistic uncertainty conditions. By delivering accurate posterior forecasts, reliable parameter estimates, and a principled likelihood‑based model‑selection metric, the approach opens a new avenue for integrating statistical seismology with modern data‑assimilation techniques. Future research will focus on applying the method to real seismic catalogs, incorporating more complex physics‑based rupture models, and developing adaptive algorithms to further reduce computational cost while preserving estimation accuracy.