Simulation study of atmosphere–ionosphere variations driven by the eruption of Hunga Tonga-Hunga Ha’apai on 15 January 2022

The volcano of Hunga Tonga-Hunga Ha ’ apai in Tonga erupted on 15 January 2022, generating severe disturbances in the atmosphere and the ionosphere. This event provided us with large amount of data of the atmosphere and the ionosphere, and various kinds of observational studies have been made. Recently several simulation studies have also been made to reproduce and understand the atmosphere–ionosphere variations driven by the volcanic eruption. Although the simulation studies have reproduced the global variations of the atmosphere and the iono-sphere successfully, phenomena related with acoustic waves have not been fully investigated. We employed an axisymmetric three-dimensional nonhydrostatic atmospheric model and the whole atmosphere–ionosphere coupled model GAIA. We found that the simulation can produce various kinds of atmospheric waves generated by the eruption, such as acoustic waves, gravity waves, Lamb waves, Pekeris waves, and TIDs concentrically propagating from the eruption site, and atmospheric oscillations with a period of a few minutes. In addition, the results indicate that the eruption generates supersonic shock waves in the volcanic region, leading to the extremely large vertical oscillations in the thermosphere and ionosphere above the volcanic eruption region.


Introduction
It is widely recognized that a volcanic eruption triggers atmospheric waves such as acoustic waves, gravity waves, and Lamb waves.Some of the atmospheric waves propagate upward and cause thermospheric and ionospheric variations.On 15th January 2022 at about 04:15 UT an unprecedentedly explosive volcanic eruption occurred at the Hunga Tonga-Hunga Ha'apai (hereafter HTHH) (20.5° S, 175.4°W), generating intense atmospheric and ionospheric disturbances.Before the volcanic eruption at HTHH, the amount and quality of atmospheric and ionospheric data associated with volcanic eruptions were not sufficient to quantitatively investigate the physical processes relevant to volcanic eruptions.The scale of the HTHH volcanic eruption was so large that the event provided enough dataset to quantitatively study the atmosphere-ionosphere variations driven by volcanic eruptions.
One of the most prominent phenomena in the atmosphere in this event is generation and propagation of Lamb waves triggered by the explosive eruption.Various numerical simulations have been performed using global atmospheric models (Amores et al. 2022;Watanabe et al. 2022;Nishikawa et al. 2022), and successfully reproduced the propagation of Lamb waves.Besides Lamb waves, Watanabe et al. (2022) detected Pekeris waves by analyzing data taken from the Himawari-8 geostationary satellite for the first time.They also performed numerical simulations and reproduced the propagation of Lamb wave and Pekeris wave.The Pekeris wave mentioned in Watanabe et al. (2022) can be interpreted as L' 1 mode of Lamb wave (Francis 1975).However, these atmospheric models did not include the upper atmosphere, and it was not clear how those waves affect the thermosphere and the ionosphere.
To investigate the response of the thermosphere and the ionosphere to the eruption at HTHH, Liu et al. (2023a) employed high-resolution Whole Atmosphere Community Climate Model with thermosphere/ionosphere extension (WACCM-X), and successfully reproduced L' 0 and L' 1 modes of Lamb waves.Their result indicates that the amplitudes of those waves grow exponentially with altitudes, causing significant variations in the ionosphere.They found that the simulated traveling ionospheric disturbances (TIDs) in the total electron contents (TECs) are in good agreement with observations.Wu et al. (2023) used the thermosphere-ionosphere electrodynamics general circulation model (TIEGCM) with output from high-resolution WACCM-X.In their model, high-resolution grids were implemented in a regional domain in the global low-resolution mesh.The model also reproduced the observed propagation of gravity waves and TIDs.Huba et al. (2023) used a coupled model of Sami3 is Also a Model of the Ionosphere (SAMI3) ionosphere/ plasmasphere model and the HIgh Altitude Mechanistic general Circulation Model whole atmosphere (HIA-MCM).Primary atmospheric gravity wave effects were taken from the Model for gravity wavE SOurces, Ray-trAcing and reConstruction model (MESORAC).They demonstrated that the eruption could produce a largescale equatorial plasma bubble (EPB).They insisted that the EPBs were caused by the atmospheric disturbance associated with the eruption, and by the effect of the dynamo electric field.Vadas et al. (2023a, b) also studied the gravity waves and traveling ionospheric disturbances (TIDs) generated by the eruption.They employed three kinds of models: MESORAC, HIAMCM, and SAMI3.They showed that at first primary gravity waves are generated by the eruption, and then secondary gravity waves are generated through forcing and heating processes in dissipation of the primary gravity waves.By comparing the simulation results with GPS/TEC and ionosonde observations, they found that the observed TIDs are caused by the secondary gravity waves, and that the horizontal phase speeds of the TIDs range 100-750 m/s, horizontal wavelengths are about 600-6000 km.They insisted that Lamb waves do not play a primary role in generating the TIDs because the relative amplitude of Lamb waves decreases as compared with that of gravity waves in the thermosphere.Miyoshi and Shinagawa (2023) used Ground-totopside Atmosphere and Ionosphere for Aeronomy (GAIA) to study atmospheric and ionospheric perturbations triggered by the eruption.They showed that ionospheric perturbations are caused by neutral wind perturbations associated with gravity waves.Gravity waves with horizontal phase speeds of 200-310 m/s are excited in the troposphere in the HTHH eruption, and propagate upward into the thermosphere.
The simulation studies mentioned above indicate that the volcanic eruption at HTHH triggers Lamb waves and gravity waves, which propagate upward and generate TIDs in the ionosphere.General features of the TIDs are reproduced well by the simulations.However, there are still some problems with current simulation studies.Most of the models are based on hydrostatic models for the atmosphere, which does not include the effect of compressibility in the vertical direction.In volcanic eruptions, explosive shock waves are generated initially, and acoustic waves propagate in all directions.Geomagnetic field observations suggest that atmospheric oscillations caused by the acoustic resonance should have occurred in association with volcanic eruptions (Iyemori et al. 2022;Yamazaki et al. 2022).Some effects of compressibility are included in the model of Vadas et al. (2023a, b), but vertically propagating acoustic waves do not clearly appear in their results.Another problem with the current numerical models is that inputs for the volcanic eruption are given to rather extensive region with a radius of more than several tens of kilometers instead of the volcanic region.Since the size of the volcanic vent is less than several kilometers (Cronin et al. 2017), it is necessary to give a large energy input in nearly a pinpoint in a simulation domain.
Although nonhydrostatic atmospheric models which include the upper atmosphere have been constructed to investigate the effects of earthquakes and tsunamis phenomena on the upper atmosphere and the ionosphere (e.g., Shinagawa et al. 2007Shinagawa et al. , 2013;;Matsumura et al. 2011Matsumura et al. , 2012;;Zettergren et al. 2017;Zettergren and Snively 2019), no nonhydrostatic atmospheric models have been constructed for the study of atmospheric disturbances driven by volcanic eruptions.
Despite a number of observational and simulation studies, it still has not been fully understood how the volcanic eruption at HTHH leads to the generation of the observed atmospheric waves and how the atmospheric waves generate the observed ionospheric variations.In this study, we present a simulation study of atmospheric and ionospheric perturbations associated with the HTHH eruption using a nonhydrostatic atmospheric model by giving nearly a pinpoint volcanic input at the eruption site.We also use a global atmosphere-ionosphere coupled model GAIA to study the response of the ionosphere to the eruption.Emphasis is placed on early stage of the atmospheric and ionospheric disturbances caused by the volcanic eruption at HTHH.

Description of the models and simulation
We use two kinds of models to study atmosphere-ionosphere variations driven by the volcanic eruption at HTHH.At first, a simulation using a nonhydrostatic atmosphere model (NHM) with an eruption input is carried out.Then, simulations using GAIA are carried out by incorporating the output velocities obtained by NHM to produce ionosphere and atmosphere variations driven by the eruption.GAIA run without the eruption input was also done to examine effects of the eruption.A brief description of the models used in this study is given in the following sections.

Nonhydrostatic atmosphere model
We used an axisymmetric three-dimensional NHM to study variations in the atmosphere caused by the volcanic eruption at HTHH.We adopted the axisymmetric spherical coordinates (Additional file 1: Figure S1), in which HTHH is located at the pole of the coordinates.Assumption of the axial symmetry seems appropriate because atmospheric perturbations caused by a volcanic eruption propagate approximately isotropically in the horizontal direction.
A two-dimensional version of the model was originally developed by Shinagawa et al. (2003) and by Shinagawa and Oyama (2006) to study thermospheric variations in the vicinity of an auroral arc.Then, the model was used to reproduce atmospheric-ionospheric variations caused by sudden sea surface uplift of the tsunami triggered by the Sumatra earthquake on 26 December 2004 (Shinagawa et al. 2007).This nonhydrostatic atmospheric model was also employed by Matsumura et al. (2011), Matsumura et al. (2012) and by Shinagawa et al. (2013) to investigate the upper atmosphere responses to atmospheric disturbances generated by the tsunami triggered by the Tohoku-oki earthquake on 11 March 2011.
In the present study, we solved the following continuity, momentum, energy equations, and equation of state, based on the single-fluid approximation, in which change of composition ratios of atmospheric species due to mutual diffusion is neglected.
Continuity equation is, where ρ is mass density, r and θ denote radial (vertical) and angular (horizontal) coordinates, respectively, as shown in Figure S1 in Additional file 1. v r and v θ are radial and horizontal velocities, respectively.To be exact, the first term in the right-hand side of Eq. ( 1) is − ∂ r 2 ∂r ρv r r 2 .However, it can be approximated by − ∂ ∂r (ρv r ) when the altitude (1) region is much smaller than the radius of the earth.This approximation is also used in the energy Eq. ( 5).In the momentum equation, Coriolis force, centrifugal force, and ion-neutral drag force are neglected because those terms are much smaller than other terms in the present simulation.The momentum equation used in this study is, where is the velocity vector, g is the gravity acceleration, p is pressure, µ is molecular viscosity.The value of µ is taken from Rees (1989).In this model, eddy viscosity is neglected.Although the eddy viscosity is small under normal atmospheric conditions, it could be significant under highly disturbed conditions driven by the volcanic eruption.This process may play an important role in the thermospheric-ionospheric variations driven by the eruption.
Assuming axial symmetry, the radial and angular components of the momentum equations are written as follows: The energy equation is, where T is neutral temperature, κ is heat conduction coefficient, Q is external heating.In the energy equation, viscous heating, and neutral-ion frictional heating are neglected.The value of κ is taken from Rees (1989).
The equation of state: where R is the gas constant. ( (4) (5) In the simulation, only perturbation quantities (ρ' , v r ' , v θ ' , T' , p') are calculated, assuming that the background quantities (ρ 0 , T 0 , p 0 ) are in equilibrium state and horizontally uniform, and that the background atmosphere is at rest, that is, v r0 = v θ0 = 0.Under those assumptions, the vertical momentum equation, the energy equation, and the equation of state in the background state are written as, Using the assumptions ( 7)-( 9), the perturbation quantities are written as follows: The radius of the earth (r 0 ) is set to 6370 km.The altitude region in the model is from 0 km altitude (the earth's surface) to 700 km altitude.The horizontal region is from θ = 0° (the eruption point of HTHH) to θ = 180° (the opposite side of the Earth) (Additional file 1: Figure S1).The vertical grid size is 2 km, and the horizontal grid size is 0.18°, which is about 20 km.The time step in this numerical simulation is 0.1 s.
The background atmosphere is obtained from NRLM-SIS 2.1 (Emmert et al. 2022), for which the time of the ( 7) (13) eruption, the location of HTHH and the F10.7 index (= 115.6) on 15 January are used as input parameters.Vertical profiles of the background neutral densities and the neutral temperature are shown in Figs.1a and b.

GAIA model
GAIA is a global atmosphere-ionosphere coupled model covering the entire atmospheric region from the ground surface to the upper thermosphere and the ionosphere (Jin et al. 2011).In Miyoshi and Shinagawa (2023), GAIA was used to study atmospheric and ionospheric perturbations triggered by the volcanic eruption at HTHH.Lamb waves and gravity waves driven by the eruption and ionospheric variations associated with those waves were reproduced.Lamb waves reproduced with GAIA are consistent with those reproduced by the model of Liu et al. (2023a).However, there are some differences between thermospheric waves in GAIA and those in other models.GAIA produces gravity waves with a phase speed of 310 m/s or less, while faster gravity waves than 310 m/s are much weaker (Miyoshi and Shinagawa 2023).Simulation studies by Vadas et al. (2023a, b) show gravity waves with several phase speeds between 100 and 600 m/s.Liu et al. (2023a) indicate that the dominant atmospheric waves are mainly Lamb waves (L' 0 and L' 1 modes) generated by the volcanic eruption, while some gravity waves with phase speeds less than 900 m/s are also produced.
Since hydrostatic atmospheric models do not include compressibility in the vertical direction, shock waves and acoustic waves generated by the eruption cannot be reproduced.To include the effect of compressibility, we incorporated the velocity perturbations obtained from the NHM simulation into the GAIA simulation, which enable us to reproduce explosive atmospheric disturbances and acoustic waves triggered by the volcanic eruption at HTHH.In the present study, the model setting of GAIA in this study is the same as one used in Miyoshi and Shinagawa (2023).

Simulation of atmospheric-ionospheric variations
We first run the NHM model with an eruption input.The eruption input is introduced as a temperature increase at the point of the HTHH volcano.Since actual temperature variation during the eruption could not be accurately determined, a temperature variation shown in Fig. 2 was given.The maximum temperature increase is assumed to be about 370 K, and the temperature returns to the background temperature in about 10 min, which was given only from the eruption point (the surface) to 2 km above the eruption point.
Previous studies of volcanic eruptions have indicated that the temperature just above a volcanic vent could reach more than 1000 K (e.g., Ogden et al. 2008a;b).Since our model does not have enough resolution to strictly reproduce actual physical processes around a volcano, we simply added the temperature increase shown in Fig. 2 at the point of HTHH, which turned out to produce fairly realistic atmospheric and ionospheric variations driven by the volcanic eruption at HTHH.However, further examination about the input will be necessary for more detailed analyses.The simulation was carried out for 24 h starting from 04:15 UT 15 January 2022.The calculated vertical and horizontal wind velocity data at every 10 s are saved for GAIA runs.
In the next step, simulations using GAIA were carried out by incorporating the velocity variation data obtained by the NHM calculation.In the GAIA simulation, perturbation velocities of NHM were added to the background neutral wind in GAIA starting from 04:15 UT.The neutral densities and temperature changes in NHM are not incorporated into GAIA.The changes in the neutral density and temperature may also affect the ionosphere, although the neutral wind variations are likely to be the main drivers of the ionospheric variations.
As was done in Miyoshi and Shinagawa (2023), meteorological reanalysis data (JRA-55) provided by Japan Meteorological Agency (Kobayashi et al. 2015;Harada et al. 2016) were incorporated below 40 km altitude by a nudging method.The JRA-55 fields are updated every 6 h, and the effect of the HTHH eruption is not included in the JRA-55 fields.The F10.7 index of the eruption day was used as a proxy for the solar EUV model.By running GAIA, global atmospheric and ionospheric variations on the eruption day were obtained.A GAIA simulation without the perturbation velocities of NHM was also carried out for the comparison with the eruption run.

Results
The NHM run starts with a sudden temperature increase at the point of the eruption described above.In our simulation, it is assumed that the heating starts at 04:15 UT.We first present the initial response of the atmosphere to the eruption obtained from NHM in "Initial responses of the atmosphere to the volcanic eruption" section.Then, the behavior of atmospheric waves generated by the eruption is described in "Propagation of atmospheric waves" section.Finally, ionospheric variations obtained by GAIA simulations are shown in "Ionospheric variations" section.

Initial responses of the atmosphere to the volcanic eruption
Figure 3 shows the initial vertical wind perturbation obtained by the NHM run.Additional file 2 shows an animation of the vertical wind perturbation.The horizontal axis is the distance from HTHH, and the vertical axis is the altitude.Red color indicates upward, and blue color indicates downward.Figure 3a is the vertical wind velocity 5 min after the eruption, revealing that shock acoustic waves are generated and propagated vertically and horizontally.The front of the wave reaches the lower thermosphere.The phase speed of the initial acoustic wave is about 450 m/s, which is supersonic.The front of the shock reaches about 500 km altitude in 10 min after the eruption (Fig. 3b), in which the vertical phase speed is about 800 m/s.The fluid velocity also reaches about 500 m/s in the thermosphere.Generation and propagation of the supersonic shock waves triggered by an explosive volcanic eruption have been investigated by a number of observational and theoretical studies, stating that supersonic shock waves can be generated when the perturbation pressure exceeds the background atmospheric pressure (e.g., Kieffer 1981;Ogden et al. 2008a, b;Koyaguchi et al. 2010;Medici et al. 2014;Watson et al. 2021).Shock waves triggered by tsunamis also have been observed in the Tohoku-oki earthquake occurred on 11 March 2011 (e.g.,  Astafyeva et al. 2011).At 04:45 UT, the shock waves reach the upper thermosphere, and produce secondary gravity waves, propagating horizontally at a speed of about 750 m/s (Fig. 3c).At 05:15 UT, the amplitude of the vertical wind decreases rapidly as the waves move away from HTHH (Fig. 3d).
Figure 4 shows horizontal velocities from 04:20 UT to 05:15 UT.Additional file 3 shows an animation of the horizontal wind perturbation.Red color indicates outward from HTHH, and blue color indicates inward.Outward velocity appears just after the eruption, and the shock waves rapidly propagate upward into the thermosphere (Fig. 4a).After 04:45 UT, gravity wave structures develop (Fig. 4b and c).While the vertical wind velocity decreases rapidly, the horizontal velocity decreases rather gradually (Fig. 4d).The gravity waves have several different phase speeds: the fastest wave moves at a speed of about 750 m/s, and some slower waves follow.The propagation speed will be discussed later in this paper.
Figure 5 shows the vertical velocity variation at altitudes of 120 km and 200 km above the eruption site between 04:15 and 05 UT.Initially, very large vertical oscillation is driven by the acoustic shock waves.After 4.6 UT, oscillation with a period of about 4 min is clearly seen.Recently Nakata et al. ( 2023) have detected oscillations of the ionosphere with a period of about 4 min associated with the eruption of HTHH using a HF Doppler sounding system.This is interpreted as resonances of acoustic mode waves propagating between the ground and the lower ionosphere.This type of oscillation was observed by geomagnetic field measurements near HTHH (Iyemori et al. 2022;Yamazaki et al. 2022).Similar oscillations in TEC or the geomagnetic field related with the acoustic resonance were also observed in past volcanic eruptions (Kanamori et al. 1994;Iyemori et al. 2005;Dautermann et al. 2009;Aoyama et al. 2016;Nakashima et al. 2016;Heki and Fujimoto 2022).
It was also found that the shape of the geomagnetic field oscillation associated with the eruption at HTHH has a "saw-tooth" form (Iyemori et al. 2022).Iyemori et al. (2022) interpreted this saw-tooth variation as

Propagation of atmospheric waves
Immediately after the volcanic eruption at HTHH, acoustic-gravity waves and Lamb waves are generated.Figure 6 shows the horizontal neutral wind velocity (Fig. 6a) and the vertical neutral wind velocity (Fig. 6b) at 11:15 UT, that is, 7 h after the eruption.To illustrate small velocities clearly, logarithms of the velocities in the unit of cm/ sec are taken.For example, "0" indicates 1 cm/s, and "2" indicates 1 m/s.Additional files 4 and 5 show animations of the horizontal and vertical wind perturbations, respectively.
In the thermosphere, large-scale gravity waves are clearly seen.The front of the fastest gravity wave has a speed of about 750 m/s, followed by several slower gravity waves.Below 100 km Lamb wave is seen at about 7800 km away from HTHH.In addition, Pekeris wave (L' 1 mode of Lamb wave) is also seen at 6000 km.Locations of Lamb wave and Pekeris wave are indicated in Fig. 6a as "L" and "P", respectively.Lamb and Pekeris waves are identified clearly as sudden surface pressure jumps (Additional file 1: Figure S2).The phase speed of Lamb wave is about 310 m/s, and that of Pekeris wave is about 240 m/s, which are approximately consistent with the results in previous simulation studies (Amores et al. 2022;Watanabe et al. 2022;Suzuki et al. 2023;Liu et al. 2023a).It is noted that both Lamb wave and Pekeris wave appear to have some influence on the neural wind even above 100 km, modifying the gravity wave structures.Global plots of the perturbation horizontal radial velocities at an altitude of 300 km are shown in Fig. 7, that is, the outward speed directed from the HTHH volcano.Since our NHM is assumed to be axisymmetric, the structures of the velocities are exactly concentric.As mentioned before, the fastest gravity wave moves at a speed of about 750 m/s, and the front of the wave reaches Japan at about 08:15 UT (4 h after the eruption).But the magnitude of the amplitude decreases rapidly, and the leading gravity wave would not affect the ionosphere above Japan very much.At about 11:15 UT gravity waves with larger amplitudes reach Japan, which could significantly affect the ionosphere.
Global plots of the variations in the vertical velocities at an altitude of 300 km are shown in Fig. 8.As shown in Fig. 5, the vertical velocity is initially very large near HTHH (Fig. 8a).The significant vertical wind variations start in Australia at about 08 UT (Fig. 8b), and in Japan at about 11 UT (Fig. 8c), which corresponds to the arrival of TIDs in those regions (e.g., Carter et al. 2023;Saito 2022;Lin et al. 2022;Shinbori et al. 2022).The vertical velocity of the gravity waves decreases rapidly as they move away from HTHH, and becomes insignificant at about 14 UT (Fig. 8d).

Ionospheric variations
To reproduce ionospheric variations associated with the volcanic eruption at HTHH, neutral wind perturbations obtained by NHM were incorporated into the GAIA model.Figure 9 shows the global distribution of variations in TEC (total electron contents), which is TEC (with eruption) minus TEC (without eruption).At 05 UT concentrically propagating TIDs are clearly seen (Fig. 9a).At 08 UT the concentric structures are still seen partially, but large-scale variations in TEC begin to develop (Fig. 9b).It seems that the concentric TIDs are covered by the large-scale TEC perturbations caused by the effect of the wind-driven dynamo electric field.At 11 UT the concentric TID structure is not clearly seen (Fig. 9c).While significant increase in TEC reaches Japan in the simulation, observations indicate more prominent TID patterns (Saito 2022;Shinbori et al. 2022;Lin et al. 2022).At 14 UT the concentric TIDs are no longer seen and global variations in TEC are seen.This is because changes in the neutral wind cover all over the world at 14 UT (Fig. 7) and because the ionospheric electric field is modified by the neutral wind-driven dynamo.
Although our simulation reproduced the concentric TIDs directly propagating from HTHH, it could not reproduce the observed concentric TIDs in geomagnetic conjugate region that are related to the directly propagating TIDs (e.g., Shinbori et al. 2022;Lin et al. 2022).Other simulation studies also could not clearly reproduce concentric TID pattern in the magnetic conjugate region (Vadas et al. 2023b;Wu et al. 2023;Liu et al. 2023a;Miyoshi and Shinagawa 2023).Figure 10 shows the global distribution of the zonal electric field, which was generated by ionospheric dynamo processes driven by the neutral wind.At 05 UT two symmetric concentric structures of the perturbed electric field are clearly seen on both sides of the magnetic equator.The concentric structure extends outward as the neutral wind perturbation structure extends outward (Fig. 8).At 08 UT, equator sides of the concentric electric field in both hemispheres collide with each other, forming an X-type electric field structure at the equatorial region (Fig. 10b).In the geomagnetic equator, the electric field is enhanced, which is likely to be caused by the increased neutral wind (Fig. 7b).
In the mid-latitude regions after 08 UT, concentric structures continue to extend, and reach Japan and Australia.The electric field variation causes E × B drift and produces concentric structure of TEC.As mentioned above, however, such structures are not seen in the simulated TEC variations (Fig. 9).This is because the contribution of the E × B drift to the electron density change is much smaller than that of the local neutral wind in the simulation.The enhanced electric field regions also show wave-like structure along the equator at 11 UT (Fig. 10c), propagating away from HTHH.The region of the electric field increased eastward moves westward, while the region of the electric field increased westward moves eastward.
Figure 11 shows altitude profiles of the electron density variations at the eruption site.From 04 UT to about 06 UT significant oscillation of the ionosphere occurs (Fig. 11a), indicating the effect of the acoustic resonance between the ground and the mesopause region.Previous studies suggest that the resonant atmospheric oscillations leak to the upper thermosphere, generating the oscillations of the ionosphere (e.g., Iyemori et al. 2022;Heki and Fujimoto 2022).After 06 UT the oscillation decreases, and the electron density gradually increases in the upper thermosphere (Fig. 11a).This is mainly due to the upward motion of the ionosphere which was caused by the upward E × B drift motion with the eastward electric field around the HTHH region (Fig. 10c).
The ionospheric observations also indicated significant decrease in TEC in HTHH region after the eruption (Aa et al. 2022a;Astafyeva et al. 2022;Sun et al. 2022b;He et al. 2023).Our simulation does not exhibit such a large decrease in electron density near HTHH, although moderate decrease occurred at altitudes of 300-400 km from 06 to 13 UT.This discrepancy is probably due to the input of the eruption.In the HTHH volcanic eruption a series of eruptions occurred (Astafyeva et al. 2022;Purkis et al. 2023) instead of a single eruption used in our simulation (Fig. 2).This long-lasting eruption might have caused large thermospheric wind perturbation and the TEC decrease near the HTHH region.
Figure 11b shows the initial part of the electron density oscillation, which is driven by the oscillatory motion in the thermosphere (Fig. 5).The period of the oscillation is about 4 min.Oscillations of the ionosphere with a period of about 4 min associated with the eruption of HTHH have been detected by a HF Doppler sounding system (Nakata et al. 2023).Geomagnetic field observations also strongly suggest the acoustic resonance driven by the eruption (Iyemori et al. 2022;Yamazaki et al. 2022;Schnepf et al. 2022).

Cause of TIDs
Recent numerical simulation studies of the atmospheric variation driven by the HTHH volcano have successfully reproduced the Lamb waves and Pekeris waves (L' 1 mode of the Lamb waves) (e.g., Watanabe et al. 2022;Liu et al. 2023a, b).The ionosphere-atmosphere coupled models have indicated that Lamb waves and/or gravity waves triggered by the eruption generate TIDs, which propagate concentrically from HTHH (Wu et al. 2023;Liu et al. 2023a;Miyoshi and Shinagawa 2023;Vadas et al 2023b).As for the large-scale variations in the ionosphere and thermosphere, our simulation result is basically the same as those studies.However, there is still an uncertainty about the triggering mechanism of the observed TIDs.While Liu et al. (2023a) suggested that the TIDs are generated by Lamb waves in the thermosphere, Vadas et al. (2023b) insisted that the observed TIDs can be explained by secondary gravity waves driven by the gravity waves triggered by the eruption rather than by the "leaked" Lamb waves in the thermosphere.Miyoshi and Shinagawa (2023) also indicated that the TIDs are generated by gravity waves in the thermosphere rather than Lamb waves.The results in the present study indicate that both Lamb wave and Pekeris wave locally exert significant influence on gravity waves in the thermosphere (Fig. 6), causing TEC variations in the region of Lamb wave and Pekeris wave.The interaction between the gravity waves and the leaked Lamb waves seems important and requires further investigation.

Acoustic waves generated by the eruption
Our simulations also demonstrated that the explosive eruption of the volcano generates the initial acoustic shock waves and subsequent acoustic resonance with a period of about 4 min.Although the oscillation has not been directly detected in the TEC observations in the HTHH region, clear signature of the magnetic field oscillation with the same period was detected (Iyemori et al. 2022;Yamazaki et al. 2022;Schnepf et al. 2022).Since oscillations in TEC were observed in previous volcanic eruption events (Kanamori et al. 1994;Dautermann et al. 2009;Nakashima et al. 2016;Heki and Fujimoto 2022), it is likely that similar oscillations of acoustic resonance should have occurred in the HTHH eruption event.

Propagation of gravity waves
Figure 12 shows the time variations in neutral wind velocities at 300 km altitude and at 175° W between latitudes of 50° S and 10° N and between 04 and 08 UT to illustrate the gravity waves propagating in the thermosphere after the eruption.The vertical wind velocity is shown in Fig. 12a.The first large vertical perturbation occurred at about 04:23 UT, which is several minutes after the eruption.Near the eruption site, prominent acoustic resonance is seen from 04:30 UT to about 06 UT as shown in Fig. 5.The gravity waves have several different phase speeds, propagating away from HTHH, and the amplitude of the vertical velocity decreases rather rapidly with time.Figure 12b shows the horizontal wind variations at 300 km altitude.As seen in Fig. 12a, the gravity waves have several different phase speeds, which are indicated by dashed lines in Fig. 12b.The wind variation begins at 04:23 UT.The fastest phase speed of the gravity waves is about 750 m/s, followed by a few slower gravity waves with a speed of 450 m/s, 300 m/s, 250 m/s, and 200 m/s.The gravity wave with a speed of 300 m/s appears to be related with Lamb waves in the lower atmosphere, and the gravity waves with a speed of 250 m/s corresponds to Pekeris waves.Those two waves propagate over long distance, while gravity waves with other speeds attenuate rapidly as they move away from the eruption site.
Miyoshi and Shinagawa (2023) also showed phase speeds of gravity waves obtained by GAIA simulations, which reproduced gravity waves propagating with phase speeds of 310 m/s, 200 m/s, and about 250 m/s.However, gravity waves propagating faster than 310 m/s did not appear in their simulations, which is different from the present study.The present simulation indicates that the faster gravity waves are secondary gravity waves generated by initial acoustic waves triggered by the eruption at HTHH.

Propagation of TIDs
Figure 13 shows the time variations in TEC at 175˚W and latitudes between 60° S and 40° N. TIDs propagating at several different speeds are seen.Near HTHH, the initial propagation speed appears to exceed 1000 m/s.However, this speed does not mean the actual horizontal movement of the structure because the initial acoustic shock wave propagates to all directions from HTHH (Fig. 3).In the region near HTHH, the front of the acoustic waves reaches 300 km altitude almost at the same time as the point just above HTHH, leading to the very fast apparent speed of propagation of the disturbance.As the TIDs move away from HTHH, the TID speeds are settled, and reach several typical propagation speeds appear: the fastest one is 750 m/s followed by TIDs at speeds of 520 m/s, 380 m/s, 310 m/s, 250 m/s and 200 m/s.
In the northern hemisphere between 10° N and 40° N, a small TEC increase occurs separately starting from about 04:30 UT.This is caused by the electric field variations transmitted from the magnetic conjugate region where significant electric field is generated by the neutral wind dynamo processes (Fig. 10).But as mentioned above the TID-like structure is not clearly seen.This TEC variation extends to the southern hemisphere, and overlap with the TEC variations propagating from HTHH.Near HTHH oscillations with a period of 4 min occur.

TIDs in geomagnetic conjugate
One of the most intriguing phenomena in the eruption event is that the TIDs are observed at about 08 UT in Japan, which is only 4 h after the eruption (Saito 2022;Shinbori et al. 2022;Lin et al. 2022).As mentioned above, the fastest gravity wave reaches Japan at about 08 UT.However, the amplitude of the wave is very small, and the horizontal wavelength is very large.Therefore, it is unlikely to generate clear TIDs.Several authors have indicated that it is caused by the electric field transmitted from the magnetic conjugate region, where TIDs are passing Australia.Although the electric field in our simulation implies the process, the TEC variations driven by the ExB drift are not significant.It seems that simulated TIDs in other studies also could not successfully reproduce such structures.For this problem, further investigation is necessary.

Uncertainties in the present study
Although our model reproduces some of the ionospheric-atmospheric phenomena associated with the eruption at HTHH, there some uncertainties in the present simulations.First, the input of the volcanic eruption is rather arbitrary because the exact physical quantities in the eruption were not obtained.In a volcanic eruption, explosive increase of pressure and temperature occurs.In addition, significant amount of water vapor is emitted from the volcano.These effects could not be treated quantitatively at present.It is also pointed out that the eruption occurred not once, but several times successively.
The present study also has some problems with the modeling method: (1) the background atmosphere in NHM is assumed to be horizontally uniform and fixed in time, which is not the same as the atmosphere in GAIA; (2) the interaction between the background wind in GAIA and eruption-generated wind in NHM is not included; (3) effects of temperature and density changes in the atmosphere on the ionosphere are not included because only velocities are incorporated in GAIA; (4) while the horizontal grid size in NHM is 20 km, that of GAIA is about 1 degree, which is not able to reproduce ionospheric structures with a scale of less than several hundred kilometers.Thus, more high-resolution and self-consistent model needs to be developed to study atmospheric and ionospheric disturbances driven by the volcanic eruption at HTHH more realistically.
Another problem is that during the eruption period, a moderate magnetic storm occurred, which might have affected the upper atmosphere (He et al. 2023).Although previous studies suggest that the effect is not very large in mid-and low-latitude regions, there might be some effects of the magnetic storm, which could modify largescale winds.

Summary
To reproduce and understand the atmosphere-ionosphere variations driven by the volcanic eruption at HTHH, we employed an axisymmetric three-dimensional nonhydrostatic atmospheric model (NHM) and the whole atmosphere-ionosphere coupled model GAIA.Variations in neutral winds calculated by NHM model were incorporated to GAIA to simulate ionospheric variations triggered by the volcanic eruption at HTHH.The following results have been obtained from the simulations: 1.The volcanic eruption at HTHH generates supersonic acoustic shock waves, leading to extremely large vertical oscillations in the thermosphere and ionosphere in the HTHH region.2. Following the initial large oscillations in the thermosphere, acoustic resonance with a period of about 4 min appears near the eruption region.This oscillation is qualitatively consistent with the observed signature of the geomagnetic field oscillation with the same period.3. The simulation reproduces various kinds of atmospheric waves generated by the eruption, such as acoustic waves, gravity waves, Lamb waves, Pekeris waves.4. Gravity waves with several different phase speeds are generated in the thermosphere.The simulated propagation speeds are roughly consistent with the observed speeds.5. TIDs concentrically propagating from the eruption site are partially reproduced by the simulation, but concentric TIDs observed in the magnetic conjugate region of the eruption site are not reproduced by the present model.6. Lamb waves and Pekeris waves modify the gravity waves to some extent in the thermosphere, and those waves could cause local TEC variations moving with those waves.7. Significant TEC decrease near the Tonga region after the eruption detected by GPS/TEC observations does not clearly appear in the simulation.
As was mentioned above, although some features of the atmosphere-thermosphere variations have been reproduced, there are still some disagreements between observations and simulations.It seems that the discrepancy is partly due to the accuracy and consistency of the model, and partly due to uncertainties with the input of the eruption.Observational studies suggest that there were multiple eruptions in the HTHH volcano (Astafyeva et al. 2022;Purkis et al. 2023), implying that the actual atmosphere-ionosphere variations are far more complicated than those triggered by a simple input.Further analysis needs to be done in the future.

Fig. 1 Fig. 2
Fig. 1 Altitude profiles of the neutral atmosphere used in NHM simulations.a Number densities, and b neutral temperature.The profiles are taken from Emmert et al. (2022) using the parameters of HTHH at the time of the volcanic eruption

Fig. 3
Fig. 3 Vertical wind perturbation obtained by the NHM simulation.a 4:20 UT, b 4:25 UT, c 4:45 UT, and d 5:15 UT.The horizontal axis is the distance from HTHH, and the vertical axis is the altitude.Red color indicates upward, and blue color indicates downward

Fig. 5
Fig. 5 Vertical neutral velocities at HTHH obtained by the NHM simulation.(Black) vertical velocity at 120 km altitude, and (red) vertical velocity at 200 km altitude.Positive sign indicates upward

Fig. 6
Fig. 6 Neutral wind velocities at 11:15 UT. a Horizontal velocity, and b vertical velocity.To illustrate small velocities clearly, logarithms of the velocities in the unit of cm/sec are taken.Positive sign indicates outward in a and upward and b, respectively.Locations of Lamb wave and Pekeris wave are indicated in a as "L" and "P", respectively

Fig. 11
Fig. 11 Altitude profiles of the electron density variations at the eruption site.a Between 4 and 20 UT, and b between 4 and 6 UT

Fig. 12
Fig. 12 Time variations in neutral wind velocities at 300 km altitude and at 175° W. a Vertical velocity, and b horizontal velocity.In a, positive sign (red) indicates upward, and negative sign (blue) indicates downward.In b, positive sign (red) indicates northward, and negative sign (blue) indicates southward.Velocities between 4 and 8 UT and between latitudes of 50° S and 10° N are shown.Dashed lines in b are estimated propagation speeds

Fig. 13
Fig. 13 Time variations in TEC variation at 175° W. TEC variations between 4 and 8 UT and between latitudes of 50° S and 10° N are shown.Dashed lines are estimated propagation speeds