Characteristics of two tsunamis generated by successive Mw 7.4 and Mw 8.1 earthquakes in Kermadec Islands on March 4, 2021

On March 4, 2021, two tsunamigenic earthquakes (Mw 7.4 and Mw 8.1) occurred successively within 2 h in 10 Kermadec Islands. We examined sea level records at tide gauges located at ~100 km to ~ 2,000 km from the epicenters, conducted Fourier and Wavelet analyses as well as numerical modelling of both tsunamis. Fourier analyses indicated that the energy of the first tsunami is mainly distributed in the period range of 5–17 min, whereas it is 8–28 min for the second tsunami. Wavelet plots showed that the oscillation of the first tsunami continued even after the arrival of the second tsunami. As the epicenters of two earthquakes are close (~ 55 km), we reconstructed the source spectrum of the second tsunami by 15 using the first tsunami as the empirical Green’s function. The main spectral peaks are 25.6 min, 16.0 min, and 9.8 min. The results are similar to those calculated using tsunami/background ratio method and also consistent with source models.

successive tsunamigenic earthquakes are very close to each other (~55 km; Figure 1) and their focal mechanisms are similar; both are thrust earthquakes. The National Emergency Management Agency, New Zealand, issued a tsunami warning for coastal areas of the North Island after the Mw 8.1 earthquake. The tsunami propagated across the Pacific Ocean and reached South America. No casualties were reported. The situation of these two consequent tsunamigenic earthquakes resembles the earthquake events also in Kermadec Islands on January 14, 1976, where two great earthquakes (Mw 7.8 and Mw 8.0) occurred 35 approximately within one hour (Power et al., 2012). The mainshock of the 1976 events (Mw 8.0) generated a moderate tsunami recorded by tide gauges.
The occurrence of two successive earthquakes provides us with a rare opportunity to study their source characteristics by different methods. The spectra of tsunami waveforms at tide gauges contain the effects of source, propagation path, and local topography (Rabinovich, 1997;Heidarzadeh et al., 2016;Corté s et al., 2017). For a single event, it is a common 40 practice to reconstruct tsunami source spectrum from tide gauge records by calculating the ratio of tsunami spectra to background spectra (Rabinovich, 1997). This is based on the assumption that the effects of transmission path on the characteristics of tsunami spectra is small compared to local topographic effects when the tsunami source is not too far away from the observation station (Miller, 1972;Rabinovich, 1997). This assumption can be readily made for two successive tsunami events in Kermadec Islands. If the sources of two earthquakes are close to each other, we can adopt the method of 45 empirical Green's function (EGF) to reconstruct tsunami source spectrum (Heidarzadeh et al., 2016). The smaller tsunami event is adopted as the EGF for the larger one. The spectral deconvolution separates the effects of propagation path and local topography around the tide gauge and gives the source spectrum of the larger tsunami. Heidarzadeh et al. (2016) successfully applied the EGF method to the 2015 (Mw 7.0) and 2013 (Mw 8.0) earthquakes in the Solomon Islands.
Here, we studied the characteristics of tsunamis generated by the two successive Kermadec Islands earthquakes and 50 calculated the source spectrum of the tsunami generated by the mainshock (hereafter, the second tsunami). Fourier analysis was applied to the sea level records at 15 tide gauges. Wavelet analysis was adopted to examine the temporal changes of the dominant spectral peaks. Finally, we used two different methods, tsunami/background ratio method and EGF method, to reconstruct tsunami source spectra. We compared two alternative methods of determining tsunami source spectra and also compared their results with USGS source models. This is a unique and rare incident and thus the data and analyses would 55 greatly help to further understand tsunami generation and propagation.

Tsunami Data
We collected sea level records of 15 tide gauges from the Intergovernmental Oceanographic Commission (http://www.ioc-sealevelmonitoring.org/list.php). These stations are located at distances between ~100 km and ~ 2,000 km 65 from the epicenters (Figure 1). The sampling rate is 60 s. First, we conducted quality control of these data, removing spikes and filling short gaps by linear interpolation. Then, we applied a second-order high-pass filter with the corner frequency of 0.00014 Hz (7,200 s) to remove the tidal components ( Figure 2) (Heidarzadeh and Satake, 2013). The time series between 17:00:00 (UTC), March 4 and 05:00:00 (UTC), March 5 was selected for analysis. However, as the Raoul Island Fishing Rock (RIFR) tide gauge ( Figure 1) lost its data communication during the mainshock, we only used its records before 70 17:41:00 (UTC) (Figure 2).

Spectral Analyses
Two types of spectral analyses were conducted to tide gauge records: Fourier analysis and wavelet (frequency-time) analysis. For Fourier analysis, a Hanning window with 40% overlaps was applied following the method of Heidarzadeh and Mulia (2021). The tsunami generated by the foreshock (hereafter, the first tsunami) was clearly recorded at eight tide gauges: 75 North Cape (NC), Great Barrier Island (GBI), East Cape (EC), Suva Viti Levu (SVL), Lenakel, Quinne, Kingston Jetty Norfolk (KJN), and RIFR, but it is difficult to accurately distinguish the arrival times at these stations. We used the 2 h segment before the arrival of the second tsunami as the input for spectral analyses. For RIFR, we utilized the 2 h segment before the time it lost connection. The second tsunami generated by the mainshock was clearly recorded at all the tide gauges except for RIFR. We selected a 2 h time segment after the tsunami arrival as the input for Fourier analysis. Additionally, 80 background spectra were calculated using the 2 h records before the arrival time of the first tsunami at the same stations. We ensured that there were no storms or other atmospheric events at the time period of the background signals, so the background spectra could exclusively reflect the frequency response of local topography. Tidal components were removed by applying a high-pass filter in a similar way to preparation of the tsunami records.
Wavelet analysis was adopted to study the temporal changes of the dominant spectral peaks and the superposition of 85 two tsunamis. We applied the wavelet package provided by Torrence and Compo (1998), which uses the Morlet wavelet mother function. We only conducted wavelet analysis to seven tide gauges that clearly recorded both tsunamis. The waveform segments for wavelet analyses are the same as those used for Fourier analyses.

Earthquake Slip Models and Tsunami Numerical Simulation
We used numerical simulation to study the propagation paths of the two tsunamis. For the initial condition, the finite 90 fault models provided by the USGS were adopted (Slip model of the first event at: https://earthquake.usgs.gov/earthquakes/eventpage/us7000dfk3/finite-fault; Slip model of the second event at: https://earthquake.usgs.gov/earthquakes/eventpage/us7000dflf/finite-fault). The source model of the first (Mw 7.4) earthquake has a rectangular dimension of 80 km (length) × 80 km (width). The strike and dip angles are 196° and 32°, respectively. The spatial intervals of sub-faults are 4 km × 4 km. The source model of the second (Mw 8.1) earthquake has a 95 rectangular dimension of 240 km (length) × 190 km (width) with strike and dip angles of 210° and 16°, respectively. The spatial intervals of sub-faults are 10 km × 10 km. We used Okada's dislocation model (Okada, 1985) to compute the seafloor deformation and used it as the initial condition for tsunami propagation modeling. The bathymetric grids were obtained from the General Bathymetric Chart of the Oceans and resampled to 0.9 arc-min. The JAGURS numerical package (Baba et al., 2015) was adopted to simulate the tsunami propagation based on linear long-wave equations. The time step for 100 numerical simulations is 1.0 s. We simulated two tsunamis separately up to 12 h after the earthquake origin time. Tsunami travel time (TTT) calculation was performed for the second tsunami by the TTT software of Geoware (2011).   (Rabinovich and Thomson, 2007) and the 2020 Aegean Sea tsunami . The largest tsunami amplitude (0.64 m) was recorded by KJN, which is also located in 120 the direction parallel to the short axis of the fault. Besides, the second tsunami has a longer wavelength than the first tsunami, leading to longer tsunami periods. At stations without clear signal of the first tsunami (e.g., Lautoka), the tsunami waveforms are dominant by the long-period components which are generated by the second tsunami. However, the stations that recorded both events show the superposition of two successive tsunamis. For example, the waveform of the second tsunami at EC is mixed up with the short-period components in the few hours after its arrival (21:00-24:00), likely caused by the oscillation 125 of the first tsunami ( Figure 2). After 24:00, the tsunami waveforms mainly present long-period components, possibly due to the dissipation of the short-period waves of the first tsunami. Similar patterns were also be observed at Ouinne (Figure 2).

Spectral Analyses
Figures 3 presents the results of Fourier analyses for two tsunamis and background signals. The gap between the tsunami and background spectra is attributed to the tsunami energy. The periods of the tsunami spectral peaks generally 130 contain the effects of tsunami source, propagation path, and local topography. According to Figure 3, the second tsunami has larger energy (spectral power) than the first tsunami due to the larger magnitude of the second earthquake. The background spectra are smoother than tsunami spectra at most stations. At most stations, the peak periods of the first tsunami are distributed in the range of 5-17 min, whereas the dominant periods range for the second tsunami is approximately 8-28 min.
In other words, the second tsunami has generated longer-period waves, which is natural as the source dimension of the 135 second earthquake (240 km × 190 km) is much larger than that of the first earthquake (80 km × 80 km).
Wavelet analyses reveal the variations of dominant tsunami peak periods over time (Figure 4). At EC, Lenakel, and KJN, the arrivals of two successive tsunamis can be clearly identified on the wavelet plots. The second tsunami has larger energy levels and longer-period waves than the first tsunami (Figure 4)

Reconstructing the Tsunami Source Spectrum
A tsunami source spectrum reveals the earthquake source characteristics without the effects of tsunami propagation path or local topography. In our study, the epicenters of two earthquakes are close to each other (~55 km; Figure 1); the earthquakes are of similar mechanism (both thrust events), and the propagation paths are similar (Figures 5a and 5b). Hence, it enables us to reconstruct tsunami source spectrum using EGF method which assumes that the smaller event acts as the 170 EGF for the larger event (Heidarzadeh et al., 2016). We computed the spectral ratio of the second tsunami to the first tsunami at seven tide gauges (Figure 5c; gray curves), and then calculated their normalized average values (Figure 5c; blue curve) through adjusting the peak energy of all spectra to that of the largest one. The source spectrum shows that the energy of the second tsunami is mainly distributed in the period range of 8-30 min, with spectral peaks at 25.6 min, 16.0 min, and 9.8 min. The period range is generally consistent with the results of Fourier analyses of each station for the second event 175 ( Figure 3). Figure 5d shows the results of the method proposed by Rabinovich (1997) which is based on dividing the tsunami spectra to that of the background to construct tsunami source spectrum. We computed the spectral ratio of the second tsunami to the background signals at all tide gauges except RIFR and calculated their normalized average. The period range of main energy (7-28 min) contains spectral peaks at 25.6 min, 16.0 min and 8.5 min, which are close to the spectral peaks calculated by EGF method. 180 Theoretically, the tsunami source period is related to earthquake rupture length and water depth (Rabinovich, 1997;Heidarzadeh and Satake, 2013;Wang et al., 2021). It can be estimated as: where is the typical size of tsunami source area (length or width), is the gravity acceleration, and ℎ is the average water depth in source area. According to the USGS model, the size of the Mw 8.1 earthquake is 240 km (length of the fault) × 190 185 km (width of the fault). The average water depth in the source area is ~ 5,000 m. Hence, the first three source periods of the short axis (width) are 28.6 min, 14.3 min, and 9.5 min. The first three source periods of the long axis (length) are 36.1 min, 18.1 min, and 12.0 min. The source periods of the short axis match well with the peaks of tsunami source spectrum, which confirms the validity of the model geometry. We note that the periods of tsunami waves are mainly influenced by the short axis (width) (Heidarzadeh and Satake, 2013). 190 In addition, the results of two methods (i.e., EGF and tsunami/background spectral ratio) show similarities in shapes and peak periods of tsunami source spectrum (Figures 5c and 5d). It is noted that the EGF method has the capability to remove both propagation-path effects and local bathymetric effects from the tide gauge records whereas the https://doi.org/10.5194/nhess-2021-369 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License. tsunami/background spectral ratio method would remove mainly local bathymetric effects. Hence, our results may imply that the effects of propagation path are negligibly small for this case, where the tide gauges are located at distances between ~100 195 km and ~ 2,000 km from the source. Both of the methods are able to reveal the source characteristics merely based on tsunami observations rather than seismological data.
https://doi.org/10.5194/nhess-2021-369 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License. to that of the first tsunami (EGF method). Blue curve is the normalized average of tsunami spectral energy at different tide gauges. (d) Spectral ratio of the second tsunami spectrum to the background signal spectrum. Green curve is the normalized average of different tide gauges.

Conclusion 205
We studied the characteristics of tsunamis generated by two earthquakes (Mw 7.4 and Mw 8.1) that occurred in Kermadec Islands successively on March 4, 2021, within ~55 km from each other and within approximately 2 h interval. We used the sea level records of 15 tide gauges. The spectra of Fourier analyses show that the dominant period bands of the first tsunami and the second tsunami are 5-17 min and 8-28 min, respectively. Two oscillation patterns with different period ranges are visible on the wavelet plot at most stations which belong to the first and the second tsunamis. We observed that 210 after the arrival of the second tsunami, the oscillation in the period range of the first tsunami still exists. We calculated the tsunami source spectrum of the larger event (i.e., the second tsunami) by two approaches: empirical Green's function (EGF) method and tsunami/background ratio method. Using the first tsunami as the EGF, spectral deconvolution indicated that energy of the second tsunami is mainly distributed in the period range of 8-30 min, with spectral peaks at 25.6 min, 16.0 min, and 9.8 min. The method of tsunami/background ratio yielded similar results to the EGF method. The source characteristics 215 were obtained merely based on tsunami data and thus these two methods could be complementary to seismological approaches in source analysis.

Data availability
The datasets were derived from sources in the public domain. The tide gauge data used in this research were provided by the Sea Level Station Monitoring Facility of the Intergovernmental Oceanographic Commission (http://www.ioc-220 sealevelmonitoring.org/list.php). We used bathymetric and topographic data of the General Bathymetric Chart of the Ocean (https://www.gebco.net/data_and_products/gridded_bathymetry_data/). We used earthquake focal mechanism catalog of the Global Centroid Moment Tensor Project (https://www.globalcmt.org/CMTsearch.html) and earthquake source model of the United States Geological Survey (https://earthquake.usgs.gov/earthquakes/eventpage/us7000dflf/executive). We used the GMT software for drafting some of the figures (Wessel and Smith, 1998). 225

Author contribution
YW is responsible for methodology, data curation, software, and writing -original draft preparation. MH is responsible for methodology, software, and writing-reviewing and editing. KS is responsible for methodology, supervision, and writingreviewing and editing. GH is responsible for data curation and software. https://doi.org/10.5194/nhess-2021-369 Preprint. Discussion started: 2 December 2021 c Author(s) 2021. CC BY 4.0 License.