Investigation of Pre-Seismic Ionospheric TEC and Acoustic-Gravity Wave Coupling Phenomena Using BDS GEO measurements: A Case Study of the 2023 Jishishan Ms6.2 Earthquake

1. Introduction

The investigation of seismo-ionospheric anomalies has garnered significant attention as a critical precursor for short-term earthquake forecasting [1,2,3]. Since the first recorded ionospheric perturbations preceding the 1964 Alaska MW9.2 earthquake, advancements in Global Navigation Satellite System (GNSS)-derived Total Electron Content (TEC) monitoring have revolutionized the detection and characterization of ionospheric anomalies associated with seismic activity [4,5,6]. Empirical studies suggest that majority of earthquakes with magnitudes exceeding Mw6.0 exhibit complex pre-seismic TEC variations [7]. Then GPS-derived TEC monitoring has become the cornerstone for pre-earthquake ionospheric anomaly detection, with well-documented cases including: the anomalous TEC attenuation preceding the 1999 Taiwan Chi-Chi Mw7.7 earthquake [8], the multi-phase TEC fluctuations observed before the 2008 Wenchuan Mw7.9 earthquake [9], the pre-seismic ionospheric anomalies preceding the 2013 Mw7.7 Pakistan earthquake [10], the pronounced positive TEC anomalies over the epicentral and magnetic conjugate regions within two weeks preceding the 2021 Yangbi–Maduo twin earthquakes [11], and the g ionospheric TEC anomalies associated with the 2024 Mw7.5 Noto Peninsula earthquake [12]. Statistical analyses further corroborate the robust correlation between seismic events and ionospheric disturbances: Liu et al. identified that approximately 80% of TEC anomalies within five days prior to Mw≥6.0 earthquakes in Taiwan were seismogenic [13], while Le et al. demonstrated a statistically significant increase in TEC anomaly frequency preceding global Mw≥6.0 earthquakes compared to background levels [14].

However, the IPPs between space satellites and ground-based receivers exhibit spatially uneven variations due to satellite motion. Consequently, the derived TEC data do not directly represent temporal variations at a fixed ionospheric location and require mathematical interpolation or modeling for spatial homogenization, introducing inherent uncertainties [15,16]. In contrast, the GEO satellites of the BeiDou Navigation Satellite System (BDS) provide a unique advantage for long-term ionospheric TEC monitoring due to their quasi-invariant orbital configurations [17]. While geostationary orbits are perturbed by Earth’s oblateness and gravitational interactions with the Sun and Moon, onboard station-keeping mechanisms minimize positional deviations, ensuring that IPP fluctuations remain negligible relative to the ionosphere’s large spatial scale (thousands of kilometers). Thus, BDS GEO satellites can be approximated as stationary relative to Earth, enabling ground stations to continuously acquire TEC data from fixed IPPs. This capability effectively establishes a network of time-continuous, position-locked monitoring nodes, significantly enhancing the spatiotemporal resolution of ionospheric observations [18,19,20].

While statistical correlations between ionospheric anomalies and seismic events have been established, significant gaps persist in our understanding of the underlying physical mechanisms governing the genesis and evolutionary processes of seismogenic ionospheric disturbances. The lithosphere-atmosphere-ionosphere coupling (LAIC) framework has emerged as a pivotal theoretical model to explain seismo-ionospheric interactions, proposing three primary coupling pathways: chemical channeling, acoustic-gravity wave (AGW) propagation, and electrostatic field modulation [21,22,23,24,25,26]. Among these, AGW-mediated coupling has gained prominence due to its traceable vertical propagation from the lithosphere to the ionosphere [27]. Seismically excited AGWs interact with ionospheric plasma through vertical energy transfer, inducing detectable TEC oscillations. For instance, Korepanov et al. employed wavelet analysis to reveal atmospheric filtering effects on AGWs, demonstrating that high-velocity, long-wavelength acoustic-gravity wave signals can propagate through the atmosphere to become captured by the ionosphere, whereas low-velocity or short-wavelength signals tend to dissipate within the atmospheric medium [28]. Notably, such high-velocity, long-wavelength acoustic-gravity wave signals are predominantly excited by surface tectonic activities such as earthquakes, rather than originating from convective weather systems. Yang et al. further validated this pathway by observing spatiotemporal coherence between stratospheric AGWs and ionospheric TEC disturbances preceding the 2016 Kumamoto earthquake [29]. Phanikumar et al. observed synchronous AGWs in VLF electromagnetic waves and mesospheric ozone anomalies prior to the 25 April 2015 Nepal Gorkha Mw7.8 earthquake [30]. Ray et al. identified dual co-seismic TIDs during the 2024 Hualien M7.4 earthquake, with velocity bifurcation (0.78 km/s vs. 0.48 km/s) potentially linked to multi-source fault rupture dynamics, highlighting the complexity of AGW-driven coupling [31]. These findings indicate that AGW coupling with the lower ionosphere not only triggers anomalous plasma oscillations but also modulates neutral atmospheric composition perturbations. The observed phenomena distinctly illustrate the vertical propagation pathway of AGWs across the troposphere-stratosphere-mesosphere (lower ionosphere) system.

Regarding the coupling mechanism between ionospheric TEC anomalies and AGWs, the E-region dynamo current hypothesis remains the most widely accepted explanation. Based on the synchronized anomalies in VLF signals and TEC variations, it is inferred that AGWs propagating to the lower ionosphere (primarily the E-region) induce wind field perturbations [32,33,34,35]. However, the precise mechanisms governing TEC-AGW synchronization remain ambiguous. Besides, current seismo-ionospheric TEC investigations predominantly employ temporally coarse datasets (e.g., global ionospheric maps with 2-hour resolution), focusing on statistical anomaly detection over pre-seismic timescales spanning days to weeks. The lack of higher-resolution spectral decomposition of ionospheric perturbations has fundamentally hindered the definitive identification of causative frequency bands underlying these anomalies. Critically, conventional GPS-derived TEC measurements inherit systematic uncertainties from their medium earth orbit (MEO) satellites. Since continuous migration of IPPs at ∼25 km/min orbital velocity introduces spatiotemporal aliasing, the derived TEC time series intrinsically conflate temporal variations and spatial plasma gradients due to IPP trajectory evolution. Besides, the orbital motion induces the amplitude attenuation in detected AGW signatures compared to stationary observations.

In light of this, this study investigates the pre-seismic ionospheric TEC anomalies and AGW coupling phenomena associated with the 18 December 2023 Gansu Jishishan Ms6.2 earthquake, utilizing high-rate TEC data from the BDS GEO satellites within the continuously operating reference stations (CORSs) of Gansu Province in China. Then we analyzed the propagation process and spatiotemporal characteristics of ionospheric disturbances on the day of the Jishishan Earthquake from multiple perspectives, including the spatiotemporal distribution, propagation velocity, time-frequency domain, disturbance amplitude, and waveform of ionospheric perturbations. It further examined the relationship between TIDs during the earthquake preparation period and weak geomagnetic storms, aiming to provide theoretical insights for studying the coupling mechanism between the ionosphere and AGWs prior to earthquakes.

2. Materials and Methods

This chapter establishes a systematic framework for ionospheric TEC retrieval, modeling, and disturbance characterization using GNSS observations. Beginning with dual-frequency code and phase measurements, it derives slant TEC (STEC) while addressing differential code biases (DCBs) through pseudorange smoothing. A Single-Layer Model (SLM) is employed to map STEC to vertical TEC (VTEC), leveraging the quasi-static geometry of GEO satellites to isolate temporal ionospheric variations. Advanced signal processing techniques—including adaptive median filtering, Butterworth bandpass filtering, and Morlet wavelet analysis—are integrated to suppress noise, extract disturbance components (DTEC), and identify dominant spectral features.

2.1. GNSS Observations and the Geometry-Free Combination

Under idealized conditions neglecting multipath effects and measurement noise, the fundamental GNSS code and carrier-phase observation equations (in length units) can be expressed as:

$$\left\{\begin{array}{c}{P}_{r,i}^s={\rho }_r^s+{\alpha }_i·I_{r,1}^s+c·B_{r,i}+c·B_i^s\\ L_{r,i}^s={\rho }_r^s-{\alpha }_i·I_{r,1}^s+{\lambda }_i·N_{r,i}^s+{\lambda }_i·b_i^s\end{array}\right.$$

(1)

where $P_{r,i}^s$ and $L_{r,i}^s$ denote the code and carrier-phase measurements from receiver $r$ to satellite $s$ at frequency $i$, respectively; ${\rho }_r^s={\varrho }_r^s+c·d{t}_r-c·d{t}^s+T_r^s$ is the frequency-independent composite term, comprising the geometric range (${\varrho }_r^s$), receiver clock-bias ($d{t}_r$), satellite clock-bias ($d{t}^s$) and tropospheric delay ($T_r^s)$; $I_{r,1}^s$ is the first-order ionospheric delay at reference frequency $f_1$ (as adopted in this study) along the line-of-sight; ${\alpha }_i={\left(f_1/f_i\right)}^2$ denotes the ionospheric scaling factor; $B_{r,i}$ and $B_i^s$ represent the receiver- and satellite-specific code hardware delays, respectively; $N_{r,i}^s$ is integer ambiguity term for carrier-phase; $b_i^s$ is phase hardware delay (typically absorbed into ambiguity parameters); $c$ is the speed of light in vacuum; and ${\lambda }_i$ is wavelength.

The cumulative nature of GNSS carrier-phase observations enables millimeter-level measurement precision (2–3 orders of magnitude higher than code), making them ideal for code smoothing to enhance usable accuracy [36,37]. This study employs dual-frequency carrier-phase smoothed code (CSC) to decouple first-order ionospheric delays and hardware bias parameters, constructing a geometry-free (GF) combination with the following observation equations:

$$\left\{\begin{array}{c}{P}_{r,ij}^s=P_{r,i}^s-P_{r,j}^s=\left({\alpha }_i-{\alpha }_j\right)·I_{r,1}^s+c·\left(DC{B}_{r,ij}+DC{B}_{ij}^s\right)\\ L_{r,ij}^s=L_{r,i}^s-L_{r,j}^s=-\left({\alpha }_i-{\alpha }_j\right)·I_{r,1}^s+\left({\lambda }_i·N_{r,i}^s-{\lambda }_j·N_{r,j}^s\right)\end{array}\right.$$

(2)

where $P_{r,ij}^s$ and $L_{r,ij}^s$ epresent the inter-frequency differenced code and phase observables, respectively. $DC{B}_{r,ij}=B_{r,i}-B_{r,j}$and $DC{B}_{ij}^s=B_i^s-B_j^s$ denote the receiver- and satellite-end DCBs.

To achieve cascaded noise attenuation, a recursive weighting algorithm is implemented, expressed as:

$${\overline{P}}_{r,ij}^s\left(t\right)={\displaystyle \frac{1}{M}}{P}_{r,ij}^s\left(t\right)+{\displaystyle \frac{M-1}{M}}\left[{\overline{P}}_{r,ij}^s\left(t-1\right)+L_{r,ij}^s\left(t\right)-L_{r,ij}^s\left(t-1\right)\right]$$

(3)

where $M$ represents the adaptive smoothing window width (typical range: 2~100), initialized with ${\overline{P}}_{r,ij}^s\left(0\right)=P_{r,ij}^s\left(0\right)$ [38,39]. The algorithm implicitly assumes quasi-linear ionospheric delay variations within the smoothing window. Under severe ionospheric disturbances (e.g., scintillation), window parameters require dynamic adjustment to balance noise suppression and temporal responsiveness.

2.2. Ionospheric TEC Extraction

In electromagnetic wave propagation theory, the geometric optics approximation applies when the radio wavelength is much smaller than the characteristic spatial scale of the ionosphere. Under this assumption, the signal propagation path from satellite to receiver can be modeled by integrating the phase velocity along the trajectory [40]. The ionospheric delay effect caused by wave propagation through the ionized medium is quantified as:

$$I_{r,1}^s={\displaystyle \frac{K}{f_1^2}}\int N_e·ds={\displaystyle \frac{K}{f_1^2}}STEC$$

(4)

where $N_e$ is the electron density (unit: electrons/m³), $K={\displaystyle \frac{e^2}{8{\pi }^2{\epsilon }_0{m}_e}}\approx 40.309m^3/s^2$ is the ionospheric constant, and $STEC$ denotes the slant total electron content along the propagation path ($1TECU={10}^{16}electrons/m^2$).

Combining Equations (1) and (4), the CSC-based STEC calculation model can be derived as

$${\overline{P}}_{r,ij}^s=\left({\alpha }_i-{\alpha }_j\right)·{\displaystyle \frac{40.309·{10}^{16}}{f_1^2}}STEC+c·\left(DC{B}_{r,ij}+DC{B}_{ij}^s\right)$$

(5)

This equation highlights the amplification effect of dual-frequency differential observations on ionospheric delays. The scaling factor, expressed as $\left({\alpha }_i-{\alpha }_j\right)·{\displaystyle \frac{40.309·{10}^{16}}{f_1^2}}$, is determined by the frequency combination properties.

2.3. Ionospheric Single Layer Modeling

The vertical distribution of ionospheric electron density exhibits a distinct stratified structure (D/E/F layers), with the F2 layer contributing approximately 80% of the TEC [41]. To establish a mapping relationship between STEC and VTEC, this study employs the single-layer model (SLM) [42,43], a widely used theoretical framework based on the following assumptions:

(1) Ionospheric free electrons are confined to an infinitely thin spherical shell at a height H above Earth’s surface.

(2) Refraction effects are neglected, simplifying the signal propagation path to a straight line within the thin layer.

(3) Horizontal ionospheric gradients and spatial inhomogeneities are ignored.

Let Earth be modeled as a sphere with radius $R$, and the thin ionospheric layer height H is set to 506.7km (based on the F2 layer’s average peak height derived from the International Reference Ionosphere (IRI) model) [44,45]. The intersection of the satellite-receiver line-of-sight with the SLM defines the IPP, whose geometry satisfies:

$$\sin \left(Z\right)={\displaystyle \frac{R}{R+H}}\sin \left(z\right)$$

(6)

where $z$ is the satellite zenith angle at the receiver’s location, and $Z$ is the equivalent zenith angle at the IPP.

Currently, five BDS GEO satellites are positioned along the equator at longitudes of 58.75°E, 80°E, 110.5°E, 140°E, and 160°E [46]. Unlike MEO satellites, each operational GEO satellite orbits at an altitude of approximately 35,787 km, maintaining a fixed position relative to the Earth. The quasi-static geometric configuration between BDS GEO satellites and ground receivers results in minimal variation in the zenith angle, leading to a nearly constant value. This unique characteristic ensures that the TEC time series derived from GEO satellites predominantly reflect temporal ionospheric variations, while spatial variability is suppressed due to the satellites’ stationary nature.

2.4. TEC Time-Series Analysis and Feature Extraction

An adaptive sliding-window nonlinear smoothing method is employed to suppress transient spike noise [47]. Mathematically, this process is defined as:

$$\begin{array}{cc}{y}_{med}\left(t\right)=median\left\{x\left(t-k\right),\cdots ,x\left(t\right),\cdots x\left(t+k\right)\right\}& \left(k\in {\mathbb{Z}}^{+}\right)\end{array}$$

(7)

where the sliding window width is an odd integer ($2k+1$). This nonlinear sorting-based algorithm eliminates >95% of impulse interference while preserving edge features of the signal.

From Equation (5), the STEC derived from satellite observations includes contributions from differential code biases (DCBs), higher-order errors and random noise. The latter two are typically negligible, but DCBs (magnitude: several TECU, comparable to ionospheric disturbances) introduce non-negligible variations. Previous studies indicate that DCBs exhibit minimal diurnal variation and can be approximated as constants [48,49], enabling their removal via time-series filtering or background subtraction.

For signal processing, the Butterworth bandpass filter is preferred due to its maximally flat passband magnitude response and monotonic stopband attenuation [50,51]. A 4th-order Butterworth filter is adopted to extract ionospheric disturbance components (DTEC), with the transfer function:

$$H\left(f\right)={\displaystyle \frac{1}{\sqrt{1+{\left(f/f_c\right)}^{2n}}}}$$

(8)

where $f_c=\left[f_{low},f_{high}\right]$ defines the cutoff frequencies and $n$ is the filter order.

To optimize $f_c$ selection for TEC disturbance extraction, a Morlet wavelet transform [52] is introduced for time-frequency joint analysis:

$$\psi \left(t\right)={\pi }^{-1/4}{e}^{-i{\omega }_{0}t}{e}^{-t^2/2}$$

(9)

$$\widehat{\psi }\left(s{\omega }_0\right)={\pi }^{-1/4}\sqrt{s}{e}^{-{\left(s\omega -{\omega }_0\right)}^2/2}$$

(10)

where the scale parameter $s$ relates to the Fourier period $T$ via:

$$T=4\pi s/\left({\omega }_0+\sqrt{2+{\omega }_0^2}\right)$$

(11)

The global wavelet power spectrum is computed by time-averaging the squared wavelet coefficients:

$$P\left(s\right)={\displaystyle \frac{1}{N}}{\sum }_{t=1}^N{\left|W\left(s,t\right)\right|}^2$$

(12)

This spectrum approximates a smoothed Fourier power spectrum constrained by the wavelet window function. By identifying regions of significant spectral power (confidence >95%), the dominant frequency components of ionospheric disturbances are determined, enabling optimized Butterworth cutoff frequency selection.

3. Results

3.1. Study Area and Data Sources

An Ms6.2 earthquake (hereinafter referred to as the Jishishan earthquake) struck Jishishan County (35.70°N, 102.79°E) in Linxia Hui Autonomous Prefecture, Gansu Province, China, at 15:59 UTC on 18 December 2023. Moment tensor inversion reveals a hybrid rupture mechanism dominated by thrust faulting with a minor strike-slip component, characterized by a rupture duration of ~8 s and a centroid depth of 12.0 km [53]. The epicenter is located within the Lajishan Fault Zone along the northeastern margin of the Tibetan Plateau, a region experiencing intense compressional deformation due to the ongoing India-Eurasia continental collision-driven tectonic stress field since the Cenozoic era.

To delineate the precursor deformation zone, we calculated the seismogenic radius ($R$) using the Dobrovolsky empirical formula [54]:

$$R={10}^{0.43M}$$

(13)

where $M$ denotes the surface wave magnitude. This yielded a seismogenic radius of 380.2 km centered on the epicenter. We extended the analysis to three standard deviations (3σ) of the radius (1,140.6 km) and compiled seismic records from the U.S. Geological Survey (USGS) earthquake catalog (https://earthquake.usgs.gov/earthquakes). Over the preceding five years, this region experienced 18 seismic events with moment magnitudes (Mw) ≥5.5 (Table 1). Notably, Mw, calculated from the energy released during rock rupture, represents the USGS’s standard magnitude measurement.

The catalog comprises 18 seismic events recorded over a six-year period (2019–2025), yielding an annual average of 3.0 events/year. Temporal variability is evident, with peak activity in 2022 (7 events) and secondary clustering in 2021 (5 events). Notably, a single day (21 May 2021) recorded four distinct events, including a Mw7.3 mainshock and three aftershocks, highlighting spatially concentrated seismicity along the Lajishan Fault Zone. Magnitudes range from Mw5.5 to Mw7.3, with 83% of events ≤ Mw6.0. The Mw7.3 event dominates the cumulative energy-release budget, accounting for ~90% of the total seismic moment.

Spatially, events cluster between 32°–38°N and 92°–103°E, correlating with the Lajishan Fault Zone and adjacent structures along the northeastern Tibetan Plateau margin. The 2023 Jishishan earthquake (Ms6.2) at 35.74°N, 102.81°E exemplifies thrust-dominated tectonics characteristic of this compressional regime. Depth analysis reveals 89% of events occur at shallow crustal levels (≤15 km), consistent with intracontinental deformation patterns, while the deepest event (15 km on 10 November 2022) suggests mid-crustal fault locking.

Statistical characteristics reveal that the current tectonic activity in this region is in a low-energy state, consistent with the strain accumulation-release cycles observed in the northeastern Tibetan Plateau since the Late Quaternary [55]. Existing studies on the Jishishan earthquake predominantly focus on source parameter inversion, Radar-derived coseismic deformation modeling, and regional geodynamic mechanisms [56,57,58]. However, significant gaps persist in understanding the coupling mechanisms between seismogenic processes and ionospheric precursors, particularly in systematically characterizing the LAIC model and resolving multi-layer energy transfer pathways.

To address this, we conduct pre-seismic ionospheric disturbance analysis using multi-frequency observations from BDS GEO satellties. While the Jishishan seismic zone (25°–46°N, 92°–113°E) hosts over 100 continuously operating GNSS reference stations, only a limited subset can track dual-frequency BDS signals. By selecting 12 Gansu Province Continuously Operating Reference Stations (CORS) (epicentral distance <100.0 km) and 4 International GNSS Service (IGS) stations, we establish a 16-station BDS GEO monitoring network (Figure 2).

In Figure 2, a red pentagram marks the epicenter of the Ms6.2 Jishishan earthquake. Red circles denote GNSS stations with dual-frequency BDS observation capabilities. Blue concentric circles delineate multi-scale spatial analysis units with 200 km radial gradients centered on the epicenter. This network can capture trans-lithospheric signals associated with stress perturbations propagating from the seismogenic zone to the ionosphere.

3.2. Analysis of Space Weather Conditions on the Day of the Seismic Event

The ionosphere, a critical component of the Sun-Earth space environment, exhibits dynamic responses governed by the interplay of solar activity, geomagnetic disturbances, and lower atmospheric processes. To investigate potential links between pre-seismic ionospheric anomalies and the Jishishan earthquake (December 18, 2023), it is imperative to systematically exclude exogenous space weather influences. Geomagnetic disturbances are screened using the Kp index to quantify global geomagnetic activity and the Dst index to evaluate ring current intensity [59,60]. Solar activity monitoring focuses on tracking sunspot numbers (SN) to determine the solar cycle phase [61], F10.7 cm solar radio flux measurements as a proxy for solar extreme ultraviolet (EUV) irradiance variability [62], and analyzing solar wind parameters, including the orientation of the interplanetary magnetic field (IMF Bz) and variations in the cross-tail electric field (Ey) [63,64]. Simultaneously, plasma environment characterization integrates vertical profiles of plasma temperature (Pt) and density (Pd) with assessments of magnetosphere-ionosphere coupling efficiency [65].

As illustrated in Figure 3 (left panels), a moderate geomagnetic storm was identified between 06:00–08:00 UT, marked by a peak Kp index of 5.3 (typical threshold: Kp ≥ 4.5) and a minimum Dst index of -69 nT (typical threshold: Dst ≤ -60 nT). Post-disturbance parameters rapidly stabilized (Kp=3.7, Dst= -43 nT), indicative of transient energy injection rather than sustained magnetospheric driving. Solar activity remained moderately elevated, with an F10.7 radio flux of 156 sfu—marginally exceeding the Solar Cycle 25 baseline (146 sfu) but below the strong activity threshold (180 sfu). Sunspot numbers (SN=138) aligned with predicted cycle maxima (125–150), confirming non-extreme solar conditions.

As shown in right panels of Figure 3, Pt reached 2.1×10⁵ K at 20:00 UT, reflecting gradual heating from coronal mass ejection (CME)-driven solar wind compression, followed by fluctuations tied to magnetopause dynamics. Pd exhibited dual enhancements at 06:00 UT and 10:00 UT, temporally consistent with coronal hole high-speed stream arrivals. A sharp southward IMF-Bz deflection (-10.2 nT) at 06:00 UT persisted for ~2 hours, enhancing magnetopause reconnection efficiency and facilitating solar wind energy transfer to the magnetosphere. Concurrently, the eastward penetration electric field (Ey) surged to 4.92 mV/m at 06:00 UT, mirroring transient equatorial electrojet intensification during substorm expansion—a hallmark of activated magnetosphere-ionosphere coupling.

Multi-parameter analysis confirms a G1-class geomagnetic storm (06:00–08:00 UT) driven by coronal hole stream-IMF-Bz coupling, with an estimated energy injection rate of 1.3×10¹¹ W [66]—significantly weaker than intense storms (>5×10¹¹ W). These conditions provide a critical baseline to differentiate space weather-driven ionospheric perturbations from potential seismogenic signals.

3.3. Time-Frequency Analysis of Ionospheric TEC in the Seismogenic Zone

To systematically decode the spectral characteristics of ionospheric TEC disturbances before and after the 2023 Jishishan Ms6.2 earthquake, we conducted a multi-scale signal processing methodology using TEC time-series data from BDS GEO satellites (PRN: C01, C02 and C03) over the seismogenic region on 18 December 2023. Initially, an adaptive sliding-window smoothing algorithm was implemented to suppress ionospheric background noise, followed by time-frequency decomposition via continuous Morlet wavelet transform.

Figure 4 illustrates the wavelet power spectral density (PSD) distributions of TEC derived from three GEO satellites, with red dashed lines marking the earthquake rupture initiation. The pre-seismic ionospheric disturbances exhibit multi-band coupling characteristics, with dominant energy concentrations in three frequency bands: (1) high-frequency band: 5–30 min (0.56–3.33 mHz); (2) mid-frequency band: 30–60 min (0.28–0.56 mHz); (3) low-frequency band: 60–90 min (0.18–0.28 mHz).

To disentangle the frequency-dependent ionospheric perturbations, we implemented a tri-band spectral decomposition framework employing fourth-order Butterworth bandpass filters. These filters demonstrate maximally flat magnitude response within designated passbands, achieving amplitude deviation ≤0.1 dB across operational frequencies through optimized pole-zero placement in the Laplace domain [67]. The cutoff frequencies were rigorously optimized through wavelet global power spectrum analysis at a 95% confidence level, defining three distinct geophysical regimes. This frequency-domain partitioning ensured physical consistency in differential TEC (DTEC) extraction while suppressing harmonic artifacts through phase-matched filter design.

Therefore, the data processing workflow comprised three stages: (1) STEC series employing a GF dual-frequency linear combination algorithm; (2) background trend removal via a sliding median filter; (3) DTEC sequence extraction through a fourth-order Butterworth bandpass filtering.

Figure 5, Figure 6, Figure 7 and Figure 8 present DTEC time-series distributions across four GNSS stations spanning critical epicentral distances: near-field (5.2 km), mid-near-field (59.7 km), mid-far-field (151.6 km) and far-field (1,220 km). In high-frequency band (0.56–3.33 mHz), significant DTEC perturbations (amplitude ≤4 TECU) were detected across all stations, yet displayed spatially incoherent phase structures. The far-field station BJF1 showed a 75% amplitude reduction compared to near-field stations, quantitatively confirming the high-frequency energy attenuation coefficient (1.85 ± 0.12 dB/100 km) through inverse square law regression [68,69]. Notably, stations within the seismogenic zone maintained compared spectral power densities, suggesting localized ionospheric coupling mechanisms dominate over geometric spreading effects in the near-field regime.

The mid-frequency DTEC (0.28–0.56 mHz) time series exhibited regular oscillatory patterns, though amplitudes were only 1/4–1/10 of those in the high-frequency band. At the surface rupture onset, all stations showed a negative or positve-phase TEC perturbation, consistent with prior studies attributing such behavior to the anisotropic nature of ionospheric charged particles under geomagnetic influence. This anisotropy is similarly observed in seismic-induced acoustic-gravity waves (AGWs) during their generation and propagation. To explain in more detail, the IPPs of C02 clustered in the southwestern quadrant of the epicenter (taking LXJS as an example, the azimuth angle is 210 ± 0.92°.). Distinct perturbations (characterized by an N-shaped waveform) initiated 1–2 hours prior to surface rupture. On the contrary, the IPPs of C03 located in the southeastern quadrant displayed dominant post-seismic perturbations. Besides, significant disturbances (~0.3 TECU) were detected at station BJF1 (1,220 km), confirming mid-frequency TIDs can propagate over 1,200 km scales (see Figure 8).

The low-frequency regime (0.18–0.28 mHz) exhibited periodic differential total electron content (DTEC) fluctuations with amplitudes of 0.1–0.2 TECU, corresponding to 50–25% of mid-frequency perturbations. A distinct phase reversal phenomenon was observed, characterized by negative-polarity responses at LXJS Station (epicentral distance: 5.2 km) contrasting with positive-polarity coherence at distal stations (>150 km). This dichotomy arises from the interplay of two wave-coupling mechanisms:

(1)

Near-field acoustic wave dominance

Proximal to the epicenter, high-wavenumber acoustic waves dominate due to their rapid energy dissipation (<10-minute relaxation timescale). These waves induce localized ionospheric turbulence through supersonic pressure perturbations, generating transient electron density depletions that manifest as negative DTEC phases.

(2)

Far-field gravity wave sustenance

At distances exceeding 100 km, long-wavelength atmospheric gravity waves prevail, maintaining positive DTEC coherence via gradual electron density enhancements. Their lower attenuation coefficients enable sustained propagation, aligning with theoretical predictions for Lamb wave modes in thermospheric waveguides [70]. These waves exhibit horizontal phase velocities of 200–300 m/s, preserving phase coherence over megameter scales [71].

The spatial partitioning of wave regimes follows frequency-dependent ($f$) attenuation laws [72,73,74,75]:

$$\left\{\begin{array}{cc}{\gamma }_a\left(f\right)\propto f^{1.8}& \left(acoustic\right)\\ {\gamma }_g\left(f\right)\propto f^{0.3}& \left(gravity\right)\end{array}\right.$$

(14)

where $\gamma$ denotes attenuation rate. Near-field stations like LXJS primarily detect acoustic modes confined to F2-layer waveguides (<80 km altitude), while far-field observations sample gravity waves ducted in the lower thermosphere (90–150 km). This transition reflects fundamental differences in lithosphere-ionosphere coupling mechanics, where geometric spreading and atmospheric stratification jointly govern wave energy redistribution.

Table 2 quantifies the spatial evolution of DTEC across frequency bands and distinct epicentral distances. In the high-frequency band (0.56–3.33 mHz), near-field stations (LXJS, LXHZ, LZSH) exhibit substantial amplitude variability (Max–Min: 8.5–9.3 TECU), with high standard deviations (STD: 1.13–1.17 TECU) indicative of localized ionospheric turbulence. Conversely, the far-field station BJF1 demonstrates a pronounced amplitude attenuation (>75% reduction), exhibiting a constrained variation range (peak-to-peak: 2.03 TECU) and substantially diminished fluctuation intensity (STD: 0.25 TECU). These characteristics are indicative of efficient energy dissipation mechanisms operative in the distal propagation regime. The mid-frequency band (0.28–0.56 mHz) displays moderate amplitude variations (Max–Min: 0.75–2.79 TECU), with LXHZ station recoding the highest perturbation (2.79 TECU). The uniformly low STD valuses (0.10–0.25 TECU) suggests coherent gravity wave-driven propagation, enabling energy transfer over 1,200 km scales. This long-range coherence is attributed to the lower attenuation rates of gravity waves compared to acoustic modes.

In the low-frequency band (0.18–0.28 mHz), amplitude ranges decline nonlinearly with distance (2.62 TECU (LXJS) to 0.47 TECU (BJF1)), reflecting stress-driven coupling oscillations between the lithosphere and ionosphere. Near-field stations (LXJS, LXHZ) maintain statistically significant residuals (STD: 0.19–0.25 TECU), contrasting with the subdued variability at far-field receiver (STD: 0.06 TECU), a dichotomy attributable to the dominance of large-scale (>500 km wavelength) atmospheric gravity waves over local acoustic disturbances.

3.4. Propagation Velocity Analysis of Ionospheric Disturbances in the Seismogenic Zone

Utilizing multi-band DTEC disturbance signatures, this study quantitatively inverts the propagation parameters of traveling ionospheric disturbances (TIDs). A novel velocity inversion framework was developed through rigorous integration of GNSS station-satellite geometric configurations, systematically addressing IPP dynamics through three principal analytical stages: derivation of DTEC temporal variations across GNSS networks, geospatial mapping of IPP positions relative to the epicentral zone, and construction of time-distance matrices for wavefront correlation analysis.

Time-distance diagrams (Figure 9, Figure 10 and Figure 11) reveal systematic wave progression patterns through linear regression of spatiotemporal wave extrema. A positive slope in the fitted line indicates proximal propagation (disturbances moving away), while a negative slope denotes distal propagation (disturbances moving toward the epicenter), where the absolute slope value represents the characteristic velocity. For BDS GEO satellite (PRN: C02), the derived velocity parameters demonstrate azimuth-dependent characteristics, with positive regression slopes (proximal propagation) dominating in the southwestern quadrant (210–230° azimuth). Quantitative analysis yields propagation velocities of 210–270 m/s (95% CI) for 30–60 min TIDs and 175–260 m/s for 60–90 min TIDs. Notably, high-frequency perturbations (5–30 min) exhibit localized transient signatures with peak amplitudes up to 4 TECU, yet lack coherent phase progression sequences across adjacent stations, suggesting non-propagating ionospheric responses to rapid energy deposition.

The results validate three fundamental atmospheric gravity wave (AGW) trapping criteria established in lithosphere-ionosphere coupling theory [28]: (1) estimated wavelengths (52–324 km) compatible with E/F-region propagation constraints, (2) propagation velocities exceed the 130 m/s threshold for vertical trapping, and (3) periodicity below the 2-hour critical threshold. Remarkably, the derived velocities (150–250 m/s) exhibits strong concordance with medium-scale traveling ionospheric disturbances (MSTIDs), implying potential energy transfer via MSTID mechanisms in the seismogenic zone.

Besides, analytical results demonstrate distinct propagation regimes: high-frequency disturbances are localized and transient, while mid/low-frequency TIDs exhibit coherent wave-like propagation. Amplitude-frequency analysis reveals an exponential decay trend in DTEC perturbations with decreasing cutoff frequencies, consistent with AGW energy dissipation mechanisms in the thermosphere-ionosphere system. Azimuthal anisotropy analysis reveals significantly higher propagation velocities in the southwest direction (225 ± 15 m/s) compared to the southeast (195 ± 20 m/s), potentially reflecting geomagnetic-ionospheric coupling effects induced by stress field adjustments in the seismogenic zone.

4. Discussion

As revealed by multi-parameter cross-validation in Section 3.2, a G1-class minor geomagnetic storm occurred between 06:00–08:00 UT on the day of the Jishishan earthquake. The driving mechanism of this event can be attributed to the synergistic interaction between a coronal hole high-speed stream and the southward component of the interplanetary magnetic field (IMF). To disentangle the coupled effects of seismogenic processes and geomagnetic storm on ionospheric perturbations, we analyzed DTEC time series spanning ±2 hours relative to the storm interval. Bandpass-filtered results from different GNSS stations are illustrated in Figure 12, Figure 13, Figure 14 and Figure 15, with the left and right subpanels correspond to C02 and C03 observations, respectively.

As shown in the figures, within the seismogenic zone (LZJS, LZHZ, and LZSH), high-frequency DTEC fluctuations (amplitudes reaching up to 4 TECU) were observed during the stress accumulation phase of the earthquake, superimposed with weak geomagnetic storm activity. While the DTEC time series exhibited quasi-periodic patterns, no statistically dominant cycles were identified (Figure 12, Figure 13 and Figure 14). The absence of coherent phase propagation sequences suggests localized transient processes (e.g., abrupt stress releases or plasma instabilities) rather than sustained external energy input. The mid-frequency DTEC components (0.28-0.56 mHz) exhibited pronounced geomagnetic-ionospheric coupling manifested through characteristic N-shaped waveforms, with C02 satellite observations (IPPs clustered southwest of the epicenter) showing higher-amplitude DTEC fluctuations compared to C03 (IPPs concentrated southeast). This azimuthal asymmetry persists under uniform geomagnetic forcing conditions, suggesting preferential stress alignment along the southwestern fault strike in the seismogenic zone. Low-frequency DTEC signals (0.18-0.28 mHz) from all near-field stations displayed characteristic quasi-sinusoidal oscillations with low intensities (<0.3 TECU). In contrast, far-field regions (Figure 15) exhibited reduced DTEC amplitudes across all frequency bands due to the limited AGW propagation range and subdued geomagnetic activity.

As established in prior studies, geomagnetic storms trigger enhanced Joule heating in the polar regions and thermospheric atmospheric expansion through solar wind energy injection and magnetosphere-ionosphere coupling processes [76,77]. The expansion propagates equatorward along magnetic field lines, leading to a global elevation of the ionospheric F-layer height and significant reorganization of electron density gradients [78,79]. Ionospheric uplift further modifies AGW dynamics, increasing horizontal phase velocities and extending wavelengths to 150–300 km scales with prolonged propagation distances [80].

Besides, during geomagnetic storms, resonant coupling between field-aligned currents and fluctuating magnetic field components in the ionosphere induces phase-synchronized modulation. Concurrently, penetration electric fields driven by interplanetary magnetic field reorientation interact with atmospheric gravity waves, amplifying quasi-periodic fluctuations in ionospheric electron density (Figure 12, Figure 13, Figure 14 and Figure 15). Compared to the effects of seismic activity alone, the superposition of geomagnetic storms leads to increased amplitude of differential total electron content (DTEC) in far-field regions, particularly for medium-to-low frequency DTEC (0.18~0.56 mHz), which exhibit enhanced quasi-sinusoidal oscillations.

The frequency-dependent interplay between seismogenic AGWs and weak geomagnetic storms is further elucidated through travel-time analysis (Figure 16, Figure 17 and Figure 18). High-frequency DTEC disturbances (5–30 min) in the seismogenic zone (stations LZJS, LZHZ, and LZSH) exhibited fragmented waveforms (Figure 16) with amplitudes comparable to non-storm conditions (Figure 9), despite superimposed geomagnetic activity. This amplitude invariance suggests high-frequency perturbations primarily originate from localized transient processes rather than sustained energy input from external geomagnetic drivers. The absence of coherent phase propagation sequences further supports the dominance of non-propagating modes in high-frequency disturbances, distinct from gravity wave-driven mechanisms observed in lower-frequency bands.

In contrast to the high-frequency components, mid-frequency DTEC signals (0.28–0.56 mHz) permit robust extraction of characteristic velocities (Figure 17). The upper velocity limit (280 m/s) aligns temporally with the seismic event, consistent with theoretical predictions for LAI coupling. However, the lower velocity bound (80 m/s) reflects significant suppression under weak-to-moderate geomagnetic activity (Figure 10), likely due to storm-induced ionospheric currents disrupting AGW coherence. This dual-range behavior highlights the modulatory role of geomagnetic storms on AGW energy transfer efficiency.

The modulation effects of weak geomagnetic storms are most pronounced in the low-frequency band (0.18–0.28 mHz). For satellites with IPPs in distinct quadrants, the lower bounds of characteristic DTEC velocities dropped below 100 m/s, indicating ubiquitous reductions in wave coherence (as shown in Figure 18). This suppression aligns with theoretical models of geomagnetic storm-AGW interactions, where perturbations in Pedersen and Hall conductivities disrupt neutral-ion momentum coupling—a critical mechanism for sustained wave energy transfer. Notably, these velocities fall below the typical range of lithosphere-driven AGWs, emphasizing the dominance of external geomagnetic forcing over tectonic processes in shaping low-frequency disturbances.

5. Conclusions

This study employs multi-frequency BDS GEO observations within the Jishishan seismic preparation zone to estimate ionospheric TEC parameters along the GNSS-satellite line-of-sight during the 2023 Ms6.2 earthquake. The DTEC sequences were extracted using a Butterworth filter, enabling a comprehensive analysis of disturbance characteristics from perspectives of amplitude, waveform, spatiotemporal distribution, propagation velocity/direction, and time-frequency domain.

This study systematically decoupled ionospheric disturbances (DTEC) across three frequency bands (0.18–0.28, 0.28–0.56, and 0.56–3.33 mHz) using optimized Butterworth filters validated by wavelet global power spectrum analysis (95% confidence). High-frequency DTEC (0.56–3.33 mHz) exhibited localized perturbations (≤4 TECU) within 300 km of the epicenter, with 75% amplitude attenuation at 1,220 km (BJF1 station), consistent with acoustic wave dominance and strong spatial decay. Mid-frequency signals (0.28–0.56 mHz) demonstrated anisotropic propagation (>1,200 km range), showing coherent “N-shaped” waveforms in southwestern sectors pre-seismicity and post-seismic disturbances in southeastern quadrants, indicative of gravity wave mediation. Low-frequency DTEC (0.18–0.28 mHz) revealed non-linear attenuation and phase reversals (negative near-field, positive far-field), suggesting lithosphere-ionosphere coupling delays during stress adjustments. Spatial evolution (Table 1) highlights frequency-dependent mechanisms: high-frequency rapid decay (localized), mid-frequency long-range transmission (gravity waves), and low-frequency phase delays (coupling oscillations). These findings provide critical constraints for distinguishing seismic precursors from ambient ionospheric noise.

Besides, this study quantifies the propagation velocities of TIDs using a GNSS geometry-based inversion model, revealing distinct frequency-dependent characteristics. High-frequency disturbances (5–30 min) exhibited localized transient responses with amplitudes up to 4 TECU but lacked coherent phase propagation. Mid- to low-frequency TIDs (30–90 min) showed anisotropic southwestward propagation (210–270 m/s for 30–60 min; 175–260 m/s for 60–120 min), consistent with AGW dynamics. Derived velocities (150–250 m/s) and wavelengths (52–324 km) align with AGW trapping criteria for E/F-region propagation. Azimuthal anisotropy analysis identified higher velocities in the southwest (225±15 m/s) than southeast (195±20 m/s), likely reflecting stress-induced geomagnetic-ionospheric coupling. Exponential amplitude decay with frequency corroborates AGW energy dissipation, while mid/low-frequency TIDs exhibited uniform proximal propagation patterns, suggesting MSTID-mediated energy transfer. These findings confirm seismogenic zone perturbations as primary drivers of multi-scale ionospheric disturbances.

Multi-parameter analysis reveals that the G1-class geomagnetic storm (06:00–08:00 UT), driven by coronal hole high-speed streams and southward IMF, superimposed frequency-dependent perturbations on seismogenic ionospheric disturbances. High-frequency DTEC (5–30 min) exhibited localized transient amplitudes ≤4 TECU during stress accumulation, with fragmented waveforms lacking coherent phase propagation, indicating dominance by abrupt stress releases rather than storm-driven energy. Mid-frequency signals (30–60 min) showed azimuthally asymmetric N-shaped waveforms (C02: southwest, 180±30 m/s; C03: southeast), correlating with seismogenic stress orientation. Low-frequency oscillations (60–90 min) maintained sinusoidal regularity (<0.3 TECU), while far-field attenuation confirmed limited AGW propagation. Despite storm-ionosphere coupling, high-frequency DTEC amplitudes remained comparable to non-storm conditions, underscoring localized seismogenic origins over external drivers. The absence of velocity amplification and coherent wavefronts further validates the predominance of crustal processes in high-frequency ionospheric perturbations.

It should be noted that our investigation into pre-seismic ionospheric TEC anomalies reveals systematic discrepancies from prior research, primarily stemming from the following methodological distinctions. Prior research predominantly emphasizes statistical anomaly characterization and phenomenological descriptions, utilizing GPS-derived TEC data with hourly resolution over extended temporal scales (days to months). In contrast, this study employs a process-driven framework, analyzing 30-second-resolution TEC data from BDS GEO satellites to investigate AGW dynamics during the immediate pre-seismic phase (0–6 hours). The geostationary BDS GEO system eliminates spatial interference caused by IPP variability inherent to medium-Earth-orbit GPS satellites, thereby enhancing observational fidelity.

Besides, Surface-originated perturbations (e.g., vertical crustal vibrations, Rayleigh waves, and potential tsunami waves) generate upward-propagating AGWs. Acoustic modes induce TEC perturbations with ~10-minute latency, while gravity waves produce oscillatory TEC anomalies after ~60 minutes. Energy conservation dictates exponential amplitude growth during atmospheric ascent, enabling ionospheric detection of even weak preseismic crustal motions. However, spatiotemporal correlations between tectonic deformation patterns (e.g., slow-slip events minutes to months pre-failure and AGW-driven TEC responses remain unresolved due to nonlinear atmospheric-ionospheric coupling.

The LAI energy transfer involves complex nonlinear interactions, obscuring causal linkages between surface deformation signatures and observed TEC disturbances. Quantifying triggering thresholds and coupling efficiencies requires advanced mechanistic models integrating multi-layer physical parameters. To resolve these challenges, we propose: Deploy integrated networks combining ionosondes, Doppler radars, and GNSS receivers to track AGW propagation dynamics across atmospheric strata with sub-minute temporal resolution. Establish cross-parameter frameworks correlating TEC anomalies with radon emissions, infrared radiation, and seismic waveforms to disentangle seismogenic signals from geomagnetic storm-modulated noise.