# Localized time-lapse elastic waveform inversion using wavefield injection and extrapolation: 2-D parametric studies

Shihao Yuan, Nobuaki Fuji, Satish Singh, Dmitry Borisov

► To cite this version:

Shihao Yuan, Nobuaki Fuji, Satish Singh, Dmitry Borisov. Localized time-lapse elastic waveform inversion using wavefield injection and extrapolation: 2-D parametric studies. Geophysical Journal International, 2017, 209, pp.1699-1717. 10.1093/gji/ggx118. insu-03748855

HAL Id: insu-03748855

https://insu.hal.science/insu-03748855

Submitted on 10 Aug 2022

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
<figure><img src="image_2.png" /></figure>

# Localized time-lapse elastic waveform inversion using wavefield injection and extrapolation: 2-D parametric studies

Shihao Yuan,<fnref n="1" /> Nobuaki Fuji,<fnref n="1" /> Satish Singh<fnref n="1" /> and Dmitry Borisov<fnref n="2" />

Accepted 2017 March 20. Received 2017 March 6; in original form 2016 December 22

## SUMMARY

We present a methodology to invert seismic data for a localized area by combining source-side wavefield injection and receiver-side extrapolation method. Despite the high resolving power of seismic full waveform inversion, the computational cost for practical scale elastic or viscoelastic waveform inversion remains a heavy burden. This can be much more severe for time-lapse surveys, which require real-time seismic imaging on a daily or weekly basis. Besides, changes of the structure during time-lapse surveys are likely to occur in a small area rather than the whole region of seismic experiments, such as oil and gas reservoir or CO₂ injection wells. We thus propose an approach that allows to image effectively and quantitatively the localized structure changes far deep from both source and receiver arrays. In our method, we perform both forward and back propagation only inside the target region. First, we look for the equivalent source expression enclosing the region of interest by using the wavefield injection method. Second, we extrapolate wavefield from physical receivers located near the Earth’s surface or on the ocean bottom to an array of virtual receivers in the subsurface by using correlation-type representation theorem. In this study, we present various 2-D elastic numerical examples of the proposed method and quantitatively evaluate errors in obtained models, in comparison to those of conventional full-model inversions. The results show that the proposed localized waveform inversion is not only efficient and robust but also accurate even under the existence of errors in both initial models and observed data.

**Key words:** Seismic monitoring and test-ban treaty verification; Body waves; Seismic tomography.

## 1 INTRODUCTION

Tomographic studies seek best models, of which synthetic data are in concordance with observed data. In general, seismic data are recorded on the Earth’s surface. Amongst many proposed tomography methodologies in seismology, full waveform inversion has emerged as one of the most promising and powerful techniques since over three decades ago (Tarantola 1984; Woodhouse & Dziewonski 1984; Mora 1988; Pratt 1999; Virieux & Operto 2009; French & Romanowicz 2014). However, waveform inversion requires full wavefield modelling. It is thus computationally expensive to invert for a large area.

Recently, there are several attempts to focus wave-equation based seismic imaging methods on a small target area in the Earth rather than the entire region. Wang et al. (2016) iteratively inverted for a small region at the uppermost mantle in the vicinity of receiver arrays using teleseismic events and wavefield injection method for the source side. For the receiver side, wave extrapolation can mitigate the computational burden of waveform tomography. Dong et al. (2009) shows after wavefield extrapolation, the computational cost of reverse-time migration in a target region is significantly

reduced. Yang et al. (2012) shows localized inversion after double-side wavefield extrapolation is computationally more efficient, and converges faster compared to a standard full-model waveform inversion. Haffinger et al. (2013) indicated that due to the reduced computational domain after wavefield extrapolation, higher frequency contents could be introduced into the waveform inversion, which naturally improves the image quality. Baring the previous studies cited above for local seismic imaging in mind, in this study, we will (i) propose a methodology that allows to locally invert seismic waveform data for the region disconnected from both source and receiver arrays, combining source-side wavefield injection method and receiver-side residual wavefield extrapolation; and (ii) examine several 2-D synthetic tests to see the efficiency, robustness, accuracy and limitations of the proposed method.

During the oil and gas production in the subsurface, it is crucial to quantitatively estimate the model changes in the reservoir fluids. Reservoir monitoring requires a huge investment. In order to maximize the return on investment, one needs to take full advantage of a large volume of time-lapse seismic data. The method that we are proposing can be particularly appealing for time-lapse surveys. We have the ‘baseline’ data set that is acquired before the

<footnote n="1">¹Institut de Physique du Globe de Paris, 1 rue Jussieu, F-75238 Paris Cedex 05, France. E-mail: nobuaki@ipgp.fr</footnote>

<footnote n="2">²Department of Geosciences, Princeton University, Guyot Hall, Princeton, NJ 08544, USA</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** A schematic illustration of the proposed localized inversion in a layered model. Red star denotes the physical source. Grey and blue triangles denote physical and virtual receivers, respectively. The dashed square represents the absorbing boundary and the solid square represents the source injection boundary.</figcaption></figure>

production and the ‘monitor’ data set that is acquired on a regular basis (daily or weekly) during or after production. We are then to interpret data difference for model difference. Unfortunately, full-model waveform inversion for every monitor data set would be still too expensive. Since the production does not change the overall substructures, there should be much redundancy in conventional full-model inversions. Furthermore, the structural changes will mainly occur in the region of fluid removal or replacement, which are far from source and receiver arrays. As the model alterations are small and deep, the corresponding changes in seismic data go only a few per cents above the background noise. For the reasons outlined above, it is indispensable to develop an efficient and accurate waveform-based method to image the local structural changes.

Fig. 1 shows how our methodology works. We have a source (array) and a receiver array. Instead of performing inversion inside the whole model region, we would like to focus on a subvolume in $V^+$, where the model changes are expected to occur. All forward and backward simulations will be conducted within the reduced absorbing boundary (the dashed square shown in Fig. 1). Such a strategy of localized waveform inversion requires two levels of wavefield reconstructions: (i) representation of seismic sources in the subdomain (wavefield injection method) and (ii) extrapolation of residual wavefield between baseline and monitor surveys recorded at the physical receivers to virtual receivers in the subdomain (residual wavefield extrapolation). Borisov et al. (2015) proposed a localized waveform inversion technique for time-lapse surveys within a target region disconnected from source arrays. However, there are few studies that extended it to a region disconnected both from source and receiver arrays. As we will show in the later sections, receiverside wavefield extrapolation combining with source-side wavefield injection (Robertsson & Chapman 2000; Capdeville et al. 2003; Singh & Royle 2010; Monteiller et al. 2013; Masson et al. 2014; Borisov et al. 2015) can naturally filter out the non-physical extrapolated events caused by limited receiver aperture (the single-sided illumination problem). We formalize localized waveform inversion

problem for time-lapse seismics and present the two levels of wavefield reconstruction methodology. To demonstrate the feasibility and robustness of our methodology, we perform a series of synthetic tests with noise and uncertainties in both model and data domain.

## 2 LOCALIZED WAVEFORM INVERSION FOR TIME-LAPSE SURVEY

In this section, we formulate localized waveform inversion for time-lapse surveys. We first start with the classical $L_2$ norm misfit function of waveforms under the framework of Tarantola (1984) and Mora (1987):

$$E(\mathbf{m}) = \frac{1}{2} |\mathbf{d} - \mathbf{u}(\mathbf{m})|^2 \quad (1)$$

with $\mathbf{d}$ and $\mathbf{u}(\mathbf{m})$ being the gathered waveform data points of observed and synthetic computed for a model $\mathbf{m}$, respectively. By evaluating the first variation of the misfit function, we obtain a Newton method expression:

$$\mathbf{A}^{\mathrm{T}} \mathbf{A} \delta \mathbf{m}_{i} = \mathbf{A}^{\mathrm{T}} (\mathbf{d} - \mathbf{u}(\mathbf{m}_{i})) \quad (2)$$

and

$$\mathbf{m}_{i+1} = \mathbf{m}_i + \delta \mathbf{m}_i \quad (3)$$

where

$$\mathbf{A} = \frac{\partial \mathbf{u}}{\partial \mathbf{m}} \quad (4)$$

which is partial derivatives (or Fréchet derivatives) of the wavefield with respect to model parameters and $\mathbf{A}^{\mathrm{T}}\mathbf{A}$ is an approximated Hessian matrix. The right-hand side of eq. (2) is the gradient direction $\partial E(\mathbf{m})/\partial \mathbf{m}$ that can be calculated by either cross-correlation of partial derivatives and waveform difference between observed and synthetics or back propagation of waveform difference in the media $\mathbf{m}_i$. Multiple iterations (eq. 3) update model to overcome the nonlinearity of this inverse problem. The details of waveform inversion can be found in many literatures, some of which are cited herein.

We use above formulation to obtain the ‘baseline’ full model $\mathbf{m}^{\mathrm{base}}$ before the production mode. We then obtain ‘monitor’ data with some deformation in the subsurface media and we would like to invert for the model alteration from the ‘baseline’ model obtained during the first step. There are mainly three strategies to monitor the changes of target regions through inversions (Asnaashari et al. 2015): (i) ‘parallel difference strategy’ independently inverts time-lapse data sets, that is, the structure changes will be calculated through the subtraction of inverted baseline model from inverted monitor models; (ii) ‘sequential difference strategy’ will use the most recently acquired former monitor model as an initial model for inversions; (iii) ‘double-difference strategy’ inverts only time-lapse data difference for a model change instead of full-data inversions. Double-difference is widely applied in earthquake location and seismic imaging (Waldhauser & Ellsworth 2000; Zhang & Thurber 2003; Zheng et al. 2011; Zhang & Huang 2013; Asnaashari et al. 2015; Yang et al. 2015) and proves to be more robust and accurate in retrieving time-lapse variations and suppressing artefacts of inversions. Double-difference strategy requires strict repeatability of seismic surveys, which is generally satisfied by the state-of-the-art acquisitions, especially in time-lapse ocean-bottom seismics (Granger et al. 2005; Ronen et al. 2005). We thus seek to minimize the difference of observed and synthetic waveform
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Snapshots of vertical component velocity wavefield using full-model and localized forward simulation. Dashed square denotes absorbing boundary and solid square denotes wavefield injection boundary. Red star represents the physical source. (a–d) Simulations on an model with no perturbation inside the injection boundary. (e–l) A perturbed model. (b,f,j) Snapshots of full-model simulations. (c,g,k) Snapshots of localized forward simulations through wavefield injection. (d,h,l) Synthetic wavefield inside the red circle. Red curves and black curves denote the waveform obtained from full-model and localized simulations, respectively.</figcaption></figure>

residuals between baseline and monitor surveys:

$$ 
\begin{aligned}
E(\mathbf{m}^{\text{mon}}) &= \frac{1}{2} |\mathbf{d}^{\text{mon}} - \mathbf{u}(\mathbf{m}^{\text{mon}})|^2 \\
&\approx \frac{1}{2} |(\mathbf{d}^{\text{mon}} - \mathbf{u}(\mathbf{m}^{\text{mon}})) + (\mathbf{d}^{\text{base}} - \mathbf{u}(\mathbf{m}^{\text{base}}))|^2 \\
&= \frac{1}{2} |\Delta \mathbf{d} - \Delta \mathbf{u}(\mathbf{m}^{\text{mon}})|^2
\end{aligned}
\quad (5) $$

with superscripts ‘base’ and ‘mon’ denoting baseline and monitor, respectively. We define the waveform residual between baseline and monitor data as

$$ \Delta \mathbf{d} = \mathbf{d}^{\text{mon}} - \mathbf{d}^{\text{base}} \quad (6) $$

$$ \Delta \mathbf{u}(\mathbf{m}^{\text{mon}}) = \mathbf{u}(\mathbf{m}^{\text{mon}}) - \mathbf{u}(\mathbf{m}^{\text{base}}). \quad (7) $$

To minimize the new cost function defined above, we assume that the baseline model is close to the real structure and that the residuals between observed and synthetic waveforms are of the same order for baseline and monitor surveys. The Newton method for this new cost function is

$$ \mathbf{A}^{\text{T}} \mathbf{A} \delta \mathbf{m}_{i}^{\text{mon}} = \mathbf{A}^{\text{T}} (\Delta \mathbf{d} - \Delta \mathbf{u}(\mathbf{m}_{i}^{\text{mon}})). \quad (8) $$

Note that A does not change from that of eq. (4).

Since we are interested in a model perturbation $\delta \mathbf{m}^{\text{mon}}$ happening inside the region far deep from the source and receiver arrays, we reconstruct $\Delta \mathbf{d}$ and $\Delta \mathbf{u}$ inside the small region of interest by

wavefield injection method and residual wavefield extrapolation method that we discuss in the following sections.

Forward modelling is performed using 2-D staggered grid elastic isotropic finite-difference method with fourth-order accuracy in space and second-order accuracy in time (Virieux 1986; Levan-der 1988). We use the convolutional perfectly matched layers (Roden & Gedney 2000; Komatitsch & Martin 2007) to absorb undesirable seismic wave reflections from the model borders or outside the target regions. The gradient is obtained from adjoint-state method (Liu & Tromp 2006; Plessix 2006) and the conjugate gradient method is used to approximate the Hessian matrix (Pok-lak 1971). Compressional wave speed ($V_p$) and shear wave speed ($V_s$) are independently inverted, while $V_s$ ($V_s \approx V_p/1.67$) and the density model is linearly linked to $V_p$ (Gardner et al. 1974).

## 3 SOURCE-SIDE WAVEFIELD INJECTION METHOD

Wavefield injection method reconstructs a source expression on the surface of target regions, which is equivalent to a physical source distribution situated outside the region. Alterman & Karal (1968) introduced wavefield injection method, originally aiming at solving source singularity problems in finite-difference computations. Robertsson & Chapman (2000) expanded this method to the localized forward modelling, which proves to be an efficient way of calculating seismic wavefield after local model alterations. Opršal et al. (2009) discussed and generalized the associated boundary
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Wavefield extrapolation examples. (a,b) represent 2-D elastic layered $V_p$ model. Red stars denote source positions. Grey and blue triangles denote physical and extrapolated receivers, respectively. Panels (c) and (d) show the comparisons of observed (black curves) and extrapolated (red waveform) wavefield of selected two shots in (a) and (b), respectively.</figcaption></figure>

condition during wavefield injection. Borisov et al. (2015) applied it to waveform inversion for 3-D media with physical receivers inside the region.

Based on the injection boundary, the modelled region is divided into two subdomains, $V^+$ and $V^-$ (Fig. 1). The wavefield in either subdomain satisfies the same elastodynamic equation and can be expressed as

$$u^+(\mathbf{x}) = u_c(\mathbf{x}) = u_b(\mathbf{x}), \quad \mathbf{x} \in V^+, \qquad (9a)$$

$$u^-(\mathbf{x}) = u_s(\mathbf{x}) = u_c(\mathbf{x}) - u_b(\mathbf{x}), \quad \mathbf{x} \in V^-, \qquad (9b)$$

where $u$ denotes the displacement wavefield. The subscript + and - denote the wavefield inside and outside the injection boundary, respectively. The subscript c represents the complete wavefield and s represents the scattered wavefield caused by the model changes within injection boundary during time-lapse surveys. $u_b$ denotes the displacement wavefield excited by the stored boundary wavefield. The stress wavefield along the injection boundary satisfies the same condition as the displacement wavefield.

Even though $u^+$ in $V^+$ and $u^-$ in $V^-$ are achieved simultaneously, they are decoupled due to the boundary condition applied on injection boundary. Specifically, we add the stored boundary data for the wavefield inside the injection boundary, while we subtract the boundary data for that outside injection boundary. Since the physical sources are extrapolated from the surface to the source injection boundary during wavefield injection, we can naturally

decrease the size of absorbing boundary (the dashed square in Fig. 1) in order to reduce the inversion area. This may prohibit the interactions between the scattered wavefield and the structures outside the absorbing boundary. However, the missing part of the wavefield is only second or higher-order scattered wavefield, whose amplitude is weaker and arrivals are later than that of the first interactions. Also, it is numerically expedient considering the associated computational cost by applying full-model absorbing boundary or the so-called exact boundary condition (van Manen et al. 2007; Willemsen et al. 2016).

In order to better illustrate the role that wavefield injection plays as the localized modelling engine during time-lapse surveys, here we perform a synthetic experiment and compare it with the full-model simulations (Fig. 2). The 2-D baseline (a perturbation-free model within the injection boundary) is 8 km in the x-direction and 4 km in the z-direction, while the reduced modelling region is 2.8 km in the x-direction and 1.4 km in the z-direction. To store the wavefield along injection boundary, we are supposed to perform the full-model forward modelling once at the beginning for each physical source. During time-lapse surveys, we assume that perturbations occur only inside the injection boundary (Figs 2e and i). Fig. 2 shows the synthetic wavefield before and after model perturbations using full-model and localized simulations. The top panel in Fig. 2 indicates the wavefield before perturbations while the middle and bottom panels indicate the wavefield after perturbations. The top and middle panels display the wave propagation at earlier time steps while the bottom panel displays the later time step wavefield.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** $V_p$ and $V_s$ models for full-model and localized inversion tests. (a,b) True time-lapse $V_p$ and $V_s$ model. The dashed square encloses the target zone where we perform the localized inversion. (c,d) True baseline $V_p$ and $V_s$ models.</figcaption></figure>

<table><caption>**Table 1.** Comparisons of $V_p$ and $V_s$ model errors for all the inversion experiments.</caption><thead><tr><th>Test number</th><th>Initial model</th><th>Wavefield extrapolation</th><th>Noise</th><th>V<sub>p</sub> error</th><th>V<sub>s</sub> error</th></tr></thead><tbody><tr><td>#1(full)</td><td>True baseline</td><td>No</td><td>No</td><td>23.8 per cent</td><td>18.4 per cent</td></tr><tr><td>#2(localized)</td><td>True baseline</td><td>Yes</td><td>No</td><td>21.9 per cent</td><td>15.3 per cent</td></tr><tr><td>#3(full)</td><td>True baseline</td><td>No</td><td>No</td><td>66.5 per cent</td><td>66.9 per cent</td></tr><tr><td>#3(localized)</td><td>True baseline</td><td>Yes</td><td>No</td><td>22.6 per cent</td><td>29.1 per cent</td></tr><tr><td>#4(localized)</td><td>True baseline</td><td>Yes</td><td>Yes</td><td>24.3 per cent</td><td>16.4 per cent</td></tr><tr><td>#5(localized)</td><td>Perturbed baseline</td><td>Yes</td><td>No</td><td>22.2 per cent</td><td>17.1 per cent</td></tr><tr><td>#6(localized)</td><td>Inverted baseline</td><td>Yes</td><td>No</td><td>22.1 per cent</td><td>17.2 per cent</td></tr></tbody></table>

At earlier stage, the scattered wavefield caused by the perturbations have not propagated long enough to interact with medium outside the target area. Thus the seismograms produced by full-model and localized simulations are exactly the same (Figs 2d and h). However, at later times, the absorbing boundary in localized forward modelling will prevent the interactions between scattered wavefield and the medium outside the target region. Therefore, neglecting the second- or higher-order wavefield when applying wavefield injection method will cause small errors at later time steps, as observed in Fig. 2(l). However, we will show in Sections 5 and 6 that this imposes limited effect on the final inversion results through various synthetic experiments.

## 4 RECEIVER-SIDE RESIDUAL WAVEFIELD EXTRAPOLATION

Extrapolation of wavefield recorded at a physical (dense) receiver array is needed in our methods in order to perform localized waveform inversion. In this study, it is identical to what is applied commonly in reverse time migration (Baysal et al. 1983; Schuster 2002). We extrapolate residuals between baseline and monitor with a truncated surface integral due to the receiver array geometry in frequency

domain (Oristaglio 1989; Wapenaar & Fokkema 2006):

$$u_n^1(\mathbf{r}_A) \sim \int_{S'} \mathrm{d}S n_j [G_{in}^{0*}(\mathbf{r}, \mathbf{r}_A) C_{ijkl}^0 u_{k,l}^1(\mathbf{r}) - u_i^1(\mathbf{r}) C_{ijkl}^{0*} G_{kn,l}^{0*}(\mathbf{r}, \mathbf{r}_A)], \quad (10)$$

where $S'$ denotes the surface covered by the physical receiver array, which cannot always be a closed boundary (cf. Fig. A1). $u_n^1(\mathbf{r})$ is the $i$th component of perturbed displacement at $\mathbf{r}$, $n_j$ denotes the normal vector, $G_{in}^0(\mathbf{r}, \mathbf{r}_A)$ denotes the $i$th component of Green function from $n$th single force at $\mathbf{r}_A$ to $\mathbf{r}$ for baseline model, whereas $C_{ijkl}$ are elastic moduli for the baseline model. The derivation of the equation above can be found in Appendix A.

This leads to the so-called single-sided illumination problems which are commonly found in migration and seismic interferometry community (Snieder et al. 2006; Wapenaar & Fokkema 2006; Löer et al. 2014). The associated spurious events are particularly severe when the subsurface is strongly inhomogeneous and receiver apertures are very limited. In this study, we choose to apply the conventional extrapolation approach. Other extrapolation techniques, such as redatuming through multidimensional deconvolution (van der Neut et al. 2011; Wapenaar et al. 2011), Marchenko redatuming
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Full-model time-lapse inversion results. (a,b) Inverted time-lapse $V_p$ and $V_s$ model. (c,d) True time-lapse perturbation in the target zone, corresponding to the difference between time-lapse and baseline model within the dashed square in panels (a) and (b). (e,f) Inverted time-lapse perturbation of $V_p$ and $V_s$ in the target zone. The error in the title is obtained through eq. (11). Panels (g) and (h) correspond to the 1-D velocity profile along the dashed lines in panels (e) and (f), respectively. Black, cyan and red curves denote the true, starting and inverted velocity respectively. For simplicity, we use the same notification in the following experiments and we will not mention them repeatedly.</figcaption></figure>

(Wapenaar et al. 2014; Meles et al. 2016), will be interesting to test against our methods in the future work.

The Green's functions $G_{in}^{0*}(\mathbf{r}, \mathbf{r}_A)$ between physical and virtual receivers, serving as back-propagator in eq. (10) are calculated and stored just once using the baseline model obtained by the first

full-model waveform inversion. We generate them by posing single source forces at $\mathbf{r}_A$. The fact that we truncate the surface integral to $S'$ means that we neglect the interactions between perturbed backpropagators (caused by time-lapse model variations) and the recorded wavefield (Ravasi & Curtis 2013).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Localized time-lapse inversion results. (a,b) Inverted time-lapse perturbation of $V_p$ and $V_s$. Panels (c) and (d) correspond to the 1-D velocity profile along the dashed lines in panels (a) and (b), respectively. Panels (e) and (f) denote the initial and final residuals of localized inversion respectively.</figcaption></figure>

The evaluation of the integral in eq. (10) requires observed data of both displacement- (or particle velocity-) and stress-/strain-types. This is practically difficult since only displacement components (or velocity, acceleration) are recorded in most current land surveys but not stress nor strains. Even though replacing stress tensor by particle velocity can still guarantee the extrapolated wavefield being kinematically correct, it would result in additional spurious events. However, in marine surveys, stress tensor can be recorded by ocean-bottom 4C acquisitions (3C displacement/velocity plus 1C fluid pressure) since normal traction will be simplified to pressure due to fluid-solid boundary condition (Ravasi & Curtis 2013). Moreover, ocean-bottom seismics have additional advantages compared to conventional marine streamer acquisitions. Receivers fixed to the ocean bottom will not only improve repeatability for time-lapse surveys but also provide shear wave information, which leads to higher resolution due to shorter wavelength. Therefore, ocean-bottom 4-D acquisitions are particularly suitable for time-lapse monitoring and elastic wavefield extrapolation.

## 4.1 Wavefield extrapolation - 2-D numerical examples

The proposed extrapolation method was tested on a 2-D elastic layered model (P-wave velocity shown in Fig. 3a). Density is linked to $V_p$ through Gardner's rule (Gardner et al. 1974). In Fig. 3(a), grey triangles and blue triangles represent physical receivers and extrapolated receivers, respectively. Extrapolated depth is 1 km. The time-lapse data difference is caused by the ellipsoid perturbation enclosed in the solid black square (Fig. 3a). Green's functions are obtained from the baseline model without the ellipsoid perturbation. Comparisons of observed (black curves) and extrapolated (red wave-form) wavefield of selected two shots are shown in Figs 3(c) and (d). We can see that the extrapolated wavefield is generally accurate both in phase and amplitude. Some spurious events come out especially when shot position is at the edge of receiver plane due to loss of stationary phase points (Snieder et al. 2006).
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Single-shot full-model and localized inversion results. (a–d) Inverted time-lapse perturbation and corresponding 1-D velocity profile of single-shot full-model inversion. (e–h) Inverted time-lapse perturbation and corresponding 1-D velocity profile of single-shot localized inversion.</figcaption></figure>

# 5 LOCALIZED WAVEFORM INVERSIONS - 2-D NUMERICAL EXAMPLES WITH LAYERED MODELS

In this section, we present numerical examples in order to examine the efficiency, robustness and accuracy of the proposed localized waveform inversion. Conventional full-model inversion is tested against our method. The grid spacing for wave propagation is 20 m

and the time step is 2 ms. There are 31 point sources at the depth of 220 m excited in the vertical direction with a spacing interval of 200 m. Source time function is a Ricker wavelet of 6.5 Hz central frequency. 351 physical receivers are deployed at the depth 240 m with an interval of 20 m. We extrapolate those physical receivers to 56 virtual receivers at the depth of 1.32 km for the subsequent localized waveform inversion. The layered baseline model (Figs 4c
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Localized inversion with white noise. (a–d) Inverted time-lapse model and corresponding 1-D velocity profile along the dashed lines in panels (a) and (b).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** One shot example of extrapolated noisy wavefield (black curves) and observed wavefield (red waveform).</figcaption></figure>

and d) is 8 km in width and 4 km in depth. The time-lapse model has a lens-shaped velocity perturbation inside $V^+$ (Figs 4a and b), which is around 6 per cent of the baseline velocity for both $V_p$ and $V_s$.

In this section, we assume that the baseline model is identical to the ‘true’ baseline model except Sections 5.5 and 5.6, where we use the unknowingly perturbed and inverted baseline model, respectively. We define the model error associated with the inverted time-lapse perturbation within the target region as

Model error (in per cent)

$$ = \frac{\int_{V^+} dV \left| |\mathbf{m}_{\text{inv}}^{\text{mon}} - \mathbf{m}_{\text{ini}}^{\text{base}}| - |\mathbf{m}_{\text{true}}^{\text{mon}} - \mathbf{m}_{\text{true}}^{\text{base}}| \right|^2}{\int_{V^+} dV \left| |\mathbf{m}_{\text{true}}^{\text{mon}} - \mathbf{m}_{\text{true}}^{\text{base}}| \right|^2} \times 100 \text{ per cent} \quad (11) $$

with subscript ‘inv’ denotes inverted model, ‘ini’ the initial model and ‘true’ the true model. Note that this error definition is directly measuring the model perturbation and not model parameters themselves, resulting in the fact that the face values are deliberately large even for the best cases. There are still few localized waveform inversion studies where model error is defined: we would like the reader to look at the relative differences in Table 1 rather than to compare the face values with other inversion publications. The comparisons of each experiment are summarized in Table. 1, where the test number corresponds to the number included in the title of following subsections. Note that all the full-model and localized inversion results are obtained after the same iterations.

## 5.1 Full-model inversion (#1)

In the first experiment, we show the conventional full-model waveform inversion results starting from the true $V_p$ and $V_s$ baseline model (Figs 4c and d). We invert for the full time-lapse model shown in Figs 5(a) and (b), with the zoom-in of the target variation (Figs 5e and f). Besides the true perturbation (Figs 5c and d), we also find some inversion artefacts outside the target region as we can see in the 1-D velocity profile (Figs 5g and h). Due to the high ill-posedness and nonlinearity of the inverse problems, it can hardly ensure that the model update is solely focused on the target time-lapse perturbation.

The more unfavourable fact is that it is unavoidable to calculate the wavefield for the whole model every iteration, although the time-lapse perturbation occurs only in a certain small target region. Therefore, considering the computational cost, it is nearly impossible to conduct full-model waveform inversion on a daily or weekly basis for practical scale problems.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Localized inversion results with unknown baseline perturbation. (a,b) Time-lapse $V_p$ and $V_s$ model with unknown perturbation at shallow depth. (c–f) Inverted time-lapse perturbation and corresponding 1-D velocity profile along the dashed lines in panels (c) and (d).</figcaption></figure>

## 5.2 Localized inversion (#2)

The injection boundary data and the Green's functions are stored beforehand through a series of full-model forward modelling. With the source- and receiver-side extrapolation as we explained in Sections 3 and 4, the inversion area could be naturally reduced to a smaller region (the dashed square shown in Figs 4a and b) where time-lapse perturbation is expected to occur. Localized inversion provide us with better results for both $V_p$ and $V_s$ (Figs 6a and b) compared to the full-model inversion (Figs 5e and f). Since model update is solely allowed within the target region, the potential inversion artefacts can be mostly prevented. This is equivalent to a preconditioning of inverse problems which leads to a more accurate gradient direction. Also, localized inversion results look sharper compared to the full-model inversion due to the reduced offset and increased moveout, thus we obtain higher sensitivity to the small-scale perturbations. More importantly, the computing time in this test is reduced by a factor of 10 compared with the full-model inversion. The achieved computational efficiency is actually proportional

to the size of the non-calculated area outside the target regions. This shows the feasibility of time-lapse elastic full waveform inversion in an efficient way for practical scale problems.

## 5.3 Single-shot inversion (#3)

In conventional waveform inversions, long-offset data sets are crucial to retrieve the long wavelength information and mitigate cycle skipping problems (Shipp & Singh 2002). However, during time-lapse surveys, we expect that the long-wavelength structures (baseline model) remain almost unchanged and then the importance of long-offset can be reduced (Plessix et al. 2010). To invert small-scale velocity perturbations occurred locally in a target region, there is probably a redundancy when applying the same acquisitions for baseline surveys. To show the capability of the proposed localized inversion in reducing field cost (the number of shots and receivers), we perform a single shot inversion for both full-model and localized inversion. The results show that the full-model inversion (Figs 7a–d)
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Extrapolated wavefield (black curves) and modelled wavefield (red waveform) using the final inverted time-lapse model.</figcaption></figure>

fails in inverting the local perturbation while the localized inversion is still reliable (Figs 7e–h).

## 5.4 Localized inversion with white noise (#4)

In this test, we would like to show the robustness and accuracy of the localized inversion in presence of random noise. Artificial Gaussian white noise is added to the recorded data set before wavefield extrapolation. The signal-to-noise ratio of extrapolated wavefield is around 6 dB (black curves in Fig. 9). Even with strong white noise, the localized inversion is still quite robust and accurate as it shows in Fig. 8.

## 5.5 Localized inversion with unknown baseline perturbation (#5)

In previous tests, there is an assumption that the baseline model (especially the parts outside the target regions) is unchanged during the time-lapse surveys. However, it is possible that small perturbation occurs during drilling or oil production. In order to test the robustness of the localized inversion in the case of slight model changes outside the target region, we introduce an unknown perturbation outside the target zone during time-lapse surveys (Figs 10a and b). Note that wavefield injection and extrapolation are still performed using the original baseline model, which does not include the shallow perturbation. The inverted perturbations are still quite accurate as is shown in Figs 10(c)–(f).

Due to the inaccurate extrapolation propagators which do not contain the effect of the unknown shallow perturbation, some non-physical events will appear in the extrapolated wavefield (Fig. 11). The non-physical events result from the wrongly extrapolated shallow reflections caused by the unknown perturbations and come out before the first arrival at the virtual receivers. We can notice that the localized inversion does not fit the spurious events through comparing the modelled wavefield and extrapolated wavefield (black curves in Fig. 11). This is because wavefield injection naturally serves as anti-causal filters and exclude the effect of the short wavelength extrapolation error existing before the physical first arrivals. In other words, the cross-correlations between the spurious extrapolated events and partial derivatives are quasi-orthogonal (eq. 2). It is interesting to refer to Kawai *et al.* (2014) who performed 3-D localized waveform inversion for the base of the mantle using sources and receivers disconnected from the Earth's surface. They

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Extended baseline model inversion. (a,b) True baseline $V_p$ and $V_s$ model. (c,d) Starting $V_p$ and $V_s$ model obtained from the smoothed true baseline model. (e,f) Inverted baseline $V_p$ and $V_s$ model.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 13. 1-D $V_p$ and $V_s$ velocity profile along the dashed lines from inverted baseline model in Figs 12(e) and (f).</figcaption></figure>

also studied the risk of over-interpretation of data for the structure outside the inversion box but it turned out that the error is small due to quasi-orthogonality of partial derivatives for points inside and outside the box.

However, in case of severe perturbations, for instance, the associated first arrivals are apparently different from those of the unperturbed baseline model, we are supposed to do the full-model inversion once again to make sure the long wavelength information is accurate enough for wavefield injection and extrapolation. We might be able to control this by looking at the kinematics of time-lapse data difference. We can qualitatively conclude that slightly perturbed baseline model will not seriously affect the accuracy of the proposed localized inversion. In fact, as long as the baseline perturbations do not dramatically affect the first arrivals, the localized inversion will be still reliable and stable.

## 5.6 Localized inversion using inverted baseline model (#6)

To acknowledge the fact that in practice, no precise information on the interfaces and medium parameters is available, the boundary data for wavefield injection and the Green's functions for wavefield extrapolation should be modelled in a more realistic baseline model. Therefore, we divide inversions into two parts in this test: (i) baseline model inversion (ii) time-lapse model inversion. Due to the lack of illumination considering the sensitivity, we need to extend the survey area to retrieve the desired baseline model (Figs 12a and b). Starting from a smoothed model (Figs 12c and d), the baseline

models of $V_p$ and $V_s$ are inverted (Figs 12e and f) with associated 1-D velocity profiles shown in Figs 13(a)–(d).

With inverted baseline models, we repeat the same procedures as in previous localized inversion experiments. The inverted local perturbation (Figs 14a–d) seems still quite accurate as compared to the localized inversion using true baseline model in Section 5.2.

# 6 LOCALIZED WAVEFORM INVERSIONS - 2-D NUMERICAL EXAMPLES WITH MODIFIED MARMOUSI MODEL

We further validate our localized inversion strategy on the modified Marmousi model (Institut Franais du Pétrole (IFP), 1988) in comparison to full-model inversions. The $V_p$ model is 7.5 km × 2.5 km (Fig. 15a). $V_s$ ($V_s \approx V_p/1.5 + 300 \text{ m s}^{-1}$) and density models are linearly linked to the $V_p$ model (Gardner et al. 1974). We generate 4 s of data recorded at 300 m depth with a regular spacing of 25 m. In total 271 receivers record two components of velocity wavefield plus stress wavefield. A Ricker wavelet with dominant frequency of 6.5 Hz is used as the source wavelet and added onto the vertical velocity component. The grid spacing for the wave propagation is 25 m and the time step is 2 ms. 55 sources are deployed at 275 m depth with the interval of 125 m. The time-lapse velocity perturbations are about 5 per cent of the average baseline model (Fig. 15a).

The 2-D $V_p$ perturbations obtained from full-model and localized inversions under true baseline model are shown in Figs 15(c) and
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Localized inversion with inverted baseline model. (a,b) Inverted time-lapse perturbation of $V_p$ and $V_s$. (c,d) corresponds to the 1-D velocity profile along the dashed lines in (a) and (b). Panels (e) and (f) denote the initial and final residuals of localized inversion respectively.</figcaption></figure>

(d), respectively. Compared to the conventional full-model inversion, the localized inversion retrieves the velocity anomalies with a higher accuracy after the same 100 iterations. The corresponding 1-D velocity profiles are shown in Fig. 16. Fig. 17 shows the vertical velocity components of initial and final residuals from one shot gather. The residuals of the major events are better fitted after localized inversion (Fig. 17d) in comparison to the final full-model inversion residuals (Fig. 17b). We specially plot the theoretical initial residuals (Fig. 17e) in case that we directly record wavefield at the same positions as the extrapolated virtual receivers. The residual difference between Figs 17(c) and (e) is shown in Fig. 17(f). This difference essentially represents the spurious events after wavefield extrapolation. However, we notice that the localized inversion does not ‘see’ those spurious events when comparing it with final residuals of the localized inversion (Fig. 17d).

We finally apply the same strategy on a smoothed Marmousi model without changing the geometry and inversion parameters (Fig. 18). Besides the reduced computational cost, the localized in-

version is also more accurate than the full-model inversion (Fig. 19) after the same 100 iterations.

## 7 DISCUSSION

In order to develop a localized waveform inversion scheme, we combine source-side wavefield injection and receiver-side extrapolation methods, where we pose some simplifications and assumptions. The closed surface integral over the source-side wavefield injection can be conducted in a rigorous numerical way, while receiver-side residual wavefield extrapolation suffers from the truncation of surface integral and sparseness of receivers. Also, it requires a challenging technology that efficiently measures strains or stresses if we would like to rigorously complete all the terms in eq. (10). During the source-side wavefield injection method, we ignored higher interactions due to local model changes between target region and the model outside the reduced absorbing boundaries. However, it turns out that it does not seriously affect the seismic images since
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Time-lapse $V_p$ model and inverted results. (a) The baseline $V_p$ model selected from the true Marmousi model and the acquisition geometry. Red star and grey triangles denote physical source and receivers, respectively. White triangles denote the extrapolated virtual receivers. Two arrows point to the locations where time-lapse perturbation occurs. (b) The zooming true time-lapse $V_p$ perturbation within the dashed square in (a). (c) Full-model inversion result. (d) Localized inversion result. The two dashed lines (#1 and #2) in (c) and (d) indicate 1-D inverted velocity profiles and perturbations which are shown in Fig. 16.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16.** 1-D inverted $V_p$ profiles and perturbations. The left two columns denote 1-D profiles from the full-model inversion result in Fig. 15(c) and the right two columns denote the 1-D profiles from the localized inversion result in Fig. 15(d). (a–d) 1-D profiles with true (black curves), starting (cyan curves) and inverted $V_p$. (e–h) 1-D true (black curves) and inverted (red curves) $V_p$ perturbations.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Vertical component of inversion residuals from one shot example under true baseline Marmousi model. (a,b) Initial and final (after 100 iterations) residuals using the full-model inversion. (c,d) Initial and final (after 100 iterations) residuals using the localized inversion. (e) Initial residuals of localized inversion using the true wavefield directly observed at virtual receiver positions without extrapolation. (f) The residuals difference between (c) and (e), which indicates that the localized inversion will not try to fit the spurious events caused the wavefield extrapolation.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** (a) The smoothed baseline $V_p$ model with the same acquisition geometry and time-lapse perturbations shown in Figs 15(a) and (b). The two dashed lines (#1 and #2) indicate 1-D velocity profiles which are shown in (b) and (c). (b,c) 1-D profiles with true (black curves) and starting (cyan curves) velocities.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 19.** Full-model and localized inversion using the smoothed Marmousi baseline model. (a) The time-lapse $V_p$ perturbations from full-model inversion. The two dashed lines (#1 and #2) indicate 1-D inverted perturbations which are shown in (b) and (c). (b,c) 1-D true (black curves) and inverted (red curves) $V_p$ perturbations from full-model inversion. (d) The time-lapse $V_p$ perturbations from localized inversion. The two dashed lines (#1 and #2) indicate 1-D inverted perturbations which are shown in (e) and (f). (e,f) 1-D true (black curves) and inverted (red curves) $V_p$ perturbations from localized inversion.</figcaption></figure>

higher-order interactions are weaker and only influence later phases than that we use for inversions. Receiver-side residual wavefield extrapolation is, in addition to the difficulty described above, assuming no change (or error) in baseline model by definition. But it shows that long-wavelength structures of the outer region is sufficient to image the interior of target regions.

We show that wavefield injection is efficient and robust enough even though part of the higher order scatter wavefield is neglected due to the reduced modelling area. With the extrapolated wavefield, albeit some spurious events due to one-side illumination problem or slightly perturbed baseline model, the localized inversion is still able to provide sufficiently accurate time-lapse images. The localized inversion, which naturally focuses on the scattered wavefield due to perturbations within target regions, will not only contribute to the reduction of computational cost to a great extent, but also to the suppression of spurious events that cause wrong model update outside target zones. We further calculated the model errors as a function of iteration index in Fig. 20, showing that localized inversions converge faster in both cases (layered and Marmousi-models). Hence, the localized inversions require fewer iterations to obtain models with the same accuracy as full-model inversions. This will further enhance the speed-up and efficiency. Although we fixed the number of iterations to be 100 for all the tests in this paper, the truncation criteria should be discussed in our future papers. The decreased computational areas will potentially improve the capability in inverting higher frequency seismic data, which will lead to a better resolution, even though it is not shown here. Based on the synthetic examples we show, our methodology is more accurate, efficient and robust for the same number of iterations.

In data domain, after wavefield extrapolation, without losing the field interaction, the overburden effect is mostly reduced so that the nonlinearity of waveform inversion and associated cycle-skipping problems are mitigated. This can partly explain why localized inversions converge faster and the final misfit is smaller compared with conventional full-model inversion. In model domain, the localized inversion is equivalent to applying a precondition on the gradient, which naturally prevents the spurious model update from the non-target regions. This strategy is specially effective if we can specify our target zones and assume that the baseline model has few perturbations. In case of large baseline changes, we are supposed to perform the full-model waveform inversion once more to achieve better baseline model for wavefield injection and extrapolation. Although we have focused on the residual waveform inversion, the method proposed here is equally validated for localized inversion of seismic data in general.

## 8 CONCLUSIONS AND PERSPECTIVES

We proposed and examined a methodology of inverting seismic waveforms for localized structures and show its feasibility through various synthetic experiments by adding errors in both model and data domain. We combine wavefield injection and extrapolation to perform the localized waveform inversion for time-lapse seismic surveys. It shows that the localized inversion can enhance the image quality of local time-lapse variation while reducing the computational cost to a large extent. The proposed strategy reveals that putting virtual sources and receivers in the vicinity of target region contribute to improving the accuracy and robustness of waveform
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Comparisons of inverted errors as a function of iteration number between full-model and localized waveform inversions. (a,b) Using layered model (Sections 5.1 and 5.2). (c,d) Using modified Marmousi model (Section 6).</figcaption></figure>

inversion even under uncertainties of the outside area. It is as if we observed through a telescope to make the distant objects look closer and thus more clearly, even without knowing what is happening between them and us.

Due to the reduced modelling size, the computing time of 2-D localized inversion in the experiments above is decreased by a factor of 10 as compared to full-model inversions. We expect that the more severe computational cost of 3-D elastic or viscoelastic full waveform inversion during practical time-lapse surveys, could be potentially reduced to a more acceptable and economic level. Although the proposed methodology is examined on 2-D cases, the extension to general 3-D cases should be straightforward.

## ACKNOWLEDGEMENTS

The research described in this paper was carried out as a part of the Paris Exploration Geophysics Group project (GPX) funded by the French National Research Agency (ANR), CGG, TOTAL and Schlumberger. A significant part of the calculations in this work was performed on S-CAPAD at Institut de Physique du Globe de Paris. The authors thank two anonymous reviewers for fruitful comments.

## REFERENCES

* Alterman, Z. & Karal, F., 1968. Propagation of elastic waves in layered media by finite difference methods, *Bull. seism. Soc. Am.*, **58**(1), 367–398.

* Asnaashari, A., Brossier, R., Garambois, S., Audebert, F., Thore, P. & Virieux, J., 2015. Time-lapse seismic imaging using regularized full-waveform inversion with a prior model: which strategy?, *Geophys. Prospect.*, **63**(1), 78–98.

* Baysal, E., Kosloff, D.D. & Sherwood, J.W., 1983. Reverse time migration, *Geophysics*, **48**(11), 1514–1524.

* Borisov, D., Singh, S.C. & 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.

* Capdeville, Y., Chaljub, E. & Montagner, J., 2003. Coupling the spectral element method with a modal solution for elastic wave propagation in global earth models, *Geophys. J. Int.*, **152**(1), 34–67.

* Dong, S., Luo, Y., Xiao, X., Chávez-Pérez, S. & Schuster, G.T., 2009. Fast 3D target-oriented reverse-time datuming, *Geophysics*, **74**(6), WCA141–WCA151.

* French, S. & Romanowicz, B., 2014. Whole mantle radially anisotropic shear–velocity structure from spectral-element waveform tomography, *Geophys. J. Int.*, **199**, 1303–1327.

* Gardner, G., Gardner, L. & Gregory, A., 1974. Formation velocity and density; the diagnostic basics for stratigraphic traps, *Geophysics*, **39**(6), 770–780.

* Granger, P.-Y., Manin, M., Boelle, J.-L., Ceragioli, E., Lefeuvre, F. & Crouzy, E., 2005. Autonomous 4C nodes used in infill areas to complement streamer data, deepwater case study, in *2005 SEG Annual Meeting, Society of Exploration Geophysicists*.

* Haffinger, P., Gisolf, A. & Van den Berg, P., 2013. Towards high resolution quantitative subsurface models by full waveform inversion, *Geophys. J. Int.*, **193**(2), 788–797.

* Kawai, K., Konishi, K., Geller, R.J. & Fuji, N., 2014. Methods for inversion of body-wave waveforms for localized three-dimensional seismic structure and an application to D'' structure beneath Central America, **197**, 495–524.

* Komatitsch, D. & Martin, R., 2007. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation, *Geophysics*, **72**(5), SM155–SM167.

* Levander, A.R., 1988. Fourth-order finite-difference P–SV seismograms, *Geophysics*, **53**(11), 1425–1436.

* Liu, Q. & Tromp, J., 2006. Finite-frequency kernels based on adjoint methods, *Bull. seism. Soc. Am.*, **96**(6), 2383–2397.

* Löer, K., Meles, G.A., Curtis, A. & Vasconcelos, I., 2014. Diffracted and pseudo-physical waves from spatially limited arrays using source–receiver interferometry (SRI), *Geophys. J. Int.*, **196**(2), 1043–1059.
Masson, Y., Cupillard, P., Capdeville, Y. & Romanowicz, B., 2014. On the numerical implementation of time-reversal mirrors for tomographic imaging, *Geophys. J. Int.*, **196**, 1580–1599.

Meles, G.A., Wapenaar, K. & Curtis, A., 2016. Reconstructing the primary reflections in seismic data by Marchenko redatuming and convolutional interferometry, *Geophysics*, **81**(2), Q15–Q26.

Monteiller, V., Chevrot, S., Komatitsch, D. & Fuji, N., 2013. A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model, *Geophys. J. Int.*, **192**(1), 230–247.

Mora, P., 1987. Nonlinear two-dimensional elastic inversion of multiooffset seismic data, *Geophysics*, **52**(9), 1211–1228.

Mora, P., 1988. Elastic wave-field inversion of reflection and transmission data, *Geophysics*, **53**(6), 750–759.

Opršal, I., Matyska, C. & Irikura, K., 2009. The source-box wave propagation hybrid methods: general formulation and implementation, *Geophys. J. Int.*, **176**(2), 555–564.

Oristaglio, M.L., 1989. An inverse scattering formula that uses all the data, *Inverse Probl.*, **5**(6), 1097–1105.

Plessix, R.-E., 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, *Geophys. J. Int.*, **167**(2), 495–503.

Plessix, R.-E., Baeten, G., de Maag, J.W., Klaassen, M., Rujie, Z. & Zhifei, T., 2010. Application of acoustic full waveform inversion to a low-frequency large-offset land data set, in *2010 SEG Annual Meeting*, Society of Exploration Geophysicists.

Polak, E., 1971. *Computational Methods in Optimization: A Unified Approach*, vol. 77, Academic Press.

Pratt, R.G., 1999. Seismic waveform inversion in the frequency domain, Part 1: Theory and verification in a physical scale model, *Geophysics*, **64**(3), 888–901.

Ravasi, M. & Curtis, A., 2013. Elastic imaging with exact wavefield extrapolation for application to ocean-bottom 4C seismic data, *Geophysics*, **78**(6), S265–S284.

Robertsson, J.O. & Chapman, C.H., 2000. An efficient method for calculating finite-difference seismograms after model alterations, *Geophysics*, **65**(3), 907–918.

Roden, J.A. & Gedney, S.D., 2000. Convolutional PML (CPML): an efficient FDTF implementation of the CFS-PML for arbitrary media, *Microw. Opt. Technol. Lett.*, **27**(5), 334–338.

Ronen, S., Comeaux, L. & Miao, X.-G., 2005. Imaging downgoing waves from ocean bottom stations, in *2005 SEG Annual Meeting*, Society of Exploration Geophysicists.

Schuster, G.T., 2002. Reverse-time migration = generalized diffraction stack migration, in *SEG Technical Program Expanded Abstracts*, vol. 21, pp. 1280–1283.

Shipp, R.M. & Singh, S.C., 2002. Two-dimensional full wavefield inversion of wide-aperture marine seismic streamer data, *Geophys. J. Int.*, **151**(2), 325–344.

Singh, S. & Royle, G., 2010. Time-lapse elastic full waveform inversion using injected grid method, in *72nd EAGE Conference and Exhibition incorporating SPE EUROPEC 2010*, EAGE.

Snieder, R., Wapenaar, K. & Larner, K., 2006. Spurious multiples in seismic interferometry of primaries, *Geophysics*, **71**(4), S111–S124.

Tarantola, A., 1984. Inversion of seismic reflection data in the acoustic approximation, *Geophysics*, **49**(8), 1259–1266.

van der Neut, J., Thorbecke, J., Mehta, K., Slob, E. & Wapenaar, K., 2011. Controlled-source interferometric redatuming by crosscorrelation and multidimensional deconvolution in elastic media, *Geophysics*, **76**(4), SA63–SA76.

van Manen, D.-J., Robertsson, J.O. & Curtis, A., 2007. Exact wave field simulation for finite-volume scattering problems, *J. acoust. Soc. Am.*, **122**(4), EL115–EL121.

Virieux, J., 1986. *P–SV* wave propagation in heterogeneous media: velocity-stress finite-difference method, *Geophysics*, **51**(4), 889–901.

Virieux, J. & Operto, S., 2009. An overview of full-waveform inversion in exploration geophysics, *Geophysics*, **74**(6), WCC1–WCC26.

Waldhauser, F. & Ellsworth, W.L., 2000. A double-difference earthquake location algorithm: method and application to the northern hayward fault, *California, Bull. seism. Soc. Am.*, **90**(6), 1353–1368.

Wang, Y. et al., 2016. The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic P waves, *Geology*, **44**, 475–478.

Wapenaar, K. & Fokkema, J., 2006. Green’s function representations for seismic interferometry, *Geophysics*, **71**(4), SI33–SI46.

Wapenaar, K., Van Der Neut, J., Ruigrok, E., Draganov, D., Hunziker, J., Slob, E., Thorbecke, J. & Snieder, R., 2011. Seismic interferometry by crosscorrelation and by multidimensional deconvolution: a systematic comparison, *Geophys. J. Int.*, **185**(3), 1335–1364.

Wapenaar, K., Thorbecke, J., Van Der Neut, J., Broggini, F., Slob, E. & Snieder, R., 2014. Marchenko imaging, *Geophysics*, **79**(3), WA39–WA57.

Willemsen, B., Malcolm, A. & Lewis, W., 2016. A numerically exact local solver applied to salt boundary inversion in seismic full-waveform inversion, *Geophys. J. Int.*, **204**(3), 1703–1720.

Woodhouse, J.H. & Dziewonski, A.M., 1984. Mapping the upper mantle: three-dimensional modeling of earth structure by inversion of seismic waveforms, *J. geophys. Res.*, **89**, 5953–5986.

Yang, D., Zheng, Y., Fehler, M. & Malcolm, A., 2012. Target-oriented time-lapse waveform inversion using virtual survey, in *2012 SEG Annual Meeting*, Society of Exploration Geophysicists.

Yang, D., Meadows, M., Inderwiesen, P., Landa, J., Malcolm, A. & Fehler, M., 2015. Double-difference waveform inversion: feasibility and robustness study with pressure data, *Geophysics*, **80**(6), M129–M141.

Zhang, H. & Thurber, C.H., 2003. Double-difference tomography: the method and its application to the Hayward Fault, California, *Bull. seism. Soc. Am.*, **93**(5), 1875–1889.

Zhang, Z. & Huang, L., 2013. Double-difference elastic-waveform inversion with prior information for time-lapse monitoring, *Geophysics*, **78**(6), R259–R273.

Zheng, Y., Barton, P. & Singh, S., 2011. Strategies for elastic full waveform inversion of time-lapse ocean bottom cable (OBC) seismic data, in *2011 SEG Annual Meeting*, Society of Exploration Geophysicists.

## APPENDIX A: REPRESENTATION THEOREM FOR A SMALL MODEL CHANGE

We review representation theorem for reconstruction of residual wavefield extrapolation. Let be $V^+$ the region where we have a model change. Note that this definition of $V^+$ is not always as same as that in Section 3. Receiver-side residual wavefield extrapolation is reconstructing the residual between baseline and monitor data ($\Delta d$ in eq. 6) at some point just outside $V^+$ from $\Delta d$ recorded at the Earth’s surface (Fig. A1). As the representation theorem writes in a compact form in frequency domain, we first recall the equation of motion for an (an)elastic medium in frequency domain

$$ \rho^0 \omega^2 v_i^0 + (C_{ijkl}^0 v_{k,l}^0)_{,j} = -g_i \quad (A1) $$

where $\omega$ is angular frequency, $\rho$ is the density, $v_i$ is the ith component of displacement, $C_{ijkl}$ is the elastic moduli tensor and has complex value when the medium is anelastic, $g_i$ is the ith component of external body force distribution. Any internal or external boundary condition can be added in a form of surface integral to this equation but we omit the discussion on this topic in this paper. The comma $A_i$ denotes spatial partial derivative with respect to the direction. The superscript $0$ denotes the seismic structure of unperturbed model. We can write the equation of motion for the perturbed model in the same manner:

$$ (\rho^0 + \rho^1)\omega^2 (u_i^* + u_i^{1*}) + [(C_{ijkl}^{0*} + C_{ijkl}^{1*}) (u_{k,l}^* + u_{k,l}^{1*})]_{,j} = -f_i^* \quad (A2) $$

where $\rho^1$ and $C_{ijkl}^{1*}$ denotes the perturbation of seismic model with respect to the initial model, $u_i^*$ is the perturbed wavefield due to the model change with a source distribution $f_i$. $A^*$ denotes the
<figure><img src="image_23.png" /><figcaption>**Figure A1.** Schematic illustration of the geometry of representation theorem discussed in Appendix A. We would like to extrapolate (residual) wavefield recorded at the surface S' to the internal point r_A by using only Green's functions inside V^-.</figcaption></figure>

complex conjugate in frequency domain or time-reversal invariance in time-domain. Assuming the medium is lossless, we use complex conjugate type (or time reversal in time domain) equation of motion for eq. (A2) so that we can develop correlation type representation theorem for backward wavefield extrapolation. When needed, it is sufficient to take off all * to obtain the representation theorem of convolution type for forward wavefield extrapolation.

We then multiply eq. (A1) by $u_i^* = (u_i^{0*} + u_i^{1*})$ and eq. (A2) by $v_i^0$, respectively, taking the difference and taking a volume integral over an arbitrary V:

$$ -\int_V dV \rho^1 \omega^2 v_i^0 u_i^* + \int_V dV \left[ u_i^* (C_{ijkl}^0 v_{k,l}^0)_{,j} - v_i^0 (C_{ijkl}^{0*} + C_{ijkl}^{1*}) u_{k,l}^* \right]_{,j} = -\int_V dV (g_i u_i^* - f_i^* v_i^0). \quad (A3) $$

This is the generalized reciprocity theorem that takes account for a model change. Unlike the classical reciprocity theorem, the volume integral over $C_{ijkl}$ does not have anymore symmetric form due to the model change. Now we set the source distribution as

$$ g_i(\mathbf{r}) = \delta_{in} \delta(\mathbf{r} - \mathbf{r}_A) \quad (A4) $$

where $\delta_{in}$ denotes Kronecker's delta and $\delta(\cdot)$ denotes Dirac's delta function. Substituting eq. (A4) into eq. (A1), taking the volume integral over V, we obtain the equation of motion for Green's function:

$$ \int_V dV \rho^0 \omega^2 G_{in}^0(\mathbf{r}, \mathbf{r}_A) + \int_V dV (C_{ijkl}^0 G_{kn,l}^0(\mathbf{r}, \mathbf{r}_A))_{,j} = - \int_V dV \delta_{in} \delta(\mathbf{r} - \mathbf{r}_A). \quad (A5) $$

We then substitute eq. (A5) into eq. (A3), arranging it, we obtain a representation theorem for a model change:

$$ \begin{aligned} u_n^*(\mathbf{r}_A) = & \int_V dV f_i^* G_{in}^0(\mathbf{r}, \mathbf{r}_A) \\ & + \int_{V \cap V^-} dV \left[ G_{in}^0(\mathbf{r}, \mathbf{r}_A) (C_{ijkl}^{0*} u_{k,l}^*(\mathbf{r}))_{,j} - u_i^*(\mathbf{r}) (C_{ijkl}^0 G_{kn,l}^0(\mathbf{r}, \mathbf{r}_A))_{,j} \right] \\ & + \int_{V \cap V^+} dV \left[ \rho^1 \omega^2 u_i^*(\mathbf{r}) G_{in}^0(\mathbf{r}, \mathbf{r}_A) + u_i^*(\mathbf{r}) (C_{ijkl}^{1*} G_{kn,l}^0(\mathbf{r}, \mathbf{r}_A))_{,j} \right] \end{aligned} \quad (A6) $$

$$ \begin{aligned} & + \int_{V \cap V^+} dV \left[ G_{in}^0(\mathbf{r}, \mathbf{r}_A) (C_{ijkl}^{0*} + C_{ijkl}^{1*}) u_{k,l}^*(\mathbf{r}) \right]_{,j} \\ & - \int_{V \cap V^-} dV \left[ u_i^*(\mathbf{r}) (C_{ijkl}^0 + C_{ijkl}^{1*}) G_{kn,l}^0(\mathbf{r}, \mathbf{r}_A) \right]_{,j} \end{aligned} \quad (A6) $$

where, again, $V^+$ is the region where the model is perturbed, that is

$$ \rho^1 \equiv 0 \quad \text{in} \ V \cap V^- \quad (A7) $$

$$ C_{ijkl}^1 \equiv 0 \quad \text{in} \ V \cap V^- \quad (A8) $$

and

$$ V^- \equiv \overline{V^+}. \quad (A9) $$

When we substitute $\rho^1 \equiv 0$ and $C_{ijkl}^1 \equiv 0$ everywhere in $V^+$ for instance, we obtain an ordinary representation theorem of $u_n^0$ for the unperturbed model:

$$ u_n^{0*}(\mathbf{r}_A) = \int_V dV f_i^* G_{in}^0(\mathbf{r}, \mathbf{r}_A) + \int_V dV \left[ G_{in}^0(\mathbf{r}, \mathbf{r}_A) (C_{ijkl}^{0*} u_{k,l}^{0*}(\mathbf{r}))_{,j} - u_i^{0*}(\mathbf{r}) (C_{ijkl}^0 G_{kn,l}^0(\mathbf{r}, \mathbf{r}_A))_{,j} \right]. \quad (A10) $$

In this paper, we are particularly interested in representation for $\mathbf{r}_A \in V^-$ for residual wavefield extrapolation. We can thus drop off integral contributions from $V \cap V^+$ in eq. (A6). We then transform volume integrals to surface integrals by integrating by parts and using Gauss' theorem:

$$ \begin{aligned} u_n^*(\mathbf{r}_A) = & \int_{V^-} dV f_i^* G_{in}^0(\mathbf{r}, \mathbf{r}_A) \\ & + \int_{\partial V^-} dS n_j \left[ G_{in}^0(\mathbf{r}, \mathbf{r}_A) C_{ijkl}^{0*} u_{k,l}^*(\mathbf{r}) - u_i^*(\mathbf{r}) C_{ijkl}^0 G_{kn,l}^0(\mathbf{r}, \mathbf{r}_A) \right] \end{aligned} \quad (A11) $$

where $n_j$ denotes the normal vector and $\partial V^-$ is a closed surface surrounding the volume $V^-$. In our case, the source is not inside the perturbation region $V^+$, we have

$$ \int_{V^-} dV f_i^* G_{in}^0(\mathbf{r}, \mathbf{r}_A) = u_n^{0*}(\mathbf{r}_A) \quad (A12) $$

which leads to

$$ u_n^{1*}(\mathbf{r}_A) = \int_{\partial V^-} dS n_j \left[ G_{in}^0(\mathbf{r}, \mathbf{r}_A) C_{ijkl}^{0*} u_{k,l}^{1*}(\mathbf{r}) - u_i^{1*}(\mathbf{r}) C_{ijkl}^0 G_{kn,l}^0(\mathbf{r}, \mathbf{r}_A) \right]. \quad (A13) $$

$\partial V^-$ has the two parts: (i) the outer surface including the Earth's surface and (ii) the internal surface of which the normal vector is opposite to that of $\partial V^+$. The latter surface integral can be then again translated into the volume integral over $V^+$ but from eq. (A3), this term will vanish if the density perturbation $\rho^1$ is 0. We thus only have to take the surface integral over the outer surface of $V^-$. When taking $V^-$ as the Earth itself except $V^+$, with a smaller contribution from further regions from the source, physical and virtual receivers. In the reality, the surface integral over the closed boundary (eq. A13) can never be realized due to the acquisition geometry and we have to truncate $\partial V^-$ to $S'$ as in Fig. A1. The formal analysis of the errors of this truncation is difficult. In Section 4, we show some numerical 'good' and 'bad' examples of extrapolated waveforms. Despite the significant errors in extrapolated waveforms as in Fig. 3, waveform inversion scheme is filtering out those spurious events and give robust models.