1 Introduction

The Hunga Tonga–Hunga Ha’apai (HTHH) volcano in the southwest Pacific Ocean has been active since 2021, and produced a series of massive eruptions on January 15, 2022 (Smithsonian Institution, 2023). The Volcanic Explosivity Index (VEI) is estimated to be 5, meaning that the total volume of ejecta was between 1 and 10 km3. Satellite images on January 15 showed that the ash plume reached an altitude of 30 km (Witze, 2022). This eruption caused disturbances in the ionosphere, atmosphere, ocean, and solid Earth. The total ionospheric election content was observed and reproduced by simulations (Astafyeva et al., 2022; D’Arcangelo et al., 2022). In the atmosphere, gravity waves of different speeds, the Lamb wave and the Pekeris wave, were detected (Suzuki et al., 2023; Watanabe et al., 2022). In the ocean, tsunamis were recorded on coastal tide and ocean bottom pressure gauges throughout the Pacific Ocean (Kubota et al., 2022). Seismic waves, both direct from the source and those coupled to atmospheric waves, were observed globally (Matoza et al., 2022).

The atmospheric Lamb wave, first proposed in 1911 (Lamb, 1911), travels at about the speed of sound (~ 310 m/s) with little frequency dispersion because it is the atmospheric boundary wave (Harkrider & Press, 1967). The Pekeris wave, an atmospheric gravity wave first proposed in 1937 (Pekeris, 1937), travels at a slower speed (200–250 m/s) and exhibits some dispersive characteristics. The vertical profiles show that the Lamb wave amplitude exponentially decays with height, while the Pekeris wave amplitude has a node at ~ 15 km and becomes largest at ~ 50 km height (Watanabe et al., 2022).

The tsunami records from this eruption show different characteristics in the near-field, such as Tongatapu Island, and in the far-field, such as Japan. The tsunami was recorded on two tide gauges in Nuku’alofa, about 70 km from the HTHH volcano, starting at 4:25 UTC with a maximum amplitude of 1.9 m at 4:46 UTC (Borrero et al., 2023). Sea level disturbances on far-field tide gauges and ocean bottom pressure gauges started earlier than the expected tsunami arrival times. In Japan, the tsunami arrivals at Chichijima and Amami were recorded at about 11 and 12 UTC, respectively, a few hours earlier than the expected tsunami arrival times (Imamura et al., 2022; Japan Meteorological Agency, 2022). Tsunami heights exceeded 1 m at Amami and Kuji, and the Japan Meteorological Agency issued a tsunami warning (Imamura et al., 2022).

The observed far-field tsunamis have been modeled by many researchers (Gusman et al., 2022; Lynett et al., 2022; Omira et al., 2022). The early arrival of the tsunami has been interpreted as being excited by the Lamb wave (Kubo et al., 2022; Kubota et al., 2022). Similar air-coupled sea waves have also been observed and modeled from the 1883 Krakatoa eruption (Harkrider & Press, 1967). To explain the observed maximum amplitude that appeared in the later phase, we need to consider the coupling effect between the atmospheric wave and the tsunami. The Proudman resonance is expected when the speed of the tsunami (long wave) is similar to the propagation speed of atmospheric waves. For the Lamb wave speed of 310 m/s, the corresponding water depth of the long wave is ~ 9800 m. On the other hand, the Pekeris wave speeds of 200–250 m/s correspond to the water depths of ~ 4000 to 6400 m. Considering that the average depth of the Pacific Ocean is about 4 km, the Pekeris wave can couple with the ocean to generate tsunami waves. The Proudman resonance effects for the Lamb wave were modeled in one-dimensional simulations with some bathymetry shapes by Sekizawa and Kohyama (2022) and also discussed by Tanioka et al. (2022). Those for the Pekeris wave were modeled in two-dimensional computations by Nishikawa et al. (2022), demonstrating the amplified tsunami waves in resonance. Mizutani and Yomogida (2023) also successfully reproduced the amplified later phases using a hybrid method of tsunami computation with the normal mode theory for atmospheric waves.

The purpose of this paper is to model the tsunamis observed both on ocean bottom pressure gauges and coastal tide gauges around the Pacific Ocean, from propagating atmospheric waves, namely, the Lamb wave and the Pekeris wave. We model the atmospheric waves without any dispersions, i.e., the propagation speed is constant along the path. We then estimate the parameters, such as the origin time and pressure wave amplitudes, and the propagation speeds of the Lamb and Pekeris waves.

2 Data

2.1 Tsunami Waveforms at DART Stations

The tsunami from the HTHH volcanic eruption was recorded at ocean bottom pressure (DART) and coastal tide gauge stations in and around the Pacific Ocean. Figures 1 and A-1 show the locations of the DARTs and tide gauges, where observed data were collected and used for analyses in this study. The U.S. DART data (Fig. 2) have increased sampling rates for the event mode: initial 15 s, following by 1 min section, and 15 min before and after the event mode. The New Zealand DART data (Fig. 3) were provided with a sampling rate of 15 s. We pre-processed the observed pressure data to extract tsunami signals by removing ocean tide signals fitted by polynomial functions which have longer periods than the tsunami signals. Then, the extracted tsunami waveforms at DARTs were resampled with a time interval of 1 min and used for inversions of the initial phases or waveform analyses, such as correlation coefficients for the later phases.

Figure 1
figure 1

Distribution of the DART and tide gauge stations in and around the Pacific Ocean from the source region of the HTHH volcano (red star). The purple squares are the New Zealand DARTs, and the other colors are the U.S. NOAA’s DARTs. The green, light-blue, blue, and purple colors correspond to the groups shown later in Figs. 2 and 3. The yellow triangles are selected tide gauges, shown later in Fig. 4, at Crescent City in North America, Honolulu, Callao in South America, and other stations in Japan. The area shown on this map corresponds to the computation area for the tsunami simulations

Full size image
Figure 2
figure 2figure 2figure 2

a Comparison of tsunami waveforms at the U.S. NOAA’s DART stations (from near source to Hawaii) shown as the green squares in Fig. 1. Time is measured from the origin time (4:16 UTC) estimated in this study. The red curves are the observed waveforms, and the light-blue ones are the assumed atmospheric pressure (Lamb and Pekeris) waves. The black curves behind the blue ones are the synthetic bottom pressure waveforms which are the sum of the assumed pressures due to the Lamb wave (amp0 = 383 hPa) and the pressure changes calculated from water levels (Sects. 3.2 and 4.1). The blue curves in front are the synthetic bottom pressures for both the Lamb and Pekeris (amp0 = − 50 hPa) waves (Sects. 3.3 and 4.2). The vertical dark-gray bars show the arrival times of the Lamb wave in the case of the propagation speed (C1) of 310 m/s. The vertical light-gray bars show the arrival times of the Pekeris wave with the best propagation speeds (C2), as listed in Tables A-1 and A-2. The vertical green bars are the expected arrival times of a tsunami (long wave) from the HTHH volcano. The thick dark-gray horizontal bars are the inversion ranges. The thick light-gray horizontal bars show the time ranges used to calculate the correlation coefficients (CC) of the observed and modeled later phases. b The DART stations (Japan–Alaska) shown as the light-blue squares in Fig. 1. c The DART stations (North–South America) shown as the blue squares in Fig. 1

Full size image
Figure 3
figure 3

Same as Fig. 2, but for the tsunami waveforms at New Zealand DART stations around the source, which are shown as the purple squares in Fig. 1

Full size image

The initial phases of the observed tsunami waveforms (red curves in Figs. 2 and 3) have positive amplitudes and correspond to theoretical arrival times (vertical dark gray bars) in a case of the atmospheric Lamb wave speed of 310 m/s. Note that the amplitudes are not the pressure values measured by the DART sensors, but the pressure changes due to atmospheric waves and/or tsunami. The peak amplitudes were several hectopascals at most DARTs, and the maximum of 16 hPa was recorded at NZG, the New Zealand DART closest to the volcano. At most of the DARTs, the first half cycles of the waves were triangular in shape with half period of ~ 20 min, but were not visible at some stations due to the coarse sampling rate. The later phases have negative onsets at the theoretical arrival times of the Pekeris wave (vertical light-gray bars, in a case of the best propagation speed for each station to be estimated in Sect. 4.2) with periods shorter than the initial phases. In the later phases, the maximum amplitudes appear either before or after the theoretical arrival times (vertical green bars) of a long wave from the HTHH volcano and show positive values at most stations but negative ones at some stations. We systematically read the arrival times (wavefront, not a crest) and the maximum amplitudes of the initial phases, the maximum positive/negative amplitudes and their appearance times of the later phases, and listed them with the station locations and distances from the HTHH volcano in Tables A-1 for the U.S. DARTs and A-2 for the New Zealand DARTs. The location data for the U.S. and New Zealand DARTs were obtained from the NOAA’s website for the real-time DART data and were provided by GNS Science, New Zealand, respectively.

2.2 Tsunami Waveforms at Tide Gauge Stations

The collected tide gauge data were also de-tided to extract the tsunami signals as described in the previous section for the DART data. The data have different sampling rates of 15 s, 30 s, or 1 min, depending on the station. The tsunami waveforms are shown for the selected far-field tide gauge stations (red curves in Fig. 4), for the near-source stations (Fig. A-2), and for the other stations in Japan (Fig. A-3). At the far-field stations, the initial phases with positive onsets appear several minutes or several tens of minutes after the theoretical arrival times (vertical dark-gray bars) of the Lamb wave (~ 310 m/s), with amplitudes ranging from a few to several tens of centimeters. In contrast to the DART data (Figs. 2 and 3), the initial peaks are followed by negative troughs of almost equal or greater amplitude. The half periods are ~ 20–30 min or longer at most stations (e.g., Kushiro, Aburatsu), but shorter at some stations (e.g., Nukualofa, Norfolk) near the source. At most stations, there are no clear arrivals of later phases, and the amplitudes increase with time. The maximum amplitudes of the later phases were recorded after the theoretical arrival times (vertical green bars) of the long wave from the HTHH volcano on most tide gauges. However, at some Japanese stations (e.g., Amami, Tosashimuzu), the maximum amplitudes appeared before the theoretical arrival times of the tsunami, but after the theoretical arrival times (vertical light-gray bars) of the Pekeris wave. The maximum amplitudes of the later phases are several centimeters to several tens of centimeters, but they exceed one meter at some stations (e.g., Amami, Crescent City). The maximum amplitudes were also negative at most stations. Similar to the DART data, we read the arrival times of the initial rises, the peak amplitudes of the initial phases, the maximum positive/negative amplitudes of the later phases, and their appearance times, and listed them in Table A-3 for the near-source (HTHH volcano) stations and Table A-4 for the far-field stations. When we placed the tide gauge locations on the computational bathymetry grid, the actual locations of some stations did not coincide with a wet grid along the coastal grid. Therefore, we show the tide gauge locations in the tables as the waveform output points on the bathymetry grid data with a grid size of 24 arc-sec (0.2 arc-min). They were selected along the coast regarding the actual tide gauge locations and used for the tsunami simulations.

Figure 4
figure 4

Comparison of tsunami waveforms on selected far-field tide gauges shown as the yellow triangles in Fig. 1. Time is measured from the origin time (4:16 UTC) estimated in this study. The red curves are the observed waveforms in meter, and the light-blue ones are the assumed atmospheric pressure (Lamb and Pekeris) waves in 100 hectopascal. The black curves behind the blue ones are the synthetic waveforms which are the water levels due to the Lamb wave (amp0 = 383 hPa). The front blue curves are the synthetic waveforms for both the Lamb and Pekeris (amp0 = − 50 hPa) waves. The vertical dark-gray bars show the arrival times of the Lamb wave in the case of the propagation speed (C1) of 310 m/s. The vertical light-gray bars show the arrival times of the Pekeris wave at the best propagation speeds (C2) of the nearby DART stations, as listed in Tables A-1 and A-2. The vertical green bars are the expected tsunami arrival times from the HTHH volcano

Full size image

3 Tsunami Modeling Method

3.1 Tsunami Computations

The effects of atmospheric pressure on tsunami propagation can be expressed by additional terms in the equations of motion (e.g., Kubota et al., 2022; Gusman et al., 2022). Assuming hydrostatic pressure, the atmospheric pressure change dP due to the Lamb or Pekeris wave can be converted to the water height change as, \({dh}_{p}=dP/\rho g\), where ρ is the density of seawater (assumed to be 1020 kg/m3), and g is the gravitational acceleration (9.8 m/s2). This equation indicates that an atmospheric pressure change of 1 hPa (102 Pa) is equivalent to a water level change of 10–2 m; hereafter, we refer to it as the “assumed atmospheric pressure (AP) due to the Lamb (or Pekeris) wave”.

We performed numerical tsunami simulations in the region shown in Fig. 1. The computation area is 240° W–70° W and 55° S–62° N. The 30 arc-sec bathymetry grid data from GEBCO 2014 (Weatherall et al., 2015) were resampled with a 24 arc-sec grid interval; thus, the numbers of the computation grid are 25,800 and 17,400 in the longitude and latitude directions, respectively. We numerically solved the linear long-wave equations (Satake, 1995) in the spherical coordinate system for 20-h tsunami propagation with a time interval of 1 s, which satisfies the stability condition. The computation time was about 13 h with GPGPU (TESLA V100 32 GB and CUDA 9.0), as used in Satake et al. (2017) and Fujii et al. (2020).

The atmospheric Lamb and Pekeris waves were modeled as propagating concentric waves with constant speeds, C1 and C2, respectively (Fig. 5). At each time step, the pressure distribution from the Lamb or Pekeris waves was used as the input for the tsunami propagation. The assumed APs are characterized by an initial amplitude (amp0) at the source (HTHH volcano) and a rise time (T). The amplitude (amp) of the APs decays with a distance as \(1/\sqrt{R{\text{sin}}\Delta }\), where R is the Earth radius (6378.137 km), Δ is the epicentral angle distance from the initial grid point. This is an asymptotic solution for the 2D nondispersive wave equation on a sphere, where the wavelength is much smaller than the radius of the propagating ring \({\text{r}}=R{\text{sin}}\Delta\), although it is not accurate at the initial stage of the propagation. In the actual simulations, we set amp0 as a unit amount at the grid point of the epicenter (the center of the HTHH volcano). For all other grid points, we set the amps calculated by \({\text{amp}}= {\text{amp}}0 /\sqrt{2\pi R{\text{sin}}\Delta /dr}\), where dr is the grid size. Since the computational grid size is 24 arc-sec (about 720 m), the amp does not diverge even at grid points immediately adjacent to the epicenter grid point. Watanabe et al. (2022) and Suzuki et al. (2023) demonstrated, theoretically and observationally, that the fast propagating (at ~ 310 m/s) Lamb wave has a positive amplitude, while the slower (200–250 m/s) Pekeris waves have negative amplitudes. The animation of the tsunami generation by the Lamb wave only shows that the positive wave in red propagates first at the Lamb wave speed, followed by the propagation of positive (red) and negative (blue) amplitudes at the tsunami speed (https://iisee.kenken.go.jp/staff/fujii/Tonga2022/tsunami.html).

Figure 5
figure 5

Schematic figures of the assumed atmospheric pressure (AP) waves propagating from the source. The left figures show the temporal change at each point of the computational grid, while the right figure shows a map view

Full size image

We added the assumed AP due to the propagation of atmospheric waves (either the Lamb wave only or both the Lamb and Pekeris waves) to the calculated water levels of the numerical tsunami simulations (Figs. 2 and 3). Since the water level change due to the Lamb or Pekeris wave is observed as a pressure change at the ocean bottom, the synthetic waveforms at DART stations (black and blue curves in Figs. 2 and 3) are the sum of the assumed AP (light-blue curves) and the associated tsunami. On the other hand, at the tide gauge stations (Figs. 4, A-2, and A-3), the synthetic waveforms (black and blue curves) show only the water level change due to the tsunami, because the pressure changes are not observed on the tide gauges.

3.2 Tsunami Waveform Inversion of the Initial Phase

We estimated the origin time, the amp0 at the source, and the speed of the propagating concentric wave (the Lamb wave), by inverting the initial phases observed at 14 U.S. and nine New Zealand DART stations. Synthetic tsunami waveforms are linear with respect to amp0; hence, amp0 can be determined from comparisons between the observed and synthetic waveforms if an origin time is fixed. We estimated the amp0 using the non-negative least squares method (Lawson & Hanson, 1974), and evaluated the degree of waveform agreement by a variance reduction (VR). We adopted three constant propagation speeds (C1 = 305, 310, and 315 m/s) of the Lamb wave without dispersion. For all three C1s, amp0s were obtained by changing the origin times between 4:00 and 4:26 UTC with a 1 min increment. We estimated the best origin time, the amp0, and the C1 that give the maximum VR among the three cases.

3.3 Waveform Analysis for the Later Phase

To simulate the tsunamis caused by the Pekeris wave, we assumed the same parameters as the Lamb wave, i.e., the origin time (4:16 UTC) and the time scale (T of 10 min as described later in Sect. 4.1), except for the propagation speeds (C2) and a negative amp0. As mentioned above, the opposite initial amplitude of the Pekeris wave from the Lamb wave has been confirmed by both theory (Watanabe et al., 2022) and observation (Suzuki et al., 2022). We varied the C2 of the Pekeris wave from 200 to 250 m/s, the range covering the wave speeds found in the previous study (Watanabe et al., 2022), with an interval of 5 m/s. We subtracted the synthetic tsunami waveforms caused by the Lamb wave from the observed waveforms and used the residual waveforms in the later phase analyses. Here, we assumed that the residual waveforms represented the later tsunami caused by the Pekeris wave.

As described in the previous section, we used VR as an indicator of waveform agreement for the initial phases when inverting the amp0 of the Lamb wave. In this section, we used a correlation coefficient (CC), which is independent of the initial amplitude (amp0) value, as an indicator of waveform reproducibility. We first estimated the best or dominant C2 for each DART station. We tentatively set amp0 = − 100 hPa and calculated the CCs between the residual waveforms and the synthetic ones for the Pekeris wave. Figure 6 shows the process of calculating the CCs for the later phase of a DART record. In this example case, the C2 of 240 m/s gives the highest CC of 0.53, so we chose this as the C2 for this DART station. After determining the best C2 for each DART station, we estimated the amp0 of the Pekeris wave by calculating the root mean squared errors (RMSE) or VRs for all the DART stations.

Figure 6
figure 6

An example of determining the best or dominant propagation speed (C2) of the Pekeris wave for a DART station. The correlation coefficients (CC) were calculated between the residual waveform (observed waveform − synthetic waveform from the Lamb wave) and synthetic waveforms from the Pekeris wave with different C2s

Full size image

4 Results

4.1 Lamb Wave Parameters

The target parameters of the atmospheric Lamb wave were the propagation speed (C1), the initial amplitude (amp0) at the source, and the origin time. Besides these parameters, we fixed the rise time (T) to 10 min to reproduce the time scales of ~ 20 min from the observed initial peaks at available DART stations. The corresponding peak width (a half wavelength) of the propagating concentric wave is ~ 370 km (C1 * 2T) for a case of C1 = 310 m/s.

From the comparison between the observed and synthetic waveforms in the case of C1 = 310 m/s, the maximum VR of about 0.6 was obtained for the origin time of 4:16 UTC. In this case, the amp0 was estimated to be 383 hPa (Fig. 7). If we slightly change the C1 to a slower value of 305 m/s or a faster value of 315 m/s (Fig. A-4), then the VRs become worse. Thus, we concluded that the best C1 is 310 m/s. The initial phases of the observed tsunami waveforms were well reproduced at the U.S. and New Zealand DART stations (see red and black curves in Figs. 2 and 3).

Figure 7
figure 7

Variance reductions (VR) and initial amplitudes (amp0) at the source versus origin times for the Lamb wave propagation speed (C1) of 310 m/s in inversions for initial phases. The filled triangles and the white circles are for VR and amp0 from inversions, respectively

Full size image

The waveform comparisons are also made for the tide gauge stations around the source region (Fig. A-2) and the selected stations (Fig. 4) in North America (Crescent City), Hawaii (Honolulu), South America (Callao), and Japan, and the other stations in Japan (Fig. A-3). The initial phases are well modeled by the synthetic tsunamis caused by the Lamb wave, while the later phases are not well reproduced, indicating that we need another source caused by the Pekeris wave. Moreover, the waveforms for the stations near the source area are very different, which means that the tsunamis recorded at these stations require some other mechanism, such as a volume change of the volcanic mountain body (Heidarzadeh et al., 2022; Pakoksung et al., 2022), rather than the propagating atmospheric waves. The amplitude of the assumed AP (amp) calculated at distant stations are several hectopascals in Japan (light-blue curves in Figs. 4 and A-3: the amplitudes are too small to see) and are consistent with the pressure wave amplitudes of ~ 2 hPa observed on barometers (Watada et al., 2023; Watanabe et al., 2022).

4.2 Propagation Speeds of Pekeris Wave

The additional propagating concentric wave corresponding to the atmospheric Pekeris wave generated the later phases of tsunami observed at the far-field DART stations, as shown by the difference between the black curves (Lamb wave only) and the blue curves (Lamb and Pekeris waves) in Fig. 2. The negative amplitude following the initial positive peaks were observed at DART stations towards Japan (51425, 52406, 52401 in Fig. 2a, 52404, 21418, 21419 in Fig. 2b), towards North America (46408, 46402, 46403, 46414, 46409 in Fig. 2b, 46411 in Fig. 2c), and towards South America (32413, 32402, 32404 in Fig. 2c).

The synthetic tsunami waveforms shown in Figs. 2 and 3 are those that best match the observations among the simulations with different propagation speeds (C2) and the amp0 of − 50 hPa. The CCs between the residual and synthetic waveforms with different C2s are shown in Fig. A-5. We picked up the best C2 for each path with the largest CC (Fig. 8), and listed them in Tables A-1 and A-2. The best C2s seem to correlate not only with distances but also with azimuths. Toward the northwest, Japan, Kamchatka, and Alaska, the faster C2s of 235–250 m/s give the better CCs, while toward the northeast, North America, and the southeast, South America, the slower C2s of 200 or 210 m/s give the better CCs. The speed variation of the Pekeris wave may be affected by the temperature and wind speed in the atmosphere, much more than the Lamb wave (Watanabe et al., 2022).

Figure 8
figure 8

The best propagation speeds (C2) of the Pekeris wave to reproduce the later phases observed at DART stations

Full size image

We calculated the CC at all the DART stations and confirmed that the CC is independent of the given amp0 (Fig. A-6) of the Pekeris wave. Then, to estimate the amp0, we calculated RMSEs and VRs (Fig. A-7) using the same waveform data as we used for the CC calculations. Although the VR values were small, we found that the waveform reproducibility was better for the Lamb + Pekeris waves than for the Lamb wave only (amp0 = 0 hPa in Fig. A-7). We also found that the amp0 of − 50 hPa was the best value. In the waveform comparisons shown so far (Figs. 2, 3, 4, A-2, and A-3), the synthetic waveforms were calculated with the amp0 of − 50 hPa.

The entire waveforms, not only the initial phases but also the later phases with their maximum amplitudes, at the far-field tide gauge stations (Figs. 4 and A-3), such as Crescent City, Honolulu, Callao, Kushiro, Chichijima, and Aburatsu, are well produced by the Lamb wave with the C1 of 310 m/s and the Pekeris wave with a suitable C2 for each station, which was determined at the DART station located offshore closest to that tide gauge station. However, on some tide gauges, such as Chichijima and Aburatsu, the synthetic tsunami waveforms of only the Lamb wave reproduce the observed ones relatively well, and the improvements with the addition of the Pekeris wave are not so significant. The later large tsunami waves are still underestimated on other tide gauges, such as Amami, Kushimoto, and Tosashimizu. For these stations, other tsunami computations using a finer bathymetry grid and considering the local effects, such as a natural oscillation of each bay, may be needed to reproduce the observed waveforms better.

5 Discussion

In this study, we successfully reproduced the observed tsunami waveforms by adding atmospheric pressure (AP) changes of the propagating concentric Lamb and Pekeris waves from the HTHH volcano. We estimated the source parameters: the initial amplitudes (amp0), the propagation speeds of the Lamb (C1) and Pekeris (C2) waves, and the origin time. In the previous studies, the Lamb wave speed and the origin time were assumed to be 300 m/s and 4:00 UTC (Kubo et al., 2022; Kubota et al., 2022), 310 m/s (Yamada et al., 2022), 317 m/s (Gusman et al., 2022), or revealed by the satellite observations to be 310 m/s (Otsuka, 2022) or 315 m/s (Shinbori et al., 2023). Our inversion of the tsunami waveforms indicated that the best C1 is 310 m/s, and the origin time of the propagating concentric Lamb wave was at 4:16 UTC. The origin time of the tsunami was estimated to be ~ 4:15 UTC from the nearby tide gauge record at Nuku’alofa on Tongatapu Island and local eyewitness accounts (Borrero et al., 2023). The USGS also estimated the eruption time to be 4:14:45 UTC, the same as the origin time of an M5.8 earthquake. The source origin time of 4:16 UTC is consistent with the earthquake origin time and the main eruption times revealed by other studies (Astafyeva et al., 2022; Diaz, 2022; Donner et al., 2023; Vergoz et al., 2022). These suggest that the eruption simultaneously generated atmospheric and seismic waves.

In the tsunami simulations to reproduce the later phases, we adopted the variation of C2 for each DART station. Our tsunami modeling showed that the C2 was faster (235–250 m/s) toward Japan, Kamchatka, and Alaska, while it was slower (200–210 m/s) toward North America and South America. These variations of the C2s suggested by the tsunami modeling are consistent with direct atmospheric observations and simulations. Watanabe et al. (2022) indicated that the speed of the Pekeris wave was 244 m/s at Yokohama, Japan, and 247 m/s at Honolulu, while it is slower (215 m/s) for the eastward propagation.

In the previous tsunami modeling studies (Gusman et al., 2022; Lynett et al., 2022), time functions of the pressure changes due to the Lamb wave were assumed to consist of positive and negative amplitudes. In this study, only positive amplitudes (amp) of pressure changes due to the Lamb wave and only negative amplitudes due to the Pekeris wave were assumed. Gusman et al. (2022) and Lynett et al. (2022) also modeled a localized tsunami source at the HTHH volcano, which was assumed to be a caldera collapse or other effects, in addition to the Lamb wave. The additional tsunami source may cause differences in the synthetic waveforms around the expected arrival time of the tsunami from the source. However, as the results of this study show, the amplification effect (Proudman resonance) of the Pekeris wave, whose propagation speed is close to that of a tsunami (long wave), is more significant in explaining the amplitudes of the later phases, especially at distant DART and tide gauge stations. In this study, the time scale and shape of propagating concentric pressure changes due to the Lamb and Pekeris waves were assumed to be the same, but they may also be sensitive to the results.

In the future, the tsunami calculation method utilized in this study has the potential to predict tsunamis resulting from volcanic eruptions, such as this 2022 HTHH event. Our method can be readily integrated into the input part of the conventional linear long-wave tsunami calculation code, and the parameters are straightforward. To provide tsunami early warning, tsunami waveforms observed near the source can be inverted using the tFISH method (Tsushima et al., 2009) or utilized for data assimilation methods (Wang et al., 2021) to forecast far-field tsunamis. This can be accomplished if pre-calculated tsunami waveforms, or Green’s functions, are prepared in advance and stored in a database. For instance, to apply the tFISH methodology, several datasets of Green’s functions can be prepared by presuming a 10-min rise time and several propagation speeds (200–315 m/s) of atmospheric waves, given a unit amount of initial amplitude at a known position of a target volcano. Then, assuming the linearity of a tsunami, we can estimate the initial amplitude (amp0) at the volcano by analyzing waveforms near the source. This enables us to forecast far-field tsunamis by multiplying the amp0 to the Green’s functions. However, predicting the later phases with slower propagation speeds may require waiting for data from offshore stations located close to the coastal area of concern. It is challenging to detect and invert the later phases observed near the source because, as evidenced in this case of the HTHH volcano, the amplitudes of the later phases do not increase significantly at the near field to the source. Additionally, the differences between the synthetic initial and later phases were small at the near-field stations, as shown by the black and blue curves in Fig. 3 for the New Zealand DARTs.

6 Conclusions

We simulated the tsunami generated by the Hunga Tonga–Hunga Ha’apai volcano eruption on January 15, 2022, adopting two propagating concentric waves corresponding to the Lamb and Pekeris waves. The tsunami waveform inversions of the initial phases observed at DART stations show the origin time of 4:16 UTC, the initial peak amplitude of 383 hPa at the volcano, and the Lamb wave propagation speed of 310 m/s, assuming the rise time of 10 min. The later phases of the tsunami can be better reproduced by adding another propagating concentric wave (Pekeris wave) with a negative initial amplitude (− 50 hPa) and propagation speeds of 200–250 m/s. The later phases recorded at DART stations around the Pacific indicate that the Pekeris wave speed is faster (235–250 m/s) toward the northwest, such as Japan, Kamchatka, and Alaska, while it is slightly slower (200–210 m/s) toward the northeast, North America, and the southeast, South America. The synthetic waveforms roughly reproduced not only the observed waveforms at DART stations but also the tsunami waveforms recorded on far-field tide gauges, including the later phases.