An Mw7.4 submarine earthquake occurred in the Kermadec Trench, northeast of New Zealand, on 18 June, 2020. This powerful earthquake triggered energetic tertiary waves (T-waves) that propagated through the South Pacific Ocean into the South Atlantic Ocean, where the T-waves were recorded by a hydrophone station near Ascension Island, 15 127 km away from the epicenter. Different T-wave arrivals were identified during the earthquake period with arrival angles deviating from the geodesic path. A three-dimensional sound propagation model has been utilized to investigate the cause of the deviation and confirm the horizontal diffraction of the T-waves at the Drake Passage.
1. Introduction
Submarine seismic activity can generate low frequency (<40 Hz) acoustic waves that propagate in the ocean with low attenuation for thousands of kilometers (Okal, 2008). These acoustic waves travel at the speed of sound in water (approximately 1.5 km/s), which is slower than P, S seismic waves, so they are known as tertiary waves (T-waves or T-phases). These waves like other low-frequency underwater acoustic waves can be profoundly influenced by bathymetric features, which can result in horizontal reflection, refraction, and diffraction (Doolittle et al., 1988; Heaney et al., 2017; Stephen et al., 2019; Oliveira and Lin, 2019). Hence, T-wave long-range propagation in the ocean can become complex, and three-dimensional (3D) underwater acoustic models are necessary for investigating its physics.
On 18 June, 2020, energetic underwater acoustic T-wave signals were recorded at the Comprehensive Nuclear-Test-Ban Treaty (CTBT) International Monitoring System (IMS) hydrophone station HA10 at Ascension Island in the South Atlantic Ocean. These T-waves were associated with a moment magnitude (Mw) 7.4 submarine earthquake with an epicenter in the Kermadec Trench, Pacific Ocean, located approximately 15 127 km from HA10 (see Fig. 1). In this work, the long-range propagation of these T-waves is investigated. First, the recorded T-wave signals at HA10 are analyzed using the Progressive Multi-Channel Correlation algorithm (DTK-GPMCC; Cansi, 1995) installed on the CTBT Organization (CTBTO) virtual Data Exploitation Centre (vDEC). Second, the long-range propagation of the T-waves from the announced epicenter to HA10 is investigated using a 3D parabolic equation (PE) model (Lin et al., 2013). This study in fact confirms the bathymetric diffraction effect at the Drake Passage suggested by Heaney et al. (2017) and the frequency dependency of this effect.
(a) Top panel: Location map of the 18 June 2020 Mw7.4 Kermadec Trench earthquake epicenter and the geodesic path (dashed white line) to the CTBT IMS H10N array (Ascension Island). Red lines present the back-azimuth paths of the T-wave arrivals detected at H10N. (b) Middle panel: Uncalibrated spectrogram of the T-wave signals at the H10N1 hydrophone. (c) Bottom panel: Time (UTC) and angles of the arrivals detected at the H10N array, and the signal waveform recorded at the H10N1 hydrophone.
(a) Top panel: Location map of the 18 June 2020 Mw7.4 Kermadec Trench earthquake epicenter and the geodesic path (dashed white line) to the CTBT IMS H10N array (Ascension Island). Red lines present the back-azimuth paths of the T-wave arrivals detected at H10N. (b) Middle panel: Uncalibrated spectrogram of the T-wave signals at the H10N1 hydrophone. (c) Bottom panel: Time (UTC) and angles of the arrivals detected at the H10N array, and the signal waveform recorded at the H10N1 hydrophone.
2. The 18 June 2020 MW7.4 South of the Kermadec Islands earthquake
As shown in the United States Geological Survey (USGS) earthquake database, an Mw 7.4 submarine earthquake occurred on 18 June, 2020, at 12:49:53.70 (hereafter all times in UTC + 00:00) in the Kermadec Trench (in the South Pacific Ocean) with hypocenter at 33.293°S 177.857°W (uncertainty ± 10 km) and depth 10.0 km (uncertainty ± 1.7 km) (USGS, 2021). The Kermadec Trench is a linear trench about 1000 km long and is one of the deepest oceanic trenches, with maximum water depths of approximately 10 km. The hypocenter was located approximately 445 km south of the Kermadec Islands. The International Seismological Centre (ISC) identified the earthquake as Event 618550810.
3. CTBT-IMS Hydrophone data
The CTBTO has been recording underwater signals using eleven hydroacoustic stations (five T-wave seismometer stations and six hydrophone stations) around the global ocean. T-wave signals from the 2020 Mw7.4 Kermadec Trench earthquake were recorded at the CTBT-IMS hydroacoustic hydrophone stations in the Pacific (HA11 and HA03) and South Atlantic (HA10) Ocean. In this work, the focus is on the analysis of the T-wave signals recorded at the HA10 station. A CTBTO's hydroacoustic hydrophone station consists of six hydrophones, suspended in the sound channel, arranged into a pair of triangular arrays except for the station at Cape Leeuwin where there are only three hydrophones. Each array has an element spacing of approximately 2 km. At the HA10 station, the two arrays (H10N and H10S) are situated on opposite sides (north and south) of Ascension Island to account for the acoustic shadow produced by the island, and the southern array H10S regularly detect T-waves from the Tonga and Kermadec Trenches (South Pacific) as shown by Heaney et al. (2017). Note that one hydrophone of the south array (H10S) during the Kermadec Trench earthquake was not recording, and data from the other two hydrophones on this array are not analyzed in this work. Only the north array (H10N) data is analyzed.
The recorded T-wave signals at H10N were first processed using the DTK-GPMCC algorithm (Cansi, 1995). A strong signal correlation was observed between the signals at all of the three H10N hydrophones (H10N1, H10N2, and H10N3), which are placed between 845 and 850 m close to the sound-channel depth as a triangular array. The center of the H10N triangular array is located at 7.838S 14.490 W, and its geodesic distance and azimuth to the announced Kermadec Trench earthquake epicenter are 15 127 km and 200.01 degrees due North (°N). Note that the Antarctic peninsula blocked the geodesic path between the announced epicenter and H10N hydrophones (see the top panel of Fig. 1), and it will be shown later in the paper that the bathymetric diffraction effect occurred at the Drake Passage and enabled the T-waves to reach to the H10N array.
On 18 June, 2020, between 15:30:00 and 15:50:00 at H10N, sound energy in the frequency range of 3 to 30 Hz was dominated by the earthquake T-waves (see the middle panel of Fig. 1). Different arrivals (see the bottom panel of Fig. 1) were identified close to the theoretical T-waves arrival time at 15:37:58, which was calculated considering water sound speed as 1500 m/s and the geodesic distance between the announced epicenter and the center of the H10N triangular array regardless the blockage of the Antarctic peninsula. Thirteen arrivals were detected between 15:30:00 and 15:50:00, with arrival azimuths varying from 203.6 to 205.6 °N (see the bottom panel of Fig. 1) and peak frequencies of detected signals ranging between 6 to 16 Hz. DTK-GPMCC algorithm was used to compute the detections in the frequency band of 2 to 40 Hz. The strongest signals were recorded between 15:39:00 and 15:47:00, with six arrivals detected in this period with arrival azimuths ranging between 203.6 and 204.4 °N.
The coupling of the in-ground earthquake to a propagating acoustic wave in the water is a very complex process. In fact, T-waves can be generated at different locations away from the epicenter (Okal, 2008). Hence, the deviation of the detected T-wave arrival angles at the H10N array could be due to T-wave generation away from the epicenter. However, the 3D propagation analysis presented in Sec. 4 will demonstrate stronger evidence showing that the bathymetric diffraction at the Drake Passage was more likely to be the mechanism enabling the T-waves to arrive at the H10N array.
4. Numerical simulations
A 3D sound propagation model is utilized based on the PE approximation and the split-step Fourier (SSF) method (Lin et al., 2013). A multi-core parallel and General-Purpose Graphics Processing Unit (GPU) version of the PE models are employed (Kushida et al., 2020). A geodesic Cartesian implementation of this 3D PE model (Oliveira and Lin, 2019) is performed to consider Earth's curvature. The properties of the bottom considered in the simulations are sound speed = 1700 m/s, density = 1500 kg/m3 and attenuation = 0.5 dB/λ, where λ is the wavelength. Also, for all the simulations, water density is considered constant of ρw = 1000 kg/m3. A bathymetric model was constructed using the 15 arc-second resolution (∼450 m) SRTM15+ global bathymetry and topography database (Tozer et al., 2019). Water salinity and temperature outputs from the∼9 km horizontal resolution global HYCOM-NCODA (Cummings, 2005) hindcast model for 13:00:00 on 18 June, 2020, was used to compute the 3D sound speed field using the equation proposed by Mackenzie (1981).
Assuming a 5-Hz T-wave source near the seafloor at the epicenter (depth 7700 m), the 3D propagation model result is presented in the left panels of Fig. 2. Frequencies below 5 Hz are not simulated as they could lead to unrealistic results due to the fact that the 3D PE model does not consider the elasticity of the bottom (Oliveira and Lin, 2019). The size of the simulation domain was 15 127 km × 2000 km × 11 km (range × cross-range × depth) and solved with a model grid where each element was 300 m × 150 m × 4 m. A range discretization smaller or equal to one wavelength, a cross-range discretization smaller than half the discretization used in range, and a vertical discretization of approximately lower or equal than one-quarter of wavelength should be used in the 3D PE model used in this work. A numerical test has been conducted to ensure that the 3D PE solution is converged. An image field was imposed to satisfy the pressure-release boundary condition at the sea surface, and artificial absorbing layers were placed around the computational domain in the vertical and transverse direction to mimic radiation boundary conditions. At each PE marching solution step, the sound pressure field in the cross-range and depth plane was computed on an NVIDIA® V100 workstation that equips with 32 giga-bytes (GB) of memory. Comparing with the azimuthally independent two-dimensional (Nx2D) model results (the right panels of Fig. 2), the 3D model results of depth-integrated sound intensity (assuming 0-dB source level) and transmission loss (TL) along the geodesic path indeed show a strong horizontal diffraction effect at the Drake Passage caused by the surrounding bathymetry. Strong focusing and shadow zones can also be observed in the 3D model but cannot be reproduced by the Nx2D simulations without horizontal coupling. The back azimuth paths from the H10N array also confirm that the detected arrivals were indeed following the bathymetric diffracted path originated at the Drake Passage, indicating a strong evidence of the horizontal diffraction effect.
3D (left) vs Nx2D (right) sound propagation model results between the epicenter and H10N for a 5 Hz source. Top panels: Depth-integrated sound intensity from a 0 dB source. Black lines present the back-azimuth paths of the T-wave arrivals detected at H10N. Green line is the geodesic path between the announced epicenter and H10N. Bottom panels: TL levels along the geodesic path from the epicenter to H10N. Black line indicates the bathymetry.
3D (left) vs Nx2D (right) sound propagation model results between the epicenter and H10N for a 5 Hz source. Top panels: Depth-integrated sound intensity from a 0 dB source. Black lines present the back-azimuth paths of the T-wave arrivals detected at H10N. Green line is the geodesic path between the announced epicenter and H10N. Bottom panels: TL levels along the geodesic path from the epicenter to H10N. Black line indicates the bathymetry.
Frequency dependency of the T-waves propagation between 5 to 15 Hz was also investigated and shown in Fig. 3. The reason for selecting this frequency band was because the T-waves were more energetic in this band at the H10N array (see Fig. 1). The same simulation domain used in the 5-Hz simulations but with a finer grid resolution (50 m×50 m×4 m) was utilized in the simulations. At each PE marching step, the pressure field on a cross-range and depth plane was discretized into a 39 996 × 4500 matrix and computed on two of the 14-core Intel® Xeon® Gold 6132 2.60 GHz CPUs and 196 GB memory. Although the computation speed of the GPU is faster (Kushida et al., 2020), the 32 GB memory limited 3D PE model computation to a 20 000 × 4500 PE marching solution matrix. Hence multi-core CPUs with larger memory size were utilized. The multi-frequency PE model results revealed that the depth-integrated sound intensity over the water column (assuming 0-dB source level) at H10N array varies from –141.7 dB to –142.9 dB (see Table 1 for details). In addition, beamforming computations considering the pressure field at 847 m depth from the 3D PE model were performed. As shown in Table 1, the predicted arrival azimuths within the considered frequency bandwidth at H10N range between 203.7 and 204.6 °N, which in fact agree with the observation shown in Fig. 1 very well.
3D sound propagation results between the epicenter and H10N for a 10 and 15 Hz source. Left panel: Depth-integrated sound intensity from a 0 dB source. Green line is the geodesic path between the announced epicenter and H10N. Right panel: TL levels along the geodesic path from the epicenter to H10N. Black line indicates the bathymetry.
3D sound propagation results between the epicenter and H10N for a 10 and 15 Hz source. Left panel: Depth-integrated sound intensity from a 0 dB source. Green line is the geodesic path between the announced epicenter and H10N. Right panel: TL levels along the geodesic path from the epicenter to H10N. Black line indicates the bathymetry.
3D sound propagation results at H10N. Depth-integrated sound intensity over the water column from a 0 dB source and arrival azimuth at 847 m depth.
Frequency (Hz) . | Depth-integrated sound intensity (dB) . | Arrival azimuth (°N) . |
---|---|---|
5 | −142.9 | 203.8 |
10 | −142.6 | 204.6 |
15 | −141.7 | 203.7 |
Frequency (Hz) . | Depth-integrated sound intensity (dB) . | Arrival azimuth (°N) . |
---|---|---|
5 | −142.9 | 203.8 |
10 | −142.6 | 204.6 |
15 | −141.7 | 203.7 |
5. Conclusions
T-waves from the 2020 Mw7.4 Kermadec Trench earthquake were recorded 15 127 km away from its announced epicenter at the CTBT IMS HA10 hydroacoustic hydrophone station. Different arrivals with azimuths varying from 203.6 to 205.6 °N were identified at the H10N array within a 20 min time window during the earthquake event. Numerical acoustic wave propagation simulations have shown that, although the Antarctic peninsula blocked the geodesic path between the epicenter and the HA10 hydrophones, the T-waves generated at the epicenter could still reach the HA10 station after being diffracted by the complicated bathymetry at the Drake Passage. Frequency dependency of this bathymetric diffraction effect has also been investigated. Modeled arrival azimuths at the H10N array for 5 to 15 Hz vary between 203.7 and 204.6 °N, which agree with the measured arrival azimuths (203.6 to 204.4) of the strongest T-waves detected at the hydrophone station. Besides that, numerical results also indicated the frequency dependency of T-wave intensity at H10N. However, due to the lack of information on T-wave source level, model comparison with field measurements can only be based on arrival azimuths. Overall, the presented numerical results confirm the importance of 3D effects on low-frequency sound passing the Drake Passage previously reported by Heaney et al. (2017) and also shed light on the frequency dependency of the phenomena. For future research directions, the authors would like to investigate the influence of geoacoustic properties in the bottom on the global scale T-wave propagation in detail, as well as T-wave generation mechanism and horizontal refraction due to sound speed gradients in the water column.
The numerical modeling shown in the paper also reveals the importance of accelerated high performance computing (HPC) to accurately predict 3D propagation of T-waves in complex large-scale ocean environments. Numerical simulations require an excessive computational domain to correctly evaluate the reflection, refraction and diffraction caused by oceanic features. In general, parallel computation with GPUs can be much faster than CPU computation. However, the GPU memory still limits the size of the computation domain. In these cases, multi-core CPU computation is a better option to perform simulations.
Acknowledgments
T.C.A.O. thanks FCT/MCTES for the financial support to CESAM (UIDP/50017/2020+UIDB/50017/2020), through national funds. Y.-T.L. would like to thank the Office of Naval Research for their support through Grant No. N00014-21-1-2416. Luís Carvalheiro is acknowledged for his help on computer code compilation on ARGUS-HPC at the University of Aveiro (UA). The authors would like to thank the Preparatory Commission for the Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO) for providing the limited access to the hydroacoustic data under the scientific project “3D long-range underwater sound modeling” (Contract No. vDEC2020-0146, contractor UA, PI: T.C.A.O.). The views expressed in this study are those of the authors and do not necessarily reflect those of the Preparatory Commission for the CTBTO.