Temporal and spatial variations in three-dimensional seismic oceanography

. Seismic oceanography is a new cross-discipline between geophysics and oceanography that uses seismic reﬂection data to image and study ocean water columns. Previous works on seismic oceanography were largely limited to two-dimensional seismic data and methods. In this work, we explore temporal and spatial variations embedded in three-dimensional (3D) seismic oceanography. From a 3D multichannel seismic survey acquired for oil and gas exploration in the Gulf of Mexico over six months period, a 3D oceanic seismic volume was derived. The 3D seismic images exhibit both temporal and spatial 5 variations of the ocean, and theoretical and empirical analyses were used to determine their contribution. Results suggest that temporal variation is largely embedded in the crossline direction, showing discontinuities due to acquisition at varied times; however, inline sections, which can be seen as snapshot representation of the water column, capture not only thermohaline structures but also temporal variation of ocean dynamics. Our ﬁndings highlight that marine 3D seismic data can provide superior oceanic images for the investigation of time-evolving mesoscale ocean dynamics.

vided very high-resolution images of the ocean structure, both vertical and, in particular, horizontal. Within the ocean water column, water density variations, fine scale (1-10 m thickness) temperature-salinity contrasts, and water turbidity fluctuations 25 result in small changes in sound speed that produce weak, but distinct, sound reflections. These reflections, 100 to 1000 times weaker than those from the solid earth below, have generally been neglected by geophysicists, whose focus is the structure of the sub-seafloor. However, studies reveal that these reflections correspond to ocean thermohaline structures (Nandi et al., 2004;Nakamura et al., 2006), and are primarily (not completely) associated with temperature gradient (Ruddick et al., 2009). This new cross-discipline, between exploration seismology and physical oceanography, has come to be known as seismic oceanogra-30 phy (Holbrook et al., 2003) and has been successfully applied to imaging mesoscale and sub-mesoscale water-column structures such as ocean fronts (Gorman et al., 2018), eddies (Song et al., 2009;Tang et al., 2014a), internal waves (Holbrook et al., 2009;Tang et al., 2014b;Buffett et al., 2017), and other thermohaline fine structures (Holbrook et al., 2003). Additionally, several theoretical studies have been derived from the application of water-column seismic images including estimation of geostrophic currents (Sheen et al., 2011;Tang et al., 2014b), wave field spectra (Fortin et al., 2016), and internal wave mixing (Dickinson 35 et al., 2017). Most of these studies used 2D seismic data and processing methods, but, in the modern exploration seismology of the oil and gas industry, 2D seismic methods are gradually upgraded to 3D seismic methods, which provide more adequate seismic images with the help of the additional dimension, leading to a more reliable interpretation (Yilmaz, 2001).
3D seismic surveys are distinguished from 2D seismic surveys by the contemporary acquisition of multiple closely spaced lines (e.g., 25 m) that provides regular data point spacing. This leads to a true data volume from which lines, planes, slices, 40 or probes can be extracted in any orientation, with nominally consistent data processing characteristics. While it is possible to acquire a dense, high-resolution grid of 2D lines, such grids are fundamentally different from a native 3D seismic survey (Lonergan and White, 1999;Yilmaz, 2001). The close spacing of 3D seismic lines and the physical illumination of the subsurface from multiple offsets and multiple azimuths remove spatial aliasing problems inherent to 2D seismic data, therefore 3D seismic oceanography has the potential to yield better spatial resolution and can be useful for studying the orientation of ocean 45 structures (Blacic and Holbrook, 2010). To date, 3D seismic oceanography is still not well developed, which can be due to multiple facts: the 3D seismic data are very expensive to acquire, there is an interdisciplinary gap between oceanography and geophysics, and most importantly, a fundamental understanding of the capability of 3D seismic oceanography is missing.
The fundamental principle of exploration seismology assumes prevalently only spatial and not temporal variations in the subsurface (at least in the survey time frame). However, this challenges the development of 3D seismic oceanography, as 50 3D seismic surveys are generally acquired over a long period (e.g., from hours to months) and the water column may evolve during the collection of 3D seismic data. Very recently, 3D seismic oceanography moved a big step forward, successfully being applied to explore the evolution of a front (Gunn et al., 2020) and the passage of an eddy (Dickinson et al., 2020). However, there are still many fundamental questions regarding 3D seismic oceanography waiting to be answered. For example, how are temporal and spatial water-column variations distributed in a single 3D seismic volume? Can we quantify them? Is the oceanic 55 3D seismic imaging meaningful in all directions? The answers to these questions will be useful to guide the development of future 3D seismic oceanography. In this study, we explore temporal and spatial variations embedded in 3D ocean seismic volume and answer the fundamental questions mentioned above. We report results of a 3D seismic oceanography study carried out at the continental slope region of the northern Gulf of Mexico, where seismic imaging of ocean structures can be useful for understanding ocean dynamics 60 and mixing. We analyzed multi-directional (inline, crossline, and depth-slice) seismic images to reveal fundamental imaging features. We further conducted theoretical and empirical analyses to understand the distribution of temporal and spatial variations. This new understanding will be useful for the interpretation of 3D seismic data as spatio-temporal evolution of ocean interior structure, which is important to promote the development of 3D seismic oceanography.

3D Seismic Data
The seismic data used in this study are a portion of a large multichannel 3D seismic survey conducted by Schlumberger West-ernGeco in the Northern Gulf of Mexico for oil and gas exploration between 2002 and 2003. Our study region ( Fig. 1) is at the continental slope of the northern Gulf of Mexico, outside the Mississippi River delta, where the water dynamics in this region is mainly dominated by the Loop Currents and the Mississippi River outflows (Coleman, 1988;Sturges and Lugo-Fernandez, 2005) 70 . Although our seismic region is a bit far north of the Loop Currents, it is still a high eddy-kinetic energy region due to eddies and jets over the continental slope (Hamilton and Lee, 2005). Seismic surveying lines run from northwest to southeast (inline azimuth, 330 • ) over a gentle continental slope (average inclination, 1.5 • ) where the seabed depth ranges from 800 m to 1300 m. Seismic data were collected using 2 airguns of 5085 cubic inches each in flip-flop configuration and 8 streamers, each accommodating 640 hydrophones. The streamer spacing was 100 m, and the hydrophone spacing was 12.5 m. We received the 75 raw shot gathers sampled at 2 milliseconds (ms) and cut at the seafloor with no processing applied apart from analog to digital anti-aliasing filtering (frequency range, 3-200 Hz). We received multiple CTD and XBT casts, including three CTD and two XBT casts inside the seismic area (see Fig. 1).

3D Seismic Data Processing
Imaging the ocean water column using 3D seismic data from the oil industry is challenging. Acquisition geometry is not op-80 timized for this target, and ocean internal reflections are inherently weak and masked by noises of different nature. However, oil industry contractors employ cutting-edge technology and therefore the obtained seismic data are generally of very high quality. The seismic data processing was started with trace editing and amplitude balancing. We have removed offsets greater than 4 km (half of the total channel number) to keep computational costs low, and to reduce the effect of long-offset deviations on the data preprocessing and velocity model building. The strong seafloor reflections were also muted to achieve a relative 85 balance of amplitude levels. The data preprocessing used here employed standard techniques for subsurface seismic imaging, described in Yilmaz (2001), adapted to work with the water-column specific issues. Typical marine noises were present in the water-column data, e.g., acquisition-related noise, direct arrivals, reflection energy returns from previous shots, etc. Environmental noises include wind, shipping activity, and inherent ocean noise. We designed a filtering strategy specifically tailored to remove or minimize as much noise as possible while preserving the weak signals of interest. Consequently, we used low-cut 5 90 Hz and high-cut 150 Hz Butterworth filtering after each processing flow in various data domains to attenuate noise. The direct waves overlapped the subtle internal reflections of the water column and complicated the imaging. We applied frequencywavenumber (FK) filtering and Radon filtering in different data domains to fruitfully attenuate unwanted linear noise, previous shot multiples, and seafloor refraction.
We processed the 3D data in a 3D fashion (which uses both inline and crossline data), rather than sequential 2D processing.

95
The spherical divergence was compensated using the conventional exponential gain function (with a power of 2) after some trials. For velocity analysis and normal move-out (NMO) correction, the data sort was changed to common-midpoint (CMP) and the semblance panels were produced in a window of 1410 to 1590 m/s. To build a reliable velocity model, sound velocities derived from CTD measurements were used to calibrate the seismic velocity model when there is a big mismatch between model and measurement. Afterward, high coherence semblances were picked manually in a distance interval of every 250 100 m in both directions to model root-mean-square (RMS) velocities for further processing. The CTD casts are high-resolution data and hence are helpful for quality control of velocity model building. In this sense, the CTD is used as guidance for the manual picking of the semblance panels. Accordingly stacking and time migration were performed using the RMS velocities to produce interpretable seismic images.
We have also tested and performed a powerful multi-parameter imaging technique known as the common-reflection-surface 105 (CRS) method for stacking and enhancing the data and also to avoid the NMO stretch effect (Yilmaz, 2001). The size of the CRS stacking surface was however chosen conservatively small, in particular in the crossline direction (50 m), to avoid stacking traces from different swaths. The CRS showed improved results compared to the conventional CMP-based imaging especially along the inline direction (Bakhtiari Rad and Macelloni, 2020). Afterward, post-stack time migration was performed using the estimated RMS velocities to correct the dipping events. Finally, the data were time-to-depth converted using the 110 estimated velocity function integrated with in situ sound velocity casts. More details of our seismic processing can be found in Bakhtiari Rad and Macelloni (2020).

3D Seismic Volume
The obtained 3D seismic volume of the ocean interior is illustrated in Fig. 2, showing the "slicing" configuration of 3D inline, crossline, and depth-line seismic images. A 3D seismic volume is generally a 3D matrix in which seismic reflections are color-115 coded according to the reflection amplitude and polarity. Each seismic reflection in seismic images represents an "interface" between two water layers with different acoustic impedances, stemming from different temperature and salinity in the water column (Ruddick et al., 2009). Because a seismic survey can be acquired in any planar directions, the notation in geographic coordinates (i.e., latitude and longitude) is not adopted. Instead, the planar position is organized in inline (the direction; down the continental slope) and crossline (perpendicular to inline; the northeast-southwest direction; along the continental slope on 120 a broad scale); see Figs. 1 and 2. The volume has a spatial resolution of 6.25 m in the inline direction, 25 m in the crossline direction, and 6-7 m in the vertical direction (assuming a sound speed of 1500 m/s and considering the dominant frequency of the airgun signal in the water column between 50 and 60 Hz). Our seismic volume shows only the portion of the water column below 200 m because the geometry of sources and receivers, optimized for deep subsurface targets (several kilometers below the seafloor), poorly image the shallowest portion of the water column (Piété et al., 2013). The signal-to-noise ratio (SNR) of 125 our inline seismic images is ranging from 8-10, using the formula of Holbrook et al. (2013).
It is important to note that for a 3D seismic oceanography study, inline sections are often acquired during short temporal intervals (typically a few hours considering the normal vessel velocity of 2.5 m/s). However, crossline sections, which are obtained interpolating the inline sections perpendicularly, may embrace a wide temporal interval (from a few hours up to several months; the average interval is about 8 hours within the linear acquisition period). We assume that water-column 130 reflectors do not move during the shot and the signal recording (the stick-slip model assumption (Klaeschen et al., 2009)). October.

Imaging Capability Analysis
To test if the 3D volume properly images the water column structures, we compared the seismic reflections with acoustic reflection coefficient derived by inverting the temperature-salinity curves collected along with the CTD profiles (Nakamura et al., 2006). The acoustic reflection coefficient for normal incident sound waves is given by (Kinsler et al., 1999): where ρ 0 , ρ 1 and c 0 , c 1 are density and sound speed; the subscripts 0 and 1 specify adjacent layers. It has been widely demonstrated that the amplitude of seismic reflections is largely related to temperature gradients in the water column (Nandi et al., 2004;Nakamura et al., 2006;Ruddick et al., 2009), but sometimes salinity may play a role up to 40% in regions prone to diffusive convection (Sallarès et al., 2009). suggests that our 3D seismic data successfully represents water-column thermohaline structures. Also, the derived reflection coefficient correlates highly (R = 0.9970) with the temperature gradient [ Fig. 3 (b)], suggesting that the seismic amplitude in our seismic images can be interpreted as the vertical temperature gradient of the water column. In addition, the depth of the 3D 150 volume seafloor reflection matches within 5-10 m error the seafloor depth from multibeam bathymetry which is derived using very high accurate sound speed measurements collected by NOAA's Ship Okeanos Explorer, confirming that our migration is accurate and our seismic images accurately represent the water column structure.

Survey Line Analysis
Once we have established the oceanic natures of the seismic reflections, we resolve the correct relationship between time and 155 space within our 3D seismic volume. We emphasize that inline sections are generally acquired during a short temporal interval been demonstrated that water-column reflectors do not move during the very short time between the shot and the recording of the reflected signal (Klaeschen et al., 2009). In 3D surveys each shot is recorded simultaneously along with parallel streamers, 160 generally referred to as swaths, which for our data set count 8 streamers. To better understand the temporal distribution of the inline/swath, we carefully investigate the acquisition time for each survey line by referring to the field acquisition records. Figure 4 shows the seismic line number, which is equivalent to the center inline number of a seismic swath, as a function of the recording date and time. Considering that 3D surveys are ideally carried out by acquiring adjacent swaths consecutively, we should observe a linear relationship between sailing line number and survey time. However, it is common that owing to 165 operational problems (e.g., weather, maritime traffic, etc.) and/or bad data acquisition some portion of the area may be covered at a different time or some lines may be reacquired. Figure 4 shows a general linear trend, in particular for the portion including Lines 1900-2140. In this interval, spatially contiguous seismic lines were collected consecutively in time, yielding a crossline surveying speed of 1.5 swath per day (equivalent to 0.6 km crossline imaging per day). We refer to this portion as linear acquisition period and we will use it in the correlation analysis discussed further. In this section, we present multi-directional images from our 3D seismic volume to improve the understanding of water-column 3D seismic images and the fundamentals of 3D seismic oceanography, and hopefully, shed a light on possible ocean dynamics that can be resolved from our seismic images. and temporal (about 4 days) separations in these images, we observe that the water-column structures varied, not only spatially (over 5 km), but more importantly, temporally (over 8 days), suggesting a strong water-column mixing process may have 180 occurred in this region. In other words, during the interval between October 20-28, 2002, over this portion (about 4 km) of the continental slope, the strong and thick thermoclines gradually weakened and thinned, and eventually almost disappeared. We interpret the water-column variation as complex and diverse internal wave variation over the continental slope of the northern Gulf of Mexico, caused by a time-varying mesoscale oceanographic process (length > 10 km, depth > 800 m, time > 10 days) which transformed the water column from highly stratified to well-mixed. The possible mesoscale processes causing this 185 transformation will be further discussed in Sec. 4.  Figure 6 shows one of the widest crossline seismic images in our seismic area. We observe significant discontinuities and a strong pattern of seismic swaths (e.g., vertical stripes). These discontinuities stem from abrupt changes of isopycnal depth (e.g., "heaving" of isopycnals) across seismic swaths, indicating the occurrence of temporal variation in the water column. In this region, the ocean dynamics that may change isopycnal depths include eddies, internal waves, internal tides, and Mississippi 190 River outflow. Swaths collected over the linear acquisition period (marked by the red line) display some continuity, while swaths collected with longer time shifts can yield significant discontinuities. For example, seismic imaging of the contrast acquisition period (marked by the blue line) looks drastically different from its neighboring parts due to a 20-day time shift (see Fig. 4).
We also notice that the discontinuity varies with depth. The discontinuity is more significant at 200-350 m [ Fig. 6(b)] than at 600-750 m [ Fig. 6(c)], suggesting that temporal variation is more significant at shallower depths. Another proof for the 195 discontinuity indicating temporal variation is the continuous seafloor, suggesting no temporal variation. Our analysis suggests that temporal variation of the water column, appearing as imaging discontinuity, is significant along the crossline direction, which lowers the quality of crossline images, and may further hinder their interpretation. alize since they offer a new view of the ocean interior. We observe little discontinuity along the inline direction, suggesting that temporal variation is weak along the inline direction. However, we observe significant swath pattern and strong discontinuity appearing along the crossline direction, and the discontinuity decreases from shallow to deep layers. This decreased discontinuity agrees with the fact that temporal variability of the top ocean is more significant than that of the deep ocean, as the ocean is mainly stirred from the top due to the wind-driven surface mixing (Talley et al., 2012). Due to the presence of 205 temporal variation, interpreting seismic depth-slice images is extremely challenging, but we may get a hint on the timescale of the ocean dynamics in Fig. 7, based on the amount of discontinuities and our 3D seismic setup (e.g., time separation between nearby seismic swaths). We suggest that in this region, the ocean dynamics have a timescale less than 8 hours at the depth of 250 m, but way longer than 8 hours at 750 m. The analysis of depth-slice images suggests that in our seismic volume, temporal variation of the water column is non-negligible in the crossline direction but negligible in the inline direction.

210
Our analysis of 3D seismic images reveals that temporal and spatial variations are both embedded in the seismic volume, and temporal variation is mostly along the crossline direction. Next, we will conduct analyses to quantitatively assess temporal and spatial variations within a 3D ocean seismic volume.

Theoretical Analysis
Here is a simple theoretical example to help us understand why temporal variation appears differently in inline and crossline 215 directions. Let us assume a water column in which temporal variation occurs at one cycle per day (1 cpd, about 10 −5 Hz) and spatial variation at one cycle per kilometer (10 −3 cpm). The values here are arbitrary but can represent mesoscale oceanographic processes such as internal waves or internal tides (Talley et al., 2012) for which seismic imaging is particularly useful.
We use single-frequency variations here for easy illustration purposes, while more complex variations can be viewed as the sum of multiple single-frequency variations. For simplification, we will refer to the total variation (sum of temporal and spatial 220 variations) in the inline direction as inline variation, and the total variation in the crossline direction as crossline variation. We assume that this water column is imaged with a common 3D seismic survey setup, having the same field geometry of our survey, and carried out swath-by-swath like the linear acquisition period. The vessel sails at a speed of 2 m/s, the seismic swath is 400 m wide, and the vessel can collect two swaths per day. Therefore, after accounting for the two-way travel time, surveying a 1-km inline section takes 10 3 s (about 0.012 days, or 0.012 temporal cycle), and surveying a 1-km crossline section 225 takes 2.5 days (2.5 temporal cycles, or 2.16 × 10 5 s).
Based on the above settings, the inline and crossline variations can be theoretically resolved [ Fig. 8 (a) and (b)], summing temporal and spatial variations in inline and crossline directions, respectively. Results suggest that in the inline direction [ Fig. 8 (a)], temporal variation is subtle, and the inline variation mostly follows spatial variation. In the crossline direction [ Fig. 8 (b)], however, temporal variation is prominent, and the crossline variation is a fair combination of both temporal and 230 spatial variations, with temporal variation being more noticeable. Furthermore, we find that this conclusion can be extended to temporal and spatial variations for all ocean dynamics. This theoretical analysis confirms that with common 3D seismic survey setups, temporal variation of mesoscale ocean dynamics is more significant in the crossline direction. In terms of seismic imaging, inline images mostly represent spatial variation of the water column, while the crossline images represent a combination of temporal and spatial variations of the water column.

235
Finally, the inline and crossline variations can be spectrally analyzed using Fourier transform, to reveal the distribution of temporal and spatial components in wavenumber domains [Figs. 8 (c) and (d)]. Here, we focused on wavenumber domains instead of frequency domains because seismic images are originally intended to resolve spatial variation rather than temporal variation. We found that the spatial component appears at the same location, e.g., at 10 −3 cpm, in both directions, while the temporal component appears at two different locations, e.g., at about 10 −5 cpm in the inline direction and 2.5×10 −3 cpm in 240 the crossline direction. In other words, in inline images, temporal variation will appear as a low-frequency component, far apart (about two orders, 100 times) from the spatial component, while in crossline images, temporal variation appears as a high-frequency component, very close (about 2.5 times, still in the same order) to the spatial component. The dashed line in the inline spectrum [ Fig. 8 (c)] illustrates the theoretical position of the temporal component, which is too weak to resolve when inline sections are too short (i.e., covering far less than one temporal cycle). These spectral analyses suggests that water-column 245 temporal variation will appear differently in inline and crossline images, posing strong high-frequency aliasing in the crossline imaging. Also, separating temporal and spatial variations in crossline imaging could be very challenging because they are spectrally close to each other.

Empirical Analysis
Our theoretical analysis provides an idea for quantitatively analyzing temporal and spatial variations based on inline and 250 crossline variations. However, real 3D ocean seismic volumes are more complicated than single-frequency-variation examples.
Here we developed a method to analyze temporal and spatial variations within our 3D seismic volume. The results will provide empirical values in our seismic volume, which can also be useful to other seismic oceanography studies with similar setups.
For a real 3D seismic volume, inline and crossline variations can be quantified with correlation, which measures the similarity between two variables or images. For 2D seismic images, x and y, the correlation is defined as: where, i and j form the 2D index, x and y are the means, and σ x and σ y are standard deviations of x and y. In practice, due to the perpendicular relationship between inline and crossline directions, the inline variation can be simply derived from correlation between crossline images [ Fig. 9 (a)], while the crossline variation can be derived from correlation between inline images [ Fig. 9 (b)]. 260 We correlate the inline and crossline image series with corresponding reference images (IL #2140, or XL #18000), respectively. The correlations are plotted as functions of inline or crossline number. The inline axis can further be converted to lateral distances (spatial domain), while the crossline axis can be converted to both lateral distances (spatial domain) and acquisition time (temporal domain). The correlation results are between -1 and 1, representing the fluctuation of image similarity with respect to the reference, with +1 being exactly similar, -1 exactly reversed, and 0 completely different. When the water-column 265 structure varies, either temporally or spatially, the similarity between seismic images varies as well. Therefore, the change of the image similarity summarizes the variation of the water column.
The correlation results allow quantitative analysis of the variation of a particular image series or at a particular depth range.
For example, we only use the linear period (see Fig. 4) for a more meaningful analysis. To best represent the variation of water column, we only use imaging of 250-750 m, excluding seafloors and low-quality shallow portion data. Furthermore, we divide 270 the water column into two parts, the upper part (250-500 m) and the lower part (500-750 m), to explore the depth dependency. Besides, the correlation function allows the examination of the correlation length, which is defined as the length to the first zero crossing, measuring the length required for two images to be completely different. The shorter the correlation length is, the faster the signal de-correlated. Finally, for periodical signals, the spectra of the correlation function can be used to study the inherent cycles within the signal.
275 Figure 9 shows the results of inline and crossline variations derived from our 3D seismic volume. Here, inline variation [ Fig.   9 (c)] is plotted as a function of space (i.e., lateral distance), while crossline variation [ Fig. 9 (d)] is plotted as a function of  variations, respectively. Note that the inline is only associated with space, while the crossline is associated with both time and space. Overall, we observe a similar descending trend in both spectra, providing us information about the spatial features in wavenumber domains. Comparing these two spectra, one can see that many peaks seen in the crossline spectrum are absent 290 in the inline spectrum. These peaks show the dominant frequencies of temporal variations. Particularly, we observe the semidiurnal tidal cycle (marked as M 2 in Fig. 10) in the crossline spectrum, which is strong evidence that the crossline variation includes temporal variations, as continental marginal waters are often heavily influenced by tidal-forcing dynamics. A significant amount of high-frequency components seen in the crossline spectrum are temporal variation, agreeing with our theory [Fig. 8 (d)] which suggests that temporal variation appears as a high-frequency aliasing. Our spectral analysis, again, confirms 295 that the crossline imaging captures temporal features of the ocean, while the inline imaging only captures the spatial features.

Discussion
The visual, theoretical, and data-driven analysis of the 3D seismic volume has provided a wealth of findings that we further discuss and summarize herein. Notice that this is a fairly new field, and we do not have a complete catalog of results and detailed concurrent oceanographic measurements that would allow unambiguous interpretation of the 3D seismic imaging in 300 light of oceanographic processes.
Fundamentally, 3D seismic oceanography resolved from common 3D seismic survey setups will be intrinsically affected by temporal variations of ocean dynamics. Ideally it is possible to design a special 3D seismic survey that can perfectly image an oceanographic process without significant temporal variation. However, 3D seismic surveys are very expensive to collect and, in most of the case, we have to rely on data acquired for other purposes. 3D seismic surveys conducted for oil and gas 305 exploration cover a very large volume of ocean water and have the right bandwidth to depict water column reflections. Those surveys often are acquired over time spanning from days to months, and we must look for the oceanographic signals embedded in the seismic data and build the best processing and analysis workflow to extract them. Mesoscale and sub-mesoscale ocean dynamics (e.g., fronts, eddies, internal waves, turbulence, etc) have frequencies 10 −6 -10 2 Hz and wavenumbers 10 −4 -10 2 cpm (Talley et al., 2012;Ruddick, 2018). As a result, wave spectra resulting from mismatched imaging inevitably produce a 310 significant overlap of temporal variations on top of spatial variation. Temporal variation in 3D imaging may vary case by case, dataset by dataset, volume by volume.
In standard 3D seismic data processing for earth interior imaging, reflections within a large volume of subsurface are utilized and account for out-of-plane acoustic scattering, to provide higher quality images than those resulting from 2D data. The general rule is to avoid summing inconsistent reflections from different seismic swaths. This is not always possible for the 315 water column due to temporal variation. Narrow aperture 3D multiparametric stacking (Bakhtiari Rad and Macelloni, 2020) or swath-by-swath common mid-point 3D stacking (Blacic and Holbrook, 2010) can both be used to create 3D seismic oceanographic images because both implement a time-dependent selection for the seismic data. However, the processing workflow for 3D seismic oceanography can be case-dependent, and we suggest always to perform an in-depth temporal analysis to avoid summing time-inconsistent seismic data.

320
Generally, we found that: Inline images are not affected by temporal variation unless the inline is very long or is imaging fast-moving water-column structures (Klaeschen et al., 2009). Each inline can be seen as a snapshot of the water column vertical thermocline structure at a given time over a given location. This temporal and spatial heterogeneity cannot be fully separated because this temporal evolution occurs at different locations. However, available data analysis methods, e.g. wavefield spectral analysis (Holbrook 325 et al., 2013) and diffusivity estimations (Fortin et al., 2017;Dickinson et al., 2017), developed for 2D seismic oceanography, can be safely extended to 3D inline images. Most importantly, because an inline section "samples" the ocean at one given time and one given place, we could have a particular sequence of inline images that adequately image both the spatial and temporal variations of oceanographic processes. This knowledge can be applied to understand time-evolving mesoscale ocean dynamics, e.g., eddies and internal waves, as the evolution of their vertical structures is very difficult to capture using traditional in-situ oceanographic measurements. Recent 3D seismic oceanography studies with the time-lapsing concept (Dickinson et al., 2020;Gunn et al., 2020) are perfect examples of the application of this temporal variability in oceanographic studies. This temporal feature makes 3D seismic oceanography an extremely powerful tool for studying the temporal evolution of mesoscale ocean dynamics.
Crossline images capture both temporal and spatial information. Temporal variation being more prominent than spatial 335 variation in the crossline direction is validated by the discontinuity in crossline images (Fig. 6), and by faster fluctuation in correlation function analysis in crossline direction (Fig. 9) and is verified by theoretical analysis [Figs. 8 (b),(d)]. In contrast, if temporal variation of an imaged structure is less prominent, it would be buried under spatial variation without showing any discontinuity in crossline imaging. For example, the seafloor appears highly continuous (does not exhibit temporal variations) in our crossline imaging (Fig. 6) because it is associated with a long geological timescale longer than years. In general, ocean 340 crossline images appear as low-quality seismic images, which might limit the imaging interpretation. This is mostly due to discontinuity induced by temporal variation, rather than fewer receivers in the crossline direction than in the inline direction (8 vs. 640). Therefore, the low quality cannot be improved through common processing techniques such as trace interpolation (Yilmaz, 2001). Instead, the crossline imaging can only be improved by optimized 3D seismic survey configurations specifically designed for seismic oceanography, e.g., faster data collection across swaths in the crossline direction. Nevertheless, crossline 345 images still provide valuable information, offering an immediate overview of the temporal distribution of the seismic survey.
In depth-slice images, temporal variation only appears in the crossline direction and is highly depth-dependent. This depth dependency comes from the timescale of ocean dynamics that occur at different depths. Due to water mixing of surface windwaves, the top portion of the water column is usually highly dynamic, while the deep water is associated with a very long timescale (Talley et al., 2012). Such a depth dependence, in terms of imaging discontinuity, can be seen in the crossline images 350 (Fig. 6) and depth-slice images (Fig. 7), which is also evidence of the presence of temporal variation in the crossline direction. This depth dependency in 3D seismic oceanography can be a useful indicator to study the timescale of ocean dynamics at different depths.
A full investigation of the ocean dynamics observed by our seismic images, which would require further analysis of concurrent in-situ measurements and satellite data, is beyond the scope of this study, and will be conducted in our future study. Here 355 we briefly discuss mesoscale ocean processes apparent in our 3D seismic images. Prior studies show the water dynamics in this region are dominated by Loop Currents and Mississippi river outflows, and possible ocean dynamics may include river plumes, internal waves, internal tides, and eddies (Coleman, 1988;Sturges and Lugo-Fernandez, 2005;Hamilton and Lee, 2005). For clear inference from our seismic images (Figs. 5-7) an ocean process should have a temporal cycle longer than 8 days (the time separation between seismic images), a length scale larger than 25 km (exceeding the maximum width of our seismic 360 images), and a depth of influence down to 800 m (generating internal wave field near the seafloor). Cross-referencing the above scales with the typical scales of possible ocean dynamics in this region, the only matched ocean dynamics are eddies. The Loop Currents create omnipresent cyclonic and anticyclonic eddies, moving northward and eastward in the Gulf of Mexico.
We suggest that our seismic images have captured the process of an eddy approaching to the northern continental slope. When eddies approach the continental margin, they start to interact with the continental slope and generate internal wave fields over the slope. This process increases the vertical mixing and reduces the stratification, leading to a change of the water column from highly stratified to well-mixed (as seen in Fig. 5).

Conclusions
This study focuses on providing fundamental understanding of spatial and temporal variations embedded in 3D seismic data collected for oil and gas exploration. A 3D seismic volume of the northern Gulf of Mexico water column has been presented, 370 and temporal and spatial variation in 3D oceanic seismic volume have been investigated for the first time using theoretical analysis and a correlation-based data analysis. Our results shows that in common 3D seismic surveys, temporal variation of mesoscale ocean dynamics can show up in the crossline direction, allowing one to view individual inline images as timelapsing snapshots. The crossline and depth-slice images are heavily affected by temporal variation, and they cannot be directly interpreted. This fundamental understanding is essential for the interpretation and analysis of 3D oceanic seismic images. It 375 allows the extension of previous 2D seismic oceanography analysis, such as wave spectral analysis and diffusivity estimations, to 3D inline images, and to study the evolution of time-varying mesoscale ocean dynamics. The correlation analysis proposed in this study is an effective tool to access the temporal and spatial information within the seismic volume. We suggest that it is important to analyze temporal variation within a 3D seismic volume before the interpretation of 3D seismic images and the analysis of wave fields in the seismic images. We believe 3D oil industry seismic data carry a wealth of oceanographic 380 information, and our future 3D seismic oceanography research will focus on further investigation of the evolution of mesoscale ocean dynamics (eddies and internal waves) in the northern Gulf of Mexico.
Data availability. The seismic data for this project were provided by Schlumberger WesternGeco. The data are available from Schlumberger Multiclient Seismic Data Library (http://www.multiclient.slb.com; area MC-14Q).
Author contributions. ZZ and LZ conducted seismic imaging analysis. PBR and LM performed 3D seismic data processing and imaging.

385
All the authors equally contribute in writing and editing the manuscript.
Competing interests. The authors declare that they have no conflict of interest.