A New Analytical Simulation Code of Acoustic-Gravity Waves of Seismic Origin and Rapid Co-Seismic Thermospheric Disturbance Energetics

Abstract

A recent study the detection of coseismic ionospheric disturbances or ionoquakes less than 400 s from the earthquake’s onset. The study also associates these rapid ionoquakes with the seismo-atmosphere–ionosphere (SAI) coupling mechanism energized by acoustic-gravity waves (AGWs) and the subsequent formation of coseismic thermospheric disturbances (CSTDs). The present study outlines a new analytical simulation code for AGWs that resolves the governing equations in the time–altitude and wavenumber domain and confirms the rapid arrival of AGWs in the thermosphere (earlier than the estimated arrival time from the ray-tracing simulation). The rapid arrivals of AGWs are associated with long wavelengths that connect to thermospheric altitudes and propagate with thermospheric sound speeds, avoiding averaging effects from the lower atmosphere. The fast simulation traces the rapid arrival of AGWs in the thermosphere and produces rapid CSTDs within 250–300 s from the earthquake’s onset. The simulation time is much shorter than the formation time of near-field CSTDs, a scenario favorable for the forecasting of CSTDs before observations of ionoquakes. In essence, the fast simulation offers an alternative tool for tracking the evolution of CSTDs.

1. Introduction

Numerous studies report the detection of coseismic ionospheric disturbances or ionoquakes [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15]. Their origin involves seismo-atmosphere–ionosphere (SAI) coupling, acoustic gravity waves (AGWs) [16,17,18] and co-seismic thermospheric disturbance (CSTD) energetics [19,20,21]. When studying the connection between SAI coupling, it is clear that AGWs significantly amplify seismic vibrations by three to four orders of magnitude at ionospheric heights. This amplification is attributed to the exponential decrease in density with height, as compared to the vertical ground velocity. The underlying reason for this AGW amplification lies in the imperative need to conserve kinetic energy. Consequently, even a minor elevation of the ground (typically in the order of millimeters per second) can induce vertical motions in the atmosphere and ionosphere, reaching velocities ranging from several tens to a few hundred meters per second at ionospheric heights [3,4,6]. It is generally assumed that ionoquakes are concentrated around the altitude of the maximum ionosphere ionization. However, recent studies by Thomas et al. (2018) [22], Astafyeva and Shults (2019) [23], and Sanchez et al. (2023) [24] reported the detection of ionoquakes at altitudes of 150–190 km, in which results are below the maximum ionization altitude of the surrounding ionosphere. Moreover, a substantial number of research findings consistently illustrate that the initial SAI coupling in proximity to the earthquake epicenter typically manifests within a timeframe of 8 min following the initiation of the seismic event [1,2,7,8,9,10,11,12,13,25].

On the other hand, a few recent studies on the 7.3-magnitude Sanriku-Oki earthquake of 9 March 2011 (often referred to as the Tohoku foreshock event) demonstrate the detection of ionoquakes in 420–480 s, i.e., earlier than the expected detection time at and above 150 km altitude from the acoustic ray-tracing estimation [22,23]. A recent study by Sanchez et al. (2023) [24] reports rapid ionoquake detection within 400 s from the onset time of the Illapel earthquake. The study also simulates SAI coupling mechanism energized by seismically triggered AGWs and interprets the rapid ionoquakes as being due to the rapid arrival of the AGWs in the thermosphere. Interestingly, Liu et al. (2011) [26] reported rapid ionoquake detection within 420 s of the Tohoku-Oki earthquake of 11 March 2011. Moreover, a numerical simulation study of AGWs, CSTDs, and ionoquakes by Kherani et al. (2012) [27] demonstrated the formation of rapid ionoquakes within 360 s of the onset time of the Tohoku-Oki earthquake.

However, the numerical simulation run time is much slower than the detection time of rapid ionoquakes, which limits their employment for the near-real-time monitoring of rapid ionoquakes. The present study aims to present a fast simulation code of AGWs, based on an analytical solution of the wave equation of AGWs. A previous study by Kherani et al. (2021) [28] presented the derivation of spatial and temporal governing equations of wave amplitudes. This work presents analytical solutions to these equations and simulates the AGWs energetics and CSTDs.

2. Analytical Simulation Code of AGWs

Appendix A.1 presents the governing equations of the SAI coupling mechanism. In a recent study, Kherani et al. (2021) [28] considered the following plane–wave solution for the wave using Equation (A1):

$$u_y(x,y,t)=u_{\mathrm{yt}}\left(t\right)u_{\mathrm{ys}}\left(y\right)e^{(i{k}_{x}x+i{k}_{y}y)}$$

$$u_x(x,y,t)=u_{\mathrm{xt}}\left(t\right)u_{\mathrm{xs}}\left(y\right)e^{(i{k}_{x}x+i{k}_{y}y)}$$

(1)

where $(x,y)$ represents epicentral distance and altitude in the Cartesian coordinate system, t represents the time, $(u_x,u_y)$ represents the horizontal and vertical amplitudes of the AGWs, $(u_{\mathrm{xs}},u_{\mathrm{ys}})$ represent the horizontal and vertical amplitudes of the AGWs that depend only on altitude, $(u_{\mathrm{xt}},u_{\mathrm{yt}})$ represent the horizontal and vertical amplitudes of the AGWs that depend only on time, and $(k_x,k_y)$ correspond to the wave parameters. Equations (A8)–(A10) of Kherani et al. (2021) [28] are a set of spatial governing equations for $(u_{\mathrm{xs}},u_{\mathrm{ys}})$ and time governing equations for $(u_{\mathrm{xt}},u_{\mathrm{yt}})$ of the following form:

$$\frac{{\mathrm{du}}_{\mathrm{ys}}}{\mathrm{dy}}=-k_0{u}_{\mathrm{ys}}$$

$$\frac{{\mathrm{du}}_{\mathrm{xs}}}{\mathrm{dy}}=-(\gamma -1)k_0{u}_{\mathrm{xs}}-(\gamma -2)k_0\frac{k_y}{k_x}\frac{u_{\mathrm{yt}}}{u_{\mathrm{xt}}}{u}_{\mathrm{ys}}$$

(2)

$$\frac{d^2{u}_{\mathrm{yt}}}{d{t}^2}=-(k_y^2{c}^2+{\Omega }_b^2)u_{\mathrm{yt}}-k_x{k}_y{c}^2\frac{u_{\mathrm{xs}}}{u_{\mathrm{ys}}}{u}_{\mathrm{xt}}$$

$$\frac{d^2{u}_{\mathrm{xt}}}{{dt}^2}=-\left(k_x^2{c}^2\right)u_{\mathrm{xt}}-k_x{k}_y{c}^2\frac{u_{\mathrm{ys}}}{u_{\mathrm{xs}}}{u}_{\mathrm{yt}}$$

(3)

where

$$k_0=\frac{\zeta }{c^2},\zeta =\frac{1}{\rho }\frac{dp}{dy},c^2=\frac{\gamma p}{\rho },\mu =\int k_{0}dy$$

$${\Omega }_b^2=\left[(\gamma -1)k_0^2-\frac{k_o}{c^2}\frac{d{c}^2}{dy}\right]c^2$$

Here, ${\Omega }_b$ is the non-isothermal non-hydrostatic Brunt–Vaisala frequency (Equation 6.7a of [29]), $\gamma$ is the ratio of the specific heats, and c is the sound speed.

We note that the spatial governing Equation (2) represents a set of first-order homogeneous differential equations with known analytical solutions, ((A3) and (A4)), as derived in Appendix A.2. The time governing Equation (3) represents coupled oscillators with known analytical solutions (A8), as derived in Appendix A.3. The study derives the general solution with viscosity in Appendix A.4 in the following form:

$$u_y=\sum u_{\mathrm{yt}}\left(t_0\right)u_{\mathrm{ys}}\left(y_o\right)exp\left(\left[\pm i\omega -\frac{1}{2}{\nu }_y^{\mathrm{col}}\right]t+i{k}_{y}y+i{k}_{x}x-\mu \right)$$

$$u_x=\sum \frac{({\omega }^2-k_y^2{c}^2-{\Omega }_b^2)}{k_x{k}_y{c}^2}{u}_{\mathrm{yt}}\left(t_0\right)u_{\mathrm{ys}}\left(y_o\right)exp\left(\left[\pm i\omega -\frac{1}{2}{\nu }_x^{\mathrm{col}}\right]t+i{k}_{y}y+i{k}_{x}x-\mu \right)$$

(4)

Here, $\omega$ is the solution of the following dispersion relation of AGWs, derived as (A5) in Appendix A.3:

$${\omega }^4-(k_y^2{c}^2+k_x^2{c}^2+{\Omega }_b^2){\omega }^2+k_x^2{c}^2{\Omega }_b^2=0$$

(5)

and ${\nu }^{col}$ relates to the viscosity dynamics of the following form:

$${\nu }_y^{\mathrm{col}}=\nu \left(k_y^2+k_x^2+2i{k}_y{k}_0-k_0^2+\frac{d{k}_0}{dy}\right)$$

and

$${\nu }_x^{\mathrm{col}}=\nu (k_y^2+k_x^2-{(\gamma -1)}^2{k}_o^2+2i{k}_y(\gamma -1)k_o)$$

The factor $e^{-\mu }$ in (4) increases with altitude due to the negative values of $(\zeta ,\mu )$. Therefore, the amplitudes $(u_y,u_x)$ of AGWs increase exponentially with altitude, a known characteristic of AGWs. With ${\Omega }_b=0$, the dispersion relation (5) leads to the acoustic wave modes $\omega =\pm {(k_x^2+k_y^2)}^{1/2}c$. From (4), this leads to the following relation between $u_x$ and $u_y$ for the non-dissipative $(\nu =0)$ AGWs:

$$u_x=\frac{\left({\omega }^2-k_y^2{c}^2\right)}{k_x{k}_y{c}^2}{u}_y\stackrel{\omega =kc}{\to }{u}_x=\frac{k_x}{k_y}{u}_y\Rightarrow \nabla \times u=0$$

Since the gravity waves are associated with shear or rotational dynamics, the irrotational velocity field implies the absence of gravity waves and only the excitation of acoustic waves, a scenario that is self-consistent with the condition ${\Omega }_b=0$. It also demonstrates the self-consistent nature of the analytical solution obtained in the present study.

Substitution of $k_y\to k_y+i(\gamma /2)k_0$ in (5) leads to the following dispersion relation:

$${\omega }^4-(k_y^2{c}^2+k_x^2{c}^2+{\Omega }_c^2+{\Omega }_b^2){\omega }^2+k_x^2{c}^2{\Omega }_b^2=0,{\Omega }_c=\frac{\gamma }{2}{k}_{0}c$$

where ${\Omega }_c$ is the acoustic cutoff frequency. Since ${\Omega }_b\ll kc,{\Omega }_b<{\Omega }_c$, the above dispersion relation can also be written as follows:

$${\omega }^4-(k_y^2{c}^2+k_x^2{c}^2+{\Omega }_c^2){\omega }^2+k_x^2{c}^2{\Omega }_b^2=0$$

(6)

which is a known dispersion relation of AGWs (Equation 6.3 of Kelley, 2009 [29]).

The present study focuses on the acoustic wave modes to understand the origin of rapid CSTDs. The condition $(k_x\ll k_y)$ relates to the solution (1) that assumes a pure plane wave solution $\exp \left(i{k}_{x}x\right)$ in the x direction, while assuming $u_{\mathrm{ys}}\left(y\right)\exp \left(i{k}_{y}y\right)$ and $u_{\mathrm{xs}}\left(y\right)\exp \left(i{k}_{y}y\right)$ in the y direction. Therefore, the wavenumber components $(k_y;k_x)$ are restricted variables as follows:

$$\frac{\pi }{N\Delta }⩽k_y⩽\frac{\pi }{\Delta }\frac{\pi }{N\Delta }⩽k_x\ll k_y$$

where $(\Delta =10\mathrm{km},\mathrm{N}=30)$ are the grid resolution and number of grids. The simulation domain covers 0–300 km in altitude and x_ep $\pm 30\Delta$ km in the epicentral distance, where $x_{ep}=0$ is the location of the epicenter. The simulation begins a few minutes before the earthquake onset time $t_{ep}$ and spans 30 min, with a time resolution of 15 s. For the selected range of $(k_y;k_x)$, the wave frequency $\omega$ is estimated using the dispersion relation (5). Furthermore, to incorporate the evanescent and resonance dynamics, the imaginary part $\left(k_{im}\right)$ of $k_y$ is estimated using the dispersion relation (6), and a factor $\exp (k_{im}\Delta )$ is multiplied by the solution (4). In the analytical form of the solution (4), the integral in $\mu$ is to be evaluated numerically with the discretization in space. Moreover, other governing Equations (A2) and (A3) are resolved numerically with the numerical code of Kherani et al. (2012) [27].

At the time $t=0$ seconds, the initial values of atmospheric variables $(\rho ;T;p)$ are derived from the NRLMSIS model [30] using the nrlmsise₋2000 python library of pypi. Figure 1A,B demonstrates the ambient atmosphere is sound speed and acoustic cut-off frequency. These profiles were calculated for the atmospheric conditions of 9 March 2011 at 02:45:20 UT, with coordinates 38.4° N, 142.8° E.

The present study considers the mechanical oscillator mechanism in which the ground vibration from seismic waves couples mechanically to the atmosphere without the loss of momentum. The continuity of vertical velocity across the Earth’s surface establishes the coupling, i.e., at the lower boundary $(y=0$ km). The continuity of the vertical velocity $\left({\mathrm{U}}_{\mathrm{SISM}}\right)$ of the ground vibration associated with the earthquake acts as the force for the generation of the AGWs, i.e.:

$$u_y(x,y=0,t)={\mathrm{U}}_{\mathrm{SISM}}(x,t)$$

$${\mathrm{U}}_{\mathrm{SISM}}(x,t)=A_{\mathrm{ep}}\left(t\right)e^{-{(x-x_{\mathrm{ep}})}^2/{\sigma }_{\mathrm{ep}}^2}$$

where $\left(A_{\mathrm{ep}}\right)$ is the ground velocity vibration (referred to as ground uplift) derived from the seismometer at or closest to the epicenter, and $(x_{\mathrm{ep}},{\sigma }_{\mathrm{ep}})$ are the epicenter and fault size of the earthquake. Figure 1C demonstrates the ${\mathrm{U}}_{\mathrm{SISM}}$ with a data sampling rate of 15 s for a typical earthquake (here the Sanriku-Oki earthquake of 9 March 2011) of magnitude 7.3.

3. Results

Figure 2A,B demonstrate the AGWs’ propagation above the epicenter. We note that from the ground uplift, numerous waves with wavefronts of different slopes or phase speeds propagate into the atmosphere. Waves with significant amplitudes of about 10–20 m/s arrive at thermospheric altitudes of about 150 km about 250–300 s from the earthquake’s onset.

Figure 3A,B demonstrate the co-seismic temperature and pressure disturbance energetics above the epicenter. We note that these co-seismic atmospheric disturbances attain a significant amplitude above 100 km altitude, i.e., co-seismic thermospheric disturbances (CSTDs) acquire large amplitudes among atmospheric disturbances. The CSTDs gain significant amplitudes in 250–300 s from the earthquake’s onset time.

Figure 2 and Figure 3 reveal that the large-amplitude AGWs and co-seismic atmospheric disturbances concentrate around the travel time $(y/c)$ of sound rays, denoted by the green circle. This suggests that the group velocities of large-amplitude waves at a given altitude are close to the sound speed at that altitude (we discuss the wave energetics in detail in the following section).

4. Discussion

The simulation results in Figure 2 and Figure 3 reveal the robust energetics of fast-propagating AGWs and their role in rapid CSTD formation. The energetics include the amplification of tiny ground vibrations by a factor of ${10}^3$ due to the amplification factor $e^{-\mu }$ in Equation (4) and the rapid propagation of AGWs with an arrival time of 250–300 s in the thermosphere.

Figure 4 demonstrates the wavefrequency $\omega$ estimation from the dispersion relation (5). A wide range of frequencies in the acoustic range are excited for vertical wavelengths $\left({\lambda }_y\right)$ ranging between 40 km to 140 km. Some of the longer waves (longer than about 70 km) have frequencies shorter than the ${\Omega }_c$ around tropopause and mesopause, and they are responsible for the resonance modes [31]. It is evident that the ground uplift triggers numerous waves with altitude-varying wave frequency and phase speeds $(\omega {\lambda }_y/2\pi )$ into the atmosphere. This is also evident from the presence of numerous slopes of wavefronts in Figure 2.

In order to understand the origin of altitude-varying phase speeds and the origin of the rapid formation of CSTDs, a simulation is carried out with average atmospheric sound speed (=422 m/s), as considered in the ray-tracing analysis. The average is obtained from ground to 200 km altitude. This particular altitude is chosen since strong CSTDs are below 200 km altitude, as evident in Figure 2 and Figure 3. Secondly, the sound speed varies rapidly with altitude below this altitude. Figure 5 presents the results from the ray-tracing simulation. We note that the propagation time of strong wavefronts follows the travel time $(y/c)$ of sound rays, i.e., all wavefronts propagate with phase speeds close to the average sound speed (drawn as green circles). Their arrival time at 200 km altitude is about 480 s from the earthquake’s onset. According to Heki and Ping (2005) [32], sound rays arrive from 480 s onward at 200 km altitude (Figure 3d of Astafyeva and Shults, 2019 [23]). The consistent ray-tracing simulation in the present study validates the AGW energetics in altitude-varying sound profiles, as presented in Figure 2 and Figure 3.

Figure 6 compares the arrival time of strong wavefronts at atmospheric heights for the original simulation results of Figure 2A and the ray-tracing simulation. We note that the original simulation produces co-seismic atmospheric disturbances earlier than the ray-tracing simulation. Moreover, CSTDs formation time ranges between 250 and 400 s in the original simulation, much earlier than the corresponding time in the ray-tracing simulation. In Figure 2, the early arriving waves have wavelengths comparable to the size of the longest atmospheric duct of about 100 km. For frequencies above the upper cutoff frequency of 4.3 mHz, the phase velocity of these waves can be larger than 430 m/s, i.e., the phase velocity is in the range of thermospheric sound speed. These waves can arrive at 200 km altitude as early as 300 s from the earthquake’s onset. Therefore, in the rapid formation of CSTDs, the high-frequency–long-wavelength AGWs participate, as also found in a simulation study of the Tohoku-Oki earthquake [27].

It is evident from Figure 6 that the AGWs follow the travel time $(y/c_s)$ trajectory of the sound ray that highlights the following propagation characteristics. Owing to the various scale heights presented in the atmosphere, numerous wavelengths are launched simultaneously from the ground at a given frequency, as evident from Figure 4. Each wavelength connects an altitude equal to the corresponding wavelength and arrives at that altitude with the sound speed of that altitude. In other words, at a given frequency, a longer wavelength arrives at a higher altitude with a faster travel time, while a shorter wavelength arrives at a lower altitude with a slower travel time. Figure 6 clearly shows this aspect. A monotonic travel time in Figure 6 corresponds to the average sound speed assumption in the ray-tracing simulation. However, in reality, an altitude-varying sound speed permits non-monotonic travel time.

In Figure 4A, the wave frequencies corresponding to the wavelengths longer than 70 km become shorter than ${\Omega }_c$ around tropopause and mesopause. Such waves become evanescent and subsequently excite the resonance modes in the atmosphere [31]. Figure 7A demonstrates the spectrogram of integrated thermospheric pressure disturbance. We note the formation of the strongest CSTDs in the 3.5–5 mHz frequency range, corresponding to the resonance modes. The present study runs another simulation that excludes wavelengths longer than 70 km, i.e., it excludes the resonance energetics. Figure 7B demonstrates the resulting spectrogram of the integrated CSTDs. We note the absence of significant amplification in the 3.5–5 mHz range, which suggests the absence of resonance energetics. Therefore, the AGWs with wavelengths longer than 70 km are responsible for the resonance energetics and strong CSTDs. Previous studies found that wavelengths comparable to or longer than the size of mesospheric ducts contribute to resonance energetics. The relative altitudes of mesopause and tropopause are the lower limits on the size of the mesospheric duct [33], which is about 60–70 km in Figure 1B. The work by Tahira (1995) [31] found resonance modes with a node at the ground and a loop at 67.4 km altitude, above which the waves at 3.7 mHz become evanescent, i.e., the lower wavelength limit of resonance modes is about 70 km, as also found in the present study. Therefore, the results of the present study are in agreement with those of previous studies.

Numerous recent studies have reported the presence of resonance modes in co-seismic ionospheric disturbances [34]. Moreover, normal mode simulation studies explain this disturbance as a product of atmosphere–ionosphere coupling energized by AGWs and CSTDs [34]. According to Astayeva and Shults (2019) [23], rapid co-seismic ionospheric disturbance detection offers a complementary scenario for achieving early warnings of tsunamis. Moreover, if their detection mode is made possible during the mainshock, it should be one of the most promising products related to near-real-time (NRT) ionospheric seismology. A recent study by Sanchez et al. (2023) [24] reported rapid ionoquakes associated with three strong earthquakes including the Sanriku-Oki earthquake. They have also presented the simulation results of rapid ionoquakes, though they did not elaborate on the fast simulation code and CSTDs energetics. The CSTDs simulation in the present study demonstrates the robust energetics of AGWs and the rapid formation of CSTDs that suggest the possibility of rapid formation of ionoquakes. However, the formation time of about 300 s is slower than the typical mainshock duration of about 30–60 s. Therefore, NRT detection of CSTDs and ionoquakes may not be possible from the measurements. However, since the simulation time of the CSTDs is about 30–60 s, the analytical simulation can forecast the CSTDs before their detection from the observations. In other words, fast simulation of rapid CSTD energetics presents numerous possibilities for their forecasting before the detection of their ionospheric counterpart, and may be a potential candidate for NRT monitoring of the earthquakes.

Studies by Thomas et al. (2018) [22], Astafyeva and Shults (2019) [11], and Sanchez et al. (2023) [24] reported the detection of ionoquakes at altitudes of 150–190 km for the Sanriku-Oki earthquake. Figure 2 of the present work finds significant amplitudes of AGWs in 150–200 km altitude regions. For this earthquake event, the ionospheric density peaks at around 280 km altitude [11]. Together with the results from the present study, it can be speculated that significant amplitude of AGWs in 150–200 km altitude regions and a significant amount of electron density at this altitude range offer favorable conditions for large-amplitude ionoquakes. The simulation of ionoquakes is beyond the scope of the present study, and we aim to present it in a separate study.

5. Summary

The study aims to simulate co-seismic atmospheric disturbances, estimate their formation time at thermospheric heights, and compare these estimations with ray-tracing estimations. The simulation examines the seismo-atmosphere coupling mechanism energized by the AGWs energetics for the formation of co-seismic thermospheric disturbances (CSTDs) from seismic vibrations on the ground. The study presents an analytical framework for AGWs’ simulation that resolves the governing equations in the time–space and wavenumber domain.

The simulation reveals robust energetics of AGWs, mainly in their rapid arrival at the thermospheric heights, as well as the rapid formation of CSTDs. The simulation produces waves with altitude-varying frequencies and phase speeds equal to the sound speed so that at thermsopheric height, their propagation is determined by the local sound speed instead of the average sound speed of 422 m/s. The waves with wavelengths longer than 70 km produce strong CSTDs within about 300 s of the earthquake’s onset, which is earlier than the ray-tracing estimation of about 480 s from the average ray speed of 422 m/s. The study also highlights the importance of wavelengths longer than 70 km for exciting resonance modes in CSTDs.

CSTDs are crucial components of the SAI coupling mechanism and the formation of co-seismic ionospheric disturbances or ionoquakes. A few recent studies have reported the detection of rapid ionoquakes in less than 400 s, which are produced by the fast arrival of AGWs at the ionospheric heights. The present study details the energetics of the fast arrival of AGWs and the formation of rapid CSTDs from a new analytical simulation code. The simulation produces these rapid CSTDs with a simulation time of 30–60 s, much faster than their time of formation. Therefore, the fast simulation code of AGWs can be employed to simulate rapid ionoquakes, the most promising observable aspect of near-real-time ionospheric seismology.