Downward continuation of wide-angle seismic data: implications for traveltime tomography uncertainty

SUMMARY

Controlled-source marine seismic experiments are key in advancing our understanding of the Earth’s subsurface structure to study tectonic, magmatic, sedimentary and fluid flow processes. Joint acquisition of wide-angle seismic (WAS) and multi-channel seismic (MCS) streamer data stands as the most robust approach for marine exploration, however effectively mapping subsurface structure remains challenging. The lack of identifiable refractions as first arrivals at short offsets in WAS data creates shallow subsurface illumination gaps up to 6–8 km offsets around Ocean Bottom Seismometers or Hydrophones. This inadequate ray coverage, more pronounced in areas with deeper water column and lower seabed velocities, limits the performance of Travel Time Tomography (TTT) techniques. Velocity determination in the sedimentary layer and reflector location are affected, and errors propagate to deeper layers. This study integrates Downward Continuation (DC) to WAS data. Similarly to our former study where DC is applied to MCS data, redatuming WAS data involves solving the acoustic wave equation backward in time. This process virtually repositions the sources to the seafloor, revealing previously masked near-seafloor refractions as first arrivals. This transformation significantly enhances ray coverage in the shallow subsurface, leading to more accurate determinations of both seismic velocity and reflector geometry. By bridging theoretical concepts with a real data application, this study demonstrates the optimization of field seismic data for improved TTT results. This methodology is particularly beneficial for deep water exploration where spatially coincident WAS and MCS are jointly inverted. In such scenarios, DC-processed WAS data provides the refracted phases key for velocity determinations, and that are typically not present in MCS data due to insufficient streamer length relative to the water column depth. Additionally, we contribute to the community by releasing our open-source, High-Performance Computing software for WAS data redatuming.

1 INTRODUCTION

Estimation of subsurface properties, including P- and S-wave velocities, density and anisotropy, is a crucial issue in geophysics. P-wave velocity (V${}_p$⁠) emerges as a fundamental parameter, aiding in rock type identification and evaluating properties such as porosity, fracturing degree, fluid content and pressure (Prada et al. 2019; Reilly et al. 2019; Arnulf et al. 2021). V${}_p$ information is also fundamental in advanced seismic imaging like pre-stack depth migration, where accurate V${}_p$ models are critical for achieving detailed images of the subsurface. The creation of accurate V${}_p$ models depends on the availability of field data and of appropriate methodologies such as traveltime and waveform tomography. Yet, the full potential of available field data is frequently undervalued.

In marine settings, seismic exploration with controlled source is key to determine the subsurface structure. Seismic velocity is typically obtained through the analysis of the seismic wavefield and the corresponding seismic phases—both reflected and refracted. Marine multichannel seismic (MCS) streamer data provides comparatively high-resolution subsurface information. However, its capacity to retrieve accurate V${}_p$ information is significantly affected by streamer length respect to target depth, particularly restricting the capture of refractions (Jimenez-Tejero et al. 2022). This limitation impacts both ray coverage and the maximum depth of reliable velocity determination. Wide-angle seismic (WAS) data complement MCS records mapping with refracted phases comparatively deeper subsurface structures, mitigating offset limitations to estimate velocity. However, WAS data are typically acquired with a limited amount of Ocean Bottom Seismometers or Hydrophones (OBS/OBH), which in crustal-scale seismic experiments are typically spaced 5–20 km. This receiver spacing results in refracted ray-coverage gaps in the shallow subsurface. In contrast, MCS provides higher coverage with receivers typically spaced as close as 6.25 m. Consequently, combining information (traveltime or waveform) from MCS and WAS data allows to enhance ray coverage at all depth levels, increasing resolution and reducing uncertainties in tomographic results and mitigating ambiguous geological interpretations.

In experiments using only MCS data, the offset and consequently the amount of refractions are constrained by the ratio of streamer length to water depth and the subsurface P-wave velocity structure (Jimenez-Tejero et al. 2022). In such scenarios, the application of redatuming techniques allows transforming the data, so that the availability of refractions can be maximized. Redatuming MCS data involves shifting the source and receiver positions to the seafloor. This transformation enables secondary refractions, initially masked in the original shot gather to emerge as first arrivals in the new virtual gather. Inverting redatumed first arrivals through traveltime tomography leads to an increased ray coverage in the shallow part of the model, thereby providing better constrained velocity models (Gras et al. 2019; Jimenez-Tejero et al. 2022; Prada et al. 2023).

While all redatuming techniques simulate data acquisition at a different surface, they diverge in their wave propagation estimation methods. These approaches fall into two main categories: model-driven and data-driven. Model-driven techniques, such as wave equation datuming (Berryhill 1979, 1984; Shtivelman & Caning 1988; Bevc 1997; Mo 1997), rely on extrapolation operators derived from the wave equation. On the other hand, data-driven methods, like cross-correlation redatuming, involve extracting operators directly from input data sets (Berkhout 1997a, b). Both approaches use various methods, including Kirchhoff summation, plane-wave methods, and finite difference techniques (Gazdag 1978; Berkhout 1981; Schuster et al. 2006; Arnulf et al. 2011; Cho et al. 2016; Gras et al. 2019; Jimenez-Tejero et al. 2022).

Just as with MCS data, challenges regarding masked secondary refractions also arise when dealing with WAS data. However, compared to the more extended use of redatuming in MCS cases, the application to WAS data remains relatively infrequent. Existing studies use redatuming of WAS data, such as wave equation datuming and wavefield extrapolation and decimation for seismic imaging (Barison et al. 2011; Giustiniani et al. 2022) and to extract V${}_s$ information from teleseismic arrivals (Zheng et al. 2019). Planert et al. (2009) employed Kirchhoff summation method following the (Berryhill 1979) approach to redatum WAS data to the seafloor, modifying software initially developed for MCS data redatuming (Arnulf et al. 2011). In their work, the redatumed OBS served as a guide for accurately picking arrival times in original OBS gathers in a marine area with strong fluctuations in reflected and refracted phases caused by water layer thickness and velocity variations. Recently, (Li et al. 2024) showed improved tomographic velocity models by performing the first arrival picking on WAS field data redatumed to the seafloor, also adapting (Arnulf et al. 2011)’s software. Additionally, Li et al. (2024) provides supplementary material testing their redatuming approach, with a focus on grid interpolation to enhance the accuracy of Downward Continuation (DC) results.

Building upon our redatuming approach for MCS data using DC (Jimenez-Tejero et al. 2022), our methodology is applied in this work to WAS data. We solve the acoustic wave equation using the finite difference method to redatume seismic data from the surface to the seafloor. Directly solving the acoustic wave equation offers an alternative to indirect methods and presents its own set of advantages. It provides a more accurate representation of complex wave propagation, enabling improved representation of scenarios with complex bathymetry, varying water velocities and lateral velocity variations. The drawback of solving the wave equation is its computational cost. Finite difference methods have greater computational cost and are slower than indirect methods like Kirchhoff migration. Fortunately, the size and RAM requirements of WAS data from OBS/OBH recordings are manageable with standard computer hardware (compared to MCS data which implies larger data sets), so this is not a limitation. Nevertheless, our code is optimized for High-Performance Computing (HPC) to address the general computational demands.

Our work extends beyond a simple application of DC to WAS data. The primary objective is to reduce velocity and depth uncertainties in tomographic models, particularly within the sedimentary layer overlying the high-velocity basement. To achieve this, we focus on strategic Travel-Time Tomography (TTT) methodologies using joint inversion of MCS and redatumed WAS data. And, as mentioned above, we also provide our open-source user-friendly HPC code to facilitate broader adoption and accelerate research

We first introduce our DC methodology for WAS data, its limitations and a comparison with MCS redatuming, including computational details. Results include synthetic tests and a field data application which shows a comparative TTT analysis with Monte Carlo uncertainty evaluation of combined MCS and WAS data. Additionally, Pre-Stack Depth Migration (PSDM) calculations for MCS-only and joint data are presented and discussed.

2 MATERIALS AND METHODS: DC APPLIED TO WAS DATA

The redatuming of WAS data via DC offers a straightforward implementation, encompassing all data without the necessity for prior phase selection. It is adaptable to irregular new datum surfaces, remains stable in steeply dipping areas, and can accommodate arbitrary V${}_p$ models. Despite these advantages, limitations are associated primarily with the finite-difference scheme and the model-dependent nature of redatuming. However, insights from DC applied to MCS data (Jimenez-Tejero et al. 2022) show that the lack of detailed velocity model information does not compromise reliability. Since propagation occurs within the water column with negligible velocity variations, employing a homogeneous yet realistic water velocity value consistently results in errors of less than 2 per cent, compared to using water velocity models derived from field data.

In this section, we first outline the main components of our DC method, which involves solving the two-way acoustic wave equation in a 2-D domain, implemented backward in time. We use benchmark synthetic homogeneous V${}_p$ models to demonstrate the impact of water column depth, as well as the near-seafloor velocity model, on the registration of refractions in WAS field data experiments. Furthermore, we discuss the main computational aspects of the DC-WAS code that we provide, and describe the data pre-processing techniques employed to enhance the quality of the results. Lastly, we describe the main aspects of the traveltime tomography method used in the analysis.

2.1 Main ingredient: two-dimensional acoustic wave equation

The tool for redatuming seismic data to the seafloor is based on the 2-D acoustic wave equation through the water layer, incorporating homogeneous, isotropic and acoustic approximations. Solving the wave propagation equation entails a finite difference scheme of sixth order in space and second order in time. To prevent artificial numerical reflections, a free surface is applied at the top, coupled with a Convolutional Perfectly Matched Layers scheme (Pasalic & McGarry 2010) at the left, right and bottom boundaries. The wave equation is defined differently outside and inside the PML layers. The acoustic wavefield propagation, represented as $p(r,t)$⁠, for a given source, $f_s(S,t)$⁠, within a velocity model, $v(r,t)$⁠, is expressed as follows:

The additional term for the PML layer, denoted as $f_p^{\text{PML}}$ is characterized with the auxiliary variables ${\psi }_k(p)$ and ${\xi }_k(p)$⁠, governing the evolution for each component (k = x, z) and time step, n, as follows:

where the parameters a${}_k$ and b${}_k$ are:

with $\mathrm{\Delta }t$ being the time step, ${\sigma }_k$ the absorption damping factor of the acoustic wave and ${\alpha }_k$ the real positive pole shifting factor (Zhang & Shen 2010).

2.2 Limitations and considerations of DC in WAS versus MCS data

Similar to the application of DC to MCS data (Jimenez-Tejero et al. 2022), DC for WAS data involves back-propagating seismic data to the seafloor through the water model, using the acoustic wave equation. For MCS data, the method requires two steps: first, redatuming receivers, and then redatuming sources as point gathers. The redatuming of WAS data is performed in a single step. For each OBS/OBH gather, the sources are propagated as a supershot backward in time, ensuring accurate waveform transformation. In this process, while the positions of the sources within each OBS/OBH gather maintain their original horizontal coordinates, their vertical locations are adjusted to conform to the seafloor bathymetry.

DC results are influenced by the underlying methodology and limitations arising from physical and experimental factors. The qualitative findings from (Jimenez-Tejero et al. 2022), extend to applying DC to WAS data. Regardless of this, we will now revisit some of the key DC methodology factors for WAS data:

Refraction content: With respect to reveal the missing refraction offset in the original data using the DC method, the length of recorded refractions is mainly governed by two key parameters: the seafloor depth and the velocity structure of the subsurface. The clear advantage of WAS data over MCS is this respect, is that the acquisition offset is not limited by the streamer length. And therefore it is true to say that applying DC to MCS is not an option in those MCS experiments where the water column depth versus streamer depth constrained the recording of refractions. However for a WAS experiment, it is always possible to recover the near seafloor refractions by applying DC to each of the OBS/OBH. In the case of WAS data, the minimum offset with visible refractions as first arrivals in a OBS/OBH gather at any marine setting is:

where D${}_{\mathrm{s}}$ and D${}_{\mathrm{r}}$ are the source and receiver depths, the angle $\varphi$ is the seafloor slope, $w=\sqrt{1-u^2}$ and $u=V_{\text{water}}/V_{\text{Earth}}$⁠. We represent eq. (7) for a flat seafloor (⁠$\varphi$=0) in Fig. 1, exploring the parameter space using benchmark models. All models maintain a constant water velocity of $V_{\text{water}}$ = 1.5 km s⁻¹ and consist of a single homogeneous subsurface layer. The y-axis indicates variations in water column depth, while the x-axis corresponds to the subsurface V${}_p$ layer. Consistent with findings from the mentioned previous study on MCS data, we observe a decrease in refractions inversely proportional to increasing water depth. Moreover, higher sub-seafloor velocities result in refractions recorded as first arrivals appearing at shorter offsets.

For real-case scenarios with non-homogeneous velocity models, this formula can still provide a useful estimate of the offset range where refracted arrivals are likely to appear in WAS data, by using an average velocity for the near-seafloor region.

Water model dependence: Water model velocity is the primary physical parameter in wave propagation. Our MCS study demonstrated that a constant velocity model outside the true range can introduce errors up to 2 per cent over water depths up to 6 km. Therefore, while our code supports XBT (Expendable bathythermographs) data input, a realistic constant water velocity is sufficient for accurate results.

Data quality and signal to noise ratio (S/N): DC results are inherently noisy due to back propagation. Optimal seismic data for DC requires an acquisition geometry with adequate grid coverage, resulting in DC results with an S/N that allows for good picking performance. Therefore, seismic data with poor acquisition geometry may require interpolation before applying DC. In this regard, the interpolation tests shown for improving the S/N ratio and the quality of the data in DC results for MCS data are equally valid for DC on WAS data (Jimenez-Tejero et al. 2022). For WAS data with insufficient intershot distances, interpolating the original data onto a denser shot grid would enhance the S/N ratio in DC results. Also, reducing the time step is a potential solution (Li et al. 2024). In our current study interpolation was not performed because the S/N ratio in the DC results on our WAS field data case allowed for reliable first-arrival picking.

Data pre-procesing: Data filtering before redatuming is essential for maximizing the quality of the results. Here, we list some useful processing techniques. OBS/OBH gathers, used as input for this software, require two types of muting to avoid noise. It is recommended to truncate the recording time just above the first multiple—the primary reflection of the seismic wave at the sea surface. In an OBS/OBH gather for the zero-offset trace, the cut-off at the first multiple is located at 3D${}_{\text{water}}$/v${}_{\text{water}}$⁠, where D${}_{\text{water}}$ refers to the depth of the water column layer. Another type of muting facilitates obtaining cleaner results by eliminating information from the trace before the first arrival. Both muting procedures can be implemented through various forms of filtering or surgical muting. Also, filters such as the Butterworth type can be applied to eliminate frequencies to amplify the signal-to-noise ratio. In our field data case, we apply Butterworth below 2 Hz and above 80 Hz.

2.3 Computational aspects

We present a new software tool for redatuming WAS data using the DC method, written in modern Fortran and tailored for HPC environments using Open MPI for efficient parallelization. This section details the computational features of the code, inheriting many aspects from the equivalent software designed for MCS data redatuming (Jimenez-Tejero et al. 2022).

HPC architecture and computational resources: MPI parallelization accounts for the number of OBS/OBH to be redatumed. Each set of shots recorded at each OBS/OBH gather is downward continued to the seafloor and propagated as a super shot through the water layer back in time. This process is independent for each OBS/OBH gather, so the number of OBS/OBH gathers is parallelized with the available process slots, that is, CPUs or cores depending on the operating system. The most efficient performance in time consumption is achieved when the code is parallelized using the same number of processes as OBS/OBH gathers. In this configuration, the computing time ranges from a few seconds to a few minutes. In terms of memory space, RAM usage is determined by the grid and water model size, which covers the shotgather length and water depth at each OBS/OBH. Consequently, RAM requirements are typically low, as OBS/OBH recordings do not necessitate high-resolution grids.

WAS data format: The WAS input field data is directly read in SU format, a modification of the standard binary SEG-Y format designed to work with Seismic Unix software (Cohen & Stockwell). Parameters like the common receiver gather location and the offset at each trace can be provided through the headers. The bathymetry is provided as a text file, and the values of the main acquisition parameters (number of receivers, nt, dt, etc) are read from a text file. The OBS/OBH gather redatumed results are also stored in SU format. For a more detailed information, we refer to the software manual for the DC-WAS code.

3 RESULTS

This paper focuses on the functionality and key role of DC-WAS data in TTT, as well as strategies for optimal TTT velocity model building using combined MCS and WAS data. The velocity models are calculated using a modified version of the code tomo2d (Begovic et al. 2017), originally developed by (Korenaga et al. 2000), to perform the joint refraction and reflection TTT (Melendez et al. 2015a, b; Merino et al. 2021). This modified version integrates source and receiver geometry of streamer MCS and WAS data for their joint inversion. The tomographic method computes synthetic traveltimes with a ray tracing algorithm based on the graph and ray-bending refinement methods (Moser 1991; Moser et al. 1992). The regularized inversion is performed using a least-square sparse matrix solver, updating the initial model iteratively to minimize the misfit between observed and synthetic traveltimes. This iterative process continues until a given Chi2 criterion, where the misfit is within the range of traveltime picking uncertainty, is satisfied.

3.1 Synthetic data example

We show as an example, synthetic tests to illustrate the influence of water depth on first-arrival traveltimes (direct rays and refractions) recorded at OBS/OBH stations, and how DC recovers hidden sub-seafloor first-arrival traveltimes that are otherwise screened in the original data. These tests involve the generation of synthetic OBS/OBH Vp records simulated using the acoustic wave equation with a 10 Hz central frequency and a 50-m shot interval. A 1-D gradient velocity model, representative of the Porcupine Basin and constrained by streamer tomography (M. Prada et al. 2019), was used in these simulations. This model featured a velocity gradient of approximately $0.8{\text{s}}^{-1}$⁠, with P-wave velocity (⁠$V_p$⁠) increasing from approximately 1.5–1.7 km s⁻¹ at the surface to 4 km s⁻¹ at the base of the sedimentary layer. While acknowledging the limitations of using the acoustic approximation and a gradient velocity model for the simulation of a synthetic gather, this approach serves the purpose to qualitatively demonstrate the impact of DC on the wavefield and its ability to recover missing near-offset refractions.

The example presented in Fig. 2 compares synthetic first arrival traveltimes of the synthetic OBS/OBH gather generated for two different water depths and the same 1-D gradient as velocity model. Fig. 2(a) illustrates the ray path obtained through the shallow-water case with a seafloor depth of 0.5 km, while Fig. 2(b) represents the deep case with a water depth of 5 km. In Fig. 2(c), both synthetic first arrival traveltimes from each scenario are compared. Notably, in the shallow-water case, the direct water wave extents to almost 3 km of offset, while in the deep-water case, it extends to around 9 km offset at each side of the OBS/OBH location.

Fig. 3 illustrates the impact of DC transformation of the OBS/OBH gather on the 5 km deep-water scenario. The DC process effectively recovers the missing near-seafloor refractions as first arrivals (panel b). Ray tracing is also shown for the OBS/OBH gather (panel c) and the ray coverage improvement after DC transformation (panel d).

3.2 Field data case and Monte Carlo analysis

We present a field data case using the 100-km-long line ATLANTIS-I, acquired during the ATLANTIS survey in 2022, southwest of Madeira Island, within the Atlantic Ocean at  31${}^{\circ }$N (Fig. 4). The data set comprises MCS data acquired with a 3-km-long streamer and WAS data acquired with 11 OBH spaced  7–8 km apart and a 100 m shot interval. Fig. 5 presents the OBH number 4 from the seismic line (Fig. 5a) with its downwards continuation to the seafloor (Fig. 5c). Figs 5(b) and (d), provide a closer zoom, highlighting the first arrival picks. Pre-processing for cleaner results for redatuming of the field data includes muting the signal up to the first arrival and from the first multiple. We do not interpolate data to thinner shot grids because the signal-to-noise ratio in the DC results is sufficient for accurate first-arrival picking. We refer to Fig. S1 in the Supporting Information, which shows the first seven OBH gathers, which are the ones redatumed for this study.

Our primary objective is to employ a joint strategy, combining MCS and WAS data, to estimate the V${}_p$ model using TTT. When considering the implementation of the joint strategy for redatuming purposes, it must be noted that the presence of refractions in MCS data is limited by the streamer length, water depth and subseafloor velocity. In this experiment, the 3-km-long streamer is too short to contain refracted arrivals given the $>$4500 m water depth (Jimenez-Tejero et al. 2022). Thus, redatuming the OBH data offers a solution to increase ray coverage, leading to improved accuracy and resolution of the V${}_p$ model in the upper 2–3 km.

The results are presented in Figs 6–8. To evaluate the use of redatumed WAS data in TTT, we invert for the sediment layer velocity structure and the geometry of the top of the basement (TOB) with three different strategies. Initially, we explore a strategy relying solely on traveltimes from MCS data, used in experiments when no OBH data are available (Fig. 6). In the second test, we implement the MCS + WAS joint strategy using the original streamer and OBH records (Fig. 7). Finally, in the third test, we extend the joint strategy by inverting WAS data redatumed to the seafloor, defined as MCS + WAS(DC) (Fig. 8).

In each of the three cases, we performed TTT following a Monte Carlo approach to compare the effectiveness of each case in reducing model parameter uncertainty (i.e. velocity and geometry of reflector). The Monte Carlo statistical analysis involves the inversion of 500 realizations. Each of them combines a randomly generated 1D velocity profile as starting model, a set of traveltimes with Gaussian noise added randomly, based on manual picking uncertainty (from 20 to 40 ms), and a randomly generated initial reflector for the TOB. Each initial velocity model was generated from randomly varying 10 per cent velocities of the reference model (Fig. S2a, Supporting Information). The reference reflector was derived by smoothing the depth-converted two-way travel times from MCS data for the top of the basement, assuming an average Vp of  2 km s⁻¹. We use the same ensemble of starting models and reflectors for the three cases. After inverting the 500 realizations, in each case, we compute the improvement factor, (⁠${\sigma }_0$-${\sigma }_f$⁠)/${\sigma }_0$⁠), which quantifies the reduction in the initial standard deviation of V${}_p$⁠. With respect to the different panels shown in Figs 6–8, the panel (a) presents the results of the average V${}_p$ velocity model of the 500 realizations, while panel (b) illustrates the V${}_p$ difference between the initial and final (average) models and panel (c) shows the improvement factor. We refer to the Supporting Information to visualize the initial V${}_p$ model and standard deviation (Fig. S2, Supporting Information).

We clarify that in Fig. 8, the redatumed OBH are receivers 1 to 7. OBH gathers from 8 to 11 were not redatumed because they already present a sufficient number of refractions because sediment thickness is comparatively larger than that beneath OBH 1 to 7 Fig. 4. This results in low standard deviations, as shown in Fig. 7(c). However, it is important to note that all OBH gathers in the same line could be systematically redatumed if desired.

4 DISCUSSION

Before analysing deeper into the results presented in Figs 6–8, we will first discuss the tools that helped us efficiently assess the validity and relative improvement of the inversion results. Table 1 provides the Traveltime misfit values for both reflection and refraction data. For the three strategies, MCS, MCS + WAS and MCS + WAS(DC), the decreasing misfit indicates that the inversion has converged to a minimum, but it does not necessarily guarantee a global minimum. While the final misfit achieves a ‘low enough’ value in any of the cases, a more reliable assessment of the inversion is obtained in this study through the Monte Carlo analysis, which provides insights into the uncertainty of the results, independent of the initial model.

As shown in Fig. 6, in the case of MCS-only inversions, the minimal reduction in initial standard deviation indicates very limited improvement, even with a decreasing misfit. Adding the WAS data (Fig. 7), especially when redatumed (Fig. 8), leads to a significant increase in the reduction of initial standard deviation in areas with significant refraction information, which indicates that the inversion process has effectively utilized the additional refraction given. These areas correspond to the raw OBHs from 8 to 11 in Fig. 7 and the redatumed OBHs from 1 to 7 in Fig. 8, respectively. In conclusion, while misfit reduction is a necessary condition for a successful inversion, the V${}_p$ difference and the Initial Standard deviation reduction between final and initial averaged models following the Monte Carlo analysis presented, is crucial for evaluating the robustness of the results. Most important aspect to clarify is that all three Monte Carlo tests (MCS, MCS + WAS, MCS + DC(WAS)) use the same ensemble of initial models and starting reflectors. Thus differences in each Monte Carlo test result primarily from differences in the set of traveltimes being inverted

4.1 Discussion of results

Results show that TTT using only reflections on MCS data provides an apparent well constrained reflector inversion (Fig. 6). However, Fig. 6(b) shows that the V${}_p$ model is similar to the initial model and that there is little reduction in standard deviation after the 500 realizations (Fig. 6c), indicating that the real V${}_p$ is poorly recovered. Even though the geometry of the reflector seems to converge towards a common solution, the final solution departs little from the initial models, with only a marginal reduction of the initial V${}_p$ standard deviation.

In contrast, when employing the joint MCS + WAS strategy, the outcome improves (Fig. 7). Notably, from OBH number 4 to number 11, there is an increase in V${}_p$ differences, and an improvement in the initial deviation occurs from OBH number 8 to 11 (Figs 7b and c). This improvement occurs in areas with thicker sedimentary basins (i.e. beneath OBH 8 to 11 Fig. 4), where the amount of refracted traveltimes through the sedimentary layer are larger. In contrast, from OBH 1 to 8, the initial deviation barely improves due to a trade-off between the number of reflections from MCS and the limited number of refractions from WAS data, rendering the V${}_p$ model unreliable (Fig. 7c). Additionally, the joint MCS + WAS inversion exhibits increased uncertainty in the location of the inverted reflector compared to the scenario using only MCS reflections. This result shows that the limited sedimentary refracted P-wave traveltimes from OBH with no DC are not sufficient to converge into a common tomographic solution in terms of velocity structure. This emphasizes the need for increasing refracted traveltime information to outweighs the tomographic information provided by reflected traveltimes that suffer from large velocity-depth trade-off.

When applying the joint strategy with WAS data redatumed to the seafloor, MCS + WAS(DC), the comparison of the initial and final models (Fig. 8b) reveals consistent reduction in the initial standard deviation (Fig. 8c) across the entire model. This suggests that the inversion converges towards a single, improved solution regardless of the initial model chosen.

Fig. 9 compares the results obtained using the three different inversion strategies at kilometers 20, 50 and 80 along the seismic line. At kilometers 20 and 50, the reflector depth is similar, while at kilometer 80 it is significantly deeper. Areas with a deeper TOB reflector naturally exhibit more refractions in the raw OBH data (Fig. S1 in the Supporting Information). Consequently, the joint MCS + WAS strategy performs reasonably well at these locations (Fig. 7). However, the joint strategy with MCS + WAS(DC) further improves the results by reducing uncertainty (Fig. 8). Overall, incorporating first arrivals through the DC redatuming strategy leads to consistently lower parameter uncertainty throughout the model. The V${}_p$ model uncertainty remains within a range of 100–150 m s⁻¹, and the reflector location uncertainty is constrained to 100 m.

Robust resolution of the velocity structure of the upper few subseafloor km is crucial for understanding processes ranging from oceanic crust formation (e.g. Arnulf et al. 2011; Grevemeyer et al. 2018; Peirce & Hobbs 2024), or the interplay between faulting, lithology and fluids (e.g. Olsen et al. 2020; Arnulf et al. 2021). Most of these settings are located in deep water environments. Our results illustrate the potential to increase the data for available resolving velocities with limited traveltime information, particularly in these deep water regions. The inversion of limited refracted traveltimes may provide velocities and geometries of reflectors biased by the initial assumption (Fig. 9). However, the implementation of DC and joint inversion strategies may mitigate this issue, reducing velocity and geometry uncertainty inherent to TTT, facilitating a more robust assessment of the shallow V${}_p$ structure, even if the receiver spacing is large.

The robust constraint provided by shallow velocity structures is essential not only for studies investigating shallow layers but also for deeper crustal-scale studies including the petrological nature of the basement of the crust (e.g. Grevemeyer et al. 2018; Merino et al. 2021) or the tectonic structure of subduction zones (Sallares et al. 2021). The joint TTT strategy presented here helps minimize uncertainties related to shallow vertical velocity gradients and the position of the TOB (Figs 9 and 10). This, in turn, should prevent the propagation of errors to deeper geological layers. This allows for an improved characterization of vertical velocity gradient of the crystalline crust, reducing uncertainties in the petrological interpretation and the location of deeper interfaces (crust–mantle boundary or subduction interplate plate boundary).

4.2 Applications beyond TTT velocity models

Pre-Stack depth migration (PSDM): While DC redatuming offers valuable advantages in jointly inverting MCS and WAS(DC) data to produce more accurate velocity models with reduced uncertainty in reflector location, caution is necessary when applying these models to tasks like PSDM calculation. A better tomographic velocity model does not always translate into a better PSDM result, especially in the presence of significant anisotropy in the sediments. Seismic waves travel at different velocities depending on the propagation direction, with horizontal waves typically exhibiting faster velocities. This results in slower apparent velocities being registered in streamer data from MCS surveys compared to WAS data, where each OBH captures rays travelling longer horizontal distances. In our field data case, the presence of significant anisotropy and the short length (3 km) of the MCS streamer prevent the velocity model obtained with DC-enhanced joint data from adequately collapsing the shallow part of the sediments in the PSDM calculation.

Fig. 11 illustrates this issue. Fig. 11(a) shows the difference between velocity models from MCS + WAS(DC) joint inversion and MCS data alone in the shallow section, highlighting the significant anisotropy, approximately 18 per cent difference between horizontal and vertical velocities. This aligns with previous studies, such as (Sallares et al. 2013), which observed a 15 per cent anisotropy at the Nicaragua convergent margin, and (Carcione et al. 2012), which reported anisotropy variations up to 30 per cent between perpendicular and parallel incidence. The MCS + WAS(DC) model exhibits faster Vp than the MCS model, as expected for marine sedimentary basins formed by subhorizontal layers. The PSDM calculations for both velocity models shown in panels (b) and (c) were performed using the raytracing and Kirchhoff Depth Migration algorithms from the Seismic Unix software (Cohen & Stockwell).

Fig. 11(c) demonstrates that the MCS + WAS(DC) velocity model, when used for PSDM processing, produces migration artifacts within the sedimentary section, characterized by a disrupted collapse of seismic reflections and a resulting discontinuous TOB horizon. To address this issue, a combined model approach is recommended (Bartolome et al. 2005). This approach utilizes the velocity model derived from MCS data for the shallow section until the TOB and the tomographic model obtained from sub-basement OBH traveltimes for the deeper section (Fig. 11b).

Full waveform inversion (FWI): The benefits of DC redatuming extend beyond improving TTT imaging. In high-resolution techniques like FWI, where success depends on the quality of the initial model to avoid nonlinear challenges like cycle-skipping, incorporating DC within the TTT framework can further enhance these TTT-derived initial models. This significantly improves the probability of FWI converging towards an accurate, high-resolution representation of the subsurface (Gras et al. 2019).

5 CONCLUSION

This work focus into the challenge of extracting valuable near-seafloor refraction information often hidden within WAS data, specially in deepwater marine environments. We provide a user-friendly DC software for processing WAS data (Jimenez-Tejero et al. 2024), adaptable for both local systems and high-performance computing clusters.

• The core of the approach lies in the DC algorithm by solving the 2-D acoustic wave equation to model complex wave phenomena as accurately as possible within the limitations of the finite-difference scheme. This direct solution, although computationally more intense, offers potentially greater accuracy over indirect methods like Kirchhoff migration, particularly in complex environments characterized by rapid lateral velocity variations, complex bathymetry, or strong velocity contrasts.

• The DC is applied to each OBS/OBH record leading to the relocation of sources to the seafloor. This process unmasks hidden secondary refractions, enhancing their identification as first arrivals in the new virtual configuration.

• Our analysis highlights the power of joint TTT with both reflections from MCS streamer data and refractions from OBS/OBH data, including arrivals revealed through DC. This approach better resolve the velocity structure of the shallow subsurface, reducing uncertainty in the sedimentary layer and consequently the position of the TOB. This mitigates the propagation of errors to deeper layers, leading to more accurate subsurface imaging.

• Building on our prior research on DC application to MCS data, this work introduces an additional effective solution by applying DC to WAS data. Through both studies and the corresponding DC softwares, we address the critical aspect of uncovering and identifying refractions regardless of the seismic dataset type. Notably, when MCS data is unavailable or unsuitable for redatuming, the DC approach remains a valid tool for processing WAS data alone or jointly with MCS, as shown in this work.

• Compared to regular TTT velocity models, DC-enhanced models provide a significant boost to the robustness of initial models used in FWI to achieve accurate high-resolution results. However, caution is necessary if the goal is applying these velocity models (from either TTT or FWI) for PSDM calculation. In scenarios with significant anisotropy and short streamer length relative to water depth, a combined model using MCS data for the shallow part of the model until the TOB and joint data for the deeper section is recommended for optimal PSDM results.

ACKNOWLEDGMENTS

This is a contribution of the “Grup de Recerca Consolidat de la Generalitat de Catalunya” Barcelona Center for Subsurface Imaging (2017 SGR 1662). The work has been supported by projects CONNECT with reference PID2021-128851OA-I00, ATLANTIS (PID2019–109559RB– I00) and THREAT (PID2021-128513OA-I00), funded by the Spanish Ministry of Science and Innovation. ICM has also had funding support of the–Severo Ochoa Centre of Excellence’ accreditation (CEX2019-000928-S), of the Spanish Research Agency (AEI). We acknowledge the use of the computational resources provided by CESGA (Galicia Supercomputing Center) for conducting the simulations reported in this paper. We thank Tim Minshull and an anonymous reviewer for their valuable feedback, which improved this manuscript.

DATA AVAILABILITY

The codes used in this study are available at GitHub repository under the GNU general public license v3.0:

• The new DC-WAS code (Jimenez-Tejero et al. 2024) placed at github repository: https://github.com/ejimeneztejero/DC_WAS.

• The 2-D code version for traveltime tomography used to obtain the results in this study (https://github.com/ejimeneztejero/TOMO2D), and our 3-D version (https://github.com/ejimeneztejero/TOMO3D_dev).

REFERENCES