Atmospheric and Ionospheric Responses to Hunga‐Tonga Volcano Eruption Simulated by WACCM‐X

High‐resolution Whole Atmosphere Community Climate Model with thermosphere/ionosphere extension is used to simulate the responses to the Hunga‐Tonga volcano eruption on 15 January 2022. Global propagation of the Lamb wave L’0 and L’1 pseudomodes are reproduced in the simulation, with the exponential growth of wave amplitudes with altitudes. The wavefront is vertical up to the lower thermosphere, and tilts outward above. These features are consistent with theoretical results. With simulated surface pressure perturbation agreeing with observations (∼100–250 Pa), thermospheric wind perturbations over 100 ms−1 are comparable with reported satellite and ground‐based observations. Traveling ionospheric disturbances in the total electron contents from the simulation show good agreement with observations, including magnitude and propagating speed and evidence of conjugacy in the first 1–2 hr after eruption. Conjugacy in E × B drift, on the other hand, is more persistent.

as L' 0 and L' 1 pseudomodes (simply referred to as modes hereafter). The former is non-dispersive propagating at acoustic phase speed, and the latter becomes nearly non-dispersive at wave periods longer than the buoy ancy period. An alternative interpretation of the large responses in the geospace system is by secondary gravity waves (Vadas et al., 2023), which were able to produce thermospheric and ionospheric signatures that agree with observations. This event provides a rare opportunity to examine the coupling of different atmosphere regions through a very fast process, and it also poses a challenge for whole atmosphere modeling (Matoza et al., 2022;Wright et al., 2022). The goal of this study is to use the newly developed high-resolution WACCM-X to simulate this event, to compare results with previous theoretical modeling results, and to quantify the wave processes that connect the various observations.
The high-resolution WACCM-X model and the simulation setup are described in Section 2. The simulation results are presented and compared with previous theoretical modeling results and observations in Section 3. A summary is given in Section 4.

Model Description
A high-resolution configuration of WACCM-X is used for this study. This configuration is based on the recent development of WACCM-X with the spectral element (SE) dynamical core. A detailed description of the model and its high-resolution configuration can be found in Liu et al. (2022). Here, we briefly describe the most important features of the model. The lower and middle atmosphere physics package in the model is the same as in WACCM version 6 (WACCM6; Gettelman et al., 2019), and the thermospheric/ionospheric physics is described in Liu et al. (2018). The SE dynamical core takes into consideration the species dependence of mean molecular weight, dry air gas constant, and specific heats. It uses a quasi-uniform cubed sphere grid , which enables simulations at higher horizontal resolutions. Molecular viscosity and thermal conductivity are now included in the horizontal direction in the SE dynamical core (those in the vertical direction are already part of the physics package). The CSLAM (conservative semi-Lagrangian multi-tracer scheme, Lauritzen et al., 2017) is used for transport by neutral wind. A regridding scheme based on Earth System Modeling Framework (Theurich et al., 2016, https://earthsystemmodeling.org/) is implemented in the model to transmit data between the physics mesh and the magnetic apex coordinate used for ionospheric calculations. The horizontal resolution of the cubed sphere used here is ∼25 km, and the vertical resolution is 0.1 scale height in most of the middle and upper atmosphere, transitioning to 0.25 scale height in the top three scale heights. There are a total of 273 levels between the Earth surface and the model top at 4 × 10 −10 hPa (∼600 km).
The simulation is initialized from a free run for January conditions . As in the free run, the solar radio flux at 10.7 cm (F10.7) is set to 120 solar flux unit, which is similar to the actual F10.7 value on 15 January 2022. The geomagnetic index Kp is set to 0.33, so the simulation does not consider the geomagnetic storm preceding the eruption on 14 January 2022 (e.g., Harding et al., 2022;Shinbori et al., 2022). This is done to isolate the responses to the eruption. The eruption is initiated with a 50 hPa surface pressure perturbation over a radius of 60 km centered at 20.5°S and 175°W, similar to what was done in the simulation by Amores et al. (2022). This is a simplified representation of the actual eruption, which was a complex sequence between UT 04:00 and UT 04:30 (Matoza et al., 2022). The model simulation is run for 44 hr after the initiation, so it covers over one passage of the Lamb wave around the globe (about 35.5 hr).

Results
Figures 1a-1d show the surface pressure perturbation (high-pass filtered with period <1 hr) from the simulation 1.5, 4, 13.5, and 18 hr after the eruption, at universal times (UTs) of 06:05, 08:35, 18:00, and 22:35 hr, respectively. These are times when the wavefront is at or near New Zealand, Hawaii, Eurasia/Africa/Atlantic, and the antipode of the eruption. While the wave phase lines form nearly perfect concentric rings at UT 06:05, the anisotropy of the wave propagation becomes more evident at later times, as seen from the shape of the wavefront and the wave amplitude variation along the phase lines. This anisotropy is likely caused by wave refraction due to changes of atmosphere temperature (thus the acoustic speed). The changes of the ring shape are in good agreement with observations (Amores et al., 2022). The amplitude of the leading Lamb mode is ∼100-250 Pa, in general agreement with barometer records at various geographical locations (  reflections are noted in the surface pressure perturbations from major mountain ranges (e.g., from the Andes at UT 18:00 hr, as seen in Figure S1 in Supporting Information S1). The reflection was seen in observations (Wright et al., 2022). The reflection from the Andes is also clearly seen in the animation of the simulation result (https:// visgallery.ucar.edu/tonga-eruption). Signature of the slower L' 1 mode is seen in the simulation. Its amplitude is between 20 and 30 Pa at UT 06:05, and drops to ∼10 Pa at UT 08:35. At UT 18:00 its signature is seen over South America, South Atlantic Ocean, Indian Ocean, and over Eurasia ( Figure S1 in Supporting Information S1). The global propagation of the wave is also seen from keogram of the surface pressure perturbation (Figures 1e and 1f). It is clear that the leading Lamb mode maintains a significant amplitude throughout the simulation, with the horizontal phase speed between 305 and 320 ms −1 . At UT 16 on January 16 (∼35.5 hr after the eruption), the wave goes around the global, though the apparent converging latitude is about 8° south of the volcano. This is probably due to wave refraction as discussed earlier. The L' 1 mode is most significant in the first 4 hr after eruption, though its phase propagation remains distinguishable till the end of 15 January. According to Francis (1973), the attenuation distance is several tens-100 Earth circumference for L' 0 mode and about half to one Earth circumference for L' 1 mode for perturbations with periods between buoyancy period and 1 hr. The global propagation of these two modes in the model simulation is thus consistent with Francis (1973).
The amplitudes of the wave perturbations grow with altitude, as can be seen in Figure 2. The amplitudes of the zonal wind, vertical wind, and temperature in the figure are obtained from the largest perturbations of these quantities (high-pass filtered with period less than 1 hr) at UT 06:05 hr. The amplitudes grow exponentially, with zonal wind from ∼1 ms −1 at the Earth surface to ∼350 ms −1 in the lower thermosphere (∼135 km), vertical wind from 0.2 to 350 ms −1 (∼300 km), temperature from 0.2 to 400 K (∼200 km). The theoretical exponential growth rates for Lamb wave and internal wave are plotted for comparison. The over two orders of magnitude increase of the zonal wind perturbation from the simulation are in good agreement with the theoretical modeling result in the presence of wave dissipation (molecular viscosity, ion drag, and surface friction; Lindzen & Blake, 1972). Its growth rate decreases above ∼135 km, similar to Lindzen and Blake (1972, Figure 11). As in the simulation, the theoretical modeling also shows larger growth rates for vertical wind and temperature perturbations than the zonal wind perturbations.
As noted in Lindzen and Blake (1972), the vertical wind perturbation is non-zero in an atmosphere with realistic background, unlike the classical Lamb wave solution in an isothermal and windless atmosphere. The horizontal wave structures of the vertical wind perturbations in the mesosphere, lower thermosphere and upper thermosphere (0.1, 2.9 × 10 −5 , and 2.6 × 10 −7 hPa) at UT 08:35 hr are shown in Figure 2. The wave amplitude increases by about one order of magnitude over each of the altitude increments. The largest wave amplitude is clearly associated with the Lamb wave (L' 0 mode), with a horizontal phase speed of about 320 ms −1 . L' 1 mode (phase speed ∼247 ms −1 ) is also evident in the plot.
The phase of the wave is also examined in Figure 2. The zonal and meridional wind perturbations (period <1 hr) at UT 06:05, normalized by p κ , are shown in longitude and latitude directions, respectively. κ is the ratio of gas constant of dry air and the specific heat of air at constant pressure. For an isothermal and windless atmosphere, the Lamb wave amplitude grows with altitude at the rate of 1/p κ . The phase of the fastest waves (Lamb waves) are perpendicular to the surface up to the lower thermosphere, and tilted outward at higher altitudes. This is consistent with the theoretical modeling results (Lindzen & Blake, 1972, Figure 4). The vertical profiles in the plot indicate the propagation distance from the epic center with the local acoustic speed for each altitude at this time. It is seen that the outer reach of the perturbations generally agree with this acoustic propagation envelope. The internal wave structure with tilted phase are also seen closer toward the epic center in the lower atmosphere, indicating slower propagating speed.
Neutral temperature and zonal wind (unfiltered) in the thermosphere are shown in Figure 3 at two different times (UT 14:05 and 18:00 hr) to demonstrate their global propagation. At UT 14:05 (9.5 hr after the eruption), global wave signatures are clearly visible in neutral temperature and zonal wind at ∼200 km: toward the east the L' 0 wavefront reaches the eastern part of North America and west coast of South America, and toward the west at Central Asia and Bay of Bengal. Over North America the temperature varies between ∼550 and 850 K, suggesting a wave amplitude of ∼150 K, and the zonal wind variation between −150 and 185 ms −1 . The location of the wavefront (at ∼78°W and 30°N) and its eastward tilt is similar to that observed by MIGHTI (Harding et al., 2022), though the observed wavefront was further eastward (60°W). To determine whether this is a systematic shift or is caused by the difference in propagating speed, the wavefront locations at UT 14:05 and 15:36 hr are compared. It is found that the longitude difference between the two wavefronts along the ICON orbit is 20° in observations (Harding et al., 2022, Figure 1) and in the simulation ( Figure S2 in Supporting Information S1). This suggests that the simulated propagating speed of the wavefront is similar to the observed speed.
The L' 1 signature is also evident in neutral temperature and zonal wind, encompassing the coastal lines of the Pacific Rim. At UT 18 (13.5 hr after the eruption), wave signatures in zonal wind remain strong, with the L' 0 mode roughly threading through the Caribbean Sea, Baltic Sea, Caspian Sea, and South Africa. The L' 1 mode is seen over the northwest coast of South America, Madagascar, and the East Siberian Sea, and aligns well with surface perturbation (Figure 1). It is noteworthy that the large wind perturbation by L' 1 mode along the northwest coast of South America is coincident with a large wind perturbation observed by the radar facility at Jicamarca Radio Observatory, and may therefore provide an interpretation of the observation.
From Figure 2 it is seen that the wind perturbation in the ionospheric E and F region altitudes can be over 100 ms −1 . Such large wind perturbations can modulate both the ionospheric dynamo and plasma transport and produce traveling ionospheric disturbances (TIDs), as shown in Figure 4. The E × B drift in the vertical direction (vertical E × B drift hereafter) displays perturbations with phase speeds of the Lamb modes, as most apparent at latitudes to the south of the epicenter in the keogram. The vertical E × B drift also has large perturbations at the magnetic conjugate latitudes in the northern hemisphere, as expected from the equipotential mapping along geomagnetic field lines. Its northward propagation is apparently strong, and symmetric with the Southern Hemisphere branch. Multiple wave components are seen more clearly in TEC perturbations (high pass filtered with period <1 hr). Large TEC perturbations up to ∼±10 TECU are seen in the first several hours after the eruption. These large perturbations are mostly seen in the northward propagation branch, when the waves traverse the equatorial ionization anomaly during daytime/dusk. Their phase speeds span a large range, up to ∼900 ms −1 -the acoustic speed in the upper thermosphere. This is consistent with the wave structure seen in Figure 2, with faster waves in the upper thermosphere. However, the very fast waves are transient and seen only within 30°-40° latitude of the epicenter. Persistent global propagation is again associated with L' 0 and L' 1 modes, with perturbation amplitude of ∼0.3 TECU (even after ∼33 hr when the L' 0 wave propagates back to the epicenter). The global propagation in TEC perturbation, as well as its alignment with the surface pressure perturbation, is also seen in the animation (https://visgallery.ucar.edu/tonga-eruption). These features, including the initial transient perturbations with large phase speeds and amplitudes, agree with observational analysis Themens et al., 2022;Zhang et al., 2022). The L' 1 mode is also identified from GNSS TEC analysis (Li et al., 2023).
Magnetic conjugacy effect is seen in TEC perturbation in Figure 4, but mainly before UT 07 hr around 20°N, with amplitude between 0.1 and 0.2 TECU (as seen from the keogram). It is weaker than the initial TEC perturbation associated with the primary waves. It appears near Hawaii shortly after the eruption (the partial ring structure to the west of Hawaii), consistent with observations (Themens et al., 2022). Since the conjugate E × B drift perturbations have similar magnitudes, the differences in the conjugate TEC perturbations should be caused by differences in transport by neutral wind perturbations (Huba et al., 2015). It is seen that as the E × B drift perturbations become smaller at later times, the conjugate TEC perturbations decrease accordingly. For example, at UT 08:55 hr (4 hr and 20 min after the eruption), the magnitude of the total E × B drift perturbations is around 10 m −1 ( Figure S3 in Supporting Information S1). This is much less than that reported in Shinbori et al. (2022). Consequently, no apparent conjugate TEC perturbation is seen at UT 08:55 ( Figure S3 in Supporting Information S1), in contrast to observational results (Lin et al., 2022;Shinbori et al., 2022). Possible causes for the discrepancy include under-estimation of perturbations to electric field and E × B at scales smaller than that resolvable by the model. Further, wave perturbations after Hunga-Tonga eruption could have triggered ionospheric bubbles Huba et al., 2022). Magnetic conjugacy effect may also manifest through bubble formation, which cannot be resolved by WACCM-X.

Summary
Even with a simplified initiation of the Hunga eruption, the high-resolution WACCM-X simulation captures salient observed features following the eruption for the Earth surface to the thermosphere and ionosphere. In particular, the global propagation of Lamb wave L' 0 and L' 1 modes are the most robust global modes, in agreement with theoretical predictions. The wave amplitudes grow exponentially with altitudes, and consequently, they cause large perturbations in the upper atmosphere, with wind perturbations reaching hundreds of ms −1 and temperatures several hundred K. The amplitude of the L' 0 is generally stronger, though L' 1 amplitudes become Figure 4. Keogram of vertical E × B drift (total) and total electron content (TEC) perturbations (period <1 hr) at 175°W (a, b). Note that the magnetic equator at 175°W is at 2°N. Vertical E × B drift and TEC perturbations (period <1 hr) at UT 05:05 hr (half hour after eruption) (c, d).
significant at higher altitudes, likely because it is an internal wave below the mesopause. The simulated Lamb L' 0 mode produces surface pressure perturbations that are comparable with barometer record (100-250 hPa), and those by L' 1 mode is an order of magnitude less. In the thermosphere, the wind perturbation caused by the two modes becomes more comparable. The phase lines of the wavefront are vertical up to the lower thermosphere and tilt outward above, consistent with previous theoretical results.
The large wind perturbations in the thermosphere and the ionospheric E and F regions modulate the ionospheric dynamo and plasma transport. Persistent TIDs with Lamb wave propagating characteristics are seen in both E × B drift and TEC globally, consistent with the neutral wind perturbations. In addition, TIDs with large amplitude and phase speeds up to 900 ms −1 appear in the first few hours after the eruption and within 30°-40° of the epicenter. E × B drift perturbations at conjugate locations are clearly seen from the simulation. However, conjugacy in TEC is evident only within the first 2 hr after the eruption, when the E × B drift perturbations are strong.
The global propagation pattern and the amplitude of the simulated surface pressure perturbation are similar to observations. In the upper atmosphere, the large zonal wind perturbations (over 150 ms −1 ) and the outward tilting of the wavefront are generally consistent with MIGHTI observations, and the propagation speed of the wavefront in the simulation is similar to the MIGHTI observations. The timing of the thermospheric wind perturbation at the northwest coast of South America due to L' 1 mode is coincident with the radar observation. Characteristics of TIDs in TEC are generally comparable with those derived from GNSS observations.

Data Availability Statement
NCAR CESM/WACCM is an open-source community model, and is available at https://doi.org/10.5065/ D67H1H0V. Model output used for this study is available through Globus (globus.org) with the shared endpoint at https://tinyurl.com/3hnwjz93. Registration for a free Globus account is required to connect through the endpoint.