CTBTO’s Data and Analysis Pertaining to the Search for the Missing Argentine Submarine ARA San Juan

The Argentine submarine ARA (Armada de la República Argentina) San Juan went missing on November 15th, 2017, with its last confirmed contact location around 600 kilometres East of Comodoro Rivadavia, Chubut, Argentina. In the days following the initiation of the multi-national search and rescue mission, the Comprehensive Nuclear-Test-Ban Treaty Organization (CTBTO) analysed data recorded by its International Monitoring System (IMS) sensor network for any signals that could support this effort. Two of the IMS hydrophone stations, HA10 in the central Atlantic Ocean and HA04 in the southern Indian Ocean, were found to have recorded an unusual hydroacoustic signal of unknown nature (a “hydroacoustic anomaly”), which originated from the vicinity of the last known location of the ARA San Juan. On December 1st, 2017, confirmation of the detection and localization capability of the IMS network was made possible by a calibration test involving a depth charge deployed by the Argentine Navy. Also this signal was detected by the two IMS hydrophone stations and its origin located to within 39 km of the announced deployment location. After extended search operations, on November 17th, 2018, the ARA San Juan was found on the seabed at approximately 900 metres depth less than 20 km from the location of the hydroacoustic anomaly provided by CTBTO. This paper presents the methodologies applied by CTBTO within months after the loss of the submarine, aimed at detecting, locating and characterizing the November 15th signal of unknown origin and the December 1st calibration explosion.


Introduction
On November 8th 2017, the Argentine submarine ARA San Juan departed from Ushuaia, the southernmost city in the world, on a cruise of approximately 11 days along a route of 2000 km in the southern Atlantic Ocean. The destination was the ARA San Juan's home base in Mar del Plata, Argentina (Fig. 1a). After 7 days of navigation, ARA San Juan contacted its home base at around 10:30 UTC on November 15th 2017 from a location approximately 600 km East of Comodoro Rivadavia, Chubut, Argentina (Fig. 1b). An extensive international search effort, led and coordinated by the Argentine Navy, was initiated and included more than a dozen surface vessels and airplanes equipped with various devices in an attempt to find the missing submarine. 1 The Comprehensive Nuclear-Test-Ban Treaty's (CTBT) International Monitoring System (IMS) consists of a world-wide network, comprising four different sensor technologies, seismic, hydroacoustic, infrasound and radionuclide, designed to detect nuclear explosions. Once completed, the network will consist of a total of 337 monitoring facilities distributed around the globe on the ground and in the oceans. At present about 90% of these monitoring stations are completed and certified as of August 30, 2019. The IMS delivers continuous data in near realtime to the International Data Centre (IDC) at the CTBTO Vienna Headquarters, Austria, for data processing and production of bulletins which are distributed to the CTBT Member States. The distribution of the world-wide sensor network is shown in Fig. 2.
The certification of hydrophone station HA04, Crozet, in June 2017 marked the completion of the hydroacoustic component of the IMS. There are in total eleven hydroacoustic stations, five of which are equipped with seismometers (T-stations, see legend in Fig. 2) to detect underwater signals converted to seismic waves near the shore. These seismometers may also contribute as auxiliary seismic stations. The other six stations comprise hydrophones deployed underwater to monitor the oceans for signs of nuclear explosions (see legend in Fig. 2). Each hydrophone station is typically composed of two triplets of hydrophones suspended in the Sound Fixing and Ranging (SOFAR) channel, one to the north and one to the south of an island except for HA01 Cape Leeuwin, Australia, which has only one triplet to the West. Placing hydrophones in the SOFAR channel ensures optimum performance of the hydro-acoustic stations in detecting signals from very long distances as a result of the low attenuation that a signal undergoes when propagating in this channel (Jensen et al. 2011).
In addition to Treaty verification, the signals recorded by the IMS network are also used for a Illustration of the approximate last known location of the ARA San Juan, indicated by the red dot, and the approximate location of the hydroacoustic anomaly detected by the CTBT IMS hydrophone stations, indicated by a yellow dot. The geodesic paths determined by estimating the back-azimuth at the hydroacoustic stations HA10 and HA04 towards the location of the hydroacoustic anomaly are indicated by the red lines in (a). b A magnified view of the area around the last known location and the estimated source location of the hydroacoustic anomaly broad range of civil and scientific applications. For example, the IMS sensors provide data for detection and real-time warning of earthquakes and tsunamis, radiation dispersal from nuclear accidents and volcanic eruptions. These data also make an important contribution to research on the Earth's core, climate change, meteorology, breaking-up of ice shelves and creation of icebergs, oceans and marine life, meteors and world-wide background radiation. 2 Shortly after the ARA San Juan was announced as missing, CTBTO undertook a comprehensive analysis of acoustic and seismic data recorded at IMS stations close to the announced search area. Additionally, also non-IMS stations were included in the analysis, in particular two local seismic stations. This analysis led to the identification and localization of a hydroacoustic anomaly found in the data recorded by the IMS hydrophone stations HA10 Ascension Island (United Kingdom), central Atlantic Ocean and HA04 Crozet Islands (France), southern Indian Ocean, on November 15th, 2017 (Fig. 1a). The event location was determined to be in the vicinity of the last known position of the submarine (Fig. 1b). The hydroacoustic anomaly was found to be an isolated event, relatively short in duration and broadband. The location of the signal, and details about its spectral and time-frequency characteristics, were provided to the Argentine authorities in a timely manner.
Two weeks after loss of contact with the ARA San Juan, the Argentine Navy conducted a controlled experiment by deploying a Mark 54 explosive depth charge on December 1st, 2017, to the north-east of the last known location of the submarine, to verify the detection and location capability of these two CTBT IMS hydrophone stations. The analysis presented in this paper deals mainly with signals recorded at HA10 Ascension Island and HA04 Crozet. The results obtained by the analysis reflect the CTBTO's effort in supporting the search activities (Nielsen et al. 2018(Nielsen et al. , 2019a.
A variety of techniques were employed, including an assessment of received power spectral levels, correlation and cepstral analysis, back-azimuth and arrival time estimation using a correlation algorithm, and analysis of direct and reflected signals from both, the November 15th event and the December 1st calibration explosion. The main purpose is to present and discuss similarities and differences between the signals originating from the November 15th event, whose origin was unknown, and those originating from the depth charge, and to provide a literature reference for the benefit of other authors who may wish to deal with the analysis of this event.
The paper is organized as follows: Sect. 2 describes the signals of the hydroacoustic anomaly of November 15th and the controlled calibration test of December 1st, recorded at HA10 and HA04, and their spectral, time-frequency and cepstral characteristics. Section 3 presents a modelling study of the propagation from the source to HA10, addressing the uncertainties in the propagation environment and the source depth, and their effect on the signal levels received at the hydrophone triplets at several thousands of km from the source. The estimation of the direction of arrival of the signals, including several reflections of the November 15th hydroacoustic anomaly from bathymetric features in the Atlantic Ocean received at HA10 and a pre-cursor associated with propagation through the seasonal Antarctic icesheet received at HA04, are discussed in Sect. 4. Section 5 presents the refinement of the event localization and the reduction of the uncertainty ellipse, commonly referred to also as coverage ellipse, using in addition to the hydrophone recordings also data from a non-IMS land based seismic station. Section 6 contains the conclusions.

Signal Detection and Association
On November 20th, 2017, the CTBTO initiated a lengthy scan for signals recorded on the two hydrophone stations HA10 and HA04 in support of the search for the missing Argentine submarine ARA San Juan. The submarine officially went missing on November 15th, 2017. The automatic processing algorithm employed at CTBTO to detect nuclear explosions on ground, in the ocean and in the atmosphere did not provide a hydroacoustic event detection automatically that could have triggered further inspection. There was no a priori indication that the hydrophone stations could detect a signal potentially associated with the loss of the ARA San Juan, and there were no indications regarding which signal characteristics could be associated with the missing submarine if any signal was emitted at all in relation to this event. Further, it could not be implied that the IDC processing algorithm could detect any signal associated to the missing submarine as the processing algorithm is designed and optimized to detect nuclear explosions whilst at the same time minimizing the number of ''false alarms''. Therefore, the algorithm necessarily tends to reject signals which cannot be associated with the primary signals of interest for CTBT monitoring.
A specialized data analysis based on the baseline IDC data analysis software was applied for a lengthy data scan in search for relevant features that might have been associable with the submarine search. Data was processed for the period November 15th-21st, 2017, as a batch process and stored to temporary result files for a manual post-processing and further inspection. The emergent results were the continuous time-bearing estimates from hydrophone time series received at stations HA10 and HA04. The geodesic path corresponding to the estimated bearing for each station was displayed to visualize possible intersections that could indicate a location of an event. These estimated locations and the associated signal arrival times were compared to the search area, and to the location and time of the last contact from the submarine to the Argentine Navy base on November 15th, as reported in the media at the time. The inspection of the data was constrained to signals that had an estimated bearing within that linked to the departure close to Ushuaia and the destination location at the Navy base in Mar del Plata, see Fig. 1. This was approximately the extent of the original search area reported by the press during those days.
The manual post-processing of the recorded signals identified an impulse-like arrival on the hydrophone station HA10 at around 14:58 UTC on November 15th, 2017. The primary signal was short in duration and had a sharp rise-time indicating a short and impulsive event happening most likely in the ocean. An arrival was also detected at H04S at around 15:20 UTC. The raw waveforms of these signals recorded on H10N, H10S and H04S are shown in Fig. 3, upper panels.
It should be noted that the H10S1 hydrophone is affected by strong electronic noise and is therefore not contributing data to IDC Operations. This electronic noise affects also the other two H10S hydrophones, as can be seen in the upper middle panel of Fig. 4. The extent to which this cross-talk affects the H10S2 and H10S3 hydrophones varies with time, with periods of high cross-talk and periods of low cross-talk. The H10S triplet will for that reason not be considered further in the data analysis presented here.
Initial back-tracing along the geodesic paths following the estimated back-azimuth from the two manually associated signals at the H10N and H04S triplet resulted in an intersection of the paths that was close to the last known location of ARA San Juan (Fig. 1a). This preliminary information from the manual data scan and associated analysis were provided to the CTBTO Lead Waveform Analysts, who were able to scrutinize and confirm the findings by using standard IDC interactive (i.e., human-assisted) analysis and review tools.
The signal received at H04S was affected by a stronger time-dispersion compared to the recordings on H10N and H10S. However, comparison of the power spectral densities shown in Fig. 5a, b showed a similar frequency content, a similar interference pattern and a similar roll-off in the early decaying part of the spectra of the two signals. These features gave further motivation to interpret the signals recorded at HA10 and HA04 as having originated from the same source. The power spectral densities of these signals are discussed further in Sect. 2.1 below.
The calibrated time-frequency spectrograms corresponding to these November 15th signals are shown in the upper panels of Fig. 4. The signals recorded at H10N and H10S are evidently short in duration and broadband, and the signal amplitudes are well above the background ambient noise level, whereas the signal recorded at H04S is clearly more affected by dispersion and hence somewhat longer in duration.
Multiple signals can be identified within a 24-min time window of the spectrograms computed for triplets H10N and H10S. These signals also exhibited similarities in terms of spectral content and timefrequency characteristics that made it possible to interpret them as horizontal reflections and refractions of the primary signal off the South American continent, islands and underwater seamounts. Also for all these signals, the spectrograms of the identified signals received at H04S confirm the increased timedispersion and a low-pass filtering compared to the HA10 recordings. These differences in the received signals are most likely caused by the propagation conditions along the paths from the event to the hydrophone triplets.
On December 1st the Argentine Navy conducted a controlled experiment to verify the detection and localization capability of the CTBT IMS hydrophone stations for impulsive events. A Mark 54 depth charge was deployed at location 45°40.0 0 S 59°24.90 0 W by an aircraft heading in the North-North-West direction. Using an equivalent TNT yield of 102 kg for the Mark 54 depth charge 3 and the observed cepstral delay of 0.45 s gives a detonation depth of approximately 30 m, using the equations published in (Willis 1963). The estimated deployment time was 20:04 UTC on December 1st at a location to the north-east of the last known position of ARA San Juan. The time series of the signal generated by the depth charge and recorded at the triplets H10N, H10S and H04S are shown in Fig. 3, lower panels. The December 1st time series are time-aligned with the signals recorded on November 15th for ease of comparison. The signal from the depth charge can be clearly seen on the three triplets and the arrival sequence between the three triplets is in striking agreement with the November 15th event.

Estimates of Received Power Spectral Levels
For comparison between the different recordings, the power spectral density (PSD) levels of the received signals from the November 15th and December 1st events were estimated for the hydrophones of stations H10N and H04S.
Welch's method (Welch 1967) was used to estimate the PSD for both signal and noise recorded on November 15th and December 1st. The start and termination time of the recorded signals were estimated by careful visual inspection of the raw time series and Welch's method was applied to the signal in this time window. A segment of similar duration as the signal, in which no arrivals could be identified, was used to estimate the PSD of the background noise. The time window for the estimation of the background noise was chosen during the time immediately preceding the signal of interest in order to represent the instantaneous noise level. The PSDs of the signal and noise were calibrated by applying the system response of the hydrophone stations. The December 1st 2017 controlled depth charge detonation recorded on the same sequence of hydrophones and triplets (lower panels) as for the November 15th 2017 event (upper panels). Note that no data are shown for the first hydrophone of triplet H10S, as this hydrophone is excluded from IDC Operations processing because of high levels of electronic noise. The extent to which this noise on H10S1 affects the detectability of signals on H10S2 and H10S3 through cross-talk varies over time, with periods of low cross talk and periods of strong cross-talk. The online version of the figure has higher resolution The calibrated PSDs averaged across the three hydrophones in each triplet for the two events recorded at the two hydrophone stations are shown in Fig. 5, where the black curve represents the signal from the event and the red curve the noise. The background noise level is slightly higher on December 1st compared to November 15th at both stations. However, the signal-to-noise ratio (SNR) is 10-20 decibels (dB) for the December 1st and up to 30 dB for the November 15th event.
It is noted that the higher received levels of the November 15th signal compared to the December 1st signal cannot be used to imply that the November 15th acoustic source was stronger, or more energetic, than the December 1st depth charge, as uncertainties about the source depth and the underwater environmental parameters along the propagation path can justify the observed difference in received levels. This aspect is further highlighted in Sect. 3 on propagation modelling.
On the other hand, the November 15th event shows clearly dominating low-frequency components and a faster roll-off at higher frequencies compared to the December 1st event. It is also noteworthy that the received levels of the November 15th signal at H04S are higher at the lower frequencies than those at H10N, although the estimated propagation distance to H04S is more than 1500 km longer. A possible explanation for this observation may be also in this case differences in the underwater propagation conditions along the two propagation paths from the same source to the two stations, which can affect the structure of the propagating signal. The increased noise level at H04S in the band from 18 to 25 Hz and in particular the singlefrequency signal at 25 Hz observed both on November 15th and December 1st are not related to the hydroacoustic anomaly or the controlled explosion. This band of relatively high noise level is almost continuously present and is generally associated with distant shipping traffic, airgun surveys or vocalizing marine mammals.
The signal PSDs from both events, the November 15th hydroacoustic anomaly and the December 1st controlled explosion, exhibit an interference pattern (or scalloping) as a function of frequency. Such an effect can generally be caused by a delayed replica of the primary signal that appears within the duration of time series considered. The delayed replica of the primary signal can be generated by two similar events separated by a short time, such as for example a reflection off the sea surface or seabed in the vicinity of the event or by an oscillating bubble pulse generated by an underwater explosion (Chapman 1985(Chapman , 1988Yamada et al. 2016). Based on the analysis of the recorded signal, it is justified to state that the December 1st depth charge most likely generated at least one bubble pulse which was observed as a delayed replica in the signal, or echo, at both hydrophone stations. However, the scalloping in the November 15th PSDs does not necessarily constitute evidence of an underwater explosion as the nature of the acoustic source is unknown a priori. The difference in frequency between the first two spectral peaks is approximately 2.8 Hz for the November 15th event and 2.0 Hz for the December 1st event.

Auto-correlation and Cepstral Processing
Auto-correlation processing is typically employed to improve the SNR and to reveal delayed replicas of a primary signal where these may be masked by the duration of the primary signal and time dispersion. The auto-correlation of a signal is simply the correlation of the signal by itself. The auto-correlation function can be calculated in the frequency domain by (1) computing the Fourier transform of the time series to obtain the complex amplitude spectrum, (2) multiplying it with the complex conjugate of the same Fourier transform and then (3) by computing the inverse Fourier transform.
The cepstrum analysis (Oppenheim and Schafer 2010) is somewhat similar to calculating the autocorrelation function and it is often used for analysing hydrophone data recorded from underwater explosions (Chapman 1985(Chapman , 1988Yamada et al. 2016). The difference between the auto-correlation and the cepstrum is that the logarithm function is applied to the amplitude squared of the spectrum in the cepstrum analysis before inverse transforming back to the time domain (or ''quefrency domain'' in cepstral analysis terminology). One of the advantages of the logarithmic function from the signal processing point of view is that it operates like an amplitude shading on the original amplitude spectrum, and subband processing can be obtained by adding logarithmic functions to the amplitude squared spectra instead of multiplying spectra. If the signal amplitude spectrum is composed of a combination of a primary signal and a delayed and scaled copy of this primary signal (echo), then the logarithmic function applied to the amplitude spectrum of the signal containing the primary signal and the echo separates. The separated part containing the echo has a specific periodicity in frequency corresponding to the inverse of the delay time between the primary signal and the echo (Oppenheim and Schafer 2010). This is what typically occurs as a result of the cepstral analysis of signals from underwater explosions where the generated gas bubble expands and contracts possibly multiple times. These multiple expansions and contractions of the gas bubble represent delayed and scaled echoes of the primary signal, the properties of which can be investigated by the cepstrum analysis.
The time resolution of the correlation or cepstrum is in general inversely proportional to the bandwidth of the signal. The time resolution provides a lower limit of which delayed replicas can be detected and separated from the primary signal. The correlation and cepstrum time-compress the time series to a bandlimited impulse and the delayed replicas can be detected although they are masked by the primary signal with longer duration than the delay time, or if the replicas are buried in ambient noise.
For completeness, both the auto-correlation function and the cepstrum were calculated for the November 15th and December 1st signals at the two hydrophone stations. The auto-correlation function (upper panel) and cepstrum (lower panel) for the November 15th signal are shown in Fig. 6a, c, respectively, and for the December 1st signal in Fig. 6b, d, respectively. The auto-correlation and cepstrum for the November 15th signal indicate an apparent delayed replica at around 0.33-0.34 s after the primary signal at both H10N and H04S. This similarity in the auto-correlation and cepstrum is a further indication that the November 15th signal was emitted from the same source, in addition to the spectral features described above in relation to the PSD's.
Such echoes detected by correlation and cepstral analyses might be indicative of a bubble pulse due to the first oscillation of a gas bubble (such as is the case with underwater explosions), a strong reflection shortly after the primary signal (e.g., from nearby bathymetric features) or two very similar underwater events happening within a short time. The delayed echoes from both the November 15th auto-correlation and cepstrum analysis reveal an apparent stronger time-dispersed replica at H04S than at H10N, which may be caused by a natural underwater propagation effect along the propagation path from the event to H04S such as low-pass filtering effects, strong ocean surface or seabed interaction.
It also appears that the delayed replica arrives slightly earlier at H04S than at H10N. Underwater explosions are isotropic sources, and it has been demonstrated that the delay of the bubble pulse from an underwater explosion is independent of the propagation distance from the event even for propagation distances differing by several thousands of Vol. 178, (2021) CTBTO's Data and Analysis Pertaining to the Search for the Missing Argentine kilometres (Chapman 1985(Chapman , 1988Yamada et al. 2016). This is in contrast with the appearance of the delayed echoes from the November 15th signal analysis recorded at HA10 and HA04, and hence these echoes do not corroborate the hypothesis that the November 15th event was an underwater explosion and indicate that the source might have been anisotropic (Blanc et al. 2018).
Correlation and cepstral analysis of controlled inwater explosive tests typically show repeated bubble pulses as the gas volume expands and recompresses possibly many times as the bubble rises towards the sea surface (Chapman 1985(Chapman , 1988. The detected signals from the December 1st controlled depth charge detonation do not appear to include successive auto-correlation or cepstral delays of bubble oscillations beyond the one presented in Fig. 6b, d. The detected delayed bubble pulse appears 0.45-0.46 s after the primary signal, and the delay time is persistent and independent of the difference in propagation distances from event to the two hydrophone triplets (Chapman 1985(Chapman , 1988Yamada et al. 2016). An example of a past study showing the second cepstral peak from an underwater explosion recorded at a hydroacoustic station and the absence of this second peak at a more distant hydroacoustic station is given in (Yamada et al. 2016). The explanation for this was simply in the lower SNR of the second bubble pulse compared to the first one.
The delayed echo is very similar between both the correlation and the cepstral analysis applied to the signals recorded on the two hydrophone triplets. It appears that the delayed echo for the November 15th event in both the correlation and cepstrum (Fig. 6a, c) is more dispersed in time for the propagation path from source location to HA04 than to HA10. This Auto-correlation function averaged over the three hydrophones at triplet H10N (black curves) and H04S (red curves) for the November 15th event (a) and the December 1st event (b). Cepstrum averaged over the three hydrophones at triplet H10N (black curves) and H04S (red curves) for the November 15th event (c) and the December 1st event (d) difference in dispersion of the delayed echo for the two propagation paths is not observed for the December 1st explosion source (Fig. 6b, d). This might indicate a difference in the nature of the underwater acoustic source of the December 1st test compared to the source of the November 15th event or changes in the propagation conditions from the estimated event locations to HA04 during November 15th to December 1st by for instance reducing the low-pass filtering effect. There exists a well-known empirical relation between the delay time of the first bubble pulse from an explosion and the explosion yield in kilograms equivalent to TNT and detonation depth (Chapman 1985(Chapman , 1988. This delay time is proportional to the yield and inversely proportional to the detonation depth. The relation is often used in analysing recorded signals from controlled explosion experiments (Chapman 1985(Chapman , 1988Willis 1963). As anticipated in Sect. 1, a fairly good agreement between the estimated delay time and observations has been obtained by applying this empirical relation for the December 1st depth charge when using the reported detonation depth and yield for the depth charge. Applying the same method to the November 15th event is of limited usefulness, as the precise nature of this acoustic source was unknown.

Propagation Modelling
In an idealized situation, the source strength for both the November 15th and the December 1st event could be determined by adding the transmission loss (TL) at the hydrophones to the received levels. However, the TL depends both on source-hydrophone geometry and the underwater environmental properties (ocean sound speed, bathymetry and seabed properties). The source depth of the November 15th event is also unknown, while the detonation depth of the depth charge deployed on December 1st is known, albeit with uncertainties on the accuracy. The underwater environmental properties can be extracted from high-resolution oceanographic databases, although also these come with the caveat of some level of uncertainty. Nevertheless, propagation modelling can be useful for understanding the different aspects which might have affected the signals presented in this study.
The two-dimensional (2D) underwater acoustic propagation model RAM (Collins 1993), based on solving the Parabolic Equation (PE) in range-dependent underwater environments, was employed to demonstrate the dependence of TL on source depth for the November 15th event. Only the 2D propagation path along the geodesic from the estimated source location to H10N was considered. The underwater environmental input for RAM was extracted from the Copernicus ocean temperature and salinity database as a function of depth (E.U. Copernicus Marine Service Information), and the water depth from the GEBCO (GEBCO Digital Atlas 2003) bathymetry database. The ocean sound speed was calculated based on the temperature and salinity profiles as a function of depth provided by the Copernicus database using Eq. 1.1 in (Jensen et al. 2011). The seabed was assumed to have sand-like properties.
The TL as a function of propagation range and depth at the acoustic frequency of 5 Hz and a source depth of 700 m is shown in Fig. 7, upper left panel. The maximum propagation range towards H10N is shown out to 6100 km while the H10N hydrophone triplet is at approximately 6050 km from the event location. In this case, the acoustic field interacts significantly with the bathymetry at short range from the acoustic source because of the downward refracting ocean sound speed profile and relatively shallow water depth. At longer ranges the field is guided by the SOFAR channel out to a distance of around 2500 km, where the acoustic field is scattered and attenuated by the interaction with the Rio Grande Rise ocean ridge. The acoustic field is then again guided back into the SOFAR channel after the Rio Grande Rise and propagates practically undisturbed until it reaches H10N at a range of 6050 km from the source location. Acoustic scattering also appears at the underwater seamount on which the H10N triplet is moored. The TL vs. range at a receiver depth of 800 m representative of the hydrophone depths at H10N is shown in Fig. 7 lower left panel, where the impact of the Rio Grande Rise and the underwater seamount at the H10N mooring location can be observed in the form of increased transmission loss. Vol. 178, (2021) CTBTO's Data and Analysis Pertaining to the Search for the Missing Argentine 2567 The TL was also computed in the band from 1 to 100 Hz in 4 Hz increments and for a source depth between 50 and 800 m in 50 m increments. The frequency averaged TL at a range of 6050 km from the source location and displayed as a function of source depth is shown in Fig. 7, right panel. The TL decreases markedly around 15 dB for a source depth from 50 to 300 m, while the decrease is more gradual below 300 m depth. This shows clearly the difficulty of estimating source strength for an unknown source depth by compensating observed received levels with a computed TL, as even one single problem parameter such as source depth can introduce large changes in the transmission loss.

Estimate of Direction of Arrival
The Multi-Channel Correlator DTK-PMCC (Cansi 1995;Cansi and Klinger 1997) was applied to the recorded time series to estimate the horizontal direction of arrival at the hydrophone triplets and to explain the multiple arrivals observed in the spectrograms (Fig. 4). Only the results from the November 15th event recorded at H10N are considered here although similar analyses were performed also for the depth charge on December 1st and triplet H04S. The DTK-PMCC performs analysis on bandpass filtered segments of the time series and cross correlates these segments pairwise for the three hydrophones of the triplet. The delay time for each of the three cross correlations and the triplet geometry are used to solve for two slowness components in the horizontal plane of the hydrophones. The back-azimuth is estimated as the arc tangent to the two slowness components by assuming a horizontally propagating plane wave.
The upper panel of Fig. 8 shows the estimated back-azimuth at H10N for a signal time window of around 15 min for multiple frequency sub-bands and time segments. Each colour corresponds to a backazimuth and is associated to an individual arrival numbered from 1 to 11. The colour is very homogeneous for each arrival in time and frequency, which indicates that the back-azimuth is unique and well determined. A total of 10 later arrivals besides the primary signal were identified by the DTK-PMCC within the 15-min time window shown here. The back-trajectory of each of the individual arrivals identified at H10N was manually traced until it intersected with an apparent obstacle (continent, island or underwater seamount). These apparent obstacles were then assumed to have scattered the trajectory towards the estimated location of the November 15th event (Fig. 8, lower panel). The travel times for all trajectories relative to the first arrival, obtained by this method and assuming a constant sound-speed value of 1485 m/s, are in close agreement with the observed individual arrival times in the upper panel.
The IDC standard processing algorithm and review tools only use one back-azimuth for each triplet to locate an event. Therefore, at least two distant triplets are necessary to perform good estimates of an event location. However, one triplet would suffice in locating an event if multiple horizontal arrival paths, such as those shown in Fig. 8 lower panel, were included in the location algorithm (Dall'Osto 2019a, b;Vergoz et al. 2019).
For the November 15th event recorded at triplet H04S, an apparent precursor was also detected, arriving 1234 s earlier than the main arrival. Spectral analysis indicated that the precursor originated from the November 15th event as both precursor and main arrival have some similarity in the interference patterns of the spectra and similar spectral roll-off. However, the precursor has a lower SNR and underwent more low-pass filtering compared to the main arrival, which most likely prevents the detection of a delayed replica in the auto-correlation or cepstrum. The correlation and cepstrum of the precursor reveal relatively high amplitude fluctuations within a time window close to zero lag time (or zero frequency) where a delayed echo is expected to be observed as for the main arrival. These fluctuations cannot be associated to and may mask the delayed echo.
The back-azimuth estimates for the precursor and main arrival were determined by applying the DTK-PMCC algorithm (Cansi 1995;Cansi and Klinger 1997) similar to the analysis shown for the H10N arrivals in Fig. 8. The DTK-PMCC result for the H04S arrivals is shown in Fig. 9 upper panel for the back-azimuth and apparent wave speed across the triplet. The apparent wave speed is simply calculated as the travel distance between the hydrophones of known geometry projected onto the direction of arrival of the assumed plane wave divided by the corresponding travel time. This results in an apparent wave speed representing the speed of sound in the ocean if the plane wave propagates horizontally. A plane wave propagating across the triplet at a vertical inclination angle will therefore appear to have an apparent wave speed that is higher than the speed of sound in the ocean. The precursor (Path 2) and main arrival (Path 1) have an apparent wave speed of around 1480 m/s, which indicates that the arrivals are sounds propagating horizontally in the water column. The estimates of the corresponding back-azimuths are homogeneous in time and frequency indicating a well-determined back-azimuth for both arrivals. The back-azimuth for Path 1 is 223.6°relative to North and represents the geodesic path from H04S used in the CTBTO location of the November 15th event. The precursor arrives along a geodesic path with a more southerly back-azimuth (Path 2) of 207.8°relative to North.
The seasonal Antarctic ice extent for November 15th 2017 was retrieved from the US National/Naval Ice Center and is shown by the green curve in Fig. 9, lower panel. The geodesic propagation path for the precursor (yellow path in Fig. 9 lower panel) and main (red path in Fig. 9 lower panel) arrival are likely to interact with the ice sheet, altering the signal propagation conditions compared to an ice-free sea surface. It is thus hypothesised that the geodesic path for the main arrival (Path 1) was partly covered by the ice sheet as indicated in Fig. 9 lower panel, and that the acoustic signal may have had significant interaction with the ice-sheet because of the typical Antarctic upward refracting ocean sound speed profile. The sound travelling along the upward refracted path may have penetrated the ice sheet, undergoing a high absorption in particular at the higher frequencies. This may explain the strong low-pass filtering of the main arrival recorded at H04S compared to the arrival at H10N for the November 15th event (Nielsen et al. 2019b).
Assuming that the precursor and main arrival were associated to the same event, the precursor must partly have propagated in a medium with higher sound speed than the water in which the main arrival has propagated, as the precursor arrived before the main arrival. Therefore, the hypothesis formulated here is that the November 15th signal coupled into the ice sheet and propagated through it. One of the possible in-ice paths that could explain the precursor would assume propagation through a relatively thin ice sheet in form of a symmetric Lamb-wave like mode (Press and Ewing 1951). The terminology ''thin'' refers to the small thickness of the ice sheet compared to the wavelength of the sound at these frequencies. The signal is assumed to then have coupled back into the ocean (Path 2 in Fig. 8 lower panel) at a location provided by the intersection of the geodesic path from H04S along the estimated Path 2 back-azimuth and the ice-sheet boundary. It is noted, that the wave traveling in the ice sheet would have had a speed of more than twice (Press and Ewing 1951;Jensen et al. 2011;Williams et al. 2018) the speed in water and could thus have arrived before the main arrival.
A simple calculation of travel time from the estimated November 15th event location to H04S, assuming a water to ice-sheet coupling and back to water geodesic signal propagation paths, was performed to test the above hypothesis. The water sound speed was assumed as constant 1480 m/s and two sound speed values of 3087 m/s and 3500 m/s in ice were applied. The first ice sound-speed is based on the speed associated with the Lamb wave hypothesis formulated above, while the second sound-speed corresponds to typical propagation speeds of a compressional elastic wave in ice. The intersection point of the geodesic propagation path from the estimated November 15th event location was then swept along the ice-sheet boundary. A geodesic propagation path was assumed from that point through the ice sheet that connected with the fixed ice-water intersection point provided by the geodesic path from H04S in the direction of the estimated back-azimuth of Path 2.
All swept paths incompatible with the time difference between the main arrival and precursor are shown as grey paths in Fig. 10. The yellow and blue curves represent paths that result in a time difference of 1234 s between the main arrival and precursor assuming the propagation speed in the ice sheet for the Lamb wave case and for the pure compressional wave case, respectively, and a speed of sound of 1480 m/s for those portions of the propagation path not covered by ice. These two paths are consistent with the hypothesis of a possible through the icesheet propagation path resulting in a precursor arriving earlier than the main arrival (Nielsen et al. 2018(Nielsen et al. , 2019a. The characteristics of the acoustic propagation of this signal under and in the ice sheet are the subject of an ongoing research effort.

Refinement of the Event Localization
The November 15th event localization obtained by using only the two hydrophone stations HA10 and HA04 is associated with large uncertainties, which can be quantified by considering the size of the 90% coverage ellipse. Other non-IMS stations located close to the estimated event location were identified to contribute to the event location and to minimize the size of the coverage ellipse.
The three-component seismic station TRQA (Tornquist on the IRIS/USGS network) and the EFI (Mount Kent on the IRIS/IDA network) detected a seismic signal compatible with an arrival time from a signal originating at the location of the hydroacoustic anomaly (13:53:12.3 UTC at TRQA). The two seismometers are situated less than 1000 km from this location. A polarization analysis similar to the method described in (Jurkevics 1988) was applied to the TRQA isolated arrival filtered in the band from 2 to 5 Hz. The SNR was determined to be too low on EFI to contribute to the event localization. The derived back azimuth, inclination and rectilinearity from a moving time window along the time series is shown in Fig. 11 lower panel. The stable azimuth (around 170°with respect to North) and inclination (around 40°with respect to vertical) and the growing rectlinearity (up to 0.66) after the arrival of the identified signal indicate the arrival of a seismic P-wave. The back azimuth of the arriving seismic P-wave is in a direction which is consistent with the back-azimuth to the location of the hydroacoustic anomaly.
The localization of the November 15th event included data recorded at the H10N and H04S triplets and the TRQA seismic station by using the standard IDC automatic localization algorithm with a following review by the waveform analyst. The automatic localization algorithm is based on an iterative, nonlinear weighted least-square minimization of origin time and azimuth residuals between observations and modelling results (IDC Processing of Seismic, Hydroacoustic and Infrasonic Data 2002). The origin time modelling results are obtained from estimated spatial (range and azimuth for each hydroacoustic station) and monthly dependent travel time tables by applying an underwater acoustic signal propagation model similar to what was used to generate Fig. 7 in Sect. 3 for all possible hydroacoustic event locations. The azimuth modelling results for all possible hydroacoustic event locations are obtained from the cross-correlation of the hydrophone signals for each HA station hydrophone triplet detecting the signal as described in Sect. 4. The event location (origin) is determined by the known station locations and the values of origin time and azimuths that minimize the residuals. The weighting is introduced by assuming a pre-defined variance of the combined measurement and modelling origin time and azimuth errors. The variance for the unknown origin time and azimuth is associated with measurement and modelling uncertainties, although the localization algorithm does not directly provide a confidence interval of the event localization based on these uncertainties. The weighting by errors introduces stability in the inversion for event origin time and azimuth rather than improving the accuracy of the event location itself (IDC Processing of Seismic, Hydroacoustic and Infrasonic Data 2002).
In the IDC standard localization algorithm, the prior measurement and modelling errors for origin time and azimuth are assumed to be known. The a-posterior residuals after the event location has been performed are not considered in estimating the event location uncertainty. Under these two assumptions, the location error is described as a coverage ellipse with a Chi squared uncertainty distribution. A Chi squared distribution assumes that the distribution of the individual unknowns is normal, i.e., follows Gaussian statistics. If the inversion for the event Figure 12 Best estimate of event localizations obtained by the IDC standard localization algorithm: December 1st test explosion (red dot) and ground truth (orange dot), November 15th event using only the two hydrophone triplets (yellow dot) and including the TRQA seismic station (green dot barely distinguishable from the yellow dot). The associated coverage ellipses are shown in red, yellow and green, respectively. The black dot indicates the reported location where the ARA San Juan was found during the night of November 16-17th, 2018 (https://www.lanacion. com.ar/politica/el-clima-adverso-le-dio-dramatismo-a-un-hallazgo-con-el-ultimo-aliento-nid2192921) location has converged, then the uncertainty of the location is provided by the coverage ellipse centred on the determined event location. The coverage ellipse in the IDC algorithm is calculated for the 90% spatial confidence interval with 90% probability that the true location is inside the coverage ellipse. The locations of both the November 15th and the December 1st events are shown in Fig. 12. The estimated December 1st depth charge detonation location (red dot) is only around 39 km from the announced drop location of the December 1st depth charge (orange dot). However, only the hydrophone triplets H10N and H04S were used in the localization, which results in a large coverage ellipse (red ellipse) of 37,212 km 2 . No seismic station was included in the depth charge analysis as the signal was not detected on any of the considered seismic stations. However, the estimated localization is impressively close to the ground truth considering the underwater acoustic signal propagation distances of more than 6000 km from the event's source to the hydrophone triplets.
Similar conclusions can be drawn for the November 15th event location using only hydrophone triplets. The coverage ellipse (yellow ellipse) is relatively large (33,659 km 2 ). However, it is reduced significantly to 748 km 2 by including the seismic station TRQA (green ellipse). The best estimate of the event location is nevertheless very similar in both cases, when using hydrophone stations only and including also the seismic station in the analysis.
A summary of the localization results is provided in Table 1. The location where the ARA San Juan was found resting on the seabed at a depth of around 900 m (black dot in Fig. 10) is at a distance of less than 20 km from the location of the November 15th, 2017, hydroacoustic anomaly provided by CTBTO. 4

Conclusions
When the submarine ARA San Juan went missing on November 15th, 2017, CTBTO provided to the Argentine authorities timely information based on acoustic signals acquired by two hydroacoustic stations, namely HA10 at Ascension Island and HA04 at the Crozet Islands. The distance between the location where the ARA San Juan made its last reported communication with the base and these hydroacoustic stations was circa 6000 km and 7700 km, respectively. The accuracy of the detection and source Table 1 Summary of event location, coverage ellipse characteristics and estimated origin time for the November 15th and December 1st 2017 hydroacoustic events https://www.lanacion.com.ar/politica/el-clima-adverso-le-dio-dramatismo-a-un-hallazgo-con-el-ultimo-aliento-nid2192921 localization was corroborated via a calibration experiment conducted by the Argentine Navy after the incident. The localization area was further refined using complementary data from the Argentine seismic station TRQA. One year later, the ARA San Juan was found on the seabed at approximately 900 m depth, less than 20 km from source location of the November 15th hydroacoustic anomaly detected by CTBTO. The analysis of the acoustic signals associated with this case is still on-going with the intention to examine all key features associated with the original signals and the various arrivals at the monitoring stations.

Acknowledgements
The CTBTO would like to acknowledge all the countries that host IMS stations. The timely availability of IMS data is crucial for the mission of the CTBT but also, in cases such as this, for humanitarian purposes. The memory of the tragic events which led to the analysis of data presented here is bound to remain with the authors of this paper and their colleagues at CTBTO, all of whom wish to express their heartfelt sympathies and condolences to those affected by the loss of the ARA San Juan.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4. 0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.