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

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.

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.¹ 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 (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 BBOBSs (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 |$\mu {\rm 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 |$H(f)$| (hereafter noted scTFFR ) can be obtained following the procedure of Bendat & Piersol (2011) which states that:

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 |$\widehat{H}(f)$| properly. A common procedure is to estimate |$\widehat{H}(f)$| in a reduced spectral range in which with coherence, defined as:

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 a 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.

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 homogenous 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 (⁠|$\lt $|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|$^{\circ }$| (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 (⁠|$\lt 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 \lt \cdots \lt 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 (⁠|$\gt 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:

where |$s(t)$| is the signal to be analysed and |$\mathrm{ 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 with values provided by the SD metric, we apply a scaling factor of 1.4826 (Rousseeuw & Croux 1993) to |$\sigma (t,N_\sigma )$|⁠. To focus on the IG-wave frequency range, the multiscale sliding MAD is computed on signals bandpassed between 10|$^{-3}$| and 5|$\times 10^{-2}$| 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_T$|⁠) 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 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 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.

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 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 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 |$\sigma$| 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).

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 10|$^{-2}$| and 3|$\times 10^{-2}$| Hz, especially around 2012 November 24 and 29. Some high-frequencies IG-wave component are also more accurately removed above 3.5|$\times 10^{-2}$| 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|$^-2$| and 3|$\times 10^{-2}$| 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 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.

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 station shows 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.

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).

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.

Footnotes

References

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 (⁠|$\lt $| 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):

where k is the wavenumber (rad |$\times$| 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⁻³) and g is the Earth’s gravitational acceleration (m s⁻²). The wavenumber k is linked to the frequency f using the dispersion equation:

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 |$\mu {\rm 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:

where |$\ast$| 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:

where f is the frequency and |$\dagger$| denotes the complex conjugate.

An optimal estimate of |$H(f)$| (written with hat |$\widehat{\cdot }$| ) 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_cp}(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):

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. Theses 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:

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 |$\widehat{H}(f)$| properly. A common procedure is to estimate |$\widehat{H}(f)$| in a reduced spectral range, with the highest level of coherence |$C_{pz_c}(f) \in [0, 1]$|⁠:

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

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 dilated 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.