# Model-based seismic data redatuming with time-domain data-assimilated wavefield reconstruction

Gaoshan Guo, Stéphane Operto

► To cite this version:

Gaoshan Guo, Stéphane Operto. Model-based seismic data redatuming with time-domain data-assimilated wavefield reconstruction. Geophysical Journal International, 2025, 242 (2), 10.1093/gji/ggaf133. hal-05157682

HAL Id: hal-05157682

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

Submitted on 11 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Model-based seismic data redatuming with time-domain data-assimilated wavefield reconstruction

Gaoshan Guo and Stéphane Operto*

University of Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, IRD, Géoazur, Valbonne, 06560, France. E-mail: stephane.operto@geoazur.unice.fr

Accepted 2025 April 3. Received 2025 March 19; in original form 2024 December 21

## SUMMARY

We present a novel model-based seismic data redatuming method based upon time-domain data-assimilated wavefield (DAW) reconstruction. Seismic redatuming refers to the positioning of sources and/or receivers from the acquisition surface to a virtual datum level along with the computation of the data set that would have been recorded by the virtual acquisition. Physically, DAWs approximate the true wavefields by least-squares time-reversal extrapolation of the recorded data in a background subsurface model. To this end, DAWs are first formulated as the sum of the background wavefields and an approximation of the scattered wavefields by the unknown contrasts between the true model and the background model. We estimate the scattered fields via the estimation of their volume scattering sources by fitting in a least-squares sense the differences between the recorded data and the simulated data in the background model, namely the restriction at receivers of the scattered fields. These underdetermined linear source problems involve two main steps: a multidimensional deconvolution of the scattered data by the data-domain Hessian, followed by the propagation of the deconvolved scattered data backward in time. Once the source problems have been solved, DAWs follow by solving the wave equation in the background model using the sum of the experimental sources and the scattering sources as extended sources. Finally, the redatumed data sets, whose accuracy depends on the precision of the background medium between the acquisition surface and the datum, are readily obtained by sampling the DAWs at arbitrary datum levels. The primary challenge in computing accurate DAWs lies in the multidimensional deconvolution of the scattered data, which requires solving a data-domain normal equation with pre-conditioned Krylov-subspace iterative methods, where each iteration requires one forward and one backward simulation. We assess the accuracy of the method and discuss its computational cost with several benchmarks representing various experimental environments (onshore with weathered layers, shallow and deep offshore).

**Key words:** Numerical modelling; Computational seismology; Wave propagation.

## 1 INTRODUCTION

Seismic data redatuming infers dataset that would be recorded at a virtual datum level in the subsurface from the recorded data by the real acquisition device (Mulder 2005; Schuster & Zhou 2006; Wapenaar et al. 2008). This method is valuable when the original acquisition is not suitable for seismic imaging, such as in the presence of irregular topography (Bevc 1997), uneven spatial sampling (Ferguson 2006), complex structures in the overburden (Bakulin & Calvert 2004), and remoteness from the target zone (Valenciano et al. 2009; Borisov et al. 2015; Ravasi et al. 2016; Yuan et al. 2017; Biondi et al. 2023). Recent applications aim at redatuming towed-streamer data on the seabed to increase the offset range over which refracted waves and diving waves from the sedimentary layers are recorded as first arrivals by removing the effect of the wave

layer from the data. Then, the velocity model used for migration can be refined by first-arrival traveltime tomography and full waveform inversion (FWI) (Huot & Singh 2018; Gras et al. 2019; Jimenez-Tejero et al. 2022; Di et al. 2023). Recently, a similar idea has also been applied to ocean bottom node (OBN) data, such that the first arrivals are separated from the coda of the water wave (Li et al. 2024) and increase illumination of the sedimentary layers below the sea bottom (Tejero et al. 2025). Another interesting application deals with the redatuming of the sources on the sea bottom for sparse OBN acquisitions to perform reflection slope tomography (Alerini 2006; Alerini et al. 2009). In this case, the issue is to estimate the slopes (the horizontal component of the slowness vector) of locally coherent events in the recorded data at the OBN positions, which are sparsely deployed. To achieve this goal, the sources are redatumed on the same surface as the OBNs, that is, the sea bottom.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Illustration of source and receiver redatuming. (a) Original acquisition surface at the top of the domain (green solid line) with source (red stars) and receiver (blue triangles) devices. (b) Redatuming of receivers (stars) to a new datum (dash green line). (c) Data rearrangement using reciprocity to set reciprocal sources (red stars with black edge) and receivers (blue triangles with black edges). (d) Redatuming of reciprocal receivers (blue triangles with black edges) to a new datum (green dash line). (e) Rearrangement using reciprocity to set source (red stars) and receiver (blue triangles) devices.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Marmousi II benchmark. (a) Marmousi II model. (b) Smoothed background model. The solid gray line at 420-m depth delineates bathymetry. Yellow dots show original shot positions at 20-m depth while pink inverted triangles show redatumed shot positions on the seabed. Every 20 shot positions are shown in the figure. The three red stars indicate OBN positions on the sea bottom.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Marmousi II benchmark. (a) Convergence curves of PMINRES (inverted triangles) and MINRES (circles) for the data-domain normal system. (b) Condition numbers of data-domain normal system with PMINRES (inverted triangles) and MINRES (circles). The OBN is horizontally located at $x=1$ km (black lines), $x=8.5$ km (blue lines) and $x=16$ km (red lines). The OBN located in the middle of the line leads to the best convergence behavior.</figcaption></figure>

Then, the slopes at a given OBN position are estimated from the densely sampled redatumed shots tied to the other OBNs using the reciprocity of Green's functions.

Typically, the prevalent methods for seismic data redatuming fall into two categories: correlation-based or model-based redatuming (Schuster & Zhou 2006). When buried receivers (in a deviated or horizontal well) exist at the new datum, correlation-based methods create virtual sources based solely on the available data, without any prior knowledge of the medium (Bakulin & Calvert 2004). Conversely, when buried sensors are unavailable, model-based redatuming relies on multidimensional deconvolution of the original data set using Kirchhoff integral or wave-equation continuation operator (Berryhill 1979) based on some prior knowledge of the overburden, specifically, the macro-model between the original datum and the new datum (Planès et al. 2018). Luiken & Van Leeuwen (2020) recast the model-based Marchenko redatuming method as a least-squares ill-posed linear inverse problem, which aims at estimating a matching filter such that the recorded reflection data are matched by the filtered direct wave sampled at the datum level. The matching filter, which is the unknown of the inverse problem, is a convolutional operator representing the impulse response of the subsurface (i.e. the Green functions connecting the receivers at the acquisition plane to the receivers at the datum plane). Therefore, the
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Marmousi II benchmark. Simulated common-OBN gathers in (a) true model and (b-c) background model. Simulations are performed with experimental sources in (a,b) and with extended sources (experimental sources + scattering sources) in (c). (d) Difference between (a) and (c). OBNs are located at (a1–c1) *x*=1 km, (a2–c2) *x*=8.5 km and (a3–c3) *x*=16 km, respectively. Data are plotted with a reduction velocity of 3.5 km/s.</figcaption></figure>

redatuming problem is recast as a multidimensional deconvolution
of the recorded data.

In this study, we recast model-based redatuming in the frame-
work of the so-called wavefield reconstruction inversion (WRI)
(van Leeuwen & Herrmann 2013, 2016; Aghamiry et al. 2019)
or extended-source FWI (ES-FWI) (Wang et al. 2016; Huang et al.
2018; Gholami et al. 2022; Operto et al. 2023b; Guo et al. 2024).
ES-FWI extends the linear regime of FWI by computing wavefields
that match the recorded data arbitrarily well, that is, without cy-
cle skipping. These wavefields can be formulated as the sum of
the background wavefields and an approximation of the scattered
wavefields by the unknown perturbation model, namely the differ-
ence between the unknown true model and the background model.
Accordingly, they can be computed by solving the wave equation in
the background model with extended sources that gather the exper-
imental sources and the estimated volume sources of the scattered
wavefields. The restriction of the scattered wavefields at the re-
ceivers is nothing else than the difference between the recorded
data and the simulated data in the background model, that is, the
data residuals at the current iteration in the FWI context. These
measured scattered data are linearly related to the scattering (or con-
trast) sources by the forward modelling operator (Green functions
sampled at the receivers) built in the background model according

to the Lippmann-Schwinger equation (Operto et al. 2023b). There-
fore, the scattering sources can be estimated in a least-squares sense
by solving an underdetermined linear scattered data fitting problem
with appropriate regularization (see e.g. Shekar & Sethi 2019; Sha-
ran et al. 2019; Aghamiry et al. 2021 for sparsity-promoting regu-
larization of source problems). We have shown that the estimation of
the scattering sources can be recast as a multidimensional predictive
deconvolution of the measured scattered data followed by the adjoint
(reversed time) simulation of the deconvolved scattered data (Op-
erto et al. 2023b). Once these scattering source problems have been
estimated, the wavefields follow by solving the wave equation with
extended sources as above mentioned. We refer these wavefields
to as data-assimilated wavefields (DAWs) (Aghamiry et al. 2020)
because they are inferred from both the wave equation and the rank-
deficient data (or observation) equation (van Leeuwen & Herrmann
2013). If we assume that we know a good approximation of the true
medium between the acquisition level and the datum level, then
sampling the DAWs at the datum level provided the redatumed data
set.

This simple redatuming approach inspired by ES-FWI was first
introduced in the frequency domain by Aghamiry *et al.* (2023).
In contrast, we formulate the method in the time domain, where
the scattering source estimation serves as a proxy for estimating
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Marmousi II benchmark. (a–c) Snapshots of simulated wavefields in (a) true model with experimental sources, and (b, c) background model with (b) experimental sources and (c) experimental + scattering sources. (a1–c1): t = 1.2 s; (a2–c2): t = 2.4 s; (a3–c3): t = 3.6 s; (a4–c4): t = 4.8 s. The yellow line denotes the bathymetry. One can readily compare the level of agreement between the true snapshots and the DA reconstructed snapshots (a versus b panels) above the bathymetry and between the true snapshots and background wavefields (a versus c panels). The OBN is horizontally located at x = 1 km.</figcaption></figure>

the wavefields that best match the true wavefields from surface measurements and the background model. Compared to the frequency domain, time-domain redatuming may offer greater robustness when accounting for the full frequency spectrum (Vargas et al. 2021) while being computationally manageable in 3-D as opposed to the frequency-domain formulation. Still, time-domain DAW reconstruction raises significant challenges, as the multidimensional predictive deconvolution of the scattered data involves solving a normal equation in the data domain using iterative linear algebra solvers. The computational cost of each iteration requires two simulations per shot and receiver (Guo et al. 2024) when both the sources and receivers are redatumed. Thus, selecting the most effective Krylov-subspace solver and pre-conditioner to ensure fast convergence is essential for developing a computationally efficient method.

In the remainder of the paper, we first introduce the theory of model-based seismic data redatuming in the framework of DAW reconstruction. In particular, we introduce the iterative solver for the accurate computation of the deconvolved scattered data. Then, we validate the algorithm with two synthetic benchmarks representative of marine and land applications (Marmousi II and Overthrust, respectively). Next, we present an application of our method on two OBN field data sets. The first corresponds to an ultra-long offset deep-offshore regional data set collected in the Nankai Trough, offshore Japan, while the second data set was recorded in an oil exploration context above the giant Gorgon gas field, offshore Australia. Last, we draw some conclusions and perspectives for future studies.

## 2 METHOD

### 2.1 Seismic data redatuming by data-assimilated wavefield reconstruction

For sake of compact notation, we write the wave equation in matrix form as

$$A(\mathbf{m}^*)\mathbf{u}_s^* = \mathbf{b}_s^*, \quad (1)$$

where A is the wave equation operator (Appendix A), m* denotes the true subsurface medium, b* is the experimental source, u* is the propagated wavefield in m*, and the subscript s labels the s^th source of the seismic survey, s = 1, ..., Ns. Moreover, the superscript * refers to measured or true quantities. We denote an approximate background model by m and we introduce the true model perturbation or contrast object δm = m* − m. Taking a first-order Taylor expansion of A around m gives

$$A(\mathbf{m}^*) = A(\mathbf{m}) + \sum_{i=1}^{M} \delta \mathbf{m}_i^* \frac{\partial A(\mathbf{m})}{\partial \mathbf{m}_i}, \quad (2)$$

where mᵢ denotes the i^th entry of the vector m, and M is the number of model parameters sampling m. The first-order development in eq. (2) is exact because the wave equation is bi-affine in u_s and m, and hence one can choose a subsurface parametrization such that the coefficients of A are linear functions of m. Inserting eq. (2) into eq. (1) leads to the following scattered field equation satisfied by u_s*:
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Marmousi II benchmark. (a) Baseline common-OBN gathers for receivers at the datum plane. (b) Redatumed common-OBN gathers on the sea bottom. (c) Difference between (a) and (b) used to check the accuracy of redatuming. The OBN is horizontally located at (a1–c1) x = 1 km, (a2–c2) x = 8.5 km and (a3–c3) x = 16 km, respectively.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Marmousi II benchmark. Direct comparison between baseline seismogram (solid black lines) and redatumed seismogram (gray lines) along with their difference (dash black lines) for different shots at new datum plane (as shown in Fig. 6). The OBN is horizontally located at x = 1 km. The offset of the seismogram is labelled at the left bottom corner of each panel.</figcaption></figure>

$$ 
\begin{aligned}
\mathbf{A}(\mathbf{m})\mathbf{u}_s^* &= \mathbf{b}_s^* - \sum_{i=1}^{M} \delta \mathbf{m}_i^* \frac{\partial \mathbf{A}(\mathbf{m})\mathbf{u}_s^*}{\partial \mathbf{m}_i} \\
&= \mathbf{b}_s^* - \sum_{i}^{M} \delta \mathbf{m}_i^* \mathbf{f}_s^{(i)} = \mathbf{b}_s^* + \delta \mathbf{b}_s^*,
\end{aligned}
\quad (3) $$

where the virtual impulse sources $f_s^{(i)} = \frac{\partial \mathbf{A}(\mathbf{m})\mathbf{u}_s^*}{\partial \mathbf{m}_i}$ are formed by sampling the incident true wavefields $\mathbf{u}_s^*$ at the position of $\mathbf{m}_i$ by the radiation pattern matrices $\partial \mathbf{A}(\mathbf{m})/\partial \mathbf{m}_i$ (Pratt et al. 1998). The weighted sum of the virtual impulse sources by the perturbation model $\delta \mathbf{m}^*$ forms the true volume scattering or contrast sources $\delta \mathbf{b}_s^*$. Furthermore, we can introduce the scattered wavefields $\delta \mathbf{u}_s^*$ by the volume sources $\delta \mathbf{b}_s^*$ as the difference between the true wavefields and the background wavefields:

$$ \delta \mathbf{u}_s^* = \mathbf{u}_s^* - \mathbf{u}_s(\mathbf{m}) = \mathbf{u}_s^* - \mathbf{A}(\mathbf{m})\mathbf{b}_s^*. \quad (4) $$

These scattered fields satisfy the wave equation in the background medium with the scattering sources $\delta \mathbf{b}_s^*$ as the right-hand side:

$$ \mathbf{A}(\mathbf{m})\delta \mathbf{u}_s^* = \delta \mathbf{b}_s^*. \quad (5) $$

The restriction of $\delta \mathbf{u}_s^*$ at receivers, that is, $\delta \mathbf{d}_s^* = \mathbf{d}_s^* - \mathbf{d}_s(\mathbf{m}) = \mathbf{d}_s^* - \mathbf{S}_s(\mathbf{m})\mathbf{b}_s^*$ represents the measured scattered data, where $\mathbf{d}_s^*$ denote the recorded common-shot gather, $\mathbf{S}_s(\mathbf{m}) = \mathbf{P}_s\mathbf{A}^{-1}(\mathbf{m})$ is the forward modelling operator, and $\mathbf{P}_s$ is the observation operator sampling wavefields at receivers (for stationary-recording system, the subscript s of P can be dropped). According to eq. (5), the measured scattered data $\delta \mathbf{d}_s^*$ and the volume contrast source $\delta \mathbf{b}_s^*$ of the true scattered wavefields are related by the rank-deficient linear operator $\mathbf{S}_s(\mathbf{m})$:

$$ \delta \mathbf{d}_s^* = \mathbf{S}_s(\mathbf{m})\delta \mathbf{b}_s^*. \quad (6) $$
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Overthrust benchmark. (a) Overthrust model. (b) Smoothed background model. Yellow inverted triangles and red stars represent original receiver and source positions at 1-m and 2-m depth, respectively. Inverted triangles and pink stars denote redatumed receiver and source positions, respectively. Two datum levels at 300-m and 500-m depth are considered. Every 10 sources and receivers are shown for the sake of clarity.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Overthrust benchmark. (a) Convergence curves of PMINRES (inverted triangles) and MINRES (circles) for the data-domain normal system. (b) Condition numbers of the data-domain normal system with PMINRES (inverted triangles) and MINRES (circles). Receiver-side redatuming (black lines) and source-side redatuming to 300 m (blue lines) and 500 m (red lines) in depth, respectively.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Overthrust benchmark. Common-shot gathers are shown. (a) Recorded data. (b, c) Simulated data in the background model without (b) and with (c) scattering sources. (d) Differences between (a) and (c) highlight the accuracy of DAWS at receivers. Data are plotted with a reduction velocity of 3.5 km/s.</figcaption></figure>

Therefore, the true volume contrast sources $\delta\mathbf{b}_s^*$ can be estimated by solving the above equation in the least-squares sense with appropriate regularization:

$$ \min_{\delta \mathbf{b}_s} \|\mathbf{S}_s(\mathbf{m})\delta \mathbf{b}_s - \delta \mathbf{d}_s^*\|_2^2 + \mu \|\delta \mathbf{b}_s\|_2^2, \quad s = 1, \dots, N_s, \quad (7) $$

where a simple damping regularization is used here and $\mu \in \mathbb{R}_+$ is a penalty parameter. The value of $\mu$ is set such that the recorded scattered data are matched as accurately as possible. Typically, a small value is used when $\mathbf{m}$ is far away from $\mathbf{m}^*$, namely, when the scattering generated by the model perturbation in $\mathbf{m}$ is significant.

This gives

$$ \begin{aligned} \delta \mathbf{b}_s^{(\text{da})} &= \mathbf{S}_s(\mathbf{m})^T (\mathbf{S}_s(\mathbf{m})\mathbf{S}_s(\mathbf{m})^T + \mu \mathbf{I})^{-1} \delta \mathbf{d}_s^* \\ &= \mathbf{S}_s(\mathbf{m})^T \mathbf{H}_{d_s}(\mathbf{m})^{-1} \delta \mathbf{d}_s^*, \\ &= \mathbf{S}(\mathbf{m})^T \delta \mathbf{d}^{(d)}. \end{aligned} \quad (8) $$

This would introduce the quantity $\delta\mathbf{d}^{(d)}$ whose computation is discussed in equation 10. Where the superscript (da) refers to the data-assimilation (DA) procedure underlying this source problem. Since the source problem is underdetermined, we formulate the Hessian $\mathbf{H}_{d_s}(\mathbf{m})$ in the data domain. The above equation shows that $\delta\mathbf{b}_s^{(da)}$
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Overthrust benchmark. Receiver-side redatuming. (a) Simulated data in the true model for receivers located on the datum planes as virtual recorded data. (b) Redatumed common-shot gathers on the receiver-side datum plane. (c) Differences between (a) and (b) for receiver-side redatuming assessment. The new datum plane is 300 m (a1–c1) and 500 m (a2–c2) in depth, respectively. Each panel is plotted with the same amplitude scale.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Overthrust benchmark. Receiver-side redatuming. Direct comparison between baseline seismogram (solid black lines) and redatumed seismogram (gray lines), along with their differences (dash black lines) for receivers at new datum plane at 300 m (a1–d1) and 500 m (a2–d2) in depth. The offset of the seismogram is labelled at the left bottom corner of each panel.</figcaption></figure>

is computed by reversed-time modelling of the weighted scattered data by the inverse data-domain Hessian. For the sake of completeness, the implementation of this adjoint-state equation is reviewed in Appendix A when the state equation is the 2-D first-order velocity-stress acoustic wave equation, which has been used in this study to perform the numerical experiments. We also review in Appendix B

the expression of the contrast source $\delta b_s^*$ for the simulation of the scattered field $\delta u_s^*$, eq. (5), with the first-order velocity-stress acoustic wave equation. Then, DAWs that best approximate the true wavefields $u_s^*$ are computed by solving the wave equation in $\mathbf{m}$ with an extended source $b_s^{(da)}$:

$$A(\mathbf{m})u_s^{(da)} = b_s^* + \delta b_s^{(da)} = b_s^{(da)}. \qquad (9)$$

The deconvolved scattered data $\delta d_s^{(d)}$ equation 8, are the solution of the following normal system:

$$H_{ds}(\mathbf{m})\delta d_s^{(d)} = \delta d_s^*, \quad s = 1, \dots, N_s. \qquad (10)$$

Therefore, the computation of the DAWs involves three main steps (Algorithm 1): The first is a multidimensional predictive deconvolution of the measured scattered data by the data-domain Hessian, eq. (10). The second is the reversed-time (adjoint) modelling of the deconvolved scattered data to build the scattering sources $\delta b_s^{(da)}$, eq. (8). The third is the forward modelling of $u_s^{(da)}$ with the extended sources $b_s^{(e)}$, eq. (9). Summing $b_s^*$ with $\delta b_s^{(da)}$ in $b_s^{(da)}$ mostly cancels the contribution of $d_s(\mathbf{m})$ at the benefit of $d_s^*$ from the measured scattered data (because $d_s(\mathbf{m})$ optimally focuses at $b_s^*$ after reversed-time modelling and hence both cancel out by subtraction), which implies that the DAWs mostly correspond to the backward/forward simulated deconvolved recorded data (Operto et al. 2023b). By virtue of wave invariance for time reversal, they match the recorded data arbitrarily well when a small value of $\mu$ is used, but they deviate from the true wavefields as they are propagated backward in time away from the receivers during the adjoint step. However, if $\mathbf{m}$ is a good approximation of $\mathbf{m}^*$ between the acquisition plane and the datum plane, then sampling $u_s^{(da)}$ on the datum plane provides accurate redatumed data set. Compared to the approach of Luiken & Van Leeuwen (2020), the unknown of the deconvolution problem is not the Green's function from the acquisition plane to the datum plane but the deconvolved recorded scattered data. Then, the redatumed data is extracted at any datum plane from the DAWs $u_s^{(da)}$. Using reciprocity, both sources and receivers can be redatumed with this approach (Fig. 1). Moreover, this redatuming strategy can be
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Overthrust benchmark. Common-receiver gathers are shown. (a) Common-receiver gathers after receiver redatuming as virtual recorded data (receiver-side redatumed data when redatumed receivers and original shots are processed as reciprocal sources and reciprocal receivers, respectively). (b, c) Simulated data in the background model (b) without and (c) with scattering sources. (d) Differences between (a) and (c) highlight the ability of DAWs to match recorded data from an inaccurate background model after receiver-side redatuming. Reciprocal sources at depth (a1–d1) 300 m and (a2–d2) 500 m and reciprocal receivers at original shot positions.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Overthrust benchmark. Source-side redatuming. Common-receiver gathers are shown. (a) Simulated data in the true model for source and receivers at the datum levels as recorded data. (b) Receiver-side and source-side redatumed data. (c) Difference between (a) and (b) for assessment of the receiver and source redatuming. The new datum plane is (a1–c1) 300 m and (a2–c2) 500 m in depth, respectively.</figcaption></figure>

applied to any kind of acquisition geometry and arrival type (reflec-
tions, diving waves,...). Also, our method handles surface multiples
as long as a free-surface boundary condition is implemented in the
forward modelling engine.

**2.2 Iterative solver for data-domain normal equation**

Previous studies approximate the inverse of the data-domain Hes-
sian in eq. (8) with scaled identity matrix (Gholami *et al.* 2022)

or matching filters (Guo et al. 2024). While these methods strike
a balance between accuracy and computational efficiency, particu-
larly considering the large number of iterations required for FWI,
their precision may be insufficient for addressing the seismic data
redatuming problem. In the context of redatuming, achieving highly
accurate DAWs is essential. As a result, we employ iterative solvers
to solve the normal equation involving the data-domain Hessian as
the left-hand side operator, eq. (10). Let us note that the scattering
sources $\delta b_s^{(\text{da})}$ might be directly computed by solving eq. (6) with
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Overthrust benchmark. Source-side redatuming. Direct comparison between baseline seismogram (solid black lines) and redatumed seismogram (gray lines), along with their difference (dash black lines) for different receivers and sources at new datum plane at 300 m (a1–e1) and 500m (a2–e2) in depth. The offset of the seismogram is labelled at the left bottom corner of each panel.</figcaption></figure>

**Algorithm 1** Model-based seismic data redatuming with DAWs

1: Input: **m**, **d*** and **b***.

2: Compute recorded scattered data $\delta\mathbf{d}^*(\mathbf{m}) = \mathbf{S}(\mathbf{m})\mathbf{b}^* - \mathbf{d}^*$.

3: Multidimensional deconvolution with MINRES solver $(\mathbf{S}(\mathbf{m})\mathbf{S}^T(\mathbf{m}) + \mu\mathbf{I})\delta\mathbf{d}^{(d)}(\mathbf{m}) = \delta\mathbf{d}^*(\mathbf{m})$.

4: Compute source extension $\delta\mathbf{b}^{(da)} = \mathbf{S}(\mathbf{m})^T\delta\mathbf{d}^{(d)}(\mathbf{m})$.

5: Compute DAWs $\mathbf{A}(\mathbf{m})\mathbf{u}^{(da)} = \mathbf{b}^* + \delta\mathbf{b}^{(da)}$.

6: Sample DAWs at the new datum.

7: Repeat the above process using reciprocity to perform source redatuming.

The Least Squares QR method (LSQR) (Paige & Saunders 1982). However, this approach is impractical due to the high dimensionality of $\delta\mathbf{b}^{(da)}$. This is the reason why we favour a data-domain implementation of the least-squares problem for $\delta\mathbf{b}^{(da)}$. The source problem remains however ill-conditioned. Guo et al. (2022) explore the conjugate gradient (CG) method but conclude that the CG method was quite unstable when handling complex data sets and strong surface multiples. We found that the minimal residual method (MINRES) (Choi et al. 2011) was providing a more robust alternative for our application. The use of MINRES in this context has been extensively discussed by Choi et al. (2011) and Chin-Lung Fong & Saunders (2012). Moreover, we apply a physics-inspired left pre-conditioner designed as a basic frequency-domain filter F to restrict the frequency band of the data residuals after each MINRES iteration:

$$ \mathbf{F}\mathbf{H}_d(\mathbf{m})\delta\mathbf{d}^{(d)}(\mathbf{m}) = \mathbf{F}\delta\mathbf{d}^*(\mathbf{m}). \quad (11) $$

The filter F performs an implicit regularization by enforcing the frequency band of the deconvolved scattered data to match that of the measured counterpart, hence preventing the onset of high-frequency artifacts. A similar approach has been proposed by Vargas et al.

(2021) to perform multidimensional deconvolution of seismic data for redatuming application. It is important to note that the multidimensional filter F is a simplistic approximation, as it does not account for the multi-dimensional effects of $\mathbf{H}_d(\mathbf{m})$ in the time-offset domain. Therefore, future works should target the development of a more effective pre-conditioner for this data-domain source problem. Designing a pre-conditioner based on asymptotic theory may be a promising area of investigation by exploiting the mathematical form of pseudo-differential operators (Symes 2008; Hou & Symes 2015; Farshad & Chauris 2021). In the following synthetic and real data tests, we illustrate the impact of the pre-conditioner F on the convergence speed by comparing the performances achieved by MINRES whether or not the pre-conditioner is used (tests labelled by MINRES and PMINRES, respectively).

## 3 NUMERICAL EXAMPLES

In the following, we validate our redatuming approach against two synthetic benchmarks. The first one corresponds to the Marmousi II model to address redatuming challenges for ocean-bottom acquisition. The second test relies on a 2-D section of the 3-D European Association of Geoscientists and Engineers (EAGE) / Society of Exploration Geophysicists (SEG) Overthrust model to assess redatuming for more challenging land applications involving highly heterogeneous media near the surface with a weathered layer. Acoustic modelling is performed using a velocity-stress finite-difference time-domain method (Virieux 1984 and Appendix A) incorporating a free-surface boundary condition at the top of the computational grid and a perfectly matched layer absorbing boundary condition on the remaining sides. Sources and receivers are positioned at arbitrary positions in the grid using Kaiser-windowed sinc functions (Hicks 2002). The stopping criterion of MINRES for the solution of the normal system (eq. 10) is provided by the relative backward error (Tournier et al. 2022):

$$ \epsilon_{\text{error}} = \frac{\| \mathbf{H}_d(\mathbf{m}) \delta \mathbf{d}^{(d)} - \delta \mathbf{d}^* \|_2^2}{\| \delta \mathbf{d}^* \|_2^2}. \quad (12) $$

We use $\epsilon_{error} = 1 \times 10^{-2}$, which guarantees that the re-simulated $\delta\mathbf{d} = \mathbf{H}_d(\mathbf{m})\delta\mathbf{d}^{(d)}$ fit accurately $\delta\mathbf{d}^*$. Moreover, we adapt dynamically the value of $\mu$ in MINRES iteration with the following expression:

$$ \mu = \alpha \frac{\| \mathbf{S}(\mathbf{m})\mathbf{S}(\mathbf{m})^T \delta \mathbf{d}^{(d)} \|_2}{\| \delta \mathbf{d}^{(d)} \|_2^2}, \quad (13) $$

where $\alpha = 1 \times 10^{-4}$. The frequency-domain filter F is a minimum-phase bandpass filter with Hanning branches with four input parameters: ($f_{low}$, $f_{hig}$, $l_{low}$, $l_{hig}$), where $f_{low}$ and $f_{hig}$ are low-cut and high-cut frequencies and $l_{low}$ and $l_{hig}$ are the lengths of the low-cut and high-cut Hanning branches. For the two following synthetic benchmarks involving noiseless data, we use a low-pass filter for F, and hence $f_{low}$ and $l_{low}$ are set to 0. For each test, we present the convergence curves of MINRES with and without pre-conditioner to demonstrate the impact of pre-conditioning. Additionally, we plot the condition numbers for each iteration of MINRES. We recall that higher condition numbers indicate ill-conditioning of the normal system, which will require more iterations to achieve convergence.

### 3.1 Marmousi II model

We first evaluate the proposed time-domain redatuming method using the Marmousi II model (Fig. 2a). The background model is
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Nankai data. Acquisition map of the SFJ-OBS survey. The shot line and the OBS line are plotted as a solid black line and a dotted red line, respectively. The yellow circles denote three OBNs processed in this study.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 17.** Nankai data. Background model for redatuming of OBN data on the sea bottom (gray line). The yellow dots denote the original shot positions, while the pink inverted triangles represent the redatumed shot positions. Every 10 shot positions are shown for the sake of clarity. Red stars indicate the positions of the three OBNs used for redatuming assessment.</figcaption></figure>

the smoothed version of the true model using the Gaussian filter with horizontal and vertical correlation lengths set to 350 m. To perform redatuming, we assume the depth of the sea bottom (set to 420 m) and the water wave velocity (set to 1500 m s⁻¹) known (Fig. 2b). We compute the recorded data (i.e. simulated data in the true model) using a Ricker wavelet with a peak frequency of 5 Hz as the source signature. The simulation length is 11.6 s. The sea-bottom fixed-spread acquisition involves a line of 851 shots (shot interval = 20 m) positioned at a depth of 20 m below the sea surface and we consider three OBN positions on the sea bed at x=1 km, 8.5 km and 16 km. We exploit the reciprocity of Green's functions for treating shooting positions as reciprocal receivers and OBNs as reciprocal sources, and we perform the redatuming of the shots on the seabed. Note that the dense shot sampling is necessary to prevent aliasing artifacts during the redatuming as it is the case for methods based on Kirchhoff integral or wave-equation continuation operators (Jimenez-Tejero et al. 2022). This dense shot sampling is generally implemented during towed-streamer and OBN field experiments.

The convergence curves of MINRES are shown in Fig. 3(a) for the three OBN positions when the pre-conditioner is used or not. For the pre-conditioner, we design the low-pass filter $F$ using $f_{\text{hig}} = 14.5$ Hz and $l_{\text{hig}} = 2.5$ Hz. PMINRES (inverted triangles) accelerates the convergence by a factor of around 2 compared to MINRES (circles). Moreover, we show in Fig. 3(b) how the pre-conditioner efficiently decreases the condition number of the data-domain operator consistently with the decrease of the number of iterations. We also note that faster convergence is obtained for the OBN located in the middle of the model. This likely results from the more limited offset range involved in the simulations. Before showing the redatuming results, we first illustrate the ability of the wavefield reconstruction process to fit the recorded data $d_s^*$ with the inaccurate background model $m_b$ shown in Fig. 2(b). The Figs 4(a), (c) and (d) show the close fit between $d_s^*$ and the computed data $d_s^{(da)} = S_s(m_b)b^{(da)}$ in $m_b$ with the extended source $b_s^{(da)}$, where it is recalled that $b_s^{(da)}$ represents the extended sources of the data-assimilated wavefields defined in eq. (9). This contrasts with the significant mismatches between $d_s^*$ and the computed data in $m_b$ with the experimental sources $b_s^*$,
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Nankai data. Common-OBN gathers. (a) Recorded data. (b, c) Simulated data in background model (b) without and (c) with scattering sources. (d) Difference between (a) and (c). Positions of OBNs are (a1–d1) x = 18.037 km, (a2–d2) x = 60.037 km, and (a3–d3) x = 94.1 km. Data are plotted with a reduction velocity of 7 km s⁻¹.</figcaption></figure>

namely $\mathbf{d}_s = S_s(\mathbf{m}_b)\mathbf{b}_s^*$ (Figs 4a and b). Additionally, several snapshots of the true wavefields, background wavefields, and DAWs are displayed in Fig. 5. Notably, the complexity of the wavefields in the shallow region is due to multiple interactions between the wavefields, the free surface and the sea bottom. Qualitatively, the DAWs

align well with the true wavefields down to the sea bottom where we want to redatum the data. However, their kinematic accuracy gradually diminishes with depth in the smoothed background model, and truncation artifacts emerge due to the limited spread of the receivers (Fig. 5c3). On the whole, the DAWS match much better the
<figure><img src="image_20.png" /><figcaption>**Figure 19.** Nankai data. Common-OBN gathers. (a) Recorded data. (b) Recorded data after redatuming of the shots at the sea bottom. Positions of OBNs are (a1–b1) *x* = 18.037 km, (a2–b2) *x* = 60.037 km and (a3–b3) *x* = 94.1 km.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20.** Nankai data. (a) Convergence curves of MINRES solver with PMINRES (inverted triangles) and MINRES (circles) for the data-domain normal system. (b) Condition numbers of the data-domain normal system with PMINRES (inverted triangles) and MINRES (circles). Positions of OBNs are *x* = 18.037 km (black lines), *x* = 60.037 km (blue lines), and *x* = 94.1 km (red lines).</figcaption></figure>

true wavefields than the background wavefields. Next, we sample the DAWs at the new datum plane (i.e. the sea bottom) (Figs 6a–c) for the three reciprocal sources. The accuracy of the redatuming is evaluated by comparing the baseline data (the simulated data in the true model sampled at the redatumed receivers on the sea bottom) with the redatumed data, namely the sampled DAWS on the sea bottom (Fig. 6). The kinematic agreement is quite good while moderate amplitude mismatches are observed in the multiscattered data between the free surface and the sea bottom (Fig. 6c). The level of data fit is further analysed in Fig. 7, which presents a direct comparison between the baseline data and the redatumed data along several vertical profiles. This comparison confirms that the discrepancies are primarily amplitude-related rather than being due to traveltime

errors. These amplitude errors likely stem from inaccuracies in the velocity contrast at the sea bottom in the background model.

## 3.2 Overthrust model

We now consider a 2-D section of the 3-D SEG/EAGE Overthrust model (Aminzadeh et al. 1996) (Fig. 8a) to demonstrate the application of redatuming in an onshore environment, where both source-side and receiver-side redatuming can be applied for static correction. We generate the background model by smoothing the true model using a Gaussian filter of horizontal and vertical correlation lengths of 300 and 200 m, respectively (Fig. 8b). The source is a Ricker wavelet with a peak frequency of 5 Hz, and redatuming is performed within a frequency bandwidth of 0 to 15 Hz. Accordingly, we design the low-pass filter F for pre-conditioning with $f_{\text{hig}} = 15$ Hz and $l_{\text{hig}} = 1.5$ Hz. The fixed-spread acquisition consists of a line of receivers spaced 75 m apart at 1-m depth, and a line of sources spaced 150 m apart at 2-m depth. We consider two datum surfaces at 300- and 500-m depth to assess the sensitivity of source and receiver redatuming to the accuracy of the background model between the acquisition surface and the datum surface.

Receiver-side redatuming is performed first. The convergence of the MINRES solver for the multidimensional deconvolution of the measured scattered data is shown in Fig. 9 (black line with circles). Again, the pre-conditioner decreases the required number of iterations to reach the stopping criterion of iteration by a factor of more than two (black line with triangles vs circles). Fitting the recorded data with the DAWs is a pre-requisite for successful redatuming. Fig. 10 shows that the fit between the recorded data and the simulated data in the background model with the extended
<figure><img src="image_22.png" /><figcaption>**Figure 21.** Gorgon data. Background model for the redatuming of the OBN data on the sea bottom (gray line). The yellow dots denote the original shot positions, while the pink inverted triangles represent the redatumed shot positions. The shot positions are shown with a reducing density factor of 10. The red stars indicate the positions of the OBNs.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 22.** Gorgon data. Common-OBN gathers of (a) recorded data and simulated data in the background model (b) without and (c) with source extension. (d) Difference between (a) and (c) used to validate the accuracy of DAWs. The OBN horizontally is located at (a1–d1) *x* = 4.674 km, (a2–d2) *x* = 9.925 km and (a3–d3) *x* = 16.3 km. Data are plotted with a reduction velocity of 3.5 km/s.</figcaption></figure>

sources is quite accurate. The DAWs are then sampled at the two new datum planes (Fig. 11). As for the Marmousi test, we assess the accuracy of the redatuming process by comparing the redatumed data (Figs 11b1 and b2) with the baseline simulated data in the true model for the virtual acquisition at the datum level (Figs 11a1 and a2). As expected, the accuracy of the redatuming decreases slightly as the datum depth increases (Fig. 12). While the traveltime match remains satisfactory due to the good kinematic accuracy of the background model, the amplitude discrepancies become more notable with increasing depth due to the accumulation of inaccuracy of DAWs as they are back-propagated in the inaccurate background model.

After receiver redatuming, we proceed with the source redatuming at the new datums using the reciprocity of Green's functions. To do this, we process the receiver-side redatumed data as recorded data and we datum the sources by processing the redatumed receivers as reciprocal sources and the original sources as reciprocal receivers. We first assess the data match between simulated data with scattering sources and the new redatumed recorded data (Fig. 13). The convergence of MINRES is accelerated by pre-conditioner as shown in Fig. 9 (blue and red lines).

Subsequently, the DAW is sampled at the new datums (Fig. 14). Notably, the kinematic agreement between the true redatumed data and the reconstructed ones remains fairly good for datum 300 m, while it is less accurate for the 500-m datum as expected (Fig. 15).

## 4 FIELD DATA APPLICATION

We next present the application of redatuming to two long-offset OBN data sets: one at the regional crustal scale from the Nankai Trough, offshore Japan and another one at the oil exploration scale from the Gorgon gas field, Australia. In both cases, we perform redatuming to move the shot positions on the bathymetry. For these real data applications, the values of $\mu$ and $\alpha$ are set to $5 \times 10^{-2}$ and $5 \times 10^{-4}$, respectively.

### 4.1 Case study: ultra-long offset OBN data from the SFJ-OBS experiment

We present an application of the time-domain redatuming method to real ultra-long offset OBS data from the Seize France Japan
<figure><img src="image_24.png" /><figcaption>**Figure 23.** Gorgon data. Common-OBN gathers of (a) recorded data and recorded data after redatuming of the shots to the sea bottom. The OBN is horizontally located at (a1–b1) *x* = 4.674 km, (a2–b2) *x* = 9.925 km and (a3–b3) *x* = 16.3 km.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 24.** Gorgon data. (a) Convergence curves of PMINRES (inverted triangles) and MINRES (circles) for the data-domain normal system. (b) Condition numbers of the data-domain normal system with PMINRES (inverted triangles) and MINRES (circles). The OBN is horizontally located at *x* = 4.674 km (black lines), *x* = 9.925 km (blue lines) and *x* = 16.3 km (red lines).</figcaption></figure>

OBS (SFJ-OBS) experiment, offshore Japan (Dessa et al. 2004; Operto et al. 2006; Górszczyk et al. 2017). The objective of the experiment is to image the deep crustal structure of the eastern Nankai trough down to Moho depths. The 2-D acquisition consists of a line of 100 OBNs spaced 1 km apart and a line of 1404 pressure sources spaced 100 m apart, both lines being perpendicular to the subduction front (Fig. 16). A kinematically accurate background model is constructed using first-arrival tomography with the correct bathymetry (Górszczyk et al. 2017) (Fig. 17). We assess redatuming for three OBNs located near the North–West end of the line (shallow water), in the middle of the line and near the South–East end of the line (deep water) (Figs 16 and 17). The data are bandpass filtered between 1.5 and 8.5 Hz. This band-pass filter is also used to build the pre-conditioner filter using $f_{\text{low}} = 1.5$ Hz, $f_{\text{hig}} = 8.5$ Hz, $l_{\text{low}} = l_{\text{hig}} = 1$ Hz. We begin by estimating the source signatures using the method of Pratt (1999) using the FWI model of Górszczyk et al. (2017) as the background model. Then, the OBNs are processed as reciprocal sources and we redatum the reciprocal receivers (i.e. shots) at the sea bottom.

Before redatuming, we check the data match between the observed data and the simulated data in the traveltime tomography model with extended sources (Fig. 18). Despite the complex anatomy of the ultra-long-offset data, we fit the recorded data fairly well both from kinematic and dynamic viewpoints (Figs 18a, c and d). For the sake of comparison, we show also the simulated data in the traveltime tomography model without source extension. In this case, we achieve a good kinematic match of the first arrivals (Figs 18a and b). However, the amplitudes of these arrivals are not matched and the later arrivals, in particular the post-critical reflections from the subducting slab, are lacking in the simulated data. Next, we sample DAWs at the seabed to generate redatumed data sets (Fig. 19). The redatuming of the shots at the seabed efficiently removes the footprint of the highly heterogeneous bathymetry, hence making the interpretation of the arrivals easier. In particular, one can see the simpler anatomy of the post-critical reflection from the subducting slab and the refraction from the upper mantle (Pn wave) beyond 35 km offset in Fig. 19(b1). We note also that high-frequency artifacts can pollute redatumed first arrivals that have been extracted from the coda of the direct water wave at short offsets, in particular for OBNs in deep water. These artifacts probably result from the relative coarse shot interval (100 m). They should be however removed efficiently by dip filtering. Overall, the quality of the redatumed data is enough to pick first and late locally coherent events arrivals and target joint first-arrival and reflection slope-tomography applications with the approach of Sambolian *et al.* (2021) and Alerini *et al.* (2009). For this application, PMINRES iterations reach the stopping criterion after approximately 70 iterations (Fig. 20a). Also, the pre-conditioner greatly decreases the condition number (Fig. 20b). Due to the complexity of the ultra long-offset data anatomy, the condition number is much higher than that of the previous numerical tests. In this case, the pre-conditioner reduces the number of MINRES iterations by 25 per cent. Additionally, the slowest convergence is obtained for the landward OBN gather (OBN 17) due to longer offsets and shallower water.

## 4.2 Case study: OBN data from the Gorgon gas field

Next, we present an application of the time-domain redatuming method on the field OBN data from the Gorgon gas field in the northwestern shelf of Australia (Operto et al. 2023a). The Gorgon field is an elongated horst block, partitioned into fault blocks by intra-horst en-echelon faults, located at 3.6–4.3 km subsea. The bathymetry varies between 200 and 1000 m and shows a sharp scarp at around 5 km distance. The shots are finely spaced 18.5 m apart. The background model, which incorporates the correct bathymetry, has been developed by reflection tomography methods (Fig. 21). The raw data is bandpass filtered between 1.5 and 24 Hz (Fig. 22a). The low-cut and high-cut frequencies of the pre-conditioner filter are 1.5 and 24 Hz with 1 and 2 Hz wide Hanning branches. Like for the Nankai benchmark, we estimate the source wavelets from the background model by least-squares frequency-domain data fitting (Pratt 1999). We perform shot redatuming at the sea bottom for three OBNs gathers. We first assess the accuracy of the DAWs by evaluating the data fit between the observed data and the simulated data (Figs 22a and c). Again, the DAWs match fairly well all the recorded arrivals both from the kinematic and dynamic viewpoints, while the simulated data in the background model with the experimental sources match only the kinematics of the first arrivals. The redatumed data are shown in Fig. 23. Compared to the results
<figure><img src="image_26.png" /><figcaption>**Figure 25.** Layer model. Snapshots of true (a), background (b), and DAW (c) wavefields in the layer model. Dash line: datum level at 3.5 km depth. Shaded area: 1-km thick layer. The source is at 4.5-km distance and at 150 m below the surface. The receiver line is at 200 m below the surface. Velocities in the homogeneous background and layer are 2 and 4 km s⁻¹, respectively. The source signature is a 5-Hz Ricker wavelet. Absorbing boundary conditions are implemented along the four edges (no surface multiples).</figcaption></figure>

of the Nankai case study, the redatumed data are free of artifacts supporting the fact that a fine sampling of the acquisition is beneficial for redatuming. The geometry of the short-offset traveltime curves clearly shows that the direct wave hyperbola has collapsed at its apex after redatuming hence revealing the diving waves from the shallow sedimentary structure. We see also that the bathymetric effects generated by the scarp are emphasized after redatuming. For this high-quality dataset, PMINRES needs around 60 iterations to reach the stopping criterion of iteration (Fig. 24). However, the number of iterations is reduced only by around 10 per cent when the pre-conditioner is used. The OBNs of the Gorgon survey are spaced 350 m apart preventing slope picking in common-shot gathers. However, this OBN sampling should be enough to allow for joint first-arrival and reflection slope tomography where shot redatuming will be beneficial for picking the first break of shallow

diving waves while allowing for the building of a sparse data set of reflection picks (Sambolian 2021).

## 5 DISCUSSION

In this study, we propose a novel seismic data redatuming approach inspired by the WRI method (van Leeuwen & Herrmann 2013). In this framework, the reconstructed DAWs are formulated as the sum of the background wavefields and the scattered wavefields by the contrasts between the true and the background subsurface media. To compute DAWs, we first estimate the volume sources of these scattered fields by fitting in the least-squares sense their restriction at receivers, namely the differences between the recorded data and the simulated data in the background model. This represents a highly
<figure><img src="image_27.png" /><figcaption>**Figure 26.** Layer model. (a) Simulated data in true model. (b) Simulated data in background model without scattering sources. (c) Same as (b) but with scattering sources. (d) Difference between (a) and (c).</figcaption></figure>

underdetermined linear source problem. The scattering sources are computed by propagating backward in time the deconvolved scattered data. Then, we compute DAWs by solving the wave equation in the background model by adding the scattering sources to the experimental sources in the right-hand side. Summing the experimental sources and the scattering sources removes the contribution of the back-propagated simulated data at the benefit of the back-propagated recorded data in the latter. This means that, to first order, the DAWs are computed by backward/forward modelling of the deconvolved recorded data in the background model. Once DAWs have been computed, redatuming simply consists in sampling them at the datum plane. The primary computational challenge of our method lies in estimating the scattering sources, which entails multidimensional deconvolution of the scattered data by the data-domain Hessian of the source problem. We propose solving approximately the normal system associated with this deconvolution process with the MINRES iterative solver, where each iteration requires one forward simulation and one adjoint simulation in the full subsurface medium. As a pre-conditioner, we use a simple frequency-domain bandpass filter to prevent boosting high-frequency noise in the deconvolved scattered data.

Our method is versatile as it can handle any kind of acquisition geometry and wave types (diving waves, refraction, reflection). As with all redatuming methods, it requires a fine sampling of the redatumed acquisition dimension to prevent aliasing artifacts. Although we assess the method in 2-D configuration, it can be readily extended to 3-D geometries. In contrast to most of the wavefield extrapolation-based methods, the redatumed sources and receivers are omnidirectional in the sense that the redatumed data account for scattered waves from the region above the new datum plane (e.g. surface multiples). However, the accuracy with which our model-driven approach accounts for the scattered waves above the datum plane depends on the accuracy of the background model above this plane. While the DAWs can fit the recorded data arbitrarily well via the data-assimilation procedure, their accuracy progressively decreases away from the receivers since they are basically computed by propagating the recorded data backward in time in the background model. All the arrivals (including multi-scattered arrivals as internal multiples) will be redatumed accurately at least from the kinematic viewpoint as long as the background medium predicts accurately their paths between the acquisition surface and the datum plane. Conversely, if the missing contrasts between the

true medium and the background medium generate reflections and multi-scattered waves, these arrivals will not be processed properly as discussed hereafter. We stress that the accuracy of the redatuming is far less sensitive to the accuracy of the background model below the datum plane because the backward/forward propagated recorded data will match the true wavefield at the datum plane by virtue of wave invariance for time reversal if the medium is accurate between the datum plane and the acquisition surface. The specification of an accurate background model above the datum plane is indeed satisfied for shot redatuming at the seabed in marine environments as illustrated in this study. For land applications, the datum levels should be shallow enough for static corrections to prevent a large accumulation of errors. In contrast, our approach will be more sensitive to the accuracy of the overburden model for target-oriented imaging where the datum level is deeper. For example, if the overburden contains salt bodies (see Vargas et al. 2021 for a benchmark), our method cannot handle accurately multi-scattered waves occurring within the salt bodies when the background model is smooth. Moreover, the primary reflections from the salt bodies, which do not propagate beneath the datum plane, will generate spurious “acausal” arrivals in the redatumed data. To illustrate these issues more quantitatively, we perform a numerical example with a homogeneous model containing a 1-km thick layer between 1.5- and 2.5-km depths (Figs 25–27). The velocities in the homogeneous medium and in the layer are 2 and 4 km s⁻¹, respectively. We compute the DAWs in a background model obtained by smoothing the true model with a vertical Gaussian filter of 200-m correlation length, and we sample them at 3.5-km depth beneath the layer. Snapshots of the true wavefields are shown in Fig. 25(a). We check in Fig. 26 that the DAWs match quite accurately the recorded data at acquisition surface. However, the DAWs poorly match the true wavefields at the datum plane (Fig. 27). First, amplitudes of the transmitted wavefields across the layer are not matched well. Moreover, spurious signals are recorded before the transmitted wavefields across the layer. Third, internal multiples from the layer are mispositioned in time and space with underestimated amplitudes. The snapshots of the DAWs allow us to understand intuitively how the scattering sources build up across time (Fig. 25c). Artificial deep virtual sources generate part of the reflections from the top and the bottom of the layer that cannot be simulated by the smooth background model (Fig. 25c2, black arrows). More precisely, they correct the inaccurate amplitudes of the simulated reflections at
<figure><img src="image_28.png" /><figcaption>Figure 27. Layer model. (a) Simulated data in true model for receivers at datum level. (b) Sampling of DAW at datum level. (c) Difference between (a) and (b).</figcaption></figure>

short angles in the smooth background model (Figs 25b4 and c4, black arrows). These fictitious upgoing wavefields from beneath the datum plane reach the latter before the downgoing wavefields triggered by the experimental sources. They generate in the redatumed data acausal artifacts associated with upgoing reflected wavefields from above the datum plane (Fig. 27). Furthermore, the snapshots show that the DAWs fail to reconstruct the internal multiples in time and space since their paths cannot be predicted by the smooth background model (Fig. 25 a6–a8, black arrows). According to this analysis, the Marchenko method is superior in its ability to account for multi-scattering occurring within the overburden (Ravasi 2017; Cui et al. 2018; Lomas & Curtis 2019; Wapenaar et al. 2021; Shoja et al. 2024). A comparative analysis of Marchenko redatuming and DAW methods deserves further investigation to find possible links and complementarity between the two approaches. For example, Hajjaj et al. (2022) estimate the focusing functions in Marchenko redatuming using a penalty objective function inspired by WRI.

Another important difference with other model-based redatuming methods is that it requires a background velocity model covering the domain in which the recorded data is propagated. This is required to compute synthetic data from which we infer the measured scattered data, namely the difference between the recorded and synthetic data. Then, it is required to perform Hessian-vector products during the multidimensional deconvolution of the scattered data by backward/forward simulations of the input vector into the full domain. This Hessian captures the correlation between all possible scattered arrivals in the subsurface by neighbour receivers (Operto et al. 2023b, their fig. 6). Finally, we need to perform one additional backward/forward simulation of the deconvolved scattered data to compute the redatumed scattered data. Let us note that the multidimensional deconvolution of the recorded scattered data followed by backward/forward simulations is equivalent to applying the data-domain resolution operator to the measured scattered data to generate the simulated redatumed scattered data $\delta\mathbf{d}^{\text{rd}}$.

$$ \delta\mathbf{d}^{\text{rd}} = \mathbf{S}^{\text{rd}}(\mathbf{m})\mathbf{S}^{\text{T}}(\mathbf{m}) (\mathbf{S}(\mathbf{m})\mathbf{S}^{\text{T}}(\mathbf{m}) + \mu\mathbf{I})^{-1} \delta\mathbf{d}^{*} = \mathbf{R}_{\text{d}}\delta\mathbf{d}^{*}, \quad (14) $$

where $\mathbf{S}^{\text{rd}}$ denotes the forward-modelling operator when the observation operator samples the wavefields at the datum level. Accounting accurately for this resolution operator requires performing the forward/backward simulations in a computational domain, the size of which is representative of the physical domain in which the

recorded data is propagated. These large-scale backward/forward simulations indeed make the method computationally expensive.

In this framework, improving the robustness and mitigating the computational burden of the underdetermined scattering source problem are compulsory, eq. (7). We regularize this source problem with a zero-order Tikhonov regularization in the frame of damped least-squares methods. Alternatively, sparsity-promoting regularization has been proposed to address source location problem (Shekar & Sethi 2019; Sharan et al. 2019; Aghamiry et al. 2021). Finally, we still need to improve the pre-conditioner of the MINRES solver that is used to perform the multidimensional deconvolution of the measured scattered data. Our pre-conditioner relies on a mono-dimensional filter in the frequency domain while it should account for the dimension of the data space. One possible approach to perform multidimensional filtering of the scattered data might be to combine multidimensional Gaussian filtering with least-squares projection on a B-spline basis (Operto et al. 2003, their appendix B).

## 6 CONCLUSIONS

We reformulate model-based seismic redatuming within the framework of least-squares time-domain wavefield reconstruction with data assimilation. Once the data-assimilated wavefields have been reconstructed, redatumed data can then be generated at any desired datum plane by sampling the reconstructed wavefields at the datum. The method can handle any kind of acquisition and wave types from above and beneath the datum level (including surface multiples). However, it requires a quite accurate knowledge of the subsurface between the surface and the datum plane in particular if this area generates reflections and multiscattering. The estimation of the data-assimilated wavefield requires performing a multidimensional deconvolution of the measured scattered data followed by a backward simulation of the deconvolved scattered data. This multidimensional deconvolution can be performed with several iterations of a Krylov-subspace linear solver. In this framework, a good pre-conditioner is a key issue to mitigate the computational burden of the method. Numerical tests and real data applications confirm the accuracy and effectiveness of the proposed method when the background medium is a good approximation of the true model above the datum plane. In this framework, obvious applications
are the redatuming of shots on the sea bottom and static corrections for marine and land applications, respectively. However, several key aspects still require improvement. Current work focuses on developing more efficient pre-conditioners for the data-domain normal system to reduce the number of iterations performed by the Krylov-subspace solver. Another challenge lies in the regularization process, where the use of zero-order Tikhonov regularization could be enhanced by adopting a more robust and effective approach. Additionally, future studies should further clarify the advantages and limitations of the proposed method in comparison to widely used redatuming techniques such as Marchenko redatuming.

## ACKNOWLEDGMENTS

We sincerely thank Editor Andrew Valentine, Assistant Editor and Reviewer Prof. Kees Wapenaar for their valuable comments. This study was funded by the WIND consortium (https://www.geoazur.fr/WIND), sponsored by AkerBP, ExxonMobil, Petrobras, Shell, and SINOPEC. The authors are grateful to the OPAL infrastructure from Observatoire de la Côte d'Azur (CRIMSON) for providing resources and support. This work was granted access to the HPC resources of IDRIS under the allocation A50050410596 made by GENCI. We also thank L. Combe (Geoazur) for her assistance during software development. Finally, we appreciate Chevron and DUG for providing a copy of the pre-processed Gorgon OBN data.

## DATA AVAILABILITY

The code package used to reproduce the results of this study is openly available at https://github.com/guogs93/Model-based-seism ic-data-redatuming-with-DAWs.git. The raw waveform data from the 2001 SFJ-OBS experiment used in this study is publicly available through JAMSTEC. However, the raw waveform data from the Gorgon gas field, provided by Chevron and DUG, is subject to restrictions and cannot be shared.

## REFERENCES

* Aghamiry, H.S., Gholami, A. & Operto, S., 2019. Improving full-waveform inversion by wavefield reconstruction with alternating direction method of multipliers, *Geophysics*, **84**(1), R139–R162.

* Aghamiry, H. S., Gholami, A., Operto, S. & Freitag, M., 2023. Efficient seismic redatuming using least-squares time-reversal, in *84th EAGE Annual Conference & Exhibition, Vol. 2023*, pp. 1–5, European Association of Geoscientists & Engineers.

* Aghamiry, H.S., Gholami, A. & Operto, S., 2020. Accurate and efficient data-assimilated wavefield reconstruction in the time domain, *Geophysics*, **85**(2), A7–A12.

* Aghamiry, H. S., Gholami, A., Operto, S. & Malcolm, A., 2021. Admm-based full-waveform inversion for microseismic imaging, *Geophys. J. Int.*, **228**(1), 259–274.

* Alerini, M., 2006. Stereotomography for nodes data, in *Extended Abstracts, 68nd Annual Meeting, (12–15 June 2006, Vienna)*, pp. B043, European Association of Geoscientists & Engineers.

* Alerini, M., Traub, B., Ravaut, C. & Duveneck, E., 2009. Prestack depth imaging of ocean-bottom node data, *Geophysics*, **74**(6), WCA57–WCA63.

* Aminzadeh, F., 1996. *3-D Salt and Overthrust Seismic Models, SEG/EAGE 3-D Modeling Series No.1*. American Association of Petroleum Geologists.

* Bakulin, A. & Calvert, R., 2004. Virtual source: new method for imaging and 4d below complex overburden, in *SEG Technical Program*

*Expanded Abstracts 2004*, pp. 2477–2480, Society of Exploration Geophysicists.

* Berryhill, J. R., 1979. Wave equation datuming, *Geophysics*, **44**, 1329–1333.

* Bevc, D., 1997. Flooding the topography: wave-equation datuming of land data with rugged acquisition topography, *Geophysics*, **62**(5), 1558–1569.

* Biondi, E., Barnier, G., Biondi, B. & Clapp, R.G., 2023. Target-oriented elastic full-waveform inversion through acoustic extended image-space redatuming, *Geophysics*, **88**(3), R269–R296.

* Borisov, D., Singh, S. & Fuji, N., 2015. An efficient method of 3-D elastic full waveform inversion using a finite-difference injection method for time-lapse imaging, *Geophys. J. Int.*, **202**(3), 1908–1922.

* Chin-Lung Fong, D. & Saunders, M., 2012. CG versus MINRES: an empirical comparison, *SQU J. Sci.*, **17**(1), 44–62.

* Choi, S.-C.T., Paige, C.C. & Saunders, M.A., 2011. MINRES-QLP: a Krylov subspace method for indefinite or singular symmetric systems, *SIAM J. Sci. Comput.*, **33**(4), 1810–1836.

* Cui, T., Vasconcelos, I., van Manen, D.-J. & Wapenaar, K., 2018. A tour of Marchenko redatuming: focusing the subsurface wavefield, *Leading Edge*, **37**(1), 67a1–67a6.

* Dessa, J.X., Operto, S., Kodaira, S., Nakanishi, A., Pascal, G., Uhira, K. & Kaneda, Y., 2004. Deep seismic imaging of the eastern Nankai trough (Japan) from multifold ocean bottom seismometer data by combined traveltime tomography and prestack depth migration, *J. geophys. Res.*, **109**(B02111). doi:10.1029/2004GL020453

* Di, H., Xie, W. & Xu, M., 2023. Upper crustal Vp/Vs ratios along the northern East Pacific Rise derived from downward-continued streamer data, *Geophys. J. Int.*, **235**(2), 1465–1478.

* Farshad, M. & Chauris, H., 2021. Accelerating the multi-parameter least-squares reverse time migration using an appropriate preconditioner, *Comput. Geosci.*, **25**(11), 2071–2092.

* Ferguson, R.J., 2006. Regularization and datuming of seismic data by weighted, damped least squares, *Geophysics*, **71**(5), U67–U76.

* Gholami, A., Aghamiry, H.S. & Operto, S., 2022. Extended full waveform inversion in the time domain by the augmented Lagrangian method, *Geophysics*, **87**(1), R63–R77.

* Górszczyk, A., Operto, S. & Malinowski, M., 2017. Toward a robust workflow for deep crustal imaging by FWI of OBS data: the eastern Nankai trough revisited, *J. Geophys. Res. Solid Earth*, **122**(6), 4601–4630.

* Gras, C., Dagnino, D., Jiménez-Tejero, C.E., Meléndez, A., Sallarès, V. & Ranero, C.R., 2019. Full-waveform inversion of short-offset, band-limited seismic data in the Alboran basin (SE Iberia), *Solid Earth*, **10**(6), 1833–1855.

* Guo, G., Operto, S., Gholami, A. & Aghamiry, H.S., 2022. A practical implementation of data-space hessian in the time-domain formulation of source extended full-waveform inversion, in *Second International Meeting for Applied Geoscience & Energy*, pp. 757–761, Society of Exploration Geophysicists and American Association of Petroleum.

* Guo, G., Operto, S., Gholami, A. & Aghamiry, H.S., 2024. Time-domain extended-source full-waveform inversion: algorithm and practical workflow, *Geophysics*, **89**(2), R73–R94.

* Hajjaj, R.F., de Ridder, S.A., Livermore, P.W. & Ravasi, M., 2022. Wavefield reconstruction inversion modelling of Marchenko focusing functions, preprint(arXiv:2210.14570).

* Hicks, G.J., 2002. Arbitrary source and receiver positioning in finite-difference schemes using Kaiser windowed sinc functions, *Geophysics*, **67**, 156–166.

* Hou, J. & Symes, W.W., 2015. An approximate inverse to the extended Born modelling operator, *Geophysics*, **80**(6), R331–R349.

* Huang, G., Nammour, R. & Symes, W.W., 2018. Source-independent extended waveform inversion based on space-time source extension: frequency-domain implementation, *Geophysics*, **83**(5), R449–R461.

* Huot, G. & Singh, S.C., 2018. Seismic evidence for fluid/gas beneath the mentawai fore-arc basin, central sumatra, *J. Geophys. Res. Solid Earth*, **123**(2), 957–976.

* Jimenez-Tejero, C.E., Ranero, C.R., Sallares, V. & Gras, C., 2022. Downward continuation of marine seismic reflection data: an under-valued tool to improve velocity models, *Geophys. J. Int.*, **230**(2), 831–848.
Jimenez-Tejero, C.E., Prada Dacasa, M., Gomez de la Peña, L., Ranero, C. & Sallares, V., 2025. Downward continuation of wide-angle seismic data: implications for traveltime tomography uncertainty, *Geophys. J. Int.*, ggaf137.

Li, L., Collier, J., Henstock, T. & Goes, S., 2024. Downward continued ocean bottom seismometer data show continued hydrothermal evolution of mature oceanic upper crust, *Geology*, **52**(9), 717–722.

Lomas, A. & Curtis, A., 2019. An introduction to Marchenko methods for imaging, *Geophysics*, **84**(2), F35–F45.

Luiken, N. & Van Leeuwen, T., 2020. Seismic wavefield redatuming with regularized multi-dimensional deconvolution, *Inverse Problems*, **36**(9), 095010. doi:10.1088/1361-6420/aba9d1

Mulder, W.A., 2005. Rigorous redatuming, *Geophys. J. Int.*, **161**(2), 401–415.

Operto, S., Lambaré, G., Podvin, P. & Thierry, P., 2003. 3-D ray-Born migration/inversion. Part 2: application to the SEG/EAGE overthrust experiment, *Geophysics*, **68**(4), 1357–1370.

Operto, S., Virieux, J., Dessa, J.X. & Pascal, G., 2006. Crustal imaging from multifold ocean bottom seismometers data by frequency-domain full-waveform tomography: application to the eastern Nankai trough, *J. geophys. Res.*, **111**(B09306), doi:10.1029/2005JB003835.

Operto, S. et al., 2023a. Is 3D frequency-domain FWI of full-azimuth/long-offset OBN data feasible? The Gorgon-data FWI case study, *The Leading Edge (special section Full Waveform Inversion)*, **42**(3), 173–183.

Operto, S. et al., 2023b. Extending the search space of FWI beyond the single-scattering born approximation, *Geophysics*, **88**(6), R671–R702.

Paige, C.C. & Saunders, M.A., 1982. LSQR: an algorithm for sparse linear equations and sparse least squares, *ACM Trans. Math. Softw.*, **8**(1), 43–71.

Planès, T., Snieder, R. & Singh, S., 2018. Model-based redatuming of seismic data: an inverse-filter approach, *Geophysics*, **83**(2), Q1–Q13.

Pratt, R.G., 1999. Seismic waveform inversion in the frequency domain, part I: theory and verification in a physical scale model, *Geophysics*, **64**, 888–901.

Pratt, R.G., Shin, C. & Hicks, G.J., 1998. Gauss–Newton and full Newton methods in frequency-space seismic waveform inversion, *Geophys. J. Int.*, **133**, 341–362.

Ravasi, M., 2017. Rayleigh-Marchenko redatuming for target-oriented, true-amplitude imaging, *Geophysics*, **82**(6), S439–S452.

Ravasi, M., Vasconcelos, I., Kritski, A., Curtis, A., Filho, C.A.d.C. & Meles, G.A., 2016. Target-oriented Marchenko imaging of a north sea field, *Geophys. J. Int.*, **205**(1), 99–104.

Sambolian, A., Gorszczyk, S., Operto, A., Ribodettì, B. & Tavakoli, F., 2021. Mitigating the ill-posedness of first-arrival traveltime tomography using slopes: application to the eastern Nankai Trough (Japan) OBS data set, *Geophysical Journal International*, **227**(2), 898–921.

Sambolian, S., 2021. Kinematically consistent slope tomography: theory and applications to velocity model building and event location, *Ph.D. thesis*, University Côte d’Azur.

Schuster, G.T. & Zhou, M., 2006. A theoretical overview of model-based and correlation-based redatuming methods, *Geophysics*, **71**(4), S1103–S1110.

Sharan, S., Wang, R. & Herrmann, F.J., 2019. Fast sparsity-promoting microseismic source estimation, *Geophys. J. Int.*, **216**, 164–181.

Shekar, B. & Sethi, H.S., 2019. Full-waveform inversion for microseismic events using sparsity constraints, *Geophysics*, **84**(2), KS1–KS12.

Shoja, A., van der Neut, J. & Wapenaar, K., 2024. Target-oriented least-squares reverse-time migration with Marchenko redatuming and double-focusing: field data application, *Geophysics*, **89**(3), S267–S274.

Symes, W.W., 2008. Approximate linearized inversion by optimal scaling of prestack depth migration, *Geophysics*, **73**(2), R23–R35.

Tournier, P.-H., Jolivet, P., Dolean, V., Aghamiry, H., Operto, S. & Riffo, S., 2022. Three-dimensional finite-difference & finite-element frequency-domain wave simulation with multi-level optimized additive Schwarz

domain-decomposition preconditioner: a tool for FWI of sparse node data sets, *Geophysics*, **87**(5), T381–T402.

Valenciano, A.A., Biondi, B.L. & Clapp, R.G., 2009. Imaging by target-oriented wave-equation inversion, *Geophysics*, **74**(6), WCA109–WCA120.

van Leeuwen, T. & Herrmann, F., 2016. A penalty method for PDE-constrained optimization in inverse problems, *Inverse Problems*, **32**(1), 1–26.

van Leeuwen, T. & Herrmann, F.J., 2013. Mitigating local minima in full-waveform inversion by expanding the search space, *Geophys. J. Int.*, **195**(1), 661–667.

Vargas, D., Vasconcelos, I., Ravasi, M. & Luiken, N., 2021. Time-domain multidimensional deconvolution: a physically reliable and stable preconditioned implementation, *Remote Sensing*, **13**(18), 3683.

Virieux, J., 1984. SH wave propagation in heterogeneous media: velocity-stress finite difference method, *Geophysics*, **49**, 1259–1266.

Wang, C., Yingst, D., Farmer, P. & Leveille, J., 2016. Full-waveform inversion with the reconstructed wavefield method, pp. 1237–1241, ed. Charles, Sicking, Society of Exploration Geophysicists.

Wapenaar, K., Draganov, D. & Robertsson, J.O., 2008. *Seismic Interferometry: History and Present Status*, Society of Exploration Geophysicists.

Wapenaar, K. et al., 2021. Marchenko redatuming, imaging, and multiple elimination and their mutual relations, *Geophysics*, **86**(5), WC117–WC140.

Yuan, S., Fuji, N., Singh, S. & Borisov, D., 2017. Localized time-lapse elastic waveform inversion using wavefield injection and extrapolation: 2-D parametric studies, *Geophys. J. Int.*, **209**(3), 1699–1717.

## APPENDIX A: IMPLEMENTATION OF THE ADJOINT-STATE EQUATION FOR SCATTERING SOURCES

In this appendix, we review the implementation of the adjoint-state equation for the scattering source $\delta\mathbf{b}^{(da)}$, eq. (8), when the state equation is the first-order velocity-stress wave equation.

### A1 Acoustic velocity-stress equation

We consider the following 2D velocity-stress acoustic wave equation:

$$ \left\{ \begin{aligned} \partial_t v_x(\mathbf{x}, t) &= bu(\mathbf{x}) \partial_x p(\mathbf{x}, t), \\ \partial_t v_z(\mathbf{x}, t) &= bu(\mathbf{x}) \partial_z p(\mathbf{x}, t), \\ \partial_t p(\mathbf{x}, t) &= \kappa(\mathbf{x}) (\partial_x v_x(\mathbf{x}, t) + \partial_z v_z(\mathbf{x}, t)) + b_p(\mathbf{x}, t), \end{aligned} \right. \quad (A1) $$

where $p(\mathbf{x}, t)$ is the pressure field, $v_x(\mathbf{x}, t)$ and $v_z(\mathbf{x}, t)$ are particle-velocity fields, $\kappa(\mathbf{x}) = \rho(\mathbf{x})c^2(\mathbf{x})$ is the bulk modulus, $c(\mathbf{x})$ is the compressional wavefield, $bu(\mathbf{x}) = 1/\rho(\mathbf{x})$ is the buoyancy, namely the inverse of the density $\rho$, $b_p(\mathbf{x}, t)$ is the pressure source, and $\mathbf{x} = (x, z)$ denotes the Cartesian coordinates.

The velocity-stress wave equation can be recast in matrix form as

$$ \partial_t \mathbf{u} = \mathbf{MD}\mathbf{u} + \mathbf{b}. \quad (A2) $$

$$ \text{with } \mathbf{u} = (v_x, v_z, p), \quad \mathbf{b} = (0, 0, b_p), \quad \mathbf{M} = \begin{pmatrix} \mathbf{B}_{2 \times 2} & \mathbf{0}_{2 \times 1} \\ \mathbf{0}_{1 \times 2} & \mathbf{C}_{1 \times 1} \end{pmatrix}, \quad \mathbf{B} = \begin{pmatrix} bu & 0 \\ 0 & bu \end{pmatrix}, \quad \mathbf{C}_{1 \times 1} = (\kappa), \quad \mathbf{D} = \begin{pmatrix} \mathbf{0}_{2 \times 2} & \mathbf{D}_{2 \times 1} \\ \mathbf{D}_{1 \times 2} & \mathbf{0}_{1 \times 1} \end{pmatrix}, \quad \mathbf{D}_{2 \times 1} = \begin{pmatrix} \partial_x \\ \partial_z \end{pmatrix}. $$

We need also to introduce initial conditions: $\mathbf{u}(0) = 0$; $\partial_t \mathbf{u}(0) = 0$.

In a more compact form,

$$ A(\mathbf{m})\mathbf{u} = \mathbf{b}, \quad (A3) $$

where $\mathbf{A} = \partial_t \mathbf{I} - \mathbf{MD}$ and $\mathbf{I}$ is the identity matrix.
## A2 Adjoint equation for first-order velocity-stress finite differences

To derive the adjoint of $\mathbf{A}$, it is convenient to define the self-adjoint (symmetric) operator $\mathbf{A}' = \mathbf{M}^{-1}\mathbf{A}$:

$$ \mathbf{A}' = \mathbf{M}^{-1}(\mathbf{m})\partial_t \mathbf{I} - \mathbf{D}, \quad (A4) $$

where $\mathbf{D}$ is symmetric unlike $\mathbf{MD}$. Let us introduce the inner product $\langle \cdot, \cdot \rangle_X$ over the subspace $\mathbf{X}$ of the state space $\mathcal{U} = \mathbf{X} \times \mathbf{T}$ where $\mathbf{X}$ and $\mathbf{T}$ denote the space and temporal dimensions spanned by $\mathcal{U}$:

$$ \int_{t=0}^{T} \langle \mathbf{v}(t), \mathbf{M}^{-1}(\mathbf{m})\partial_t \mathbf{u}(t) - \mathbf{D}\mathbf{u}(t) \rangle_X dt, \quad (A5) $$

To derive the adjoint operator, we perform an integration by part with respect to time and space:

* With respect to time

$$ \begin{aligned} \int_0^T \langle \mathbf{v}(t), \partial_t \mathbf{u}(t) \rangle_X &= [\langle \mathbf{v}(t), \mathbf{u}(t) \rangle]_0^T - \int_0^T \langle \partial_t \mathbf{v}(t), \mathbf{u}(t) \rangle_X \\ &= \langle \mathbf{v}(T), \mathbf{u}(T) \rangle_X - \int_0^T \langle \partial_t \mathbf{v}(t), \mathbf{u}(t) \rangle_X. \end{aligned} $$

* With respect to space

$$ \int \partial_x \mathbf{u}(t) \mathbf{v}(t) dx = - \int \mathbf{u}(t) \partial_x \mathbf{v}(t) dx, \quad (A6) $$

which means that $\mathbf{D}^T = -\mathbf{D}$.

Substituting these expressions in eq. (A5) gives

$$ \begin{aligned} \int_{t=0}^{T} \langle \mathbf{v}(t), (\mathbf{M}^{-1}(\mathbf{m})\partial_t - \mathbf{D}) \mathbf{u}(t) \rangle_X dt \\ = \int_{t=0}^{T} \langle (-\mathbf{M}^{-1}(\mathbf{m})\partial_t + \mathbf{D}) \mathbf{v}(t), \mathbf{u}(t) \rangle_X dt. \end{aligned} \quad (A7) $$

with a final condition $\mathbf{v}(T) = 0$.

Plugging the expression of this adjoint operator in the equation satisfied by $\delta\mathbf{b}^{(da)}$, eq. (8), gives

$$ - \mathbf{M}^{-1}(\mathbf{m}) \partial_t \delta \mathbf{b}^e(t) = - \mathbf{D} \delta \mathbf{b}^e(t) + \mathbf{M}^{-1}(\mathbf{m}) \mathbf{P}^T \delta \mathbf{d}^{(d)}(t) \quad (A8) $$

with $\delta\mathbf{b}^{(da)}(T) = 0$ and $\mathbf{P}^T$ is a prolongation operator which extends with zeroes the deconvolved scattered data in the full computational domain. For numerical solving with an explicit time integration scheme, we move back the subsurface properties on the right-hand side

$$ - \partial_t \delta \mathbf{b}^e(t) = - \mathbf{M}(\mathbf{m}) \mathbf{D} \delta \mathbf{b}^e(t) + \mathbf{P}^T \delta \mathbf{d}^{(d)}(t) \quad (A9) $$

We perform a change of variable $\delta\mathbf{b}'(t) = \delta\mathbf{b}^{(da)}(T-t)$ to transform the final condition into an initial condition:

$$ - \partial_t \delta \mathbf{b}'(T-t) = - \mathbf{M}(\mathbf{m}) \mathbf{D} \delta \mathbf{b}'(T-t) + \mathbf{P}^T \delta \mathbf{d}^*(t), \quad (A10) $$

followed by a second change of variable $t' = T - t$:

$$ \partial_{t'} \delta \mathbf{b}'(t') = - \mathbf{M}(\mathbf{m}) \delta \mathbf{D} \delta \mathbf{b}'(t') + \mathbf{P}^T \delta \mathbf{d}^*(T-t'). \quad (A11) $$

with $\delta\mathbf{b}^{(da)}(t') = 0$.

The differential operator $\mathbf{D}$ is the same as that used for the wave equation, eq. (A2). Therefore, the same finite-difference wave modelling code can used to compute the DAWs and the scattering

sources. Comparing the wave equation, eq. (A2), and the adjoint equation, eq. (A11) shows that one just needs to reverse in time the scattered data on the right-hand side and reverse the sign of the time step $dt$ during time integration to implement the adjoint equation satisfied by the scattering sources $\delta\mathbf{b}^{(da)}$.

## APPENDIX B: LIPPMANN-SCHWINGER EQUATION FOR SCATTERED-FIELD SIMULATION

In this appendix, we review the Lippmann-Schwinger equation when the wave-equation operator relies on the first-order velocity-stress formulation.

We start from the symmetric form of the wave equation (Appendix A):

$$ \mathbf{A}'(\mathbf{m})\mathbf{u} = \mathbf{M}^{-1}\mathbf{b}. \quad (B1) $$

For the sake of simplicity, we limit ourselves to the case where the density is constant and equal to 1 and $\mathbf{m}$ is the squared slowness:

$$ \mathbf{M}^{-1}(\mathbf{m}) = \begin{pmatrix} \mathbf{I}_{2 \times 2} & \mathbf{0}_{2 \times 1} \\ \mathbf{0}_{1 \times 2} & \mathbf{m} \end{pmatrix}, \quad (B2) $$

where $\mathbf{m} = 1/\mathbf{c}^2$. The radiation pattern matrix for the $i^{th}$ entry of $\mathbf{m}$ reduces to

$$ \frac{\partial \mathbf{A}'(\mathbf{m})}{\partial \mathbf{m}_i} = \begin{pmatrix} \mathbf{0}_{2 \times 2} & \mathbf{0}_{2 \times 1} \\ \mathbf{0}_{1 \times 2} & \mathbf{R}_i \partial_t \mathbf{I}_{1 \times 1} \end{pmatrix}, \quad (B3) $$

where $\mathbf{R}_i$ is a sampling or detection operator such that $(\mathbf{R}_i\mathbf{m})_j = \mathbf{m}_j$ if $i = j$ and 0 otherwise. Accordingly, the $i^{th}$ virtual (pressure) source is given by

$$ \mathbf{f}^{(i)} = \begin{pmatrix} \mathbf{0}_{2 \times 1} \\ \mathbf{R}_i \partial_t \mathbf{p}^* \end{pmatrix}. \quad (B4) $$

The scattered-field Lippmann-Schwinger equation is given by

$$ \mathbf{A}'(\mathbf{m})\delta\mathbf{u} = -\sum_{i=1}^{M} (\mathbf{0} \quad \mathbf{R}_i \delta\mathbf{m}_i)\mathbf{f}^{(i)}, \quad (B5) $$

where $\delta\mathbf{m}_i = 1/\mathbf{c}_i^2 - 1/\mathbf{c}_i^2$. Moving back the subsurface properties in the right-hand side gives

$$ (\partial_t - \mathbf{MD}) \delta \mathbf{u}^* = - \sum_{i=1}^{M} \begin{pmatrix} \mathbf{0}_{2 \times 2} & \mathbf{0}_{2 \times 1} \\ \mathbf{0}_{1 \times 2} & \mathbf{R}_i \left( \frac{\mathbf{e}^2}{\mathbf{e}^2 + 1} - 1 \right) \end{pmatrix} \partial_t \mathbf{u}^* = \delta \mathbf{b}^*. \quad (B6) $$

The above tedious notations mean that the pressure volume scattering source on the right-hand side is the sum of the impulse sources formed by the product between the time derivative of the true wave-field sampled at the position of the diffractor point $i$ and the contrast $\mathbf{c}_i^2/\mathbf{c}_i^{*2} - 1$.