The 2008 M_S8.0 Wenchuan earthquake was one of the most destructive continental earthquakes in recent decades. It caused massive human casualties and severe economic losses. Nearly eighteen years after the mainshock, seismic activity within its aftershock zone remains elevated and has not yet returned to the background level. This raises two key scientific questions: Are recent earthquakes in the Wenchuan aftershock zone aftershocks or background events? If they are aftershocks, how long will the aftershock sequence persist? Herein, we address these questions via a statistics-based approach relying on long-term seismic catalogs instead of purely physical modeling.

We first compiled the 1970–2024 seismic catalog for the Sichuan–Yunnan region, which covers the Longmenshan fault zone, the seismogenic belt of the Wenchuan earthquake. Earthquakes with M≥4.0 form the primary analytical dataset owing to complete catalog coverage spanning the pre- and post-mainshock periods; events with M≥3.0 serve as an auxiliary dataset to test result robustness. We adopt the nearest-neighbor （NN） algorithm to separate aftershocks from background earthquakes by quantifying temporal, spatial and magnitude correlations that characterize earthquake clustering. Compared with conventional declustering algorithms （e.g., the Gardner–Knopoff time-space window method and the epidemic-type aftershock sequence （ETAS）model）, the NN method bears prominent strengths: it requires no preset empirical parameters, with classification thresholds objectively derived from observational data. Such merits render the NN framework particularly suitable for resolving long-lived aftershock sequences, especially within continental tectonic settings where aftershock activity can sustain for decades. Additionally, the method is less susceptible to catalog incompleteness, a critical advantage for practical observational seismic data.

Prior to NN implementation, two core parameters—the fractal dimension （d_f） and Gutenberg–Richter b-value—are quantified. The fractal dimension constrains the spatial clustering of earthquakes and mirrors the geometric complexity of fault systems, whereas the b-value describes the relative frequency of small versus large earthquakes according to the Gutenberg–Richter scaling law. We compute the fractal dimension via the correlation integral technique and estimate the b-value using maximum-likelihood inversion. With calibrated parameters and catalog data, we calculate the nearest-neighbor distance （NND） for each earthquake and examine its statistical distribution. The resultant NND distribution across the Sichuan–Yunnan region features a distinct bimodal pattern: the low-NND peak corresponds to clustered aftershocks with strong inter-event correlation, while the high-NND peak represents loosely correlated background seismicity. Using this bimodal signature, we first set a preliminary NND cutoff for aftershock discrimination. Focusing on the Wenchuan aftershock domain, we further optimize this threshold under the widely accepted constraint of stationary background seismicity before and after the mainshock, yielding an optimal cutoff and stable declustering outputs. Another methodological improvement of this work is the treatment of the Wenchuan mainshock as an extended finite fault source rather than a point source. Given its several-hundred-kilometer rupture along the Longmenshan fault, we discretize the mainshock into multiple sub-events distributed along the rupture trace. This modification mitigates distance-calculation biases for distal aftershocks far from the epicenter and improves the reliability of aftershock/background discrimination.

Following declustering, we compute the temporal evolution of background seismicity rate and quantify the temporal variation in the aftershock fraction of total seismicity. The aftershock proportion decreases monotonically with time: it peaked at about 97% immediately after the mainshock and has fallen to approximately 75% in recent years. Aftershocks still dominate contemporary seismicity across the study area, demonstrating persistent post-Wenchuan perturbation on regional tectonic stress. We adopt the Omori–Utsu law, which governs the temporal decay of aftershock occurrence rate, to constrain the total lifespan of the Wenchuan aftershock sequence. By fitting observed seismicity decay to the Omori–Utsu formulation and extrapolating toward the background baseline rate, we derive the timescale required for regional seismicity to recover to pre-mainshock background levels, i.e., the total duration of the aftershock sequence. We further incorporate parametric uncertainties via a 100 000-trial Monte Carlo simulation to enhance estimate reliability. Our results suggest that the Wenchuan aftershock sequence will continue for 45–160 years, with bounds varying with magnitude cutoff and confidence interval selection; the wide interval originates from uncertainties in raw catalog data and statistical model assumptions. Our estimates are generally consistent with prior physics-driven simulations, and minor discrepancies likely arise from divergent assumptions on tectonic loading rates and postseismic processes including viscoelastic relaxation, fault afterslip and pore-fluid migration. Sensitivity tests against study-region boundary, magnitude threshold and parametric settings show that predicted aftershock durations fluctuate within an approximate10-year window, verifying the robustness of our conclusions.

In summary, the NN-based declustering approach proves effective for aftershock identification and long-term aftershock-duration evaluation. The Wenchuan aftershock sequence will persist for several decades to over a century, with essential implications for regional seismic hazard assessment and disaster risk governance.