Three-dimensional global scale underwater sound modeling: The T-phase wave propagation of a Southern Mid-Atlantic Ridge earthquake

Underwater infrasound (f < 20 Hz) can be generated in the ocean by mechanical energy transfer from the Earth's crust (e.g., earthquakes or volcanoes) and by energy transfer occurring at the water surface (e.g., the interaction of opposing gravity waves, ice-quakes, or localized pressure changes), and it can travel thousands of kilometers in the ocean. The infrasound generated by submarine earthquakes is known as T-phases or T wave (tertiary phase wave) (Ewing et al., 1952; Okal, 2008), and it has been successfully applied to estimate submarine earthquake rupture model parameters (Okal, 2003; de Groot-Hedlin, 2005; Tolstoy and Bohnenstiehl, 2005). In this regard, it is considered that T wave energy is lower for normal and reverse fault events than strike-slip events (Dziak, 2001). T waves travel at the speed of sound in water (approximately 1.5 km/s) and slower than P, S seismic waves and Rayleigh surface waves which propagate at approximately 5.0–8.0, 3.0–4.8, and 3.0–4.0 km/s, respectively. Compared to seismic and surface waves, T waves correspond to the high-frequency part of the earthquake source spectra. In this study, a three-dimensional (3D) geodesic Cartesian parabolic equation model is utilized to study the global propagation of T-phase waves.

As low-frequency sound generated by submarine earthquakes travels much faster than tsunamis, it was anticipated that appropriate measurement and analysis of underwater infrasound triggered by ocean bottom movements would enhance current early tsunami detection systems (Abdolali et al., 2015; Oliveira and Kadri, 2016). However, the reality is more complicated than that since the infrasound propagation to a certain distance from an epicenter can be influenced by a number of physical oceanography and marine geological effects. These propagation effects can plausibly explain the decorrelation between T-phases amplitudes and tsunami generation reported by several studies investigating the relation between tsunamis and T waves triggered by tsunami and tsunamigenic earthquakes (Okal, 2008; Ewing et al., 1950; Hiyoshi et al., 1992; Walker and Bernard, 1993; Okal, 2003; Okal et al., 2003).

The understanding of long-range sound propagation plays a vital role in analyzing and characterizing sound from submarine earthquakes. In this regard, we must consider that underwater sound propagation can be profoundly influenced by three-dimensional (3D) effects (Fawcett, 1993; Luo and Henrik, 2009; Colin et al., 2013; Lin et al., 2015; Heaney and Campbell, 2016; Duda et al., 2011). A variety of geological and physical oceanographic features can cause horizontal refraction, reflection, and diffraction on long-range underwater sound propagation (Heaney et al., 2013; Heaney and Campbell, 2016; Heaney et al., 2017). In these cases, 3D underwater sound propagation models are required for accurately predicting the sound pressure field. Several 3D oceanic acoustic propagation models have been developed over the past decades (Kuperman et al., 1991; Luo and Henrik, 2009; Lin et al., 2015; Porter, 2016; Calazan and Rodriguez, 2017). However, solving the underwater sound propagation accurately for fully 3D environments involves significant scientific challenges and still leads to very high computational costs. A parabolic equation sound propagation model is utilized in this study, and this type of sound propagation model is focused on its computational efficiency and accuracy.

In this work, we aim to improve the understanding of long-range propagation from deep to shallow water of infrasound triggered by submarine earthquakes. More specifically, we aim to answer the following question: What are the 3D effects on long-range submarine earthquake sound propagation? Our study focuses on the high-frequency part of the earthquake source spectra (5–20 Hz), the T-phase wave spectrum, with a future research plan on localization of T-phase wave generation/excitation locations. In order to reach our goal we analyze and model the propagation of T waves generated by a moment magnitude (M_w) 4.8 Southern Mid-Atlantic Ridge earthquake that occurred on August 31, 2006, at approximately 1030 km South of Ascension Island. Data recorded off the coast of Ascension Island by the deep water hydrophones of the Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO) network and off the US east coast by the hydrophones of the Shallow Water 2006 (SW06) experiment (Tang et al., 2007; Newhall et al., 2007) are analyzed. The data were recorded at approximately 8655 and 1050 km from the epicenter of the earthquake and showed discrepancies between expected and measured earthquake sound arrival angles, which suggests the plausible occurrence of horizontal sound reflection and/or refraction in the long waveguide journey. The arrival angle discrepancies can also indicate that the T waves were likely not excited at the epicenter, so the arrival angles of the T waves at the hydrophones deviated from the predicted geodesic azimuths. Nonetheless, our study here is to focus on the 3D propagation effect through numerical modeling and examine the possibility of horizontal reflection and/or refraction causing the arrival angle deviation. There is indeed a likelihood that both of the effects co-exist and jointly contribute to what was observed. A 3D and N × 2D parabolic equations model (Lin, 2013; Lin et al., 2013), considering realistic environmental parameters, is utilized. An assessment of the capability and robustness of the numerical propagation is carried out.

This paper is composed of five sections. In Sec. II, we present the Southern Mid-Atlantic Ridge earthquake studied in this paper and the earthquake T-phase signal received on the CTBTO hydrophone station near Ascension Island and the hydrophone array deployed in the SW06 experiment. Next, we brief the three-dimensional parabolic-equation sound propagation and ocean environment models employed in this study. Numerical results are presented in Sec. IV. Finally, the paper is concluded in Sec. V.

A moment magnitude (M_w) 4.8 submarine earthquake occurred on August 31, 2006, at 05:44:03 UTC in the Southern Mid-Atlantic Ridge with hypocenter at 17.230°S 15.336°W and depth 10.0 km (United States Geological Survey, 2019). The epicenter was located at approximately 1030 km South of Ascension Island. The Mid-Atlantic Ridge, the most extensive chain of mountains on Earth, together with the rift valley along its crest constitute a complex deepwater bathymetry environment. Figure 1 shows the location of the epicenter and the geodesic paths to CTBTO and SW06 hydrophone stations, which also recorded the event. Both geodesic paths indicate that the earthquake T-phase signal should cross different parts of the Mid-Atlantic Ridge in its path to the hydrophone stations.

Measurements were acquired from a hydrophone station located near Ascension Island, installed and maintained by the CTBTO. The Ascension Island station is one of six hydrophone stations of the CTBTO International Monitoring System (IMS). The station consists of six hydrophones, moored in the sound channel at approximately 850 m depth, arranged into a pair of triangular arrays (see Fig. 1 and Table I). Each array has an element spacing of approximately 2 km. The two arrays are situated on opposite sides of the island, one to the north and one to the south of the island, in order to account for the acoustic shadow produced by the island. The Northern array (HA10N) is approximately 20 km from the island, whereas the southern array (HA10S) is approximately 110 km from the island. The signals recorded at each hydrophone are low-pass filtered and sampled at 250 Hz. Buried fiber optic cables carry the digitized signals to a station on Ascension Island, where they are transmitted via satellite to the CTBTO headquarters in Vienna for near real-time monitoring.

The Southern Mid-Atlantic Ridge earthquake on August 31, 2006 was detected at hydrophones H10S1 (Southern array) and H10N1 (Northern array) at 05:54:14-05:55:55 and 05:55:40-05:56:42 UTC, respectively. Figure 2 presents, for the earthquake T-phase signal recorded period, the spectrograms at each hydrophone of the Southern and Northern arrays. From these spectrograms, we can observe that energy in the 1–30 Hz frequency band decreases from the Southern to the Northern array indicating the sound shadow effect induced by the Island. The spectrograms also show the dispersive nature of propagation with the highest frequency part of the signal terminating before the lower frequencies.

The geodesic path between the epicenter and HA10N suggests that the arrival angle of the earthquake sound at the array should be 184.9° due North. However, a beamforming of horizontal line array signals showed that the arrival angle was 190° due North.

The SW06 experiment was conducted during the summer of 2006 on the edge of the continental shelf off New Jersey (Tang et al., 2007; Newhall et al., 2007). This large interdisciplinary experiment had two main components for studying acoustics and oceanography. An acoustic receiving system of vertical and horizontal hydrophone arrays was deployed during the experiment. The vertical array consisted of 16 hydrophones covering three-fourths of the water column (80 m deep) from 13 to 75 m in depth, and the horizontal array consisted of 32 hydrophones spaced at 15 m intervals and deployed on the bottom. The nominal array orientation was North to South, and a long baseline transponder system was deployed to monitor the array position (Newhall et al., 2007).

The Southern Mid-Atlantic Ridge earthquake of August 31, 2006 was detected at the vertical line array during the SW06 experiment at 07:20:15-07:21:45 UTC. Spectrograms for sound recorded at hydrophones channels 2, 8, and 15, at 21.00, 43.50, and 77.75 m water depth, respectively, are shown in Fig. 3. The geodesic path between the epicenter and the vertical array (see Fig. 1) suggests that the arrival angle of the earthquake sound at the array should be 124.33° due North. However, data recording during SW06 experiments indicate that the arrival angle must have been 122.2° due North. The possible causes of this 2° discrepancy could be related to the accuracy of array location and/or the occurrence of horizontal sound refraction/reflection.

A 3D sound propagation model using the parabolic equation (PE) approximation and the split-step Fourier (SSF) method in a Cartesian system (Lin and Duda, 2012) is applied in this study. The model solves the approximated Helmholtz wave equation by taking a square root of the propagation operator and utilizing the SSF solution marching scheme. In addition, this model uses a density-reduced pressure variable (Bergman, 1946) to handle the density variation across the seafloor interface. To maintain the accuracy of the square root Helmholtz operator, the model utilizes the wide-angle approximation proposed by Feit and Fleck (1978).

In order to identify 3D propagation effects that cannot be detected by 2D models, N × 2D (Perkins and Baer, 1982) numerical simulations are also carried out. This type of simulation uses a 2D model with the same PE technique used in the 3D simulations but limiting the transverse coupling of sound energy and only considering its variation in the radial direction.

For the application of global scale propagation, we need to map the Earth's surface onto the Cartesian coordinates of the PE model. This mapping consists of two steps. We first let the x axis of the PE grids (the solution marching axis) to align with the geodesic path between source and receiver. We then, at each marching step, determine the transverse geodesic grids (the y axis) along the cross-range perpendicular to the x axis. Since this model grid mapping is based on geodesics, we refer this model to “geodesic Cartesian PE.” Figure 4 shows the horizontal model grid utilized for sound propagation modeling from the epicenter to the SW06 site. As for the vertical axis, the flat-earth transformation (Aki and Richards, 2002) is applied to sound speed profiles and bathymetry to account for Earth's curvature.

A bathymetric model was constructed using the 1/60° resolution database ETOPO1 Global Relief Model from the National Geophysical Data Center (2008). Water salinity and temperature models, from which a 3D sound speed model was derived (Mackenzie, 1981), were characterized from the 1/12° resolution global HYCOM-NCODA (Cummings, 2005) hindcast model. Figure 5 shows salinity and temperature at 100 m depth given by the HYCOM-NCODA hindcast model for August 31, 2006. The flat-earth transformation (Aki and Richards, 2002) was applied to the sound speed field and water depth models in order to account for Earth's curvature. Essentially, this consists of changing the sound speed based on the latitude, longitude and depth of each grid point leading to an increase of sound speed at deeper water depths.

For the 3D simulations, the grid size in the marching direction was set to be 50 and 25 m for 5 and 20 Hz, respectively. This grid size setting was decided from a series of numerical convergence test. The same grid sizes were used for transverse and marching directions. The bathymetry and sound speed were updated through interpolation every 500 and 2500 m in the marching direction, respectively. The grid size in the water depth was considered 2 m. Two 3D sound propagation models were set up: (a) epicenter to CTBTO Network and (b) epicenter to SW06. The computation domains were a 1050 × 100 × 14 and 8655 × 500 × 14 km (x × y × z, marching × transverse × vertical direction) volume for CTBTO and SW06 model, respectively. An image field to satisfy the pressure-release boundary condition at the sea surface was applied. Artificial absorbing layers were placed at the ends of the vertical and transverse direction to mimic radiation boundary conditions. Simulations were run on an Intel(R) Xeon(R) CPU E5-2687 W v4 @ 3.00 GHz, 48 CPUs with 512 GB RAM memory.

In order to identify 3D propagation effects, N × 2D (Perkins and Baer, 1982) numerical simulations were also carried out. This type of simulation uses a 2D model with the same PE technique used in the 3D simulation but limiting the transverse coupling of sound energy and only considering its variation in the radial direction.

The bathymetric and sound speed models along the geodesic for epicenter to CTBTO Network and SW06 models are shown in Fig. 6. In this figure, y = 0 corresponds to the geodesic from the epicenter to hydrophones at Ascension Island CTBTO Network and SW06 Network.

The properties of the bottom considered in all the simulations are sound speed c_b = 1800 m/s, density ρ = 1500 kg/m³ and attenuation β_b = 0.1 dB/λ. Also, for all the simulations, water density is considered to be constant of ρ_w = 1000 kg/m³. Although sound propagation in the ocean bottom is out of the scope of this study, bottom properties can influence infrasound propagation in the water column. Due to this fact, we investigated the influence first, of three different bottom sound speeds (c_b = 1600, 1700, and 1800 m/s), and second, of an elastic bottom on the acoustic mode group velocity. Figure 7 compares the group velocity dispersion curves of the first mode for the four ocean bottom types analyzed considering a h = 4 km water depth ocean, water density ρ_w = 1000 kg/m³ and water sound speed c = 1500 m/s. A homogeneous fluid layer with a soft top and a homogeneous, higher velocity fluid bottom (Frisk, 1994) with bottom density ρ_b = 1500 kg/m³ is considered to study the influence of c_b. In the elastic bottom case (Ellis and Chapman, 1985), solid density is ρ_s = 2750 kg/m³, pressure wave velocity is c_p = 6300 m/s and shear-wave velocity is c_s = 3550 m/s. As we can see in Fig. 7, for frequencies higher than 5 Hz the group velocity is not too much influenced by the bottom properties. However, for frequencies lower than 5 Hz an accurate group velocity representation can only be obtained if the elastic properties of the bottom are considered.

For a 0-dB point source transmitting 5 Hz sound and placed on the bottom at 3300 m water depth, a comparison between 3D and N × 2D depth-integrated sound intensity results from the epicenter to Ascension Island CTBTO Network is shown in Fig. 8(a). Strong focusing, reflection, refraction and diffraction effects can be seen in the 3D model at different ranges. However, these effects are not reproduced by the N × 2D model. These 3D effects are clearly induced by the strong topographic changes along the propagation paths crossing the Mid-Atlantic Ridge. On the other hand, a shadow zone along the geodesic (y = 0) after 500 km range cannot be reproduced by the N × 2D model leading to an underestimation of depth-integrated TL up to 40 dB. Differences between the two models are also very evident in the shadow zone induced by Ascension Island due to the diffraction effect, where the North array of CTBTO Network is located (x = 1050 km, y = 0 km). Significant differences of TL can also be observed after 500 km range at different water depths when analyzing the sound intensity profile along the geodesic [see Fig. 8(b)].

3D effects are also found along the propagation paths from the epicenter to SW06 experiments for 5 Hz (see Fig. 9) due to the presence of the Mid-Atlantic Ridge. Differences between the 3D and N × 2D models are very evident at the acoustic shadow on the backside of Bermuda Islands [located at range 7600 km in Fig. 9(a)]. In some points of this shadow zone, the N × 2D model over-predicts the depth-integrated TL given by the 3D model by more than 60 dB. It should be noticed that also significant differences between both models along the shelf break and continental shelf shoreward are also observed.

Figure 10 presents the beamforming output for 5 and 20 Hz along the geodesic between epicenter and CTBTO network. Beamforming computations considered the pressure field at 846 m depth from the 3D PE model. A beam aperture of ±1500 and ±5000 m centered at y = 0 was considered for 20 and 5 Hz, respectively. Then, at each x position (over range) a beamforming analysis is performed to compute horizontal arrival angles. The apertures chosen for beamforming computation take into account that (i) aperture will affect beamforming results since larger aperture provide better resolution and (ii) possible effects of finite beamwidth. A cross-range pressure field resolution of half wavelength was used to avoid potential aliasing in beamforming computations. In Fig. 10 the zero degrees on the plot indicates the geodesic path, and as it can be seen, the arrival angles vary with the geodesic distance. Observing this figure, and overlapping the bathymetry with the beamforming outputs, one can find that the large angle deviation areas match with the Mid-Atlantic Ridge and Ascension Island both for 5 and 20 Hz.

Figure 11 gives the arrival horizontal angles at 76 m depth between the epicenter and SW06 experiments also for 5 and 20 Hz. Also in here, the large angle deviation areas match with the Mid-Atlantic Ridge. Moreover, large angle deviations can also be observed along the Shelf break and continental shelf close to the location of SW06 experiments. Hence, large angle deviation areas occur where 3D effects induced by bathymetric features are expected to be stronger. Therefore, beamforming results confirm the previous observations based on analysis of TL results.

Figure 12 shows the beam power patterns at CTBTO and SW06 locations obtained with the 3D model for 5 and 20 Hz. Level differences in Fig. 12 are because different array apertures were used for 5 and 20 Hz. The horizontal arrival angles (relative to the geodesic path) are higher at CTBTO station than SW06. At CTBTO the arrival angles are −2.99° and −13.75° for 5 and 20 Hz, respectively. The difference in the arrival angles for 5 and 20 Hz could be due to wavelength and vertical mode angle dependency. At SW06 the arrival angles decrease to −0.96° and −0.45° for 5 and 20 Hz, respectively. The broadband arrival angle measured at CTBTO (–5.1°) is in the interval of the modeled arrival angles for 5 and 20 Hz. However, at SW06 the measured broadband arrival angle (2.13°) is slightly higher than the arrival angles modeled.

Our study investigated the 3D effects induced by the Mid-Atlantic Ridge and Atlantic Islands on long-range propagation of low-frequency underwater sound from submarine earthquakes. We used a 3D geodesic Cartesian PE model to propagate infrasound T-phase waves (5 to 20 Hz) triggered by the Southern Mid-Atlantic Ridge earthquake of August 31, 2006. Model results lead us to conclude that sound propagation was extremely influenced by the complex 3D reflection effects induced by the ridge. In fact, results suggest that the ridge induced horizontal focusing, reflection, refraction, and diffraction in different parts of the geodesic path between the epicenter and CTBTO and SW06 hydrophone stations. The influence on sound propagation could be observed close to the sound source in the model, where the epicenter was located and the T-phase waves were triggered as well in the far field thousands of km away from the epicenter. Close to the source the ridge influence can be observed in the sound scattering from the source to the waveguide. In the far field, even though the bathymetry does not reach the channel axis at 1000 m, the rises of the sea floor still extend into the SOFAR channel, which is bounded by the conjugate depth at 4000 m, and influence the propagating sound over and around the ridge.

Our study also investigated the 3D effects of acoustic shadowing and diffraction induced by Atlantic islands on long-range submarine earthquake sound propagation. The acoustic shadow produced by Ascension Island during the earthquake, observed by the two CTBTO arrays situated on opposite sides of the island, was represented by the 3D model. On the contrary, the N × 2D PE model failed to provide an accurate prediction for the sound field around the island because it does not account for 3D diffraction around the island.

Although physical oceanographic features can cause horizontal refraction on long-range underwater sound propagation (Heaney and Campbell, 2016), numerical results suggested that observed discrepancies between expected and measured arrival angles could mainly be induced by the Mid-Atlantic Ridge, islands, shelf break, and continental shelf.

The broadband arrival angle measured at CTBTO is in the interval of the horizontal arrival angles given by the 3D model for 5 and 20 Hz. However, at SW06 the measured broadband arrival angle deviates about 3° from the angles given by the 3D model. This deviation could be due to uncertainties on SW06 horizontal array orientation, which could range between 1° to 2° from its reported orientation (Newhall et al., 2007). Moreover, as shown by the numerical results, the shelf break and shelf could have a significant impact on sound horizontal reflection induced by geological features. Nevertheless, the global 1/60° resolution bathymetry dataset used in this study cannot represent accurately such geological features in shallow waters.

In the numerical simulations, the sound source was placed at the ocean bottom and at the earthquake epicenter, which assumes that the T-phase wave was triggered by the submarine earthquake just above the hypocenter. However, T-phase waves are ultimately a seismic problem and the assumption that T-phase source is located at the earthquake epicenter may not be realistic (de Groot-Hedlin and Orcutt, 2001). Consequently, uncertainties on the exact location of the source of T-phase wave could lead to errors on the calculations of arrival angles at SW06 and CTBTO in the numerical simulations; however, the 3D effects observed in the models are promising, so it is believe that the horizontal reflection and refraction did contribute to the observed arrival deviations.

Due to the lack of knowledge on the exact sound source function, it becomes challenging to compare broadband model results with field data. In this sense, future works should focus on broadband underwater submarine sound propagation based on more well-known sources. This should help to improve the use of 3D sound propagation models for the localization of T-phase wave excitation locations.

In this work, our focus is at the bathymetric effects on 3D global scale sound propagation, especially in the presence of the mid-Atlantic ridge, also with oceanographic influence. Although this work does not consider bottom elasticity, this should be considered in future work. For that, 2D elastic parabolic equation solutions (Frank et al., 2013; Frank et al., 2015) can be adapted to improve the treatment of 3D T-phase generation in the global scale propagation problem.

The authors acknowledge David Bradley and Stephen Nichols for the access and support analyzing CTBTO data. The first author acknowledges the Postdoctoral Scholar Program at the Woods Hole Oceanographic Institution, with funding provided by the University of Haifa. The funding supports for the second author were from the Office of Naval Research, USA, under Grant No. N00014-17-1-2692.