# The First Near-Real-Time Compatible Numerical Method for Co-Seismic Ionospheric Disturbances Simulation

S. A. Sanchez, E. A. Kherani, E. Astafyeva, E. R. Paula

► To cite this version:

S. A. Sanchez, E. A. Kherani, E. Astafyeva, E. R. Paula. The First Near-Real-Time Compatible Numerical Method for Co-Seismic Ionospheric Disturbances Simulation. Geophysical Research Letters, 2025, 52, 10.1029/2025GL115208. insu-05080906

HAL Id: insu-05080906

https://insu.hal.science/insu-05080906v1

Submitted on 23 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geophysical Research Letters®

## RESEARCH LETTER

10.1029/2025GL115208

**Key Points:**

* New fast simulation of Acoustic Gravity Waves and co-seismic ionospheric disturbances (CSID, or ionoquakes)

* In less than a minute, the simulation reproduces near-field CSID, that is, about 5–7 min before they are detected in observational data

* This method can be used in near-real-time (NRT) to confirm detection of CSID and for future NRT ionosphere-based assessment of earthquakes

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

S. A. Sanchez,
saul.juarez@inpe.br;
saulsjf@gmail.com

**Citation:**

Sanchez, S. A., Kherani, E. A., Astafyeva, E., & de Paula, E. R. (2025). The first near-real-time compatible numerical method for co-seismic ionospheric disturbances simulation. *Geophysical Research Letters, 52*, e2025GL115208. https://doi.org/10.1029/2025GL115208

Received 5 FEB 2025
Accepted 1 MAR 2025

© 2025. The Author(s).
This is an open access article under the
terms of the Creative Commons
Attribution License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.

<figure><img src="image_2.png" /></figure>

## The First Near-Real-Time Compatible Numerical Method for Co-Seismic Ionospheric Disturbances Simulation

S. A. Sanchez¹ ID, E. A. Kherani¹ ID, E. Astafyeva² ID, and E. R. de Paula¹ ID

¹Instituto Nacional de Pesquisas Espaciais, São José dos Campos, Brazil, ²Université Paris Cité, Institut de physique du globe de Paris (IPGP), CNRS UMR, Paris, France

**Abstract** Numerous studies demonstrated a huge potential for ionospheric total electron content (TEC) measurements to be included in earthquake and tsunami risk assessments. However, up to now, only a few methodologies or tools can be used in near-real-time (NRT) for these purposes. For the first time, this study presents a method allowing for fast simulations of Co-Seismic Ionospheric Disturbances (CSID) or ionoquakes. The new analytical simulation code models the propagation of acoustic-gravity waves (AGWs) generated by a co-seismic uplift of the ground/seafloor, by resolving the governing equations in the time-altitude-horizontal domain. The method models the near-field CSID in about 60 s of simulation time, that is, before the waves are detected in the ionospheric TEC measurements. The developed method is, therefore, among the first most promising products in Ionospheric Seismology that can be used for NRT applications, such as tsunami early warning.

**Plain Language Summary** This study presents a new method for fast simulation of disturbances generated by earthquakes that are often referred to as Co-Seismic Ionospheric Disturbances (CSID). Using a combined analytical and numerical simulation code, our tool models acoustic gravity waves generated by seismic ground or seafloor uplift. The model can simulate these CSID in approximately 60 s. The method has been validated on significant earthquakes such as the Mw7.3 2011 Sanriku-Oki and the Mw7.5 2024 Noto-Peninsula earthquakes, showing fast and accurate results compared to real observations. This advancement enables near-real-time use for early tsunami detection, providing a valuable tool to assess risks associated with earthquakes and tsunamis.

## 1. Introduction

Earthquakes and tsunamis are the most disastrous natural hazards; every year, they cause numerous casualties and significant economic loss. The ability to detect and warn of the arrival of tsunamis as early as possible is crucial to minimize the impact on human life. The tsunami early warning systems depend on a rapid and accurate estimation of a tsunami potential of an earthquake, which is based on the amplitude of the co-seismic uplift of the Earth's crust and the extent of the source, which is difficult to reliably obtain in near-real-time (NRT). Conventional seismological and seismo-geodetic approaches still fail to correctly estimate the tsunamigenic potential for large earthquakes (LaBrecque et al., 2019). In response to this need, the ionospheric seismology offers a promising alternative. Earthquakes are known to produce disturbances in the ionosphere, referred to as co-seismic ionospheric disturbances (CSID) or ionoquakes (Afraimovich et al., 2001; Astafyeva, 2019; Bagiya et al., 2023; Cahyadi & Heki, 2015; Calais & Minster, 1995; Heki et al., 2022; Jin et al., 2015; Occhipinti et al., 2013; Rolland et al., 2013; Sanchez et al., 2023; Zedek et al., 2021). The CSID manifest themselves as acoustic-gravity, acoustic, and shock-acoustic waves in the ionospheric radio sounding measurements (Astafyeva et al., 2011; Chum et al., 2012; Heki & Ping, 2005; Rolland et al., 2011). The CSID can be detected as early as 480 s after the earthquake onset time (Astafyeva & Shults, 2019; Cahyadi & Heki, 2015). The CSID parameters depend on parameters of earthquakes and tsunamis (e.g., Astafyeva, Shalimov, et al., 2013; Astafyeva, Rolland, et al., 2013; Kamogawa et al., 2016; Manta et al., 2020; Rakoto et al., 2018; Ravanelli et al., 2021). All this opens a possibility for estimation of earthquake parameters and, therefore, for assessment of the tsunamigenic potential of earthquakes from the ionosphere in NRT.

Going toward this development, the first step is a rapid and reliable detection of CSID in NRT (e.g., Brissaud & Astafyeva, 2022; Fuso & Ravanelli, 2024; Maletckii & Astafyeva, 2021, 2022; Martire et al., 2023; Ravanelli et al., 2024; Ten et al., 2025). The second step is confirming the seismic origin of the detected disturbances that
can be done by estimating the spatio-temporal parameters of CSID in NRT (Maletckii & Astafyeva, 2024), and/or by fast modeling. However, up to now, CSID modeling by using numerical simulations remains time-consuming (e.g., Inchin et al., 2021; Inchin et al., 2024; Kherani et al., 2012; Meng et al., 2015; Mikesell et al., 2019; Nozuka et al., 2024), and, therefore, they are not applicable in NRT. This work presents the first method for rapid analytical simulation of CSID that could be applicable in NRT.

## 2. Methods

### 2.1. Earthquake Events

The Sanriku-Oki earthquake of Mw7.3, occurred on 9 March 2011, at 02:45:20 UT at a depth of 32 km, with an epicenter at 38.435°N, 142.842°E (https://earthquake.usgs.gov/earthquakes/eventpage/usp000hvjh/executive), due to shallow thrust faulting near the subduction zone interface between the Pacific and North American plates.

The Noto-Peninsula earthquake of Mw7.5, occurred on 1 January 2024, at 07:10:09 UT at a depth of 10 km, with an epicenter at 37.487°N, 137.271°E (https://earthquake.usgs.gov/earthquakes/eventpage/us6000m0xl/executive). It occurred from shallow reverse faulting in the Earth's crust. Both events generated noticeable CSID (Astafyeva & Shults, 2019; Heki, 2024; Thomas et al., 2018).

### 2.2. Simulation Input for the Source

We employed the ObsPy library from PyPI to process the recorded ground motions in count format available at the National Research Institute for Earth Science and Disaster Resilience (NIED)-F-net (https://www.fnet.bosai.go.jp/fnet/event/joho.php?LANG=en). The module downloads the data in this format and estimates the vertical ground velocity or uplift ($U_{sism}$) by applying instrumental response corrections. These corrections include the removal of the frequency response of the seismometers, as well as the effects of any amplifiers, analog and digital filters, and digitization. For both earthquakes, we used a band-pass filter on the seismograms before deconvolution, defined by a tuple specifying the four cut-off frequencies: 0.02, 0.12, 10, and 20 Hz (Hosseini & Sigloch, 2017).

Figures 1a and 1b demonstrate the $U_{sism}$ at the epicenter for the Sanriku-Oki and Noto-Peninsula earthquakes, respectively. The oscillations in the seismic vibrations in Figure 1a show a higher frequency compared to Figure 1b, even though the two-time series have a resolution of 7.5 s. It can be seen that the overall energy release of the seismic event captured in Figure 1a is more concentrated or intense, depending on the initial peak and rapid decrease in amplitude. The seismic event in Figure 1b shows a longer or sustained energy release, based on the longer duration of the vibrations before subsiding. This suggests that the response in generating acoustic-gravity waves (AGWs) differs, which also affects the generation of CSID.

Because the location and the magnitude of the uplift that is, the epicentral distribution of $U_{sism}$ is not known in NRT, in this first version of our model, we construct the simulation input from two main observational parameters: (a) the $U_{sism}$ vibrations, obtained from seismometers following the method described by Sanchez and Kherani (2024) and placing the strongest vibration at the epicenter; (b) the fault size Gaussian variation on either side from the epicenter up to the fault size. The fault size is taken from the USGS models. In Figures 2c and 2d, the $U_{sism}$ variation in time and epicentral distance for the two earthquake events are shown.

### 2.3. CSID Simulation

The simulation is based on the Seismo-Atmosphere-Ionosphere coupling mechanism and energetics that include the generation of AGWs from the seismic uplift at the ground and the subsequent generation of CSID. Our approach combines analytical and numerical methods: AGWs are modeled analytically, while the TEC simulation is performed numerically. Balachandran (1968) and Tahira (1995) had obtained an analytical solution of amplitude of non-viscous AGWs of the following form:

$$u_y^{BC} = e^{i\omega t - ik_x x} (ae^{ik_y y} + be^{-ik_y y})e^{-k_0 y} \quad (1)$$

$$\omega^4 - (k_x^2 s^2 + k_y^2 s^2 + \Omega_c^2)\omega^2 + k_x^2 s^2 \Omega_b^2 = 0 \quad (2)$$

$$s^2 = \frac{\gamma P}{\rho}, \quad \Omega_c = \frac{\gamma}{2}k_0 s, \quad k_0 = -\frac{1}{2} \frac{\gamma g}{s^2}, \quad \Omega_b = \left[ (\gamma - 1)k_0^2 - \frac{k_0}{s^2} \frac{ds^2}{dy} \right]^{1/2} s$$
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a–b) Show the ground uplift, in terms of vertical velocity $U_{\text{sism}}$ at the epicenter for the Sanriku-Oki and Noto earthquakes, respectively. (c–d) Show $U_{\text{sism}}$ distribution with epicentral distance for the two earthquakes, respectively. (e) Shows the sound speed profiles from the NRLMSISE00 model and (f) shows the ionospheric electron number density profiles from the IRICORE model.</figcaption></figure>

here, (t, y, x) represent the time, altitude and epicentral distance, (ω, kx, ky) are the wave frequency and wave-numbers in the altitude-horizontal domain, (a, b) are the constants associated with the boundary conditions, (Ωc, Ωb) are the acoustic cut-off and Brunt-Vaisala frequencies, (γ, s) are the specific heat ratio and sound speed.
Similar to the above solution, Sanchez and Kherani (2024) have derived the vertical and horizontal amplitude of viscous-AGWs as shown in solutions 3 and 4, respectively. The viscosity model is the same as those used by Campos (1983), Kherani et al. (2012, 2016):

$$u_y(t, y, x) = \sum_{k_x, k_y} U_{\text{sism}} \exp(-i\omega t + ik_y y + ik_x x) \exp\left[-\frac{k_x^2 + k_y^2}{2}\nu t\right] e^{-\mu} \quad (3)$$

$$u_x(t, y, x) = \sum \frac{\omega^2 - k_y^2 s^2 - \Omega_b^2}{k_x k_y s^2} u_y \quad (4)$$

$$k_0 = \frac{1}{s^2} \frac{1}{\rho} \frac{dp}{dy}, \quad \mu = \int k_0 dy,$$

here, $U_{\text{sism}}$ is the ground uplift associated with the seismic vibration, and $\nu$ is kinematic viscosity.

A few important differences between the solution $u_y^{BC}$ in Equation 1 and $u_y$ in Equation 3 are that the solution Equation 3 considers the non-hydrostatic atmosphere and retains the integral of $k_0$ with altitude via $\mu$ that is, it accounts for the non-local effects of atmospheric density and temperature variations with altitude. In Equation 3, the integral $\mu$ is the non-dimensional scaling factor that arises from the fact that atmospheric flow tends to scale with the characteristic scale, that is, the scale height $(1/k_0)$ (equation 3 of Tahira (1995)). Since the atmosphere density/temperature have different scale heights at different altitudes, the factor ($\mu$) allows the waves to have all possible wavelengths ranging between the lowest scale height to the largest scale height (Sanchez & Kherani, 2024; Tahira, 1995). The factor $e^{-\mu}$ increases exponentially with altitude and is responsible for the amplitude growth of AGWs by 3–4 order at thermospheric altitudes from the tiny ground uplift $U_{\text{sism}}$. Solution 4 maintains $\nabla \times u = 0$ for a pure acoustic wave that is, for $\omega = \text{ks}$, $\Omega_b = 0$.

The simulated TEC and CSID amplitude ($\Delta TEC_{\text{sim}}$) are estimated with the following methods:

$$\text{TEC}_{\text{sim}}(t, y, x) = \int_{150\text{km}}^{300\text{km}} n(x, y, t) dy, \quad (5)$$

$$\Delta \text{TEC}_{\text{sim}}(t, y, x) = \text{TEC}_{\text{sim}}(t, y, x) - \text{TEC}_{\text{sim}}(t_p, y, x) \quad (6)$$

here (n) is the ionospheric density, (t, t_p) are the current and previous times, respectively. The ionospheric motions and currents due to the wind amplitude ($u_y$, $u_x$) of AGWs cause the electron density change governed by the following continuity equation:

$$\frac{\partial n}{\partial t} + \nabla \cdot \left[ n \bar{\eta} \cdot \left( \vec{E} + \vec{u} \times \vec{B}_0 + \vec{u} \cdot \vec{B}_0 \frac{\vec{B}_0}{|\vec{B}_0|} \right) \right] = 0 \quad (7)$$

where $(\bar{\eta}, \vec{B}_0, \vec{E})$ are the electron mobility tensor, geomagnetic field strength, and electric field caused by current due to the wind amplitude ($\vec{u}$) of AGWs. A detailed set of equations and numerical methods are described by Kherani et al. (2012) and outlined in the Text S1 in Supporting Information S1.

The simulation domain covers 0–300 km (150–300 km) altitudes in the atmosphere (ionosphere) and ±500 km in epicentral distance with $x = 0$ at the location of the epicenter. The simulation begins a few minutes before the earthquake onset time ($t_{ep}$) and spans for 30 min with a time resolution of 7.5 s. In the analytical form of the solution Equation 3, the integral in $\mu$ is to be evaluated numerically with discretization in space. For a chosen range of wavenumbers ($\frac{\pi}{N\Delta} \le k_y \le \frac{\pi}{\Delta}$, $\frac{\pi}{N\Delta} \le k_x \ll k_y$), where ($\Delta = 10$ km, N = 30) are the grid resolution and number of grids in altitude, simulation estimates ($\omega, u_y, u_x$) from Equations 2–4 and $\Delta TEC$ from Equation 6. At $t = 0$ seconds, the atmospheric ($\rho, T, p$) and ionospheric parameters are inferred from the NRLMSISE00 and IRICORE (Python library of PyPI) models. Figures 1e and 1f demonstrate the atmospheric sound and ionospheric density profiles above the epicenter for the Sanriku-oki and Noto-Peninsula earthquakes.
# 3. Results and Discussions

## 3.1. Simulation Results

Figures 2a and 2c demonstrate the AGWs amplitude during their upward propagation above the epicenters of the 2011 Sanriku-oki and the 2024 Noto-Peninsula earthquakes, respectively. In both cases, from the ground uplift, numerous waves with wavefronts of different slopes or phase speeds propagate into the atmosphere. The AGWs attain the largest amplitude of about 50 m/s in mesosphere-lower-thermosphere, owing to the amplification factor $e^{-\mu}$ and ducting mechanism (Kherani et al., 2011). AGWs also attain a secondary maximum (denoted by the black circle) in the upper thermosphere. In the case of the Sanriku-Oki earthquake, the secondary maximum is reached at 180 km altitude at 02:52 UT, that is, the waves with notable amplitude arrive in the ionosphere in about 400 s from the onset time of 02:45:20 UT. Similarly, in the case of the Noto earthquake, notable amplitude waves reach an altitude of 200 km (the secondary maxima denoted by the black circle) at 07:18 UT, which means that they take about 490 s from the onset time of 07:10:09 UT to reach the ionosphere (Figure 2c).

Figures 2b and 2d show the electron density disturbances at different ionospheric altitudes above the epicenter of the Sanriku-Oki and the Noto-Peninsula earthquakes, respectively. Figure 2b shows the formation of CSID of significant amplitude from 02:52 UT onwards. We note that the CSID attains the largest amplitude at around 210 km altitude, which is highlighted in black. Therefore, the most rapid CSID in TEC measurements should originate at about 210 km altitude of a majority of CSID, close to the origin altitude of 190 km reported by Astafyeva and Shults (2019). In Figure 2d, we note the formation of CSID from 07:18 UT onwards above the epicenter of the Noto earthquake. The strongest CSID (highlighted in black) originates at around 190 km altitude, as evident in Figure 2d.

## 3.2. Comparison With the Observations

### 3.2.1. Mw7.3 Sanriku-Oki Earthquake of 9 March 2011

Figure 3a presents CSID signatures in TEC measurements along the epicentral distance, as registered by five GNSS-stations and the GPS satellite G07. These satellite-receiver pairs are the same as those considered in previous studies (Astafyeva & Shults, 2019; Thomas et al., 2018). We note the CSID detection at 2:52-2:54 UT at an epicentral distance of 0.5°–1.75°, similar to those reported from previous studies. The green line with a slope of ~1 km/s connects the peak of CSID and suggests the CSID propagation speed of ~1 km/s which is in agreement with the speed reported by Astafyeva and Shults (2019).

Comparison is the simulation results with the observations is shown in Figure 3b, and in Figures S1a and S2a in Supporting Information S1. We note similar waveform characteristics of the first CSID oscillation from observation and simulation. The differences between the two are attributed partly to the non-seismic contributions in the TEC data, which are absent in the simulations, and partly to the linear nature of the analytical simulation that excludes the secondary wave generation. Also, the TEC estimation method from simulated electron density differs from the TEC estimation from the observations. The observed TEC is calculated by combining carrier phase and code measurements from dual-frequency GNSS signals (Hofmann-Wellenhof et al., 2008), representing the number of electrons along the line of sight between a satellite and a ground-based receiver. In contrast, the simulated TEC is obtained by Equation 5. Figure 3c, Figures S1b and S2b in Supporting Information S1 present the analysis of observed and simulated CSID detection time, considering the first peak of the CSID, as suggested by Astafyeva and Shults (2019). The simulation estimates the detection time to be around 480 s at the epicentral distance of about 0.5°, confirming the observed detection time by Astafyeva and Shults (2019). Considering the detection time method of Thomas et al. (2018), Figure S1c in Supporting Information S1 presents the CSID detection time defined as the beginning time of the first CSID oscillation since the onset of the earthquake. The observed detection time for station 0940 PRN G07, at the epicentral distance of 0.5°, is approximately 420 s, a value close to that reported by Thomas et al. (2018) for the same receiver-satellite pair.

### 3.2.2. Mw7.5 Noto-Peninsula Earthquake of 01 January 2024

Figure 4a presents TEC signatures of CSID generated by the Noto-Peninsula earthquake recorded by five GNSS-stations and satellite G16. The first CSID are observed at an epicentral distance of 0.5° about 550 s after the earthquake onset time, considering the first peak.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Observation and Simulation comparison for the Sanriku-Oki earthquake: (a) The observed ΔTEC disturbances that is, CSID from five receiver-G07 pairs where numbers denote receiver stations at the top of each curve, (b) The ΔTEC disturbances that is, CSID from observation (red curve) and simulation (black curve), (c) The CSID detection time estimation from observation (red-square) and simulation (black circle).</figcaption></figure>

Figure 4b and Figures S3a and S4 in Supporting Information S1 present the simulation results and observations for CSID recorded by different satellites. Similar features can be seen in the waveform of the first oscillation of the CSID in both observations and simulations. Figure 4c shows the estimated detection time of CSID, considering the first peak of the CSID, both observed and simulated. The simulation estimates the detection time of the CSID at less than 550 s at an epicentral distance of approximately 0.5°, confirming the observed detection time. Additionally, in Figure S3b in Supporting Information S1 with a greater number of measurements, the simulation estimates the detection time of about 500 s considering the first maximum peak, and the detection time of about 410 s in Figure S3c in Supporting Information S1, considering the detection time method of Thomas et al. (2018), in both cases just above the epicenter. To reinforce the validity of our simulations, in addition to GPS data, we present the results for GLONASS (R09) and Galileo (E31) satellites (Figures S5 and S6 in Supporting Information S1), respectively. These results provide additional confirmation of the robustness of our findings.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Observation and Simulation comparison for the Noto earthquake: (a) The observed ΔTEC disturbances that is, CSID from five receiver-G16 pairs where numbers denote receiver stations at the top of each curve, (b) The ΔTEC disturbances that is, CSID from observation (red curve) and simulation (black curve), (c) The CSID detection time estimation from observation (red-square) and simulation (black circle).</figcaption></figure>

As shown at Figure 1f, during the Sanriku-Oki event, the F region of the ionosphere (i.e., the maximum electron density), is located at a slightly higher altitude than that during the Noto earthquake. This lifts the altitude of the strongest CSID to higher altitudes despite the lower altitudes of the strongest AGWs on the day of Sanriku-Oki earthquake. Therefore, the amplitude of AGWs and the ionospheric layer height collectively decide the altitude of the strongest CSID.

### 3.2.3. Rapid Coupling Energetics From Fast AGWs

The rapid simulation results in Figures 2–4 reveal a good agreement between the observations and the simulations. The fast-propagating AGWs play a crucial role in rapid CSID formation. The dynamics of the coupling include the amplification of tiny ground vibrations by a factor of about $10^3$ due to the amplification factor $e^{-\mu}$ in Equation 1.

Based on AGWs simulation, the present study validates the observed CSID detection time in 420–470 s, earlier than the “typical” 480–600 s after an earthquake. Therefore, fast AGWs energetics present numerous possibilities for rapid CSID formation. The first arrival time of the AGWs at an altitude closely follows the acoustic ray travel
time (shown as green dots in Figures 2a and 2c). Therefore, the AGWs arrive at an altitude with phase speeds comparable to the sound speed at that altitude. As known, numerous waves with different phase speeds propagate into the atmosphere generated by the ground seismic uplift (e.g., (Sanchez et al., 2023; Sanchez & Kherani, 2024)). It is owing to the various scale heights and ducts of varying sizes present in the sound speed profile that allows numerous wavelengths and phase speeds at a given frequency to sustain in the atmosphere. The longer wavelengths connect to the thermospheric altitudes and arrive at those altitudes with the thermospheric sound speed, that is, they arrive at thermospheric altitudes earlier than expected from integrated or averaged sound speed. We note that such earlier arrivals are in agreement with much more advanced simulation tools (Inchin et al., 2024; Nozuka et al., 2024).

### 3.3. Seismic Source Approximation

Despite a simplified source on the ground, the present study demonstrated a close agreement between observations and simulations. This means that the two important inputs, the time history of the $U_{SISM}$ and the fault size used to define the source in the present study, are adequate to reproduce realistic CSIDs even before the CSID are actually detected in the ionosphere. Other detailed information such as the distribution of the ground vibration in fault regions and shake maps might enhance the degree of agreement in observation and simulation. However, the present study does not rely on detailed source energetics while reproducing CSID.

For future NRT applications, we plan to reduce dependency on source parameters that are often unavailable in NRT in seismic measurements. Instead, we are planning to only use $U_{SISM}$ that can be estimated in NRT either from seismic stations, or, eventually, from ground-based high-rate (HR) GNSS receivers Melgar and Hayes (2019).

### 3.4. The Use of the Simulation Tool for Future NRT Tsunami Early Warning

When an earthquake occurs, its magnitude and epicenter location can be estimated from seismic data in about 2–5 min after the earthquake onset time. However, these parameters are not enough to correctly assess the tsunamigenic potential of an earthquake. The most crucial parameters for tsunami early warning are the co-seismic displacement of the ground or the sea floor, as well as its location and the horizontal dimensions. Recent advances in seismo-geodesy, with a combined use of seismic data and HR- GNSS data of co-seismic ground displacements, allow to estimate earthquake size 5–15 min after the earthquake onset time (e.g., Kawamoto et al., 2016, 2017; Melgar & Bock, 2015; Melgar & Hayes, 2019; Ohno et al., 2021). However, the results of seismo-geodetic solutions largely depend on the position of GNSS-stations with respect to a seismic fault’s location and dimensions. Besides, HR-GNSS coverage remains limited in offshore subduction zones, where real-time near-field data are often unavailable (Melgar & Hayes, 2019). Therefore, the NRT assessment of tsunamigenic potential remains quite challenging, especially for large earthquakes with magnitudes above 8 (e.g., LaBrecque et al. (2019)). In response to this need, it has been suggested that the tsunami potential can be estimated from the ionosphere, owing to the close link between the earthquake source parameters and the CSID (e.g., Astafyeva et al., 2011; Astafyeva, Shalimov, et al., 2013; Astafyeva, Rolland, et al., 2013; Bagiya et al., 2020; Manta et al., 2020; Sunil et al., 2021).

Once the CSID are detected in NRT, we need to confirm that a detected disturbance is related to the earthquake, since many other disturbances of other origins can propagate through the ionosphere at any time. The first attempt to identify the source of ionospheric disturbance has been made recently by (Maletckii & Astafyeva, 2024). In addition to the source identification, CSID can be rapidly confirmed by using fast modeling tools, which did not exist until now. This work presents the first tool allowing to rapidly simulate CSID in the near-epicentral area, which is the first step toward NRT assessment of tsunamigenic potential of earthquakes from the ionosphere. The steps after will include estimation of an earthquake's magnitude from the ionosphere, and ionospheric imaging of seismic faults, which are the subjects of our future works.

## 4. Summary

This study presents a new innovative combined analytical and numerical code to simulate co-seismic ionospheric disturbances (CSID) or ionoquakes in NRT. This approach uses analytical code to model the propagation of acoustic-gravity waves generated by an earthquake and numerical simulation to calculate ionospheric density, providing results in around 60 s, that is, much faster than the CSID detection time by GNSS sounding. The
approach was validated with the 2011 Sanriku-Oki and 2024 Noto-Peninsula earthquakes. Both produced detectable CSID, and the simulation replicated similar waveform characteristics of the first oscillation. Other existing numerical simulation tools can reproduce CSID very well; however, they are slower for NRT applications. For the Sanriku-Oki (Noto) earthquake, the earliest CSID emerge at 420 (470) seconds from the earthquake onset and attain the peak at 470 (540) seconds. The CSID detection time from observation and simulation are in fairly good agreement. The fast simulation capability positions the method as a significant advance for NRT applications. The rapid CSID detection and simulations could potentially be used for early tsunami warnings at distant locations. At near-field locations, their formations from simulation can be forecasted before their detection by GNSS. Since key input parameters, including HR-GNSS data, seismic data, the background atmospheric and ionospheric conditions, and positions of GNSS line-of-sight, can be quickly obtained after an earthquake, our newly developed method can, indeed, be used in NRT.

## Data Availability Statement

The GNSS data for the 2011 Sanriku and 2024 Noto-Peninsula earthquakes are available from the GeoSpatial Authority of Japan (GSI, https://terras.gsi.go.jp/) via https://www.gsi.go.jp/ENGLISH/geonet_english.html. The seismic data can be downloaded from the National Research Institute for Earth Science and Disaster Resilience (NIED)-F-net (https://www.fnet.bosai.go.jp/fnet/event/joho.php?LANG=en) administers the seismic data.

## References

Afraimovich, E. L., Perevalova, N. P., Plotnikov, A. V., & Uralov, A. M. (2001). The shock-acoustic waves generated by earthquakes. *Annales Geophysicae*, *19*(4), 395–409. https://doi.org/10.5194/angeo-19-395-2001

Astafyeva, E. (2019). Ionospheric detection of natural hazards. *Reviews of Geophysics*, *57*(4), 1265–1288. https://doi.org/10.1029/2019RG000668

Astafyeva, E., Lognonné, P., & Rolland, L. (2011). First ionospheric images of the seismic fault slip on the example of the Tohoku-Oki earthquake. *Geophysical Research Letters*, *38*(22), L22104. https://doi.org/10.1029/2011GL049623

Astafyeva, E., Rolland, L., Lognonné, P., Khelfi, K., & Yahagi, T. (2013). Parameters of seismic source as deduced from 1HZ ionospheric GPS data: Case study of the 2011 Tohoku-oki event. *Journal of Geophysical Research: Space Physics*, *118*(9), 5942–5950. https://doi.org/10.1002/jgra.50556

Astafyeva, E., Shalimov, S., Olshanskaya, E., & Lognonné, P. (2013). Ionospheric response to earthquakes of different magnitudes: Larger quakes perturb the ionosphere stronger and longer. *Geophysical Research Letters*, *40*(9), 1675–1681. https://doi.org/10.1002/grl.50398

Astafyeva, E., & Shults, K. (2019). Ionospheric GNSS imagery of seismic source: Possibilities, difficulties, and challenges. *Journal of Geophysical Research: Space Physics*, *124*(1), 534–543. https://doi.org/10.1029/2018JA026107

Bagiya, M. S., Heki, K., & Gahalaut, V. K. (2023). Anisotropy of the near-field coseismic ionospheric perturbation amplitudes reflecting the source process: The 2023 February Turkey earthquakes. *Geophysical Research Letters*, *50*(14), e2023GL013931. https://doi.org/10.1029/2023GL013931

Bagiya, M. S., Thomas, D., Astafyeva, E., Bletery, Q., Lognonné, P., & Ramesh, D. S. (2020). The ionospheric view of the 2011 Tohoku-Oki earthquake seismic source: The first 60 seconds of the rupture. *Scientific Reports*, *10*(1), 5232. https://doi.org/10.1038/s41598-020-61749-x

Balachandran, N. K. (1968). Acoustic-gravity wave propagation in a temperature- and wind-stratified atmosphere. *Journal of the Atmospheric Sciences*, *25*(5), 818–826. https://doi.org/10.1175/1520-0469(1968)025(0818:AGWPIA)2.0:CO;2

Brissaud, Q., & Astafyeva, E. (2022). Near-real-time detection of co-seismic ionospheric disturbances using machine learning. *Geophysical Journal International*, *230*(3), 2117–2130. https://doi.org/10.1093/gji/gga167

Cahyadi, M. N., & Heki, K. (2015). Coseismic ionospheric disturbance of the large strike-slip earthquakes in North Sumatra in 2012: Mw dependence of the disturbance amplitudes. *Geophysical Journal International*, *200*(1), 116–129. https://doi.org/10.1093/gji/ggu343

Calais, E., & Minster, J. B. (1995). GPS detection of ionospheric perturbations following the January 17, 1994, Northridge earthquake. *Geophysical Research Letters*, *22*(9), 1045–1048. https://doi.org/10.1029/95GL00168

Campos, L. M. B. C. (1983). On three-dimensional acoustic-gravity waves in model non-isothermal atmospheres. *Wave Motion*, *5*(1), 1–14. https://doi.org/10.1016/0165-2125(83)90002-1

Chum, J., Hruska, F., Zednik, J., & Lastovicka, J. (2012). Ionospheric disturbances (infrasound waves) over the Czech Republic excited by the 2011 Tohoku earthquake. *Journal of Geophysical Research*, *117*(A8), A08319. https://doi.org/10.1029/2012JA017767

Fuso, F., & Ravanelli, M. (2024). Probing the ionospheric effects of the 2020 Aegean Sea earthquake: Leveraging GNSS observations for tsunami early warning in the Mediterranean. *Journal of Geophysical Research: Space Physics*, *129*(12), e2024JA032946. https://doi.org/10.1029/2024JA032946

Heki, K. (2024). Decoding multiple source signatures in coseismic ionospheric disturbances of the 2024 January Mw7.5 Noto-Peninsula earthquake, Central Japan. *Earth and Planetary Science Letters*, *640*, 118796. https://doi.org/10.1016/j.epsl.2024.118796

Heki, K., Bagiya, M. S., & Takasaka, Y. (2022). Slow fault slip signatures in coseismic ionospheric disturbances. *Geophysical Research Letters*, *49*(24), e2022GL101064. https://doi.org/10.1029/2022GL101064

Heki, K., & Ping, J. (2005). Directivity and apparent velocity of the coseismic ionospheric disturbances observed with a dense GPS array. *Earth and Planetary Science Letters*, *236*(3), 845–855. https://doi.org/10.1016/j.epsl.2005.06.010

Hofmann-Wellenhof, B., Lichtenegger, H., & Wasle, E. (2008). GNSS - global navigation satellite systems. Springer-Verlag Vienna. https://doi.org/10.1007/978-3-211-73017-1

Hosseini, K., & Sigloch, K. (2017). Obspydmt: A python toolbox for retrieving and processing large seismological data sets. *Solid Earth*, *8*(5), 1047–1070. https://doi.org/10.5194/se-8-1047-2017

## Acknowledgments

This work was supported by the joint French-Brazilian project "IONO-DIET" funded by the French National Research Agency (ANR, Grant ANR-22-CE49-0011) and by the Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) process 2022/03502-6. We thank the Brazilian Ministry of Science, Technology and Innovation and the Brazilian Space Agency. S. A. Sanchez thanks the FAPESP for financial support through the process 2023/07807-9 and Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) process 88887.351785/2019-00. E. A. Kherani thanks the financial support of the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) through Grant 307496/2015-5. E. Astafyeva thanks the financial support of the French Space Agency (CNES), project "RealDetect." E. R. de Paula acknowledges the support of CNPq 202531/2019-0 and INCT GNSS-NavAer Project under grants CNPq 465648/2014-2 and FAPESP 2017/50115-0.
Inchin, P. A., Kaneko, Y., Gabriel, A.-A., Ulrich, T., Martire, L., Komjathy, A., et al. (2024). Three-dimensional numerical modeling of coseismic atmospheric dynamics and ionospheric responses in slant total electron content observations. *Geophysical Research Letters*, 51(21), e2023GL108068. https://doi.org/10.1029/2023GL108068

Inchin, P. A., Snively, J. B., Kaneko, Y., Zettergren, M. D., & Komjathy, A. (2021). Inerring the evolution of a large earthquake from its acoustic impacts on the ionosphere. *AGU Advances*, 2(2), e2020AV000260. https://doi.org/10.1029/2020AV000260

Jin, S., Occhipinti, G., & Jin, R. (2015). Gnss ionospheric seismology: Recent observation evidences and characteristics. *Earth-Science Reviews*, 147, 54–64. https://doi.org/10.1016/j.earscirev.2015.05.003

Kamogawa, M., Orihara, Y., Tsurudome, C., Tomida, Y., Kanaya, T., Ikeda, D., et al. (2016). A possible space-based tsunami early warning system using observations of the tsunami ionospheric hole. *Scientific Reports*, 6(1), 37989. https://doi.org/10.1038/srep37989

Kawamoto, S., Hiyama, Y., Ohta, Y., & Nishimura, T. (2016). First result from the GEONET real-time analysis system (REGARD): The case of the 2016 Kumamoto earthquakes. *Earth Planets and Space*, 68, 1–12. https://doi.org/10.1186/s40623-016-0564-4

Kawamoto, S., Ohta, Y., Hiyama, Y., Todoriki, M., Nishimura, T., Furuya, T., et al. (2017). REGARD: A new GNSS-based real-time finite fault modeling system for GEONET. *Journal of Geophysical Research: Solid Earth*, 122(2), 1324–1349. https://doi.org/10.1002/2016JB013485

Kherani, E. A., Abdu, M. A., Fritts, D. C., & de Paula, E. R. (2011). The acoustic gravity wave induced disturbances in the equatorial ionosphere. In M. A. Abdu & D. Pancheva (Eds.), *Aeronomy of the earth’s atmosphere and ionosphere* (pp. 141–162). Springer Netherlands. https://doi.org/10.1007/978-94-007-0326-1_10

Kherani, E. A., Lognonné, P., Hébert, H., Rolland, L., Astafyeva, E., Occhipinti, G., et al. (2012). Modelling of the total electronic content and magnetic field anomalies generated by the 2011 Tohoku-Oki tsunami and associated acoustic-gravity waves. *Geophysical Journal International*, 191(3), 1049–1066. https://doi.org/10.1111/j.1365-246X.2012.05617.x

Kherani, E. A., Rolland, L., Lognonné, P., Sladen, A., Klausner, V., & de Paula, E. (2016). Traveling ionospheric disturbances propagating ahead of the Tohoku-oki tsunami: A case study. *Geophysical Journal International*, 204(2), 1148–1158. https://doi.org/10.1093/gji/ggv500

LaBrecque, J., Rundle, J., Bawden, G., Surface, E., & Area, I. F. (2019). Global navigation satellite system enhancement for tsunami early warning systems. *Global Assessment Report on Disaster Risk Reduction*.

Maletckii, B., & Astafyeva, E. (2021). Determining spatio-temporal characteristics of coseismic travelling ionospheric disturbances (CTID) in near real-time. *Scientific Reports*, 11(1), 20783. https://doi.org/10.1038/s41598-021-99906-5

Maletckii, B., & Astafyeva, E. (2022). Near-real-time analysis of the ionospheric response to the 15 January 2022 Hunga Tonga-Hunga Ha'apai volcanic eruption. *Journal of Geophysical Research: Space Physics*, 127(10), e2022JA030735. https://doi.org/10.1029/2022JA030735

Maletckii, B., & Astafyeva, E. (2024). Near-real-time identification of the source of ionospheric disturbances. *Journal of Geophysical Research: Space Physics*, 129(11), e2024JA032664. https://doi.org/10.1029/2024JA032664

Manta, F., Occhipinti, G., Feng, L., & Hill, E. M. (2020). Rapid identification of tsunamigenic earthquakes using GNSS ionospheric sounding. *Scientific Reports*, 10(1), 11054. https://doi.org/10.1038/s41598-020-68097-w

Martire, L., Krishnamoorthy, S., Vergados, P., Romans, L. J., Szilágyi, B., Meng, X., et al. (2023). The GUARDIAN system-a GNSS upper atmospheric real-time disaster information and alert network. *GPS Solutions*, 27(1), 32. https://doi.org/10.1007/s10291-022-01365-6

Melgar, D., & Bock, Y. (2015). Kinematic earthquake source inversion and tsunami runup prediction with regional geophysical data. *Journal of Geophysical Research: Solid Earth*, 120(5), 3324–3349. https://doi.org/10.1002/2014JB011832

Melgar, D., & Hayes, P. G. (2019). Characterizing large earthquakes before rupture is complete. *Science Advances*, 5(5), eaav2032. https://doi.org/10.1126/sciadv.aav2032

Meng, X., Komjathy, A., Verkhoglyadova, O. P., Yang, Y.-M., Deng, Y., & Mannucci, A. J. (2015). A new physics-based modeling approach for tsunami-ionosphere coupling. *Geophysical Research Letters*, 42(12), 4736–4744. https://doi.org/10.1002/2015GL064610

Mikesell, T. D., Rolland, L. M., Lee, R. F., Zedek, F., Coïsson, P., & Dessa, J.-X. (2019). IonoSeis: A package to model coseismic ionospheric disturbances. *Atmosphere*, 10(8), 2073–4433. https://doi.org/10.3390/atmos10080443

Nozuka, Y., Inchin, P. A., Kaneko, Y., Sabatini, R., & Snively, J. B. (2024). Earthquake source impacts on the generation and propagation of seismic infrasound to the upper atmosphere. *Geophysical Journal International*, 238(1), 537–556. https://doi.org/10.1093/gji/ggae170

Occhipinti, G., Rolland, L., Lognonné, P., & Watada, S. (2013). From Sumatra 2004 to Tohoku-Oki 2011: The systematic GPS detection of the ionospheric signature induced by tsunamigenic earthquakes. *Journal of Geophysical Research: Space Physics*, 118(6), 3626–3636. https://doi.org/10.1002/jgra.50322

Ohno, K., Ohta, Y., Kawamoto, S., Abe, S., Hino, R., Koshimura, S., et al. (2021). Real-time automatic uncertainty estimation of coseismic single rectangular fault model using GNSS data. *Earth Planets and Space*, 73, 1–18. https://doi.org/10.1186/s40623-021-01425-0

Rakoto, V., Lognonné, P., Rolland, L., & Coïsson, P. (2018). Tsunami wave height estimation from GPS-derived ionospheric data. *Journal of Geophysical Research: Space Physics*, 123(5), 4329–4348. https://doi.org/10.1002/2017JA024654

Ravanelli, M., Constantinou, V., Liu, H., & Bortnik, J. (2024). Exploring AI progress in GNSS remote sensing: A deep learning based framework for real-time detection of earthquake and tsunami induced ionospheric perturbations. *Radio Science*, 59(9), e2024RS008016. https://doi.org/10.1029/2024RS008016

Ravanelli, M., Occhipinti, G., Savastano, G., Komjathy, A., Shume, E. B., & Crespi, M. (2021). GNSS total variometric approach: First demonstration of a tool for real-time tsunami genesis estimation. *Scientific Reports*, 11(1), 3114. https://doi.org/10.1038/s41598-021-82532-6

Rolland, L. M., Lognonné, P., & Munekane, H. (2011). Detection and modeling of Rayleigh wave induced patterns in the ionosphere. *Journal of Geophysical Research*, 116(A5). https://doi.org/10.1029/2010JA016060

Rolland, L. M., Vergnolle, M., Nocquet, J.-M., Sladen, A., Dessa, J.-X., Tavakoli, F., et al. (2013). Discriminating the tectonic and non-tectonic contributions in the ionospheric signature of the 2011 Mw7.1, dip-slip Van earthquake, Eastern Turkey. *Geophysical Research Letters*, 40(11), 2518–2522. https://doi.org/10.1002/grl.50544

Sanchez, S. A., & Kherani, E. A. (2024). A new analytical simulation code of acoustic-gravity waves of seismic origin and rapid co-seismic thermospheric disturbances energetics. *Atmosphere*, 15(5), 592. https://doi.org/10.3390/atmos15050592

Sanchez, S. A., Kherani, E. A., Astafyeva, E., & de Paula, E. R. (2023). Rapid detection of co-seismic ionospheric disturbances associated with the 2015 Illapel, the 2014 Iquique and the 2011 Sanriku-Oki earthquakes. *Journal of Geophysical Research: Space Physics*, 128(9), e2022JA031231. https://doi.org/10.1029/2022JA031231

Sunil, A. S., Bagiya, M. S., Bletery, Q., & Ramesh, D. S. (2021). Association of ionospheric signatures to various tectonic parameters during moderate to large magnitude earthquakes: Case study. *Journal of Geophysical Research: Space Physics*, 126(3), e2020JA028709. https://doi.org/10.1029/2020JA028709

Tahira, M. (1995). Acoustic resonance of the atmospheric at 3.7 Hz. *Journal of the Atmospheric Sciences*, 52(15), 2670–2674. https://doi.org/10.1175/1520-0469(1995)052(2670:AROTAA)2.0.CO;2

Ten, A., Sorokin, A., Shestakov, N., Ohzono, M., & Titkov, N. (2025). Detecting covolcanic ionospheric disturbances using GNSS data and a machine learning algorithm. *Advances in Space Research*, 75(1), 1052–1065. https://doi.org/10.1016/j.asr.2024.10.030
Thomas, D., Bagiya, M. S., Sunil, P. S., Rolland, L., Sunil, A. S., Mikesell, T. D., et al. (2018). Revelation of early detection of co-seismic ionospheric perturbations in GPS-TEC from realistic modelling approach: Case study. *Scientific Reports*, 8(1), 12105. https://doi.org/10.1038/s41598-018-30476-9

Zedek, F., Rolland, L. M., Mikesell, T. D., Sladen, A., Delouis, B., Twardzik, C., & Coïsson, P. (2021). Locating surface deformation induced by earthquakes using GPS, GLONASS and Galileo ionospheric sounding from a single station. *Advances in Space Research*, 68(8), 3403–3416. https://doi.org/10.1016/j.asr.2021.06.011