Upper-atmosphere responses to the 2022 Hunga Tonga–Hunga Ha′apai volcanic eruption via acoustic gravity waves and air–sea interaction

A multi-group of strong atmospheric waves (wave packet nos. 1–5) over China associated with the 2022 Hunga Tonga–Hunga Ha′apai (HTHH) volcano eruptions were observed in the mesopause region using a ground-based airglow imager network. The horizontal phase speed of wave packet nos. 1 and 2 is approximately 309 and 236 m s⁻¹, respectively, which is consistent with Lamb wave L0 mode and L1 mode from theoretical predictions. The amplitude of the Lamb wave L1 mode is larger than that of the L0 mode. The wave fronts of Lamb wave L0 and L1 below the lower thermosphere are vertical, while the wave fronts of L0 mode tilt forward above the lower atmosphere, exhibiting internal wave characteristics which show good agreement with the theoretical results. Two types of tsunamis were simulated; one type of tsunami is induced by the atmospheric-pressure wave (TIAPW), and the other type of tsunami is directly induced by the Tonga volcano eruption (TITVE). From backward ray-tracing analysis, the TIAPW and TITVE were likely the sources of wave packet nos. 3 and 4–5, respectively. The scale of tsunamis near the coast is very consistent with the atmospheric AGWs observed by the airglow network. The atmospheric gravity waves (AGWs) triggered by TITVE propagate nearly 3000 km inland with the support of a duct. The atmospheric-pressure wave can directly affect the upper atmosphere and can also be coupled with the upper atmosphere through the indirect way of generating a tsunami and, subsequently, tsunami-generating AGWs, which will provide a new understanding of the coupling between ocean and atmosphere.

1 Introduction

Hunga Tonga–Hunga Ha′apai (HTHH) volcano, which erupted at 04:14:45 UT on 15 January 2022, produced the largest volcanic eruption in terms of the energy release of a single event since the Krakatoa volcanic eruption (Symons, 1888) in 1883. This volcanic eruption triggered broad-spectrum atmospheric disturbances (Adam, 2022; Duncombe, 2022; Wright et al., 2022), including Lamb waves (Zhang et al., 2022), acoustic waves, gravity waves (GWs) (Liu et al., 2022), and shock waves (Astafyeva et al., 2022). In addition, the traveling ionospheric disturbances (TIDs) caused by this volcanic eruption have also been reported (Themens et al., 2022; Lin et al., 2022).

Lamb waves are external waves propagating along Earth's surface at the speed of sound (Beer, 1974). They are non-dispersive or nearly non-dispersive (Francis, 1973) and can propagate horizontally over long distances. A Lamb wave mainly occupies the troposphere, and its perturbation pressure decays exponentially with height (Yeh and Liu, 1974). The Lamb waves excited by the Tonga volcano eruptions went around the Earth several times (Amores et al., 2022; Duncombe, 2022). Sepúlveda et al. (2023) found that the wind field strongly affects the morphology and propagation of a Lamb wave. Liu et al. (2023) reproduced the Lamb wave L0 and L1 modes consistently with theoretical predictions (Francis, 1973) using the high-resolution Whole Atmosphere Community Climate Model with the thermosphere/ionosphere extension (WACCM-X). Li et al. (2023) identified the Lamb wave L1 mode using phase-leveling amplitude technology based on a global navigation satellite system (GNSS) total electron content (TEC). Poblet et al. (2023) reported that the strong perturbations in the meteor radar horizontal wind field over South America are caused by the Lamb wave L1 mode associated with the 2022 HTHH volcano eruption.

Acoustic gravity waves (AGWs) are mechanical waves in compressible fluids in a gravity field (Gossard and Hooke, 1975). If the frequencies are much larger than the buoyancy frequency, then AGWs tend towards an acoustic wave mode; when the frequency is much smaller than the buoyancy frequency, then the fluid can be considered incompressible, and the AGWs tend towards internal GWs mode. The term “acoustic gravity waves” is usually used when restoring forces due to both gravity and compressibility is important. AGWs are known to play a significant role in the coupling between the atmosphere or ionosphere and the ocean (Press and Harkrider, 1962; Harkrider and Press, 1967; Donn and Balachandran, 1981; Azeem et al., 2017). Atmospheric-pressure waves are mechanical waves that are related to the density of the atmosphere. Compression and expansion are the high-pressure and low-pressure regions of motion in a medium.

The 2022 HTHH volcano eruption triggered tsunamis that affected the whole world (Carvajal et al., 2022; Ghent and Crowell, 2022). Conventional tsunamis are typically generated by localized sea surface displacements caused by sources such as earthquakes and volcanoes, similar to the tsunamis directly induced by the 2022 Tonga volcano eruption (TITVE). Another type of tsunami is induced by an atmospheric-pressure wave (TIAPW) (Kubota et al., 2022; Gusman et al., 2022). Tsunamis can generate upward-propagating AGWs through the water–air interface and propagate to the thermosphere/ionosphere (Hines, 1972; Peltier and Hines, 1976; Hickey et al., 2009, 2010; Occhipinti et al., 2013; Vadas et al., 2015; Laughman et al., 2017; Nishikawa et al., 2022; Pradipta et al., 2023). Using the red line airglow imager, Makela et al. (2011) detected an airglow disturbance in Hawaii that arrived 1 h earlier than the tsunami generated by the 11 March 2011 Tohoku earthquake. Also using the red line airglow, Smith et al. (2015) observed a tsunami and GWs almost simultaneously in Chile. Inchin et al. (2020) used a three-dimensional (3-D) numerical model to simulate the atmospheric AGWs generated by a tsunami. They found that bathymetry variations significantly affected the tsunamis and the AGWs excited by tsunamis, leading to their nonlinear evolution process. More recently, Inchin et al. (2022) performed the numerical simulations of mesopause airglow radiation fluctuations induced by tsunami-generated AGWs and found that large-scale tsunamis can cause detectable and quantitative disturbances of mesopause airglow through AGWs.

As far as we know, the research on the impact of tsunami-induced atmospheric AGWs on the atmosphere and ionosphere shown above is all related to conventional tsunamis. There are only two rare studies on the ground-based airglow observations of AGWs caused by this conventional tsunami, and both are limited to red line observations (Makela et al., 2011; Smith et al., 2015). However, the observation of tsunami-induced AGWs in the mesopause region observed by ground-based airglow imaging has never been reported. In this study, we first reported the propagation characteristics of the AGWs generated by the tsunamis triggered by the 2022 HTHH volcano eruptions in the mesopause region using the ground-based airglow imager observation network. We then focus on the coupling process of atmospheric-pressure waves triggering tsunamis, and then we look at tsunamis generating atmospheric AGWs through the air–water–air coupling process in the far-field area of the 2022 HTHH volcano eruption.

2 Data and methods

2.1 Multi-layer airglow imager network

A multi-layer airglow observation network (Xu et al., 2021) was built to study atmospheric disturbances excited by severe weather events, such as thunderstorms (Xu et al., 2015), typhoons (Li et al., 2022), and volcanic activities. Figure 1 shows the distribution of the multi-layer airglow observation network station. The multi-layer airglow observation network mainly includes a hydroxyl radical (OH) airglow network, which has been used to observe the airglow layer at the height of 87 km; an atomic oxygen emission (OI) airglow network has been used to observe the airglow layer at the height of 250 km. In addition, there were 557 nm airglow and Na airglow imagers installed at some stations, such as Xinglong station (40.4° N, 117.6° E) and Lhasa (29.7° N, 91.0° E). The airglow network can provide observations with a high temporal and spatial resolution. The time resolution of an OH airglow imager is 1 min, while the resolution of the OI 557 nm and OI 630 nm airglow imager is 3 min, respectively. The spatial resolution of the airglow imager at the airglow layer is not uniform. The resolutions of OH, OI 557 nm, and OI 630 nm airglow in the zenith direction are 0.27, 0.29, and 0.77 km, respectively, while in the zenith angle of 60° the resolutions are 1.01 (OH), 1.11 (OI 557 nm), and 2.65 km (OI 630 nm), respectively.

Figure 1 The distribution of airglow network stations, along with the large circle centered on the Tonga volcano and its radius length, is marked in the figure.

2.2 Spectral analysis of atmospheric wave parameters

The airglow image was calibrated with the help of standard star map (Garcia et al., 1997) and projected into geospatial space. The background radiation is removed by the time differential (TD) method (Swenson and Mende, 1994) to highlight atmospheric fluctuations. The atmospheric wave parameters (horizontal wavelength λ_h, observed horizontal phase speed c, and the relative intensity perturbation $I^{\prime }/I$) are extracted from the spectral analysis method. Figure 2c presents the two-dimensional cross-spectrum obtained from Fig. 2a and b. Zonal (k_x) and meridional (k_y) wave numbers are determined from the peak position of the spectra. The horizontal wavelengths λ_h are obtained from the expression of ${\mathit{\lambda }}_{\mathrm{h}}=\mathrm{2}\mathit{\pi }/\sqrt{k_x^{\mathrm{2}}+k_y^{\mathrm{2}}}$. The observed speeds c are calculated from the phase (φ) (Fig. 2d) at the maximum peak of the cross-spectrum as $c=\frac{\mathit{\phi }}{\mathrm{2}\mathit{\pi }}\frac{{\mathit{\lambda }}_{\mathrm{h}}}{\mathrm{\Delta }t}$, where Δt is the time interval between the two TD images. The amplitudes of intensity perturbations were calculated by integrating the power surrounding the central peaks of the power spectrum. To eliminate noise, the energy of the wave spectrum should be greater than 10 % of the total spectrum (Tang et al., 2005).

Figure 2 The time difference images in panels (a) and (b) have been obtained from the Xinglong OH airglow imager on the night of 15 February 2022. Each image is projected on an area of 900 km×900 km. The (c) cross-spectrum and (d) phase obtained from the yellow box area in panels (a) and (b) using 2-D fast Fourier transform.

2.3 Tsunami simulation model

Tonga submarine volcano erupted on 15 January 2022 and generated tsunamis that were detected around the globe, particularly affecting the Pacific region. In this study, two types of tsunamis were simulated, namely conventional tsunami simulations and atmospheric-pressure, wave-induced tsunami simulations. The linear shallow-water equations in the spherical coordinate system are used to simulate the tsunamis from the localized source and atmospheric-pressure wave. The continuity equation of a linear shallow-water wave model with spherical coordinates is as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\eta }}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{R\sin \mathit{\theta }}}\left[{\displaystyle \frac{\partial \left(ud\right)}{\partial \mathit{\phi }}}+\sin \mathit{\theta }{\displaystyle \frac{\partial \left(vd\right)}{\partial \mathit{\theta }}}\right]=\mathrm{0},\end{array}$$

where η is the free-surface elevation (m), d is the water depth (m), R is the Earth's radius (6 371 000 m), φ is longitude, and θ is colatitude.

While the momentum equations of the linear shallow-water wave model are

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial u}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{R\sin \mathit{\theta }}}\left[g{\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathit{\phi }}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial p}{\partial \mathit{\phi }}}\right]+fv=\mathrm{0},\text{(3)}& {\displaystyle \frac{\partial v}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{R}}\left[g{\displaystyle \frac{\partial \mathit{\eta }}{\partial \mathit{\theta }}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial p}{\partial \mathit{\theta }}}\right]-fu=\mathrm{0},\end{array}$$

where u is the velocity along the lines of longitude (m s⁻¹), v is the velocity along the lines of latitude, g is the gravitational acceleration (9.81 m s⁻²), p is the atmospheric pressure (Pa), ρ is the seawater density (1026 kg m⁻³), and f is the Coriolis coefficient. For the atmospheric-pressure wave-induced tsunami simulation, the moving pressure term is used as an input to the tsunami simulation momentum equation. The atmospheric-pressure wave model is based on the Eq. (1) in Gusman et al. (2022).

For the tsunami simulations from a localized source, a B-spline function (Koketsu and Higashi, 1992) is used below to represent the circular water uplift source at the volcano.

$$\begin{array}{}\text{(4)}& {\displaystyle }\begin{array}{rl}{\displaystyle }f(x,y)={\sum }_{i=\mathrm{0}}^{\mathrm{3}}{\sum }_{j=\mathrm{0}}^{\mathrm{3}}& {\displaystyle }{c}_{k+i,l+j}{B}_{\mathrm{4}-i}\left({\displaystyle \frac{x-x_k}{h}}\right)\\ {\displaystyle }& {\displaystyle }\times B_{\mathrm{4}-j}\left({\displaystyle \frac{y-y{}_l}{h}}\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(5)}& B_i\left(r\right)=\left\{\begin{array}{ll}{r}^{\mathrm{3}}/\mathrm{6},& i=\mathrm{1}\\ (-\mathrm{3}{r}^{\mathrm{3}}+\mathrm{3}{r}^{\mathrm{2}}+\mathrm{3}r+\mathrm{1})/\mathrm{6},& i=\mathrm{2}\\ (\mathrm{3}{r}^{\mathrm{3}}-\mathrm{6}{r}^{\mathrm{2}}+\mathrm{4})/\mathrm{6},& i=\mathrm{3}\\ (-r^{\mathrm{3}}+\mathrm{3}{r}^{\mathrm{2}}-\mathrm{3}r+\mathrm{1})/\mathrm{6},& i=\mathrm{4}.\end{array}\right.\end{array}$$

x_k and x_l stand for the coordinates of the knots along the x and y axes, h is the characteristic diameter of water uplift, r is the great-circle distance from the volcano eruption center, and $c_{\mathrm{1},\mathrm{1}}=\mathrm{1}$ and the other $c_{k+i,l+j}=\mathrm{0}$. In this study, the modeling domain covers the Pacific Ocean and some parts of the Indian Ocean and the Caribbean with a grid size of 5 arcmin. For detailed tsunami simulation algorithms, please refer to Gusman et al. (2022).

The models for the 2022 HTHH volcanic eruption used in this study was estimated and validated with observations at offshore DART stations around the Pacific Ocean in a previous study (Figs. 3 and 7 in Gusman et al., 2022).

2.4 Ray-tracing method

The following ray-tracing equations (Lighthill, 1978) describe the propagation path of AGWs.

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}{x}_i}{\mathrm{d}t}}={\displaystyle \frac{\partial \mathit{\omega }}{\partial k_i}}=c_{g_i},\text{(7)}& {\displaystyle \frac{\mathrm{d}{k}_i}{\mathrm{d}t}}=-{\displaystyle \frac{\partial \mathit{\omega }}{\partial x_i}},\end{array}$$

where x_i, k_i, $c_{g_i}$, ($i=\mathrm{1},\mathrm{2},\mathrm{3}$), and ω are the position vector, wavenumber vector, group speed, and intrinsic frequency, respectively.

Using the dispersion relation of acoustic gravity wave (Yeh and Liu, 1974), we can assess the vertical propagation state of AGWs. The dispersion relation is as follows:

$$\begin{array}{}\text{(8)}& m^{\mathrm{2}}={\displaystyle \frac{{\mathit{\omega }}^{\mathrm{2}}}{c_{\mathrm{s}}^{\mathrm{2}}}}\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\omega }}_{\mathrm{a}}^{\mathrm{2}}}{{\mathit{\omega }}^{\mathrm{2}}}}\right)-k^{\mathrm{2}}\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\omega }}_{\mathrm{b}}^{\mathrm{2}}}{{\mathit{\omega }}^{\mathrm{2}}}}\right),\end{array}$$

where m is the vertical wave number, k is the horizontal wave number, c_s the local speed of sound, $\mathit{\omega }=k(c-u)$ is the intrinsic frequency, and u is the background wind speed in the direction of wave propagation from meteor radar observations and ERA5 (Hersbach et al., 2020). ${\mathit{\omega }}_{\mathrm{a}}^{\mathrm{2}}=\frac{g}{T}\frac{\mathrm{d}T}{\mathrm{d}z}+\frac{\mathit{\gamma }g}{\mathrm{4}H}$ is the acoustic cutoff frequency, ${\mathit{\omega }}_{\mathrm{b}}^{\mathrm{2}}=\frac{g}{T}\frac{\mathrm{d}T}{\mathrm{d}z}+\frac{(\mathit{\gamma }-\mathrm{1})g}{\mathit{\gamma }H}$ is the buoyancy frequency, g is the gravitational acceleration, and T is temperature from the Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) instrument on the Thermosphere, Ionosphere, Mesosphere Energetics and Dynamics (TIMED) satellite. When ω>ω_a or ω<ω_b, m²>0, AGWs can propagate freely, whereas when ${\mathit{\omega }}_{\mathrm{b}}<\mathit{\omega }<{\mathit{\omega }}_{\mathrm{a}}$, m²<0, the wave is evanescent.

Figure 3 Five strong group atmospheric waves associated with the Tonga volcano eruptions were observed in the mesopause region by the ground-based airglow network. Different colored triangles correspond to each wave event sampling point, while red, blue, green, yellow, and cyan correspond to wave packet nos. 1, 2, 3, 4, and 5, respectively. The red time markers in this figure and the following figure represent the lapse time since the volcano eruption.

3 Results and discussion

3.1 Upper-atmospheric airglow responses to HTHH volcanic eruption via Lamb waves

Five groups of atmospheric waves (wave packet nos. 1–5) were observed in the mesopause region by the ground-based airglow network. Refer to the Supplement (https://doi.org/10.5446/66190; Li, 2024c) for the detailed wave propagation status. To eliminate random disturbances, we also made videos of 2 d before and after the volcanic eruption (https://doi.org/10.5446/s_1689​​​​​​​; Li, 2024b). From the videos, it can be seen that the OH airglow layer was very calm during this period. Figure 3 shows wave packet no. 1 observed by the airglow imager network (top panels). Wave packet no. 1 entered the view of the airglow network approximately 8 h after the HTHH volcanic eruption (left panel of the top row). A total of 3 h after wave packet no. 1 entered the field of view, wave packet no. 2 was observed by the airglow network. The leading front of wave packet no. 2 has an uninterrupted continuous front, which almost covers the whole Chinese mainland (middle panels). Interestingly, we observed AGWs accompanying wave packet no. 2 (hereafter wave packet no. 3) over the northwestern region of the Yellow Sea (left panel of middle row). Wave packet no. 2 always keeps a stable state in the process of propagation and maintains a regular front when propagating over Lhasa station (29.7° N, 91.0° E). Wave packet no. 4 exhibits strong instability characteristics during propagation. Compared to the continuous leading front of wave packet no. 2, the fronts of wave packet nos. 4 and 5 are separated (bottom panels). We also found that wave packet no. 5 propagates more than 3000 km inland (propagating to the area west of 90° E).

Figure 4 shows the distribution of wave parameters for multi-group of atmospheric waves (wave packet nos. 1–5) from cross-spectral analysis. The phase speed of wave packet no. 1 leading front is approximately 309 m s⁻¹. Wave packet no. 2 displays a slightly slower phase speed, with an average phase speed of 236 m s⁻¹. The horizontal phase speeds of the group with wave packet nos. 3–5 are mainly distributed in the range of 200 to 215 m s⁻¹, which is smaller than that of wave packet nos. 1–2. The horizontal wavelengths of these five grouped wave packets are mainly distributed in 80–105 km, while the observation periods are relatively small and mainly concentrated in 5.7–7.2 min. For amplitude, the average amplitude of the Lamb wave L1 mode (5.4 %) is higher than that of the Lamb wave L0 mode (3.2 %). Wave packet nos. 3, 4, and 5 have relatively small amplitudes, mainly distributed between 0.85 % and 1.25 %.

Figure 4 Distribution of (a) horizontal wave wavelength, (b) phase speed, (c) period, and (d) amplitude parameters for a multi-group of atmospheric waves (wave packet nos. 1–5). The calculation of wave packet parameters comes from the average value of the wave passing through the sampling points in Fig. 3.

Figure 5 (a) OI 630 nm airglow observation at 13:13:18 UT. OH airglow network observations when (b) wave packet no. 1 and (c) wave packet no. 2 pass through the zenith direction of Xinglong station at 13:13:34 and at 16:32:16 UT, respectively. (d, e) Himawari-8 6.2 µm brightness temperature at 13:10:00 UT. (f) The surface time series of surface pressure obtained from Xinglong station. The red line represents the time derivative of the pressure. The sudden change in air pressure at 13:15 UT indicates the arrival time of Lamb wave L0. A small disturbance of air pressure that occurs at 16:33 UT indicates the arrival time of Lamb wave L1. The yellow stars represent the location of Xinglong station.

The HTHH volcano eruption produced Lamb waves that propagated around the globe (Wright et al., 2022), causing sudden changes in surface pressure (Omira et al., 2022; Takahashi et al., 2023). Figure 5f shows the surface air pressure data of Xinglong station (40.4° N, 117.6° E). At 13:15 UT on 15 January 2022, the air pressure dropped sharply from 920 to 917.7 Pa, indicating that a Lamb wave had arrived at the surface of Xinglong station at 13:15 UT. A small disturbance of air pressure occurred at 16:33 UT. Figure 5e and d present the Himawari-8 6.2 µm brightness temperature at 13:10:00 UT (Otsuka, 2022). It can be seen that the leading front of Lamb wave L0 mode happens to pass through the zenith direction of Xinglong station. The time when wave packet no. 1 (Fig. 5b) and wave packet no. 2 (Fig. 5c) reach the zenith direction of Xinglong station from the OH airglow observation is 13:13:34 and 16:32:16 UT, which matches the time for surface pressure disturbances quite well. The phase speed of wave packet no. 1 leading front (∼309 m s⁻¹) is very close to the speed of surface Lamb wave (L0 mode). From the Fig. 5, it can be seen that the phase of the Lamb wave L0 mode is almost vertical from the ground to the stratosphere and then to the mesosphere. Wave packet no. 2, with a slower phase speed (∼236 m s⁻¹), is consistent with the Lamb wave L1 mode in theoretical predictions (Francis, 1973) and simulations from the WACCM-X model (Liu et al., 2023). However, at almost the same time, the wave front observed in the thermosphere (Video Supplement; https://doi.org/10.5446/66280; Li, 2024a) with a slightly faster phase speed of 342 m s⁻¹ is nearly 550 km ahead of the wave front in the mesopause region in the horizontal propagation direction and ahead of time by approximately 30 min (Fig. 5a). This is in good agreement with theoretical and modeling results (Fig. 4 in Lindzen and Blake, 1972; Fig. 2 in Liu et al., 2023), which show that the wave fronts of the Lamb wave below the lower thermosphere are vertical and tilt forward above. As for Lamb wave L1 mode, the ground and mesopause region provide waveguide surfaces, resulting in maximum wave energy between the two layers, while the phase does not change with height (Francis, 1973).

Figure 6 The solid red lines indicate leading wave front of wave packet no. 2. The solid yellow lines mark wave packet no. 3; the lines are clearly not parallel to the wave fronts of wave packet no. 2.

As for why the observed Lamb wave L0 shape in the OH airglow layer is not a strong leading wave with much weaker trailing waves, this may be caused by the following factors. It is seen from model simulations that the wave amplitudes of the L0 and L1 modes are not uniform at the wave front. This non-uniformity becomes more pronounced in the upper atmosphere (e.g., Fig. 2 in Liu et al., 2023), probably as a result of the large variation in the background atmosphere propagation conditions. It is thus possible that over certain regions the trailing waves become comparable to the leading wave. It is also possible for the leading wave to gradually dissipate energy and become invisible during propagation by generating trailing waves. In addition, due to the smaller field of view of the airglow imager compared to satellite observations, some structures may be related to local fine structures, especially in the middle and upper layers, where many internal waves have significant amplitudes, which may be relatively more significant than Lamb waves.

As mentioned above, the amplitude of the Lamb wave L1 mode in the mesopause region is greater than that of the L0 mode, which may be due to the fact that the L1 mode is an internal wave below the mesopause (Liu et al., 2023). For an isothermal atmosphere, the Lamb wave L0 mode amplitude grows with altitude z as $e^{\mathit{\kappa }z/H}$, where H is the scale height, $\mathit{\kappa }=(\mathit{\gamma }-\mathrm{1})/\mathit{\gamma }$, and γ is the ratio of specific heats (∼1.4). However, the amplitude of the internal GWs varies as $e^{z/\mathrm{2}H}$. The amplitude of internal waves increases with height at a rate greater than that of surface modes.

Poblet et al. (2023) reported an observation of the Lamb wave L1 mode in the horizontal wind field of meteor radar, but they do not see the Lamb wave L0 mode and argue that the L0 mode is likely a higher-frequency wave and was averaged out. Stober et al. (2023, 2024) found that the anomalous peak signal in the meteor radar wind field cannot be completely determined to be caused by the Lamb wave generated by the Tonga volcanic eruption. On the one hand, meteor radar observations may have filtered out high-frequency Lamb waves. On the other hand, even if Lamb waves are observed in the upper atmosphere, there is still debate over whether they propagate directly to the upper atmosphere or through multi-step vertical coupling process described by Becker and Vadas (2018), Vadas and Becker (2018), and Vadas et al. (2018, 2023).

Figure 6 shows the time sequence of propagation image of wave packet no. 3. We found that with the propagation of wave packet no. 2, there is an AGW (wave packet no. 3) with a certain angle between its phase plane (solid yellow line) and the phase plane of wave packet no. 2. This implies that the source of wave packet no. 3 is different from that of wave packet no. 2. The horizontal wavelength of wave packet no. 3 near the coast is 84 ± 5 km.

Figure 7 Snapshots of simulated tsunamis induced by the atmospheric-pressure wave (left panels) and tsunamis directly induced by the Tonga volcano eruption (right panels).

3.2 Simulation of tsunami induced by HTHH volcano eruption

The 2022 HTHH volcano eruption triggered global atmospheric-pressure waves. The simulated atmospheric-pressure waves propagate at an approximate constant speed of 317 m s⁻¹, and the amplitude decreases with the distance from the volcano (Gusman et al., 2022). Figure 7 shows snapshots of the TIAPW and TITVE simulation results. The leading TIAPW excited by the pressure disturbances travels at the same speed as the atmospheric-pressure wave and is followed by subsequent sea waves generated earlier in the atmospheric-pressure wave propagation which thereafter travels at the conventional tsunami propagation speed. Under a given pressure gradient, the discharge flux in the deep sea is much greater than that in shallow water. A deep-bathymetric feature such as the Kermadec–Tonga Trench can more effectively generate tsunami waves. The wave train following the leading wave traveling over the trench appears to be larger than those traveling in other directions. The propagation speed of TITVE from the shallow-water (long) wave approximation is $v=\sqrt{g{H}_{\mathrm{0}}}$ (Salmon, 2014), where g is the gravitational acceleration, and H₀ is the ocean depth. For seawater with a general depth of 4 km, the speed of shallow-water wave is about 200 m s⁻¹. Therefore, the TIAPW is significantly faster than the TITVE. The amplitude of TITVE is greater than that of tsunamis generated by atmospheric-pressure waves. The wave train following the leading wave of TITVE exhibits finer structures with scales smaller than that of TIAPW. We found that the TIAPW arrived along the coast of the Chinese mainland about 4–5 h earlier than the TITVE.

3.3 Upper-atmosphere responses to HTHH volcanic eruption via air–sea interaction

Figure 8 shows the simulation results of TIAPW and TITVE near the coast of the Chinese mainland 11 h (15:15 UT) and 15 h (19:15 UT) after the volcanic eruption, respectively. Air pressure waves are not very efficient at directly exciting tsunamis in shallow water due to the weaker air–sea coupling (Gusman et al., 2022; Yamada et al., 2022). The Yellow Sea is quite shallow, so the amplitude of the leading of TIAPW is very small there. The leading wave is followed by subsequent waves with larger amplitudes, which propagate in the same direction as the leading wave but at the conventional tsunami speed (Gusman et al., 2022). We found that the TIAPW and TITVE on the continental shelf have shorter wavelengths compared with those in the deep ocean. When the tsunamis approached the coast of China, three groups of AGWs (wave packet nos. 3, 4, and 5) were observed by the airglow network. The time at which the AGWs entered the view of the airglow network was very close to the time when the Tonga tsunamis reached the coast of the Chinese mainland. Wave packet no. 3 entered the airglow network at 15:30 UT, and wave packet nos. 4–5 entered the airglow network at 19:40 UT. This strongly suggests that wave packets detected by the airglow network are correlated to the tsunamis near the coast. We found that as the tsunamis approached the coast of China, they diffracted between Taiwan and the Philippines and became discontinuous. And wave packet nos. 4 and 5 that we observed were also discontinuous, which further confirms the correlation between wave packet nos. 4–5 and discontinuous tsunamis. We estimate that the average wavelength of TIAPW near the coast of the Yellow Sea is approximately 82 ± 4 km, which is very consistent with the horizontal wavelengths of the atmospheric AGWs observed by the airglow network, as mentioned above (84 ± 5 km), while the average wavelengths of TITVE near the coast of the Yellow Sea and South Sea are 95 ± 5 km and 86 ± 5 km, respectively.

Figure 8 Simulated tsunamis induced by the atmospheric-pressure waves (a) and tsunamis directly induced by the Tonga volcano eruption (b) near the coast of the Chinese mainland. The marked time represents the time after the volcanic eruption.

Figure 9 (a) Ascending and descending SABER/TIMED satellite tracks over the Chinese mainland. The background is a representative ocean depth map. (b) Square of vertical wave number m² profiles, with the solid black line profile derived from the ascending track no. 2 (marked by the red circle), the dotted line profile derived from the ascending track no. 1 north (marked by the red circle), the dashed line profile derived from the ascending track no. 1 south (marked by the red triangle), and the dashed–dotted line profile derived from the average of ascending track no. 1 and descending track no. 1 (marked by the red circle) from the SABER/TIMED measurement locations shown in panel (a). The red line represents the OH 1.6 µm emission intensity obtained by SABER/TIMED.

Figure 9a shows three TIMED satellite tracks with descending track no. 1 along the coast of China, ascending track no. 1 located east of the Korean Peninsula, and ascending track no. 2 located in inland China. Figure 9b shows the square of vertical wave number m² profile (black) derived from the average temperature from the limb viewing of the sounding of the atmosphere using SABER/TIMED measurement locations marked by the red circles and triangles in Fig. 9a. We take the average temperature of ascending track no. 1; descending track no. 1 serves as the background temperature for wave packet no. 3, and ascending track no. 1 serves as the background temperature of wave packet nos. 4–5 when they propagate in the coastal vicinity. We take ascending track no. 2 as the background temperature of wave packet nos. 4–5 when they propagate in inland China. The peak height of the OH airglow layer is 87 km. We found that the propagation of wave packet no. 3 (dashed–dotted line) is in a state of free propagation in the coastal vicinity.

Figure 10 The background field used for the ray-tracing analysis for the TIAPW event. (a) Saber temperature (red) comes from the average temperature of ascending track no. 1 and descending track no. 1 in Fig. 9, and the ERA5 temperature (black) comes from the average of 15:00 and 16:00 UT. (b) Meteor zonal wind field (red) and ERA5 zonal wind field (black). (c) Meteor meridional wind field (red) and ERA5 meridional wind field (black). The two red and black lines in panels (b) and (c) are, respectively, from 15:00 and 16:00 UT. The green lines represent the average of two lines. The meteor radar wind field is from Beijing station.

Figure 11 (a) Backward ray-tracing results of wave packet no. 3 observed by the OH airglow network. The red triangles and red crosses represent the trace start and termination points, respectively. (b) Simulated tsunamis induced by the atmospheric-pressure wave (TIAPW) corresponding to the dotted rectangular area in panel (a). (c) Ray paths of the wave starting from the seven sampling points in panel (a).

Figure 10 shows the background field used for ray-tracing analysis for the TIAPW event. The temperature comes from TIMED/SABER and ERA5, and wind data comes from meteor radar and ERA5. Meteor radar wind field is from Beijing station (40.3° N, 116.2° E). Figure 11 shows the results of ray tracing for wave packet no. 3. We find that the source location of AGWs over the coast of the Chinese mainland falls near-coast where the tsunami occurred.

A tsunami simulation shows that the surface wave height along the coast of the Chinese mainland is of the order of 2 cm. There have been theoretical (Peltier and Hines, 1976) and observational (Grawe and Makela, 2015, 2017) studies on the relationship between the amplitude of tsunamis and GWs. Peltier and Hines (1976) found that a tsunami amplitude of ±1 cm at sea level can cause a vertical motion of the ionospheric E layer and F layer ±100 m. More direct observational evidence is that Grawe and Makela (2017) provided an airglow observation of tsunami-generated ionospheric signatures over Hawaii caused by the 16 September 2015 Illapel earthquake. They found that vertical disturbances on the sea surface not exceeding 2 cm (Fig. 3b in Grave and Makela, 2017) can create detectable signatures in the ionosphere (Fig. 1 in Grave and Makela, 2017). Therefore, we suggest that the waves with larger amplitudes following the leading of TIAPW interact with the atmosphere after arriving at the coast of the Chinese mainland to generate the upward-propagating AGW packet.

Figure 12 Similar to Fig. 10 but for ray-tracing analysis for the TITVE events. The SABER temperature field in panel (a) comes from ascending track no. 1 (21:17:50, 21:18:33, 21:19:43, and 21:20:43 UT) in Fig. 9, and the meteor radar wind fields in panels (b) and (c) come from Beijing station. The SABER temperature field in panel (d) is from ascending track no. 1 (21:12:51, 21:14:01, and 21:14:44 UT) in Fig. 9, and the meteor radar wind fields in panels (e) and (f) are from Ledong station.

According to the theory of AGW dispersion, the AGW propagating obliquely has the following approximate relationship: $\sin \left(\mathit{\phi }\right)\sim T_{\mathrm{B}}/T$; φ is the oblique propagation angle, T_B is the buoyancy period, and T is the intrinsic period. Azeem et al. (2017) found disturbances in the ionosphere excited by the 2011 Tohoku tsunamis when they reached the West Coast of the United States. They concluded that the fluctuations observed in TEC satisfy the AGW dispersion relation, and the period and horizontal wavelength of the TEC disturbances increased with distance from the West Coast of the USA.

From the airglow network observations, we found that wave packet nos. 4–5 excited by the tsunamis continue to propagate over the mainland more than 3000 km from the coast. If the AGWs observed by the airglow network propagate freely rather than being constrained by a duct, we will obtain the propagation characteristics similar to that observed by Azeem et al. (2017) in the ionosphere from TEC observations. T_B is about 5 min from the SABER/TIMED observation. The period of wave packet no. 3 is between 5.5 and 8.5 min. The minimum propagation angle φ​​​​​​​ equals 35°, and the corresponding maximum propagation distance L is 125 km from $L\sim H_{\text{oh}}/\tan \left(\mathit{\phi }\right)$ estimation, where H_oh=87 km is the height of OH airglow layer. However, our observation does not satisfy the free oblique propagation dispersion theory of AGWs. In addition, we did not find that the GW horizontal wavelength increased with the distance from the shore, as predicted by the theory of AGW oblique propagation. Therefore, the AGWs excited by the tsunami we observed in the mesopause region may be modulated by a duct.

We did find a duct structure between 80 and 93 km (solid black line in Fig. 9b), while wave packet no. 3 was in a state of free propagation when it propagated around the coastal vicinity of the Chinese mainland (dotted line and dashed line). The duct almost includes the whole OH airglow layer. Therefore, we believe that AGWs generated by TITVE may enter the duct in the process of propagation over the Chinese mainland. The duct structure over the Chinese mainland can explain that the GWs generated by the tsunamis can propagate thousands of kilometers inland.

Figure 13 (a) Backward ray-tracing results of the fourth and fifth groups of GWs observed by the OH airglow network. The red triangles and red crosses represent the trace start and termination points, respectively. Panels (b) and (c) show the simulated tsunami directly induced by the Tonga volcano eruption (TITVE), corresponding to the dotted rectangular area in panel (a). (c) Ray paths of the wave starting from the seven sampling points in panel (a).

Figure 13 shows the results of ray tracing for wave packet nos. 4–5. The background field used for ray-tracing analysis for wave packet nos. 4–5 is from Fig. 12. The meteor radar wind field is from Ledong station (18.3° N, 109.4° E). The horizontal wavelength of wave packet nos. 4 and 5 is observed near the coast by the OH airglow network at approximately 89 ± 6 km and 80 ± 4 km. We find that the source location of AGWs over the coast of the Chinese mainland falls in the near-tsunami area, while the location of AGW ray termination over the inland is around 80 km (positions B6 and B7 in Fig. 13d), which indicates that the wave meets the evanescent layer (Wrasse et al., 2006). This is consistent with the duct structure obtained through dispersion relation. Therefore, we suggest that TITVE interact with the atmosphere after arriving at the coast of the Chinese mainland to generate the upward-propagating AGW packet. After reaching the mesopause region, this wave packet enters the wave duct structure in the horizontal propagation process, and this wave duct supports wave packet no. 5 to propagate more than 3000 km inland China.

4 Conclusions

Strong atmospheric disturbances, including Lamb waves, acoustic waves, and GWs, were triggered by the 2022 HTHH volcano eruption. The HTHH submarine volcanic eruption also triggered an unusual tsunami, which can generate atmospheric gravity waves (Fig. 14). We observed five strong groups of atmospheric waves associated with the HTHH volcano eruption from the ground-based airglow network observations.

Figure 14 The Tonga volcano eruptions triggered two types of tsunamis: one type of tsunami is induced by the atmospheric-pressure wave (TIAPW), and the other type tsunami is directly induced by the Tonga volcano eruption (TITVE). The acoustic gravity waves (AGWs) caused by tsunamis can propagate to the mesopause region.

The phase speed of the wave packet no. 1 leading front is approximately 309 m s⁻¹, which is observed almost simultaneously with the surface Lamb wave L0 mode. The high-frequency wave trains following the wave packet no. 1 leading front observed by the northern OH airglow imager network may also be related to the dissipation of the leading waves. Wave packet no. 2, with an average phase speed of 236 m s⁻¹, may be considered Lamb wave L1 mode, which exhibits internal GW behavior. Wave packet no. 3 and wave packet nos. 4–5 are generated by TIAPW and TITVE from backward ray-tracing analysis. The horizontal phase speed distribution range of wave packet nos. 3–5 is 200 to 215 m s⁻¹, which is smaller than that of wave packet nos. 1–2. For amplitude, the average amplitude of the Lamb wave L1 mode (5.4 %) is higher than that of the Lamb wave L0 mode (3.2 %), while wave packet nos. 3, 4, and 5 have relatively small amplitudes, mainly distributed between 0.85 % and 1.25 %. The horizontal wavelengths of the atmospheric AGWs observed by the airglow network are very consistent with those of the tsunami near the coast. This is the first time that we observed the AGWs in the mesopause region triggered by the tsunamis using optical detection equipment. It is also the first time to report atmospheric gravity waves excited by TIAPW.

When the wave excited by TITVE propagates far away from the coast, the characteristics of AGWs are not consistent with the dispersion of free-propagation AGWs. We find these wave packets are controlled by the duct, which can support the propagation of these GWs for thousands of kilometers after the tsunami was stopped at the coast. Therefore, tsunamis can have a significant impact on the upper atmosphere over inland areas far from the ocean through AGWs.

The 2022 HTHH volcano eruption formed a complex coupling relationship in the land–ocean–atmosphere system (Fig. 14). First, the heat released by the eruption has a direct impact on the ocean, causing temperature changes in the surrounding waters. This can lead to changes in the marine environment, affecting the behavior, distribution, and ecosystem structure of organisms.

Meanwhile, volcanoes release gases such as carbon dioxide and sulfur dioxide. Carbon dioxide is one of the greenhouse gases that can cause an increase in the Earth's temperature, leading to global warming. Sulfur dioxide can cause sulfuric acid mist in the atmosphere, which affects the reflectivity and temperature of the atmosphere and thus affects the global climate.

Moreover, the 2022 HTHH volcano eruptions also trigger atmospheric waves and tsunamis. The surface atmospheric-pressure wave generated by the 2022 HTHH volcano eruption can affect the upper atmosphere. The conventional tsunami triggered by the Tonga volcano generated AGWs. The atmospheric-pressure wave from the eruption generated a fast tsunami never before observed by tsunami observation networks. When the tsunamis reach the coast, their speeds decrease but their amplitudes increase, and the AGWs generated by them will also affect the upper atmosphere. These AGWs play an important coupling role between the ocean and the atmosphere by affecting the density and pressure distribution of the atmosphere during propagation, leading to changes in the wind field and affecting global atmospheric circulation. This study exhibits special dynamic coupling process between the air and sea via acoustic gravity waves (Fig. 14). This indirect impact on the upper atmosphere provides a new perspective for us to study the coupling between the ocean and the atmosphere and a key opportunity to improve the air–sea coupling model, thereby enhancing our future ability to make tsunami warning forecasts.