# Automatic Analysis of Hydroacoustic Signals Related to the Activity of the Fani Maoré Submarine Volcano

Pierre-yves Raumer, Sara Bazin, Romain Safran, Dorian Cazau, Jean-yves Royer

► To cite this version:

Pierre-yves Raumer, Sara Bazin, Romain Safran, Dorian Cazau, Jean-yves Royer. Automatic Analysis of Hydroacoustic Signals Related to the Activity of the Fani Maoré Submarine Volcano. Geochemistry, Geophysics, Geosystems, 2025, 26 (12), 10.1029/2025GC012572. hal-05443182

HAL Id: hal-05443182

https://hal.univ-brest.fr/hal-05443182v1

Submitted on 6 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Geochemistry, Geophysics, Geosystems®

## RESEARCH ARTICLE

10.1029/2025GC012572

**Key Points:**

* New method enables to automatically catalog hydroacoustic P-phases and volcanic H-waves recorded by four moored hydrophones

* Demonstration of successful seismicity cataloging around the Fani Maoré volcano in the Mozambique channel, from continuous hydroacoustic monitoring

* Analysis of the hydroacoustic catalog reveals evolving b-values in two distinct seismic clusters

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

P.-Y. Raumer,
pierreyves.raumer@gmail.com

**Citation:**

Raumer, P.-Y., Bazin, S., Safran, R., Cazau, D., & Royer, J.-Y. (2025). Automatic analysis of hydroacoustic signals related to the activity of the Fani Maoré submarine volcano. *Geochemistry, Geophysics, Geosystems*, 26, e2025GC012572. https://doi.org/10.1029/2025GC012572

Received 10 JUL 2025

Accepted 15 NOV 2025

<figure><img src="image_2.png" /></figure>

## Automatic Analysis of Hydroacoustic Signals Related to the Activity of the Fani Maoré Submarine Volcano

Pierre-Yves Raumer¹,² ID, Sara Bazin¹ ID, Romain Safran¹ ID, Dorian Cazau², and Jean-Yves Royer¹ ID

¹Geo-Ocean, University of Brest, CNRS, Ifremer, IUEM, Plouzané, France, ²Laboratoire des Sciences et Techniques de l'Information, de la Communication et de la Connaissance (Lab-STICC), École Nationale Supérieure de Techniques Avancées (ENSTA), Brest, France

**Abstract** Due to the properties of sound propagation underwater, including the existence of the SOFAR channel, sounds in the ocean can propagate over large distances with little attenuation. This makes passive acoustics a relevant method for monitoring natural events such as ice calving, earthquakes, and underwater volcanic eruptions. However, the lack of automated techniques (e.g., such as beamforming for small-aperture array) for large (>50 km) hydrophone arrays makes the data analysis very time-consuming; indeed, analysts must recognize signals such as earthquake T-phases on several hydrophones, manually pick and associate their arrival times to build event catalogs. This lengthy iterative process by trial and error often leads to underutilization of the recorded data and exposes the resulting catalogs to the analyst's subjectivity. Here we propose a pipeline that fully automates the creation of event catalogs from continuous hydroacoustic data in two steps: (a) detection of events at each recording station and (b) association of detections at multiple stations to trilaterate the source. In order to evaluate this approach, we tested it on real data recorded by the Mayotte Hydrophone Network (MAHY), which monitors the activity of the new underwater volcano off Mayotte Island in the Mozambique Channel. Its seismicity and volcanic activity were analyzed with the catalogs obtained.

**Plain Language Summary** Sounds in the ocean can propagate over very large distances. As a result, sensors for underwater sounds, called hydrophones, are very relevant for monitoring natural sources of “noisy” events and, in particular, for detecting small earthquakes over large distances. However, manual analysis of such data is cumbersome and time-consuming. We propose a new suite of tools to partially automate this task. It automatically detects events in different hydrophone records, then gathers detections corresponding to the same event, and finally obtains a location and time origin of the source that generated the detected sounds. We tested this processing pipeline on real data to demonstrate its relevance for building earthquake catalogs from hydroacoustic data.

## 1. Introduction

The vast field of passive hydroacoustics enables one to monitor many sound sources with a single sensor, the hydrophone. The underwater soundscape results from a combination of biological, anthropogenic and geophysical sound sources. The latter category of sounds includes, in addition to weather-related noise (e.g., Gros-Martial et al., 2025), earthquake generated sounds (e.g., Royer et al., 2015), seawater-magma interactions (e.g., Bazin et al., 2023, 2025), landslides (e.g., Ingale et al., 2025), or hydrothermal activity (e.g., Smith & Barclay, 2024). T-waves, hydroacoustic waves related to earthquakes, mainly result from the conversion of seismic waves into acoustic waves in the area of its epicenter; when captured by the Sound Fixing and Ranging (SOFAR) channel, T-waves can travel over very long distances (>1,000 km) with little attenuation. Consequently, hydrophones are sensitive to magnitudes as small as 3.0–3.5 on a near-ocean scale range (Dziak et al., 2004; Ingale et al., 2021), and provide insights about sources as diverse as major subduction earthquakes (e.g., De Groot-Hedlin, 2005), strike-slip events along transform faults (e.g., Parnell-Turner et al., 2022), normal faulting events (e.g., Olive & Escartín, 2016) or magmato-tectonic swarms along spreading ridges (e.g., Ingale et al., 2021; Tsang-Hin-Sun et al., 2016). Direct arrivals of seismic phases can also convert into acoustic waves at the foot of moored hydrophones, if the event is close enough or has a high enough magnitude (Dziak et al., 2004); then, a hydrophone can quasi be used as a seismometer to perform some basic tomography. In such case, the direct P- and S-wave arrivals converted in acoustic waves can be termed as P- and S-phases, to disambiguate them from

© 2025 The Author(s). Geochemistry, Geophysics, Geosystems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
T-phases. The energy content of these three types of phases is typically distributed in the 2–60 Hz frequency-band.

More recently, some hydrophones (moored or fixed to Ocean Bottom Seismometers (OBS)) recorded hydroacoustic waves corresponding to energetic impulsive signals in the 2–120 Hz frequency-band. Based on their source co-location with active lava-flows, they seem to be related to lava-seawater interactions (Bazin et al., 2023; Le Saout et al., 2020; Tan et al., 2016); however, the mechanism which generates them is still debated (Saurel et al., 2022). Steam burst could involve phase transitions, but at a 3,000 m depth, pressure prevents the formation of gas bubbles. Supercritical water or volcanic gases such as CO₂ may however exist and lead to volumetric changes (Schipper & White, 2010). Cavitation-like phenomena could also cause a local pressure variability. Implosions, explosions and physical state changes may then produce acoustic signals. Lava quenching, which does not require compressible matter (Murch et al., 2022), may be another possible acoustic source (Bazin et al., 2025). Recently, Wang et al. (2025) extensively investigated the matter and finally opted for the implosion of bubbles generated by trapped water or hydrated sediments. These signals will be called H-waves as they are directly produced in the water column.

The analysis of T-wave signals has long proven to be relevant for studying low-magnitude events at oceanic basin scale, but most of the times relies on manual processing (e.g., Fox et al., 2001; Ingale et al., 2021; Tsang-Hin-Sun et al., 2016). As in standard seismology, analysts must go through a “picking” stage of wave arrival-times, then an association stage by trial-and-error to combine them for the trilateration stage to determine the location and origin-time of the events, to ultimately build a seismic catalog ready for analysis. Such approach is cumbersome, user-dependent, and often hinders the full exploitation of available data. Although it is possible to automate the detection and/or classification of geophysical signals in a hydroacoustic data set (e.g., Raumer et al., 2023, 2024; Sukhovich et al., 2014), the steps of association and trilateration remain necessary to obtain a complete catalog. Methods used to analyze acoustic data from small-aperture arrays (e.g., beamforming, cross-correlation) require a high coherence between signals from different stations, and thus are not applicable for large-aperture arrays (>50 km). For such arrays, this paper proposes a fully automated processing pipeline implementing detection, association and trilateration steps, and tests it on hydroacoustic data from an array composed of four hydrophones forming a square of 70 km × 70 km, called the MAHY network (Figure 1). It was set up to monitor the seismic and volcanic activity of the Fani Maoré underwater volcano.

Fani Maoré is a newborn underwater volcano 50 km off Mayotte Island, the easternmost island of the Comoros Archipelago, West of Madagascar. Following an intense seismic activity in May 2018, a scientific expedition conducted in 2019 led to the discovery of a new edifice approximately 820 m high (Feuillet et al., 2021). It also led to the completion of the seismologic network on Mayotte Island and to the establishment of the REVOSIMA observatory (REseau VOlcanologique et SIsmologique de MAyotte), which, offshore, includes OBS's and four moored hydrophones, the latter forming the MAHY network (Saurel et al., 2022). Since then, land-based seismometers and OBS's detected many Volcano-Tectonic (VT), Long-Period (LP) or Very-Long-Period (VLP) events, which spread over two main clusters: an island proximal cluster speculated to correspond to a deep magmatic reservoir and a distal cluster corresponding to the narrow plumbing system that fed the new volcanic edifice (Lavayssière & Retailleau, 2023; Retailleau et al., 2023). Jacques et al. (2024) proposed a piston model where a piston-like structure delineated by the proximal earthquake swarm sustained the driving pressure of magma withdrawal, leading to the eruption. West of the main edifice, near the proximal cluster, a complex region named the Horseshoe area showed high gas (CO₂, CH₄) concentration (Heumann et al., 2025). The MAHY hydrophones, installed in October 2020, dated the end of the lava outpouring to 4 December 2020 when all H-waves stopped, marking the end of the final eruption phase (Lavayssière et al., 2024). All the manually detected and located H-waves originated from the Tiktak area, a circular lava field, 5–10 km North-West of the main volcano edifice. Their Source Levels (SL) reached approximately 190 dB, comparable to a 300 cubic-inch airgun-like source.

The seismo-volcanic activity after 4 December 2020, however, still needs a systematic and comprehensive investigation to verify the absence of H-waves in the longer term. Moreover, given that the local seismicity is well constrained by the OBS network (Saurel et al., 2022), the full potential of the available hydroacoustic data for seismic monitoring has not been thoroughly assessed. As an application to our new automated approach, we propose to address these two issues. This paper will first present the new hydroacoustic automatic analysis
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Location of the MAHY network of autonomous hydrophones (AuH). Yellow stars are the positions of the four AuH's. Bathymetry compiled from MAYOBS surveys (Rinnert et al., 2019). The red triangle corresponds to the new submarine volcano named “Fani Maoré.” The white crosstar is the epicenter of the earthquake whose signal is shown in Figure 2, as detected by MAHY01.</figcaption></figure>

methodology, and second test its performance on the MAHY data, and evaluate its relevance to explore the
volcano-tectonic underwater soundscape in this region.

**1.1. Data and Signals of Interest**

The MAHY network consists of four AuH, sampling at 240 Hz, and moored at ~1,300 m depth in the SOFAR channel in order to enhance the detection capability of hydroacoustic waves. They were deployed in October 2020 and were recovered after 6 months, and since then, redeployed at the same sites on a quasi-yearly basis (Figure 1). In the recovery process, a comparison with a GPS clock enables us to estimate the sensor clock drift and to correct the data time-base assuming this drift is linear with time, so that all hydrophones are synchronized. In 1 year, the clock drift is typically less than 2 s. Some battery outages however prevented this correction upon a few recoveries. Stations are named MAHY*ij, with *i* being the deployment number and *j* the sensor location among the four moorings. MAHY*ij designates all the deployments at location *j* regardless of the recording period, and MAHY*i*, all the stations of deployment number *i*. Data from MAHY0 to MAHY4 were collected from October 2020 to October 2024, during the following oceanographic campaigns: MAYOBS15 (Rinnert et al., 2020), MAYOBS18 (Rinnert et al., 2021a), MAYOBS21 (Rinnert et al., 2021b), MAYOBS23 (Jorry et al., 2022), MAYOBS25 (Thinon & Lebas, 2023) and MAYOBS30 (Komorowski & Paquet, 2024). It must be noted that during MAHY2, from late 2021 to mid-2022, MAHY24 mooring was not redeployed due to lack of time, resulting in a network with only three stations.

The inspection of spectrograms reveals many signals that can be grouped in two categories: seismic-related waves lasting ~1 min and H-waves lasting only ~5 s (Figure 2). The MAHY AuH's being close (<100 km) to the epicenters, one could expect the recorded waves to include P-, S-, T-phases in addition to secondary phases as found in similar situations (Dziak et al., 2004; Slack et al., 1999). Typical phases observed in the waveforms are
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Spectrograms (top) and waveforms (bottom) of an earthquake (left) of magnitude *M*<sub>L</sub> = 3.9 detected by MAHY01 and of a H-wave (right) as detected by MAHY04. Wave arrivals: P phase, Pw Surface-reflected P-phase, Pww Bottom reflected Pw, and S phase. The bottom plot is an enlargement of the earthquake signal. The black dashes mark the time origin of the earthquake as given by the reference catalog, blue dashes show the P-phase arrivals, green dashes that of S-phase arrival, and the golden rectangle represents the time-interval where T-phases would be expected. On the bottom plot, P-phases are picked at 10.88, 12.51, and 14.63 s, and the S-phase at 16.44 s. The spectrogram color scale is log-normalized.</figcaption></figure>

shown in Figure 2. There are three P-wave arrivals, marked by blue dashed lines, with two secondary P-phases interpreted as surface and ocean-bottom reflections of the first arrival. These P-phase arrival times are consistent with a conversion from seismic to hydroacoustic just below the AuH, with a subsequent vertical propagation in the water column. The event epicenter is shown in Figure 1 and its hypocenter was 40 km deep. The travel-time calculation is expanded in the Text S1 in Supporting Information S1. Our interpretation discarded the possibility of these arrivals being T-phases. In the given geometry (Figure 3), the yellow rectangle in Figure 2 outlines the time interval where T-phase arrivals would be expected; the first P-arrivals are clearly ahead of time due to the fast travel of P-waves in the crust (~5,000–8,000 ms⁻¹) compared to the travel time of T-waves in the water column (~1,490 ms⁻¹). Similar predictions were made for all the examples that we inspected. In the second signal category, one can observe a clear impulsive signal or H-wave followed by a secondary arrival due to a sea-surface reflection (Figure 2). Thus, in the absence of clearly distinguishable T-phases at this range, our analysis focused on the detection of P-phases and H-waves. Figure 3 shows the celerity models in the solid and water media and the ray-paths of the two signals of interest. In the water-column, the sound-speed was derived from the WOA23 global database at this location (Reagan et al., 2024); in the crust, we used the 1-D regional model determined by Lavayssière et al. (2022).

## 2. The Developed Pipeline

Our overall aim is to automate the construction of an event catalog from a set of raw hydroacoustic data. The most common approach consists of at least three steps as follows: the analyst first browses simultaneously records from the available hydrophones, displayed as spectrograms and/or time-signals, then picks signals that likely originate from the same event, and finally attempts to localize its source in a grid with a least-squares algorithm. The event
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Predicted ray-paths of P- and H-waves, respectively originating from an earthquake at 40 km in the crust, and from the seafloor. Left: Hydroacoustic (Reagan et al., 2024) and seismic celerity (Lavayssière et al., 2022) profiles, both linearly interpolated. Right: Propagation of (purple) a H-wave as modeled by Bellhop (Porter, 2011) and (blue) a P-phase as modeled by TauP (Crotwell et al., 1999) and a linear acoustic path. The bathymetry profile (black line), the P-wave source (bottom star) position and the AuH (black dot for MAHY04) position correspond to the experimental setting and data shown in Figure 2. The two less attenuated rays of H-waves are displayed in dark purple while highly attenuated rays are visible in the background. Note the larger vertical exaggeration (x10) for the water column relative to the crustal section.</figcaption></figure>

association is tested by trial-and-error based on the trilateration results: whether the trilateration algorithm is able or not to converge to a solution, or the misfit is acceptable or not, or the location is plausible or not. Once a plausible location/association is found, the process is generally reiterated by refining the arrival-time picking in zoomed-in time windows. These three steps can be named detection, association and filtering. For the detection part, we first adapted an already developed automatic scheme, mainly designed to detect T-waves (Raumer et al., 2024), then we automated the association part, including a filtering part. As in any automatic method, an analyst must still closely examine the final results, but will be left to focus only on limited problematic recording periods or questionable events.

## 2.1. Detection

### 2.1.1. Definition

Let detection be the process of finding, given some acoustic data of a single station, the precise time of arrival of a particular hydroacoustic event. Given a station *s*, the expected output takes the shape of a set of detection times $T_s$, where each $T_s$ is a date.

### 2.1.2. Method and Evaluation

This work focused on two types of geophysical signals: direct crustal ray-path from earthquake (P-phases) and volcanic H-waves. To detect them in a noisy soundscape and take advantage of a time-frequency representation, an automatic method based on spectrograms is used.

The process started with the published *TiSSNet* model (Raumer et al., 2024), a Fully Convolutional Network (FCN) developed for such detection purpose and trained using a large-scale AuH network in the Indian Ocean, OHASISBIO (Royer, 2009). The ambient sound around Mayotte Island appears different from the one present in the training data of *TiSSNet*, encompassing more low-frequency anthropic signals such as ship traffic and seismic
<table><caption>**Table 1**
Summary of the Training Results of TiSSNet</caption><thead><tr><th rowspan="2">Model</th><th colspan="3">Training samples</th><th colspan="2">Evaluation</th></tr><tr><th>Positive</th><th>Negative</th><th>F1</th><th>Samples</th><th>Precision</th></tr></thead><tbody><tr><td>P-TiSSNet (1st)</td><td>909</td><td>4,595</td><td><strong>77.0%</strong></td><td>1,000</td><td><strong>87.4%</strong></td></tr><tr><td>H-TiSSNet (1st)</td><td>5,102</td><td>4,734</td><td><strong>88.7%</strong></td><td>1,000</td><td><strong>45.0%</strong></td></tr><tr><td>H-TiSSNet (2nd)</td><td>6,701</td><td>13,634</td><td><strong>77.1%</strong></td><td>1,000</td><td><strong>79.2%</strong></td></tr></tbody></table>

Note. P-TiSSNet corresponds to the model used to detect P-phases and is obtained by fine-tuning the published TiSSNet model using MAHY data. H-TiSSNet corresponds to the model used to detect H-waves and is obtained in a similar way, requiring two iterations of fine tuning starting from the annotations of Lavayssière et al. (2024). The F1 score, intended to account for both accuracy and recall, is defined by the formula $\frac{2TP}{2TP+FP+FN}$, where TP, FP, FN stand for True Positive, False Positive and False Negative. The bold values highlight the metric score of each training. Only the precision was evaluated on the evaluation part, this being defined as $\frac{TP}{TP+FP}$. (1st) and (2nd) refer to the iterations, the process being continued until an evaluation precision of 75% is reached.

airgun surveys. Moreover, TiSSNet was mainly trained as a T-phase detector, while the two signals targeted here have specific signatures that may not be fully captured by TiSSNet. To improve its detection capacity for these new classes, an active fine-tuning approach was adopted. The uncertainty sampling method (e.g., Mosqueira-Rey et al., 2023) was applied on the TiSSNet probability outputs to create a new training data set from the MAHY data and iteratively fine-tuned TiSSNet on these new classes. The same methodology was applied for our two target classes, independently of each other, thus enabling the development of two new detection models called “H-TiSSNet” (fine-tuned for H-phase waves) and “P-TiSSNet” (fine-tuned for P-phase waves). In more details, TiSSNet outputs, for each time bin of the input spectrograms, an estimate of the probability *p* that this time bin is part of a geophysical acoustic signal. To actually obtain a binary classification, a threshold *χ* is required such that results with *p* > *χ* are considered as positive and negative otherwise. With a threshold of *χ* = 0.5, uncertain positive classifications correspond to probabilities slightly above 0.5, while uncertain negative classifications are beneath. A total of 5,000 TiSSNet detections were sampled with estimated probabilities lying between 0.4 and 0.6, and were manually either discarded or labeled as “positive” or “negative.” In the case of

H-waves, which turned out to be more ambiguous than P-phases, any impulsive-shaped signal which did not show a periodicity like air-gun shots for seismic exploration was classified as positive.

The resulting annotations were then used as a training data set to fine-tune TiSSNet. A 5-fold cross-validation with 80/20% split between the training and validation data sets was used. To finish, the precision of the fine-tuned model was evaluated by manually checking 1,000 more detections above the $\chi = 0.5$ threshold. Given the scarcity of the target signals in the data set, blind random sampling of 200-s windows is unlikely to contain P-phases or H-waves. This would lead to a particularly low and unreliable false negative rate estimate. This is why only the model precision is estimated for evaluation, leaving the recall to be estimated from cross-validation in the training data set. This was iteratively repeated, using up to 10,000 new samples for training and validation in further iterations. Iterations have been continued until the evaluation precision reached the arbitrary value of 75% with $\chi = 0.5$. This value has been deliberately kept fairly low to enable the association to carry out more detailed filtering of detections based on source locations at a later stage. While P-phases were globally well-detected by the original TiSSNet model, because of their similarity to T-phases, only a few H-waves were actually detected. In the case of H-TiSSNet, to bootstrap the process, the available manual annotations from Lavayssière et al. (2024) were added to the first training loop.

The training metrics are summarized in Table 1. After one iteration, P-TiSSNet reached the expected precision of 75%, whereas H-TiSSNet required an extra iteration to reach this score. As observed during the evaluation and data set building for both detectors, most false positives were either seismic airgun-shots, unidentified anthropic noises or marine mammal vocalizations.

## 2.2. Association Task

### 2.2.1. Problem

Let association be the process of grouping detections from several stations that correspond to the same source event. A similar problem arises in seismology, leading to the development of methods such as GAMMA (Zhu et al., 2022). However, our problem differs from the previous one in at least two respects: first, we expect at most one detection per source event and per station, and second, greater ambiguity is likely to occur due to the wider physical detection range of hydrophones. Furthermore, unlike some networks that use closely spaced hydrophone arrays (approximately 2 km apart) to take advantage of signal coherence, thereby enabling cross-correlation (Vergoz et al., 2021), this coherence is likely to be lost due to the distance between our receivers (~70 km apart).

The problem can be formulated as follows. Let $s \in S$ be a station, and let the detection step described in Section 2.1 result in a corresponding set of detection times $T_s$. An association *a* is defined as a set
$$a = \{t_s \in T_s | s \in S_a \subseteq S\},$$

where $S_a$ contains no duplicate stations. In other words, each association contains at most one detection $t$ from each station.

Given the full set of detection times $T = \{T_s | s \in S\}$, an association algorithm is expected to produce a set $A$ of physically plausible associations. An association $a \in A$ is considered plausible if there is a source position (in a 2D lat/lon plan) $(x, y)$ in the study area and an origin time $t_0$ such that

$$t_0 + \theta_{(x,y),s_i} \approx t_{s_i} \quad \forall t_{s_i} \in a,$$

where $\theta_{(x,y),s_i}$ is the sound propagation time from the position $(x,y)$ to station $s_i$.

Let us now consider the computational complexity of the problem. For a detection $t_{s_i}$ from station $s_i$, and assuming each other station $s_j \in S \setminus \{s_i\}$ has at least $n$ candidate detections, the number of associations containing $t_{s_i}$ is at least $(n + 1)^{|S| - 1}$. This leads to a lower bound complexity of $O(n^{|S|})$ for the total number of candidate associations.

However, the set of valid associations is constrained: a detection $t_{s_j}$ from station $s_j$ is a valid candidate to pair with $t_{s_i}$ if and only if

$$|t_{s_i} - t_{s_j}| = |\theta_{(x,y),s_i} - \theta_{(x,y),s_j}| < \theta_{s_i,s_j} \quad (\text{triangle inequality}) \quad (1)$$

where the left-hand side is the absolute value of the so-called Time Difference of Arrival (TDoA), and $\theta_{s_i,s_j}$ is the sound propagation time between stations $s_i$ and $s_j$. Assuming a uniform spatio-temporal distribution of sources, and that each source event is detected by all stations, the number of candidate detections verifying Equation 1 grows linearly with $\theta_{s_i,s_j}$. The value $\theta_{s_i,s_j}$ is assumed to grow linearly with the distance between $s_i$ and $s_j$.

Altogether, assuming that $\theta_{s_i,s_j} \ge \delta$ for all $s_i, s_j \in S$, with $s_i \ne s_j$, a naive association search algorithm has a lower bound complexity of $O(\delta^{|S|})$.

In practice, hydrophone networks can contain more than 10 stations spanning several thousand kilometers, like the OHASISBIO network (Royer, 2009). This corresponds to inter-station propagation times of up to 1 hr. Within such time windows, seismic swarms may generate hundreds of T-phases, leading to an infeasible large number of potential associations. To address this issue with an algorithm applicable to such large networks, we propose a scalable approach based on spatial and temporal progressive pruning, called TDoA Acoustic Progressive Algorithm for Association (TAPAAs).

### 2.2.2. TAPAAs

For each pair of stations $(s_i, s_j), i \ne j$, a geodesic grid of expected TDoA is computed for any source position $(x, y)$ in the study area. A non-rectangular grid is designed so that the orthodromic distance between two adjacent cells remains constant regardless of their latitude. This is important for large-scale station networks covering a wide range of latitudes.

The actual TDoA of an event within a given cell may differ from the predicted value for that cell. This is due both to the source not being centered in the cell, and to uncertainties in detection. To account for these effects, a total TDoA tolerance $\Delta$ is defined. It combines a geometrical tolerance $\Delta_g$ for the source position offset within the cell, with a pick-time tolerance $\Delta_p$, and a clock drift tolerance $\Delta_t$.

A cell $(x, y)$ should account for any detection of coordinates $(x', y')$ in the square defined by the corners $(x - \frac{l}{2}, y - \frac{l}{2})$ and $(x + \frac{l}{2}, y + \frac{l}{2})$. In addition, the acoustic model used to compute the TDoA values often implies an uncertainty in the sound-speed $c$ of upper bound $d_c$. Let $\delta_{i,j,x',y',c}$ be the predicted TDoA between $(s_i, s_j)$ in a cell $(x, y)$ with a sound-speed $c$, then the geometrical tolerance $\Delta_g$ can be empirically estimated by:
$$ \Delta_g = \max_{\substack{x' \in \left[x - \frac{l}{2}x + \frac{l}{2}\right] \\ y' \in \left[y - \frac{l}{2}y + \frac{l}{2}\right] \\ c' \in [c - d_c, c + d_c]}} \left| \delta_{s_i, s_j, x', y', c'} - \delta_{s_i, s_j, x, y, c} \right| $$

The pick time tolerance $\Delta_p$ is defined as the upper bound of the difference between used detection time and actual detection time. Similarly, a clock drift upper bound $\Delta_t$ is defined. The total tolerance of TDoA $\Delta$ is

$$ \Delta = \Delta_g + 2 * \Delta_p + 2 * \Delta_t $$

The effects of non-punctual sources are ignored and exact hydrophone coordinates are assumed. Given a pair of observed detections $t_{s_i}$, $t_{s_j}$, $(x, y)$ is considered a valid cell if and only if

$$ |t_{s_j} - t_{s_i} - \delta_{s_i, s_j, x, y, c}| \le \Delta \quad (2) $$

The values $\Delta$ are estimated for each cell of the grid.

We now consider detections from TiSSNet. Each detection $t$ is given with an estimated probability $p_t$. Let $T = \{t_{s_i} | s_i \in S, p_{t_i} > p_1\}$ with $p_1$ a meta-parameter serving as probability threshold. Similarly, $T_{s_i} = \{t_{s_i} | p_{t_i} > p_2\}$ with $p_2$ another metaparameter such that $p_2 < p_1$. Finally, let $G$ be the coordinates of the grid cells. The Algorithm 1 and its sub-routine Algorithm 2 gives a more detailed explanation of TAPAAs.

### Algorithm 1. TAPAAs

1: **for** $t_{s_i} \in T$ **do**

2: &nbsp;&nbsp;&nbsp;&nbsp;a $\leftarrow \{t_{s_i}\}$

3: &nbsp;&nbsp;&nbsp;&nbsp;valid $\leftarrow \{(x, y) | (x, y) \in G\}$

4: &nbsp;&nbsp;&nbsp;&nbsp;TAPAAs- SUB (a, valid, $\{s_i\}$)

5: **end for**

### Algorithm 2. TAPAAs-Sub

1: **Require:** a, valid, $S_{used}$

2: &nbsp;&nbsp;&nbsp;&nbsp;S$_{remain}$ $\leftarrow$ S $\setminus$ $S_{used}$

3: &nbsp;&nbsp;&nbsp;&nbsp;S$_i$ $\leftarrow$ first value of S$_{remain}$

4: &nbsp;&nbsp;&nbsp;&nbsp;TAPAAs- SUB (a, valid, $S_{used} \cup \{s_i\}$)

5: &nbsp;&nbsp;&nbsp;&nbsp;candidates $\leftarrow \{t_{s_i} \in T_{s_i} | \forall t_{s_j} \in a, |t_{s_i} - t_{s_j}| < \theta_{s_i, s_j}\}$

▷ Try a without $s_i$

▷ Equation 1

6: **for** $t_{s_i} \in$ candidates **do**

7: &nbsp;&nbsp;&nbsp;&nbsp;valid' $\leftarrow$ valid $\cap$ $\{(x, y) \in G | \text{Eq. 2 satisfied for } t_{s_i} \text{ and each } t_{s_j} \in a\}$ ▷

Equation 2

8: &nbsp;&nbsp;&nbsp;&nbsp;a' $\leftarrow$ a $\cup$ $\{t_{s_i}\}$

9: &nbsp;&nbsp;&nbsp;&nbsp;Save a', valid'

10: &nbsp;&nbsp;&nbsp;&nbsp;TAPAAs- SUB (a', valid', $S_{used} \cup \{s_i\}$)

11: **end if**

12: **end for**

Associations that include a detection that is already searched in the main part of the algorithm are pruned, as these associations have already been tested. Additionally, for any detection $t$, given the associations $A$ that contain it, $a \in A \Rightarrow |a| \ge |a'| - n \quad \forall a' \in A$, with $n$ being a meta-parameter used to save memory. It is set to two for the
rest of this work. Finally, only associations of size three or greater are saved, as it is not possible to localize a
source using less than three stations with no prior knowledge. The algorithm produces redundant associations;
that is, for two sets of associations $a_1$ and $a_2$, it is possible that $a_1 \cap a_2 \neq \emptyset$. Posterior filtering steps are expected to
resolve such redundancies.

# 3. Application of the Processing Pipeline for Submarine Volcano Monitoring

## 3.1. Methodology

P-TiSSNet was first applied on the MAHY data, selecting detections with a probability output of more than 0.05. The resulting detections had, at most, the time-detection accuracy of TiSSNet, which is ~0.53 s (due to the resolution of spectrograms). Since P-phases mainly travel as seismic waves (at ~5.5 km/s), this value represents an inaccuracy of at least 3 km, which is not satisfactory for our purpose. To address this issue, given that P-phases usually have a very steep start in the time domain, a re-picking step was applied before the association. To this effect, a 6 s waveform centered on the detected signal was band-pass filtered between 4 and 60 Hz, then squared and averaged in a 0.1 s sliding window; the first occurrence of a value greater than the fifth of the maximum was then selected as the P-phase arrival time. An example is given Figure S2 in Supporting Information S1.

T-phases are generally supposed to originate from a single point on the ocean-bottom, corresponding to the epicenter for shallow earthquakes, enabling a 2D location. However, P-phases directly come from the hypocenter and the TDOA at different stations will greatly depend on its depth. A 3D version of TAPAAs, 3D-TAPAAs, was then used. A rectangular grid between the corners (13.4°S, 45.25°E) and (12.4°S, 46.25°E) with a cell side of 1.11 km was created, enclosing the MAHY network. Vertically, the grid was extended down to -100 km with a resolution of 1 km, given that local seismicity occurs at a depth of approximately 40 km (Retailleau et al., 2023). Overall, the 3D grid comprises about a million cells. The propagation time from each cell to the ocean-bottom beneath each hydrophone was estimated using the 1D seismic velocity model from Lavayssière et al. (2022) and TauP (Crotwell et al., 1999). The acoustic propagation time between the seafloor and each hydrophone was then estimated, assuming that the acoustic waves are generated at the foot of the AuH mooring line, by dividing the length of the line by a fixed sound-speed of 1,490 ms⁻¹. To finish, we used as metaparameters $p_1 = 0.3$, $p_2 = 0.15$ and $\Delta_p = 1$ s. To evaluate the accuracy of the cataloging pipeline, the TAPAAs catalog was compared to the catalog of earthquakes from Lavayssière and Retailleau (2023) until March 2021 and then from the REVOSIMA observatory data (Retailleau et al., 2022) until October 2024. These catalogs were compiled from data recorded by land-based seismometers, processed manually for the first and automatically for the second. A first iteration of association was applied to fix the clock drift of the sensors, using uncorrected detection dates, that is which suffer from clock drift. $\Delta_t = 2$ s was fixed. The resulting associations were compared to the reference catalog, considering an association matches a reference event if it occurs within 15 s of it. For each deployment MAHYij, the difference between the expected and actual times of arrival as a function of time was determined. The result is expected to be well-fitted by a linear regression, whose slope $d_{MAHY\ ij}$ corresponds to the clock drift of MAHYij and intercept $\Delta_{MAHY\ *j}$ corresponds to some propagation uncertainties resulting, among other causes, from an inaccurate seismic velocity model. Linear regressions were then performed imposing a similar value for all stations with the same position $\Delta_{MAHY\ *j}$. The obtained values were used to linearly correct all the P-phase and H-wave detection times from both the clock drift and the intercept values, before a second iteration of association with $\Delta_t = 0.1$ s was used for P-phase analysis.

Once the second iteration of P-phase association was run, a final filtering was applied. All associations locating an event outside the MAHY network defined by [12.95, 12.65]°S and [45.27, 45.67]°E were initially discarded. To eliminate any duplicates, any association of size three that was a subset of an association of size four was discarded, and for all associations sharing a detection, the associations with the lowest TAPAAs residues were selected.

The Received Levels (RL) of each event on each AuH were then computed as the root mean square of the signal in the 30 s following the P-phase detection times. We compared those values with earthquake magnitudes $M_L$ (local scale) for each station using the association matched with reference events. The distribution of the RL values of a single AuH was then plotted for the generated catalog in a Gutenberg-Richter fashion. This enabled to derive a slope, comparable to a b-value, that we name $b_{RL}$ to avoid ambiguity. To compute this value in presence of a nonlinear tail, we use the truncated Gutenberg-Richter technique of Kagan (2002), assigning a fixed level of
completeness $RL_c$ and corner level $RL_T$. We performed this analysis for the two seismic clusters independently, defining the proximal group as being west of the 45.5°E meridian and the distal group as being east of it. While the obtained $RL$ values completely ignore the transmission losses, we assume the losses are similar within the same cluster. Within a cluster, the transmission losses are thus expressed as a constant offset, and do not affect $b_{RL}$. This causes the obtained $b_{RL}$ values independently computed in the two clusters to be consistent and comparable, while a $b_{RL}$ value computed on the whole catalog may be unreliable. The evolution of $b_{RL}$ during the recording time was estimated using the method of Jacobs and McNutt (2010). A sliding window of 400 events with an overlap of 50% was chosen, and a succession of $b_{RL}$ over the data set were estimated.

For H-wave association, the 2D-TAPAAs algorithm was applied with the same meta-parameters as for 3D-TAPAAs, with a cell size of 750 m, and metaparameters set to $p_1 = 0.5$ and $p_2 = 0.35$. The propagation times were estimated using great-circle distances and a fixed velocity of $1,490 \text{ ms}^{-1}$, consistent with the local SOFAR velocity as shown in Figure S1 in Supporting Information S1. The catalog of Lavayssière et al. (2024) was used as ground-truth, but it must be reminded that these annotations were used for H-TiSSNet training. Before the comparison, the same filtering process as for P-phase associations was used for H-waves, in the rectangle bounded in [12.9, 12.8]°S and [45.6, 45.75]°E. Finally, the $RL$ were computed as peak amplitudes in a 2 s window centered on the signals, and converted to SL using a cylindrical transmission loss correction such as that described in Royer et al. (2015), except that the vertical component of the ray path was neglected.

## 3.2. Results

The application of P-TiSSNet and H-TiSSNet resulted in 227,012 and 2,891,241 automatic detections within the 4 years of continuous data, as shown in Figure 4. While the number of P-phases rapidly decreased, following the trends shown by the reference seismic catalog, H-waves remained frequent during the whole recording period. This is due to the existence of many impulsive signal sources, resulting in false positives that may not be discardable at this stage.

The first iteration of P-phase association using 3D-TAPAAs resulted in 28,194 associations of size three or more, of which 10,744 remained after geographical filtering and duplicate removal. Only 2,847 of them matched an earthquake from the reference catalog. This represents 42.7% of the 6,666 events in the reference catalog happening in the period of study. The differences between expected and observed arrival times allowed us to estimate the clock drifts of each hydrophone, which ranged from -0.065 to 0.182 ppm, with intercept from 0.08 to 1.64 s. Figure S3 in Supporting Information S1 shows an example of the evolution of this difference over the recording time of MAHY21. For 13 deployments, out of 19, the clock drift was also successfully estimated from the post-deployment GPS synchronization (upon deployment all hydrophones clocks are synchronized with the GPS clock). The two methods show a good general agreement (Table S1 in Supporting Information S1).

After correcting all the AuH detection times resulting from P-TiSSNet and H-TiSSNet, the more restrictive iteration of TAPAAs led to 26,018 P-phase associations of size three or more, of which 11,606 remained after geographical filtering and duplicate removal. Among them, 3,215 matched an earthquake in the reference catalog, which represent 48.2% of the 6,666 events in the reference catalog happening in the period of study. Figure 5 shows the time distribution of the association catalog and compares it to the time distribution of the reference catalog. One can notice expected decreases due to AuH absence. The sudden increase in the number of reference events in March 2021 is an artifact, as it corresponds to the start of the application of automatic methods to create the REVOSIMA catalog.

Figure 6 shows the spatial distribution of the P-phase associations resulting from 3D-TAPAAs application. We observe a striking concordance between the catalog resulting from our automatic processing pipeline for continuous acoustic data (AuH) and the reference catalog obtained with the regional seismological network. The AuH network is also able to distinguish the seismic clusters: the proximal and distal one with its elongated shape.

Regarding the H-waves, TAPAAs produced 10,880 associations of size three or more, of which 1,093 remained after geographical filtering and duplicate removal. Most of the prefiltering associations correspond to false positive and are distributed randomly in space. On the other hand, the overall spatial distribution of associations clearly highlight two groups, as shown in Figure 6, which is likely caused by true positives. This makes the impact of geographic filtering particularly significant. Once the geographical filter applied, a total of 470 events in the catalog of Lavayssière et al. (2024) were detected out of 585, that is 80.3%. The temporal distribution of
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Weekly chronograms of the number of automatic detections in function. Left: P-phase detections. Right: H-wave detections. The light, four first rows correspond to detections obtained from fine-tuned *TiSSNet* model on each of the four MAHY stations. The darker last row corresponds to the filtered associations resulting from the final application of *TAPAA*s. Gray rectangles account for periods of data gaps due to instrument recovery or failure. The Σ value corresponds to the total number of events in each subfigure.</figcaption></figure>

associations is shown in Figure 4 and accurately matches the distribution shown by the reference catalog. A closer look at the first 2 months of the distribution of these signals is given as Figure S4 in Supporting Information S1 and highlights the two lava vents that have been observed by the same work. Their spatial distribution is shown in Figure 6; the H-waves cluster occurs in the Tiktak area as previously observed by Lavayssière et al. (2024). A second cluster of automatic detections lies near MAHY*4. To understand the cause of such clustering, we examined all events comprised between [12.60, 12.45]°S and [45.45, 45.60]°E. As shown in Figure S5 in Supporting Information S1, the 807 selected events mostly occurred during working hours, and are therefore likely of anthropic origin. This, combined with our knowledge of mining on the island of Mayotte, suggests that these are blasts in a quarry. The reason why they are located near MAHY*4, however, remains unknown and may
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Top: Time distribution of the automatic associations computed with a sliding-window of 10 days and 25% overlap. The red line on 01/05/21 separates Phase I from Phase II, characterized by two different levels of seismic activity detected by the AuH. Gray rectangles show periods with limited AuH availability. Bottom: Time distribution of earthquakes from the reference catalog computed with the same sliding-window. The bisque rectangle marks the time before the automatic cataloging method adopted by REVOSIMA analysts. The considered reference catalog was limited here to the AuH availability periods. Dashed gray vertical lines mark earthquakes with $M_L > 4$, and black ones, two earthquakes felt by the population of Mayotte Island.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Spatial distribution of the associations resulting from TAPAA s. Blue and green cells represent P-phase automatic detection and earthquake epicenters while purple cells represent H-waves automatic locations. The blue/green delimitation represents the proximal/distal cluster delimitation and was fixed to 45.5°E. The corresponding earthquake and H-wave source ground truths are represented with contours. The outer contours correspond to 10 earthquake occurrences, while the black inner contour is for 30 earthquakes. Grid cells with less than 200 H-waves or 1,000 P-phase associations were left blank. “Petite Terre,” contoured in gray, is the closest island from the Mayotte Archipelago. The red triangle represents Fani Maoré main edifice.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Distributions of errors for the automatic locations with respect to the reference seismic catalog. An error of 0.2° in position corresponds to ~1 km. Red, dashed lines and their labels correspond to the 5% and 95% quantiles.</figcaption></figure>

be due to an efficient seismic to acoustic conversion on the island submarine slopes. Figure S6 in Supporting Information S1 shows an example of such a signal corresponding to a day of known blasting activity in a quarry.

The differences between the chosen P-phase association locations and the reference catalog are shown in Figure 7, selecting the closest association in space in case of ambiguity. It has to be noted that origin times from the catalog of Lavayssière and Retailleau (2023) have a precision of 1 s, thus increasing the standard deviation in the time difference. Interestingly, 10,360 refined associations were not associated with any event from the reference catalog. It is speculated that most of them correspond to badly constrained events that were discarded while compiling the reference earthquake catalog.

The left-hand part of Figure 8 shows *RL* values as a function of *M*<sub>L</sub> values from the reference catalog, for all stations MAHY*3. MAHY*3 was chosen because they provide the best linear fit with local magnitudes, as shown later. The obtained relationship for each station is given in Table S2 in Supporting Information S1. Similarly, Figure S7 in Supporting Information S1 displays the same plot for all stations and for three different *RL* computation techniques, showing the technique exposed here led to the best linear fit. The right-hand part of Figure 8 displays the Gutenberg-Richter distribution, relative to *RL* values, for the 9,905 MAHY*3 P-phase associations. Similar curves, relative to *M*<sub>L</sub> values, are plotted for the whole reference catalog, and its subset matching the AuH associations.

The Gutenberg-Richter curve of the TAPAA*s* association resulting catalog (blue curve) shows a linear behavior from a conservative value *RL*<sub>*T*</sub> = 140 dB up to the completeness level *RL*<sub>*c*</sub> = 113 dB, which may be compared to *M*<sub>*L*</sub> = 1.35 in the local magnitude scale. The slope of the power-law distribution yields a *RL* b-value *b*<sub>*RL*</sub> = 0.052 for the 9,905 MAHY*3 associations. The usual ratio between small earthquakes and large ones in a tectonic context generally implies a seismic b-value of 1, which would correspond to a *b*<sub>*RL*</sub> of 0.070 here. Although, b-values are usually larger than 1 in volcanic regions, we observe that *b*<sub>*RL*</sub> is lower than expected. For comparison, the *b*-value obtained from the reference catalog (black curve) using local magnitudes, with a magnitude of completeness of *m*<sub>*b*</sub> = 1 and no corner magnitude, is 0.569 while it is 0.653 for the subset matching the AuH
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Left: comparison of $M_L$ values with the estimated RL Received Levels for all MAHY\*3 detections belonging to associations matching an event in the reference seismic catalog. Right: Gutenberg-Richter relationship linking a value of RL or $M_L$ to the number of events with a greater magnitude. The orange curve corresponds to the whole reference catalog and the gray curve to the subset matching AuH associations, both plotted against their $M_L$ values. The blue curve is the distribution of the whole MAHY\*3 association catalog (TAPAAs), plotted against RL values. Numbers in parentheses are the number of events in each curve. Matching events yield the correspondence between RL and $M_L$. The blue vertical dashed lines show the level of completeness $RL_c$ and the corner level $RL_T$ used for the Gutenberg-Richter fit of the blue curve.</figcaption></figure>

automatic detections. The $b_{RL}$ of the Gutenberg-Richter relation was further examined as a function of year (Figure 9). It shows that, considering the whole association catalog, $b_{RL}$ slowly decreases over the years similarly as the number of events per year. We can take advantage of our 4-year long catalog obtained in a consistent way with our new processing suite to further investigate the evolution of the weak seismicity over-time and its effect on $b_{RL}$. Applying a sliding event-window to estimate the $b_{RL}$, as shown in the right part of Figure 9, confirms the previous statement and shows a globally decreasing $b_{RL}$ over time. The behavior of $b_{RL}$ can be visually separated in two distinct phases (red line in the Figure). From the beginning of the recordings until September 2021, the proximal cluster data show a high and variable $b_{RL}$ value, with a mean of 0.090 and a standard deviation of 0.012, falling to 0.075 with a slightly lower standard deviation of 0.010. In comparison, $b_{RL}$ in the distal cluster falls from 0.040 to 0.029, and is more than twice smaller and more stable than in the proximal cluster. It is worth noting that the Gutenberg-Richter distribution for the distal cluster seismicity is less linear than that for the proximal cluster.

## 4. Discussion

The absence of T-phases in the inspected AuH signals is the first unexpected result. The combination of the great depth of events with the local seismic-to-acoustic conversion setup may have hindered the excitation of T-phases. Okal (2008) predicts that deep earthquakes should result in emergent T-phases, as they result from the integration of many body-wave conversions. In the MAHY framework, we guess that the AuHs are too close to integrate enough such waves, resulting in a low-energy diluted T-phase. This near-field hypothesis in the presence of moderately deep earthquakes needs to be tested by modeling both seismic and hydroacoustic waves using tools such as spectral elements (Komatitsch et al., 2024).

The results from the association pipeline, in particular the automatic localization of hypocenters, closely matches the reference catalog. TAPAAs was successfully extended to three dimensions. To our knowledge, this is the first attempt to localize hypocenters with moored AuHs. Until now, they were mainly used to monitor regional seismic activity in 2-D from T-phases, assimilating acoustic radiators with epicenters, for shallow earthquakes. T-phase monitoring enhances the detection of low-magnitude events with a completeness of the order of $m_b = 3.0$ (Dziak et al., 2004) or 3.4 (Ingale et al., 2021) depending on the array dimension. This work demonstrates the capability of a hydrophone network, on a scale of approximately 100 km or less, to estimate earthquake depths and to lower the local magnitude of completeness to ~1.35. However, it must be emphasized that the RL-$M_L$ linear conversion only has a determination coefficient of 0.76. However, even if it applies fairly well to the proximal and distal
<figure><img src="image_11.png" /><figcaption>**Figure 9.** RL Gutenberg-Richter distributions and *b* values by year for the whole TAPAAs association catalog (top), and for the catalogs corresponding to the proximal (middle) and distal (bottom) clusters (as delineated in Figure 6). Left: RL Gutenberg-Richter distributions by yearly catalogs (in different colors); numbers in parenthesis are the number of events per catalog. Right: *b*<sub>RL</sub> of the Gutenberg-Richter distributions calculated in sliding windows of 400 events with a 50% overlap. Vertical error bars show the 95% confidence interval of the regression. The mean *μ* and standard deviation *σ* of *b*<sub>RL</sub> are given for the two Phases of seismic activity detected by the AuH (separated by the vertical red line). The horizontal dashed line at 0.070 represents the equivalent seismic magnitude *b*-value of 1 converted to RL using the relationship in Figure 8. All RL values are from MAHY*3 data.</figcaption></figure>

clusters, it would be difficult to generalize it to other contexts where transmission losses in the crust might be
different. Still, such detection capabilities are of interest for marine seismic observatories without OBS. This calls
for further studies with AuH deployments at similar scales. It must be noted, though, that nearly 50% of the events
in the reference catalog were missed by the AuH's, likely because of magnitudes lower than the magnitude of
completeness or the randomness of the detection capabilities in a noisy acoustic environment. For example,
Mayotte Island generates an important marine traffic, producing a significant low-frequency noise. Conversely,
more than 50% of the events located with AuH were not present in the reference catalog, which may reflect its
incompleteness—or, less likely, false positives on the AuH side. Regarding clock drift corrections, a local
reference catalog is not necessarily needed to calculate it; distant teleseismic P-waves can be used instead (e.g.,
Safran et al., 2025; Sukhovich et al., 2021).

The discrepancies between the automatic AuH catalog and the seismic reference catalog are quantified in terms of
location and time differences in Figure 7. Latitude and longitude difference distributions are almost centered
gaussian, while depth and time errors reveal a bias. This might be caused by an overrated seismic velocity which would shift the poorly constrained depths downward. The use of a 3D-crustal velocity model could account for inhomogeneities in the crust. S-phases, and other secondary phases, may also be of interest and precisely detected with tools similar to PhaseNet (Zhu & Beroza, 2018) and stochastically associated with tools such as GAMMA (Zhu et al., 2022). Moreover, the double difference method (Zhang, 2003) would eliminate inaccuracies resulting from clock drift or local velocity heterogeneities, and should produce better relative locations.

The RL sound levels of received P-phases correlate well with the local magnitudes of the sources. However, the relation linking both values is not trivially linear, as shown by the obtained R² value of 0.76. Being aware of this limitation, M_L values obtained by hydrophones may be relevant for large statistical analysis like in the present work. The associated catalog was complete enough to estimate b_RL, comparable to standard b-values, and allowed us to observe a decrease in b_RL over time. Jacobs and McNutt (2010) noted higher b-values as Etna approached its eruptive phase and hypothesized that this was due to the presence of magma at a shallower depth. A less expected result, however, is the important difference between the two clusters, with a high b_RL for the proximal cluster events, which means a high amount of small events despite the rapid decrease of activity during the recording period. The distal cluster, on the other hand, shows a smaller and more stable value of b_RL after a rapid decrease. This demonstrates that the two systems are different. The proximal cluster has been proposed to delimit the shape of a ~10 × 10 × 5 km reservoir acting like a piston structure (Jacques et al., 2024) while the distal cluster is related to a dike propagation towards the new volcanic edifice (Lavayssière & Retailleau, 2023). Consequently, the proximal cluster is probably long-lasting when the distal one is more transient. The b-value is often thought to be linked to the stress regime, degassing or magma transfer (Firetto Carlino et al., 2025; Jacobs & McNutt, 2010). Two types of monitoring were maintained to track magmatic fluid emissions in the area: a submarine glider has been measuring CO₂ and CH₄ in the Horseshoe area near the proximal cluster since 09/2021 (Heumann et al., 2025) and a land station has been measuring CO₂ concentration in the air near the airport on the east coast of the island since 11/2020 (Bénard et al., 2023; REVOSIMA, 2025). It is likely that degassing fluctuates in response to magmatic activity and it consequently modifies the seismogenic volume within the proximal cluster. We propose that the variations we observe in the seismic b-values over a time-scale of several months (Figure 9) are related to degassing fluctuations near the proximal cluster. In this fractured and fluid-rich medium, seismicity occurs preferentially in low-magnitude earthquakes, resulting in a higher b-value than in the distal cluster. Those values of b_RL may however be affected by some saturation-like phenomena: the relation linking seismic moments to acoustic levels are not well known and may present nonlinearities causing high or low moments to be underrepresented by a Gutenberg-Richter relation drawn using received levels. Our observations highlight the potential of AuH for magnitude and b-value analyses, and calls for attempts to better estimate magnitudes from acoustic levels. It is important to note that P-phases were used in this work, and that T-phases are likely to yield much different results.

H-waves are harder to discriminate, leading to many false positive detections that actually look similar to the expected signals. Figure 4 shows the high amount of false positive resulting from the detection process, with many detections after the presumed end of the eruption and incompatible with the expected locations. However, the last histogram in the figure shows the distribution of the associations from TAPAAs looks like that of an eruption cluster. In particular, Figure S6 in Supporting Information S1 shows the distribution of spatially filtered associations closely matches the reference catalog and thus leading to the same conclusion as that of Lavayssière et al. (2024), namely the end of the eruption on 4 December 2020. This highlights the high potential of TAPAAs to filter noisy data, making it a relevant tool for statistical analysis of large data sets.

The application of the newly developed pipeline should facilitate volcanic monitoring and, consequently, risk analysis by the REVOSIMA observatory. The tuned *H-TiSSNet* algorithm will be applied on the upcoming and future MAHY data sets with the objective to rapidly detect new lava flows. *H-TiSSNet*, *P-TiSSNet* and their processing by *TAPAAs* were all run on a single laptop, with a total run time of nearly 10 hr for the whole MAHY data set (4 quasi-continuous years from 4 hydrophones). This represents an average of 25 s per day of data, or 2h30 to process 1 year of MAHY data after recovery.

**5. Conclusions**

In this study, we presented a new algorithm to associate automatic detections identified on continuous hydroacoustic data. The proposed pipeline allows to build catalogs of geophysical events generated by two different
processes—earthquakes and lava-seawater interactions—and displaying two different signatures—acoustic P-phases and H-phases. This allowed us to establish a multi-year overview of the seismic activity (P-phases) and of the occurrence of submarine lava flows (H-waves) as recorded by a network of four hydrophones moored in the SOFAR channel. Its purpose was to monitor the submarine seismo-volcanic activity off-shore Mayotte in the Northern Mozambique Channel, as part of the REVOSIMA observatory. Its low level of completeness has made it possible to catalog a sufficiently large number of earthquakes to allow for a realistic estimate of regional b-values. Our results show a change in the regime ~6 months after the end of the eruption: the initial phase I presented high $b_{RL}$ with large temporal variability while phase II presented stable and low $b_{RL}$. Moreover, the proximal cluster displayed a much higher $b_{RL}$, potentially associated with magmatic and degassing processes, than the distal cluster, which may be linked to fracturion associated with the diking event. In addition, the new pipeline will enable the design of a fast and lightweight tool for an automatic analysis of impulsive events in the event of renewed lava flows.

## Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

## Data Availability Statement

All marine operations of REVOSIMA are performed as part of the MAYOBS set of cruises (https://doi.org/10.18142/291) and contribute to the French Service National d'Observation en Volcanologie (SNOV). The MAHY data is available as FDSN network 1T (https://doi.org/10.15778/RESIF.1T2018). Data used to train *TiSSNet* partly originate from annotations of Lavayssière et al. (2024) and are given together with the trained model weights and the results (Raumer, 2025). The code used in this work was based on the hydroacoustic toolbox (Raumer & Safran, 2025). The reference earthquake catalog used in this work was made from the one available at IPGP website (https://ws.ipgp.fr/fdsnws/event/1, for more information see the website http://volobsis.ipgp.fr/revosima/) and the supplementary materials of Lavayssière and Retailleau (2023). In addition to GEBCO data, Mayotte bathymetric data may be found at https://sextant.ifremer.fr/Donnees/Catalogue.

## Acknowledgments

The hydrophones used to collect the data were funded by the *Contrats de projet Etat-Région Bretagne (CPER ODO and O3DO)* and are now part of the Epos-France Mobile Marine Seismology action (smm.epos-france.fr), which is partially funded by the MARMOR project (ANR 21-ESRE-0020). We thank the captains and crews of R/V Marion Dufresne and Pourquoi Pas? for the successful deployments and recoveries of the MAHY network, and Jonathan Tanrin and Mickael Beauverger for instrument operations. We also thank Lise Retailleau for the fruitful discussions regarding the public earthquake catalogs published by REVOSIMA, and the REVOSIMA itself. Pierre-Yves Raumer was supported by a PhD fellowship from the University of Brest and from the Regional Council of Brittany (ARED), through the Interdisciplinary Graduate School for the Blue Planet (ISblue), co-funded by ANR (ANR-17-EURE-0015). We would also like to thank the anonymous reviewer for their very helpful comments.

## References

* Bazin, S., Raumer, P.-Y., Lavayssière, A., Fauvel, H., Safran, R., Ingale, V. V., et al. (2025). Monitoring underwater eruptions using networks of hydrophones in the SOFAR channel. In *Proceedings of UACE 2025*.

* Bazin, S., Royer, J.-Y., Dubost, F., Paquet, F., Loubrieu, B., Lavayssière, A., et al. (2023). Initial results from a hydroacoustic network to monitor submarine lava flows near Mayotte Island. *Comptes Rendus Geoscience*, **354**(S2), 257-273. https://doi.org/10.5802/crgeos.119

* Bénard, B., Di Muro, A., Liuzzo, M., Gurrieri, S., Boissier, P., Brunet, C., et al. (2023). Seasonal environmental controls on soil CO₂ dynamics at a high CO₂ flux sites (Piton de la Fournaise and Mayotte Volcanoes). *Journal of Geophysical Research: Biogeosciences*, **128**(6), e2023JG007409. https://doi.org/10.1029/2023JG007409

* Crotwell, H. P., Owens, T. J., & Ritsema, J. (1999). The TauP Toolkit: Flexible seismic travel-time and ray-path utilities. *Seismological Research Letters*, **70**(2), 154-160. https://doi.org/10.1785/gssrl.70.2.154

* De Groot-Hedlin, C. D. (2005). Estimation of the rupture length and velocity of the Great Sumatra earthquake of Dec 26, 2004 using hydroacoustic signals. *Geophysical Research Letters*, **32**(11), 2005GL022695. https://doi.org/10.1029/2005GL022695

* Dziak, R. P., Bohnenstiehl, D., Matsumoto, H., Fox, C., Smith, D., Tolstoy, M., & Fowler, M. (2004). P- and T-wave detection thresholds, Pn velocity estimate, and detection of lower mantle and core P-waves on ocean sound-channel hydrophones at the mid-Atlantic ridge. *Bulletin of the Seismological Society of America*, **94**(2), 665-677. https://doi.org/10.1785/0120030156

* Feuillet, N., Jorry, S., Crawford, W. C., Deplus, C., Thinin, I., Jacques, E., et al. (2021). Birth of a large volcanic edifice offshore Mayotte via lithosphere-scale dyke intrusion. *Nature Geoscience*, **14**(10), 787-795. https://doi.org/10.1038/s41561-021-00809-x

* Firetto Carloino, M., Lo Bue, R., Taroni, M., Coltelli, M., Cannavò, F., Cocina, O., et al. (2025). Earthquake frequency-magnitude distribution at Mount Etna sheds light on magma ascent in the volcano's plumbing system. *Science Advances*, **11**(41), eadx9873. https://doi.org/10.1126/sciadv.adx9873

* Fox, C. G., Matsumoto, H., & Lau, T. A. (2001). Monitoring Pacific Ocean seismicity from an autonomous hydrophone array. *Journal of Geophysical Research*, **106**(B3), 4183-4206. https://doi.org/10.1029/2000JB900404

* Gros-Martial, A., Delaigue, L., Cauchy, P., Pensieri, S., Bozzano, R., Bonnel, J., et al. (2025). Benchmarking models for ocean wind speed estimation based on passive acoustic monitoring. In *Oceans 2025 brest* (pp. 1-10). https://doi.org/10.1109/OCEANS58557.2025.11104747

* Heumann, A., Margirier, F., Rinnert, E., Lherminier, P., Scalabrin, C., Géli, L., et al. (2025). A 30-month data set of glider physico-chemical data off Mayotte Island near the Fani Maoré volcano. *Earth System Science Data*, **17**(9), 4535-4554. https://doi.org/10.5194/essd-17-4535-2025

* Ingale, V. V., Bazin, S., & Royer, J.-Y. (2021). Hydroacoustic observations of two contrasted seismic swarms along the southwest Indian Ridge in 2018. *Geosciences*, **11**(6), 225. https://doi.org/10.3390/geosciences11060225

* Ingale, V. V., Parnell-Turner, R., Fan, W., Talling, P. J., & Neasham, J. (2025). Hydroacoustic signals recorded by CTBTO network suggest a possible submarine landslide in Trou Sans Fond Canyon, Offshore Ivory Coast, March 2024. *Seismological Research Letters*, **96**(4), 2253-2263. https://doi.org/10.1785/0220240448
Jacobs, K. M., & McNutt, S. R. (2010). Using seismic b-values to interpret seismicity rates and physical processes during the preeruptive earthquake swarm at Augustine Volcano 2005-2006. In *Chapter 3 in The 2006 eruption of Augustine Volcano, Alaska (Professional Paper No. 1769)*. U.S. Geological Survey. (Series: Professional Paper).

Jacques, E., Hoste-Colomer, R., Feuillet, N., Lemoine, A., Van Der Woerd, J., Crawford, W., et al. (2024). Ring faulting and piston collapse in the mantle sustained the largest submarine eruption ever documented. *Earth and Planetary Science Letters*, **647**, 119026. https://doi.org/10.1016/j.epsl.2024.119026

Jorry, S., Paquet, F., & Lebas, E. (2022). MAYOBS23 cruise, RV Marion Dufresne. https://doi.org/10.17600/18002494

Kagan, Y. Y. (2002). Seismic moment distribution revisited: I. Statistical results: Seismic moment distribution: I. *Geophysical Journal International*, **148**(3), 520–541. https://doi.org/10.1046/j.1365-246x.2002.01594.x

Komatitsch, D., Tromp, J., Garg, R., Gharti, H. N., Nagaso, M., Oral, E., et al. (2024). SPECFEM/specfem3d: SPECFEM3D v4.1.1. Zenodo. https://doi.org/10.5281/zenodo.10823181

Komorowski, J.-C., & Paquet, F. (2024). MD246 / MAYOBS30 cruise, RV Marion Dufresne. https://doi.org/10.17600/18003592

Lavayssière, A., Bazin, S., & Royer, J.-Y. (2024). Hydroacoustic monitoring of Mayotte submarine volcano during its eruptive phase. *Geo-science*, **14**(6), 170. https://doi.org/10.3390/geosciences14060170

Lavayssière, A., Crawford, W. C., Saurel, J.-M., Satriano, C., Feuillet, N., Jacques, E., & Komorowski, J.-C. (2022). A new 1D velocity model and absolute locations image the Mayotte seismo-volcanic region. *Journal of Volcanology and Geothermal Research*, **421**, 107440. https://doi.org/10.1016/j.jvolgeores.2021.107440

Lavayssière, A., & Retailleau, L. (2023). Capturing Mayotte's deep magmatic plumbing system and its spatiotemporal evolution with volcanotectonic seismicity. *Volcanica*, **6**(2), 331–344. https://doi.org/10.30909/vol.06.02.331344

Le Saout, M., Bohnenstiehl, D. R., Paduan, J. B., & Clague, D. A. (2020). Quantification of eruption dynamics on the north rift at axial seamount, Juan de Fuca Ridge. *Geochemistry, Geophysics, Geosystems*, **21**(9), e2020GC009136. https://doi.org/10.1029/2020GC009136

Mosqueira-Rey, E., Hernández-Pereira, E., Alonso-Ríos, D., Bobes-Bascarán, J., & Fernández-Leal, A. (2023). Human-in-the-loop machine learning: A state of the art. *Artificial Intelligence Review*, **56**(4), 3005–3054. https://doi.org/10.1007/s10462-022-10246-w

Murch, A. P., Portner, R. A., Rubin, K. H., & Clague, D. A. (2022). Deep-subaqueous impulsive volcanism at West Mata seamount, Tonga. *Earth and Planetary Science Letters*, **578**, 117328. https://doi.org/10.1016/j.epsl.2021.117328

Okal, E. A. (2008). The generation of T waves by earthquakes. *Advances in Geophysics*, **49**, 1–65. https://doi.org/10.1016/S0065-2687(07)49001-X

Olive, J., & Escartín, J. (2016). Dependence of seismic coupling on normal fault style along the Northern Mid-Atlantic Ridge. *Geochemistry, Geophysics, Geosystems*, **17**(10), 4128–4152. https://doi.org/10.1002/2016GC006460

Parnell-Turner, R., Smith, D. K., & Dziak, R. P. (2022). Hydroacoustic monitoring of seafloor spreading and transform faulting in the equatorial Atlantic Ocean. *Journal of Geophysical Research: Solid Earth*, **127**(7), e2022JB024008. https://doi.org/10.1029/2022JB024008

Porter, M. B. (2011). The BELLHOP manual and user's guide: Preliminary draft (Tech. Rep.).

Raumer, P.-Y. (2025). PYLRR/toolbox-MAHY: 2025. https://doi.org/10.5281/zenodo.17531175

Raumer, P.-Y., Bazin, S., Cazau, D., Ingale, V. V., Lavayssière, A., & Royer, J.-Y. (2023). Application of machine learning to hydro-acoustic seismic and magmatic events detections. https://doi.org/10.5194/egusphere-egu23-7028

Raumer, P.-Y., Bazin, S., Cazau, D., Ingale, V. V., Royer, J.-Y., & Lavayssière, A. (2024). An open source hydroacoustic benchmarking framework for geophonic signal detection. *Seismica*, **3**(2). https://doi.org/10.26443/seismica.v3i2.1344

Raumer, P.-Y., & Safran, R. (2025). GEO-OCEAN-hydroacoustic/toolbox: 2025. https://doi.org/10.5281/zenodo.17531167

Reagan, J. R., Boyer, T. P., García, H. E., Locarini, R. A., Baranova, O. K., Bouchard, C., et al. (2024). World ocean atlas 2023.

Retailleau, L., Saurel, J.-M., Laporte, M., Lavayssière, A., Ferrazzini, V., Zhu, W., et al. (2023). Automatic detection for a comprehensive view of Mayotte seismicity. *Comptes Rendus. Géoscience*, **354**(S2), 153–170. https://doi.org/10.5802/crgeos.133

Retailleau, L., Saurel, J.-M., Zhu, W., Satriano, C., Beroza, G. C., Issartel, S., et al. (2022). A wrapper to use a machine-learning-based algorithm for earthquake monitoring. *Seismological Research Letters*, **93**(3), 1673–1682. https://doi.org/10.1785/0220210279

REVOSIMA. (2025). Mensual bulletin may 2025 [REVOSIMA Mensual Bulletins]. Retrieved from https://www.ipgp.fr/actualites-du-revosima/

Rinnert, E., Lebas, E., Komorowski, J.-C., Paquet, F., Jorry, S., Feuillet, N., et al. (2019). MAYOBS. Sismer. Retrieved 2025-06-20, from Retrieved from https://campagnes.flotteoceanographique.fr/series/291/doi:10.18142/291

Rinnert, E., Thinon, I., & Feuillet, N. (2020). MD 228 / MAYOBS15 cruise, RV Marion Dufresne. https://doi.org/10.17600/18001745

Rinnert, E., Thinon, I., & Lebas, E. (2021a). MAYOBS18 cruise, RV Pourquoi pas? https://doi.org/10.17600/18001984

Rinnert, E., Thinon, I., & Lebas, E. (2021b). MAYOBS21 cruise, RV Marion Dufresne. https://doi.org/10.17600/18001986

Royer, J.-Y. (2009). OHA-SIS-BIO—Observatoire hydroacoustique. Sismer. Retrieved from https://campagnes.flotteoceanographique.fr/series/229/doi:10.18142/229

Royer, J.-Y., Chateau, R., Dziak, R., & Bohnenstiehl, D. (2015). Seafloor seismicity, Antarctic ice-sounds, cetecean vocalizations and long-term ambient sound in the Indian Ocean basin. *Geophysical Journal International*, **202**(2), 748–762. https://doi.org/10.1093/gji/ggv178

Safran, R., Sukhovich, A., Perrot, J., Raumer, P.-Y., Bazin, S., Briais, A., & Royer, J.-Y. (2025). Seismo-acoustic analysis for accurate clock-drift correction of autonomous long-term hydrophones. In *Oceans 2025 brest* (pp. 1–7). https://doi.org/10.1109/OCEANS58557.2025.11104407

Saurel, J.-M., Retailleau, L., Deplus, C., Loubrieu, B., Pierre, D., Frangieh, M., et al. (2022). Combining hydro-acoustic sources and bathymetric differences to track the vent evolution of the Mayotte eruption, Mozambique Channel. *Frontiers in Earth Science*, **10**, 983051. https://doi.org/10.3389/feart.2022.983051

Schipper, C. I., & White, J. D. L. (2010). No depth limit to hydrovolcanic limu o Pele: Analysis of limu from Lō'ihi Seamount, Hawai'i. *Bulletin of Volcanology*, **72**(2), 149–164. https://doi.org/10.1007/s00445-009-0315-5

Slack, P. D., Fox, C. G., & Dziak, R. P. (1999). *P* wave detection thresholds, *Pn* velocity estimates, and *T* wave location uncertainty from oceanic hydrophones. *Journal of Geophysical Research*, **104**(B6), 13061–13072. https://doi.org/10.1029/1999JB900112

Smith, B., & Barclay, D. (2024). The sound of hydrothermal vents. In G. Bayrakci & F. Klingelhoefer (Eds.), *Geophysical monograph series* (1st ed., pp. 97–118). Wiley. https://doi.org/10.1002/9781119750925.ch7

Sukhovich, A., Irisson, J., Perrot, J., & Nolet, G. (2014). Automatic recognition of *T* and teleseismic *P* waves by statistical analysis of their spectra: An application to continuous records of moored hydrophones. *Journal of Geophysical Research: Solid Earth*, **119**(8), 6469–6485. https://doi.org/10.1002/2013jb010936

Sukhovich, A., Perrot, J., & Royer, J.-Y. (2021). Using teleseismic *P*-wave arrivals to calibrate the clock drift of autonomous underwater hydrophones. *Bulletin of the Seismological Society of America*, **111**(1), 21–35. https://doi.org/10.1785/0120200217

Tan, Y. J., Tolstoy, M., Waldhauser, F., & Wilcock, W. S. D. (2016). Dynamics of a seafloor-spreading episode at the East Pacific Rise. *Nature*, **540**(7632), 261–265. https://doi.org/10.1038/nature20116

Thinon, I., & Lebas, E. (2023). MD242 / MAYOBS25 cruise, RV Marion Dufresne. https://doi.org/10.17600/18003404
Tsang-Hin-Sun, E., Royer, J.-Y., & Perrot, J. (2016). Seismicity and active accretion processes at the ultraslow-spreading Southwest and intermediate-spreading Southeast Indian ridges from hydroacoustic data. *Geophysical Journal International*, *206*(2), 1232–1245. https://doi.org/10.1093/gji/ggw201

Vergoz, J., Cansi, Y., Cano, Y., & Gaillard, P. (2021). Analysis of hydroacoustic signals associated to the loss of the Argentinian ARA San Juan submarine. *Pure and Applied Geophysics*, *178*(7), 2527–2556. https://doi.org/10.1007/s00024-020-02625-7

Wang, P., Tan, Y. J., Bohnenstiehl, D. R., Wilcock, W. S., Tolstoy, M., Waldhauser, F., et al. (2025). Source mechanism of impulsive seafloor events that track submarine lava flows. *Science Advances*, *11*(34), eadk3942. https://doi.org/10.1126/sciadv.adk3942

Zhang, H. (2003). Double-difference tomography: The method and its application to the Hayward Fault, California. *Bulletin of the Seismological Society of America*, *93*(5), 1875–1889. https://doi.org/10.1785/0120020190

Zhu, W., & Beroza, G. C. (2018). PhaseNet: A deep-neural-network-based seismic arrival time picking method. *Geophysical Journal International*. https://doi.org/10.1093/gji/ggy423

Zhu, W., McBrearty, I. W., Mousavi, S. M., Ellsworth, W. L., & Beroza, G. C. (2022). Earthquake phase association using a Bayesian Gaussian mixture model. *Journal of Geophysical Research: Solid Earth*, *127*(5), e2021JB023249. https://doi.org/10.1029/2021JB023249