# Three-dimensional elastic full waveform inversion in a marine environment using multicomponent ocean-bottom cables: a synthetic study

Dmitry Borisov, Satish C. Singh

► To cite this version:

Dmitry Borisov, Satish C. Singh. Three-dimensional elastic full waveform inversion in a marine environment using multicomponent ocean-bottom cables: a synthetic study. Geophysical Journal International, 2015, 201, pp.1215-1234. 10.1093/gji/ggv048. insu-03579970

HAL Id: insu-03579970

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

Submitted on 18 Feb 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>

# Three-dimensional elastic full waveform inversion in a marine environment using multicomponent ocean-bottom cables: a synthetic study

Dmitry Borisov and Satish C. Singh

Institut de Physique du Globe de Paris (IPGP), 1 Rue Jussieu, F-75005 Paris, France. E-mail: borisov@ipgp.fr

Accepted 2015 January 28. Received 2015 January 27; in original form 2014 April 12

## SUMMARY

Full waveform inversion (FWI) is a powerful tool used to quantify the elastic properties of the subsurface from seismic data. Because of very high computational cost, the technique has so far been used for either 2-D full elastic or 3-D acoustic media while the extension to 3-D elastic media to a realistic model size is still a challenging task. However, the Earth being 3-D, elastic and highly heterogeneous, one would require a full 3-D elastic wave equation for accurate modelling of amplitudes and phases within the inversion process. The acoustic approximation could significantly impact the final waveform inversion results, mainly due to the amplitude variation with offset effect. This effect becomes extremely important in the presence of strong contrasts in S-wave velocity and density, specifically when long-offset reflection data are included for waveform inversion. Recent increase in computer power allows for more efficient parallel computing using thousands of processors simultaneously thus making 3-D elastic waveform inversion feasible today. In this paper we consider a synthetic study based on a 3-D elastic medium for inversion of both *P*- and *S*-wave velocities using multicomponent, ocean-bottom cable seismic data. Both the forward modelling part and the inversion part are carried out in the time domain. The inverse problem is parametrized in terms of *P*- and *S*-wave velocities, while the density, being difficult to reconstruct, is not inverted and is linked to the *P*-wave velocity. Among several synthetic examples, a successful experiment on a small part of a 3-D SEG/EAGE overthrust model is presented, demonstrating the feasibility of inverting and accurately quantifying both *P*- and *S*-wave velocities. The resolution analysis of the waveform inversion is tested using a checkerboard model. Our results show that 3-D elastic FWI of sparsely spaced sources can retrieve *P*- and *S*-wave velocities accurately.

**Key words:** Inverse theory; Computational seismology; Wave propagation.

## 1 INTRODUCTION

Full waveform inversion (FWI) is a powerful tool used to determine quantitative images of the earth subsurface from seismic data. The technique was introduced as a local optimization problem (Lailly 1983; Tarantola 1984) where wave propagation was carried out using a finite-difference (FD) method followed by least-squares minimization of the misfit between observed and modelled data using an adjoint state technique (e.g. Plessix 2006).

There have been many debates concerning the appropriate level of physics incorporated within FWI process and the question is still open. Generally *P*-wave velocity is inverted using either acoustic or elastic approximations. Inversion of additional parameters, such as *S*-wave velocity, density, attenuation, anisotropy, could lead to increasing ill-posedness of the inverse problem due to the increase in the degrees of freedom and the difference in sensitivity

to different parameters (Virieux & Operto 2009). Moreover, the multiparameter inversion could lead to a dramatic increase in computational requirements, especially for a 3-D medium. Therefore, it is currently widely accepted to apply the acoustic approximation for seismic modelling. This leads to a single-parameter inversion, that is the *P*-wave velocity (*Vp*) only, assuming the density is constant. A few authors (Vigh et al. 2010; Warner et al. 2012) have shown that the acoustic wave equation is sufficient for FWI while inverting mainly turning rays, provided the purpose of FWI is to generate *Vp* models for depth migration. Warner et al. (2012) also suggest that attenuation and anisotropy are much more important than elasticity. However, waveform inversion accounts for the full information contained in the seismic data, and the integrated modelling engine should theoretically honour as much as possible the physics of wave propagation. Given the fact that the Earth is 3-D, elastic and highly heterogeneous, a full 3-D elastic wave equation is required for accurate modelling and for preserving correct amplitudes within the
inversion process, particularly when combined near and far offset data are used.

The acoustic approximation cannot be used widely, particularly because of the amplitude variation with offset (AVO) effects in the presence of strong contrast in S-wave velocity (Vs) and density (Barnes & Charara 2009). The AVO effect becomes very important at large offsets that are also used in retrieving medium scale-length features of velocity model within the FWI process (Shipp & Singh 2002). Arnulf et al. (2014) show that the acoustic approximation for an elastic medium can lead to an error of ±500 m s⁻¹ while inverting only turning rays because of erroneous partitioning of energy at the seafloor. On the other hand, the use of elastic FWI could provide valuable high-resolution estimates of the Vp/Vs ratios and density (Shipp & Singh 2002). These elastic parameters could act as a lithology indicator (a low Vp/Vs ratio often indicates the presence of sands) and could give valuable information about fluid content, for example help distinguish hydrocarbon from water (Sears et al. 2008; Lu et al. 2013).

Because seismic modelling in a 3-D elastic medium is computationally quite expensive, elastic FWI has been used mainly for 2-D media in marine (Sears et al. 2008, 2010) and land (Brossier et al. 2009) environments. However, the effect of using a 2-D approximation in modelling could negatively impact the final inversion results because a simple 3-D to 2-D correction is generally applied to take 3-D wave propagation from a point source into account but that in the presence of complex geological structures such a correction is not valid. For 3-D media, an acoustic approximation was applied to synthetic data sets (Ben-Hadj-Ali et al. 2008). It has also been used in real data sets to invert long offset refracted arrivals for marine data sets (Plessix 2009; Sirgue et al. 2010) and in land environments (Plessix et al. 2012). Since computer power has increased by three orders of magnitude in terms of FIOPS (Floating-point Operations Per Second) in the last decade, and assuming that this rate keeps up in the future, one would expect 3-D elastic FWI to be widely implemented during the upcoming decade.

In this study, we present our implementation of 3-D FWI for an elastic isotropic medium in the time distance domain. We first briefly discuss the theoretical and numerical aspects of wave propagation and an FWI scheme formulated in the time domain for a 3-D elastic medium. A key aspect in the inversion algorithm is an efficient forward modelling engine, which is achieved using high-performance parallel computation on a distributed-memory platform. We then present several synthetic examples from a small part of a 3-D SEG/EAGE overthrust model in a marine environment using multicomponent OBC (ocean bottom cable) data. The importance of including elasticity as well as the effect of 3-D wave propagation within waveform inversion is demonstrated. We also present a successful application of the technique to image shallow structures. Finally, a checkerboard test is performed to investigate the resolution characteristics of the inverted P- and S-wave velocity models as a function of depth.

## 2 3-D ELASTIC FWI

Waveform inversion provides high-resolution results that have a large potential for different applications. Compared to conventional seismic processing and imaging methods, waveform inversion is a data driven approach. Therefore, there is no need to identify and pick the events, as is required for travel time tomography, for instance. FWI allows to take into account the full information contained in seismic data, exploiting both travel times and waveform attributes

such as amplitude, phase and frequency content. All types of waves could be theoretically involved in the optimization process, including pre- and post-critical reflections, diving waves and multiples. Currently the results of FWI are mainly used as an improved input velocity model for pre-stack depth migration. However it can also be used as a quantitative seismic image for geological interpretation (Jaiswal et al. 2008; Sears et al. 2010; Arnulf et al. 2012), lithology estimation, reservoir monitoring (Queißer & Singh 2013) and pore pressure estimation (Roberts et al. 2008).

The FWI algorithm implemented in this study is based on the theoretical framework of Tarantola (1986) and Mora (1987). The waveform inversion algorithm can be divided in two major parts, summarized briefly in the following subsections: the forward problem (modelling) and the inverse problem.

### 2.1 Forward problem

The forward problem represents an important part within a waveform inversion scheme because most of the computing time is generally spent in forward modelling. Successfulness of the inversion strongly depends on the balance between accuracy, efficiency and the level of physics included in the forward modelling solver.

The forward problem of a seismic wave-field in a 3-D isotropic elastic medium involves simulating wave propagation through the earth and can be described by the following set of first-order equations in the velocity-stress formulation:

$$ 
\begin{aligned}
\partial_t v_x &= 1/\rho (\partial_x \sigma_{xx} + \partial_y \sigma_{xy} + \partial_z \sigma_{xz} + f_x), \\
\partial_t v_y &= 1/\rho (\partial_x \sigma_{xy} + \partial_y \sigma_{yy} + \partial_z \sigma_{yz} + f_y), \\
\partial_t v_z &= 1/\rho (\partial_x \sigma_{xz} + \partial_y \sigma_{yz} + \partial_z \sigma_{zz} + f_z),
\end{aligned} 
\quad (1) $$

$$ 
\begin{aligned}
\partial_t \sigma_{xx} &= (\lambda + 2\mu)\partial_x v_x + (\partial_y v_y + \partial_z v_z), \\
\partial_t \sigma_{yy} &= (\lambda + 2\mu)\partial_y v_y + (\partial_x v_x + \partial_z v_z), \\
\partial_t \sigma_{zz} &= (\lambda + 2\mu)\partial_z v_z + (\partial_x v_x + \partial_y v_y), \\
\partial_t \sigma_{xy} &= \mu (\partial_y v_x + \partial_x v_y), \\
\partial_t \sigma_{yz} &= \mu (\partial_z v_y + \partial_y v_z), \\
\partial_t \sigma_{xz} &= \mu (\partial_z v_x + \partial_x v_z),
\end{aligned} 
\quad (2) $$

where v is the velocity vector, σ is the stress tensor, ρ is the density, λ and μ are Lamé parameters, {f_s, f_z, f_y} are the body force components, and the symbols ∂_t, ∂_x, ∂_y and ∂_z are short representations of the derivatives ∂/∂_t, ∂/∂_x, ∂/∂_y and ∂/∂_z, respectively. This coupled system of equations gives three equations for particle velocities {v_x, v_y, v_z}, three equations for normal stresses {σ_{xx}, σ_{yy}, σ_{zz}} and three equations for shear stresses {σ_{xz}, σ_{xy}, σ_{yz}} wavefield components. A FD solution using a staggered grid scheme, second order in time and fourth order in space [O(2, 4)], introduced by Levander (1988), allows for an accurate modelling of P and S waves in the presence of a fluid-solid interface, such as the seafloor. The stability condition is maintained using the relationship between time step (Δt) and grid spacing (Δx) (Yao & Margrave 2000)

$$ \Delta t \le \frac{6\Delta x}{7\sqrt{3}Vp_{\max}}, \quad (3) $$

where Vp_{max} is a maximum value of P-wave velocity. Since both P and S waves are generated and propagate in the media, the maximum frequency (f_{max}) that can be accurately modelled is limited by a
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic illustration of the spatial 3-D domain decomposition technique used for efficient parallel computing. Initial 3-D volume (left-hand panel) divided into many blocks (right-hand panel). Each computer processor performs calculations within small blocks separately and then interchanges results with neighbourhood processors.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Efficiency of our algorithm for wave propagation in a 3-D elastic medium as function of number of cores as compared to a linear increase in number of cores. The model dimensions are 12.8 × 3.2 × 12.8 km in X, Y and Z, respectively, with 20 m for grid spacing. Black bold line is computing time for different number of cores; dashed line is a linear relationship (100 per cent efficiency).</figcaption></figure>

minimum of 5 gridpoints per shortest wavelength and depends on the minimum shear velocity ($Vs_{min}$)

$$f_{\max} < \frac{Vs_{min}}{5\Delta x}. \qquad (4)$$

The passage from 2-D to 3-D modelling within FWI obviously affects the total running time significantly due to the increased number of gridpoints in the second horizontal direction of the model and the increased number of source simulations. The difference could be as big as ~3–5 orders of magnitude. One can reduce the computational cost using an acoustic approximation (Plessix 2009; Sirgue et al. 2010), but as discussed above, the use of an acoustic approximation within FWI could hamper the algorithm to converge to the local minimum and in some cases could lead to erroneous

results (Barnes & Charara 2009; Solano et al. 2013; Arnulf et al. 2014; Borisov et al. 2014).

The increase in computational cost of 3-D modelling in elastic media in comparison to acoustic modelling mainly comes from three factors. (1) First, in the acoustic medium one of the Lamé parameters ($μ$) is equal to zero, which leads to neglecting all three shear stress components. The acoustic formulation can be further simplified replacing the three remaining normal stress components by single pressure component. These simplifications reduce the total number of equations required to be solved from nine for elastic to four for acoustic. (2) Secondly, a bigger grid spacing can be used for acoustic modelling for the same frequency band, because $Vp_{min}$ in eq. (4) is used and $Vp_{min} > 1.42 \times Vs_{min}$. Therefore, the total number of points within a 3-D grid would be significantly smaller for an acoustic medium as compared to that for an elastic medium. (3) Finally, the stability condition (eq. 3) also allows larger intervals for time discretization, thus reducing the total number of time steps. In total, the numerical modelling in a 3-D elastic medium as compared to acoustic media, could be 1–3 orders of magnitude computationally more expensive, depending of course on the experiment and the particular model characteristics (e.g. maximum and minimum velocity values of P and S waves). Therefore, a significant increase in computing resources is required to perform numerous forward modelling within an elastic waveform inversion, which makes 3-D elastic FWI for seismic imaging a real challenge.

The implementation of non-reflecting artificial boundaries plays an important role in seismic modelling. Unsplit convolutional perfectly matched layers (C-PML) are used in our study to absorb undesirable reflections from model boundaries (Komatitsch & Martin 2007). It allows for an accurate wavefield absorption using layers of only 10-gridpoint in thickness, even at grazing incidence angles, that is when wave is propagating nearly parallel to the absorbing boundary. This occurs when a seismic source is placed close to the model boundary where the classical PML method based on wavefield splitting (Collino & Tsogka 2001) produces erroneous reflection of large amplitude (Komatitsch & Martin 2007). Moreover, this type of boundary conditions has proven to be efficient for the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** P-wave velocity from the SEG/EAGE overthrust model, modified by adding water layer of 200 m on top of the model: 3-D view (a), vertical slices through the middle of the model (b), and horizontal slices at 1.8 km depth (c). Inline and cross-line correspond to *X* and *Y* directions.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Wavefield computed using a 5-Hz Ricker source wavelet after 1.6 s (a) and 4.0 s (b) of the source injection in the model. The source position is indicated as the yellow star, OBC-receivers are marked by white dots.</figcaption></figure>

elastic wave equation from a numerical point of view. For simplification, in most examples presented below, the boundary conditions were implemented at all model edges including the top surface, and thus multiples were not taken into account in this study. For a more realistic modelling, a free surface boundary condition should be considered at the top of the grid. Robertsson (1996) compared various numerical implementations of the free surface conditions for FD modelling and proposed a method, which is relatively simple to implement.

Although wavefield modelling in 3-D elastic media using the FD method was formulated in the 20th century (Randall 1989; Graves 1996), one would require extensive parallelization for the FWI. Our algorithm is heavily parallelized on a distributed memory system

using MPI (message passing interface) and a domain decomposition technique, allowing for efficient and accurate calculation of the 3-D wavefield (Minkoff 2002; Toxopeus et al. 2002). For parallelization the whole 3-D volume is decomposed into sub-volumes in both horizontal (*X*, *Y*) and vertical (*Z*) directions, constituting $N_x \times N_y \times N_z$ cubes of equal size. At each time step, each processor independently updates the wavefield within its portion of the grid and exchanges the wavefield information at the interface with the neighbourhood subvolume (Fig. 1). We calculated the efficiency of the algorithm using a fixed size problem by evaluating the time of wave propagation within a fixed 3-D volume for an increasing number of used computer processors. In Fig. 2 the bold line represents the computing time for different numbers of cores, while the dashed
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison of OBC records from the overthrust SEG/EAGE model: 3-D acoustic modelling (a), 3-D elastic modelling (b), conversion of data to the 2-D data after implementation of lateral filtering in the manner of Wapenaar (1992) (c) and 2-D elastic modelling (d).</figcaption></figure>

line shows an ideal case, where the computing time decreases linearly with the number of cores and therefore efficiency equals to 100 per cent. The algorithm shows excellent speedup values for a small number of processor cores (less than 16) and then gradually decreases for a large number of cores (up to 2048). Fig. 2 shows that our algorithm has a typical behaviour of parallel efficiency, calculated using the domain decomposition approach on the fixed-sized model. In fact, the total volume of communication between subdomains becomes more and more considerable with the increasing number of processors. The upgoing communication fills the available limited computer memory bandwidth and thus decreases the total efficiency. However, the simultaneous use of 2048 cores still seems reasonable.

To demonstrate our implementation of the forward modelling algorithm described above, we simulated wave propagation on the reduced 3-D SEG/EAGE overthrust model with an additional water layer on the top to reproduce marine environment. The P-wave velocity model (Fig. 3) was used to create both S-wave velocity and density models. The S-wave velocity was derived using a multiplication factor, which gradually varies with depth from 1.5 to 1.8 and a minimum value equal to 1300 m s⁻¹. The values for density were derived using the Gardner relationship (Gardner *et al.* 1974). The model covers an area of 12.8 km × 3.2 km × 12.8 km (X-, Z- and Y-directions, respectively) and it is discretized with 20 m cubic cells, representing a uniform mesh of 640 × 160 × 640

gridpoints with a total number of 65.536 × 10⁶ nodes. A convolutional perfectly matched layer was used at all boundaries with 10-gridpoints length. A Ricker wavelet with a 5 Hz dominant frequency was used as a source, which was injected in the water layer at 10 m depth with horizontal coordinates at 0.8 km, and 6.4 km in X- and Y- directions. Ocean bottom cables (OBC) were placed on the seafloor at 200 m water depth. An OBC is 10.76 km long, contains 270 receivers at 40 m interval. Fig. 4 shows a wavefield propagated in the model after 1.6 and 4.0 s, respectively. Normally, wave propagation for a single shot would require almost an hour of computation on a single processor computer, but the computational time was reduced to only a few seconds using our efficient parallel implementation on 1024 processor cores. We have also performed wave propagation using an acoustic approximation assuming that S-wave velocity in the medium is zero. Fig. 5 shows OBC records for vertical components for different types of media: (i) 3-D acoustic data, (ii) 3-D elastic data, (iii) 3-D elastic data converted to 2-D data after application of a lateral filter (Wapenaar 1992) and (iv) a 2-D elastic simulation. We clearly observe a significant difference between acoustic and elastic data, as the shear wave energy is fully absent in the acoustic case. Moreover, as mentioned earlier, the amplitudes for most of the events (non-zero-angle of incidence) are also different. It is interesting to note that elastic data sets calculated from 3-D and 2-D models look very similar. However, it will be shown in the next section that the restriction to a
<figure><img src="image_8.png" /><figcaption>Figure 6. 3-D model of the true P-wave velocity used in FWI example 1, taken from the 3-D SEG/EAGE overthrust model shown in Fig. 3.</figcaption></figure>

2-D geometry could significantly reduce the quality of the inversion results.

## 2.2 Inverse problem

The FWI scheme implemented in this study is described more elaborately by Shipp & Singh (2002) and Sears et al. (2010). The inversion estimates detailed P- and S-wave velocity models by minimizing the difference between synthetically calculated seismic data ($d_{syn}$) and observed data ($d_{obs}$) in a least-square sense ($L_2$-norm):

$$s = \sum_{\text{shots}} \int_{0}^{T} dt \sum_{\text{recs}} [d_{\text{syn}} - d_{\text{obs}}]^2, \qquad (5)$$

where $s$ is a misfit function. Eq. (5) gives a quantitative measure of the misfit, which is obtained by summing the differences for all receivers over the total recording time $T$, for each shot gather. We do not include the regularisation term, which is commonly used in traveltime tomography (Hobro et al. 2003). This is because the resolution of the waveform inversion of seismic reflection data is a quarter of the dominant wavelength, and this will be smeared out due to regularisation. Besides the least-square criterion used here, other norms, such as the $L_1$-norm criterion (Djikpéssé & Tarantola 1999) or the logarithmic type norm using the Cauchy criterion (Crase et al.

1990) could be used in the case of low signal-to-noise ratio but here we use the $L_2$ norm for noise-free synthetic data examples.

The inversion is formulated as a local optimization problem using a conjugate gradient algorithm in the time domain. As compared to the frequency domain approach, the time domain approach provides straightforward and efficient implementation of the inversion algorithm using parallel computing facilities, and affordable memory requirements for 3-D elastic FWI. As shown by Tarantola (1984), the cross-correlation of the forward modelled wavefield, $\bar{u}$ with the back-propagated residuals (difference between observed and calculated synthetic data), $\bar{\psi}$ gives the gradient ($g$) of the misfit function

$$g = \sum_{\text{shots}} \int_{0}^{T} dt (\bar{u} \cdot \bar{\psi}). \qquad (6)$$

Eq. (6) is also known as adjoint technique (e.g. Plessix 2006 for review) and provides an estimation of the update for the Lamé parameters ($\delta\lambda, \delta\mu$) and density ($\delta\rho$). Since the choice of parameters has a strong influence on inversion, although theoretically equivalent, the gradients for Lamé parameters are converted to the gradient for P-wave and S-wave velocities ($\delta V_p, \delta V_s$) in the manner of Mora (1987) thus improving inversion convergence (Debski & Tarantola 1995):

$$\delta V_p = 2\rho V_p \delta\lambda, \qquad (7)$$

$$\delta V_s = -4\rho V_s \delta\lambda + 2\rho V_s \delta\mu, \qquad (8)$$

Density is a parameter that is poorly resolved by the inversion process (Tarantola 1986). Therefore, it is not inverted for in this study and is updated using empirical relationship with P-wave velocity (Gardner et al. 1974; Hamilton 1978) in the manner of Shipp & Singh (2002). At the final stage of iteration $n$, the new model ($m_{n+1}$) is calculated by updating the current model ($m_n$)

$$m_{n+1} = m_n - \alpha_n \gamma_n, \qquad (9)$$

where $\gamma_n$ is a pre-conditioned gradient direction and $\alpha_n$ is an optimal step-length. A positive scalar $\alpha$ parameter is calculated using a linear assumption during each iteration (Pica et al. 1990):

$$\alpha = \varepsilon \frac{(d_{\text{pert}} - d_{\text{syn}})^T (d_{\text{syn}} - d_{\text{obs}})}{(d_{\text{pert}} - d_{\text{syn}})^T (d_{\text{pert}} - d_{\text{syn}})}, \qquad (10)$$

<figure><img src="image_9.png" /><figcaption>Figure 7. Vertical slices, extending through a middle of the cross-line axis of the model shown in Fig. 6. Velocity models for true P wave (a), true S wave (b), initial P wave (c) and initial S wave (d). White bold lines indicate the C-PML boundary used for modelling.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Schematic representation of geometry for different approaches in FWI example 1: 3-D propagation combined with 3-D FWI (a), 3-D propagation combined with 2-D FWI (b) and 2-D propagation (local 3-D) combined with 2-D FWI (c). The red stars indicate the source locations and the black dots indicate the position of OBC-receivers.</figcaption></figure>

where $d_{\text{pert}}$ is data generated from the current model $m_n$, perturbed by a small known amount $\epsilon$ (e.g. 1 per cent) in the direction of previously calculated gradient. This constrains the inversion to remain in the vicinity of the current model and makes the linear assumption acceptable. To improve the convergence rate, a line-search could be performed, which would require the solution of at least one more

additional forward calculation and increase the computation time significantly. Therefore, eq. (10) was used for step-length calculation in our study.

At each iteration, the simultaneous inversion of both P- and S-wave velocities requires four forward modelling computations for each shot namely: (1) calculation of synthetics for misfit estimation, (2) wavefield reconstruction, (3) back-propagation of residuals and (4) step-length calculation. Step (2) helps to avoid the storage of the whole 3-D snapshots $\bar{u}$, used for the gradient calculation in eq. (6), which is performed in the reverse order at each time step. For example, a 3-D model with model dimensions of 400 × 400 × 200 and 4000 time steps would require ~9.2 TB of computer memory space for storage of all wavefield components. For the current level of computer power, this volume of information could not be easily stored in random access memory (RAM). Therefore, in this case the overall computing time will be significantly increased due to the required storage on an external disk and consequently increase the reading/writing time of such a considerable volume. Since the time sampling is very small, the wavefield is oversampled. However, the storage of the entire wavefield at every time step could be avoided using a reconstruction technique (Symes 2007). In our implementation, we use the time reversibility of the FD scheme where the wavefield is re-established from the final condition and the information at the boundaries is stored for every time step during the first calculation of synthetic seismograms. By performing one additional FD-simulation, our wavefield reconstruction allows to use fast RAM and therefore reduce the memory requirement by a factor of ~10² in our 3-D inversion examples. This approach enables simple and effective calculation of the cross-correlation between back-propagated residuals and forward modelled wavefields.

In this study, we consider a multiparameter inversion with the goal of inverting P- and S-wave velocity models. In addition to the simultaneous inversion approach used here, where both $Vp$ and $Vs$ are inverted together at the same iteration, one can use either alternate or hierarchical approaches (Freudenreich & Singh 2000). In the alternate approach, the $Vp$ is inverted in the first iteration and then the $Vs$ is inverted in the second iteration, and the process is repeated until the misfit function is considered to be sufficiently reduced. In the hierarchical approach, each elastic parameter is inverted in order of their sensitivity during several sequential inversion iterations while keeping other parameters fixed. For example, the P-wave velocity is inverted first followed by the S-wave velocity and so on. A hybrid approach is a combination of alternate and hierarchical approaches. In our 1-D, 2-D and 3-D inversion tests, the simultaneous and alternate approaches have shown a better convergence. However, we prefer to use the simultaneous approach, which requires less forward modelling steps and hence is more efficient. To adjust the sensitivity of different parameters during the simultaneous inversion (multiparameter inversion) of P- and S-wave velocities, a gradient preconditioning is required. An ideal pre-conditioner would be the inverse of the Hessian matrix, which would provide an exact solution in one single iteration (Mora 1987):

$$H = [D^* C_d^{-1} D + C_m^{-1}]^{-1}, \quad (11)$$

where D is a Fréchet derivative matrix, $C_d$ and $C_m$ are covariance matrices for data and model, respectively, the asterisk indicates a transpose of matrix. However, the estimation of eq. (11) leads to prohibitively expensive calculations. The gradient could be normalized using some reference values, such as their a priori variance (Mora 1987; Pratt et al. 1996) or average background velocity (Mora 1988). We applied an approach adopted by Freudenreich (2002), where a form of covariance matrix containing cross-coupling terms
<table><caption>**Table 1.** Comparison of the CPU time expressed as a total computational time multiplied by the number of computer cores for four different tests in the FWI example 1.</caption><thead><tr><th></th><th>FWI mode</th><th>Total computing time (hr)</th><th>Number of computer cores</th><th>CPU time (hr)</th></tr></thead><tbody><tr><td rowspan="3">Elastic</td><td>(1) 3-D propagation and 3-D inversion</td><td>41.6</td><td>128</td><td>5324.8</td></tr><tr><td>(2) 3-D propagation and 2-D inversion</td><td>5.9</td><td>128</td><td>755.2</td></tr><tr><td>(3) 2-D propagation and 2-D inversion</td><td>6.1</td><td>16</td><td>97.6</td></tr><tr><td>Acoustic</td><td>(4) 3-D propagation and 3-D inversion</td><td>Expected: 10.5</td><td>128</td><td>1344.0</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Comparison of final results for vertical (inline) slices through middle of the model from FWI example 1. Velocities of *P* and *S* waves obtained from test 1 (a and b), test 2 (c and d), test 3 (e and f) and test 4 (g). White bold lines indicate the C-PML boundary.</figcaption></figure>

between *P*- and *S*-wave velocities was used. In this approach the covariance matrix is not explicitly derived, but the effect of the model covariance on the model parameters is calculated. To improve the convergence, an additional gradient preconditioning was applied to boost lower gradient values (Shipp & Singh 2002), in the form of linear scaling with depth.

Although we use the FWI in the time domain, the FWI could be equally applied in the frequency domain (Pratt & Worthington 1990;

Pratt & Shipp 1999). Even though the two approaches are theoretically equivalent for the same frequency range, there is a difference between the two approaches in practice, as each of the strategies has its own advantages and drawbacks (Virieux & Operto 2009). The inversion of discrete frequencies provides a natural workflow for multiscale imaging (Bunks et al. 1995), by successively inverting higher frequencies (Ravaut et al. 2004; Operto et al. 2006; Prieux et al. 2010). However, the implementation of the frequency
<table><caption>**Table 2.** Comparison of relative differences in percentage between the true model and the resulting (elastic/acoustic) FWI models for *P*- and *S*-wave velocities, calculated within a slice extending through the middle of cross-line axis.</caption><thead><tr><th>Models</th><th>Difference for Vp with true model (per cent)</th><th>Difference for Vs with true model (per cent)</th></tr></thead><tbody><tr><td>Initial model</td><td>100</td><td>100</td></tr><tr><td rowspan="3">Elastic FWI</td><td>(1) 3-D propagation and 3-D inversion</td><td>29.2</td><td>28.3</td></tr><tr><td>(2) 3-D propagation and 2-D inversion</td><td>43.2</td><td>32.7</td></tr><tr><td>(3) 2-D propagation and 2-D inversion</td><td>64.3</td><td>129.2</td></tr><tr><td>Acoustic FWI</td><td>(4) 3-D propagation and 3-D inversion</td><td>65.8</td><td>—</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Vertical component of the observed data recorded at one line of receivers passing through the middle of the cross-line axis in example 1.</figcaption></figure>

domain approach for FWI in 3-D has been difficult due to very large memory requirements associated with inversion of the Helmholtz operator (Yingst et al. 2011). Moreover, the inversion of a set of sparse frequencies tends to suffer from ringing artefacts associated with the Gibbs phenomena (Brenders et al. 2012), and this becomes a crucial problem in inverting reflection data and therefore, the applications of the frequency domain FWI have been limited to refraction arrivals using an acoustic approximation. For this reason, the time domain approach is more suitable for inverting reflection data on current computer architectures (Vigh et al. 2009). Moreover, the time domain implementation provides a natural way for time-windowing of data that allows the selection of specific arrivals and consequently decreases the nonlinearity of the inverse problem. In order to address the problem of the memory requirement, some authors have performed 3-D acoustic modelling in the time domain and then performed the inversion in the frequency domain (Sirgue et al. 2008), but such an approach has been limited to refraction arrivals only. In the time domain approach, the number of shots and the number of time steps become crucial factors, but the implementation of the FWI is straightforward as the method can be easily parallelized and run effectively on current computer clusters.

In the following synthetic examples, we assumed that the source function is known. Furthermore, we invert all frequencies simultaneously, and thus we do not implement a multiscale approach (Bunks et al. 1995). However, we do invert both reflection and refraction arrivals. Due to limitation of computer resources, relatively

small models were used in the following FWI examples. Since modelling is scalable, a larger model can easily be tested on large computers.

# 3 SYNTHETIC EXAMPLES

## 3.1 Example 1: FWI of the 3-D SEG/EAGE overthrust model with high contrast layers

In the first example, a subvolume was taken from the SEG/EAGE overthrust model with a complex 3-D structure (Fig. 3). The dimensions of the reduced volume are 4.8 km × 4.8 km × 2.0 km (240 × 240 × 100 gridpoints including boundary conditions) and the total number of gridpoints is 5.76 × 10⁶. We focus our attention on two 150–250-m-thick high velocity layers and a fault, that starts from the right-hand side of the model at about 0.7 km depth and plunges along the inline direction (Figs 6, 7a and b). Using our 3-D elastic propagator, we have generated 2.16 s of data recorded at the OBC-geophones placed just below the seabed at 200 m depth at regular spacing of 40 m. In total 10 201 (101 × 101) receivers recorded the three velocity wavefield components ($v_x$, $v_y$, $v_z$). A Ricker wavelet with a dominant frequency of 5.5 Hz was used as the source wavelet. The grid spacing for wavefield computation was 20 m and the time step was 1.8 ms, small enough to avoid dispersion. In total forty-nine sources (7 × 7) spaced at 600 m were used to generate the observed data. C-PML of 10 gridpoints or 200 m in thickness were used at all boundaries. The starting model was created from the true models by applying a smoothing in both horizontal (inline and cross-line) (X, Y) and vertical (depth) directions (Z). Fig. 7 shows vertical inline velocity slices taken at the middle of the cross-line axis of the true and starting velocities for both P and S waves. Subsequently, using the 3-D elastic synthetic data, we have performed four different FWI experiments in order to assess different approximations used within FWI. These are as follows:

(1) 3-D elastic FWI for *P*- and *S*-wave velocities using 3-D elastic wave propagator;

(2) 2-D elastic FWI for *P*- and *S*-wave velocities using 3-D elastic wave propagator;

(3) 2-D elastic FWI for *P*- and *S*-wave velocities using 2-D elastic wave propagator;

(4) 3-D acoustic FWI for *P*-wave velocity using 3-D acoustic wave propagator.

In the first test, the model dimensions and the geometry of the recording system are identical to that used in generating the observed data (Fig. 8a). For the second test, the 2-D inversion was
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Vertical component of data from test (1) conducted in FWI example 1: initial synthetics (a), final synthetics (b), initial residuals (c) and final residuals (d). All plots are shown at the same scale.</figcaption></figure>

performed along one source-receiver line in the middle of the cross-line axis (Fig. 8b), which has seven shots and 101 receivers. Although the total number of shots was reduced from 49 to seven, the wavefield could propagate in the entire volume using the 3-D wavefield propagator and the velocity was updated in the whole 3-D volume, but the results are shown only along the 2-D profile. In the third test, the 2-D propagation signifies that the effective size of the cross-line dimension (without counting 0.2-km-thick C-PML layers and 0.1-km-thick layers of tapering from each side) was reduced from 4.4 km to only 0.2 km (Fig. 8c). In such a way, the FWI does not take into account the lateral variations in the model while preserving correct amplitudes due to geometrical spreading. This could be considered as 2.5-D FWI, where the medium is close to 2-D geometry (very thin offline direction) while a 3-D point source is used instead of line source. Finally, in the fourth test, the acoustic approximation for the wavefield modelling is applied while keeping the acquisition geometry and the model dimensions identical to those used in the first test, that is full 3-D.

To perform numerical simulations, eight nodes with 16 cores each were used on a local IPGP cluster, providing a total of 128 MPI processors. Each MPI processor had 4 Gb of RAM. In comparison with the first test, which lasted almost two days, all the remaining tests were less computationally expensive (Table 1). In tests (2) and (3), the total computation time decreased proportionally to the number of shots, that is seven instead of 49. Moreover, in the third test, the number of parameters for inversion decreased considerably due to the reduced cross-line dimension (0.8 km instead of 4.8 km). For the acoustic FWI, we have used zero shear wave velocity, without

changing any specific algorithm or modifying the number of parameters, i.e. the grid spacing was the same as for the other tests. Therefore, the total computing time for tests (1) and (4) was of the same order. However, as mentioned earlier, in practice the application of the acoustic approximation could reduce the computing time by a few orders of magnitude, and the difference in the final calculation time will depend on the ratio between the minimum values for *P*- and *S*-wave velocities. Here, in this example, we have utilized a model with a hard water bottom layer. Therefore, the ratio between the smallest value of primary velocity in the water layer and the smallest value of the shear wave velocity in the layer below the water column is equal to about 1.15 ($Vp_{min}$ (in water)/$Vs_{min}$ (at seafloor) = 1500/1300 ≈ 1.15), which is negligible. However, after adopting the grid spacing and the time step, the expected total computing time for acoustic FWI could be reduced almost by a factor of 4.

A comparison of the final results from different inversion tests is shown in Fig. 9. The inversion results from the first test show good agreement between the true and the final FWI models. The inversion has recovered the main features, including the two main high-velocity layers and the fault in both for *P*-wave (Fig. 9a) and *S*-wave (Fig. 9b) velocity images. Even though the shot spacing was 600 m, the 3-D acquisition geometry has allowed to recover the main features. It is interesting to note that the *P*-wave velocity is poorly resolved as compared to the *S*-wave velocity, which is due to shorter wavelengths for *S* waves. In the inverted *S*-wave image, the layer boundaries are sharp and the fault is well delineated. The results from the second test are also quite satisfactory, even though the
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Vertical component of data from test (2) conducted in FWI example 1: initial synthetics (a), final synthetics (b), initial residuals (c) and final residuals (d). All plots are shown at the same scale.</figcaption></figure>

layers and the fault are less clear in *P*- and *S*-wave velocity images (Figs 9c and d), which is due to the reduced number of sources and receivers used for establishing the gradient. The results from the third test are poorer than the previous two tests. The two main high-velocity layers are still recognizable in the final FWI models in both *P*- and *S*-wave velocity images (Figs 9e and f), however the results are strongly contaminated by erroneous updates in the shallow part of the model and the faults are no longer detectable. These artefacts are introduced into the model by FWI because the 3-D effect has not been properly accounted for both in modelling and inversion, which results in inappropriate model updates. The degradation of the results is expected to be even worse for a pure 2-D experiment, where the cross-line dimension is completely absent and an inaccurate 2-D geometrical spreading effect would introduce both amplitude and phase artefacts. Finally, the quality of the *P*-wave velocity image from the fourth test (Fig. 9g) is also poorer than those from tests (1) and (2). Fig. 9(g) shows that the two main high velocity layers and the fault are hardly identifiable. This is because the FWI has tried to fit elastically modelled observed data with acoustically modelled data. It is well known from the Zoeppritz equation that the resulting amplitudes of the seismic response for elastic and acoustic data are different. Furthermore, the mode-converted waves present in the real data were not included in the acoustic modelling. Using 1-D and 2-D synthetic examples, it has been shown (Barnes & Charara 2009; Solano et al. 2013; Borisov et al. 2014) that the errors in acoustic waveform inversion of elastic data will be more in the case of strong contrasts in the elastic parameters, especially for the long offset where the AVO effect is generally much more important. In

such cases, waveform inversion could introduce strong erroneous artefacts into the velocity updates and even diverge quite far from the correct model. From the comparison of the inversion results shown in Fig. 9, we can conclude that tests (2), (3) and (4) while being computationally attractive produce less accurate velocity updates than test (1), that is full 3-D elastic FWI.

A quantitative comparison of the accuracy of FWI results is presented in Table 2, where the differences in percent between the true model and the final FWI models are shown at each gridpoint in the vertical slice, extending through the middle of the cross-line axis. Zero percent will correspond to a perfect match between the inverted and true models, starting from a 100 per cent difference between the true model and the initial model. From this table one could conclude that tests (1) and (2) provide results of superior quality than the results from tests (3) and (4). In the last two examples, the difference in the *P*-wave velocity is more than 60 per cent. The situation is even worse for the *S*-wave velocity in the third test, where the difference is >100 per cent. As we can see, there are no updates for the *S*-wave velocity in the case of acoustic inversion. It should be noted that the difference presented here is valid for the examples shown in this paper, and should not be generalized.

The vertical component of the observed data for one line of receivers placed in the middle of the cross-line axis is shown on Fig. 10. The data are very complex and contain several reflections, direct and refracted arrivals, for both *P* and *S* waves. Note that one of the reflected *P* waves (blue dashed line) has negative moveout due to a dip in the left-side model structure (Figs 7a and b).
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Vertical component of data from test (3) conducted in FWI example 1: initial synthetics (a), final synthetics (b), initial residuals (c) and final residuals (d). All plots are shown at the same scale.</figcaption></figure>

Comparisons of data for FWI tests (1)–(4) are shown in Figs 11–14, respectively, where a vertical component for the initial synthetics (a), final synthetics (b), initial residual (c) and final residuals (d) are shown at the same scale. From visual comparison between the observed data and the initial synthetics we can conclude that they look similar and therefore the initial smooth model is an accurate choice as a starting model, though the amplitude of the initial residuals (Fig. 11c) is quite strong. After 30 iterations, the synthetic data match well with the observed data for tests (1) and (2) (Figs 11b and 12b), which is confirmed by the final residuals (Figs 11d and 12d) where the amplitudes have been reduced significantly. On the other hand, in test (3) the amplitude of final residuals is only slightly reduced (Fig. 13d), while in test (4) it has remained almost at the same level (Fig. 14d) as in the initial residual. We suggest that the low misfit reduction in the last inversion test could be due to the acoustic approximation used that is not capable of reproducing observed seismic records correctly. It is possible that the results of acoustic inversion could improve by applying some data pre-processing steps, such as trace-normalization to reduce the AVO effect (Dessa *et al.* 2004; Ravaut *et al.* 2004), or by removing shear wave energy using a FK filter.

The above example suggests that 3-D elastic FWI is capable of recovering both P- and S-wave velocities from a sparsely sampled 3-D seismic data. The inversion of 2-D line data using the 3-D elastic inversion can also recover these parameters, but with slightly less accuracy. On the other hand, 2-D FWI of a 2-D profile from 3-D data does not lead to satisfactory results.

## 3.2 Example 2: FWI of the 3-D SEG/EAGE overthrust model with shallow channel

In the second example, the inversion was tested on another part of the SEG/EAGE overthrust synthetic model (Fig. 3). This time we focused our attention on a channel, which is located at 1 km depth (Fig. 15a). Since the first approach consisting of full 3-D elastic propagation and 3-D waveform inversion provides the best results, we tested only this approach in this example. The model dimensions are 6.4 km × 6.4 km in horizontal directions and 1.6 km in depth (320 × 320 × 80 gridpoints including boundary conditions). The grid spacing of 20 m makes a total number of gridpoints of about $8.2 \times 10^6$. A Ricker source wavelet with a dominant frequency of 7 Hz was used and the sampling interval was 3 ms. Again, minimum S-wave velocity in this example equal to 1300 m s⁻¹. In total, 111 sources (11 × 11) sparsely spaced at 500 m were used to generate the synthetic data. 19 821 (141 × 141) receivers regularly spaced at 40 m intervals were placed close to the surface. Each receiver measured three velocity field components. C-PML were used at all boundaries with 10-gridpoints length. The starting model was a smoothed version of the true model in all three directions (X, Y and Z).

The inversion results show good agreement between the true and final FWI models by recovering the main features for both P- and S-wave velocities over the whole depth-scale (Figs 16 and 17). The channel is accurately recovered as well as the background velocity in the upper layers. As previously shown, the S-wave velocity is better resolved than the P-wave velocity. Some footprints from sources are visible on the FWI results for both P- and S-wave
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Vertical component of data from test (4) conducted in FWI example 1: initial synthetics (a), final synthetics (b), initial residuals (c) and final residuals (d). All plots are shown at the same scale.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** 3-D model of a true (a) and starting (b) P-wave velocity used in FWI example 2, taken from the 3-D SEG/EAGE overthrust model shown in Fig. 3.</figcaption></figure>

velocities. These artefacts could be reduced by applying a smooth-
ing filter to the gradient. Fig. 18 shows observed data (a), initial
(c) and final synthetics (d). From visual comparison between the
observed data and the initial synthetics, we can see that they look
very similar and, therefore the initial smooth velocity model is a
good choice as a starting model, in spite of large initial amplitude
residual (Fig. 18e). After 30 iterations of the inversion, the total
misfit was reduced by 95 per cent (Fig. 18b) and the amplitude of
final residuals has decreased significantly (Fig. 18f). To perform
this numerical experiment, 16 quadri-core bi-processor nodes were
used providing 128 MPI processors. Each MPI processor used 1 Gb
of RAM, therefore a total of 128 Gb of RAM was allocated. The
FWI results of 30 iterations were obtained after 3 d of continuous
calculations performed on a local IPGP cluster.

**3.3 Example 3: Checkerboard test**

The third example of 3-D elastic waveform inversion consists of
a checkerboard test, which was performed in order to assess the
resolution of the 3-D elastic FWI scheme as a function of depth for
a given geometry of data acquisition. The model dimensions are
slightly smaller in horizontal directions than the previous model ex-
amples. It is 4.8 km × 4.8 km × 1.6 km (240 × 240 × 80 gridpoints
including boundary conditions). The background *P*- and *S*-wave
velocity models consist of a water layer overlying a constant ve-
locity layer. Checkerboard models for *P* and *S* waves were created
by perturbing the background velocity model with two different
grid sizes: The large size cubes are 800 m (40 gridpoints) in each
<figure><img src="image_18.png" /></figure>

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

<figure><img src="image_20.png" /><figcaption>**Figure 16.** 3-D elastic FWI results for *P*-wave velocity from example 2. (a) Vertical slices through the middle of the true (top panel), starting (middle panel) and inverted (bottom panel) model. (b, c) Horizontal slices at 0.7 km and at 1 km depths, respectively, of the true (left-hand panel), starting (middle panel) and inverted model (right-hand panel).</figcaption></figure>

direction with a velocity perturbation of ±3 and ±4 per cent for *P*- and *S*-wave velocities, respectively, whereas the small size cubes are 160 m (eight gridpoints) in each dimension with a velocity perturbation of ±2 and ±3 per cent for *P* and *S* waves, respectively, super-

imposed over the large size grid perturbation. Eighty-one sources (9 × 9) spaced at 500 m were recorded by 10 201 (101 × 101) equally spaced receivers. A Ricker wavelet with a dominant frequency of 7.5 Hz was the source wavelet. The minimum *S*-wave
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 17.** 3-D elastic FWI results for *S*-wave velocity from example 2. For more explanation cf. Fig. 16.</figcaption></figure>

velocity in this example is 1250 m s⁻¹. The background velocity was taken as a starting model. Again the FWI results show good recovery for both *P*-wave (Fig. 19) and *S*-wave (Fig. 20) velocities. The horizontal slices show that FWI has correctly recovered the amplitude and the form of the perturbations almost in the whole

section, although the resolution has decreased at the bottom of the model and close to the model edges due to poor ray coverage. The half wavelength of the *P* wave is ~160 m and close to the size of the small cubes, whereas the half wavelength of the *S* wave is ~110 m. Therefore, as in the previous example, it could be
<figure><img src="image_23.png" /><figcaption>**Figure 18.** Vertical component of inversion data ($v_z$) from one source, placed in the middle of the model (panels a, c, d, e and f) and misfit (panel b) from FWI example 2. Shot gather from observed (a), initial (c) and final (d) FWI models. (e) Initial residual, and (f) final residuals.</figcaption></figure>

observed that the inverted S-wave velocity model is better resolved than the P-wave velocity model. Models with different relative amplitude perturbations and grid sizes were tested and the 3-D FWI was able to recover models of half of the wavelength particularly in the upper part of the model due to a good ray coverage.

# 4 DISCUSSIONS AND CONCLUSIONS

We have developed a 3-D elastic FWI scheme in the time domain and successfully applied it to three synthetic marine examples; a model

with hard velocity contrasts, an example with a channel structure and checkerboard structures of different sizes. Although our test was performed on a relatively small 3-D model size with sparse source spacing, similar results are expected for a larger scale model, provided adequate computer resources are available. In the first experiment, we have compared several different cases and showed that the full 3-D elastic FWI performs much better than using acoustic and 2-D approximations. It was shown that 2-D restriction introduces significant artefacts, particularly when sparse sources are used. We have also shown that acoustic approximation could impact significantly on final inversion results due to amplitude variation with the
<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 19.** 3-D elastic FWI results for *P*-wave velocity from the checkerboard test. (a) Vertical slices through the middle of the true (top panel), starting (middle panel) and inverted (bottom panel) model. (b, c) Horizontal slices at 0.7 km and at 1.3 km depths of the true (left-hand panel), inverted model (middle panel) with the zoom area (right-hand panel).</figcaption></figure>

offset (AVO) effect, especially in the presence of strong contrast in *S*-wave velocity and density. The second example shows that 3-D elastic waveform inversion could provide excellent results, and can retrieve small-scale features such as shallow channel structures. In

the third example, a checkerboard test was carried out. The results confirm that the resolution of half-wavelength could be achieved in the upper part of the models and generally the inversion of *S*-wave velocity provides a better resolution due to the corresponding
<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 20.** 3-D elastic FWI results for S-wave velocity from the checkerboard test. For more explanation cf. Fig. 19.</figcaption></figure>

shorter wavelength. The retrieved elastic parameters could be used to characterize reservoirs on a fine scale thus providing valuable information about the fluid content and lithology.

The successfulness of our inversion results was achieved by a combination of several factors, a perfectly known source signature, the absence of any noise, and an accurate starting model

ensuring the convergence. We recognize that in practice these limiting factors could reduce the quality of the final FWI models, but our objective here was to show that the 3-D elastic waveform inversion could be performed on present-day computers and show excellent recovery even in the presence of sparse sources.
The results could be further improved using a larger number of sources and finer grid spacing. In the case of small *S*-wave velocity contrast, a 3-D acoustic FWI might be more efficient and may lead to good inversion results. In the case of a 2-D earth model, 2-D elastic inversion may provide good results. However, in the presence of 3-D elastic complex structures, 3-D elastic FWI would be required to fully explain the data and retrieve the model accurately. Here, we have assumed that the density is related to *P*-wave velocity. In reality, one might consider inverting density as well, but the sensitivity of density on data is rather poor, and hence inverting all three parameters, *P*- and *S*-wave velocities and density, may lead to a non-unique solution. In that case, it might be wiser to use a hierarchical inversion approach.

These synthetic tests show that our FWI algorithm is efficient and ready to be implemented on real data sets. However, the application of full 3-D elastic waveform inversion is still limited to relatively small subsets of surveys. At this stage several recently developed techniques might improve the efficiency of the FWI. Borisov & Singh (2013) have showed that the use of local inversion strategy based on a grid injection method (Robertsson & Chapman 2000) could dramatically reduce total computational time for 3-D elastic waveform inversion of time-lapse seismic data. The use of simultaneous sources could lead to a significant reduction in computation time (Neelamani et al. 2010). The anisotropy and attenuation were not considered in our inversion strategy. They could play important role in the convergence of algorithm, but this would increase the computation time further.

## ACKNOWLEDGEMENTS

We thank LITHOS consortium for founding this project. We are also very grateful to James Martin and Tim Sears (BP) for fruitful discussions. This work was granted access to the HPC resources of TGCC under the allocation t2013047092 made by GENCI (Grand Equipment National de Calcul Intensif). A significant part of the calculations in this work was performed on the local cluster at the Institut de Physique du Globe de Paris (IPGP). We would like to thank Johan Robertsson and an anonymous reviewer for constructive comments.

## REFERENCES

* Arnulf, A.F., Harding, A.J., Singh, S.C., Kent, G.M. & Crawford, W., 2012. Fine-scale velocity structure of upper oceanic crust from full waveform inversion of downward continued seismic reflection data at the Lucky Strike Volcano, Mid-Atlantic Ridge, Geophys. Res. Lett., 39(8), doi:10.1029/2012GL051064.

* Arnulf, A.F., Harding, A.J., Singh, S.C., Kent, G.M. & Crawford, W.C., 2014. Nature of upper crust beneath the Lucky Strike volcano using elastic full waveform inversion of streamer data, Geophys. J. Int., 196, 1471–1491.

* Barnes, C. & Charara, M., 2009. The domain of applicability of acoustic full-waveform inversion for marine seismic data, Geophysics, 74(6), WCC91–WCC103.

* Ben-Hadj-Ali, H., Operto, S. & Virieux, J., 2008. Velocity model building by 3D frequency-domain, full-waveform inversion of wide-aperture seismic data, Geophysics, 73(5), VE101–VE117.

* Borisov, D. & Singh, S.C., 2013. An efficient 3D elastic full waveform inversion of time-lapse seismic data using grid injection method, in Proceeding of the 84th SEG Annual International Meeting, Expanded Abstracts, pp. 954–958.

* Borisov, D., Stopin, A. & Plessix, R.E., 2014. Acoustic pseudo-density full waveform inversion in the presence of hard thin beds, in Proceeding of the 76th EAGE Conference & Exhibition, Extended Abstracts.

* Brenders, A.J., Albertin, U. & Mika, J., 2012. Comparison of 3D time-and frequency-domain waveform inversion: benefits and insights of a broadband, discrete-frequency strategy, in Proceeding of the 83rd SEG Annual International Meeting, Expanded Abstracts, pp. 1–5.

* Brossier, R., Operto, S. & Virieux, J., 2009. Seismic imaging of complex onshore structures by 2D elastic frequency-domain full-waveform inversion, Geophysics, 74(6), WCC105–WCC118.

* Bunks, C., Saleck, F.M., Zaleski, S. & Chavent, G., 1995. Multiscale seismic waveform inversion, Geophysics, 60(5), 1457–1473.

* Collino, F. & Tsogka, C., 2001. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media, Geophysics, 66(1), 294–307.

* Crase, E., Pica, A., Noble, M., McDonald, J. & Tarantola, A., 1990. Robust elastic nonlinear waveform inversion: application to real data, Geophysics, 55, 527–538.

* Debski, W. & Tarantola, A., 1995. Information on elastic parameters obtained from the amplitudes of reflected waves, Geophysics, 60, 1426–1436.

* Dessa, J.X., Operto, S., Kodaira, S., Nakanishi, A., Pascal, G., Virieux, J. & Kaneda, Y., 2004. Multiscale seismic imaging of the eastern Nankai trough by full waveform inversion, Geophys. Res. Lett., 31(18), doi:10.1029/2004GL020453.

* Djikpéssé, H.A. & Tarantola, A., 1999. Multiparameter $l_1$ norm waveform fitting: interpretation of Gulf of Mexico reflection seismograms, Geophysics, 64(4), 1023–1035.

* Freudenreich, Y. & Singh, C.S., 2000. Full waveform inversion for seismic data—frequency versus time domain, in Proceedings of the 62nd EAGE Conference & Exhibition, Extended Abstracts.

* Freudenreich, Y.P., 2002. *P*- and *S*-wave velocity estimation from full wave-field inversion of wide-aperture seismic data, PhD thesis, University of Cambridge.

* Gardner, G.H.F., Gardner, L.W. & Gregory, A.R., 1974. Formation velocity and density—the diagnostic basics for stratigraphic traps, Geophysics, 39, 770–780.

* Graves, R.W., 1996. Simulating seismic wave propagation in 3D elastic media, Geophysics, 86, 1091–1106.

* Hamilton, E.L., 1978. Sound velocity-density relations in sea-floor sediments and rocks, J. acoust. Soc. Am., 65, 909–922.

* Hobro, J.W., Singh, S.C. & Minshull, T.A., 2003. Three-dimensional tomographic inversion of combined reflection and refraction seismic traveltime data, Geophys. J. Int., 152(1), 79–93.

* Jaiswal, P., Zelt, C.A., Bally, A.W. & Dasgupta, R., 2008. 2-D traveltime and waveform inversion for improved seismic imaging: Naga thrust and fold belt, India, Geophys. J. Int., 173, 642–658.

* Komatitsch, D. & Martin, R., 2007. An un-split convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation, Geophysics, 72(5), 155–167.

* Lailly, P., 1983. The seismic inverse problem as a sequence of before stack migrations, in Proceedings of the Conference on Inverse Scattering, Theory and Application, Society of Industrial and Applied Mathematics, Expanded Abstracts, pp. 206–220.

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

* Lu, R., Lazaratos, S., Wang, K., Cha, Y.H., Chikichev, I. & Prosser, R., 2013. High-resolution elastic FWI for reservoir characterization, in Proceedings of the 75th EAGE Conference & Exhibition, incorporating SPE EUROPEC 2013, doi:10.3997/2214-4609.20130113.

* Minkoff, S., 2002. Spatial parallelism of a 3D finite difference velocity-stress elastic wave propagation code, SIAM. J. Scient. Comput., 24(1), 1–19.

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

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

* Neelamani, R., Krohn, C.E., Krebs, J.R., Romberg, J.K., Deffenbaugh, M. & Anderson, J.E., 2010. Efficient seismic forward modeling using simultaneous random sources and sparsity, Geophysics, 75(6), WB15–WB27.
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, *J. geophys.
Res.*, **111**, doi:10.1029/2005JB003835.

Pica, A., Diet, J.P. & Tarantola, A., 1990. Nonlinear inversion of seismic reflection data in a laterally invariant medium, *Geophysics*, **55**(3), 284–292.

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

Plessix, R.É., 2009. Three-dimensional frequency-domain full-waveform inversion with an iterative solver, *Geophysics*, **74**(6), WCC149–WCC157.

Plessix, R.É., Baeten, G., de Maag, J.W., ten Kroode, F. & Rujie, Z., 2012. Full waveform inversion and distance separated simultaneous sweeping: a study with a land seismic data set, *Geophys. Prospect.*, **60**(4), 733–747.

Pratt, R.G. & Shipp, R.M., 1999. Seismic waveform inversion in the frequency domain, Part 2: fault delineation in sediments using crosshole data, *Geophysics*, **64**, 902–914.

Pratt, R.G. & Worthington, M.H., 1990. Inverse theory applied to multi-source cross-hole tomography. Part 1: acoustic wave-equation method, *Geophys. Prospect.*, **38**(3), 287–310.

Pratt, R.G., Song, Z. & Warner, M., 1996. Two-dimensional velocity models from wide-angle seismic data by wavefield inversion, *Geophys. J. Int.*, **124**, 323–340.

Prieux, V., Brossier, R., Kommendal, J.R., Barkved, O.I., Operto, S. & Virieux, J., 2010. Application of 2D acoustic frequency-domain full-waveform inversion to OBC wide-aperture data from the Valhall field, in *Proceedings of the 80th SEG Annual International Meeting. Expanded Abstracts*, Vol. 29, pp. 920–924.

Queißer, M. & Singh, S.C., 2013. Localizing CO₂—seismic images versus P-wave velocities from waveform inversion, *Geophysics*, **78**, 1–16.

Randall, C.J., 1989. Absorbing boundary condition for the elastic wave equation: velocity-stress formulation, *Geophysics*, **54**, 1141–1152.

Ravaut, C., Operto, S., Improta, L., Virieux, J., Herrero, A. & Dell'Aversana, P., 2004. Multiscale imaging of complex structures from multifold wide-aperture seismic data by frequency-domain full-waveform tomography: application to a thrust belt, *Geophys. J. Int.*, **159**(3), 1032–1056.

Roberts, M.A., Singh, S. & Hornby, B.E., 2008. Investigation into the use of 2D elastic waveform inversion from look-ahead walk-away VSP surveys, *Geophys. Prospect.*, **56**(6), 883–895.

Robertsson, J.O.A., 1996. A numerical free-surface condition for elastic/
viscoelastic finite-difference modeling in the presence of topography,
*Geophysics*, **61**, 1921–1934.

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

Sears, T.J., Singh, S.C. & Barton, P.J., 2008. Elastic full waveform inversion of multi-component OBC seismic data, *Geophys. Prospect.*, **56**(6), 843–862.

Sears, T.J., Barton, P.J. & Singh, S.C., 2010. Elastic full waveform inversion of multicomponent ocean-bottom cable seismic data: Application to Alba Field, U.K. North Sea, *Geophysics*, **75**(6), R109–R119.

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

Sirgue, L., Etgen, J.T. & Albertin, U., 2008. 3D frequency domain waveform inversion using time domain finite difference methods, in *Proceedings of the 70th EAGE Conference & Exhibition, Extended Abstracts*, F022.

Sirgue, L., Barkved, O.I., Dellinger, J., Etgen, J., Albertin, U. & Kommedal, J.H., 2010. Full waveform inversion: the next leap forward in imaging at Valhall, *First Break*, **28**(4), 65–70.

Solano, C.P., Stopin, A. & Plessix, R.E., 2013. Synthetic study of elastic effects on acoustic full waveform inversion, in *Proceedings of the 75th EAGE Conference & Exhibition Incorporating, SPE EUROPEC 2013, Extended Abstracts*, doi:10.3997/2214-4609.20130294.

Symes, W., 2007. Reverse time migration with optimal checkpointing, *Geo-physics*, **72**(5), SM213–SM221.

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

Tarantola, A., 1986. A strategy for nonlinear elastic inversion of seismic reflection data, *Geophysics*, **51**, 1893–1903.

Toxopeus, G., Fokkema, J.T., Groenenboom, J. & Lin, H.X., 2002. 3-D parallel elastic finite difference wave modeling, in *Proceedings of the 8th Annual Conference of the Advanced School for Computing and Imaging*, eds Deprettere, E.F. et al., Lochem, the Netherlands, pp. 231–237.

Vigh, D.V., Starr, W.E.S. & Dingwall, K.K., 2009. 3D prestack time domain full waveform inversion, in *Proceedings of the 71st EAGE Conference & Exhibition*, Vienna, Extended Abstract, U043.

Vigh, D.V., Starr, W.E.S. & Elapavuluri, P., 2010. Acoustic waveform in-
version applicability on elastic data, in *Proceedings of the 72nd EAGE
Conference & Exhibition*, Barcelona, Extended Abstract, A019.

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

Wapenaar, C.P.A., Verschuur, E. & Herrmann, P., 1992. Amplitude prepro-
cessing of single- and multi-component seismic data, *Geophysics*, **57**,
1178–1188.

Warner, M., Morgan, J., Umpleby, A., Stekl, I. & Guasch, L., 2012. Which physics for full-wavefield seismic inversion?, in *Proceedings of the 74th EAGE Conference & Exhibition, Copenhagen, Denmark, 4–7 June 2012*, Extended Abstract.

Yao, Z. & Margrave, G.F., 2000. Elastic wavefield modelling in 3D by the fourth-order staggered finite-difference technique, *CREWES Res. Rep.*, **12**, 157–167.

Yingst, D.V., Wang, C., Park, J., Bloor, R., Leveille, J. & Farmer, P., 2011. Application of time domain and single frequency waveform inversion to real data, in *Proceedings of the 73rd EAGE Conference & Exhibition*, Vienna, Austria, 23–26 May 2011, Extended Abstracts.