Source triangulation utilizing three-dimensional arrivals: Application to the search for the ARA San Juan submarine

On Nov. 15, 2017, an intense signal was detected on an underwater acoustic hydrophone network, coincident with the disappearance of the Argentine navy submarine ARA San Juan. The purpose-built acoustic network is part of the larger global International Monitoring System (IMS) operated and maintained by the Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO) to verify compliance with the treaty. This continuous, real-time acoustic monitoring of the worlds' oceans provided crucial evidence to find the missing submarine. Analysis of the signal detected at two stations, HA10 at Ascension Island in the Mid-Atlantic and the other HA04 at the Crozet Islands (Ile de la Possession) in the Southern Indian Ocean, provided a location where a search effort ensued (see Fig. 1). A significant refinement, established by timing of a land-based seismic detection, reduced this search area to a circle of 20 km radius.¹ Unfortunately, even though the initial search effort encompassed the location where the ARA San Juan was eventually found in 2018, it was unsuccessful at locating the submarine.

In September of the following year, the search effort was restarted by the seabed exploration company Ocean Infinity. After searching the high priority area established by the CTBTO findings,¹ the author was contacted to provide a confirmation of this location. One striking feature in these recordings are three-dimensional (3-D) propagation features observed at HA10,¹ manifest as delayed arrivals with different azimuth than the primary, main arrival. The highest amplitude of these delayed arrivals corresponds to acoustic energy bathymetric refracted,² redirected seaward though multiple bottom interactions as the wavefronts propagated up the continental slope.³ While this observation ruled out possible locations on the shelf (as sources there do not produce refracted arrivals), the timing of these signals, as will be shown, is of great utility for source triangulation. Exactly a year and a day after ARA San Juan's disappearance, the submarine was discovered following this extensive 2+ month effort to acoustically map the ocean floor.⁴ 

In an attempt to refine the location and confirm the CTBTO analysis, an air-dropped depth-charge was deployed on Dec. 1, 2017 by the Argentine military in the same region [46° S, 60° W]. As will be shown, deployment of the depth charge as a known source was a good idea. Not only did it confirm the signal characteristics of an impulsive signal that propagated 8000 km range, it allowed for verification of modeled propagation speed. Still, when only the main arrival is considered at two stations, there is ambiguity in when or where the event occurred. To resolve ambiguity in the San Juan event origin, CTBTO analysis included detection on a land-based seismic station (TRQA) roughly 120 km north of the origin.¹ This third detection decouples event time and location, however, this same improvement could be made with an additional hydroacoustic detection. Detection on a third hydroacoustic station is not necessary, as distinct propagation paths of the observed 3-D arrivals add “virtual” stations which can satisfy the three detection triangulation requirement.

Virtual stations are built by forming the representative 3-D energy pathlines, connecting the origin of the sound to the hydrophones. Three-dimensional propagation paths are modeled with both a horizontal modal-ray description^5,6 and energy pathlines based on an adiabatic-mode parabolic equation (PE) solution.⁷ Modal-rays and energy pathlines are effectively equivalent in deep-water propagation. The PE method is favored for modeling propagation over a realistic continental slope, as it provides an unambiguous formulation of acoustic intensity and requires no distinction between a reflection or refraction process.⁸ While the analysis is limited here to the study of the first propagating mode of 10-Hz sound, this reduced model captures the relevant aspects global propagation “involving significant bottom interaction and azimuthal coupling.”⁷ 

Triangulation is accomplished here in two parts: (1) a directional assessment of sound detected on the multiple hydrophones at each station, and (2) triangulation through the finding the intersection of contours of equal time (or isochrones) representing the possible prior location of signal wavefronts. Geosynchronous timing of this hydroacoustic network allows for great precision in triangulation. Error in localization is ultimately due to uncertainty in precisely when an event occurred, along with under-constrained estimates of the oceanography between the source and hydrophones.

With event time and location defined, the measured propagation speed relates directly to the ocean temperature structure (i.e., ocean tomography). Here, propagation speeds are estimated from the range and depth dependent sound speed (see Fig. 2) derived from a data product developed by the consortium for “Estimating the Circulation and Climate of the Ocean” (ECCO2).⁹ Complimented with direct sampling and space-based oceanographic data, the network and its 20-yr hydroacoustic record is a trove of data for study of trends in ocean temperature.^10,11 While the focus of this study is on additional information garnered from 3-D arrivals needed for triangulation, it lays the framework to infer the origin and time of sources of opportunity with potential for ocean thermometry using historical and future CTBT data to improve our basic understanding of the under-sampled, deep ocean temperature structure.

The balance of this paper is organized as follows: Section II details the hydroacoustic stations that recorded the San Juan and depth-charge event, and characterizes the time and frequency dependent features of the signals. Section III covers the procedure for source triangulation, and addresses uncertainties based on the propagation assumptions. Section IV implements a model to reconstruct one of the observed 3-D acoustic propagation paths arising from refraction by the continental slope, and uses this to triangulate the origin of the two acoustic events. Section V concludes the paper with a short summary and discussion on how 3-D propagation features can reduce uncertainties in localization.

There are six hydroacoustic (HA) stations integrated into the CTBTO IMS to support verification of member states' compliance with the CTBT. An HA station consists of either one or two hydrophone triplets, an array of three hydrophones on an equilateral triangle with ∼2 km baseline separation. This layout provides an unambiguous measurement of sound direction. When combined with geo-synchronized timing of arrivals, multiple stations are used to triangulate the origin of sound sources.

The hydrophones are suspended as close as possible to the axis of the sound channel, near the minimum of sound speed where decreasing ocean temperature is balanced by increasing density. This minimum forms the axis of a sound channel, where sound can propagate without reflection from the seafloor or sea surface. This allows sound to propagate great distances without the scattering and loss mechanisms from incurred by boundary reflections. Although the sound channel structure is not present near the polar region, sound energy still does not interact with the ocean floor where most loss occurs. Here hydrophone triplets are suspended closer to mid-water, at depths ∼500 m, to provide optimal sensitivity.

The two hydroacoustic stations that sensed the sound events near [46° S, 60° W] were HA10 at Ascension Island [8° S, 14° 30′ W] in the mid-Atlantic Ocean, and HA04 at the Crozet Islands [46° 30′ S, 51° 45′ E] in the Southern Indian Ocean. Both HA10 and HA04 consist of a pair of hydrophone triplets, separated by ∼200 km, with one triplet located to the north and the other to the south of the island base station providing power and communications (see Fig. 3). During the two-week period Nov. 15 to Dec. 1, 2017, five hydrophones were in operation at HA10 (the complete triplet to the north—H10N, and only two of the three hydrophones at the south—H10S). Only the southern triplet of the HA04 station (H04S) was operational during the Nov. 15th event, although both were operational for the Dec. 1st depth-charge detonation.

Perhaps what is most revealing to the location of these two events are the additional arrivals appearing after the main arrivals at HA10 (Ascension Island). Like the main arrival, these additional 3-D arrivals all exhibit the repetitive spectral feature associated with the source signature. In fact, ten of such late arrivals are documented at HA10,¹ although the focus here will be on the most significant of these that occurs ∼90 s after the direct arrival. This arrival corresponds to acoustic energy refracted from the continental slope, and will be used in triangulation of the sound events. While the HA04 recordings do not have distinct, delayed 3-D arrivals, they exhibit another fascinating propagation feature: an advanced, faster than seawater arrival (see Fig. 4 at 14:59) presumed to be related to propagation through Antarctic ice.¹ 

On Nov. 15, 2017, a distinct signal well above the ambient noise level was received at HA10 (peak arrival for H10S at 14:57:57 Coordinated Universal Time (UTC) and for H10N at 14:59:08 UTC) and at HA04 (peak arrival at 15:19:26 UTC) (see Fig. 4). The southern station of HA10 was the first station to detect the signal, and measured the greatest peak pressure, a level of 145 dB re 1 μPa. This impulse signal is well-dispersed, with 95% of the energy arriving within 16 s at the south station and within 21 s for the north station. A distinct peak in this signal occurs 4.79 s after the onset for the south station, and after 5.32 s for the north station. This peak is likely associated with the mode-1 arrival, while the earlier arrivals are higher order modes which in deep-water propagation have faster group velocities. The difference in dispersion at the north station is due to the additional 100 km of propagation distance, as supported by the ∼70 s delay relative to the south station.

The signal arriving at HA04 is much more dispersed, over 35 s, and does not have a distinct peak. A defining feature of the signal is the mode-1 dispersion, readily observed in a time-frequency analysis. Due to the slowing group velocity with frequency, the peak of the signal starts at 3 Hz and follows a linear upsweep in frequency, at a rate of 0.4 Hz per second.¹³ This feature corresponds to mode-1 of the surface duct formed by the frigid polar waters. In this linear gradient sound-speed profile, the group velocity of mode-1 at lower frequencies is faster (extends deeper into the ocean where sound speed is greater) than at higher frequencies. Remember this dispersion feature, as it is an important aspect of matching the propagation time of mode-1 in the signal.¹⁴ 

The frequency content appears to be limited to $≳3 Hz$⁠, with a peak arrival occurring at ∼10 Hz (see Fig. 5). Since the dispersion of the signals at the two stations look drastically different, the periodic nature of the frequency content was paramount to associating these arrivals to the same event. The periodicity to the spectral levels, Δf, is presumed to be caused by second impulse, or doublet, delayed from the first by 1/Δf. This feature is consistent at the three locations (H10N1, H10S3, and H04S2), and upon auto-correlation or cepstral analysis¹ (i.e., delay spectrum) is revealed to be 335 ms ± 15 ms. Note that the frequency peaks and valleys at each station are not identical, but both spectra show Δf = 2.9 Hz.

On Dec. 1, 2017, two weeks after the incident, an air dropped Mk-54 depth-charge [45 40 S, 59 24 W] was detonated within the general search area. Again, a distinct signal well above the ambient noise level (see Fig. 6) was received at HA10 (first arrival for H10S at 21:10:45 UTC and for H10N at 21:11:57 UTC) and at HA04 (first arrival at 21:32:42 UTC). This detection confirmed the existence of propagation channels to both stations, and the richness of 3-D arrivals at HA10. Dispersion characteristics of this signal are nearly identical to that of the Nov. 15th acoustic event, albeit the signal is less intense. Again, energy concentration was observed in the early part of the HA10 arrivals, while the HA04 arrival is more evenly dispersed. Unlike the San Juan signal, the HA04 energy is more evenly distributed in modes, which makes identifying a clear mode-1 “feature” difficult. This is a potential source of error caused by “feature picking” inaccuracy.

Similarly, the Dec 1 depth charge also exhibits a distinct repetitive spectral feature from which to associate azimuthally distributed 3-D arrivals. As the source signature of explosions is well documented, the repetitive feature is believed to be caused by the doublet impulse formed by the detonation and subsequent bubble collapse. This delay can be predicted based on the charge parameters as¹⁵ 

where W is charge weight (NEW, or net equivalent weight of TNT), z is the charge depth, both in units MKS, and the constant 10.1 represents atmospheric pressure in meters of seawater. The parameters of this Mk-54 charge (an ordnance now obsolete) are known: charge weight, 250 lb Torpex multiplied by a 1.3 conversion factor to NEW,¹⁶ W = 147 kg NEW; planned detonation depth, z = 38.1 m; and notional detonation time 20:04:00 UTC. Equation (1) predicts a bubble period of 440 ms, consistent with the observed value 455 ms ± 2 ms (Δf = 2.2 Hz) considering minor errors in conversion factor, charge weight and detonation depth.

Note that although a doublet signature was observed in the Nov. 15th signal associated with the ARA San Juan tragedy, there is no absolute evidence showing its natural cause.

The IMS hydroacoustic network is designed to localize the origin of intense sound anywhere in the oceans of the world. The triad of sensors at each HA station can form a back-azimuth to the source bearing, as input for triangulation. Furthermore, the geosynchronous measurement of arrivals on multiple stations allows for extrapolation of contours of equal travel time, called isochrones, which intersect at the origin of the source. Combined, these inputs produce an area of probable source location (an error ellipse) which can be reduced with detections on additional stations.¹⁷ 

The bearing of the acoustic arrivals at each station can be formed either via complex spectra, or through delay corresponding to the peak in cross-correlation between two hydrophone stations. The high signal-to-noise ratio (SNR) of the two events proves no distinct advantage of either method; both produce comparable estimates of arrival angle. The arrival angle, as computed from the delay Δt, and the relative location of the two hydrophones follows the relation

where $ϕ$ is the incident bearing, $ϕ0$ is the broadside alignment of the two hydrophones, r₀ is the horizontal separation of the hydrophones, and the ± represents the left–right ambiguity of the estimate. When reversed, these angles form back-azimuths that can be traced along their geodesic to the source location (see Fig. 7).

Ideally, the intersection of these two paths represents the source location such that the timing of the signal (T) satisfies,

where $dk¯r/dω$ is the reciprocal average group velocity over the propagation path, which has length L. Source location can be determined without back-azimuths by finding the isochron intersection, which for three stations is distinct and establishes an event time (the received arrival time less the propagation time, T). Computation of isochrones with Eq. (3) requires a global map of group velocity, which again depends on the latitude and longitude (and temporally dependent) sound speed profile. As input, a global model of ocean temperature, salinity and depth are translated to depth dependent sound-speed profiles via empirical equations.¹⁸ The sound-speed profiles along the paths from the two stations show a strong sound channel in the mid-latitudes, sharply transitioning at the circumpolar front to surface duct propagation in the southern latitudes (see Fig. 2). After accounting for the curvature of the Earth through an Earth flattening transformation, the sound speed profiles are then translated to a group velocity to produce the global group velocity map for a particular mode and frequency.

Each pair of hydrophones provides a back-azimuth: three from H10N, one from H10S, and three from H04S, while each hydrophone (eight total for Nov. 15, 2017) provides an arrival time to calculate isochrones. Separation of stations in both in range and azimuth is desirable, as this makes intersection points of back-azimuths and isochrones are distinct. To form isochrones, travel time is integrated along back-azimuths. A representative group velocity of the signal at 10 Hz, taken as the averaged group velocity of mode-1 along the propagation path (see Fig. 2), is 1482 m/s at HA10 and 1455 m/s at HA04. A uniform geoacoustic parameterization representative of the continental slope,^19,20 is implemented, characterized by a sound speed 1.14 times, and density two times greater than seawater.

Based on a notional event time of 20:04:00 (arrivals outlined in Table I) the isochron intersection is 50 km from the actual source location. The intersection occurs along a line in the east–west direction, for an uncertain arrival time. This orientation is due to the relative location of the stations. If we allow for the notional detonation time to be inaccurate, a correction of +31 s brings the isochron intersection to within 3 km. To correct for this offset, the modeled propagation speed can be adjusted (e.g., through modifying ocean temperature) although it is worth noting that this correction may not be physically justifiable.

In order to triangulate with the additional 3-D arrivals, a model of the horizontally refracted propagation path (its energy pathline) is required. One approximate representation of the energy pathline is through horizontally refracted modal-eigenrays. As refraction is caused by gradients in phase speed, e.g., steepness of contour levels in Fig. 1, at 10-Hz refraction is resulting primarily from changes in ocean depth. As such, modal-rays are a great tool to determine if sound can reach a particular place in the ocean, but estimating the intensity of refracted arrivals (for example, by examining ray density or through forming a ray-tube) proves difficult. Specifically, as the refraction process involves multiple aspects of acoustic propagation, e.g., forward scattering, diffraction, etc., interpretation of scattered rays from a realistic bathymetry is not straightforward. With these considerations, and to convince the reader that this horizontally refracted propagation path is intense enough to detect on the CTBT hydroacoustic stations, a wave-based numerical code is implemented to identify distinct refracted wavefronts.

Energy pathlines, built from a full-wave simulation of the acoustic intensity vector field, provide an accurate ray-like depiction of the propagation path a wavefront follows.²¹ Here, the vector field is constructed from gridded pressure field simulations via finite-difference approximation to the pressure gradient.²² Results of a “broad-band” simulation are summed (i.e., Fourier synthesis representative of an impulse) to add time-dependence and, given sufficient bandwidth, provides time separation between direct and refracted arrivals to form energy pathlines for each discrete arrival. Here, triangulation of the San Juan and depth-charge events is based on pathlines constructed from hybrid single mode-parabolic equation (PE) simulation of the horizontal propagation of mode-1 at 10 Hz.

First, mode-1 properties: horizontal wavenumber, k_r, and group velocity, at 10 Hz are computed over the modeling region (e.g., using a mode-code Kraken¹²). For simplicity, a half-space bottom model is used for these computations (density ratio at the interface of 2, and sound speed ratio of 1.14, and attenuation of 0.25 dB/λ), noting that group velocity of the refracted path will be most sensitive to this characterization. Next, the horizontal wavenumber (k_r) is used as input to the gridded environment of a PE propagator, specifically phase velocity (2π/Re{k_r}) as the propagation speed and the imaginary part of k_r as attenuation.²² While this analysis implements the ocean state estimate from ECCO2,⁹ commensurate results can be obtained from tabulated seasonal ocean atlases as the low-frequency propagation is only effected by slowly varying temperature structure.

Figure 8 shows a snapshot in time of the broad-band intensity field, constructed by a summation of the complex narrow-band fields from 9.7 to 10.3 Hz computed in 1/100 Hz increments. Based on this sampling, the concentric rings are wavefronts representing multiple instances of the same source (at the Dec. 1 depth-charge location) repeating every 100 s. Roughly 350 km north of the source, the wavefronts begin to exhibit bathymetric refraction, turning down-slope due to the relative sharp gradient in phase speed as the water shallows approaching the continental shelf (contours in Fig. 8 are at 25 m/s increments of phase speed, between 1475 and 2000 m/s). Still, some sound propagates up onto the shelf, where the propagation is subject to substantial bottom loss corresponding to a large imaginary part of k_r. While this on-shelf sound decays relatively quickly compared to deep-water propagation, the reverse situation—detection of on-shelf sources in deep water—is possible from such intense signals.²³ 

The vectors shown in Fig. 8 identify the intensity direction of the refracted and direct wavefronts, with those corresponding to the refracted 200 s advanced from the others. Energy pathlines connect these vectors to the source, constructed by starting at their end (arrows) stepping the field backwards in time applying a phase factor of e^iωt with t defined in fixed time steps. Depending on which propagation feature is of interest, a recipe to trace a energy pathline is: take an initial point on a wavefront, compute the horizontal intensity, and then back-propagate (step backwards in time) along that direction at the local group velocity. This process is repeated incrementally until the pathline traces the wavefront back to the source, smoothly tracing a path through the bathymetric refraction process. The propagation time is integrated along the pathline following Eq. (3).

Application of this numerical technique to an entire ocean basin is computationally expensive, and some reduction of the field is desirable to efficiently handle 3-D propagation effects. Pathlines are straight and well characterized as rays in deeper water, and application of a full-wave solution is only necessary in the vicinity of the continental shelf, where modes approach cut-off depth. Conceptually, these pathlines form a connection to a piece-wise ray description along which isochrones of refracted paths can be formed. Thus, to maintain efficiency when propagating through blue-water over thousands of kilometers, rays are launched to continue pathlines after the refraction process.¹³ 

Pathlines in Fig. 8 are incorporated into the triangulation by continuing propagation to HA10 along the geodesics (labeled 1 and 2 in Fig. 8), and in the opposite direction extrapolated beyond the origin. Isochrones are formed for each arrival perpendicular to the respective pathline. Figure 9 shows the triangulation isochron intersection based on the extracted signal parameters in Table I. The triangulation falls within 3 km of the documented depth charge location (star symbol) predicting an event time of 20:04:31, and similarly for the San Juan (square symbol) predicting an event time 13:51:16. Uncertainty in triangulation is ultimately due error in the modeled group velocity, or precision in feature picking the mode arrival time within the recorded signals. Uncertainty in triangulation for ±1 m/s error in group velocity (or 3 s error in feature picking) corresponds to ±5 km accuracy in triangulation. Factors contributing to this uncertainty include the oceanographic model input, as well as characterization of geoacoustic properties the continental slope sediments. Sediment characterization significantly effects the group velocity where sound interacts with the seafloor, i.e., ocean depth less than ∼1000 m for mode-1, 10-Hz sound.

This paper detailed the underwater acoustic recordings of a massive acoustic event, made at distances 6000 and 8000 km, relating to the loss of the submarine ARA San Juan. These detections provided crucial evidence to define a priority search area,¹ where eventually the wreckage was discovered.⁴ This study was initially motivated to support the search effort and a qualitative assessment was provided based on the observed 3-D arrivals: limit the search area to depths greater than the 250 m contour following the continental slope. The San Juan was discovered on the continental slope [45°56.4′ W, 59°48.6′ S] along the 926 m depth contour. Development of a 3-D triangulation technique continued after the San Juan was discovered, and that technique is applied here to recordings of that event and a Mk54 depth-charge detonated two weeks after the disappearance in the same general location.

Precision in triangulation improves by including more stations, and at a minimum three arrivals are needed to decouple event time from location. In light of the limited number of hydroacoustic stations in the CTBT network (six total, worldwide), additional “virtual stations” are formed by modeling the propagation paths of 3-D arrivals recorded at a hydroacoustic station. While all source locations may not be advantageous to produce detectable 3-D arrivals, especially for low level signals, the intensity of a nuclear detonation will most certainty produce detectable 3-D arrivals, regardless of location.²⁴ With such 3-D analysis, the CTBT network station density is sufficient to accomplish its mission.

For the triangulation of the events detailed here, one additional station is needed to supplement the main arrival received at two hydroacoustic stations (such a seismic stations¹). Here a virtual station was formed from the first and most intense delayed 3-D arrival, corresponding to sound that propagated up the continental slope and refracted back toward the mid-Atlantic station (HA10). Analysis here was simplified to consider only the mode-1 arrival at 10 Hz, chosen based on its distinct arrival time and at low enough frequency where oceanic refraction from horizontal temperature gradients is insignificant.⁵ Now with the San Juan found, quantitative assessment of the triangulation based on 3-D arrivals is possible. For triangulation, the parameters controlling uncertainty are feature picking, i.e., identifying the modal arrival time within the recorded signals, and the estimates of modal group velocity for each path. Triangulation of both signals was based on modern maps of ocean temperature produced estimates accurate to within 5 km, corresponding to +/–3 s accuracy in feature picking (or 1 m/s error in the group velocity).

From a scientific perspective, recordings of impulsive signals with a precisely known origin (time and location) provide great insight into oceanographic processes—they can be used for ocean tomography (or thermometry, as temperature is the primary driver of sound speed variability). Through matching observed and modeled signals, an iterative process can refine a constrained ocean state estimate by minimizing the error. While in practice known sources are used for tomography, the 3-D triangulation technique presented here provides a framework to do historical thermometry with the CTBT network. With this goal in mind, and complimented with direct sampling and space-based oceanographic measurements, historical and future analysis will necessarily improve our understanding of the deep ocean temperature structure.

This work was supported by ONR code 32, and made possible by the virtual Data Exploration Center (vDEC) maintained by CTBTO. Thanks to colleagues at CTBTO, who deserve commendation for their excellent work to quickly respond, analyze, and document these events. Their report in the Proceedings of the European Geophysical Union¹ is an excellent summary of their effort. And thanks to Andy Sherrell, for trusting the precision of the acoustic evidence which led to eventual discovery of the ARA San Juan. Deepest sympathies go out to the families and loved ones of the submarinistas lost to this tragedy.