Alexey Sukhovich

Julie Perrot

Jean-Yves Royer

J.-Y Royer

Using teleseismic P-wave arrivals to calibrate the clock drift of autonomous underwater hydrophones

Keywords: Autonomous underwater hydrophone, clock drift, teleseismic P -wave

Networks of Autonomous Underwater Hydrophones (AUHs) are successfully employed for monitoring the low-level seismicity of mid-oceanic ridges by detecting hydroacoustic phases known as T -waves. For a precise localization of a seismic event from T -wave arrival times, all AUHs must be synchronized. To this effect, at the beginning of the experiment all instrument clocks are set to GPS time, which serves as a common reference. However, during the experiment the instrument clock often deviates from GPS time, and since the amount of deviation differs from one instrument to another, the synchronization of the AUHs deteriorates as the experiment progresses in time. Just after the instrument recovery, the time difference (called "skew") between the instrument and the GPS clocks is measured. Assuming that the skew varies linearly with time, the correction of a time series for the clock drift is a straightforward procedure. When the final skew cannot be determined, correcting for the clock drift is not possible and any event localization becomes problematic. In this paper, we demonstrate that the clock drift rate (assumed to be time-independent) can be successfully estimated from arrival times of teleseismic P -waves, commonly recorded by AUHs. Using a ray-tracing code, and accounting for the uncertainties in event hypocenter locations, origin times and the Earth seismic-velocity model, confidence intervals of the estimated drift rates are deduced. The validity of the approach is tested on data from two AUHs with known clock-drifts. Our results show that a reliable estimation is possible for skews as small as 4 s per two years (corresponding to a drift rate of about 5.5 ms • day -1 ). This method can also be applied to correct data of other recording instruments subject to internal-clock drift, such as ocean bottom seismometers, when the skew is unknown.

INTRODUCTION

As demonstrated by several studies, a network of Autonomous Underwater Hydrophones (AUHs) is an efficient tool to monitor the seismicity of mid-oceanic ridges (e.g. [START_REF] Fox | Monitoring Pacific Ocean seismicity from an autonomous hydrophone array[END_REF][START_REF] Smith | Hydroacoustic monitoring of seismicity at the slow-spreading Mid-Atlantic Ridge[END_REF][START_REF] Simão | Regional seismicity of the Mid-Atlantic Ridge: observations from autonomous hydrophone arrays[END_REF][START_REF] Goslin | Spatiotemporal distribution of the seismicity along the Mid-Atlantic Ridge north of the Azores from hydroacoustic data: Insights into seismogenic processes in a ridge-hot spot context[END_REF][START_REF] Tsang-Hin-Sun | Seismicity and active accretion processes at the ultraslow-spreading Southwest and intermediate-spreading Southeast Indian ridges from hydroacoustic data[END_REF]. The seismic activity of spreading ridges consists of many low-magnitude earthquakes most of which are not detected by land-based seismic networks (e.g. [START_REF] Bergman | Earthquake swarms on the Mid-Atlantic Ridge: Products of magmatism or extensional tectonics?[END_REF][START_REF] Bohnenstiehl | Aftershock sequences in the mid-ocean ridge environment: an analysis using hydroacoustic data[END_REF].

However, these weak events are capable of producing low-frequency acoustic waves (5-40 Hz), called "Tertiary" or T -waves, of an appreciable amplitude [START_REF] Williams | Hydroacoustic events located at the intersection of the Atlantis (30 • N) and Kane (23 • 40 N) Transform Faults with the Mid-Atlantic Ridge[END_REF][START_REF] Okal | The generation of T waves by earthquakes[END_REF].

The Sound Fixing and Ranging (SOFAR) channel, acting as a waveguide, allows T -waves to propagate in the ocean over long distances with very little attenuation. For the most efficient detection of T -waves, an AUH is typically positioned at the depth of the SOFAR channel axis, where it continuously records acoustic pressure variation. The Laboratoire Geosciences Ocean (LGO) currently maintains two AUH networks, HYDROMOMAR and OHASISBIO, deployed in the Atlantic and Indian oceans, respectively [START_REF] Perrot | HYDROMOMAR, Hydroacoustic observatory of the Mid-Atlantic Ridge[END_REF][START_REF] Royer | OHASISBIO -Hydroacoustic observatory for the seismicity and biodiversity in the Indian Ocean[END_REF][START_REF] Tsang-Hin-Sun | Seismicity and active accretion processes at the ultraslow-spreading Southwest and intermediate-spreading Southeast Indian ridges from hydroacoustic data[END_REF][START_REF] Giusti | The August 2010 earthquake swarm at North FAMOUS-FAMOUS segments, Mid-Atlantic Ridge: geophysical evidence of dike intrusion[END_REF]. T -wave arrivals recorded by several AUHs are subsequently used to deduce the location and the origin time of their source events. For a precise event localization, it is crucial that the time series of all instruments are synchronized. At the beginning of an experiment, the synchronization is ensured by synchronizing each AUH clock (which, in our case, is based on a Temperature-Compensated Crystal Oscillator) with the GPS time shortly before the instrument deployment. Once deployed at depth, the instrument clock triggers the data sampling and keeps track of the timing of the acoustic records. The instrument clock is however subject to a drift; in addition, drift rates vary from one instrument to another. As the experiment progresses in time, the synchronization between instruments is lost. Typically, AUHs are deployed for long time periods varying from one to two years and the skew can sometimes reach a substantial value. As soon as possible after a recovery of an instrument, the skew between the AUH and GPS clocks is measured with a microsecond precision, which allows us to correct each time series and thus restore the synchronization of the hydrophone array. The correction is done by assuming the skew to be a linear function of time (or equivalently, that the drift rate is time-independent). This plausible assumption is indirectly confirmed by the success of event localizations in all our hydroacoustic experiments. The same approach for the clock-drift correction is routinely used in experiments involving data acquisition with networks of ocean bottom seismometers (e.g. [START_REF] Geissler | Focal mechanisms for sub-crustal earthquakes in the Gulf of Cadiz from a dense OBS deployment[END_REF]).

However, in some cases measuring the final skew is not possible. For example, in July 2016, all five AUHs of the HYDROMOMAR network stopped functioning before they were recovered (due to logistical issues, the instruments were deployed longer than expected which lead to a premature exhaustion of the batteries). Without a clock-drift correction, event localization from T -waves becomes too imprecise to be useful for a seismicity analysis. Several authors showed that the clock drift can be estimated from interstation cross-correlations of the ambient noise (e.g. [START_REF] Stehly | Traveltime measurements from noise correlation: stability and detection of instrumental time-shifts[END_REF][START_REF] Hannemann | Measuring of clock drift rates and static time offsets of ocean bottom stations by means of ambient noise[END_REF][START_REF] Gouédard | Correction of Ocean-Bottom Seismometer Instrumental Clock Errors Using Ambient Seismic Noise[END_REF][START_REF] Hable | Clock errors in land and ocean bottom seismograms: high-accuracy estimates from multiple-component noise crosscorrelations[END_REF]. In a different approach, we propose to estimate the clock drift by using time markers naturally provided by the Earth: teleseismic P -wave arrivals detected by AUHs (see Table 1). Acoustic signals, generated by teleseismic P -waves refracting at the ocean bottom-water interface, are commonly observed (e.g. [START_REF] Slack | P wave detection thresholds, Pn velocity estimates, and T wave location uncertainty from oceanic hydrophones[END_REF][START_REF] Dziak | 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[END_REF][START_REF] Simons | On the potential of recording earthquakes for global seismic tomography by low-cost autonomous instruments in the oceans[END_REF][START_REF] Sukhovich | Seismic monitoring in the oceans by autonomous floats[END_REF]. In general, source events for such signals can be easily identified with the help of land-based seismic catalogs such as the one provided by the National Earthquake Information Center (NEIC, see Data and Resources Section), which we have used in this study. Using a ray propagation code, one can calculate the travel time between the hypocenter and the foot of the AUH mooring.

By adding the time that an acoustic wave takes to reach the AUH, the P -wave arrival time according to the GPS clock can be predicted. Finally, the drift rate of the AUH clock can be estimated by comparing expected and observed arrival times of many P -waves, which, ideally, are uniformly distributed over the entire duration of the experiment. Teleseismic P -waves (i.e. P -waves produced by distant events located at ∆ > 50 • away from the AUH), as compared to Clock drift from P-wave arrivals 5 seismic waves produced by nearby events, are better suited for such analysis since they traverse larger distances inside the globe, and the effect of inhomogeneities, responsible for deviation of actual travel times from calculated ones, is averaged out and becomes less important. Under the assumption of a time-independent drift rate, the difference between expected P -wave arrival times and those observed in the AUH record before correction should be a linear function of time.

To validate the method, we applied it to data from two AUHs of the HYDROMOMAR network, namely M2 and M7 (Figure 1), with known (measured) final skews. Both instruments operated for two years and were recovered in June 2018. The M7 hydrophone had an unusually high clock drift with a final skew of about 28 s, while that of the M2 hydrophone was approximately 7 times smaller (about 3.8 s). Thus, the analysis of the M2 data also tests the sensitivity limit of the method. Uncertainties in the predicted P -wave arrival times are estimated by a direct numerical approach and take into account uncertainties in the hypocenter depths, epicenter locations and origin times of the source events as well as the uncertainties in the Earth seismic-velocity model (used by the ray-tracing code). For both AUHs, we find that the actual drift rate (deduced from the measured skew) is located within the 95% confidence interval of the estimated drift rate. We also investigate the method performance as a function of the dataset size by considering subsets of P -wave arrivals to represent a possible scatter of estimated values. Our results demonstrate that, provided that a sufficient number of teleseismic P -wave arrivals are available, the method can reliably discern clock drifts of as small as 2 s per year.

In what follows, we first present the mathematical details of the method and then discuss the results.

METHOD Fitting

The term "clock drift" designates the deviation of a given clock relative to some high-precision reference clock, whose own deviation from the absolute time can be regarded as negligibly small.

In our case, we consider the drift of the AUH clock relative the GPS clock during the data acquisition period. At the beginning of the experiment the two clocks are synchronized, and the time difference (or "skew") ∆t between the time readings of two clocks is thus zero. We indicate by t AUH and t GPS the time intervals elapsed since the synchronization instance according to the AUH and the GPS clocks, respectively. At any moment of the record ∆t = t AUH -t GPS . Assuming that the skew is a linear function of time, it can be expressed as :

∆t = γ t GPS (1)
where γ indicates a constant drift rate (note that at the moment of synchronization, t AUH = t GPS = 0 and ∆t = 0). γ can be obtained from the final skew ∆t f measured at the end of the experiment :

γ = ∆t f t f,GPS (2)
Since ∆t f t f,GPS , γ is frequently expressed in either parts per million (ppm) or ms • day -1 .

Clock drift from P-wave arrivals 7

One can rewrite (1) to describe the clock drift in terms of t AUH and t GPS :

γ = t AUH t GPS (3)
where

γ = γ + 1 (4)
In the absence of clock drift, γ = 1 and γ = 0. If the instrument's clock runs faster or slower than the GPS clock, one finds γ > 1 (γ > 0) or γ < 1 (γ < 0), respectively.

Suppose that a P -wave, produced by a known event, is observed in the AUH record at some moment t AUH = t obs . At the same time, using a global seismic-velocity model, one can predict the time t pred at which the P -wave should have appeared in the record according to the GPS clock :

t pred = t or + t earth + t water (5)
where t or is the GPS time of the event origin (here, the time elapsed between the synchronization instance and the event origin), t earth is the P -wave travel time from the event hypocenter to the point on the ocean bottom at which the AUH mooring anchor is located, and t water is the time for a refracted acoustic wave to propagate from the ocean bottom to the hydrophone (note that the design of our moorings prevents the direct transmission of P -waves or their propagation through the mooring line).

If the AUH clock drifts, t obs and t pred will not be the same, and based on equation (3), their ratio will yield the value of γ :

γ = t obs t pred (6)
and subsequently, that of γ . Since both t obs and t pred are subject to uncertainties, a single P -wave arrival is not enough to reliably estimate γ. One needs to consider a sufficiently large number of P -wave arrivals, ideally distributed over the entire duration of the hydrophone record, and perform an optimization procedure. We estimate γ using a least-squares fit, which, in general, minimizes the quantity :

χ 2 = N i=1 (y i -y i ) 2 σ 2 i ( 7
)
where y i are experimental data points of some variable y, y i are its predicted values, N is the number of data points, and σ i is an uncertainty of a data point y i [START_REF] Bevington | Data reduction and error analysis for the physical sciences[END_REF]. The predictions y i are obtained from some assumed functional dependence of y on an

Clock drift from P-wave arrivals 9 independent variable x and, in general, K fitting parameters α k :

y i = f (α k , x i ) (k = 1, . . . , K)
To estimate the "goodness-of-fit", we use the quantity χ 2 ν :

χ 2 ν = χ 2 ν = χ 2 N -K (8)
classically known as "reduced chi square" with ν being the number of degrees of freedom [START_REF] Bevington | Data reduction and error analysis for the physical sciences[END_REF].

Our goal is to modify the time axis of the AUH clock so that all observed arrival times t i,obs agree with predicted arrival times t i,pred in the least-squares sense. Using equation ( 6), the functional dependence of t pred on t obs is :

f (γ, t obs ) = t obs γ (9)
This dependence is that of a straight line with a slope b = 1/γ and passing through the origin (remember that at the beginning of the experiment the instrument clock is synchronized with the GPS clock and thus, necessarily, t GPS = 0 when t AUH = 0). It is important to note that predicted arrival times t i,pred are based on a global seismic-velocity model which does not necessarily match identically the actual seismic-velocity profile at the hydrophone site. To take into account the site effect and to estimate correctly the slope b, a common time correction δt should be applied to all predicted arrival times : t i,pred → t i,pred + δt (this correction can be either positive or negative, depending on the differences between the local and global velocity profiles).

The value of δt is a priori unknown but can be estimated simultaneously with the slope b during the same fitting procedure from the intercept of the best-fitting line. Identifying f (γ, t obs ) as y i and t i,pred + δt as y i , equation ( 7) becomes :

χ 2 = N i=1 b t i,obs + a -t i,pred 2 σ 2 i ( 10
)
where a = -δt and N is the number of observed P -wave arrivals. Equation ( 10) represents a linear least-squares fit with two fitting parameters. It should be noted that the overall offset a can be a sum of several offsets, as will be shown in the Results section. Both γ and γ can be readily obtained from b (equation (4)).

The travel time t water was estimated from the bottom-hydrophone distance and the vertical sound-speed profile. In this work, we systematically took the length of the mooring line d l as the bottom-hydrophone distance, effectively assuming that the AUH was suspended in the water column directly above the mooring anchor. The validity of this assumption is discussed at the end of this section. Taking into account that both hydrophones were deployed in flat areas (within a radius of at least 5 km around the target position, see Figure 1) and the important acoustic impedance contrast between the sea-water and the oceanic crust, the seismic waves are expected to refract into the water column at angles very close to normal. Assuming also that the sound speed varies only with depth and not laterally, the acoustic waves are also expected to propagate

Clock drift from P-wave arrivals 11 over nearly vertical paths. The precise calculation of t water in case of a depth-dependent sound speed involves an integration over the sound wave propagation path. We found however that due to the short propagation distance, with a good degree of precision, one can use a simplified expression t water = d l cmean , where c mean is taken as the average of the two extreme values of the sound-speed profile (values at the ocean bottom and the SOFAR channel axis, where the AUH is positioned). By the same reason, equally well can be neglected a possible uncertainty in the acoustic travel time introduced by the uncertainty in the sound-speed profile. The time t water is thus treated as a constant with no associated uncertainty. To calculate t water , we used sound-speed profiles from the Generalized Digital Environmental Model produced by US Naval Oceanographic Office [START_REF] Countryman | The Navy Ocean-Temperature Temporal-Variability model, in "Challenges of Our Changing Global Environment[END_REF].

Estimation of uncertainties

According to equation ( 10), the calculation of χ 2 requires a knowledge of the uncertainty σ i for each fitted point (in what follows, index i is omitted to simplify the notation). Both t obs and t pred are subject to its own error, which we discuss below.

Uncertainty in observed arrival times

The uncertainty in t obs is user-estimated and represents the precision with which an operator performing the time-picking is able to identify the onset of a P -wave. It is not a true statistical error but rather a numerical criterion of the quality of the time-pick. We thus call this uncertainty "Quality Criterion" (QC). Obviously, its value mostly depends on the complexity of the signal onset and the signal-to-noise ratio, but also on the operator's experience. For the most complicated P -wave arrivals, QC may amount to 2 s. When the signal onset is clear and the time-pick is evident, we ascribed to QC a value of 0.05 s (equivalent to approximately 12 sampling intervals). Two examples of seismic phases with low and high QC are presented in Figure 2. As explained later, the QC was used to make a decision whether a given arrival should or should not be included in the fit. In the fit itself, no uncertainties were ascribed to t obs .

Uncertainty in predicted arrival times

The predicted arrival time t pred is a sum of the event origin time, the P -wave travel time in the Earth and the acoustic-wave travel time in the water column (equation ( 5)). The uncertainty σ th in t pred is thus a combination of the uncertainties in the event origin time σ or and the travel time σ prop . Its calculation is presented after the discussion of σ or and σ prop .

Assuming an event with a hypocenter located at depth d 0 , and epicenter's latitude and longitude θ 0 and φ 0 , respectively, several factors contribute to the σ prop :

-the uncertainty in the focal depth σ d , which characterizes the spread around the point d 0 of the possible depths of the hypocenter along the radial line passing through the point (θ 0 , φ 0 ) ; -the uncertainty in the hypocenter position σ p , which, at any given hypocenter depth, characterizes the spread around the point (θ 0 , φ 0 ) of possible hypocenter positions on the surface of a sphere whose radius is equal to the Earth radius minus the hypocenter depth ;

-the uncertainty in the Earth seismic-velocity model σ vm .

Due to the complexity of travel-time calculations, an analytical derivation of an error propagation formula to calculate σ prop directly from the above uncertainties is a complicated task.

Instead, we resorted to a simpler numerical approach. We first generated a collection of points Π, randomly distributed around the point (d 0 , θ 0 , φ 0 ), each point representing a possible hypocenter position. Second, we generated an ensemble of velocity models randomly distributed around the reference model (ak135 in our case) to reflect the uncertainty in our knowledge of the seismic Clock drift from P-wave arrivals 13 velocities in the Earth. For each pair "point in Π/velocity model", we then calculated a travel time to the AUH. This resulted in an ensemble of travel times t tr , from which the uncertainty σ prop was estimated.

In addition to the event's d 0 , θ 0 and φ 0 , the NEIC seismic catalog generally reports two errors, named depthError and horizontalError (both measured in km) and defined as the largest projections of the three principal axes of the error ellipsoid on a vertical axis and a horizontal plane, respectively. We used the value of depthError as σ d to generate an ensemble of points D which are randomly distributed according to the normal law with the mean d 0 and the standard deviation σ d .

The uncertainty σ p was assumed to be depth-independent and equal to horizontalError.

For each point d j in D (j is a dummy index indicating a particular depth value), the haversine formula was used to make a transition from σ p to uncertainties in latitude σ lat,j and longitude σ lon,j . Next, for each d j , a random ensemble P j was generated, with θ and φ distributed around θ 0 and φ 0 according to the normal law with standard deviations σ lat,j and σ lon,j , respectively. Points P j , combined for all depths d j , form the collection of points Π. In a more rigorous approach to generate the collection Π, one might attempt taking into account the directivity of the error ellipse axes (not reported by the NEIC catalog). However, the simpler approach adopted in this paper is sufficient to demonstrate the effectiveness of the method.

As an example, Figure 3 illustrates the procedure of estimating the spread in angular distance ∆ between the M2 hydrophone and an event, which happened on September 8, 2017 at 04:49:19.18 UTC (Table 1). According to the NEIC catalog, the event's depth d 0 was 47.39 km and its epicenter, located at θ 0 = 15.02 Calculating the travel times from all points in Π to the AUH location would result in a collection of travel times distributed with some standard deviation σ prop which is a function of σ d and σ p . The remaining source of uncertainty contributing to the σ prop is due to the uncertainty σ vm in the velocity model used for calculation of travel times. To take into account σ vm , a simple and straightforward way is to calculate travel times between points in Π and the AUH location using many different velocity models. An ensemble of such velocity models should be representative of the uncertainty in the principal velocity model and its realization is discussed later. This will result in an even larger collection of t earth , whose standard deviation will be a function of all three uncertainties σ d , σ p and σ vm .

Travel time calculations were performed with the ObsPy Python toolbox (see Data and Resources Section). As a principal model, we used the ak135 model introduced by [START_REF] Kennett | Constraints on seismic velocities in the Earth from traveltimes[END_REF]. According to these authors, with the exception of the inner core and the D region, the P -wave speed is constrained within ±0.01 km/s with a reasonable level of confidence. Upon examination of Figure 9 in [START_REF] Kennett | Constraints on seismic velocities in the Earth from traveltimes[END_REF], it seems that one can assign a higher uncertainty of ±0.015 km/s to the speed of P -waves in both the inner core and the D region. We generated 265 velocity models uniformly distributed around ak135 within the specified limits.

Another requirement stipulated that all velocity gradients in any generated velocity model did not deviate by more than 35% of the corresponding gradients in ak135. S1 andS2 of Supplemental material, for some events σ vm can provide a significant contribution to the final error σ prop . Therefore, for each event we systematically applied all three uncertainties σ d , σ p and σ vm to estimate σ prop .

In a final step, the uncertainty σ prop was combined with the uncertainty in origin time σ or , the latter being retrieved from the catalog of the International Seismological Centre (ISC, see Data and Resources Section) for the selected NEIC solutions. For all events, we found that the resulting distributions of travel times follow closely the normal law (see Figure 5(e) showing the September 8, 2017 event as an example). We assumed the event origin time t or to be also normally distributed with the variance σ 2 or . The sum of two normal random variables of variances σ 2 1 and σ 2 2 also being a normal random variable of variance σ 2 equal to the sum of both variances (see for example [START_REF] Mood | Introduction to the Theory of Statistics, McGraw-Hill international editions: Statistics series[END_REF]), the error in an event origin time can be accounted for by

considering σ th = σ 2 prop + σ 2 or .
Tables 2 and3 report the values of QC and all uncertainties for each observed P -wave arrival.

Arrivals for which QC was larger than the corresponding σ th were rejected and not used in the fit. For the remaining arrivals, σ th was used as σ in equation ( 10).

In this work, the location of the AUH corresponded to that of the mooring anchor whose location was known within a few meters. The bottom-hydrophone distance was taken to be equal to the length of the mooring line. It may happen that, in presence of deep oceanic currents, the mooring line will slant in the current direction and thus both AUH elevation and position will be different from their assumed values. A quantitative estimation of vertical and horizontal offsets introduced by currents is a complicated task which should be performed separately for each hydrophone location since in general the currents vary from one site to another. Moreover, in a given area, the offsets most likely vary with time as the buoy carrying the hydrophone moves around some equilibrium position. In a similar deployment of our AUH in the Indian Ocean, we were able to measure a horizontal offset of approximately 500 m for a 3560 m long line which is equivalent to a slant angle of about 8 • . The M2 and M7 hydrophones had shorter mooring lines, 2750 and 2450 m, respectively. Using simple geometric arguments, for the same slant angle, we estimated for the M2 and M7 hydrophones a horizontal offset of about 380 and 340 m and a vertical offset of about 27 and 24 m, respectively. The mean speed of sound c mean being about 1508 m • s -1 at both sites, the error in t water induced by neglecting the vertical offset would be less than 0.02 s, one order of magnitude smaller than the smallest estimated uncertainty σ th .

Similarly, the horizontal error in the location of the AUH (< 400 m) should not affect significantly the value of σ th compared to the combined effect of the uncertainties σ d and σ p . Therefore, neglecting possible depth and position uncertainties in the hydrophone position seemed to be justified.

Clock drift from P-wave arrivals 17

RESULTS

Using the NEIC seismic catalog, we first identified events of magnitude 6 or higher that we expected to be strong enough to produce acoustic signals detectable by the HYDROMOMAR network. Next, the origin times of the events, listed by the catalog in UTC time, were corrected for the leap second offset to be expressed in GPS time. Using the ObsPy toolbox, arrival times of all possible seismic phases were predicted and then matched with the acoustic signals (when observed). For the best clock-drift estimation, only signals with sufficiently high signal-to-noise ratio were retained. It is worth noting that not all of the selected events produced useful acoustic signals in the time series of both instruments (possibly, due to site effects and/or different ambient noise levels at the hydrophone locations on the day of detection). Jointly, the M2 and M7 hydrophones recorded 27 events, listed in Table 1 and shown in Figure 6.

Tables 2 and3 provide more details on the detected events, including the seismic phases expected to arrive at the hydrophone locations. Seismic phase names follow the current ObsPy convention according to which PKPab and PKPbc are both labeled as PKP while PKPdf is labeled as PKIKP. In case of triplications, it was not obvious to establish an unambiguous correspondence between an observed acoustic signal and its corresponding seismic phase for the phases predicted to arrive close in time. To avoid such ambiguity, we preferred to combine the seismic phases whose travel times were separated by less than 0.5 s, to produce a joint collection of travel times. For example, for the July 29, 2016 event, the combined phases were PKIKP and PKiKP, while for the August 31, 2016 and September 24, 2016 events, two PKP phases were combined (Table 2). For these events, the arrival time t pred was taken as the average of the arrival times predicted with the ak135 model for each of the combined phases.

In the ak135 model, the Earth is approximated by a solid sphere with an uppermost 20-kmthick layer in which the P -wave velocity is equal to 5.8 km • s -1 . In all travel-time calculations, ray paths originated at the earthquake hypocenter and ended at the surface of the sphere, which in our case corresponds to the ocean surface. Thus, all travel times calculated with ak135 for the P -waves arriving at the hydrophone site will exceed the actual travel times by the amount approximately equal to the ratio of the water depth at the hydrophone site to the P -wave velocity (5.8 km • s -1 ). The exact value of this offset obviously depends on the ray path and may vary slightly from one event to another, but to a good degree of precision can be considered the same for all arrivals. The presence of this time offset, as well as any eventual offset due to a site effect, is automatically taken into account by the parameter a in fit (10).

For the M7 hydrophone, the measured clock drift γ m was found to be of 0.437 ppm or 37.72 ms • day -1 (the corresponding final skew ∆t f , rounded to two decimal places, is 28.10 s;

note that its actual value is known with a microsecond accuracy). Figure 7(a) presents the predicted t pred versus observed t obs arrival times. At the scale of the entire experiment, the effect of the observed clock drift is not discernible, as illustrated by the dashed line of slope 1 representing the case of no clock-drift (i.e. t pred = t obs ). To make the presence of the clock drift evident, Figure 7(b) displays the differences (t obs -t pred ) versus the event number. Since γ m is positive, one expects that for each event t obs > t pred and their difference should be increasing linearly as the experiment progresses in time (equation ( 1)). This is indeed the case as can be seen in Figure 7(b). By fitting the data shown in Figure 7(a) (with the exception of two excluded points not consistent with the QC criterion), the estimated drift rate γ e was found to be 0.467 ppm (fit M7-1). The estimated γ e and observed γ m are in good agreement, the latter falling within 95% confidence interval of the former. A reasonably low value of χ 2 ν of 1.49 is consistent with the hypothesis of a constant drift rate. The value of the parameter a is 0.52 s and agrees well with a quick estimate of 0.59 s, calculated as a ratio of the water depth at the M7 site (3.4 km) to the P -wave velocity (5.8 km/s) assumed by ak135 in the corresponding layer.

To test the sensitivity limit of the method, we performed a similar analysis on the data Clock drift from P-wave arrivals 19 from the M2 hydrophone, for which the skew measured at the end of the experiment was almost 10 times smaller than that of the M7 hydrophone. More precisely, γ m was found to be 0.059 ppm or 5.12 ms • day -1 (the corresponding ∆t f , rounded to two decimal places, is 3.82 s).

Figures 8(a) and 8(b) present t pred versus t obs arrival times and their differences, respectively.

Two points inconsistent with the QC criterion (Table 2), were not used in the first fit indicated as fit M2-1. The fit yields γ e = 0.083 ppm and a 95% confidence interval bracketing γ m . The consistency between a linear model and the observed data is again supported by a χ 2 ν of 0.41.

The parameter a is equal to 0.57 s, again, consistent with a quick estimate of 0.62 s (the M2 site depth is 3.6 km).

The agreement between estimated and measured clock drifts can be improved even further if one takes into account that the ak135 model is developed from arrivals of seismic waves generated by earthquakes mostly located in the oceans and recorded by land-based seismological stations [START_REF] Kennett | Constraints on seismic velocities in the Earth from traveltimes[END_REF]. One might expect that in our case the ak135 would best describe P -waves from earthquakes originating in the continental crust and detected by a receiver (AUH) located in the ocean, since this situation is somehow similar to that of [START_REF] Kennett | Constraints on seismic velocities in the Earth from traveltimes[END_REF] except that receivers (at sea) and sources (on land) are switched. There is an equivalent global seismic-velocity model called ak135-f whose version for an oceanic lithosphere (see Data and Resources Section) mainly differs from the ak135 model only in the range between 0 and 35 km (water layer, sediment layer, oceanic crust and upper mantle instead of a 35-km-thick continental crust in ak135). To make sure that ak135 is applicable to all events used in the fit, we excluded, in fits M2-2 and M7-2, all oceanic events (i.e. events whose epicenters are in ocean)

shallower than 35 km (Table 1), since they certainly did not occur in a continental crust. With these events excluded, we find γ e of 0.434 ppm (fit M7-2, Figure 7(c)) and 0.059 ppm (fit M2-2, Figure 8(c)) for the M7 and M2 hydrophones, respectively. The values of the offset a, -0.93 s and 0.33 s for the M7 and M2 hydrophones, respectively, can be used to approximately quantify the site effects by subtracting from the estimated offset a the offset due to the ak135 model (which amounts, as explained above, to 0.59 s and 0.62 s for the M7 and M2 hydrophones, respectively).

One obtains the site contribution of -1.52 s for the M7 hydrophone and -0.29 s for the M2 hydrophone. This means that at both hydrophone locations, the P -waves arrive systematically later as predicted by the ak135 model (remember that in (10), a = -δt). This explanation seems to be quite reasonable considering that both hydrophones are located in the area influenced by a low-velocity anomaly due to the presence of the Azores hotspot [START_REF] Yang | Upper mantle structure beneath the Azores hotspot from finite-frequency seismic tomography[END_REF][START_REF] Silveira | Azores hotspot signature in the upper mantle[END_REF], with the M7 hydrophone actually positioned in the presumed hotspot location [START_REF] Silveira | Azores hotspot signature in the upper mantle[END_REF]). The amount of slowing down of P -waves is consistent with the results of Silveira et al. ( 2010) who found an evidence for the existence of a low P -wave velocity profile in the first 150 km of the upper mantle beneath the Azores. The site effect is less obvious in fits M7-1 and M2-1, possibly due to the presence of the arrivals of the P -waves originating in the oceanic crust. Table 4 summarizes the results of all four fits. More information on each fit can be found in Supplemental Material.

Figure 9 is a different way to present the results of the fits. For both hydrophones, it compares the observed time differences (t obs -t pred ) with those expected from equation (1) (after having taken into account the time offset a) at moments t pred based on measured γ m and estimated γ e drift rates.

To check how well the method performs when the sample size is decreased, we investigated for each hydrophone the spread of γ red values obtained from fits performed on reduced sets containing approximately 80% and 60% of the arrivals used in fits M2-2 and M7-2. Subsets were formed by choosing all possible combinations of n arrivals (irrespective of their order within a given combination) out of a larger collection of N arrivals (the total number of combinations C is given by N ! (N -n)! n!

). For the M2 hydrophone this resulted in selecting n = 12 and n = 9 arrivals out of N = 15 arrivals and for the M7 hydrophone in selecting n = 10 and n = 8 arrivals out of N = 13 arrivals. Figure 10 presents the histograms of obtained γ red values while

Clock drift from P-wave arrivals 21 the numerical results are presented in Table 5. In all cases, the mean value γ mean remains close to γ e found by fitting all arrivals. To measure the scatter of γ red values, for each collection shown in Figure 10 we calculated the ratio of the standard deviation σ n to γ mean . This quantity can be considered as an analog of the coefficient of variation σ/µ employed in statistics to characterize the dispersion of a probability distribution. Only in one case the ratio exceeds 20%, implying a reasonably narrow spreading of the estimated values γ red around the measured value γ m (this spread decreases when number of arrivals used in the fit increases). Finally, in each case we considered an interval containing 95% of γ red values. This interval can be viewed as an analog of a 95% confidence interval of a true random variable. As can be seen from Figure 10, for both hydrophones the measured value γ m always falls within the corresponding 95% interval.

CONCLUSION

We demonstrated that the drift rate of an AUH internal clock can be successfully estimated using arrivals of teleseismic P -waves. We employed a ray-tracing code in combination with uncertainties in the location, depth and origin time of the events as well as uncertainties in the

Earth seismic-velocity model, to place precision limits on estimated drift rates. By applying the method to data from two hydrophones with known clock drifts, we found that the drift rate resulting in a skew as small as 4 s per two years (about 5.5 ms • day -1 ) can be estimated with a good precision. Unlike techniques involving cross-correlations of the ambient noise (e.g. [START_REF] Gouédard | Correction of Ocean-Bottom Seismometer Instrumental Clock Errors Using Ambient Seismic Noise[END_REF][START_REF] Hable | Clock errors in land and ocean bottom seismograms: high-accuracy estimates from multiple-component noise crosscorrelations[END_REF]), our method cannot discern time variations in the drift rate, which we assumed to be constant over the entire duration of the experiment. However, the proposed approach is significantly less computationally intensive and simpler to implement. At the same time, the assumption of a time-invariant drift rate seems quite reasonable as evidenced by the low values of χ 2 ν obtained in all fits. The narrow spread of drift rates estimated from reduced subsets of arrival times demonstrates that the method is not very sensitive to the size of the dataset. In view of designing a practical tool for estimating the clock drift, selecting only events occurring in the continental crust and all other events with hypocenters located deeper than 35 km and applying the ak135 model provides a simple recipe and good results. A further improvement could be a calculation of travel times with a model which combines a continental structure on the source side, an oceanic structure on the receiver side and a common part for depths > 35 km. Additionally, the uncertainties in the observed arrival times t obs could be quantified, e.g. as recently proposed by [START_REF] Simon | Multiscale Estimation of Event Arrival Times and Their Uncertainties in Hydroacoustic Records from Autonomous Oceanic Floats[END_REF]. Our method is not limited to the field of hydroacoustics, and can be applied to any autonomous recording instruments, such as ocean bottom seismometers.

Clock drift from P-wave arrivals 23

Data and resources

The AUH data used in this study was acquired by and belongs to the LGO. The NEIC seismic catalog (https://earthquake.usgs.gov/earthquakes/search/, last accessed February 2020) was used to identify the source events of the observed P -wave arrivals. The full list of the networks contributing to the NEIC catalog can be found at https://earthquake. For each event, Table reports its time, depth, angular distance to the hydrophone, generated seismic phases and their travel times (as predicted by the ray-tracing code), horizontal, depth and origin time uncertainties, estimated uncertainty σprop in the travel time, estimated uncertainty σ th in the arrival time, the time-pick quality criterion QC and the uncertainty σ, used in equation ( 10). For the arrivals satisfying the condition QC < σ th , σ is equal to σ th . Otherwise, the arrival is rejected (not used in the fit) and its uncertainty σ is omitted. For each event, Table reports its time, depth, angular distance to the hydrophone, generated seismic phases and their travel times (as predicted by the ray-tracing code), horizontal, depth and origin time uncertainties, estimated uncertainty σprop in the travel time, estimated uncertainty σ th in the arrival time, the time-pick quality criterion QC and the uncertainty σ, used in equation ( 10).

For the arrivals satisfying the condition QC < σ th , σ is equal to σ th . Otherwise, the arrival is rejected (not used in the fit) and its uncertainty σ is omitted.

Clock drift from P-wave arrivals 33 The horizontal axes are identical column-wise. Clock drift from P-wave arrivals 39 4).

Clock drift from P-wave arrivals 47

  • N and φ 0 = 93.90 • W, was separated by the angular distance ∆ 0 of 57.35 • from the hydrophone location. The catalog reports the values of 6.2 and 3.7 km for the horizontalError and depthError, respectively. A histogram of the corresponding ensemble D is shown in Figure 3(a). Figures 3(b) and 3(c) show histograms of latitudes and longitudes of points in Π (all depths combined), respectively. Finally, Figure 3(d) displays a histogram of angular distances ∆ separating each point in Π and the mooring anchor of the M2 hydrophone. For each event, ensemble D contained 100 points, while each ensemble P j contained 900 points. All depths combined, each event's hypocenter was represented by a collection of 90000 points.

  Figure 5(c)). The histograms of the resulting collections of standard deviations are shown in Figures 5(b),(d) and (f). The mean of each collection is taken as the final value of σ prop while the number of significant decimal places in σ prop is estimated from the standard deviation of each collection. Comparison of Figures 5(b) and (d) shows that for this event the effect of the uncertainty in the velocity models on the width of the travel time distribution is significantly smaller than the effect of the uncertainty in the hypocenter location. Nevertheless, taking into account both uncertainties results in a slight increase of σ prop (compare the values displayed in Figures 5(b) and (f)). Moreover, as shown in TablesS1 and S2of Supplemental material, for

  (color online) Bathymetric maps showing the deployment areas of (a) the M2 and (b) the M7 hydrophones. Bathymetric data are extracted from the ETOPO1 global relief model (Amante and Eakins 2009). Isobaths are separated by 100 m. 2 Examples of P -waves detected by the M7 hydrophone with (a) high and (b) low quality criterion (QC) (see Table 3). Top row panels show the entire acoustic signals while the bottom row panels present an arbitrary zoom on the P -wave onset. The time axes correspond to the AUH time. The time origins were shifted with respect to the synchronization instant to make the reading of time axes easier. (a)-(b) Vertical lines indicate the predicted arrival times t pred . Although measured in the GPS time, t pred are indicated to make the presence of the clock drift visually evident. In its absence, t pred would coincide with (or be very close to) the observed onsets. (c)-(d) Vertical bars indicate the user time-picks while the width of the gray vertical bands shows the assigned QC value. The sampling frequency is 240 Hz. 3 Histograms illustrating the spread in (a) depths, (b) latitudes, (c) longitudes and (d) angular distances of all generated hypocenters used to estimate the uncertainty in the arrival time of the P phase produced by the September 8, 2017 event at 04:49:19.18 UTC. According to the NEIC catalog, the event depth d 0 was 47.39 km, the latitude θ 0 and the longitude φ 0 of its epicenter were 15.02 • N and 93.90 • W, respectively, corresponding to an angular distance ∆ 0 of 57.35 • to the M2 hydrophone location. 4 (color online) Generated seismic velocity models compared to the principal model ak135. (a) The ak135 model. (b)-(d) The three panels compare the ak135 model (thicker line with filled circles) and the generated models (thinner lines with empty circles) by zooming at specific depths, indicated in (a) by horizontal dashed lines. Left panels : Histograms of travel times for the P phase propagating between the hypocenter of the September 8, 2017 event and the M2 hydrophone. The travel times are calculated under different conditions : (a) using the ak135 model and the hypocenters with depths and angular distances shown in Figure 3; (c) using the velocity models shown in Figure 4 and the hypocentre position listed in the NEIC catalog (d 0 = 47.39 km and ∆ 0 = 57.35 • ); (e) using all possible pairs "hypocenter/velocity model" formed from the collections of hypocenters and velocity models employed in the first two cases. Vertical dotted lines indicate the travel time predicted by the ak135 model for the hypocenter at (d 0 , ∆ 0 ). Solid line in (e) is the normal probability density function found by fitting the data. Right panels : Histograms of the standard deviations σ prop obtained (see text) from the ensemble of the travel times shown in the corresponding left panel. All histograms are normalized to have an area of 1.6 Global map of the teleseismic events (circles,

(

  solid line) drift rates for the events used in (a) fit M2-2 and (b) fit M7-2. The error bars represent uncertainties σ on t pred and the gray zone is bordered by limiting values of the 95% confidence interval on γ e (as reported in Table4).10 Histograms of the γ red obtained from the fits based on subsets of the arrivals used in fits M2-2 and M7-2. Each subset contains (a)-(b) 80% of total number of arrivals (12 and 8 in case of the M2 and M7 hydrophones, respectively) and (c)-(d) 60% (11 and 7 in case of the M2 and M7 hydrophones, respectively). In each panel, the vertical line indicates the value of the measured γ m while the darker gray color outlines an interval occupied by 95% of γ red values. All histograms are normalized to have an area of 1.

Figure 1 .

 1 Figure 1. (color online) Bathymetric maps showing the deployment areas of (a) the M2 and (b) the M7 hydrophones. Bathymetric data are extracted from the ETOPO1 global relief model (Amante and Eakins 2009). Isobaths are separated by 100 m.

Figure 2 .

 2 Figure 2. Examples of P -waves detected by the M7 hydrophone with (a) high and (b) low quality criterion (QC) (see Table 3). Top row panels show the entire acoustic signals while the bottom row panels present an arbitrary zoom on the P -wave onset. The time axes correspond to the AUH time. The time origins were shifted with respect to the synchronization instant to make the reading of time axes easier. (a)-(b) Vertical lines indicate the predicted arrival times t pred . Although measured in the GPS time, t pred are indicated to make the presence of the clock drift visually evident. In its absence, t pred would coincide with (or be very close to) the observed onsets. (c)-(d) Vertical bars indicate the user time-picks while the width of the gray vertical bands shows the assigned QC value. The sampling frequency is 240 Hz.

Figure 3 .

 3 Figure 3. Histograms illustrating the spread in (a) depths, (b) latitudes, (c) longitudes and (d) angular distances of all generated hypocenters used to estimate the uncertainty in the arrival time of the P phase produced by the September 8, 2017 event at 04:49:19.18 UTC. According to the NEIC catalog, the event depth d 0 was 47.39 km, the latitude θ 0 and the longitude φ 0 of its epicenter were 15.02 • N and 93.90 • W, respectively, corresponding to an angular distance ∆ 0 of 57.35 • to the M2 hydrophone location.

Figure 4 .

 4 Figure 4. (color online) Generated seismic velocity models compared to the principal model ak135. (a) The ak135 model. (b)-(d) The three panels compare the ak135 model (thicker line with filled circles) and the generated models (thinner lines with empty circles) by zooming at specific depths, indicated in (a) by horizontal dashed lines.

Figure 9 .

 9 Figure 9. Comparison of the observed time differences (t obs -t pred ) with those expected from equation (1) at moments t pred in case of measured γ m (dashed line) and estimated γ e (solid line) drift rates for the events used in (a) fit M2-2 and (b) fit M7-2. The error bars represent uncertainties σ on t pred and the gray zone is bordered by limiting values of the 95% confidence interval on γ e (as reported inTable 4).

Figure 10 .

 10 Figure 10. Histograms of the γ red obtained from the fits based on subsets of the arrivals used in fits M2-2 and M7-2. Each subset contains (a)-(b) 80% of total number of arrivals (12 and 8 in case of the M2 and M7 hydrophones, respectively) and (c)-(d) 60% (11 and 7 in case of the M2 and M7 hydrophones, respectively). In each panel, the vertical line indicates the value of the measured γ m while the darker gray color outlines an interval occupied by 95% of γ red values. All histograms are normalized to have an area of 1. The horizontal axes are identical column-wise.

Table 1 :

 1 usgs.gov/data/comcat/contributor/ (last accessed February 2020). The ISC catalog (http://www.isc.ac.uk, last accessed February 2020;[START_REF] Bondár | Improved location procedures at the International Seismological Centre[END_REF] was used to obtain the uncertainties on the event origin times. Maps in Figures1 and 6were made using the Generic Mapping Tools version 4.5.18 of the paper. The authors also thank : Mickaël Beauverger, the chief manager of the AUH park at LGO, for his significant contribution to all AUH experiments; captains and crews of the research vessels for their help and seamanship during the HYDROMOMAR campaigns in 2016 and 2018; and last but not least, ObsPy and GMT developers for their time and effort devoted to the design and improvements of these useful open-source tools. Detected teleseismic events.

	Clock drift from P-wave arrivals 27

(www.soest.hawaii.edu/gmt;

[START_REF] Wessel | Generic mapping tools: Improved version released[END_REF]

. Theoretical travel times were computed using the package obspy.taup which is a part of a publicly available ObsPy tool box (https://docs.obspy.org/packages/ last accessed August 2019;

[START_REF] Beyreuther | ObsPy: A Python Toolbox for Seismology[END_REF]

. The oceanic version of the ak135-f model is available at http://rses.anu.edu.au/seismology/ak135/ak135f.html (last accessed September 2020). The fitting was performed using Curve Fitting Toolbox of MATLAB. The toolbox allows one to define a model and reports both values and confidence intervals of the fitting parameters. Supplemental material to this paper provides information on how, for each event, uncertainties in the seismic-velocity model and in hypocenter position contribute to the final uncertainty in travel time. It also includes extra details on the fits. improvement the time series of the M2 and/or M7 hydrophones. The instruments were deployed on June 17 and June 13, 2016 and recovered on July 2 and June 28, 2018, respectively. Column "Oceanic" identifies the events whose epicenters are oceanic and depths are shallower than 35 km.

Table 2 :

 2 Events detected by the M2 hydrophone (39.42 • N, 34.11 • W, site depth : 3600 m).

	Event time (UTC)	Depth	∆ ( • )	Phase	Travel	σp	σ d	σor	σprop	σ th	QC	σ
		(km)			time (s)	(km)	(km)	(s)	(s)	(s)	(s)	(s)
	2016-07-29T21:18:24.740Z	196.00 122.34	PKIKP	1111.57	7.50	1.80	2.80	0.30	2.82	0.05	2.82
				PKiKP	1111.77
	2016-08-31T03:11:34.420Z	476.00 143.88	PKP	1119.28	8.10	1.90	1.27	0.35	1.32	0.05	1.32
				PKP	1119.32
				PKIKP	1120.34
				PKiKP	1124.00
	2016-09-24T21:28:41.700Z	596.40 143.51	PKP	1106.38	9.20	3.50	1.67	0.49	1.74	0.40	1.74
				PKP	1106.43
				PKIKP	1107.55
				PKiKP	1111.17
	2016-11-24T18:43:47.710Z	10.00	55.47	P	574.89	5.60	1.60	2.07	0.44	2.12	0.30	2.12
	2017-01-22T04:30:22.960Z	135.00 145.93	PKIKP	1162.68	7.30	1.80	2.48	0.27	2.49	0.50	2.49
				PKP	1163.60
				PKP	1164.30
				PKiKP	1166.86
	2017-02-21T14:09:04.320Z	595.98 64.70	P	580.46	9.60	3.00	1.83	0.66	1.94	0.20	1.94
	2017-04-03T17:40:18.560Z	29.00	82.92	P	742.01	7.50	1.80	2.25	0.45	2.30	0.30	2.30
	2017-04-24T21:38:30.820Z	28.00	80.22	P	727.88	4.90	1.80	1.48	0.38	1.53	1.00	1.53
	2017-06-02T22:24:47.440Z	5.00	84.46	P	753.68	6.70	1.70	1.84	0.43	1.89	0.70	1.89
	2017-06-14T07:29:04.390Z	93.00	55.97	P	568.62	2.80	2.20	4.34	0.28	4.35	0.40	4.35
	2017-06-22T12:31:03.490Z	38.12	55.91	P	574.00	6.80	3.50	2.96	0.61	3.02	0.20	3.02
	2017-07-17T23:34:13.740Z	10.00	84.45	P	752.80	7.00	1.80	1.37	0.45	1.44	2.00	-
	2017-07-20T22:31:11.260Z	7.00	47.60	P	516.42	4.30	1.70	1.59	0.41	1.64	2.00	-
	2017-09-08T04:49:19.180Z	47.39	57.35	P	583.26	6.20	3.70	2.47	0.55	2.53	0.40	2.53
	2017-09-19T18:14:38.090Z	48.00	58.82	P	593.46	4.50	1.80	2.32	0.35	2.35	0.05	2.35
	2017-11-12T18:18:17.180Z	19.00	61.96	P	618.55	4.90	1.70	2.62	0.41	2.65	0.20	2.65
	2018-01-10T02:51:33.290Z	19.00	47.95	P	517.26	5.80	1.70	1.65	0.49	1.72	0.30	1.72
	2018-01-14T09:18:45.540Z	39.00	66.71	P	646.91	7.00	1.80	2.43	0.45	2.47	0.30	2.47
	2018-01-23T09:31:40.890Z	14.06	70.19	P	672.23	6.00	2.80	1.57	0.54	1.66	0.05	1.66
	2018-02-16T23:39:39.280Z	22.00	59.69	P	602.66	4.30	1.80	2.49	0.39	2.52	0.30	2.52
	2018-03-26T09:51:00.430Z	40.00	145.87	PKIKP	1174.28	6.50	0.80	2.39	0.20	2.40	0.50	2.40
				PKP	1174.92
				PKP	1175.40

Table 3 :

 3 Events detected by the M7 hydrophone (33.45 • N, 32.40 • W, site depth : 3400 m).

		Event time (UTC)	Depth	∆ ( • )	Phase	Travel	σp	σ d	σor	σprop	σ th	QC	σ (s)
			(km)			time (s)	(km)	(km)	(s)	(s)	(s)	(s)
	1	2016-07-29T21:18:24.740Z	196.00 128.26	PKIKP	1122.90	7.5	1.8	2.80	0.30	2.82	0.60	2.82
					PKiKP	1123.54
	2	2016-08-19T07:32:22.710Z	10.00	88.38	P	772.12	8.4	1.7	1.59	0.48	1.66	0.40	1.66
	3	2016-08-29T04:29:57.860Z	10.00	36.03	P	421.37	8.4	1.7	1.76	0.70	1.90	5.00	-
	4	2016-11-25T14:24:30.710Z	17.00	80.62	P	731.70	5.9	1.7	1.78	0.42	1.83	0.80	1.83
	5	2016-12-25T14:22:27.010Z	38.00	85.39	P	753.33	5.4	1.9	1.59	0.36	1.63	0.05	1.63
	6	2017-02-21T14:09:04.320Z	595.98 60.45	P	553.13	9.6	3.0	1.83	0.67	1.95	0.05	1.95
	7	2017-04-03T17:40:18.560Z	29.00	78.23	P	716.88	7.5	1.8	2.25	0.47	2.30	0.50	2.30
	8	2017-04-24T21:38:30.820Z	28.00	75.94	P	704.15	4.9	1.8	1.48	0.40	1.53	0.50	1.53
	9	2017-06-14T07:29:04.390Z	93.00	56.69	P	573.72	2.8	2.2	4.34	0.32	4.35	0.20	4.35
	10	2017-07-20T22:31:11.260Z	7.00	48.31	P	521.95	4.3	1.7	1.59	0.42	1.64	0.05	1.64
	11	2017-09-08T04:49:19.180Z	47.39	58.21	P	589.27	6.2	3.7	2.47	0.58	2.54	1.50	2.54
	12	2017-09-19T18:14:38.090Z	48.00	60.33	P	603.86	4.5	1.8	2.32	0.37	2.35	0.10	2.35
	13	2017-11-12T18:18:17.180Z	19.00	63.18	P	626.70	4.9	1.7	2.62	0.41	2.65	0.40	2.65
	14	2018-01-10T02:51:33.290Z	19.00	48.35	P	520.31	5.8	1.7	1.65	0.49	1.72	3.00	-
	15	2018-01-14T09:18:45.540Z	39.00	63.44	P	625.71	7.0	1.8	2.43	0.47	2.48	0.05	2.48
	16	2018-01-23T09:31:40.890Z	14.06	75.97	P	706.44	6.0	2.8	1.57	0.58	1.67	0.05	1.67
	17	2018-02-16T23:39:39.280Z	22.00	60.95	P	611.29	4.3	1.8	2.49	0.39	2.52	0.50	2.52
	18	2018-04-02T13:40:34.840Z	559.00 61.22	P	561.06	6.7	1.9	2.01	0.47	2.06	0.15	2.06

Table 4 :

 4 Results of the fits.For each fit, Table reports the estimated drift rate γ e , its 95% confidence interval (CI), estimated value of the time offset a and the goodness-of-fit χ 2 ν .

	AUH fit id	γ e	γ e 95% CI interval a (s) χ 2 ν
	M2-1	0.083	[0.048 0.118]	0.57 0.41
	M2
	M2-2	0.059	[0.021 0.096]	0.33 0.31
	M7-1	0.467	[0.382 0.551]	0.52 1.49
	M7
	M7-2	0.434	[0.368 0.501]	-0.93 0.54

Table 5 :

 5 Estimation of the drift-rate spread.For each hydrophone, subsets containing 80% and 60% of the recorded arrivals were considered. Tablereportsthe number of available arrivals N , the number of arrivals n composing a subset and the number of fitted subsets C given by the number of combinations formed by choosing (irrespective of the order within a given combination) n items from N items. As explained in the text, each fit yields a different γ red . From the resulting collection of γ red (Figure10), a mean value γ mean , an interval occupied by 95% of γ red values and a ratio of the standard deviation σn to γ mean are estimated in each case and reported in Table.

	AUH N	n	C	γ mean 95% interval	σ n /γ mean (%)
		12 455 0.060 [0.045 0.074]	12.08
	M2 15
		9 5005 0.060 [0.035 0.087]	22.12
		10 286 0.435 [0.402 0.466]	3.86
	M7 13
		8 1287 0.435 [0.392 0.490]	6.04

Table 1

 1 Predicted t pred versus observed t obs arrival times. The dashed line indicates the case of no clock-drift (i.e. t pred = t obs ). The uncertainties in t pred (Table3) are not indicated as they are not visible at the scale of the experiment. (b) Observed time differences (t obs -t pred ). Filled circles indicate events labeled as oceanic in Table1. (c) Residuals (t obs /γ e -t pred ), obtained in the fit M7-2, compared with the corresponding uncertainties σ (Table3) represented by bars (each bar takes an interval from -σ to σ). The absolute value of most residuals is smaller than the corresponding σ, which is reflected by the small value of χ 2 ν . Arrival numbering corresponds to that of Table 3. Oceanic events shallower than 35 km and events not passing the quality criterion are excluded from the fit. Estimation of γ for the M2 hydrophone. (a) Predicted t pred versus observed t obs arrival times. The dashed line indicates the case of no clock-drift (i.e. t pred = t obs ). The uncertainties in t pred (Table 2) are not indicated as they are not visible at the scale of the experiment. (b) Observed time differences (t obs -t pred ). Filled circles indicate events labeled as oceanic in Table 1. (c) Residuals (t obs /γ e -t pred ), obtained in the fit M2-2, are compared with the corresponding uncertainties σ (Table 2) represented by bars (each bar takes an interval from -σ to σ). The absolute value of almost all residuals is smaller than the corresponding σ, which is reflected by the small value of χ 2 ν . Arrival numbering corresponds to that of Table 2. Oceanic events shallower than 35 km and events not passing the quality criterion are excluded from the fit. 9 Comparison of the observed time differences (t obs -t pred ) with those expected from equation (1) at moments t pred in case of measured γ m (dashed line) and estimated γ e

	) recorded by the M2 (black

Acknowledgments

Alexey Sukhovich received financial support from a CNRS "Chaire d'Excellence" award. The HYDROMOMAR experiments were financially supported by the French "Agence Nationale de la Recherche" ("Laboratoire d'Excellence" LabexMER program, ANR-10-LABX-19-01), the French government ("Investissements d'Avenir" grant), INSU-CNRS via the EMSO-Açores program and the French Oceanographic Fleet. The authors wish to thank two anonymous Reviewers and Editors for a constructive criticism and valuable suggestions, which undoubtedly led to an

LIST OF TABLES

Detected teleseismic events.

2 Events detected by the M2 hydrophone (39.42 • N, 34.11 • W, site depth : 3600 m).

3 Events detected by the M7 hydrophone (33.45 • N, 32.40 • W, site depth : 3400 m).

4

Results of the fits. The dashed line indicates the case of no clock-drift (i.e. t pred = t obs ). The uncertainties in t pred (Table 3)

are not indicated as they are not visible at the scale of the experiment. (b) Observed time differences (t obs -t pred ). Filled circles indicate events labeled as oceanic in Table 1. (c) Residuals (t obs /γ e -t pred ), obtained in the fit M7-2, compared with the corresponding uncertainties σ (Table 3) represented by bars (each bar takes an interval from -σ to σ). The absolute value of most residuals is smaller than the corresponding σ, which is reflected by the small value of χ 2 ν . Arrival numbering corresponds to that of Table 3. Oceanic events shallower than 35 km and events not passing the quality criterion are excluded from the fit.

Clock drift from P-wave arrivals 45 The dashed line indicates the case of no clock-drift (i.e. t pred = t obs ). The uncertainties in t pred (Table 2) are not indicated as they are not visible at the scale of the experiment. (b) Observed time differences (t obs -t pred ). Filled circles indicate events labeled as oceanic in Table 1. (c) Residuals (t obs /γ e -t pred ), obtained in the fit M2-2, are compared with the corresponding uncertainties σ (Table 2) represented by bars (each bar takes an interval from -σ to σ). The absolute value of almost all residuals is smaller than the corresponding σ, which is reflected by the small value of χ 2 ν . Arrival numbering corresponds to that of Table 2. Oceanic events shallower than 35 km and events not passing the quality criterion are excluded from the fit.