Observations from the seafloor: ultra-low-frequency ambient ocean-bottom nodal seismology at the Amendment field

SUMMARY

Large-scale ocean-bottom node (OBN) arrays of 1000s of multicomponent instruments deployed over 1000s of square kilometres have been used successfully for active-source seismic exploration activities including full waveform inversion (FWI) at exploration frequencies above about 2.0 Hz. The analysis of concurrently recorded lower-frequency ambient wavefield data, though, is only just beginning. A key long-term objective of such ambient wavefield analyses is to exploit the sensitivity of sub-2.0 Hz energy to build long-wavelength initial elastic models, thus facilitating FWI applications. However, doing so requires a more detailed understanding of ambient wavefield information recorded on the seafloor, the types, frequency structure and effective source distribution of recorded surface-wave modes, the near-seafloor elastic model structure, and the sensitivity of recorded wave modes to subsurface model structure. To this end, we present a wavefield analysis of low- and ultra-low-frequency ambient data (defined as <1.0 and <0.1 Hz, respectively) acquired on 2712 OBN stations in the Amendment Phase 1 survey covering 2750 km² of the Gulf of Mexico. After applying ambient data conditioning prior to cross-correlation and seismic cross-coherence interferometry workflows, we demonstrate that the resulting virtual shot gather (VSG) volumes contain evidence for surface-wave and guided P-wave mode propagation between the 0.01 and 1.0 Hz that remains coherent to distances of at least 80 km. Evidence for surface-wave scattering from near-surface salt-body structure between 0.35 and 0.85 Hz is also present in a wide spatial distribution of VSG data. Finally, the interferometric VSG volumes clearly show waveform repetition at 20 s intervals in sub-0.3 Hz surface-wave arrivals, a periodicity consistent with the mean active-source shot interval. This suggests that the dominant contribution of surface-wave energy acquired in this VSG frequency band is likely predominantly related to air-gun excitation rather than by naturally occurring energy sources. Overall, these observations may have important consequences for the early stages of initial model building for elastic FWI analysis.

1 INTRODUCTION

Energy from compressed air guns recorded on towed streamer or ocean-bottom nodal (OBN) arrays is commonly used for active-source seismic data acquisition for exploration activities. Typical active-source processing strategies work to isolate energy of individual shots (or separated from other contributions in simultaneous source acquisition), which are used in subsequent seismic data processing, velocity model building, and migration imaging activities. While air-gun energy sources have long been an industry staple and can provide energy rich in frequencies above roughly 2.0 Hz, they face significant technical limitations in generating sub-2.0 Hz energy at magnitudes sufficient for high-end velocity model building. In particular, full waveform inversion (FWI) requires starting earth models that are sufficiently accurate to enable the simulation of waveforms to within half a wavelength of recorded data. Not satisfying this criterion causes cycle-skipping phenomena that can lead to the FWI optimization processes not converging to the global minimum. This has motivated much research in the development of expanding the lower-frequency bandwidth of energy sources and the acquisition of longer source–receiver offset data, both of which demonstrably improve the stability of FWI analyses (Pérez Solano & Plessix 2023).

A potential alternative source of low-frequency (i.e. sub-1.0 Hz) information is ambient seismic wavefield energy. Naturally occurring energy, generated by swell-induced ocean gravity (and potentially infragravity) waves with typical dominant periods between 1 and 25 s, is known to transfer wavefield energy into the subsurface in the |$10^{-3}-10^{0}$| Hz frequency band (Longuet-Higgins 1950; Webb 1998; Bromirski et al. 2005). This energy propagates predominantly as surface waves (i.e. Scholte, Love and in some circumstances, leaky Rayleigh) at and below the seafloor at velocities controlled by the elastic properties of the solid medium, the acoustic properties of sea water and with frequency-dependent wavefield magnitudes that generally decay away from the fluid-solid interface. Ambient seismic waveforms have been observed throughout the world by, among others, the longstanding Ocean Seismic Network of continuously recording seismometers typically buried in shallow boreholes at 0.1 km depth below the seafloor (Stephen 1998). These high-quality recordings with broadband sensitivity at frequencies between 10⁻³ and 10³ Hz have greatly improved the seismology community’s temporal understanding of microseism phenomena.

The growth in deployments of large-scale OBN arrays consisting of 1000s of recording stations at fixed seafloor locations for up to three months presents an opportunity to greatly improve the spatio-temporal understanding of marine ambient wavefield phenomena. OBN instruments generally consist of four-component (4-C) sensors with a triaxial geophone measuring one vertical and two horizontal components embedded in the solid medium measuring vector particle velocity and a single hydrophone sensor situated in the fluid layer measuring the scalar pressure field. OBN geophones and hydrophone sensors usually have frequency corners between 2 and 15 Hz and thus (are thought to) become decreasingly sensitive to ambient wavefield energy when progressing to decreasingly lower frequencies in the sub-1.0 Hz band. OBN recordings at frequencies lower than the stated geophone and hydrophone corners also are subject to increasing magnitude and phase distortions with decreasing frequency. This fact is commonly assumed to make high-fidelity individual station observations challenging without applying careful instrumentation corrections. Thus, an outstanding question is to what degree are ambient seismic wavefield data recorded on 4-C OBN stations useful for low- and ultra-low-frequency seismic investigation (respectively defined herein as |$\lt 1.0$| and |$\lt 0.1$| Hz)?

Over the past decade, numerous researchers have investigated the OBN array response to low-frequency maritime ambient seismic wavefield energy, usually in the context of seismic interferometry. Olofsson (2010) investigated low-frequency ambient wavefield energy in the 1–10 Hz frequency band recorded for a five-day period on the Astero OBN array located 70 km offshore of the Norwegian coast. [Table 1 presents a list of notable ambient seismic wavefield interferometry applied to ocean-bottom node (OBN) or cable (OBC) arrays.] Bussat & Kugler (2011) applied ambient seismic wavefield interferometry to the same North Sea data set and generated virtual shot gathers (VSGs) that demonstrate the recovery of usable surface-wave waveforms to as low as 0.1 Hz. That work subsequently used the recovered VSG data to constrain the shear-wave velocity structure to 4.0 km depth. de Ridder & Dellinger (2011) demonstrated the use of ambient noise eikonal tomography results for near-seabed imaging at the North Sea Valhall field. de Ridder & Biondi (2015) presented a further case study of ambient seismic noise tomography at the North Sea Ekofisk field. Girard et al. (2023) applied a prestack ambient processing workflow (Girard & Shragge 2020) to data acquired in the 0.3–1.6 Hz frequency band on a Gulf of Mexico (GoM) OBN array. The ensuing interferometric results clearly demonstrated the ability to recover surface-wave information in both Scholte and Love modes between 0.3 and 0.8 Hz as well as waveform sensitivity to large-scale salt structure. The various wave modes measured on 4-C OBNs and their demonstrated sensitivity to complex geologic structure is evidence of the value of using multicomponent data when conducting marine seismic experiments in complex elastic media.

While these studies successfully demonstrate that ambient wavefield data acquired on the ocean bottom can be processed to generate coherent wave propagation across recording arrays and that the resulting waveforms can be used in seismic imaging and inversion investigations, a number of important questions remain about the limits of this style of ambient seismic OBN investigation: (1) can coherent ambient waveforms be recovered by seismic interferometry on dense OBN arrays significantly larger than the typically reported 100–400 km|$^2$| area with sampling sparser than the 4–16 stations per km²? (2) are conventional 4-C OBN instruments capable of recovering usable coherent ambient wavefield information at low- and ultra-low frequencies? (3) does ambient wavefield information extracted at these frequencies on conventional instruments coherently propagate over distances ranging up to 80 km? and (4) what are the implications for long-wavelength elastic model building [e.g. ambient FWI (Sager et al. 2018; de Ridder & Maddison 2018)]?

To examine these questions we use continuous low-frequency ambient wavefield recordings from the Amendment Phase 1 OBN survey, which covers 2750 km² of the Mississippi Canyon and Atwater Valley regions of the GoM. We note that while the large-scale active-source OBN survey was not originally intended for ambient wavefield investigations, the deployment featured a 35-d period when 2750 OBN stations were continuously and simultaneously recording. This extended period of synchronous data acquisition greatly facilitates the extraction and analysis of low-frequency energy, data conditioning prior to cross-correlation and seismic interferometry for estimating VSGs with virtual source–receiver offsets reaching over 80 km in length. Thus, a key objective of this work is presenting the results and observations of applying an ambient pre-correlation data conditioning workflow and seismic interferometry to this unique OBN data set.

The paper begins with an overview and characteristics of the Amendment Phase 1 OBN data set and a description of the ambient pre-correlation data conditioning and cross-coherence seismic interferometry workflows applied to estimate the VSGs. We then present results in terms of frequency decomposed VSGs and illustrate observed surface-wave propagation at the seafloor for the vertical and pressure (Z and P) components including a repeating waveform with a 20 s period that we attribute to air-gun contributions excited with the same periodicity. After investigating the observed wave modes and spatial heterogeneity of observations as a function of absolute offset, we present observations of significant surface-wave scattering from subsurface structure. The paper concludes with a discussion on the potential benefits and inversion opportunities provided by ultra-low-frequency ambient wavefield observations on large-scale OBN arrays.

2 AMENDMENT OBN DATA SET

The Amendment Phase 1 OBN data set was acquired by TGS for an approximately 122-d period in 2019 (Roende et al. 2020). A total of 2750 4-C OBNs were deployed over a roughly 80 km by 40 km area in an staggered grid pattern with 1.0 km mean station spacing in both the inline and crossline directions (see Fig. 1). The survey covers water depths ranging from 0.60 to 2.07 km. The deployed ZXPLR nodal hardware used a chip-scale atomic clock for timing accuracy and 3 Hz hydrophone and 15 Hz triaxial geophone as sensing elements. Sensor hardware settings included 0 and 6 dB preamp gains on the hydrophone and geophones, respectively.

Continuous OBN field records received at Colorado School of Mines from TGS were partitioned into 30-min recordings, with separate files for each of the four components. Because this experiment focused on low-frequency data, continuous waveform data were low-passed with a 4.0 Hz high corner frequency and subsequently subsampled to 0.060 s by TGS personnel prior to being written to disk. The work presented here only analyses the Z- and P-component data; however, we note that the horizontal components are likely useful for complementary identification of different wave modes and phenomena contained within the data set.

After inspecting the recordings from each receiver, we identified a 35-d period when all receivers were concurrently recording. We also discovered that some receivers either were deployed later or ended recording earlier than the vast majority of the survey in that time window. In addition, we noted several nodes that had unidentifiable polarity changes; while these may have been corrected in the active-source survey, it is challenging to identify the corrections needed with ambient records. Due to these uncertainties, we removed the 38 affected OBNs from our analysis, resulting in 2712 OBNs being used for the ensuing experiments.

During the OBN acquisition period an active-source survey was being conducted by three source acquisition vessels, each using two air-gun arrays. There were approximately 2.06 million air-gun shots at nominally 20 s intervals, which were designed to generate frequencies to as low as 1.7 Hz. Because we are predominantly examining energy in the sub-1.0 Hz frequency bands, we did not expect significant overlap from the active sources in these bandpassed and subsampled recordings.

After analysing the ambient data and identifying ambient wavefield characteristics that suggested subsurface influence, we were provided with a P-wave velocity (⁠|$V_P$|⁠) model by TGS. This FWI model was derived from the higher-frequency active-source P-component data using acoustic FWI (Huang et al. 2020). TGS personnel subsequently downsampled the high-resolution model to a uniform 0.1 km spacing in all three dimensions. No velocity information was used in processing the ambient records or calculating the VSGs; however, we use this information for independent corroboration of observed data features.

3 AMBIENT PROCESSING WORKFLOW

This section presents the pre-correlation ambient wavefield data workflow applied to the low-frequency Amendment OBN data set. We note that there are numerous approaches and open-source packages that can be used for ambient wavefield processing (see, e.g. Prieto et al. 2011; Lecocq et al. 2014; Jiang et al. 2020). Here, we follow the approach outlined in Girard & Shragge (2020) that is applied within the open-source Madagascar data processing framework (Fomel et al. 2013). This section highlights the four main workflow steps applied in this study: (1) time-header synchronization; (2) data window selection; (3) time debursting and (4) frequency debursting. The applied workflow has been largely adapted with only minor changes to that presented in Girard & Shragge (2020); readers interested in additional procedural details are referred to this work.

3.1 Time-header synchronization

When using long-time seismic recordings for seismic interferometry, it is imperative that the timing of each sensor is consistent throughout the survey because the interferometry process will otherwise generate incorrect correlation-lag information (and consequently incoherent VSGs). Because acquisition information indicated that the clock errors were less than a single 0.06 s time sample, the raw data were not corrected for clock drift. Therefore, the first step was to ensure that each trace has identical start and stop times to facilitate correlation of traces and recovery of wave fronts propagating across the array with the correct correlation lag for each receiver pair. However, some nodes prematurely stopped recording (due to battery or other mechanical failure) and were therefore removed from the data set in favour of a longer global recording window. Because this affected less than 0.5 per cent of the nodes, we decided that excluding them would not be detrimental to the interferometric analysis.

We ensured that each OBN record was windowed to the same length with the same origin time by examining the header information; fortunately this required no modification from the original records. The data set was organized in 30-min windows upon delivery, a structure that was maintained when generating ‘common-ambient-window’ gathers (i.e. the equivalent of a common shot gather in active-source seismic acquisition). This window duration was deemed sufficient to recover ultra-low frequency information when generating VSGs and was therefore left unchanged during the ambient pre-correlation data processing.

3.2 Data window selection

Not all ambient seismic data are valuable for identifying low-frequency information through seismic interferometry analysis (see, e.g. Prieto et al. 2011). Some records will be dominated by singular large-amplitude events (e.g. earthquakes, OBN deployment ROVs) that dominate the interferometric stacking process. In other cases, individual OBN stations will be influenced by nearby energy sources that are too faint to be detected at other receiver stations and are thus not correlatable across the array and offer no value for interferometric analysis. Therefore, it is important to judiciously select a set of optimal windows with ‘appropriate statistics’ such that calculated VSGs are more likely to highlight weak ambient energy propagating through the earth.

This work approaches data selection of ambient records through a multistep procedure that is due to the non-stationary nature of unwanted signals contained within the data volumes. For example, an impulsive high-amplitude event compared to a moderate-amplitude but longer-duration event can have different effects on the overall interferometric stack quality. One way to address this challenge is to normalize the windows to all have the same energy level (e.g. Draganov et al. 2007), but in this survey there are indeed earthquakes, a hurricane and possibly other high amplitude events that occur during recording time, which can cause unbalanced directional correlations. Here, we use a selection process that aims to eliminate statistically anomalous high-amplitude windows (Nakata et al. 2015). The window-selection step involves removing windows with residual high root-mean-square (rms) energy amplitudes (Issa et al. 2017). To do so, we computed short- versus long-term averages (McEvilly & Majer 1982) to prioritize high-energy windows for removal. Based on this information, we defined a global magnitude threshold (70 per cent) using the pressure component from every OBN (though this could be done using any individual or combination of components) for eliminating windows with abnormally large rms energy values. We do this because of calendar variations in environmental conditions (e.g. effects of severe weather disturbances, distal earthquakes) that cause some windows to exhibit relatively high unwanted signal levels. The resulting recording time used for the remainder of the experiment is 588.0 hr or equally 24.5 d.

3.3 Time debursting

There are different scenarios that can cause individual channels to have a strong ‘burst-like’ energy disturbance on ambient recordings, and including these unwanted energy sources in long-time stacks can skew the statistical convergence of interferometric analyses. Here, we remove burst-like data from individual seismic time-series data using an L|$_1$| iteratively reduced least-squares (IRLS) debursting approach (Claerbout 2014). (We term this ‘time debursting’ to differentiate it from the ‘frequency debursting’ procedure discussed below.) This time-debursting operation addresses residual high-energy events remaining after performing the data masking and window selection steps. To apply this filter, we calculate the envelope for each trace individually and choose a preservation threshold (70 per cent) based on the window rms energy and the highest residual spike amplitude remaining after the data selection step. Waveforms with magnitudes greater than this threshold are reduced to the selected threshold value without hard clipping, while those with lower amplitudes are unaffected by the debursting process. We assert that this data conditioning approach is a judicious alternative to hard clipping or sign-bit normalization, which can introduce frequency-domain artifacts via discontinuous particle accelerations.

3.4 Frequency debursting

Other types of unwanted signals can cause frequency-domain spiking when they are stacked over long periods (e.g. electromechanical signal of repeated turning lights on and off). To address these types of unwanted signals sources, we apply a similar process to the previous processing step to mitigate strong monochromatic (or narrow-band) energy, which manifests as ringing in time-domain VSGs. To do this, we modify the time-debursting method of Claerbout (2014) to operate on Fourier magnitude spectra while leaving the corresponding phase spectra untouched. This filter again down-weights monochromatic energy of significantly greater levels than the background magnitude spectrum to a user-specified level. We then combine the untouched phase and filtered magnitude spectra and apply an inverse Fourier transform to complete the ambient pre-correlation data processing workflow. By removing strong, localized frequency-domain energy we aim to minimize coherent monochromatic noise in processed time-domain ambient data while preserving phase-component information important for reconstructing empirical Green’s function kinematics.

There are a number of different approaches that can be used to address spike-like structure in frequency-domain data. While notch filtering is commonly applied to remove various types of monochromatic noise in active-source seismic experiments (Linville 1994), designing non-stationary notch filters for each window and characteristic frequency makes automation challenging if not impractical. In addition, notch filters can introduce artefacts by affecting phase information. Our approach differs from notch filtering in that it neither requires prior knowledge of the frequency structure nor designs a suite of filters to remove energy at specific peak frequencies. This frequency-debursting technique leaves filtered spectral magnitudes at levels commensurate with those of nearby Fourier components, which is unlikely to be true with notch filtering applications. The parameter defining how much to filter spiky amplitudes was chosen through parameter testing during processing on representative VSG examples.

3.5 Processing QC check

To visualize the effects of the processing sequence detailed above, Fig. 2 presents representative spectrograms taken immediately after the selection process (Fig. 2a) and after applying the full data processing workflow (Fig. 2b) for a station located in the middle of the Amendment OBN array. The first spectrogram exhibits strong vertical banding with otherwise limited energy below 0.4 Hz and between 2600 and 4000 min. Relative to the first spectrogram, the fully processed version now clearly shows significant sub-1.0 Hz energy with coherent energy appearing down to 0.001 Hz, or a period of 1000 s. The first 1300 and final 1400 min of representative recording time also show ‘Dirac combing’ effects that manifest as horizontal lines between 0.05 and 0.20 Hz (see Section 5).

3.6 Cross-coherence interferometry

After data processing, there are several available techniques of interferometry available to reconstruct VSGs from ambient seismic wavefields. Aki (1957) introduces the concept of using autocorrelations to identify wave-mode characteristics to analyse stationary waves in the Earth. Claerbout (1968) arguably develops a precursor to seismic interferometry for a 1-D earth as a method to retrieve a seismic impulse response (Green’s function) by cross-correlating a signal with itself (auto-correlating), thereby creating a ‘virtual’ source for a 1-D medium. Interferometric VSGs also can be generated with an improved cross-correlation-plus-stack workflow that extracts the empirical Green’s function response (Wapenaar 2004). The spectral balance of a VSG can be improved through deconvolution (Wapenaar et al. 2011) or cross-coherence (Nakata et al. 2011) processing, which allows for a choice of smaller regularization parameter and remains stable because amplitude is not explicitly preserved (Prieto et al. 2009).

We calculate VSGs using a CUDA-based code for cross-coherence interferometric calculations (Girard et al. 2023) that first transfers wavefield traces for each window from the CPU to the GPU, computes the forward Fourier transforms over the time axis and then calculates the (symmetric) cross-coherence VSG |$I_{ij}$| contribution between the ith and jth OBN component, |$U_i$| and |$U_j$|⁠, according to:

where |$\overline{U_i}$| represents the complex conjugate of the wavefield |$U_i$|⁠; |${\bf x_A} = (x_A,y_A)$| and |${\bf x_B} = (x_B,y_B)$| are the coordinates of OBN stations A and B; |$\omega$| is angular frequency; m is the window index; M is the total number of windows; the wavefield magnitude is given by |$|U_i| = \sqrt{\left(\Re (U_i)\right)^2 + \left(\Im (U_i)\right)^2}$| and |$\epsilon =0.05$| is a small positive real constant (i.e. after trace normalization) used for spectral-whitening operation (Wapenaar et al. 2011). The sum of the cross-coherence calculation stacks each window after calculating |$I_{ij}({\bf x_A},{\bf x_B},\omega )$|⁠. The GPU code then applies an inverse Fourier transform over the frequency axis to recover the |$I_{ij}({\bf x_A},{\bf x_B},\tau )$|⁠, where variable |$\tau$| is the two-sided temporal correlation lag. Compared to other preprocessing workflow steps, the cross-coherence calculation is extremely fast, taking approximately 12 min on a single NVidia V100 GPU card to compute a single VSG at one virtual-source location for all 1176 windows each with 30 000 samples and 2712 receivers.

4 RESULTS AND OBSERVATIONS

We next present the results of applying the pre-correlation ambient processing and seismic interferometry workflow to generate the output VSG volumes. We first discuss the numerical procedure for generating wavefield propagation images and then analyse the frequency composition and wave propagation embedded within the Z- and P-component VSG volumes. Next, we highlight observations of signal waveforms with an approximately 20 s periodicity asserted to be associated with the same periodicity of active-source air-gun source interval. Finally, we illustrate observations of strong surface-wave scattering from subsurface velocity structure interpreted to be associated with a shallow salt body.

4.1 Gridding and image generation

After interferometric processing, the resulting VSGs were output at each OBN location (see Fig. 1) as a 3-D data cube with the following axes: time lag |$\tau$|⁠, and the virtual shot |${\bf x_A}$| and receiver |${\bf x_B}$| locations. To minimize unused space in the following images, figures appearing below are presented in an inline-crossline ‘deployment’ coordinate system |${\bf x}=[x_i,x_c]$| that is rotated approximately 45° clockwise from geographic north. The first visualization step involved bandpassing individual traces to different frequency bands of interest for the ensuing frequency-decomposition analysis. We then used the OBN geometry to grid the narrow-band 3-D VSG volume for each virtual shot at 0.5 km intervals in both the inline and crossline direction. These intervals were approximately half the mean (staggered) station spacing, which allowed for more accurate spatial nearest-neighbour binning though at the cost of introducing ‘holes’ within the gridded 3-D volumes. For visualization purposes, we next applied uniaxial triangular convolutional smoothing along the inline and crossline directions to infill holes and reduce binning-based high spatial-wavenumber content prior to applying sinc interpolation to spatially map the data to a uniformly sampled 0.25 km grid. While this procedure proved sufficient for figure generation, more advanced spatial data gridding and interpolation techniques is recommended before using the interpolated 3-D VSG data for any follow-on imaging and inversion work.

4.2 Vertical-component VSG analysis

The first set of examples presents a frequency-decomposition analysis of the Z-component VSGs. Fig. 3 shows a representative Z-component VSG time slice extracted at time lag |$\tau =$|10.0 s from an OBN located at inline and crossline coordinates |$[x_i,x_c]=[46.0,16.5]$| km. This volume has been filtered in eight frequency bands within the |$10^{-2}-10^0$| Hz range of interest to highlight different wave phenomena (Shen et al. 2012). Overall, the panels exhibit a remarkable coherency across the illustrated spectral range.

Fig. 3(a) presents the wavefield estimate in the lowest frequency band between 0.008 and 0.04 Hz and shows a near azimuthally symmetric waveforms about the OBN location to a radial distance of about 35 km. Progressing through the next two frequency bands of 0.008–0.075 Hz (Fig. 3b) and 0.01–0.15 Hz (Fig. 3c), the radial expression of the VSG data reduces to respectively about 30 and 25 km. The next three frequency bands of 0.075–0.275 Hz (Fig. 3d), 0.17–0.45 Hz (Fig. 3e) and 0.30–0.65 Hz (Fig. 3f) exhibit increasingly compact waveforms. The final two frequency bands of 0.35–0.83 Hz (Fig. 3g) and 0.50–1.00 Hz (Fig. 3h) continue this trend, though are getting close to the limit where spatial aliasing becomes evident due to the OBN sampling geometry. Overall, the frequency decomposition analysis suggests that the processed long-time VSGs contain coherent information at distances approaching a minimum of 40 km over the 10|$^{-2}-10^0$| Hz frequency range.

Fig. 4 presents an example of wave propagation contained within a VSG volume. The eight wavefield snapshots are extracted from the 0.35–0.825 Hz frequency band shown in Fig. 3(g) and start at 4.0 s (Fig. 4a) and increase with 3.0 s increments to a 25.0 s maximum (Fig. 4h). The wave front expands nearly circularly with increasing lag time; however, various outward ‘kinks’ suggest faster propagation in some locations likely due to the presence of shallow salt bodies.

4.3 Pressure-component VSG analysis

Fig. 5 presents a complementary analysis to that shown in Fig. 3, but now for the P-component data at 10.0 s for a VSG with a virtual source located at |$[x_i,x_c]=[58.0,19.5]$| km. The panels again show eight different narrow frequency bands between 0.008 and 1.0 Hz. Similar to the Z-component data examples, the waveforms appear coherent in all bands with the lowest frequency bands exhibiting a broader expression than at higher frequencies. In addition, the wavefield is starting to appear spatially aliased in Fig. 5(h). Finally, we note that Figs 5(f)–(h) show secondary scattering radiating outward from a point centred at |$[x_i,x_c]=[65.0,16.8]$| km. These observations are discussed further in Section 4.5.

Fig. 6 presents a wave-propagation example in the 0.35–0.85 Hz frequency band extracted from P-component VSG data for the same virtual-source location as in Fig. 5. Figs 6(a)–(h), respectively, present wavefield time slices starting at 4.0 s at 3.0 s increments to a maximum of 25.0 s. We note that as the energy propagates, the wavefield especially toward low inline coordinates becomes increasingly dispersive. Finally, Fig. 6(b) shows the onset of the aforementioned surface-wave scattering that is more evident in Figs 6(c)–(e).

4.4 VSG spatial heterogeneity

The broad radial symmetry of the observed VSGs wavefields suggests that an alternate way to visualize and analyse the observed waveforms is to bin the 3-D VSGs using the following coordinates: correlation lag |$\tau$|⁠, source–receiver absolute offset r and source–receiver azimuth |$\phi$|⁠. One can then compute the average stack over the azimuthal coordinate |$\phi$| to extract a 2-D mean radial |$\tau -r$| VSG panel. The main purpose of the stacking is thus two-fold: (1) improve the overall signal-to-noise ratio of the VSG observations and (2) reduce the data dimensionality from 3-D to 2-D for more effective visual presentation.

Fig. 7 presents representative Z-component radial |$\tau -r$| panels for a VSG from a deep-water OBN located at |$[x_i,x_c]=[3.0,17.0]$| km and 1723 m water depth that have been bandpassed to the same frequency ranges as those presented in Figs 3 and 5. The panels depict coherent energy over almost the full 80 km of absolute offset range. Interestingly, Figs 7(b)–(e) show a clear multiple-like pattern repeating at approximately 20 s intervals. The other panels exhibit weak-to-no repetition, suggesting that the signal band predominantly falls between 0.05 and 0.5 Hz. This observation is further analysed in Section 4.6. A further observation in Figs 7(f)–(h) is that the dominant arrivals start to develop ‘shingling’ behaviour, as indicated by the discontinuous surface-wave arrivals. This behaviour is most prominently observed between 0.4 and 1.0 Hz, and is consistent with Abrahams et al. (2023), which describes the separation of the phase and group velocities for the oceanic Rayleigh wave mode.

Fig. 8 presents the same absolute-offset stack but for P-component observations. Again, propagated energy is visible over nearly the full 80 km of absolute offset. Repetitive signals with a 20 s recurrence interval are visible, though only in Figs 8(b)–(e). There is also a notable change in the character of the observed waveforms between Figs 8(e) and (f) suggesting a transition between different wave-mode types between 0.3 and 0.4 Hz. The waveforms at frequencies higher than this transition range clearly exhibit dispersive behaviour.

Given the bathymetric variations observed in the survey area, an interesting question is whether there are visible differences between VSG volumes for OBNs situated in deep versus shallow water. To address this question, Figs 9 and 10 present results for OBN stations located in deep and shallow water, respectively. Fig. 9(a) presents a phase-velocity-frequency (PVF) plot calculated using Z-component VSG data for a deepwater OBN located at coordinate |$[x_i,x_c]=[3.0,17.0]$| km and at 1.72 km water depth. The image depicts two distinct wave modes with different characteristics. The first mode in frequency bands between 0.15 and 0.45 Hz exhibits a 1.38 km s⁻¹ phase velocity that increases moderately at about 0.08 Hz to 1.40 km s⁻¹. Fig. 9(b) illustrates the Z-component |$r-\tau$| panel filtered to a 0.05–0.28 Hz frequency band, which depicts a strong linear arrival at the expected 1.4 km s⁻¹ moveout. The second mode in Fig. 9(a) falling between 0.45 and 1.25 Hz is significantly more dispersive and exhibits phase velocities in the 2.4–1.7 km s⁻¹ range. Fig. 9(c), which presents the Z-component |$r-\tau$| panel filtered between 0.35 and 0.83 Hz, shows waveforms with a similar 1.4 km s⁻¹ moveout as the interpreted surface waves in Fig. 9(b); however, there is a ‘shingling’ effect comprised of waveforms with shallower dips indicating elements of dispersive wave propagation that are in agreement with PVF-visualized moveouts. Fig. 9(d) presents the P-component PVF plot for the same OBN station and again depicts two types of waveforms. The first is interpreted as a low-frequency surface wave between 0.05 and 0.25 Hz that is significantly lower in magnitude in the P-component relative to Z-component recordings. Unlike in Fig. 9(a), there is no evidence for surface waves in P-component waveforms between 0.25 and 0.45 Hz in this image. Fig. 9(e) presents the P-component |$r-\tau$| panel filtered in the 0.05–0.28 Hz frequency band. The main arrival is very similar to that observed on the Z component (Fig. 9b). The second set of waveforms are interpreted as a guided P-wave package with a dominant lower modes and perhaps as many as three visible higher-order modes. These guided compressional waveforms are generated by multiply reflected/refracted waves trapped within the water column (Hokstad 2004; Shi et al. 2023). Fig. 9(f) presents the P-component |$r-\tau$| panel filtered between 0.35 and 0.83 Hz. The main arrival is distinct from that observed on the Z-component (Fig. 9c) and exhibits dispersive mode interference.

We repeat the analysis presented in Fig. 9 but for data acquired for a virtual source at a shallow-water OBN station located at coordinate |$[{x_i,x_c}]=[79.0,14.0]$| km at 0.96 km water depth (see Fig. 10).

Compared to the PVF plot in Fig. 9(a), the Z-component PVF panel presented in Fig. 10(a) highlights an interpreted surface wave mode largely falling within the same frequency band. The relatively reduced definition may be due to the shorter 45 km distance over which coherent arrival is observed (Fig. 10b) as well as more complex, shallow bathymetry to the northwest of the array. We also note that the dispersive waveforms clearly evident in Fig. 9(c) are now largely absent from Fig. 10(c). The P-component PVF panel in Fig. 10(d) again depicts a low-frequency wave mode at 1.4 km s⁻¹ in the 0.05–0.25 Hz band, though at lower magnitudes and for shorter absolute offsets than that observed in the deep-water example (Fig. 10e). In addition, unlike the deepwater example, only a single interpreted guided P-wave mode is visible in Fig. 10(d).

4.5 Air-gun source contributions

One observation discussed above is the presence of multiple-like energy in VSG data recurring at approximately 20 s intervals in both the vertical and pressure component data (see, e.g. Figs 7b–f and 8b–f). This energy is present in both shallow- and deep-water settings and thus is unlikely associated with bathymetric variations. Thus, an interesting question is what is the cause of these repeating 20 s periodic signals? To investigate this question, we examine the available tabulated active-source shot-timing records. First, we parse out sequential shot-timing data for each of two air-gun arrays on the three boats used to acquired active-source data. We then compute the time difference between shots by differencing the successive shot times. The resulting times are then stacked and binned at 1.0 s intervals to generate the shot-delay histogram presented in Fig. 11. Based on these data, we observe that the dominant shot interval falls between 18 and 22 s with 20 s having the highest recurrence count. Overall, we are not aware of any additional anthropogenic or natural signal occurring with such a dominant 20 s periodicity.

The implication of this observation is that the air-gun array is responsible for generating the interpreted surface waves within the 0.01–0.30 Hz frequency band. Previous research examines the generation of lower-frequency (i.e. 1–5 Hz) surface waves (specifically Scholte waves) excited by explosions or vertical impact sources on or near the seafloor and recorded by ocean-bottom seismometers (e.g. Essen 1980; Schirmer 1980; Rauch 1986; Gimpel 1987; Stoll 1991; Ewing et al. 1992; Stoll & Bautista 1994; Krone 1997). Other work demonstrates that air-gun sources in relatively shallow-water settings can generate surface waves (particularly Scholte waves) of measurable amplitude though at frequencies generally above 2 Hz (e.g. Ritzwoller & Levshin 2002; Klein 2003; Bohlen et al. 2004; Kugler et al. 2007). However, we are unaware of any previous reports of surface waves generated by air-gun arrays in the 0.01–0.30 Hz frequency band. We stress that, for any individual active-source shot record, the generated surface-wave energy would be below the noise floor; however, the consistency of air-gun waveforms over the approximately 2.06 million shots, with the majority of shots falling outside any given interferometric pair, combined with the cross-coherence plus stack processing enables surface-wave energy to become the dominant observable signal in VSGs in this frequency band.

4.6 Surface-wave scattering

A further interesting observation is the surface-wave scattering noted above that persists across a wide range of VSGs and suggests the presence of a point-like diffraction scatterer located in the vicinity of |$[x_{i},x_{c}]=[65.0,16.8]$| km. Fig. 12 presents representative time-slice panels extracted from VSGs in the frequency band (0.58–0.83 Hz) for eight different VSG locations approximately situated at the four cardinal and four intercardinal orientations. The centre panel shows the |$V_P$| model extracted at 2.5 km depth, where the blue and red colours represent sediments and salt, respectively. Each VSG panel shows both portions of the clipped outward-propagating direct surface wave as well as a scattered wavefield emanating from the scattered location. The crosshairs in the panels indicate the approximate centre of the circular scattering. To better illustrate the subsurface feature potentially generating the scattering, Fig. 13 presents a cube view showing the inline and crossline |$V_P$| cross-sections at |$[x_{i},x_{c}]=[65.0,16.8]$| km. This observation suggests that the scattering may be related to a nearby shallow salt ‘pinnacle’ located approximately 1.0 km below the seafloor.

To investigate whether surface-wave scattering in the 0.5–0.8 Hz frequency range might be expected at 1.0–2.0 km depths, we model theoretical 1-D Rayleigh-wave sensitivity kernels using the disba software package (Luu 2021) for the generic sediment |$V_P$| and |$V_S$| profiles shown in Fig. 14(a). Fig. 14(b) presents 0.2–2.0 Hz Rayleigh-wave sensitivity kernels scaled by the period for visualization purposes. The green 0.5 Hz curve shows that the sensitivity peak is at approximately 0.5 km depth and ranges from a null at the seafloor to near zero by 2.5 km depth. Similarly, the orange 1.0 Hz curve peaks at 0.25 km depth and ranges between the null at 0 km and is nearly zero by 1.5 km depth. Thus, the assertion that surface-wave scattering from a salt body located at approximately 1.0 km below the seafloor is not inconsistent with the Rayleigh-wave sensitivity kernels at the given frequency range.

5 DISCUSSION

This section addresses a number of key questions posed in the Introduction: (1) can coherent waveforms be recovered by seismic interferometry on arrays with station spacing of approximately 1 km or larger? (2) can conventional 4-C OBN instruments recover usable ultra-low frequency content? (3) does ultra-low-frequency ambient wavefield energy coherently propagate across larger arrays? and (4) how do these observations affect the prospectus for low-frequency elastic model building?

5.1 VSGs and OBN array sparseness considerations

The sensor spacings used in typical exploration-scale OBN array deployments (e.g. 300–500 m, see Table 1) usually are sufficient for non-aliased spatial sampling of low-frequency wavefields like those observed in the VSG data presented in this experiment. The Amendment Phase 1 OBN survey, though considered to be a sparse array when deployed at 1 km nominal spacing, remains sufficient to recover long-wavelength wavefield information. For ultra-low-frequency data used to create the Amendment VSGs, we still effectively oversample surface waves for wavelengths ranging from ten kilometres to a few tens of kilometres. However, at frequencies nearing 1.0 Hz (and as low as 0.75 Hz for 1.4 km s⁻¹ surface waves), this array nears the Nyquist sampling criteria of two samples per wavelength required for unaliased recovery of the corresponding spatial wavelengths.

For geological scenarios where surface-wave phase velocities are significantly lower than those shown herein (i.e. 1.4 km s⁻¹), one could easily encounter spatial aliasing at frequencies around 1.0 Hz. For example, de Ridder & Biondi (2015) presents an example at the Ekofisk field in the North Sea where observed surface-wave phase velocities fall between 0.4 and 0.6 km s⁻¹ for frequencies ranging between 1.3 and 0.4 Hz. Thus, station sampling around 0.3–0.5 km would be needed to adequately sample the shorter wavelength surface-wave contributions in that particular geological environment.

Overall, we consider OBN spatial sampling may be an important factor when aiming to use VSG wavefield information falling near the spatial Nyquist value for FWI model building activities.

5.2 VSGs and geophone/hydrophone corner frequency

One of the more surprising observations is that the geophone and hydrophone sensors were able to recover coherent wavefield information two-to-three decades below the stated cut-off frequencies—especially considering that no frequency-dependent instrument phase and magnitude corrections were applied as part of this processing. We speculate that there are two key associated factors: the use of interferometric cross-coherence processing and the statistics of long-term stacking. First, we note that the interferometric process uses cross-coherence plus stacking in the mth of M windows of a wavefield recorded at station A with magnitude |$A_m=A_m({\bf x_A},\omega )$| and phase |$\phi ^A_m=\phi ^A_m({\bf x_A},\omega )$| [i.e. |$U_i({\bf x_A},\omega ,m)=A_m{\rm e}^{i\phi ^A_m}$|] with that at station B with magnitude |$B_m=B_m({\bf x_B},\omega )$| and phase |$\phi ^B_m=\phi ^B_m({\bf x_B},\omega ,)$| [i.e. |$U_j({\bf x_B},\omega ,m)=B_m{\rm e}^{i\phi ^B_m}$|]. To investigate the phase component of the interferometric processing, we insert these expressions into the numerator of eq. (1) to obtain:

where the phase difference is due to evaluation of the complex conjugate. Let us now consider a model where the phase component of each signal window is represented by the sum of three elements: (1) the true window-independent wavefield phases |$\psi _A(\omega )$| and |$\psi _B(\omega )$|⁠; (2) the deterministic instrument phase error |$\gamma =\gamma (\omega )$| and (3) random zero-mean noise terms |$\epsilon ^A_m$| and |$\epsilon ^B_m$| usually assumed to arise due to a random Gaussian process. Inserting |$\phi ^A_m=\psi _A+\gamma +\epsilon ^A_m$| and |$\phi ^B_m=\psi _B+\gamma +\epsilon ^B_m$| into eq. (2) yields

which has no explicit dependence on instrument phase error |$\gamma _\phi$|⁠. Moreover, as M approaches ‘large’ (e.g. over a 35-d acquisition period), one assumes that the net contribution of the zero-mean Gaussian error difference |$\epsilon ^B_m-\epsilon ^A_m$| ideally becomes negligible through repeated stacking. Thus, we expect the phase response of calculated VSG data to be sufficiently accurate due to the use of the long-term interferometric cross-coherence-plus-stacking process—even at frequencies much lower than the stated sensing element cut-off values.

We note that a similar analysis on wavefield magnitude spectra is made challenging by the myriad preprocessing steps applied before and during VSG generation. In fact, one of the reasons why we performed narrow-band filtering of VSG volumes in this work is due to the significant variations in magnitude spectra. Specifically, VSG data between |$10^{-3}-10^{-1}$| Hz decades are significantly weaker than those in the |$10^{-1}-10^{0}$| Hz decade. Thus, we stress the importance of performing narrow-band frequency decomposition when analysing broadband ambient VSG energy contributions due to the complexities of handling the variable amplitude scales.

5.3 Ultra-low-frequency ambient wavefield coherence

A similarly notable finding from the VSGs for the Amendment data set is a demonstration of coherent wavefields propagating at ultra-low frequencies with very long associated wavelengths. This is largely due to the much larger aperture of the Amendment array (80 km by 40 km) than those listed in Table 1, which makes it possible to identify propagating surface waves with wavelengths in the tens of kilometre range. In addition to larger aperture, the frequencies recorded in the VSG volumes likely are generated by active-source air-gun excitation with a regular 20 s shooting interval (or equivalently 0.05 Hz). While ocean waves and swell are known to generate energy that transfers into the subsurface in the |$10^{-3}-10^{0}$| Hz band and is generally recognized as a key source of observed low-frequency energy in VSGs, this work presents strong evidence for measured active-source contributions at these low frequencies generated by 2.06 million repeated air-gun excitations.

5.4 Model building prospectus

The effective surface-wave propagation velocity may be considered as a weighted average of the elastic model properties over the depth range where the associated sensitivity kernel exhibits meaningful values (Ekström et al. 2009). At the ultra-low-frequencies shown in Fig. 14(c) (i.e. 0.03–0.20 Hz), Rayleigh waves are sensitive to depths exceeding 10 km. This suggests that the surface-wave modes observed in Amendment VSG data in the 0.05–0.20 Hz range are likely useful for constraining the long-wavelength 3-D elastic model components at the 0–10 km depth range most important for seismic exploration. In addition, any secondary scatterers present in the various frequency bands of observation (see, e.g. Fig. 12) could be used to identify locations of anomalous short-wavelength geological formations (e.g. salt pinnacles) that can be used to further constrain elastic model building analyses.

A corresponding challenge in elastic model building, though, is the need to correctly identify the different types of wave modes present in the data. Initial coupled acoustic-elastic modelling efforts indicate that a wide variety of factors (e.g. source characteristics; observation frequency; bathymetry; presence of guided P-wave modes; background S-wave velocity gradients; and presence or absence of shallow, fast salt canopy of variable thickness) combine to contribute to a large range of possible forward modeling outcomes. While surface-wave-mode (in particular Rayleigh and Scholte) sensitivities are indeed related (Bagheri et al. 2015), they do exhibit distinct characteristics that if incorrectly identified can lead to erroneous inversion-based velocity-model results. A further confirmatory forward modeling effort is currently under way to assist with wave-mode identification; however, this topic remains beyond the scope of this work. The ambient forward modelling framework can also be used to test and validate the ambient wavefield response to model structure including scattering from shallow salt bodies like that asserted in Section 4.6.

6 CONCLUSIONS

This paper presents the results of an ambient wavefield study using low- and ultra-low-frequency data acquired on the large-scale Gulf of Mexico Amendment Phase 1 OBN array. We demonstrate that combining pre-correlation ambient data processing and cross-coherence interferometry workflows leads to the recovery of coherent surface-wave arrivals from as low as 0.008 Hz to about 1.0 Hz. Stacking virtual shot gather (VSG) data over azimuths leads to lag-offset panels that show strong coherency of wavefield arrivals to distances up to (and likely exceeding) 80 km. Phase-velocity-frequency plots suggest the presence of interpreted low-frequency surface wave-mode arrivals below 0.4 Hz in both the Z- and P-component data. We highlight the presence of surface-wave scattering from a shallow salt-body pinnacle that appears in numerous VSGs located at numerous azimuths with respect to the scattering point. Finally, we present evidence that air-gun energy stacked over long periods is measurable on OBN arrays at sub-0.3 Hz frequencies. This assertion is based on the observed 20 s periodicity of waveforms, which is consistent with the mean 20 s active-source shooting interval. This suggests that the dominant generator of ‘ambient’ wavefield energy during the Amendment ambient data acquisition is likely the excitation of active-source air-gun arrays rather than naturally occurring microseism energy. Overall, these findings suggest that ultra-low-frequency seismic energy acquired on standard OBN hardware, after appropriate preprocessing, can generate high-quality, coherent, and interpretable VSG volumes. Moreover, the resulting VSG waveforms show a broad sensitivity to subsurface velocity structure and, thus, may provide a potential pathway forward for generating elastic starting models for FWI analyses.

AUTHOR CONTRIBUTIONS

A. J. Girard (Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Validation, Visualization, Writing – original draft, Writing – review & editing); J. Shragge (Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Validation, Visualization, Writing – original draft, Writing – review & editing); M. Danilouchkine (Data curation, Formal analysis, Investigation, Supervision, Validation, Writing – review & editing); C. Udengaard (Formal analysis, Investigation, Validation, Writing – review & editing) and S. Gerritsen (Formal analysis, Funding acquisition, Investigation, Project administration, Supervision, Writing – review & editing).

ACKNOWLEDGEMENTS

We acknowledge Shell Global Solutions International B.V. for financial support, TGS for the data, and both for the opportunity to publish this work. We also thank our Shell and TGS colleagues for their intellectual contributions over the duration of the project. We also recognize the Center for Wave Phenomena and the Geophysics Department at the Colorado School of Mines, as well as the Wendian HPC cluster and Mines Research Computing for support of this project. Lastly, we thank Deyan Draganov and one other anonymous reviewer for their valuable input during the review process for this paper.

DATA AVAILABILITY

Data associated with this research are confidential and cannot be released.

REFERENCES