Multi-scale Simulation of Subsequent Tsunami Waves in Japan Excited by Air Pressure Waves Due to the 2022 Tonga Volcanic Eruption

The 2022 Hunga Tonga-Hunga Ha’apai eruption generated tsunamis that propagated across the Pacific Ocean. Along the coast of Japan, nearshore amplification led to amplitudes of nearly 1 m at some locations, with varying peak tsunami occurrence times. The leading tsunami wave can generally be reproduced by Lamb waves, which are a type of air-pressure wave generated by an eruption. However, subsequent tsunamis that occurred several hours after the leading wave tended to be larger for unknown reasons. This study performs multi-scale numerical simulations to investigate subsequent tsunami waves in the vicinity of Japan induced by air pressure waves caused by the eruption. The atmospheric pressure field was created using a dispersion relation of atmospheric gravity wave and tuned by physical parameters based on observational records. The tsunami simulations used the adaptive mesh refinement method, incorporating detailed bathymetry and topography to solve the tsunami at various spatial scales. The simulations effectively reproduced the tsunami waveforms observed at numerous coastal locations, and results indicate that the factors contributing to the maximum tsunami amplitude differ by region. In particular, bay resonance plays a major role in determining the maximum amplitude at many sites along the east coast of Japan. However, large tsunami amplification at some west coast locations was not replicated, probably because it was caused by amplification during oceanic wave propagation rather than meteorological factors. These findings enhance our understanding of meteotsunami complexity and help distinguish tsunami amplification factors.


Introduction
On January 15, 2022, a massive eruption occurred at the Hunga Tonga-Hunga Ha'apai volcano (hereafter referred to as 'the volcano') in Tonga, an island nation in the South Pacific.While the eruption had been ongoing since December 2021, the most significant event occurred at approximately 4:14 UTC on this date (Purkis et al., 2023).This eruption was one of the most intense volcanic events in recent history (Terry et al., 2022), releasing an enormous amount of energy and ejecting a substantial volume of volcanic ash and gas into the atmosphere.It generated a series of tsunami waves that propagated across the Pacific Ocean, impacted numerous coastal regions and caused considerable damage (Borrero et al., 2023).
Tsunami waves after the major eruption have been observed in various areas along the Pacific coast, including New Zealand (Gusman et al., 2022), Australia (Munaibari et al., 2023), Mexico (Ortiz-Huerta & Ortiz, 2022;Ramı ´rez-Herrera et al., 2022), Chile (Carvajal et al., 2022) and Japan (Imamura et al., 2022;Wang et al., 2023).The tsunami waves traveled over the Pacific Ocean, caused widespread damage and exhibited very complex behavior worldwide.The far-reaching impacts of the tsunami underscore the vulnerability of coastal communities to such natural disasters.
In Japan, tsunamis reaching amplitudes of 0.5-1.0m were observed over a wide area.No serious tsunami inundation damage was reported; however, an amplified tsunami caused ships to capsize in harbors (Imamura et al., 2022).The maximum amplitude was delayed by several hours or more after the arrival of the leading tsunami.For example, the maximum amplitudes were recorded over 6 h after the leading tsunami at some gauges such as Kushiro, Nemuro and Urakawa.In addition, the leading tsunami reached the coast of Japan several hours earlier than expected.The sudden and unexpected tsunami caught many residents off guard, highlighting the need for improved early warning systems and coastal resilience measures.Wang et al. (2023) analyzed records of this tsunami event from high-frequency radars, offshore bottom pressure gauges and tidal gauges on the east coast of Japan, respectively, and attempted to develop a prompt forecast network of tsunami waveforms using data assimilation.Their results showed that it is difficult to predict tsunamis using a traditional method because the mechanism of wave generation is different from that of general tsunamis.
The eruption of the volcano generated a complex set of pressure, Lamb and air gravity waves (Amores et al., 2022;Yuen et al., 2022).These waves, which traveled at different speeds, contributed to the generation and propagation of the tsunami waves observed in the Pacific.Atmospheric pressure generated by volcanic eruptions may excite the sea surface, resulting in a type of tsunami that occurs less frequently than other natural disasters.The last confirmed event of this type was the 1883 eruption of Krakatau in Indonesia (e.g., Nomanbhoy & Satake, 1995;Paris et al., 2014).These tsunami events are classified as meteotsunami (or meteorological tsunami; e.g., Monserrat et al., 2006;Rabinovich & Monserrat, 1996;Rabinovich 2020).The interaction between air pressure changes and the ocean surface can lead to the formation of waves that can be amplified and travel large distances, resulting in significant impacts on coastal regions.When the celerity of long waves in the ocean matches the speed of the air pressure waves, the ocean wave can be continuously excited and largely amplified.This physical phenomenon is well known as Proudman resonance (Hibiya & Kajiura, 1982;Vilibic ´et al., 2008).
The Tongan tsunami event excited by Lamb waves of preceding pressure waves has been the subject of several studies on global (e.g., Kubota et al., 2022;Lynett et al., 2022) and regional scales (e.g., Gusman et al., 2022;Tanioka et al., 2022).These studies were fairly successful in reproducing tsunami waveforms synchronized to the arrival time of the Lamb wave and showed that the Lamb wave was dominant in the leading tsunami behavior.The relationship between pressure-induced ocean waves and sea bathymetry was also investigated (Ho et al., 2023;Liu & Higuera, 2022).The occurrence of Proudman resonance has also been suggested.However, there is currently no simulation capable of reproducing a large tsunami amplification in the range of 1 m on the Japanese coast using realistic bathymetry and pressure wave data.Suzuki et al. (2023) performed a tsunami simulation around Japan using realistic bathymetry and pressure waves that were numerically solved.Their study reproduced resonance in the open ocean but did not lead to significant amplification along the coast.The difficulty in simulating large amplified tsunamis lies in the diversity of their amplification sources.Proudman resonance is most likely to occur in deep, flat, open ocean regions, where ocean wave speeds and atmospheric pressure velocities correspond to each other.The horizontal scale of this phenomenon is on the order of a few kilometers.However, to reproduce resonance in bays, it is necessary to capture tsunami behavior at the smaller scale of each bay.Tsunami simulations that include phenomena at various spatial scales are thus required to understand in detail the factors that cause tsunami amplification.
Complicating this tsunami phenomenon is the diversity of tsunami wave sources.During this tsunami event, tsunamis other than meteotsunami were also observed near the volcano (Borrero et al., 2023).In the islands of Oceania near the volcano, the topographic change caused by the caldera collapse may have contributed to the tsunami run-up.In addition, although the Lamb wave and air gravity waves have different velocities, they overlap each other.Thus, the tsunami behavior near volcanoes was complex and simultaneously influenced by various factors.On the other hand, in Japan, the effect of tsunami caused by the caldera collapse was almost negligible, and the Lamb wave and gravity waves were separately observed.Therefore, elucidating the mechanism of tsunami amplification in Japan can contribute to understanding the characteristics of pressure wave of this eruption and global tsunami phenomenon.
This study aims to provide a comprehensive multi-scale numerical simulation of tsunami waves in Japan excited by air pressure waves caused by the 2022 Tonga volcanic eruption.Particular emphasis was placed on the waves after the initial arrival of the tsunami.The atmospheric pressure was modeled by considering the dispersion relation of air gravity waves.The representation of pressure waves passing around Japan was tuned using pressure observation records.High-resolution computationally demanding simulations over a wide coastal area were conducted using an efficient simulation method.Note that we did not consider tsunamis triggered by caldera collapse, as a preliminary study indicated that such tsunamis did not affect the Japanese coast, and previous studies (Heidarzadeh et al., 2022;Terry et al., 2022) have also shown that their effect is negligible in the North Pacific region.
The remainder of this article is structured as follows.Section 2 describes the curation and processing of observational data, pressure wave modeling and methods for creating pressure data and provides an overview of numerical models for tsunami simulations, simulation conditions and bathymetry data.Section 3 describes the results of the meteotsunami simulation and discusses their comparability to observational records.Section 4 uses the simulation results to explain the amplification sources along the coast of Japan.Finally, conclusions from the study are presented in the last section.

Observation Data and Processing
In this study, offshore and coastal waveform observation records were used to validate the results of the tsunami simulations.Data were sourced from offshore ocean buoys forming part of the Deep-ocean Assessment and Reporting of Tsunamis (DART) project.DART is a real-time tsunami detection system managed by the United States National Oceanic and Atmospheric Administration (NOAA) and provides critical information to protect coastal communities from the impacts of tsunamis.The DART program has a long history of successful operations (e.g., Mungov et al., 2013;Rabinovich & Eble ´, 2015), and its buoys deployed throughout the world's oceans provide a wealth of information for researchers.During the 2022 Tongan tsunami, five DART buoys (21,418, 21,420, 52,401, 52,402 and 52,404) were operational in the vicinity of Japan.To extract tsunami signals, periodic components [ 100 min were filtered out.
The Japan Meteorological Agency (JMA) tidal or tsunami gauges were used for coastal observation records.The JMA operates an extensive network of tidal gauges distributed throughout Japan, which are essential for monitoring sea level changes and predicting coastal hazards, such as storm surges and tsunamis.This study used data from 40 of these gauges on the Pacific Ocean side, as shown in Fig. 1.The gauges where the maximum amplitudes were below a threshold of 5 cm or not significant were omitted.As shown in the figure, the gauges for evaluation form a dense series along the Japanese coastal areas.They were generally numbered from west to east, except for gauge number 30, which was located on an isolated island (Chichijima).
Table 1 provides detailed information on the 40 JMA gauges, showing coordinates, peak amplitude during the event and the time of its occurrence (time elapsed after the main explosion at UTC 04:14 on January 15, 2022).This relative time measure is used in the following.The maximum tsunami amplitude for each gauge was extracted from the observed waveform by applying a bandpass filter with corner frequencies of 1/300-1/6000 Hz (5-100 min) to remove unrelated signals.Note that the maximum tsunami amplitudes used here do not necessarily coincide with those published by JMA because periodic components \ 5 min were filtered out.The reason for filtering out the shorter side of the period is the horizontal resolution of the tsunami simulation.Wave components with a period of a few minutes cannot be reproduced unless the resolution is quite high.The tsunami amplitudes and occurrence times differed between gauges by several hours or more.This occurred even between nearby sites, indicating the complexity of this event and the diversity of the factors causing major tsunamis.
In addition, JMA atmospheric pressure observation data were used to tune the simulated spatiotemporal distributions of atmospheric pressure.Time series data from four gauges (Naze, Nishinoomote, Hachijojima and Kamaishi) were used to target only the eruption-induced pressure changes.These gauges face the Pacific Ocean, are located on small islands and are less affected by topography.Their eruption-induced signals were, therefore, relatively clear among the atmospheric pressure gauges near Japan.Similar to the ocean wave observations, the long-period components were filtered out for comparison purposes.
Spectral and wavelet analyses were used to analyze the observed waveforms, perform preprocessing and compare the characteristics of the generated atmospheric pressure as shown in the following sections.These analyses are important tools for observing wave characteristics in the frequency domain and have been used in numerous tsunami cases.In particular, wavelet analysis is an essential method for representing the time variation of the signal magnitude at each frequency and is most useful for events such as the 2022 Tonga tsunami, where the characteristics of the air pressure as a wave source change with time (e.g., Wang et al., 2022).

Pressure Wave Generation
The numerical tsunami model used a realistic spatiotemporal distribution of pressure sourced from observational records.The pressure waves caused by a volcanic eruption can be broadly classified into two types: Lamb waves and acoustic gravity waves (e.g., Kulichkov et al., 2022;Wright et al., 2022).Lamb waves propagate horizontally at the speed of sound across the Earth's surface exhibiting no vertical oscillation in their velocity.Acoustic gravity waves propagate vertically and horizontally and result from buoyancy forces within the atmosphere, causing  periodic oscillations in air pressure.The 2022 Tonga volcanic eruption produced both types of waves, and they have been detected and analyzed in many studies striving to understand the impact of the eruption on the atmosphere.The maximum pressure anomaly for Lamb waves near Japan has been reported as approximately 2 hPa (Imamura et al., 2022;Nishikawa et al., 2022).As these waves propagate at the speed of sound and have no dispersion, the observed waveforms can be reproduced relatively easily by simple amplitude tuning.This study modeled Lamb waves in the same manner as Gusman et al. (2022) and Kubota et al. (2022), assuming that the pressure amplitude of peak P p ðrÞ and trough P t ðrÞ was proportional to r À1=2 , where r is the distance from the volcano (km).Our model also employed separate amplitude parameters for the peak and trough of the wave; P p ðrÞ ¼ 169= ffiffi r p (hPa) and P t ðrÞ ¼ À109= ffiffi r p (hPa).These amplitude parameters were adjusted to approximately 2 hPa for the pressure peaks and 1 hPa for the troughs.The speed of sound was assumed to be 310 m/s.Acoustic gravity waves have dispersive characteristics and require more complex modeling.The dispersion effect and smaller amplitude also make it difficult to determine the intensity of each individual wave at each wavelength from the collection of observed waveforms.We did not perform numerical simulations of pressure waves but simply superimposed parametric pressure waves based on the dispersion relation of the internal gravity waves.The assumptions behind this relation are briefly described below.The modeling is identical to that in Wright et al. (2022).
The basic linear theory can often be used to describe atmospheric gravity waves in terms of minor deviations from a vertically varying and stably stratified background state.The driving force behind Gray-shaded areas represent the range of acoustic and atmospheric gravity waves, respectively the gravity waves is the buoyancy generated by the adiabatic displacement of air parcels.First, the horizontal-vertical (x-z) plane is considered (as in Nishikawa et al., 2022).The Coriolis force is neglected, and the buoyancy frequency N is considered constant.Although the buoyancy frequency N and speed of sound c s generally depend on altitude above sea level, they are considered constant here to simplify the solution form.This assumption is useful for understanding the characteristics of gravity waves in a dispersion relation.Let the solution of the traveling gravity waves be in the form of where w is the vertical velocity and ŵ is the amplitude of the vertical velocity fluctuations.vector k ¼ ðk; mÞ is the wavenumber of the disturbance and x is the frequency.The characteristics of the solution are conveniently illustrated in the case of an isothermal atmosphere, and the velocity decreases exponentially with height: These assumptions and the linearized forms of the fundamental equations [e.g., Gill (1982) and Fritts and Alexander (2003)] yield the dispersion relation This equation is fourth order in x and represents both acoustic and gravity waves.Solving x 2 yields the following equation: The positive root represents acoustic waves, and the negative root represents gravity waves.
Figure 2 illustrates Eq. ( 4), by showing x as a function of k.The Lamb wave characteristics are also shown in a separate relation x 2 ¼ k 2 c 2 s using a hydrostatic approximation.For clarity, the vertical and horizontal scales are non-dimensionalized on the representative scales N and H s , where H s is Earth's scale height.N and H s in this figure are assumed to have typical values of 1:18 Â 10 À2 s À1 and 8:4 Â 10 3 m, respectively.This x-k plane is divided into two regions: x [ kc s represents sound waves, and x\kc s represents gravity waves.In the gravity wave region, the frequency asymptotically approaches N as k increases and reaches its maximum when m ¼ 0.
By fine-tuning different N and c s parameters and changing the number of wave superpositions and the amplitude of each wave, various spatiotemporal pressure patterns can be generated using the same equations.Following the dispersion relations, several patterns of composite waves were generated by superimposing waves of various wavelengths to reproduce pressure waveforms consistent with observed records in Japan.Prior to tuning the parameters, from a comparison between the time series of tidal gauges and atmospheric pressure gauges, ocean waves are amplified at many locations immediately after the passage of the atmospheric gravity waves.Meanwhile, it is evident from Fig. 2 and Eq. ( 4) that the atmospheric acoustic wave always travels faster than the atmospheric gravity waves.For these reasons, acoustic waves were not considered in the study.For simplicity, the amplitude of gravity waves was determined by tuning the parameter to the À1=2 power of the distance from the volcano for each wavelength, as in the case of Lamb waves.
Four patterns of spatiotemporal pressure fields were generated.Table 2 shows an outline of the generated pressure conditions.One was composed only of the Lamb wave (referred to as Pressure L), and the other three consisted of the Lamb and gravity waves (referred to as Pressures A, B and C).The components of the Lamb wave are the same for all pressure conditions.For gravity waves, Pressure A, B and C have different intentions.Obviously, not all possible scenarios could be considered, but we set up fewer patterns so that a broad effect could be tested.The details of the conditions and their rationale are described below.First, we recognized from the spectral analysis of atmospheric pressure observations that the pressure waves in a wide range of frequencies arrived at Japan.However, the source and spatio-temporal characteristics of atmospheric pressure in the whole frequency domain are poorly understood.Besides, the contribution of each frequency band to the tsunami amplification can be identified by creating pressure distributions with a specific frequency band separately.Based on this idea, several different patterns were deliberately created the following patterns were considered: Pressure A should include as wide a frequency component as possible in a short time window and have a similar time series to the observed ones at JMA atmospheric pressure gauges; Pressures B and C should test whether pressure wave trains with a small amplitude (0.5 hPa) can cause resonance at bays, assuming that the pressure waves continue to excite the ocean surface; Pressures B and C should have different frequency bands.
Figure 3 shows the horizontal distribution of generated pressure fields on the sea surface.These pressure fields are shown for the same time point in the four panels, where green circles represent JMA pressure observation gauges.Patterns of Pressures A-C have different characteristics.At Pressure A, waves of various wavelengths radiated simultaneously from the volcano, each with a different amplitude.When these waves reached Japan, they exhibited a more complex time series because of the different phase velocities of each wave.This pattern is a sea level pressure wave pattern generated to match the waveforms at the JMA pressure gauges at Naze and Hachijojima.The number of waves, wavelength and wave amplitude were adapted to agree with the time series at these gauges.
In contrast, Pressures B and C have long wave groups, as shown in the lower panels of Fig. 3.This was intended to reproduce long pressure oscillations with small amplitudes, assuming that waves were continuously radiated.The spatial pattern of the wave group in these patterns was more monotonous than that at Pressure A, implying that the wavelength variations of the gravity waves were less diverse.Wavelengths of wave groups differed between Pressures B and C, intending to produce different pressure wave velocities.
Figure 4 shows the air pressure waveforms and wavelet analyses for Pressures A, B and C at the Naze gauge.The elapsed time was referenced to the main eruption time, which was also used in the following figures.Pressure L was omitted because the leading Lamb wave section did not differ between Pressures A to C (top panel), which was generally consistent with the observations, but those waveforms differed after the Lamb wave passed.The amplitude following the Lamb wave was small (\ 0.5 hPa) and exhibited different frequency characteristics.Pressure A contained a wide frequency range of 10-to 50-min components several hours after the subsequent incident waves (9-12 h).The relatively high intensities of these gravity wave frequency bands correspond to many waves with horizontal velocities of 100-250 m/ s; this broadly includes the velocity ranges that are estimated to cause Proudman resonance (Kubota et al., 2022;Suzuki et al., 2023;Tanioka et al., 2022).Pressures B and C each cover a different frequency range for an extended duration, with B being centered mainly on 20-40 min and C on 10-30 min.
The pressure dataset developed in this study is publicly available.Comparisons of the created pressure waveforms with observations at Nishinoomote, Hachijojima and Kamaishi are provided in the Supplementary Information.Following the creation of the atmospheric pressure waveforms, four patterns of water surface waves were simulated based on the spatiotemporal distribution of sea surface pressure, including one in which only Lamb waves were considered.These are described in the next section.

Numerical Tsunami Simulation
Numerical simulations of the meteotsunami event were conducted under the pressure wave conditions generated above.Figure 5a shows the entire computational domain and the horizontal distribution of Figure 5b presents an enlarged map of Japan.On the way from Tonga to the coast of Japan, pressure and ocean waves pass through several ocean basins and ridges.The Japan Trench and Izu-Ogasawara Trench, which run parallel to the east coast of Japan, are bordered by various ocean basins and ridges.This suggests that different types of ocean waves were excited offshore along the eastern and western coasts.The depth of the offshore basins ranged from approximately 4000 to 6000 m, corresponding to a possible surface wave celerity of 200-240 m/s.It is unclear from which region pressure-induced waves started that contributed to the maximum tsunami to the coasts of Japan.Therefore, we solved the ocean wave propagation in the entire path from the volcano to Japan.
The bathymetry data for this domain were obtained from the General Bathymetric Chart of Oceans (GEBCO) Dataset 2022.The GEBCO2022 data have a horizontal resolution of 15 arcseconds, which is generally sufficient for offshore tsunami simulations.However, as this is considered insufficient for solving phenomena that require small spatial scales, such as resonance in bays, we also obtained high-resolution bathymetry data.Bathymetry data provided by the Central Disaster Management Council (CDMC), Cabinet Office of Japan, and data from the Japan Hydrographic Association were used to solve with a resolution of 5 arcseconds.There were 25 bathymetric data domains at this resolution, capturing the detailed topography and bathymetry around all 40 coastal gauges, as shown in Fig. 1.The CDMC bathymetry data are properly converted from the Cartesian coordinates (originally provided) to the geographical coordinates using the Gauss-Kru ¨ger Projection (Kawase, 2013).The bathymetry of these domains is shown in the Supplemental Information.Note that, as described below, tsunami behavior was not necessarily resolved at this high resolution.
The numerical package GeoClaw (Mandli & Dawson, 2014), which solves nonlinear shallow water equations, was used for the tsunami simulations.Geoclaw employs a finite-volume method with adaptive mesh refinement (AMR) and has been used in previous simulations of this tsunami event (e.g., Omira et al., 2022;Yuen et al., 2022) and other meteotsunamis (Kim & Omira, 2021).In simulations employing the commonly used nested grid system, the resolution and integration time intervals were fixed for each computational domain.This type of system would usually incur a high computational cost for tsunami simulations from the Tongan volcano to Japanese coastal areas because of the intervening distance of [ 7000 km.However, the use of AMR and variable resolutions as well as integration time intervals make such simulations feasible without exceptionally high-performance computers.The original GeoClaw package supports the simulation of ocean surface waves induced by pressure and wind fields by assuming a typhoon/hurricane model that can be expressed using simple parameters.In the present study, the package was extended to solve for water waves with arbitrary pressure fields.
In the AMR method, a finite number of resolution levels are set at the beginning of the computation, and the resolution of each level is determined in advance.During computation, area resolutions were adapted as required.Five resolution levels were used to optimize the computational cost: 12 arcminutes, 4 arcminutes, 1 arcminute, 15 arcseconds and 5 arcseconds.GEBCO bathymetry data were used for the fourth level (resolutions up to 15 arcseconds), while the 5 arcsecond level was based on higher resolution datasets.These resolution levels varied during each simulation, with sea surface height as the reference value.The highest resolution was applied to the surrounding computational domain when the sea surface height exceeded 2 cm.This criterion enabled the simulation to capture the tsunami wavefront.Even if a tsunami is\2 cm, waves in the vicinity can be resolved at the second highest resolution level.Because the horizontal scale of tsunami behavior in b Figure 4 Pressure waveforms and wavelet analysis at Naze pressure gauge, located on the Amami Islands at the southwestern edge of Japan.The top panel shows the comparison of generated with observed pressure waves.The horizontal axis represents the elapsed time after the main eruption.The other four panels show the wavelet (ft) diagrams for the observed data and Pressure A, B and C, respectively such waters (e.g., DART buoys) can be several hundred meters, the waves do not necessarily have to be resolved at high resolution.
This study performed tsunami simulations under four different pressure conditions (Fig. 3).The effect of sea level excitation due to spatiotemporal changes in air pressure caused by the explosion was considered.Computations were performed from the time of the start of the explosion to 16 h later, which covered the tsunami peaks at all 40 coastal gauges.A spherical coordinate system was adopted to solve tsunami behavior over a wide area.The simulations also accounted for Manning's bottom friction and, in contrast to atmospheric waves, for Coriolis forcing.Wind forcing and variations in atmospheric pressure, other than eruptive factors, were not considered.

Tsunami Offshore Japan
We first confirmed that tsunami waves with appropriate amplitudes have been recorded in the open ocean near Japan.Figure 6 shows the simulated and observed waveforms for the five DART buoys.In deep oceans, surface waves excited by atmospheric pressure scatter quickly, which results in small surface amplitudes.Under Lamb-wave-only conditions, the oscillations were almost inactive after the arrival of the initial pressure wave, and the subsequent waves were not reproduced.However, conditions in which both Lamb and air gravity waves were considered reproduced the subsequent highfrequency component with a maximum amplitude error of 1 cm.Components of much higher frequency were not the focus of this study because they are  difficult to reproduce under simplified pressure conditions.Overall, the observed waveforms from buoys closer to Japan were well reproduced because the atmospheric pressure conditions were tuned to these coast-proximate observations.It was noted that at three DART buoys of 21,418, 52,401 and 52,402, located in the open ocean southeast/northeast of Japan, the initial pulses arrived earlier than observed.This can be attributed to the effect of contrary background winds (Kulichkov et al., 2022), which were not considered in the generated pressure conditions.Therefore, a phase shift (6 or 12 min) was incorporated in the following comparisons to accommodate this phase difference on the eastern coast.We found that the incident water wave in the simulation had sufficient amplitudes at distances of tens or hundreds of kilometers offshore from the Japanese coast.The amplitudes of the subsequent waves for Pressures B and C were almost identical to those for Pressure A because of identical pressure amplitudes.In the following sections, the simulated tsunami waveforms are compared with those recorded at the JMA gauges along the coast of Japan.

Tsunami Along the Coasts of Japan
We examined how deep-ocean tsunami waves with a amplitude of 5 cm may be amplified in the coastal areas.Figure 7 shows the spatial distributions of the simulated maximum tsunami amplitude at 0-16 h after the main eruption for the four pressure conditions.Circles filled with grayscale represent JMA coastal gauges and the deviation between the simulated and observed amplitudes at these gauges.A negative value indicates that the simulated amplitude was lower than the observed amplitude.The Lamb wave-only condition in the upper left panel (Pressure L) shows a 10-cm tsunami along the eastern coast of Japan, while no significant amplification was observed in the other areas.Simulated tsunami amplitudes deviated from the observed maximum tsunami amplitudes at most points along the coast in this scenario, indicating that Lamb wave-only forcing could not reproduce this tsunami.
In contrast, Pressure A (top right panel), in which the forcing includes gravity waves, shows extensive amplification even in offshore areas.The bathymetry distribution (Fig. 5b) indicates that the large tsunami areas corresponded to the ridge lines.Differences between observed amplitudes were also significantly smaller than under Pressure L. However, the amplitudes of the three DART buoys in deep waters all were \ 5 cm.This result is consistent with the observations shown in Fig. 6 and indicates that the tsunami was greatly amplified during propagation to shallower areas.
Under Pressure B, the amplification area was narrower; instead, amplification was observed on a local spatial scale very close to the shore, and errors in simulated amplitudes were smaller on the east coast.This suggests that the maximum tsunami at these locations was caused by a bay-scale phenomenon rather than driven by offshore processes.
Pressure C had pressure properties similar to Pressure B, but the spatial distribution of the maximum tsunami was very different, owing to the difference in frequency bands.This indicates that the frequency band covered by Pressure C (10-20 min) did not cause bay-scale amplification as under Pressure B.
Figure 8 shows a comparison of simulated (Pressures L and A) and observed waveforms at the eight coastal gauges.The gauges that recorded particularly large tsunamis and the gauges with later peak occurrence times are displayed.The differences between the L and A scenarios indicate the effect of air gravity waves on sea surface elevation.The time of the leading wave (first peak) at each gauge corresponds to the time of the Lamb wave arrival, as confirmed in other studies (e.g., Tanioka et al., 2022;Yamada et al., 2022).
Overall, the Lamb wave-only forcing scenario reproduced the leading three or four peaks and deviated from the observed waveforms a few hours after the Lamb wave arrival.Subsequently, the water wave excited by the Lamb wave decayed, and the oscillation became significantly smaller.However, the simulation of Pressure A forcing, which included atmospheric gravity waves, showed that the subsequent waves were amplified and reproduced the observed waveforms well.In particular, the following sites reproduced the observed waveforms well, including the maximum amplitude: 4 Tanegashima, 18 Shirahama, 30 Chichijima and 37 Kuji.At 30  Chichijima, an amplitude of 80 cm was reproduced, which was assumed to be due to the effect of atmospheric gravity waves.As shown in Fig. 7, in Pressure A the tsunami amplified along the ridges and islands to a significant extent.
At Gauges 3 Amami and 9 Muroto, the amplitudes of the subsequent waves under Pressure A were larger than those under Pressure L but deviated significantly from the observed amplitudes.Even under the same pressure conditions, there was a large difference in reproducibility between 3 Amami and 4 Tanegashima, which are geographically close to each other.This result also suggests a predominance of tsunami amplification at the bay scale.In addition, at most gauges, the simulation amplitudes decreased approximately 6 h after the initial peak, and the reproducibility dropped significantly.This indicates that to reproduce a long-lasting tsunami as observed, a different type of amplification mechanism is needed than that present under Pressure A.
Figure 9 shows the simulated (Pressure B and C forcing) and observed waveforms from the same eight coastal gauges as in Fig. 8.These air pressures focused on the continuous forcing of the trailing waves, as described above.Because of the forcing of pressure waves with particular frequency bands, large oscillations were observed at several gauges.Under Pressure A, 12 h after the eruption, the oscillations were significantly smaller than the observed waveforms at 3 Tenagashima, 28 Muroto and 30 Chichijima.However, under Pressure B, resonances developed, and the amplitudes were equal to or larger than those observed.Differences in the occurrence of resonance between conditions B and C were due to the difference in the resonance period at each gauge.At gauge 37 Kuji, the simulated amplitude was larger than the observed amplitude because the specific frequency component was more intense, as shown in the time series and wavelet analyses (Fig. 4).It is also noteworthy that at gauges 3 Amami and 9 Muroto, no scenarios showed large amplification (conforming to the observational data), despite the long duration of the pressure waves with a wide frequency band.The frequency components of the air pressure waves contained frequencies that were predominant at the observed water level; however, no amplification was observed at these gauges.Thus, the pressure forcing with frequency components corresponding to the resonance period did not necessarily reproduce the amplitude of the observed water level.
The above results demonstrate that the observed maximum tsunami wave levels were successfully reproduced based on comparison with recorded data from a subset of gauges.The scenario that bestreproduced observations varied from gauge to gauge among the three pressure wave conditions incorporating Lamb and air gravity waves.At some gauges, Pressure A best reproduced the observed waveform because the maximum tsunami occurred within a few hours of the first arrival of the air gravity waves.However, where the maximum tsunami occurred later, Pressures B and C, which had longer radiation times, were better suited to reproducing the observed waveform.These results demonstrate the diversity of the amplification factors for this large tsunami.
Figure 10 shows a summary of the numerical tsunami simulations and the comparison with the observations.The observed and simulated maximum tsunami amplitudes for all 40 gauges are shown as boxes and bars.Bar colors represent the pressure conditions (A, B, or C) with the least error versus the observed data.Simulated results considering only Lamb waves are indicated by circles for reference purposes.This figure compares only the maximum values (from 0 to 16 h) and does not consider a comparison of peak occurrence times.Overall, the maximum amplitudes were reproduced well for gauges 27-40 in Eastern Japan, and Pressure B here exhibited the least error in the simulation.In contrast, there were only a few gauges with high reproducibility along the western coast of Japan, and Pressure A provided the lowest error results here.At gauges with maximum amplitude \ 20 cm, errors were almost below 5 cm.Particularly large tsunamis were observed at gauges 3 Amami, 8 Tosashimizu and 19 Kushimoto.Here, maximum simulated amplitudes differed from observed amplitudes by [ 30 cm.However, the simulated amplitudes under these conditions (A-C) were several tens of centimeters larger than those under the Lamb wave conditions alone.While tsunamis [ 50 cm are difficult to reproduce, we have shown that they may be reproduced at gauges such as 30 Chichijima by modeling air gravity waves to create waveforms similar to those observed.

Discussion
Numerical meteotsunami simulations were performed using air pressure forcing data modeled by a simple linear theory for air waves.In the simulation considering only Lamb waves in the pressure wave, the tsunami was relatively larger along the eastern coast of Japan and smaller along the western coast.A pressure wave \ 2 hPa reproduced the tsunami waveform approximately 1 h after arrival; however, its contribution to subsequent waves was negligible.This result is consistent with those of previous studies (e.g., Kubota et al., 2022;Nishikawa et al., 2022;Yamada et al., 2022).
The simulation results showed that the maximum amplitudes of air gravity waves, in addition to the Lamb wave, were successfully reproduced, particularly along the eastern coast of Japan.The east coast of Japan has a complex shape with numerous bays ranging from several kilometers to several tens of kilometers on a horizontal scale.Their geometry and the pressure waves were probably caused by the bay resonances at various locations.
To support the conclusion about the occurrence of bay resonance, a further simulation was performed using topography with a deliberately reduced horizontal resolution.While the original simulation used a maximum horizontal resolution of 5 arcseconds (approximately 130 m in the mid-latitude zone), this variant was performed at 1 arcminute (1/12 resolution).Figure 11  Especially in the eastern part of Gauges 27, 33, 34 and 37, the high-resolution simulation increases the amplitude.This amplification is more significant under the Pressure B condition.This is consistent with the intention of radiating long duration pressure waves to generate resonance in bays.
Figure 12 shows a comparison of the simulated waveforms between high-and low-resolution models for Pressure B. The high-resolution simulation reproduced the amplitudes to some extent, although the times of the peak occurrences (phases) were different.However, the low-resolution simulation showed a significant decrease in the amplitudes of subsequent waves.The phases of the waveforms also differed between resolutions, although the same air pressure conditions were used.The difference in tsunami amplitude between the high and low resolutions was larger after an elapsed time of 12 h when the pressure wave was weak.This indicates that the difference was not simply due to the effect of the detailed coastal depth profiles and that reproducing the substantial amplification seen in this tsunami does require a resolution sufficient to capture the bay resonance.
Using Gauge 37 Kuji as an example, the maximum amplitude at different resolutions is focused on.Figure 13 shows the spatial distributions of the maximum amplitude in the two simulations and their difference.The left and middle panels indicate the high-and low-resolution simulations, respectively.In the high-resolution simulation, Kuji Bay, a semicircular-shaped bay, was entirely amplified.On the other hand, the low-resolution simulation showed amplification only at the inner area of the bay.This result suggests the occurrence of bay resonance in the high-resolution simulation.Moreover, in the neighboring bays to the south of Kuji Bay, the differences in the maximum amplitude by resolution were small.This is probably because the resonance period of this bay, typically characterized by bay size and depth, is longer than that of Kuji Bay, which is not consistent with Pressure B. Therefore, the presence or absence of bay resonance is considered to be the main cause of the difference between the two simulations.
In addition to the resonance aspects, the effect of offshore tsunami amplification was quantitatively verified.The spatial distribution of the maximum tsunami amplitude (Fig. 7) shows that the amplitude increased from the Japan Trench toward the eastern coast of Japan, probably because of wave shoaling effects.In contrast, the maximum amplitude offshore of the trench axis in the deep ocean was approximately 5 cm and never exceeded 10 cm, which is consistent with the findings of Suzuki et al. (2023).Even though the tsunami excited offshore was small, Pressures A and B reproduced the amplitudes of the subsequent tsunamis.This suggests that phenomena occurring in the offshore basin, such as Proudman resonance, were not the primary cause of tsunami amplification along the eastern coast of Japan.This result was in contrast to South America, where the contribution of the Proundman resonance was large (Ren et al., 2023).We also found that even weak air pressure waves (with an amplitude of 0.5 hPa) could amplify the tsunami to [ 50 cm in the coastal area because of the bay resonances.
Amplification by island chains and ocean ridges was prominent in the spatial distribution of the maximum amplitude (Fig. 7).These bathymetries excite waves and constrain their energy based on the angle and frequency of the incoming tsunami (Buchwald & Adams, 1968;Buchwald & Longuet-Higgins, 1969).The simulated results under Pressure A forcing showed strong amplification in shallow areas and reproduced observed tsunami waveforms fairly well.In particular, at Gauge 30 Chichiijma, tsunamis were amplified over the entire area surrounding the Ogasawara Islands.This suggests that amplification caused by topographic features is one of the primary causes of large tsunamis.In addition to the modeled tsunami event, there have been other events in which tsunamis were amplified on Chichijima because of the topography (Koshimura et al., 1999a, b).
These excitations due to bay resonance and seafloor topography are highly dependent on frequency characteristics.During the eruption, the pressure waves were estimated to have radiated over a very wide range (Kulichkov et al., 2022;Omira et al., 2022;Tarumi & Yoshizawa, 2023).Generally, a sharper pulse results in a wider spectral width (as in the Fourier transform of Dirac's delta function).The modeled eruption event also represented an abrupt change in atmospheric conditions, which resulted in the formation of pulse waves and the generation of strong pressure waves over a wide frequency range (e.g., Nishikawa et al., 2022).This extended spectral width may have generated pressure waves prone to resonance in any region.
Although the waveforms were well reproduced at several gauges, the amplifications at Gauges 3 Amami, 8 Tosashimizu and 9 Muroto, for example, could not be explained under any pressure conditions.To demonstrate that the reason for the deviation is not simply the horizontal resolution of bathymetry, an additional simulation was performed for the Muroto region as an example.Bathymetry data with a higher resolution of 1 arcsecond around Gauge 9 Muroto (from the CDMC dataset) were added to the simulation.Thus, a total of six levels of resolution were used.The simulation was compared with the original 5 arcsecond resolution simulation to confirm the difference in tsunami amplitude due to the resolution.
Figure 14a shows the comparison in waveforms between these two simulations for Pressure A condition.The blue line represents the original simulation, identical to that in Fig. 8.Although the additional simulation was performed at 1 arcsecond resolution, the tsunami amplitude did not increase.Figure 14b shows the spatial distribution of the difference in the maximum amplitude between the two simulations.The only regions where the amplitude increased with higher resolution were those with the shape of a bay.In other regions, including Muroto, the amplitude increased approximately several centimeters even at the shore.It is possible to perform even more high-resolution simulations.However, we think that a finer resolution will only be effective in reproducing amplitudes with shorter periods and will not lead to a fundamental solution to the large discrepancy from the observed amplitudes.Similar results where increasing the resolution did not improve the amplification were obtained for Amami and Tosashimizu.
To investigate the reason for the large deviation of simulated amplitude, wavelet analysis of the observed tsunami waveforms at these two gauges was performed.Figure 15 shows the results of these analyses.Amplitudes observed at both gauges exceeded 70 cm; however, they showed different frequency characteristics.
For Amami (Fig. 15a), the predominant component was concentrated between 10.5 and 11.5 relative hours in the time domain and amplified in a wide range of 10-30 min in the frequency domain.By contrast, for Tosashimizu (Fig. 15b), the predominant component occurred at 10.5-16.0h relative time, and its dominant period was relatively narrow, ranging from 20 to 30 min.A comparison with the wavelet analysis of air pressure provided information about the excited waves.Naze, an air pressure gauge, and Amami, a tidal gauge, were located on the same island (Amami Oshima), which would suggest that both should record an almost identical pressure waveform.Wavelet analyses of the pressure waveforms (Fig. 4) show that Pressure A features an intense component in the corresponding period of predominance at Amami.Wavelet analysis of Pressure B indicates that a component corresponding to the predominant period at Tosashimizu was incident for a long period.Pressures A and B were higher than the actual observed pressure waveforms; however, the observed tsunami simulation waveform was significantly amplified, whereas the simulated waveform was not.One reason for this discrepancy is the systematic biases in bathymetry and topography data.Another reason may be the existence of secondary amplification factors in the oceanic tsunami propagation process, such as wave trapping due to island topography (Longuet-Higgins, 1967).Amplifications at these sites may thus not have been caused directly by the atmospheric pressure wave.
The assumption used in constructing both the atmospheric and tsunami models in this study results in some limitations of their application to complex and realistic phenomena.Modeling atmospheric pressure waves based on the vertical structure of the atmosphere may improve the simulated tsunami results.Coriolis force and background winds were disregarded in the pressure model.In addition, buoyancy frequency and sound velocity, which vary with the temperature structure, were assumed to be constant.These parameters were tuned to match the waveforms observed in Japan.These assumptions are the same as those used in Wright et al. (2022) and Gusman et al. (2022).However, such linear superposition and simplification of parameters ignore the possible propagation paths of atmospheric gravity waves, which may be a source of error in tsunami waveform reconstruction.
Furthermore, under the simple assumptions of this study, it was difficult to match the amplitudes with the observations when considering a wide frequency band and long-lasting waves with features of both Pressure A and Pressure B. Acoustic waves, which constitute another type of gravity waves in a linear wave solution, were also neglected, although it has been suggested that they contribute to sea-surface excitation (Nishikawa et al., 2022).Numerical simulations based on a sophisticated atmospheric model that implements a three-dimensional structure of the atmosphere (e.g., Watanabe et al., 2022) can address these issues and may lead to a better understanding of tsunami events.The tsunami model was based on nonlinear long waves (as in Omira et al., 2022) and did not consider dispersion.The wavelength and period were shorter than those of tsunamis generated by typical subduction earthquakes (Lynett et al., 2022); incorporating dispersion (e.g., Baba et al., 2017;Kirby et al., 2013) may be effective in reproducing such tsunamis.In this tsunami event, the significance of the dispersion is found to be negligibly small in the far field in deeper waters (Ren et al., 2023).The significance for shallower and shorter wavelength components should also be investigated in the future.
Furthermore, the underlying bathymetry data used in the model may not be precise, which in turn may affect the reproducibility of tsunamis.In particular, the resonance phenomena prominent in the modeled event were closely related to water depth.It would be useful to validate the bay resonance period obtained from the numerical experiments in this study by comparing it with the resonance period empirically obtained from other tsunami events (e.g., Rabinovich, 1997;Ren et al., 2021) or with a theoretical method (e.g., Lee, 1971).In addition to the resonance phenomena, a detailed investigation of the characteristics of long-wave trapping and the amplification conditions of specific frequencies at each location on the islands may help understand the amplification factors of tsunamis in this region.

Conclusions
In this study, we conducted simulations of tsunamis on the Japanese coast caused by the 2022 Hunga Tonga-Hunga Ha'apai volcano eruption on the Japanese coast, where the tsunamis were significantly amplified in many areas by air pressure waves from the eruption.Linearized wave theory was used to generate atmospheric pressure data, atmospheric Lamb waves and gravity waves were modeled, and dispersion relations were considered for gravity waves.A spatiotemporal pressure wave model was developed to reproduce atmospheric pressure fluctuations near Japan by superpositioning waves with various combinations of wavelengths and amplitude parameters.The numerical tsunami simulations successfully reproduced the maximum tsunami amplitudes at various locations along the Japanese coast, while this proved difficult to do using only Lamb waves or gravity waves of a single wavelength.The major conclusions from this study are as follows: • The effect of surface waves excited by the Lamb wave dominated the surface water fluctuations in Japan immediately after the Lamb wave passed until approximately 1 h later.Following this, the influence decreased, and their contribution to tsunami amplification at the time of the maximum tsunami was very small at most sites.• Tsunami simulation results under conditions using composites of Lamb waves and air gravity waves reproduced tsunami waveforms and maximum tsunami amplitudes well at many locations.The gravity waves likely amplified the tsunami in shallow areas such as ocean ridges and island chains.• At several locations, bay resonance appears to have been the primary cause of tsunami amplification, especially on the east coast of Japan.In contrast, the maximum amplitudes in offshore basin regions did not exceed 10 cm, and the tsunami was amplified as it moved from the trench axis toward the coast due to shoaling effects.Amplification in deep offshore areas, such as by Proudman resonance, was likely not a primary amplification factor, although it may have been contributory.• At several locations in western and southwestern Japan, the amplification could not be reproduced even though the pressure data included a pressure component equivalent to the predominant period of the observed tsunami waveforms.Even with the additional high-resolution simulation (1 arcsecond), its contribution to the tsunami amplitude was small.This suggests that the cause of the tsunami amplification at these locations was not direct atmospheric pressure excitation but rather an amplification of a specific frequency band generated during wave propagation in the ocean.
These results were obtained by modeling tsunami behavior at a sufficiently fine scale over a wide area along the coast of Japan.The AMR method allows for detailed and efficient tsunami simulations in these areas.Models of this resolution are clearly useful since the numerical results show diverse amplification factors, and the major amplification factors differed even at gauges close to each other.The methodology and simulations used in this study are readily applicable to other regions.If separation of the tsunami component caused by the caldera collapse and the meteotsunami component is desired, these could be performed in regions closer to the volcano.Further clarification of the causes of tsunami amplification requires an investigation of the amplification characteristics of long waves in island chains and ridge topography, pressure wave data based on a detailed atmospheric model and the verification of bathymetry data bias.

Figure 1
Figure 1 Location of JMA tidal or tsunami gauges around the Pacific coast of Japan.Red circles represent JMA tidal or tsunami gauges.Number labels represent gauge numbers used for identification in this study.Detailed gauge locations and observed tsunami amplitudes following the modeled eruption are shown in Table1.The two yellow squares represent DART buoys off the coast of Japan

Figure 2
Figure 2 Dispersion relation of air acoustic and gravity waves.The vertical and horizontal axes are non-dimensionalized on the representative scales.Gray-shaded areas represent the range of acoustic and atmospheric gravity waves, respectively

b
Figure 3 Snapshots of spatial distribution of air pressure field generated in the simulation.Four patterns of pressure waves were created and are shown as sub-panels at 8 h after the eruption.The red triangle in the South Pacific represents the volcano from which the pressure waves radiate bathymetry.The longitudinal and latitudinal extent of the domain was from 115 to 200 E and from 55 S to 55 N. Yellow boxes indicate the DART buoys available during the tsunami.The waveforms recorded at these DART buoys were used to validate offshore results.

b 5
Figure 6Comparison of simulated and observed pressure-forced waves

Figure 7
Figure 7Horizontal distributions of simulated maximum amplitude for the modeled pressure conditions.Circles filled with grayscale shading represent JMA coastal gauges and the deviation between the simulated and observed amplitudes at those gauges.Negative value indicates a simulated amplitude lower than observed amplitude.Yellow boxes represent offshore DART buoys.Note that the number of JMA gauges displayed has been reduced for visibility

Figure 8
Figure 8 Comparison between simulated and observed tsunami waveforms at coastal gauges.Black lines represent observed waveforms.Red and blue lines represent simulated waveforms for Pressure L and A conditions, respectively

Figure 9
Figure 9 As in Fig. 8 but for Pressure B and C conditions

Figure 10
Figure 10 Comparison of maximum amplitude between simulated and observed results.Only simulated results for the condition with the least error in maximum amplitude among the three pressure conditions are shown.Conditions resulting in the least error are indicated by the color of the bars.Black boxes represent observed maximum amplitude.Simulated results for the Lamb wave-only condition (Pressure L) are plotted as circles for reference

Figure 11
Figure 11 As in Fig. 10 but for all simulated results.Horizontal hatched bars indicate the added amplitude of the high-resolution (5 arcseconds) simulations from the low-resolution (1 arcminute) simulations

Figure 12
Figure 12Comparison of waveforms between the high-resolution (5 arcseconds) and low-resolution (1 arcminute) simulations for Pressure B forcing condition.Five coastal gauges on the eastern coast of Japan were selected to evaluate dependency on horizontal resolution.The orange lines represent the high-resolution results (identical to Fig.9for Pressure B but different gauges), and the blue lines represent the low-resolution results shows the difference of amplitude between high-and low-resolution results.This figure is generally the corresponding figure to Fig. 10 but shows the simulated results for Pressures A through C. High-resolution results are indicated as stacked bars hatched with lines, while low-resolution results are indicated as bars filled with plain colors.

Figure 13
Figure 13 Spatial distributions of the maximum amplitude with the high (left) and low (middle) resolution simulations and the difference (right).The yellow circle of each panel represents the location of Gauge 37 Kuji

b
Figure 15 Wavelet analysis of observed tsunami waveforms for gauge a Amami and b Tosashimizu, located in southwestern Japan.The upper panel shows the tsunami waveform and the lower panel shows the wavelet (f-t) diagram

Table 1 .
The two yellow squares represent DART buoys off the coast of Japan

Table 1
Properties of JMA tidal gauges and maximum amplitudes

Table 2
Outline of pressure conditions generated in this study