Meteotsunami observed in Japan following the Hunga Tonga eruption in 2022 investigated using a one-dimensional shallow-water model

On January 15, 2022, the volcano Hunga Tonga about 8000-km away from Japan explosively erupted. Following the eruption, tsunami-like sea-level fluctuations were observed in Japan, much earlier than expected based on the oceanic long-wave propagation from Tonga. By contrast, atmospheric pressure disturbance presumably due to the eruption was also observed about 30 minutes before the sea-level change. Therefore, the observed sea-level fluctuations can be considered as meteotsunamis forced by the pressure perturbation rather than tectonically forced by the eruption, but the mechanism is not yet fully understood. This study attempts to understand the nature of this meteotsunami by using a simple one-dimensional shallow-water model. The results show that the time and amplitude of the observed sea-level changes are consistent with the simulated sea-level changes forced by the atmospheric forcing. A set of experiments with different bathymetry profiles also reveals the importance of amplification due to near-Proudman resonance over deep basins and the shoaling effect over the continental slope, while extremely deep and narrow topography such as trenches is of second-order importance.


Introduction
At a volcano called Hunga Tonga-Hunga Ha'apai located in the Tonga Islands, a large-scale volcanic eruption occurred at around 1pm on January 15, 2022 (we hereafter use Japan Standard Time as local time). In response to this eruption, the Japan Meteorological Agency (JMA) announced at around 7 pm that sea-level changes were expected in Japan, but they expected no concern about large fluctuations that could be deadly for human lives or induce societal damages. This announcement was based on the evidence that the magnitude of tsunami events had been limited in other Pacific regions closer to Tonga.
What happened, however, was different from what we had expected. Figure 1a shows the time series of observed sea levels at Amami (28°N, 129°E), Ayukawa (38°N, 141°E), and Tokyo (35°N, 139°E) (Fig. 2). These time series exhibit prominent sea-level fluctuations after 8:30-9 pm, which means that it has taken ~7.5 hours for the information of the eruption to travel about 8,000 km from Tonga to Japan. If this information had been conveyed by oceanic tsunami waves, whose phase velocity is calculated as where g is the gravitational acceleration and H is the equivalent depth, then H would be about 9,000 m. Considering the ocean depth between Japan and Tonga is typically about 5,000 m (Fig. 2), this equivalent depth is too deep. The sea-level fluctuations were observed too early to be due to tsunami waves directly excited by the eruption. In fact, later in the evening, JMA announced tsunami warnings in Japan. The observed sea-level fluctuations were earlier and larger than expected, and the mechanism is yet to be fully explained.
As a promising mechanism for understanding similar sea-level fluctuations, Press and Harkinder (1966) reported that atmospheric waves generated by the notable volcanic eruption of Krakatoa in 1883 caused large sea-level fluctuations. Their following work, Harkrider and Press (1967), numerically investigated multiple modes of air-sea coupled waves excited by an air pulse produced by Krakatoa. They showed that, even if the oceanic equivalent depth H is fixed, atmospheric free waves with phase velocities close to transfer energy to the ocean efficiently. In fact, this kind of sea-level fluctuations, which is often referred to as meteotsunamis or Abiki, has been extensively studied in contexts unrelated to volcanic eruptions (e.g., Hibiya and Kajiura 1982;Fukuzawa and Hibiya 2020;Kubota et al. 2021). The main mechanism of meteotsunami is called the Proudman resonance (Proudman 1929), a physical process where pressure perturbations of atmospheric waves amplify oceanic waves that have close phase speeds.
In this study, we hypothesize that the sea-level fluctuations that followed the recent volcanic eruption in Tonga were caused by atmospheric waves, and present theoretical considerations on whether this hypothesis is consistent with observed evidence. Because the numerical simulations in this study are idealized, we cannot immediately conclude that the presented physical process was of first-order importance in the real world. Nevertheless, we show that a plausible explanation can be given if atmospheric waves worked as a main forcing. In particular, we pose the following scientific questions:  Where on Earth did the Proudman resonance contribute to the sea-level fluctuations most efficiently?
 Could the existence of deep regions such as the ocean trench have been important?
 How could the coastal shallow region and the seafloor topography enhance and/or sustain sealevel fluctuations?
This paper is organized as follows. In section 2, data and methods are presented. In section 3, results obtained from idealized models are presented. In section 4, we will present conclusions.

Observations
The sea-level data is from the Realtime Tidal Data archived by Hydrographic and Oceanographic department of Japan Coast Guard available at https://www1.kaiho.mlit.go.jp/TIDE/gauge/ (Accessed: 11:13 am, January 20). The time span used in this study is from 6 pm, January 15, to the following noon, and the temporal resolution is 5 minutes. The observed air pressure data is from the Soratena meteorological stations owned by Weathernews Inc., which is distributed for free for academic use. The detailed description of the data is presented at https://global.weathernews.com/news/16551/. The time span used in this study is from 7 pm, January 15, to the following midnight, and the temporal resolution is 1 minute.

Simple shallow-water model
To investigate the oceanic response to atmospheric pressure waves, we have conducted numerical experiments with a simple onedimensional model. A set of fundamental equations of the model is the linearized, non-rotating shallow-water equations: where η(t,x) and u(t,x) are the sea-level change from its mean state and the current velocity, respectively, at any time t and location x on the Cartesian coordinates with 1-km resolution. ρ and g are the density of seawater (1027 kg m -3 ) and gravitational acceleration (9.81 m s -2 ), respectively. We have integrated Eq. (1) by the leap-frog scheme from the rest state over the 10,000-km long domain (-2,000 km ≤ x ≤ 8,000 km) as illustrated in Fig. 3. The boundary condition is u = 0 at x = -2,000 and 8,000 km, although it may overemphasize the amplitude of sea-level fluctuations at the boundaries. We consider the boundary at x = 8,000 km as the coast of Japan and focus on only the domain 0≤x≤8,000 km. Around the opposite boundary at x = -2,000 km, strong linear damping is applied to attenuate reflected waves.
The atmospheric pressure forcing P is imposed as an external forcing and assumed to be a period of sinusoidal function as where the amplitude P0, wavelength L and propagation speed V are prescribed as 2 hPa, 500 km and 310 m s -1 , respectively, in accordance with observations in Japan (see section 3.1). At t=0, the pressure disturbance center is located 8000-km away from the coast (x = 0), which is approximately equal to the distance between Tonga and Japan. In reality, both atmospheric and oceanic waves propagate on a sphere, rather than on a plane with Cartesian coordinates. Nevertheless, we consider that the difference in coordinate system does not seriously affect the results, since the effect of ocean topography relatively close to Japan is important as we see later.
We have designed six experiments with different bathymetry profiles as indicated in Figs. 4 and 5 (black lines): realistic bathymetry (RB), idealized bathymetry (IB), 6,000-m flat floor (FF60), 4,500-m flat floor (FF45), no trench (NoT), and no shelf (NoS). The inhomogeneity in affects wave behaviors through the volume flux convergence in Eq. (1a). In the RB experiment, H(x) has been taken from the real bathymetry along the great circle between Tonga and Japan ( Fig. 2) based on the General Bathymetric Chart of the Oceans (GEBCO)'s 2021 gridded data (https://doi.org/10.5285/c6612cbe-50b3-0cff-e053-6c86abc09f8f). The bathymetry is simplified in the IB experiment to represent deepening ocean basin, trench, and steep continental slope. The two flat floor experiments, FF60 and FF45, have flat seafloors of 6,000 and 4,500 m, respectively. The NoT and NoS experiments exclude domains deeper than 6,000 m around the trench and shallower than 2,000 m near the coast, respectively, from the idealized bathymetry.

Observed evidence
The sea-level fluctuations started around 8:30-9 pm, and lasted more than half a day. Figure 1a shows the sea-level fluctuations at three representative sites (Fig. 2). At Amami, the sea-level fluctuations started weakly, and a few hours later, exhibited one of the largest fluctuations in Japan. Sea-level fluctuations are observed even in the inland sea such as Tokyo, though they are more moderate than those observed at stations closer to the open ocean. As shown in the Introduction, such sea-level fluctuations are too early to be tectonically-induced tsunamis.
The pressure fluctuations with amplitude of ~2hPa also arrived in Japan around 8:30 pm, but started about 30 minutes earlier than the sea-level fluctuations. Figure 1b shows air pressure fluctuations at three stations near the locations shown in Fig. 1a. The eruption occurred at around 1 pm, so the pressure fluctuations traveled for less than 7.5 hours over 8,000 km from Tonga to Japan. Hence, the phase velocity of the atmospheric waves is estimated to be about 310 m s -1 , which virtually represents the speed of sound, consistent with satellite-based estimation by Otsuka (2022). In addition, since the pressure fluctuations took about 30 minutes to rise and fall, the wavelength is estimated to be about 500 km. In the RB experiment, many waves propagating in both directions are extensively generated (Fig. 4a). Only the rightward first wave propagates at the speed of the pressure disturbance V throughout most of the basin, while all other waves at the phase speed of long wave c = (free waves). The rightward free waves emerge not only at t = 0 but also on the way the first wave propagates.

Simulated sea-level fluctuations
If the realistic topography is simplified as in the IB experiment, only a pair of the first wave and the following free wave propagating rightward from x = 0 is generated (Fig. 4b). Comparing the sealevel anomalies at the coast (x = 8,000 km) in both experiments (Figs. 6a and 6b), the IB experiment reproduces well the first wave seen in the RB experiment. The first wave arrives about 30 minutes later than the pressure wave and about 4 hours earlier than the free wave propagating throughout 8,000 km, both of which are consistent with the observations (Fig. 1). The first wave is synchronized with the pressure wave while travelling over the basin. When it reaches the continental slope, however, the effect of volume flux convergence due to the steep slope becomes significantly large, which acts to amplify the sea-level displacement, rather than propagating it forward [Eq. (1a)]. The wave thus slows down and begins to propagate as a free wave separated from the pressure wave (Fig. 5b), resulting in the late arrival behind the pressure wave at the coast. The separation of the first wave from atmospheric forcing over shoaling slopes was also demonstrated in numerical experiments by Vilibić (2008) and Chen and Niu (2018).
In contrast to the similarity in the first wave between the RB and IB experiments, the latter shows much weaker subsequent sea-level oscillations with partially reflected waves over the continental slope (Figs. 5b and 6b). This lack of following waves suggests that the oscillatory behavior is a manifestation of the effect of relatively small-scale topography. In the RB experiment, when the first wave rides up on the shoals along the way from Tonga to Japan and cannot catch up the atmospheric forcing, the first wave is partially converted into the subsequent free waves repeatedly (Figs. 4a and 5a). However, since such shoals represented in the RB experiment mostly reflects isolated seamounts rather than ridges and waves can propagate around them in the real world, the effect may be exaggerated in this one-dimensional model. Nevertheless, the effect of relatively large-scale topography, such as island arcs and ridges, may have some impacts.
If the ocean basin is flattened, the first wave amplitude is nearly constant while propagating over the basin after separating from the free wave near x = 0 (Figs. 4c and 4d). According to previous studies (Proudman 1929;Vilibić 2008;Williams et al. 2021), the sea-level displacement η over the uniform depth can be written as a superposition of a forced wave, a rightward free wave and a leftward free wave as , ! 1 ,-. / 2 1 ,-% / 2 1 % ,-. 3 η0 = -P/ρg is the equivalent static sea-level change in response to a steady pressure anomaly (~2 cm for 2 hPa), and Fr = V/c is the Froude number. At H = 6,000 m, the forced-wave amplitude, η0/(1-Fr 2 ), is thus estimated to be ~3.2 cm, which is consistent with to that of the simulated first wave (Fig. 4c). The amplitude is larger than η0, indicating the near-Proudman resonance. In the FF45 case, the forced wave is weaker with the amplitude of ~1.7 cm (Fig. 4d), again consistent with the theoretical estimation. Interestingly, despite the varying depth in the IB experiment, the first wave amplitude seems to be determined by the forced wave in Eq. (3) given the local water depth. The first wave gradually grows while propagating across the ocean basin (Fig. 4b), and the amplitude at x = 7,000 km is almost equivalent to that in the FF60 experiment (Fig. 4c). Unlike the flat floor cases, the first wave in the Fig. 5. Same as in Fig. 4 but for the domain near the coast (7,500 km ≤ x ≤ 8,000 km). IB experiment over the basin is not purely forced wave, since it is not completely in phase with the pressure forcing (Fig. S1). The net positive work done by the atmosphere due to the slight leading of the first wave to the forcing can be supporting, but not sufficient, evidence of the amplification. We speculate that the slope in the IB experiment is gentle enough for the amplification to work slowly, so Eq. (3) can locally determine the first wave amplitude.
On the basis of Eq.
(3), one may expect that the trench can act to amplify the wave significantly. However, Figs. 6b and 6e reveal that the trench has little effect on the amplitude. This insensitivity to the trench can be attributed to the narrowness of the trench relative to the wavelength of pressure forcing.
In all of the experiments, the first wave undergoes significant amplification near the coast (Figs. 5a-e). This amplification is mainly due to the shoaling effect on long waves: deceleration of the phase speed shortens their wavelength and grows their height to conserve the energy flux. This effect makes the wave amplitude inversely proportional to the fourth root of the depth, which is referred to as the Green's law, together with the contribution from narrowing intervals of wave rays (Lamb 1932). This law is widely used to analyze not only meteotsunamis but also tectonically-induced tsunamis (e.g., Sandambata et al. 2018). Shallowing from the 6,000-m to 20-m depth results in the amplification by a factor of ∜(6,000 m/20 m) = ~4.2. We have decomposed the sea-level fluctuations in the NoT experiment into rightward and leftward components (Fig. 7). Figure 7a and 7b indicate the partial reflection over the continental slope (7,800 km < x < 7,900 km). The amplitude of the leftward reflected wave is ~1/6 of that of the rightward first wave. Green's law applied to the simulated first wave in rightward component at x = 7,800 km (Fig 7c, green line) reasonably explains the simulated amplitude (Fig 7c, orange line). However, on the coastal side of the continental slope (7,950 km < x), Green's law slightly overestimates it. This overestimation is presumably due to the partial reflection as pointed out by Hibiya and Kajiura (1982). Indeed, the amplification of waves becomes weaker in the NoS experiment (Figs. 5f and 6f).

Conclusions
After the Hunga-Tonga eruption in 2022, tsunami-like sea-level fluctuations were observed in Japan following the atmospheric pressure disturbance presumably generated by the eruption. To explain the nature of the observed sea-level fluctuations as simply as possible, we have performed six numerical experiments with different depth profiles by using a one-dimensional shallow-water model. Based on the simulations, we have reached the following conclusions about this tsunami-like waves:  The tsunami-like waves were forced by atmospheric pressure waves propagating faster than oceanic long waves, resulting in the arrival of the waves 2 to 3 hours typically earlier than predicted.
 The first wave propagated over the ocean basin synchronized with the atmospheric forcing and was converted to the free waves over the continental slope. Thus, the arrival of the first wave at the coast is delayed by a few tens of minutes behind the pressure wave.
 The first wave is amplified by the near-Proudman resonance over the 6,000-m class basin and shoaling effect (Green's law) near the coast.
 The trench is too narrow for the Proudman resonance to be effective, so it has little effect on the first wave amplitude.
The detailed mechanisms involved in the forced-to-free waves conversion and partial reflection over varying bathymetry should be elucidated in future studies. We also note that the model we used oversimplifies the real world and cannot explain all aspects of the observed phenomenon, although we consider that it captures the essence. Since the one-dimensional model does not have coastal topography in the direction perpendicular to wave propagation, we cannot evaluate amplification mechanisms though energy flux convergence into shoals, harbor resonance and interaction with edge waves, which are important in ordinary meteotsunamis (e.g., Hibiya and Kajiura 1982;Monserrat et al. 2006;Fukuzawa and Hibiya 2020). Moreover, oceanic waves directly forced by the volcano eruption and atmospheric waves slower than the speed of sound may have also arrived several hours later, as pointed out by Harkrider and Press (1967) for the Krakatoa eruption. Numerical simulations with more realistic configurations will shed light on these issues in future studies.