# Revisiting the OBS seafloor compliance signal removal with a stationarity and stacking-based approach: the BRUIT-FM toolbox

Simon Rebeyrol, Stephan Ker, Laurent Duval, Wayne C. Crawford

► To cite this version:

Simon Rebeyrol, Stephan Ker, Laurent Duval, Wayne C. Crawford. Revisiting the OBS seafloor compliance signal removal with a stationarity and stacking-based approach: the BRUIT-FM toolbox. Geophysical Journal International, 2024, 239 (1), pp.386-401. 10.1093/gji/ggae265. hal-04731639

HAL Id: hal-04731639

https://hal.science/hal-04731639v1

Submitted on 14 Oct 2024

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 Attribution 4.0 International License
# Revisiting the OBS seafloor compliance signal removal with a stationarity and stacking-based approach: the BRUIT-FM toolbox

Simon Rebeyrol<sup>1</sup>,<sup>1</sup> Stéphan Ker<sup>1</sup>,<sup>1</sup> Laurent Duval<sup>2,3</sup> and Wayne C. Crawford<sup>4</sup>

Accepted 2024 July 19. Received 2024 July 7; in original form 2023 August 8

## SUMMARY

This study focuses on improving the seafloor compliance noise removal method, which relies on estimates of the compliance transfer function frequency response (the deformation of the seafloor under long-period pressure waves). We first propose a new multiscale deviation analysis of broad-band ocean-bottom seismometer data to evaluate stationarity properties that are key to the subsequent analysis. We then propose a new approach to removing the compliance noise from the vertical channel data, by stacking daily estimated transfer function frequency responses over a period of time. We also propose an automated transient event detection and data selection method based on a cross-correlation criterion. As an example, we apply the method to data from the Cascadia Initiative (network 7D2011). We find that the spectral extent of long-period forcing waves varies significantly over time so that standard daily transfer function calculation techniques poorly estimate the transfer function frequency response at the lowest frequencies, resulting in poor denoising performance. The proposed method more accurately removes noise at these lower frequencies, especially where coherence is low, reducing the mean deviation of the signal in our test case by 27 per cent or more. We also show that our calculated transfer functions can be transferred across time periods. The method should allow better estimates of seafloor compliance and help to remove compliance noise at stations with low pressure-acceleration coherence. Our results can be reproduced using the BRUIT-FM Python toolbox, available at https://gitlab.ifremer.fr/anr-bruitfm/bruit-fm-toolbox.

**Key words:** Fourier analysis; Seismic noise; Surface waves and free oscillations.

## 1 INTRODUCTION

The vertical component of broad-band ocean-bottom seismometers (BBOBS) record a seafloor compliance signal at frequencies below approximately 100 mHz (Crawford et al. 1998). The compliance signal results from the vertical motion of the seabed in response to pressure variations of long-period ocean waves. These pressure variations decay quasi-exponentially with depth, with a decay constant that depends on the wavelength, which itself depends on the frequency (Webb et al. 1991; Aucan & Ardhuin 2013). Thus, the maximum frequency of the compliance signal decreases from 120 mHz at a water depth of 100 m to 20 mHz at a water depth of 4000 m (Bell et al. 2015). The compliance signal can be inverted to assess the sub-seabed structure as a function of depth and has been used to study melt in the oceanic crust and uppermost mantle, gas

hydrates and sub-basalt sediments (Crawford 1994, 2004; Willoughby et al. 2008). However, the compliance signal also masks long-period earthquake signals, and needs to be removed for better analysis of these signals (Crawford & Webb 2000).

The removal of compliance noise has been facilitated by public-domain code packages such as ATaCR (Janiszewski et al. 2019) and TisKit.<fnref n="1" /> The standard technique to remove compliance noise, first proposed by Webb & Crawford (1999) and Crawford & Webb (2000), is based on estimating the seafloor compliance transfer function frequency response (scTFFR) following the methodology of Bendat & Piersol (2011). This approach uses power spectral and cross-power spectral density (PSD and CPSD) estimates

<footnote n="1">¹Geo-Ocean UMR 6538 CNRS – Ifremer – UBO – UBS, 1625 route de Sainte-Anne CS 10070 - 29280 Plouzané, France.</footnote>

<footnote>E-mails: simon.rebeyrol@murena.io; stephan.ker@ifremer.fr</footnote>

<footnote n="2">²ESIEE-Paris, Université Gustave Eiffel, F-77454 Marne-la-Vallée, France</footnote>

<footnote n="3">³IFP Energies nouvelles, 1-4 avenue du Bois-Préau, Rueil-Malmaison F-92852, France</footnote>

<footnote n="4">⁴Géosciences marines, Institut de Physique du Globe de Paris, Paris, France</footnote>

<footnote n="1">¹https://github.com/WayneCrawford/tiskitpy</footnote>
(Welch 1967) of pressure and vertical component signals, which are considered as stationary random processes. Several key parameters and methods are used when calculating the compliance, including: the size of the time windows used to compute the PSDs, techniques to avoid transient signals, rejecting biased or noisy PSDs and defining a coherence threshold to determine the frequency range over which the removal is realized. The key parameters are defined in Table A1 and further explained in Section 2.

To 'assure' the signal stationarity assumed by the transfer function frequency response calculation, most studies limit the overall estimation window to one day and the length of individual time windows to between 17 min and 3 hr (1000–10 000 s). This methodology requires careful removal of transient events, which have a significant influence on signal stationarity. Another approach, proposed by McNamara & Buland (2004) and Yang et al. (2012) reduces the overall estimate window to 1 hr, sliced into 13 overlapping segments. While the resulting frequency resolution does not cover the compliance frequency band, this parametrization allows many more PSD estimations, reducing further the noise and allowing easy PSD-level rejection of transient events.

The values of key parameters, like the length of the time windows, are generally defined by users based on their experience, need and bias and vary significantly (Bell et al. 2015; Janiszewski et al. 2019, 2022; Crawford & Webb 2000). The choice of these parameters may strongly influence the calculated compliance transfer function and the efficiency of compliance noise removal, but no study of this influence exists.

We study here the robustness of the compliance transfer function estimation and the efficiency of compliance noise removal. After a quick overview of the theory of transfer functions (Section 2), we describe data selection and conditioning (Section 3), then propose a multiscale deviation analysis to determine the parameter values for which the stationarity assumption is valid (Section 4). Section 5 presents an improved compliance noise removal method, using the parameters calculated in the previous section and presenting the important automated steps of transient event detection and transfer function stacking. The results are shown in Section 6, followed by a discussion and conclusion.

# 2 REMINDER OF BACKGROUND THEORY AND APPLICATION

## 2.1 Seafloor compliance signal removal

For more than two decades, highly sensitive BBOBS (Webb 1998) have been used to record subtle seismological signals below 0.2 Hz, including microseism and infra-gravity (IG) waves. IG waves are generated by nonlinear interactions between ocean surface waves, mostly near coastlines, and propagate into the deep ocean with little attenuation (Webb et al. 1991). The seafloor deforms by up to 10 µm vertically under the pressure fluctuations generated by IG waves (Crawford et al. 1998), which is enough to be recorded by the BBOBS. Seafloor compliance is defined as the ratio of the differential seafloor vertical displacement to the seafloor pressure.

We summarize below the procedure to estimate and remove seafloor compliance from the data. To simplify comparison with other methods, we use data with little to no current-induced noise and we do not apply noise correction. The most commonly used current-noise reduction process is analogous to the compliance noise removal: our comparison is more direct if we focus only

on compliance. In the frequency range of the IG waves, the seafloor compliance relation is considered linear and characterized by its impulse response $h(t)$. This relation links the clean (compliance-free) vertical acceleration $z(t)$, the recorded acceleration $z_c(t)$ and the pressure $p(t)$ signals. An optimal estimate of the seafloor compliance frequency response $\hat{H}(f)$ (hereafter noted scTFFR) can be obtained following the procedure of Bendat & Piersol (2011) which states that:

$$ \hat{H}(f) = \frac{G_{z_c} p(f)}{G_p(f)}, \qquad (1) $$

where $G_{z_c}$ is the cross-spectrum for signals $z_c(t)$ and $p(t)$ and $G_p(f)$ is the individual auto-power spectrum of $p(t)$. In practice, one often improves this estimate by splitting the pressure and acceleration signals (of total length $N_T$ seconds) into K windows of length $N_w$, with an overlap of $N_o$ seconds. Each signal is apodized by a suitable taper window. These are the key parameters used in the BRUIT-FM (Bruits des Fonds Marins) toolbox and are summarized in Table A1. We then calculate scTFFR Welch's modified periodogram (Welch 1967), by averaging windowed cross and power spectra.

The length $N_w$ of each windowed segment is a critical compromise. A longer $N_w$ provides higher frequency resolution and lower Signal-to-Noise Ratio (SNR), whereas a shorter $N_w$ may be needed to ensure that the signal is stationary within each segment. There are also several practical condition mismatches to account for. First, $z(t)$ and $p(t)$ may not be always be fully uncorrelated, for example, a low-magnitude/frequency earthquake could be recorded by both pressure sensors and seismometers. A selection step is usually conducted to remove data contaminated by known earthquakes (e.g. Bell et al. 2015; An et al. 2020). Second, $z_c(t)$ and $p(t)$ may not carry enough information to estimate $\hat{H}(f)$ properly. A common procedure is to estimate $\hat{H}(f)$ in a reduced spectral range in which with coherence, defined as:

$$ C_{pz_c}(f) = \frac{|G_{pz_c}(f)|}{\sqrt{G_{pp}(f)G_{z_c z_c}(f)}} \in [0, 1] \qquad (2) $$

is relatively high. Appendix A2 describes the complete methodology used in the BRUIT-FM toolbox.

## 2.2 Signal stationarity

As stated above, the derivation of scTFFR strongly relies on stationarity assumptions. Second-order stationarity is a strong requirement that must be carefully verified, especially for year-long or multiyear data sets. To mitigate this, most of the past studies processed at most 1 d of data at a time, imposing a reduced window length and relatively low-frequency resolution (Crawford et al. 1991, 1998; Crawford 2004; Bell et al. 2015; Crawford & Webb 2000). For example, Crawford (2004) estimated the scTFFR using 29 hr of data and 1024-s data windows. More recently, An et al. (2020) used an 10 000-s data windows (and an unspecified overall time period), and Janiszewski et al. (2022) used 7200-s windows and one day of data, averaging the resulting PSDs over 25 d. These parameter values are summarized in Table 1.

To our knowledge, no stationarity checks were conducted prior to removing seafloor compliance noise in the above studies. A stationarity check assesses variations in the signal deviation over different length-scales. We use this assessment to determine the range of values of the key parameters $N_w$ and $N_T$ that respect the stationarity assumption.
<table><caption>**Table 1.** Key parameter values of the PSDs estimation used in past studies. $N_w$ is the window length and $N_T$ is the length of the signal to be analysed.</caption><thead><tr><th>Studies</th><th>N<sub>w</sub> (s)</th><th>N<sub>T</sub></th></tr></thead><tbody><tr><td>Webb & Crawford (1999)</td><td>1024</td><td>24 hr</td></tr><tr><td>Crawford (2004)</td><td>1024</td><td>29 hr</td></tr><tr><td>Yang et al. (2012)</td><td>100</td><td>1 hr</td></tr><tr><td>Bell et al. (2015)</td><td>2000</td><td>11 hr 7 min</td></tr><tr><td>Doran & Laske (2019)</td><td>3600</td><td>Not specified</td></tr><tr><td>An et al. (2020)</td><td>10 000</td><td>Not specified</td></tr><tr><td>Janiszewski et al. (2022)</td><td>7200</td><td>25 × 24 hr</td></tr></tbody></table>

<table><caption>**Table 2.** Selected station information. We estimated the tilt angles using the TisKit library. The stations are ABALONES-4x4, the seismometers are Trillium Compact and the pressure sensors are differential pressure gauges.</caption><thead><tr><th>Station</th><th>Depth (m)</th><th>Start date</th><th>End date</th><th>Tilt angle (°)</th></tr></thead><tbody><tr><td>J09B</td><td>252</td><td>2012-09-02</td><td>2013-06-21</td><td>0.23</td></tr><tr><td>M11B</td><td>1109</td><td>2012-09-02</td><td>2013-06-18</td><td>-0.74</td></tr><tr><td>G02B</td><td>1920</td><td>2012-09-03</td><td>2013-06-19</td><td>-0.17</td></tr></tbody></table>

# 3 DATA SELECTION AND CONDITIONING

## 3.1 Data set and instrument selection

To evaluate our technique, we use data from the Cascadia Initiative (FDSN network 7D2011, IRIS OBSIP 2011; Toomey et al. 2014), because of its large range of water depths (from 100 to 5000 m), large sets of homogeneous instruments, and extensive previous analysis (e.g. Bell et al. 2015; Janiszewski et al. 2019; Stone et al. 2018; Tian & Ritzwoller 2017; Hilmo & Wilcock 2020). We wanted to use homogeneous instruments because instrument design is the first factor responsible for dissimilarities between PSDs at a site (Janiszewski et al. 2022). The instrument response of DPGs can slightly vary between two units, which can lead to dissimilarities in the PSDs between two stations (Doran et al. 2019). Such variations will have no impact on the denoising process. To concentrate on the compliance signal, we selected data which was minimally contaminated by tilt-induced seafloor current noise. We chose the BBOBS stations J09B, M11B and G02B, which share the same instrument design and period of acquisition and for which An et al. (2020) found no significant tilt-induced noise. These stations are respectively located at water depths of 252, 1019 and 1920 m. Their features are summarized in Table 2.

We checked the tilt noise on the stations by analysing the coherence between the seismometer channels on a day with low IG-wave energy (Fig. 1). The coherence values in the compliance frequency band are low (<0.4) at stations M11B and G02B, whereas coherence is as high as 0.85 in a narrow band at J09B. While this is not ideal for a pure test of compliance calculation, it allows us to observe the effect of tilt noise on the compliance calculation. The estimated tilt angles for these stations are all less than 0.17° (Table 2).

## 3.2 Data conditioning

We downloaded 1-d data files using the International Federation of Digital Seismograph Networks (FDSN) application programming interface, then concatenated the data to obtain 1-yr gapless data files. The data conditioning workflow consisted of anomaly detection, downsampling and instrument response removal. The time

spanned by the detected anomalies are set to zeros and tapered to prevent from any unstabilities during processing. Since the signal of interest lies below 50 mHz and the data were sampled at 50 Hz, we downsampled the data to 1 Hz after applying an anti-aliasing filter. We evaluated the instrument response in the frequency domain, using a water level of -60 dB, then removed the response using a deconvolution in the Fourier domain: details are available in Appendix B.

# 4 MULTISCALE DATA ANALYSIS

Signal stationarity is the main pre-requisite for frequency response transfer function-based seafloor compliance signal removal. The computation of the PSDs requires the signal to be stationary both over the period of each $N_w$-second spectral window and along the overall $N_T$ signal analysis period. Short-duration changes (< $N_w$) in signal statistics will make the Fourier transforms of each spectral window significantly different from one another. In this case, averaging their squared absolute values, as in eq. (A7), will lead to inaccurate PSD estimates because of the mix between the transient and the background signals at the frequencies where the transient signal occurs. Medium-duration ($N_w < \dots < N_T$) changes in signal statistics can also affect the averaging step at the lowest frequencies because of residual trends that are not well removed. Long-duration changes (> $N_T$) in signal statistics should not affect the estimation of the PSDs. In addition to, and sometimes in conflict with, the stationarity assumption, the durations used to calculate the PSDs must allow sampling of the longest IG-wave periods.

We propose an approach for assessing the time varying properties of pressure and vertical acceleration signals in order to control the stationarity properties of these signals and to optimize the selection of the window length ($N_w$) and the total signal analysis period ($N_T$). We propose a multiscale sliding deviation analysis, to assess second-order statistical variations across time over a variety of scales. The scales $N_\sigma$ are closely related to $N_w$ and $N_T$, as they define the length of each sliding segment on which the deviation is computed. They are therefore chosen to span a range covering the extreme possible values of $N_w$ and $N_T$.

We set a medium scale of 1 hr and other scales as power of two series up to 64 hr and down to 3 min and 45 s. These scales cover all the window and signal lengths that were used in previous studies, plus shorter lengths to account for transient events. The shortest duration (3 min and 45 s) is on the order of the maximum IG-wave period expected to be seen on the vertical component channel (Webb et al. 1991). The choice of the deviation metric is important, because different metrics can be sensitive to different parts of the signal. The signals are the combination of ‘background’ signals, such as IG waves and noise, and ‘foreground’ signals, such as transient events. Since the transient events are generally of higher amplitude than the background signal, and therefore easy to detect, we focus on the more subtle time-varying statistics of the background signal. The standard deviation (SD) metric is unsuited to reveal subtle changes in the background signals, because it is sensitive to outliers such as transient events. We choose the median absolute deviation (MAD), which is much less sensitive to transient events. The MAD is defined as:

$$ \sigma(t, N_{\sigma}) = \text{med}_{t-N_{\sigma}/2}^{t+N_{\sigma}/2} \left( |s(t) - \text{med}_{t-N_{\sigma}/2}^{t+N_{\sigma}/2}(s(t))| \right), \quad (3) $$

where $s(t)$ is the signal to be analysed and $\text{med}_{t-N_\sigma/2}^{t+N_\sigma/2}$ is the median value calculated between $t - N_\sigma/2$ and $t + N_\sigma/2$. For consistency
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Coherence between seismometer components on the M11B station on 2012 November 21 (N<sub>w</sub> =1800 s, N<sub>0</sub> =900 s, a 4π discrete prolate spheroidal window). The coherences are computed between the BHZ and BH1 + BH2 channels of station (a) M11B, (b) G02B and J09B. The straight line indicates the coherence significance at 95 per cent.</figcaption></figure>

with values provided by the SD metric, we apply a scaling factor of 1.4826 (Rousseeuw & Croux 1993) to σ(t, Nσ). To focus on the IG-wave frequency range, the multiscale sliding MAD is computed on signals bandpassed between 10⁻³ and 5×10⁻² Hz.

We analyse two 1-month time periods: one located at the end of autumn (2012 November 21 to December 21) and one in spring (2013 May 17 to June 17). We refer to them hereafter as the autumn and spring data sets. We will fully describe the analysis of M11B, located at intermediate depth, and summarize our analysis of stations G02B and J11B below and in the Supporting Information.

## 4.1 Station M11B

The multiscale deviations of station M11B are shown in Fig. 2. The pressure component (BDH) reveals a stronger IG-wave activity in autumn than in spring. Consequently, the amount of high amplitude signal on the pressure and vertical components is greater and the signal-to-noise ratio (SNR) is higher in autumn.

The differences in deviation across the scales are insignificant from half an hour to 32 hr for all dates and channels. For scales of

32 hr or longer, the deviations are smoother across time, especially for the 64 hr scale. Most of the deviations are stable for up to 24 hr (corresponding to the spacing between the dotted white lines in Fig. 2), but not more. Hence, the maximum period of analysis (N<sub>T</sub>) should not be longer than a day.

For scales shorter than half an hour, the days-long trend (stable deviations over 24 hr or longer), is still present, but three types of variations appear. The first type (type 1), are high amplitude, short and sparse variations, which are only present on the vertical component. They are most likely transient seismic events. The second type (type 2) seems to be periodic, suggesting that they are related to variations at the longest IG-wave periods. The third type (type 3) are high-amplitude variations, seen on both pressure and vertical components. The type 3 variations are better seen in smaller subsets (Fig. 3). Fig. 3 subsets highlight transition periods during which changes in signal statistics occur and exhibit the three different types. Type 2 variations only occur during strong IG-wave periods, supporting the assumption that they are due to IG waves. Type 3 variations are strong and not scale-dependent, which could be due to a temporary increase in meteorological activity. During
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Multiscale sliding MAD of station 7D.M11B. The deviations were computed on bandpassed time-series between 1 and 50 mHz. Left: deviations values for the autumn data set. Right: deviation values for the spring data set. The dotted white lines are spaced every 24 hr.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Multiscale sliding MAD subsets of station 7D.M11B. The deviations were computed on bandpassed time-series between 1 and 50 mHz. Left: deviations values during the autumn data set. Right: deviation values during the spring data set. The dotted white lines are spaced every 2 hr. The red ellipse shows an example of type 1 variations, the green ellipses show type 2 variations and the orange ellipses show type 3 variations.</figcaption></figure>

low-amplitude IG-wave periods, only type 1 (transient) events are seen, meaning that the background signal can be considered stationary.

The majority of the 2-hr periods (separated by white dotted lines in Fig. 3) contains significant variations of the deviation metric for scales from 0.0625 hr (225 s) to 1 hr. These periods are therefore not suitable to use as $N_w$ (window length) in the spectral analysis because the signal would not be stationary. A half-hour (1800 s) or quarter-hour (900 s) period is short enough to ensure that there would not be enough non-stationary windows in the PSD averaging step to significantly decrease the PSD estimate accuracy. Since we expect a maximum IG-wave period of 200 s, 15-min spectral windows (900 s) can resolve the IG-wave spectrum.

## 4.2 Station G02B

Station G02B is deployed at a deeper site (1920 m) than station M11B, so the IG-wave signal should be of lower amplitude and with a lower cut-off frequency. The station's multiscale sliding MAD values (Figs S1 and S2 of the Supporting Information) also show stronger IG-wave activity in autumn than in spring and more frequent high MAD values in autumn. The IG-wave amplitudes are lower than for station M11B but the MAD values are similar, indicating that the IG-wave signal power depends strongly on depth but that its overall variability does not. The result is similar for the vertical component. The parameter analysis for G02B is similar to that for M11B but, since the IG-wave signal power and the cut-off frequency are lower, the short scales do not show
<figure><img src="image_5.png" /><figcaption>**Figure 4.** M11B station day-to-day vertical-pressure coherence for the (a) autumn and (b) spring periods.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** G02B station day-to-day vertical-pressure coherence for the (a) autumn and (b) spring periods. The upper limit of the frequency scale is slightly different between (a) and (b).</figcaption></figure>

peak values during strong IG-wave periods. This smooths the type 3 variations and relaxes the signal stationarity, allowing us to use longer $N_w$.

## 4.3 Station J09B

Station J09B is deployed in very shallow water (252 m), so its IG-wave signal and cut-off frequency should be higher than at the other sites. The station's multiscale sliding MAD values (Figs S3 and S4 of the Supporting Information) correlate visually with those at the other stations for the autumn data set, but not for spring. During the spring period, the vertical channel records periods of high energy with no correlated IG signal (e.g. between 2013 May 17 and 21). This should lower the coherence between the two channels and make the transfer function frequency response harder to estimate during

the spring season. The parameters chosen for station M11B are also valid for J09B.

## 5 ROBUST COMPLIANCE NOISE REMOVAL

In this section, we propose a new procedure to more accurately estimate the scTFFR and to improve compliance noise removal from the vertical component. We pay particular attention to expanding the frequency range over which compliance scTFFR is estimated, which pays off especially at the lowest frequencies. We first carefully select the optimum parameter values for the PSD estimation ($N_w$, $N_o$ and $N_T$) and the taper window ($N_w$ and $N_T$). We then present an automated methodology for the window selection. Finally, we propose a month-stacking procedure for the daily-estimated scTFFR to increase its stability.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** J09B station day-to-day vertical-pressure coherence for the (a) autumn and (b) spring periods.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** M11B Monthly estimated scTFFRs magnitude and unwrapped phase values for the autumn (left) and spring (right) data sets.</figcaption></figure>

In order to focus on the improvements obtained using this method, we do not remove tilt, nor known earthquakes, before our analysis. These operations may be separately applied before using our method.

## 5.1 Optimum parameter value selection

The spectral analysis relies on the values of $N_w$, $N_T$ and $N_o$. As we determined in Section 4, $N_T$ should not be greater than 32 hr. We choose a conservative value of 24 hr. The $N_w$ parameter directly impacts the spectral resolution and uncertainty: as $N_w$ increases, the frequency resolution increases but the number of windows to average decreases. Furthermore, the longer $N_w$, the more samples will be discarded due to transients, but the better will be compliance calculation/removal at the lowest frequencies. Our analysis (Section 4) indicates that $N_w$ should not be greater than half an hour (1800 s). The expected maximum IG-wave period is approximately 200 s and so $N_w$ should not be less than 400 s, to ensure that these

long-period IG waves will be correctly sampled. Given these two limits, we select $N_w = 900$ s.

The taper window type does not have a significant influence on the scTFFR estimation. We choose the discrete prolate spheroidal window (Slepian 1978), with a standardized half bandwidth equal to four samples (Crawford & Webb 2000). We set $N_o = N_w/2$, to equally cover all of the samples in the signal without significant redundancy.

## 5.2 Spectral window selection

### 5.2.1 Transient event detection

The signals may contain transient events that are non-stationary in the IG-wave frequency range. A detection stage is required to eliminate these events. We focus of detecting transient signals that are stronger than the background IG waves, under the assumption that transient events that are less powerful than the background
<figure><img src="image_9.png" /><figcaption>**Figure 8.** G02B monthly estimated scTFFRs magnitude and unwrapped phase values for the autumn (left) and spring (right) data sets.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** J09B monthly estimated scTFFRs magnitude and unwrapped phase values for the autumn (left) and spring (right) data sets.</figcaption></figure>

IG-wave signal do not significantly affect the of the overall signal stationarity. To detect transient events, we compute deviation values over a sliding window. We use the SD algorithm here, to focus on foreground signals. The higher amplitude transient signals in the IG-wave frequency range last from 60 to 1000 s, so we choose a 225 s (1/16 hr) window, which is the lowest scale evaluated in Section 4. To account for any long-term changes in the signal amplitude, the sliding SD is divided into 2-hr segments, with the median value for each segment multiplied by a factor *a* to determine the event threshold. The factor *a* value should be set between 2 (tight) and 5 (very loose). Any sliding SD value above this threshold is labelled as a transient event and any window containing transient events is excluded from spectral analysis.

### 5.2.2 Correlation criterion

We add a supplementary criterion for the automated window selection, based on the correlation coefficient between the pressure and vertical component signals. The final coherence between the

vertical and pressure component signals is strongly related to their correlation. Some spectral windows may have low correlation in the IG-wave frequency range, bringing down the coherence values. The following methodology allows us to eliminate these windows from spectral analysis. First, the signal windows are filtered between 3.5 and 50 mHz. Then, the unbiased correlation coefficient is computed between each pair of windows. If the absolute value of the correlation coefficient is less than 0.5, the windows are discarded.

## 5.3 Transfer function frequency response stacking

In past studies, the scTFFR and the noise-free vertical component signal were estimated on a day-to-day basis in the frequency range where the coherence value is above the 95 per cent significance value. We propose to extend compliance estimation to a monthly basis, assuming that the seafloor compliance is stable across time.

The monthly scTFFR estimate, hereafter referred to as the monthly stacked scTFFR, is the average of the 24-hr scTFFR estimates, stacked in the largest frequency range found during
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Day-to-day PSDs of the original vertical acceleration signal of station M11B for the (a) autumn and (d) spring data sets and for the vertical signal acceleration corrected from the compliance component using a day-to-day estimated scTFFR for the (b) autumn and (e) spring data sets and the monthly stacked scTFFR for the (c) autumn and (f) spring data sets. The corrections were carried out in the frequency range defined by the day-to-day coherence significance threshold for the day-to-day correction. For the monthly stacked correction, it was carried out using the best frequency range found of the season related day-to-day correction. The PSDs have been estimated using a 4-pi spheroidal window with $N_w = 900$ s, $N_o = 450$ s and a median average.</figcaption></figure>

the month. The denoising process is then applied to each day's data using the monthly stacked scTFFR. By applying the monthly stack scTFFR to each day's data, we take advantage of the monthly scTFFR's lower uncertainty and constant frequency range, resulting in a more efficient and stable noise reduction.

The assumption that compliance is stationary over a month's time could be false if the subsurface changes rapidly, such as over a volcano in the period surrounding an eruption (Doran & Crawford 2020), or if the subsurface structure is strongly 2-D and the IG waves change directions over multiday periods (Crawford et al. 1998). This can be evaluated by comparing the formal uncertainty of the stacked estimates with the variance between the daily values, with a significantly higher variance indicating that compliance is not stationary. In this case, the monthly estimate could be replaced by a shorter, N-day estimate that satisfies stationarity conditions.

## 5.4 Performance assessment

We use two methodologies to quantify how much of the IG-wave signal is removed from the vertical component. First we calculate the deviation of the daily signal and evaluate how much it decreases and stabilizes around what should be a compliance-free background vertical acceleration value. Since the background signal is of interest, we use the σ MAD metric. Second, we compare the performance

of the proposed approach to that of the reference method described by Crawford & Webb (2000), hereafter referred to as the 'classical method'. To ensure an unbiased comparison of the two approaches, the classical method was applied to the same preprocessed data, including quality control, and used the same parameter values ($N_w$, $N_o$ and window type).

# 6 RESULTS

## 6.1 Coherences and transfer function frequency responses

We processed the autumn and spring data'sets of the three stations using the classical and proposed denoising methodologies. We computed daily coherences between the vertical and pressure component signals, using the classical method. On the M11B station, there is a clear difference in coherence between the autumn and spring seasons (Fig. 4). The frequency band of the autumn coherence is wider and the values higher than in the spring, indicating that the IG-wave amplitudes are greater in autumn. This is consistent with the higher level of deviation in autumn (Fig. 3). This is also true at the G02B station, except that the maximum coherence frequency bands are the same for both seasons (Fig. 5). At station J09B, the coherence levels are approximately the same in both seasons (Fig. 6).
<figure><img src="image_12.png" /><figcaption>Figure 11. Day-to-day PSDs of the original vertical acceleration signal of station G02B. More details can be found in the caption of Fig. 10.</figcaption></figure>

More windows are rejected in spring than in autumn at stations M11B and G02B (Figs 4 and 5), meaning that more transient events are detected during spring and/or the correlation between the two channels is more often below the threshold value. Both of these are probable consequences of the lower IG-wave amplitude in spring.

The monthly stacked scTFFR amplitude for each station is not significantly different between spring and autumn (Figs 7–9). The phase values lie around $\pi$ radians, consistent with the physics of seafloor compliance. Uncertainties in the DPG instrument response may affect the estimated compliance (Godin et al. 2013) but will not affect the IG noise removal.

## 6.2 Performance

We now analyse the performance of the IG-wave component correction, first for the M11B station and then for the G02B and J09B stations. We first apply the classical correction method using the day-to-day scTFFR to both autumn and spring data sets in the frequency range defined by the 95 per cent significance coherence threshold for the given day. We then apply our proposed approach for each day using the monthly stacked scTFFR of the related season in the best frequency range found by the classical method during the corresponding month. Since the instrument response has been removed on the year-long data set, the resolution of the low-frequency content is much improved than if the instrument response was removed for only one day. Therefore, complex large-scale trends are

present after the instrument response removal when considering only 1-d-long segments. To stabilize the noise removal process in the lower frequencies, we found that polynomials with an order of at most 10 are capable of removing these trends. Consequently, each day-long pressure and acceleration signals are detrended using an appropriate 10th-order polynomial. Moreover, to ensure Fourier operation stability, the edges of the signals are tapered using a 2-min-long half-shaped taper window. Since the amount of tapered data is minimal compared to the whole signal, the shape of the window does not impact the following Fourier operation. The PSDs of the original and compliance-corrected vertical component signals are shown in Figs 10–12 for stations M11B, G02B and J09B. They have been estimated using a 4-pi spheroidal window with $N_w = 900$ s and $N_o = 450$ s. We use the 4-pi spheroidal window to conform to Crawford & Webb (2000). The frequency range covered by the IG-wave component variations are significant during both seasons for station M11B, which is consistent with the coherence seen in Fig. 4. Due to the lack of correction and the less stable daily scTFFR under $8 \times 10^{-3}$ Hz, low-frequency IG-wave components remain. This is seen during both seasons for station M11B, but the reduced frequency range is more pronounced during spring. These remaining IG-wave components are greatly reduced in the M11B signal corrected by the monthly stacked scTFFR and the PSDs of the signals corrected by the monthly stacked scTFFR are more consistent across time. This is especially true for the autumn season. The correction seems to remove slightly more IG-wave component using the monthly stacked scTFFR than the day-to-day one between
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Day-to-day PSDs of the original vertical acceleration signal of station J09B. More details can be found in the caption of Fig. 10.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** Day-to-day PSDs of the original and corrected vertical acceleration signal for (a) 2012 December 19 and (b) 2013 May 29 for station M11B.</figcaption></figure>

10⁻² and 3×10⁻² Hz, especially around 2012 November 24 and 29. Some high-frequencies IG-wave component are also more accurately removed above 3.5×10⁻² Hz, around 2013 May 21 for instance. As we could have expected considering the coherence level and the amount of rejected windows, the corrected spring data set signals show much noisier PSD levels between 10⁻² and 3×10⁻² Hz, even when using the monthly stacked scTFFR. The noisier PSD comes from 2013 May 24. This is mainly due to the occurrence of a magnitude 8.4 earthquake at 5:44:48 which is recorded by the

M11B station. Although its decay is long enough to affect the median averaging for the visualization PSD estimation, the IG-wave component was successfully removed. This is due to both the automated event detection procedure and the window cross-correlation selection. Because the event detection method divides the signal in 2 hr segments, only the part of the earthquake signal that has the higher magnitude level is labelled as an event. But, the cross-correlation coefficient value between the bandpassed pressure and vertical acceleration windows containing the earthquake is not high
<figure><img src="image_15.png" /><figcaption>**Figure 14.** 20-min bandpassed segments of corrected vertical acceleration signal (top), extracted IG-wave acceleration signal (middle) and original negative pressure signals (bottom) for 2012 December 19 (left) and 2013 May 29 (right) for station M11B.</figcaption></figure>

<table><caption>**Table 3.** Mean deviation values (MAD) for 30 d in the autumn data set (from 2012 November 21 to December 21) and 31 d in the spring data set (from 2013 May 17 to June 17). The *o*, *d* and *s* indices stand respectively for original, day-to-day corrected and monthly stacked corrected vertical component data. The overline symbol $\bar{x}$ stands for the average over all days in the corresponding periods. The best values are shown in bold.</caption><thead><tr><th>Data set</th><th>σ<sub>o</sub> (m s<sup>-2</sup>)</th><th>σ<sub>d</sub> (m s<sup>-2</sup>)</th><th>σ<sub>s</sub> (m s<sup>-2</sup>)</th></tr></thead><tbody><tr><td>Autumn</td><td>1.03 × 10<sup>-7</sup></td><td>8.90 × 10<sup>-9</sup></td><td><strong>6.17 × 10<sup>-9</sup></strong></td></tr><tr><td>Spring</td><td>4.87 × 10<sup>-8</sup></td><td>1.15 × 10<sup>-8</sup></td><td><strong>8.39 × 10<sup>-9</sup></strong></td></tr></tbody></table>

enough. Consequently, the rest of the event is not kept during the spectral analysis, which ensured an accurate scTFFR estimate. The same analysis can be drawn for the station G02B (Fig. 11), the main difference between the stations G02B and M11B is the magnitude of the original PSD, which is greatly lower for the station G02B due to its deeper location. The lower frequency IG wave remaining components are still present after the compliance correction using the classical method and they are accurately removed by the proposed method. The corrected PSDs for station J09B (Fig. 12) are however not different between the proposed and the classical method. Some days even show a poorer corrected PSD using the monthly stacked scTFFR (2012 November 28). We assume that this is due to the shallower location of this station. The background signal is certainly not only populated by the IG waves but also by microseismicity which amplitude can change across the month and affect the monthly stacked scTFFR estimation. We hereafter conduct a more detailed analysis and only focus on station M11B.

To illustrate the improvement of the proposed approach, more details of the compliance signal removal results are shown for 2012 December 19 and 2013 May 29. The PSDs of the original and corrected signals for these days are shown in Fig. 13. While the correction results in the PSDs is equivalent for both methods above $10^{-2}$ Hz, the proposed approach using the monthly stacked scTFFR outperforms the classical method under $10^{-2}$ Hz, especially for 2012 December 19. This improvement can also be seen in the time domain. Fig. 14 shows two 20-min segments of the corrected acceleration signals, the extracted IG waves acceleration component and the original negative pressure signals. The acceleration signals corrected using the monthly stacked scTFFR do not have long-period components unlike the signals corrected using the classical day-to-day scTFFR. The acceleration signal corrected with the day-to-day scTFFR is then of higher amplitude than the one corrected with the monthly stacked scTFFR. Finally, the mean deviation value (MAD) of the original and vertical acceleration signal are given in Table 3. The best performance is achieved when using the monthly stacked scTFFR with a reduction of the MAD of 94 and 82.8 per cent, respectively, for the autumn and spring data set, while the reduction is of 92.2 and 79.7 per cent when using the daily scTFFR, respectively, for the autumn and spring data set. The mean deviation of the compliance signal corrected using the daily scTFFR is reduced by 22.8 and 15.2 per cent by using the monthly stacked scTFFR, respectively, for the autumn and spring data set. The decrease of deviation is less important for the spring data set due to the fact that the magnitude of the IG waves during this season is lower than during the autumn season.
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Day-to-day PSDs of the original vertical acceleration signal of station M11B for the (a) spring season and for the (b) vertical signal acceleration corrected from the compliance component using monthly spring stacked scTFFR and (c) the monthly autumn stacked scTFFR. The PSDs have been estimated using a 4-pi spheroidal window with $N_w = 900$ s, $N_o = 450$ s and a median average.</figcaption></figure>

## 6.3 Limitations

This study was carried out on only three stations with a focus on the M11B station whose depth is low enough to let a substantial amount of IG waves to reach the seafloor without overlapping with the microseismic frequency range. A study conducted on the entire Cascadia Initiative network would overcome this uncertainty. Moreover, we restricted our studies to $N_w = 900$ s which is prior knowledge from the multiscale deviation analysis. We believe that an optimum value of $N_w$ can be found between 900 and 1800 s as 1800 s is seen as the limit of stationarity. Also, the deviation information can be further used during the compliance signal removal process as a stationarity index, making possible to detect non-stationary windows and discard them. It should also be noted that the performance are blindly assessed, meaning no true vertical component noiseless signal are available to precisely quantify the noise removal process. Further investigations are required to quantify precisely the denoising performance with synthetic data. Such data can be generated using eqs (A1) and (A2) for the IG-wave modelling and a seafloor compliance model (Crawford 2004). This study also focus on high-coherence pressure versus vertical component data sets. It should be interesting to quantify the denoising performance data sets showing more frequent lower coherence values.

## 6.4 Transferability of the monthly stacked scTFFR

It is of interest to evaluate the transferability of the monthly stacked scTFFR from the autumn season to the spring season to check if the mantle crust has evolved between these seasons (Doran & Crawford 2020; Crawford et al. 1998; Webb & Crawford 1999) and if we can use a better estimated monthly stacked scTFFR on another data set. To this end, the corrected PSDs of station M11B and G02B during the spring season, using both monthly autumn and spring

stacked TFFR, are shown in Figs 15 and 16. On the one hand, for the M11B station, both scTFFRs share the same frequency range and yield to fairly identical corrected PSDs. On the other hand, the G02B stations show a wider coherent frequency range when computing the stacked scTFFR during the autumns than during the spring seasons. This is due to IG waves having more frequency content during the autumn season. Removing the compliance signal from the spring data set using the autumn scTFFR slightly improves the SNR below 5 mHz. Two conclusions can be drawn from this result. The first conclusion would be that the compliance response of the sediment remains unchanged between 2012 autumn and 2013 spring seasons. The second one is that a different amount of IG waves accross the seasons can lead to slightly different robustness of the transfer function estimate and prevents from completely removing the compliance signal. Even if the robustness in the transfer function estimate is not systematic, it is relevant to remove the compliance signal using the scTFFR estimated during the period having the strongest IG waves that cover a larger frequency range.

## 7 CONCLUSION

This study focused on improving the seafloor compliance noise removal method, which is based on the frequency response transfer function between a BBOBS's vertical and pressure channels. We first conducted a multiscale deviation analysis to assess the validity of the critical assumption, implicit in the calculation of the scTFFRs, that the two signals are stationary. This analysis allows us to calculate the range of values of the key parameters, in particular the window length and the overall analysis duration, over which the stationarity assumption is valid. Using data from broad-band seafloor stations from the Cascadia Initiative, we found that the upper limit for the analysis duration is 32 hr and the valid range of window lengths is 15-30 min.
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Day-to-day PSDs of the original vertical acceleration signal of station G02B for the (a) spring season and for the (b) vertical signal acceleration corrected from the compliance component using monthly spring stacked scTFFR and the (c) monthly autumn stacked scTFFR.</figcaption></figure>

We propose a multistage approach to improve estimates of seafloor compliance and to remove its related noise. An automated window selection step discards windows containing transient events or whose vertical-pressure cross-correlation coefficient is below a threshold value. Next, by stacking daily scTFFR estimates, we increase the SNR of the corrected signal. These long-term compliance estimates also allow us to uniformly remove compliance noise over all of the days in the period, including days with low pressure-vertical coherence. When applied to the Cascadia Initiative stations M11B, G02B and J09B, our method visibly improves noise removal at the lowest compliance frequencies. We also demonstrated that the compliance was stationary over the year measurements in our test data and that, in this case, compliance estimated during a period with good pressure-acceleration coherence can be used to remove compliance noise from other periods.

## ACKNOWLEDGMENTS

The authors would like to deeply thank Martin Schimmel from Geosciences Barcelona (GEO3BCN - CSIC), Barcelona and Éleonore Stutzmann from Institut de Physique du Globe de Paris, as well as all the members of the ANR BRUIT-FM project for their valuable expertise and time. This work was funded by the BRUIT-FM project sponsored by the French National Research Agency (ANR-21-CE01-0031).

## SUPPORTING INFORMATION

Supplementary data are available at GJIRAS online.

supplementary.tar.gz

Please note: Oxford University Press are not responsible for the content or functionality of any supporting materials supplied by

the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

## DATA AND SOFTWARE AVAILABILITY

The Cascadia Initiative station data are open and available for download from the IRIS website. The BRUIT-FM Python toolbox developed for the purpose of this study is available at https://gitlab.ifremer.fr/anr-bruitfm/bruit-fm-toolbox.

## REFERENCES

* An, C., Wei, S.S., Cai, C. & Yue, H., 2020. Frequency limit for the pressure compliance correction of ocean-bottom seismic data, *Seismol. Res. Lett.*, **91**(2A), 967–976.

* Aucan, J. & Ardhuin, F., 2013. Infragravity waves in the deep ocean: an upward revision, *Geophys. Res. Lett.*, **40**(13), 3435–3439.

* Bell, S.W., Forsyth, D.W. & Ruan, Y., 2015. Removing noise from the vertical component records of ocean-bottom seismometers: results from year one of the cascadia initiative, *Bull. seism. Soc. Am.*, **105**(1), 300–313.

* Bendat, J.S. & Piersol, A.G., 2011. *Random data: Analysis and Measurement Procedures*, Chapter 6: Single-Input/Output Relationships, p. 407, Wiley-Interscience, New York, 4th edn.

* Biésel, F., 1952. Equations générales au second ordre de la houle irrégulière, *La Houille Blanche*, **38**(3), 372–376.

* Crawford, W., 1994. *Determination of Oceanic Crustal Shear Velocity Structure from Seafloor Compliance Measurements*, Ph.D. thesis, University of California, San Diego.

* Crawford, W.C., 2004. The sensitivity of seafloor compliance measurements to sub-basalt sediments, *Geophys. J. Int.*, **157**(3), 1130–1145.

* Crawford, W.C. & Webb, S.C., 2000. Identifying and removing tilt noise from low-frequency (<0.1 Hz) seafloor vertical seismic data, *Bull. seism. Soc. Am.*, **90**(4), 952–963.
Crawford, W.C., Webb, S.C. & Hildebrand, J.A., 1991. Seafloor compliance observed by long-period pressure and displacement measurements, *J geophys. Res.*, **96**(B10), 16151, doi:10.1029/91JB01577.

Crawford, W.C., Webb, S.C. & Hildebrand, J.A., 1998. Estimating shear velocities in the oceanic crust from compliance measurements by two-dimensional finite difference modeling, *J. geophys. Res.: Solid Earth*, **103**(B5), 9895–9916.

Crawford, W.C., Satish C. Singh, T.H. & Smallwood, J.R., 2005. Applications of seafloor compliance measurements in the faroes-shetland basin, in Ziska, Thomas. *Proceedings from Faroe Islands Exploration Conference 2004*, Vol. **43**, pp. 32–43. Faroe University Press.

Doran, A.K. & Crawford, W.C., 2020. Continuous evolution of oceanic crustal structure following an eruption at axial seamount, juan de fuca ridge, *Geology*, **48**(5), 452–456.

Doran, A.K. & Laske, G., 2019. Seismic structure of marine sediments and upper oceanic crust surrounding hawaii, *J. geophys. Res.: Solid Earth*, **124**(2), 2038–2056.

Doran, A.K., Rapa, M., Laske, G., Babcock, J. & McPeak, S., 2019. Calibration of differential pressure gauges through in situ testing, *Earth Space Sci.*, **6**(12), 2663–2670.

Dougherty, E.R., 1992. *An Introduction to Morphological Image Processing*, Tutorial Texts in Optical Engineering, SPIE Optical Engineering Press, Bellingham, Washington, 2nd edn., Literaturverz. S. 144–160.

Godin, O.A., Zabotin, N.A., Sheehan, A.F., Yang, Z. & Collins, J.A., 2013. Power spectra of infragravity waves in a deep ocean, *Geophys. Res. Lett.*, **40**, 1–7.

Hilmo, R. & Wilcock, W.S.D., 2020. Physical sources of high-frequency seismic noise on cascadia initiative ocean bottom seismometers, *Geochem. Geophys. Geosyst.*, **21**(10), doi:10.1029/2020GC009085.

IRIS OBSIP, 2011. Cascadia initiative community experiment - obs component. https://doi.org/10.7914/SN/7D_2011

Janiszewski, H.A., Daherty, J.B., Abers, G.A., Gao, H. & Eilon, Z.C., 2019. Amphibious surface-wave phase-velocity measurements of the Cascadia subduction zone, *Geophys. J. Int.*, **217**(3), 1929–1948.

Janiszewski, H.A., Eilon, Z., Russell, J.B., Brunsvik, B., Daherty, J.B., Mosher, S.G., Hawley, W.B. & Coats, S., 2022. Broadband ocean bottom seismometer noise properties, *Geophys. J. Int.*, doi:10.1093/gji/ggz051

McNamara, D.E. & Buland, R.P., 2004. Ambient noise levels in the continental United States, *Bull. seism. Soc. Am.*, **94**(4), 1517–1527.

Rousseeuw, P.J. & Croux, C., 1993. Alternatives to the median absolute deviation, *J. Am. Stat. Assoc.*, **88**(424), 1273–1283.

Slepian, D., 1978. Prolate spheroidal wave functions, fourier analysis, and uncertainty-V: the discrete case, *Bell Syst. Tech. J.*, **57**(5), 1371–1430.

Stone, I., Vidale, J.E., Han, S. & Roland, E., 2018. Catalog of offshore seismicity in cascadia: Insights into the regional distribution of microseismicity and its relation to subduction processes, *J. geophys. Res.: Solid Earth*, **123**(1), 641–652.

Thompson, R. O. R.Y., 1979. Coherence significance levels, *J. Atmos. Sci.*, **36**(10), 2020–2021.

Tian, Y. & Ritzwoller, M.H., 2017. Improving ambient noise cross-correlations in the noisy ocean bottom environment of the Juan de Fuca plate, *Geophys. J. Int.*, **210**(3), 1787–1805.

Toomey, D.R. et al., 2014. The cascadia initiative: a sea change in seismo-logical studies of subduction zones, *Oceanography*, **27**(2), 138–150.

Webb, S.C., 1998. Broadband seismology and noise under the ocean, *Rev. Geophys.*, **36**(1), 105–142.

Webb, S.C. & Crawford, W.C., 1999. Long-period seafloor seismology and deformation under ocean waves, *Bull. seism. Soc. Am.*, **89**(6), 1535–1542.

Webb, S.C., Zhang, X. & Crawford, W., 1991. Infragravity waves in the deep ocean, *J. geophys. Res.*, **96**(B10), 16151, doi:10.1029/91JB01577.

Welch, P.D., 1967. The use of fast fourier transform for the estimation of power spectra: A method based on time averaging over short, modified periodograms, *IEEE Trans. Audio Electroacoust.*, **15**(2), 70–73.

Willoughby, E.C., Latychev, K., Edwards, R.N., Schwalenberg, K. & Hyn-dman, R.D., 2008. Seafloor compliance imaging of marine gas hydrate deposits and cold vent structures, *J. geophys. Res.*, **113**, B07107. doi: 10.1029/2005JB004136.

Yang, Z., Sheehan, A.F., Collins, J.A. & Laske, G., 2012. The character of seafloor ambient noise recorded offshore New Zealand: Results from the MOANA ocean bottom seismic experiment, *Geochem. Geo-phys. Geosyst.*, **13**(10), n/a–n/a, doi:10.1029/2012GC004201.

# APPENDIX A: BACKGROUND THEORY

## A1 Infra-gravity wave signal and seafloor compliance

For more than two decades, highly sensitive BBOBSs (Webb 1998) have been used to record subtle seismological signals below 0.2 Hz, including microseism and IG waves. These waves are generated by nonlinear interactions between ocean surface waves, mostly near coastlines, and propagate into the deep ocean with little attenuation (Webb et al. 1991). IG waves also include waves directly generated by ocean swells (Biésel 1952). Their seafloor signal drops dramatically at frequencies corresponding to wavelengths that are less than the water depth (Bell et al. 2015; Crawford et al. 1998), making them clearly recognizable in pressure PSDs. In shallow water (< 500 m), IG waves may occupy some of the same frequency range as microseisms, making them more difficult to recognize in pressure PSDs. This frequency-depth dependence is easily seen in the formulation (Webb et al. 1991):

$$P_b(k) = \frac{P_0}{\cosh(kD)} = \frac{a\rho g}{\cosh(kD)}, \quad (A1)$$

where $k$ is the wavenumber ($\text{rad} \times \text{m}^{-1}$), $P_0$ is the pressure magnitude (Pa) at the top of the water column, $D$ is the depth beneath the sea surface, $\rho$ is the seawater surface density (kg m$^{-3}$) and $g$ is the Earth's gravitational acceleration (m s$^{-2}$). The wavenumber $k$ is linked to the frequency $f$ using the dispersion equation:

$$(2\pi f)^2 = gk \tanh(kD). \quad (A2)$$

The average height of IG waves is about 10 mm and the wavelengths of IG waves detectable at the seafloor cover a range from $D$ to $\approx 300\sqrt{(g \times D)}$ (2–40 km for a 2 km water depth) (Crawford et al. 2005).

The seafloor deforms under the pressure fluctuations generated by IG waves (Crawford et al. 1998). The vertical deformation does not exceed 10 µm but is enough to be recorded by the BBOBS. Seafloor compliance is defined as the ratio of the differential seafloor vertical displacement to the seafloor pressure.

## A2 Conditions for estimating the transfer function frequency response

We summarize here the conditions under which we can relate signals to seafloor compliance through frequency-domain estimates. We assume no seafloor current-induced noise. In the frequency range of the IG waves, the seafloor compliance relation is considered linear and characterized by its impulse response $h(t)$. This relation links the clean (compliance-free) vertical acceleration $z(t)$, the recorded acceleration $z_c(t)$ and the pressure $p(t)$ signals through:

$$z_c(t) = z(t) + p(t) * h(t), \quad (A3)$$

where * denotes the convolution product.

We can rewrite eq. (A3) to express the clean signal in the spectral domain, using capital letters to represent a spectral representation of signals, possibly obtained over a limited time range and windowed:

$$Z(f) = Z_c(f) - P^\dagger(f)H(f), \quad (A4)$$

where $f$ is the frequency and $\dagger$ denotes the complex conjugate.
<table><caption>**Table A1.** Description of BRUIT-FM toolbox parameters.</caption><thead><tr><td>Parameter</td><td>Description</td></tr></thead><tbody><tr><td>N<sub>w</sub></td><td>Length (seconds) of the analysis window</td></tr><tr><td>N<sub>o</sub></td><td>Length (seconds) of the overlap between two windows</td></tr><tr><td>N<sub>T</sub></td><td>Total duration (seconds) of the signals to be analysed</td></tr><tr><td>Taper type</td><td>Shape of the taper window, for example, Hann, Hamming and discrete prolate spheroidal</td></tr></tbody></table>

An optimal estimate of $H(f)$ (written with hat $\hat{\ }$) can be obtained, assuming that all three signals are second-order stationary. We introduce the notion of cross-spectrum for signals $x$ and $y$: $G_{xy}(f) = X^{\dagger}(f)Y(f)$. We obtain $G_{z_c p}(f)$ and the individual (auto-)power spectra $G_p(f)$ and $G_{z_c}(f)$. If the compliance system is time-invariant, and $z(t)$ and $p(t)$ are uncorrelated, an optimal scTFFR is estimated as the Wiener filter (Bendat & Piersol 2011):

$$ \hat{H}(f) = \frac{G_{z_c p}(f)}{G_p(f)}. \qquad (A5) $$

In practice, one often improves this estimate by splitting the pressure and acceleration signals (of total length $N_T$ seconds) into $K$ windows of length $N_w$, with an overlap of $N_o$ seconds. These parameters are summarized in Table A1. Each signal is apodized by a suitable taper window. These are the key parameters used in the BRUIT-FM toolbox and are summarized in Table A1.

Then, through Welch's modified periodogram (Welch 1967), by averaging windowed cross and power spectra denoted by $l$, their window index, a novel estimate of scTFFR is obtained with:

$$ \hat{G}_{z_c p}(f) = \frac{1}{L} \sum_{l=1}^{L} Z_{cl}^{\dagger}(f) P_l(f), \qquad (A6) $$

$$ \hat{G}_p(f) = \frac{1}{L} \sum_{l=1}^{L} P_l^{\dagger}(f) P_l(f). \qquad (A7) $$

The length $N_w$ of windowed segments drives the frequency resolution and the SNR of the estimates. In a nutshell, the larger the segment, the better the frequency resolution and the lower the SNR. A minimal frequency resolution (larger $N_w$ s) is required to accurately evaluate the IG waves bandwidth. However, second-order stationarity is not granted in practice. To ensure that the signal variance does not significantly change over time, smaller $N_w$ s are preferred. A compromise must then be found.

Several practical condition mismatches can be acknowledged. First, $z(t)$ and $p(t)$ may not be fully uncorrelated in practice, for example, a low-magnitude/frequency earthquake can be recorded by both pressure instruments and seismometers. A selection step is usually conducted to remove the data contaminated by any earthquake signals from the spectral analysis. Nevertheless, correcting the vertical acceleration signal while earthquake signals are present in the pressure signal leads to distortions (Bell et al. 2015; An et al. 2020). Secondly, $z_c(t)$ and $p(t)$ may not carry enough information to estimate $\hat{H}(f)$ properly. A common procedure is to estimate $\hat{H}(f)$ in a reduced spectral range, with the highest level of coherence $C_{pz_c}(f) \in [0, 1]$:

$$ C_{pz_c}(f) = \frac{|G_{pz_c}(f)|}{\sqrt{G_{pp}(f)G_{z_c z_c}(f)}}. \qquad (A8) $$

One can also define the statistical coherence significance (Thompson 1979):

$$ c_{\alpha} = 1 - (1 - \alpha)^{\frac{1}{N-1}}, \qquad (A9) $$

where $\alpha$ denote the confidence sought in the statistical result, also known as p-value. The spectral range is then usually defined for any coherence value above a 95 per cent significance threshold. This can restrain the process to a limited frequency range of the IG-wave spectrum or even make it impossible to estimate the scTFFR if the coherence level is too low.

## APPENDIX B: DATA CONDITIONING

They come as raw data designated as $r_c(t)$, where $c$ stands for the channel number, and need to be preprocessed before trying to extract any features or information from it.

### B.1 Anomaly detection

The first pre-processing stage is the detection and tagging of anomalies like saturated amplitudes or zero-valued data. The invalid data time spans are diluted by 60 s using a binary morphological mathematics tool (Dougherty 1992) and convolved with a 60-s taper window. Given the few amount of anomaly in the data, the shape of the taper window should not have any impact on the signal, hence we used the classic Hamming window. The signal is multiplied by the resulting convolved mask to smooth the discontinuities. This operation also tapers the edges of the signal, which improves the stability of the forthcoming filtering operation.

### B.2 Downsampling

Since the signal of interest lies below 50 mHz and the data are usually sampled at 50 Hz, we downsample them to 1 Hz to reduce data storage consumption and computation cost. Prior to downsampling, an anti-aliasing filter is applied. We need a filter with a flat response in the passband. A nonlinear phase shift is not a problem since the filter can be applied both forward and backward, cancelling the phase component of the filter. Two filter designs are well suited for this task: the Butterworth filter, which has a smooth transition between the passband and the stopband, and the type 2 Chebyshev filter which has a sharper gain transition between the passband and the stopband, and has ripples in the stopband. Both designs complete the task well and their produced results would be close to each other. We use a fifth-order low-pass Butterworth design with a cut-off frequency of 0.4 Hz. The filter is performed on an extended version of the original signal, whose edges have been mirrored over an eighth of the signal to mitigate stability issues. The filtered data are then linearly interpolated in the time domain to the targeted sampling frequency.

### B.3 Instrument response removal

The downsampled instrument response is evaluated in the frequency domain for the month-long signal. A water level of -60 dB is applied and a deconvolution is performed in the Fourier domain.