# On the adjoint state method for the gradient computation in full waveform inversion: a complete mathematical derivation for the (visco-)elastodynamics approximation

Ludovic Métivier, Romain Brossier

► To cite this version:

Ludovic Métivier, Romain Brossier. On the adjoint state method for the gradient computation in full waveform inversion: a complete mathematical derivation for the (visco-)elastodynamics approximation. Geophysical Journal International, 2025, 240 (2), pp.942-966. 10.1093/gji/ggae421. hal-05336533

HAL Id: hal-05336533

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

Submitted on 29 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# On the adjoint state method for the gradient computation in full waveform inversion: a complete mathematical derivation for the (visco-)elastodynamics approximation

L. Métivier<fnref n="1" />,<fnref n="2" /> and R. Brossier<fnref n="1" />

Accepted 2024 November 21. Received 2024 November 18; in original form 2024 July 4

## SUMMARY

High-resolution seismic imaging at all scales using full waveform inversion is now routinely used in the industry and in the academy. One key element for the success of this approach is a numerical method, named adjoint state method, originally designed for optimization problems constrained by partial differential equations, a category to which full waveform inversion belongs. This method provides an efficient way to compute the gradient of the full waveform inversion misfit function, which is the most computationally demanding task in the implementation of full waveform inversion. While well known, the complete and rigorous mathematical derivation of the adjoint state method for full waveform inversion remains missing in the scientific bibliography. The aim of this study is to remedy this lack. The derivation is performed in general settings, that is in the elastodynamics approximation, with and without considering viscosity. Through the calculus, the mechanism of the adjoint state strategy makes clear the connection between the incident and adjoint fields, especially regarding their initial and boundary conditions. The impact of introducing the viscosity is carefully analysed. The resulting gradient formulas are analysed and shown to be consistent with already published ones. The generic approach which is adopted also makes it possible to derive misfit function gradients with respect to other quantities than the subsurface mechanical parameters, for instance with respect to the initial or the boundary conditions, which could be of interest for specific applications where the reconstructed parameters are not only volumetric mechanical parameters but boundary parameters or initial field values.

**Key words:** Inverse theory; Computational seismology; Waveform inversion.

## 1 INTRODUCTION

Full waveform inversion (FWI) is now a widely adopted high-resolution seismic imaging tool, in the industry, for imaging the crust and exploration purposes (Plessix & Perkins 2010; Sirgue et al. 2010; Warner et al. 2013; Stopin et al. 2014; Vigh et al. 2014; Operto et al. 2015; Raknes et al. 2015; Solano & Plessix 2019), and in the academy, for global Earth (Bozdağ et al. 2016; Karaoğlu & Romanowicz 2018; Lei et al. 2020; Thrastarson et al. 2022), regional (Yuan et al. 2014; Fichtner & Villaseñor 2015; Lu et al. 2020; Nouibat et al. 2023) and deep crust imaging (Górszczyk et al. 2017) of typically subduction, seismogenic or volcanic zones. Recently, FWI even tends to be applied at smaller scales: from near surface scales (Koehn et al. 2018; Irnaka et al. 2022) to medical imaging (Guasch et al. 2020; Marty et al. 2021). A review on FWI and its application can be found for instance in Virieux et al. (2017).

Unlike tomography methods based on the interpretation of traveltimes, which provide low-resolution estimation of the seismic wave velocities, FWI aims at interpreting the whole recorded signal, that is phases and amplitudes. The benefit is the ability to reconstruct quantitative estimations of the subsurface mechanical parameters, potentially beyond seismic wave velocities (density, attenuation and anisotropy) and with much higher resolution, up to half the minimum propagated wavelength following diffraction tomography analysis (Devaney 1984; Wu & Toksöz 1987; Sirgue & Pratt 2004).

The FWI problem is formulated as the minimization of the mismatch between the seismic data acquired on the field, and what is commonly called synthetic seismic data, which are computed from the solution of partial differential equations modelling the propaga-

<footnote n="1">¹ISTerre, Univ. Grenoble Alpes, F-38000 Grenoble, France. E-mail: ludovic.metivier@univ-grenoble-alpes.fr</footnote>

<footnote n="2">²LJK, CNRS, Univ. Grenoble Alpes, F-38000 Grenoble, France</footnote>
tion of waves in the subsurface. This minimization, through local optimization strategy, requires the repeated solution of these partial differential equations. In 3-D, with frequency content leading to the propagation of several tens to several hundreds of wavelength in all directions (typical problem size for FWI applications), solving such a problem is already a computationally intensive task, requiring the use of parallel computing. Calculating the solution of the FWI problem is thus extremely intense in terms of computational effort. This has notably been one of the obstacle to the generalization of the use of FWI when the method was conceived in the 80s (Lailly 1983; Tarantola 1984): at that time, the available computational resources were not sufficient to apply FWI on realistic scale 3-D problem. This paradigm has changed since the end of years 2000's and the first applications of FWI to 3-D field data in the acoustic approximation at exploration scales (Sirgue et al. 2010) and in elastic approximation at regional scales (Fichtner et al. 2009; Tape et al. 2009). However, FWI remains a computationally heavy task, especially when considering FWI application in the (visco-)elastic approximation with large number of propagated wavelengths (Zhu et al. 2012; Komatitsch et al. 2016; Trinh et al. 2019; Tromp; Cao et al. 2022).

Interestingly, a necessary condition to make FWI applicable, from a computational point of view, is a mathematical method arising from the theory of optimal control, developed by the French mathematician Jacques-Louis Lions. This method is named the adjoint-state method (Lions 1968). It is designed for the problem of estimating, given data acquired on a physical system described by partial differential equations, a parameter of these partial differential equations. This is a category of problems FWI belongs to. The adjoint state strategy provides a method to compute the gradient of the misfit function for such kind of problems. In a nutshell, the gradient is expressed as the product of a wavefield solution of the partial differential equations constraining the optimization problem, and an adjoint wavefield, solution of the adjoint system to these partial differential equations. In the context of FWI, it shows that the gradient of the misfit function corresponds to the sum over each source position of the zero-lag time-correlation between the wavefield emitted by the source, and the wavefield resulting from the time-reversed propagation of the residuals, that is the difference between the observed and synthetic data, from the receivers. This result has been put forward initially by Chavent (1974). Its analogy and its link with the concept of exploding receiver and imaging through migration (Claerbout 1985) has been further explored by Lailly (1983). A review on the application of the adjoint state strategy to seismic imaging has been proposed by Plessix (2006).

This result is key to the implementation of FWI because it reduces the computational cost to evaluate the misfit function gradient to the computation of two wavefields per source position. A direct computation of the gradient would, in comparison, require the evaluation of a Jacobian matrix which, in the best case, demands the solution of a number of wavefields equal to the number of receivers for each source position. The adjoint state method is thus the cornerstone of any implementation of FWI algorithm. Nonetheless, it appears to us that the complete and rigorous mathematical derivation of the application of this method in the general context of elastodynamics equations has not yet been proposed. What is usually presented is a formal derivation, where the adjoint equations appear somehow artificially, especially regarding their initial/final conditions and their boundary conditions (Tromp et al. 2005; Plessix 2006; Fichtner & van Driel 2014; Yang et al. 2016). Basically, it is acknowledged that the adjoint of the wave equation operator is also a wave equation operator “backpropagating in time” (Plessix 2006). The system is complemented by hand with ad-hoc final and boundary conditions.

To remedy this lack, we propose in this study a systematic and rigorous mathematical derivation of the adjoint state equations, for elastodynamics equations, first in the purely elastic case, then considering the addition of viscosity. We show that the initial and boundary conditions need to be incorporated in the definition of the Lagrangian operator so that in the derivation of the adjoint state equations, the corresponding adjoint final and boundary conditions also appear. The related derivation is a little bit tedious and lengthy, however we feel that, given the importance of the adjoint state technique for FWI, a widespread cutting-edge technology for seismic imaging at all scales, it is necessary that a complete and rigorous derivation is presented and exist in the literature. Beyond the derivation of the adjoint state equation, we also provide the associated gradient formulas for the density, the stiffness tensor coefficients, and the quality factor in the visco-elastic case. We show that the expressions we obtain are compatible with the ones existing in the literature. The generic equations that are considered can be specialized to any configuration from isotropic elasticity to the triclinic anisotropy. The generic aspect of the development also makes it easier to see the impact of changes in the partial differential equations, their initial conditions, and their boundary conditions in the adjoint state system. Finally, a byproduct of this general derivation is that the gradient of the misfit function with respect to parameters related to the initial condition or the boundary condition are naturally obtained.

The structure of the study is as follows: We recall the general full waveform inversion formalism and the formal derivation of the gradient through the adjoint state strategy in Section 2. We then present the complete computation of the adjoint state equations in the elastodynamics approximation in Section 3. The generalization to visco-elasticity with standard linear solid mechanisms is presented in Section 4. A discussion is proposed in Section 5. We conclude in Section 6.

## 2 FULL WAVEFORM INVERSION: NOTATIONS AND FORMALISM

We first introduce the observed seismic data as a collection of real-valued functions $d_{\text{obs,s}}(\mathbf{x}_r, t)$, depending on the receiver position $\mathbf{x}_r$ and time $t$, indexed by the source number $s$. We introduce $N_s$ and $N_r$ the total number of source and receiver positions, respectively. The total recording time is denoted by $T$.
FWI consists in solving the following optimization problem

$$ \min_{m \in \mathcal{M}} J(m), \quad J(m) = J(d_{\text{cal}}[m], d_{\text{obs}}), \qquad (1) $$

where $J$ is a positive function

$$ J: d_1, d_2 \rightarrow J(d_1, d_2) \\ \mathcal{D} \times \mathcal{D} \rightarrow \mathbb{R}^+, \qquad (2) $$

with $\mathcal{D}$ the data space. In eq. (1), $d_{\text{cal}}[m]$ corresponds to the calculated data, parametrized by $m(\mathbf{x})$, a subsurface model living in the model space $\mathcal{M} \subset L^2(\Omega)$. In the remainder of the study, $\Omega \subset \mathbb{R}^d$ is the domain of interest, in which we reconstruct the subsurface mechanical parameters $m(\mathbf{x})$, with $d$ the physical dimension. Usually $d = 3$, sometimes $d = 2$ in specific scenario (synthetic tests, acquisition on a 2-D line).

The calculated data is obtained through the solution of partial differential equations describing the wave propagation in the domain of interest, a two steps process which for now we denote with general notations as

$$ A(m)u_s = b_s, \quad d_{\text{cal},s}(\mathbf{x}_r, t) = R u_s. \qquad (3) $$

In eq. (3), $A(m)$ is a wave propagation operator, involving partial differential operators in space and time, and depending on the subsurface model $m(\mathbf{x})$. The source term $b_s(\mathbf{x}, t)$ corresponds to the source $s$, $u_s(\mathbf{x}, t)$ is the incident field propagating from $b_s(\mathbf{x}, t)$, and $R$ is an extraction operator which can be written as a convolution with a Dirac delta function $\delta(\mathbf{x} - \mathbf{x}_r)$ as

$$ d_{\text{cal},s}(\mathbf{x}_r, t) = (R u_s)(\mathbf{x}_r, t) = \int_{\Omega} \delta(\mathbf{x} - \mathbf{x}_r) u_s(\mathbf{x}, t) d\mathbf{x}, \quad r = 1, \dots, N_r. \qquad (4) $$

The minimization problem (2) is solved through local optimization. Starting from an initial guess $m_0$, a sequence of subsurface models $m_k$ is built using the following recurrence

$$ m_{k+1} = m_k + \alpha_k \Delta m_k, \qquad (5) $$

where $\alpha_k$ is a linesearch parameter satisfying the Wolfe's conditions and $\Delta m_k$ is a descent direction (Nocedal & Wright 2006). The state-of-the art technique consists in relying on the $l$-BFGS quasi-Newton method, which computes $\Delta m_k$ as

$$ \Delta m_k = -Q_k \nabla J(m_k), \qquad (6) $$

where $Q_k$ is an approximation of the inverse Hessian operator $(\nabla^2 J(m_k))^{-1}$, built from $l$ previously stored gradient $[\nabla J(m_k), \nabla J(m_{k-1}), \dots, \nabla J(m_{k-l+1})]$. The $l$-BFGS approximation formula provides a low-rank representation of the inverse Hessian operator based on successive rank-2 updates of an initial estimation which can be taken as the identity. Besides, it provides a matrix-free formalism where $Q_k$ is never assembled. Only its action on the gradient $\nabla J(m_k)$ is evaluated through a double recursion formula (Nocedal 1980; Nocedal & Wright 2006). The $l$-BFGS method can be combined efficiently with a pre-conditioner, that is a matrix $P$ approximating the inverse Hessian operator such that

$$ P \approx (\nabla^2 J(m_k))^{-1}. \qquad (7) $$

In this case, $P$ serves as the initial estimation of $(\nabla^2 J(m_k))^{-1}$ which is further updated by the $l$-BFGS formula (Métivier & Brossier 2016).

The main topic of this study is how the gradient $\nabla J(m)$ can be computed thanks to the adjoint state method (Plessix 2006). In this section, we recall the formal result. Following the adjoint-state technique, the gradient is obtained by the summation over all the sources of the zero-lag time-correlation between the incident wavefield $u_s(\mathbf{x}, t)$ and the adjoint wavefield $\lambda_s(\mathbf{x}, t)$, such that

$$ \nabla J(m) = \sum_{s=1}^{N_s} \int_{0}^{T} \left\langle \frac{\partial A}{\partial m} u_s(\mathbf{x}, t), \lambda_s(\mathbf{x}, t) \right\rangle dt, \qquad (8) $$

where $\lambda_s$ is the solution of the adjoint wave equation

$$ A(m)^{\dagger} \lambda_s = R^{\dagger} \frac{\partial J(d_{\text{cal},s}, d_{\text{obs},s})}{\partial d_{\text{cal},s}}. \qquad (9) $$

In eq. (9), the subscript $\dagger$ means “adjoint operator”. The source, for the adjoint field, is related to the first-order derivatives of the misfit function with respect to the calculated data, and is localized at the receiver positions through the operator $R^\dagger$. When a least-squares misfit function is used for $J$, this first-order derivatives simply give the difference between observed and calculated data, that is the residuals. Also, the notation $\langle \dots \rangle$ denotes a scalar product in the wavefield space. The importance of this result is that it provides a way to compute the gradient of the FWI misfit function at the cost of solving two wave propagation problems per source. It is key to any efficient implementation of FWI algorithm, and to the application of FWI in real 3-D configurations.

The purpose of this study is to make a complete mathematical derivation of this result, in the general case where the operator $A(m)$ corresponds to the (visco-)elastodynamics equations.
# 3 THE ADJOINT STATE METHOD FOR ELASTIC FWI

## 3.1 Incident field equations in the elastic approximation

We consider a hexagonal domain $\Omega$ such that

$$ \Omega = \Omega_x \times \Omega_y \times \Omega_z \subset \mathbb{R}^3, \text{ with } \Omega_x = [0, L_x], \Omega_y = [0, L_y], \Omega_z = [0, L_z], \quad (10) $$

and a time interval $[0, T] \subset \mathbb{R}^+$. In the remainder of the study, the spatial position is denoted by $\mathbf{x} = (x, y, z) \in \Omega$ and the time by $t \in [0, T]$. The general 3-D elastodynamics equations are derived from the second Newton's law, relating the particle acceleration to the external forces, and the Hooke's law, relating the stress tensor to the strain tensor. This corresponds to the following hyperbolic system of partial differential equations:

$$ \left\{ \begin{array}{l} \rho \partial_t v - D(\nabla)\sigma = F_v, \\ \partial_t \sigma - C D(\nabla)^T v = F_\sigma, \end{array} \right. \quad (11) $$

where $v(\mathbf{x}, t) \in L^2(\Omega \times [0, T])^3$ is the particle velocity,

$$ v(\mathbf{x}, t) = [v_x \ v_y \ v_z]^T, \quad (12) $$

$\sigma(\mathbf{x}, t) \in L^2(\Omega \times [0, T])^6$ is the stress tensor, represented as a 6-component vector using Voigt notation

$$ \sigma(\mathbf{x}, t) = [\sigma_{xx} \ \sigma_{yy} \ \sigma_{zz} \ \sigma_{yz} \ \sigma_{xz} \ \sigma_{xy}]^T := [\sigma_1 \ \sigma_2 \ \sigma_3 \ \sigma_4 \ \sigma_5 \ \sigma_6]^T, \quad (13) $$

and $D(\nabla)$ is the matrix of spatial differential operators

$$ D(\nabla) = \begin{pmatrix} \partial_x & 0 & 0 & 0 & \partial_y \\ 0 & \partial_y & 0 & \partial_z & 0 & \partial_x \\ 0 & 0 & \partial_z & \partial_y & \partial_x & 0 \end{pmatrix}. \quad (14) $$

The parameters describing the subsurface elastic properties are the density $\rho(\mathbf{x}) \in C(\Omega, R)$ and the stiffness matrix $C(\mathbf{x}) \in C(\Omega, S_6^+)$ where $S_6^+$ is the space of symmetric positive 6 × 6 matrix, such that

$$ C(\mathbf{x}) = \begin{pmatrix} C_{11} & C_{12} & C_{13} & C_{14} & C_{15} & C_{16} \\ C_{22} & C_{23} & C_{24} & C_{25} & C_{26} \\ C_{33} & C_{34} & C_{35} & C_{36} \\ C_{44} & C_{45} & C_{46} \\ C_{55} & C_{56} \\ C_{66} \end{pmatrix}. \quad (15) $$

The representation of the stiffness tensor $C$ as a second-order tensor (a matrix) is due to the use of Voigt notation. The source is represented by the term $F_v \in L^2(\Omega \times [0, T])^3$ which gathers the components of a force source, and the term $F_\sigma \in L^2(\Omega \times [0, T])^6$ which gathers the components of a moment tensor source.

In a more compact form, the elastic system can be recast as

$$ S \partial_t w - A(\nabla)w = F, \quad (16) $$

where

$$ w = \begin{bmatrix} v \\ \sigma \end{bmatrix}, \quad S = \begin{pmatrix} \rho I_3 & 0 \\ 0 & C^{-1} \end{pmatrix}, \quad A(\nabla) = \begin{pmatrix} 0 & D(\nabla) \\ D(\nabla)^T & 0 \end{pmatrix}, \quad F = \begin{bmatrix} F_v \\ C^{-1} F_\sigma \end{bmatrix}. \quad (17) $$

In (17), $I_3$ denotes the rank 3 identity matrix and $C^{-1}$ is the compliance matrix<fnref n="1" />. In a more conventional hyperbolic form, this system is equivalent to

$$ S \partial_t w - A_x \partial_x w - A_y \partial_y w - A_z \partial_z w = F, \quad (18) $$

where

$$ A_x = \begin{pmatrix} 0_3 & D_x \\ D_x^T & 0_6 \end{pmatrix}, \quad A_y = \begin{pmatrix} 0_3 & D_y \\ D_y^T & 0_6 \end{pmatrix}, \quad A_z = \begin{pmatrix} 0_3 & D_z \\ D_z^T & 0_6 \end{pmatrix}, \quad (19) $$

with

$$ D_x = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 1 & 0 \end{pmatrix}, \quad D_y = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \end{pmatrix}, \quad D_z = \begin{pmatrix} 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 \end{pmatrix}, \quad (20) $$

<footnote n="1">¹By definition, $C$ is symmetric positive definite therefore it admits an inverse which is also symmetric positive definite</footnote>
and where $0_3, 0_6$ denote the rank 3 and rank 6 zero square matrices, respectively. The matrices $A_x, A_y$, and $A_z$ are symmetric, which is inherited from the symmetry of $A(\nabla)$. As $S$ is symmetric (by construction), this indicates that the elastodynamics system is symmetrizable with $S$ as a symmetrizer. Such systems are known to be energy conservative (Kreiss & Lorenz 2004). This is consistent with the analysis performed for instance in Yang et al. (2016).

For the incident field, it is usual to consider zero initial conditions $w(\mathbf{x}, 0) = 0$: the medium is at rest before it is excited by the source $F$. For the boundary conditions, a medium with infinite extension has to be simulated. Perfectly matched layers (PML) techniques have become the standard for modelling electromagnetic waves and acoustic waves in such (semi)-infinite media (Bérenger 1994). There is no such consensus when it comes to elastodynamics equations, the main reason being the instability of the PML for anisotropic media (Becache et al. 2003; Métivier et al. 2014). For general elastodynamics equations, it is thus usual to consider a combination of sponge layers (Cerjan et al. 1985) and first-order absorbing boundary conditions, as is done for instance in the spectral-element based full waveform modelling and inversion code named SEM46 (Trinh et al. 2019; Cao et al. 2022). This is the strategy we consider in this study.

From a mathematical standpoint, the sponge layer techniques consists in adding a zero-order term in the first-order hyperbolic system (Métivier et al. 2014). The original domain $\Omega$ is surrounded by a layer such that the computation domain becomes $\tilde{\Omega} \supset \Omega$. The wavefield is damped in this layer so that it reaches zero at the external boundary. The hyperbolic system becomes

$$S\partial_t w - A_x \partial_x w - A_y \partial_y w - A_z \partial_z w + Bw = F, \quad (21)$$

where the zero-order term $B$ is such that

$$B = b(\mathbf{x})I_9, \quad b(\mathbf{x}) \in C^\infty(\tilde{\Omega}), \quad \forall \mathbf{x} \in \Omega, \quad b(\mathbf{x}) = 0. \quad (22)$$

The coefficient $b(\mathbf{x})$ is non-zero only in the layer $\tilde{\Omega} \setminus \Omega$, and grows smoothly as $\mathbf{x}$ gets further from the border between the domain of interest $\Omega$ and the layer. Typically exponential or cosine functions are used. For convenience, we also introduce the 3-by-3 upper diagonal block of $B$ as $B_v$ and the 6-by-6 lower diagonal block of $B$ as $B_\sigma$ such that in developed form, the incident field equations write

$$\begin{cases} \rho \partial_t v - D(\nabla)\sigma + B_v v = F_v, \\ C^{-1} \partial_t \sigma - D(\nabla)^T v + B_\sigma \sigma = C^{-1} F_\sigma. \end{cases} \quad (23)$$

For the sake of simplicity, we now keep the full domain (with sponge layers) denoted by $\Omega$ and keep in mind that a part of it is “non-physical”, corresponding to sponge layers. The damping of the wavefield by sponge layers is complemented by transparent boundary conditions at the external boundary of the layer, as formulated in Lysmer & Kuhlemeyer (1969), Komatitsch (1997), and Kouroussis et al. (2011). This is an adaptation of the Clayton–Enquist strategy to the elastodynamics case. The normal stress is related to the normal velocity, while the shear stress is related to the tangential velocity. If we consider the border at $x = 0$ for instance, this transparent boundary condition can be written as

$$\left\{ \begin{aligned} \sigma_{xx}(0, y, z, t) &= \rho V_P(0, y, z) v_x(0, y, z, t), \\ \sigma_{xy}(0, y, z, t) &= \rho V_S(0, y, z) v_y(0, y, z, t), \\ \sigma_{xz}(0, y, z, t) &= \rho V_S(0, y, z) v_z(0, y, z, t), \end{aligned} \right. \quad (24)$$

while at $x = L_x$ we have

$$\left\{ \begin{aligned} \sigma_{xx}(L_x, y, z, t) &= -\rho V_P(L_x, y, z) v_x(L_x, y, z, t), \\ \sigma_{xy}(L_x, y, z, t) &= -\rho V_S(L_x, y, z) v_y(L_x, y, z, t), \\ \sigma_{xz}(L_x, y, z, t) &= -\rho V_S(L_x, y, z) v_z(L_x, y, z, t), \end{aligned} \right. \quad (25)$$

where $V_P(\mathbf{x})$ and $V_S(\mathbf{x})$ are respectively the $P$-wave and $S$-wave velocities, derived from the stiffness tensor $C(\mathbf{x})$. The same type of boundary conditions is applied on all sides of the medium, except on top of the domain, at the surface $z = 0$, where a free surface condition is imposed (zero traction). Using matrix operators, the boundary conditions applied on the six faces of $\Omega$ can be written

$$\left\{ \begin{aligned} (D_x \sigma)(0, y, z, t) - M_x(0, y, z, t)v(0, y, z, t) &= 0, \\ (D_x \sigma)(L_x, y, z, t) + M_x(L_x, y, z, t)v(L_x, y, z, t) &= 0, \\ (D_y \sigma)(x, 0, z, t) - M_y(x, 0, z, t)v(x, 0, z, t) &= 0, \\ (D_y \sigma)(x, L_y, z, t) + M_y(x, L_y, z, t)v(x, L_y, z, t) &= 0, \\ (D_z \sigma)(x, y, 0, t) &= 0, \\ (D_z \sigma)(x, y, L_z, t) + M_z(x, y, L_z, t)v(x, y, L_z, t) &= 0, \end{aligned} \right. \quad (26)$$

where

$$M_x = \begin{pmatrix} \rho V_P & 0 & 0 \\ 0 & \rho V_S & 0 \\ 0 & 0 & \rho V_S \end{pmatrix}, \quad M_y = \begin{pmatrix} \rho V_S & 0 & 0 \\ 0 & \rho V_P & 0 \\ 0 & 0 & \rho V_S \end{pmatrix}, \quad M_z = \begin{pmatrix} \rho V_S & 0 & 0 \\ 0 & \rho V_S & 0 \\ 0 & 0 & \rho V_P \end{pmatrix}. \quad (27)$$
We obtain a more compact expression of these boundary conditions by introducing the operators $B_x^\pm, B_y^\pm, B_z^\pm$ as

$$B_i^\pm = [\pm M_i D_i], \quad i \in \{x, y\}, \quad B_z^- = [0_3 D_z], \quad B_z^+ = [M_z D_z], \qquad (28)$$

which give the following expression for the boundary conditions

$$\left\{
\begin{array}{l}
(B_x^- w)(0, y, z, t) = 0, \\
(B_x^+ w)(L_x, y, z, t) = 0, \\
(B_y^- w)(x, 0, z, t) = 0, \\
(B_y^+ w)(x, L_y, z, t) = 0, \\
(B_z^- w)(x, y, 0, t) = 0, \\
(B_z^+ w)(x, y, L_z, t) = 0.
\end{array}
\right.
\qquad (29)$$

Note that in practice, to implement a free surface boundary condition on the $z = 0$ plane the damping coefficients $b(x)$ are adapted so that no damping is introduced on the top part of the model.

We end up with the following incident field equations, first in condensed form

$$\left\{
\begin{array}{l}
S \partial_t w - A_x \partial_x w - A_y \partial_y w - A_z \partial_z w + B w = F, \\
w(\mathbf{x}, t) = 0, \\
(B_x^- w)(0, y, z, t) = 0, \\
(B_x^+ w)(L_x, y, z, t) = 0, \\
(B_y^- w)(x, 0, z, t) = 0, \\
(B_y^+ w)(x, L_y, z, t) = 0, \\
(B_z^- w)(x, y, 0, t) = 0, \\
(B_z^+ w)(x, y, L_z, t) = 0.
\end{array}
\right.
\qquad (30)$$

and in developed form

$$\left\{
\begin{array}{l}
\rho \partial_t v - D(\nabla)\sigma + B_v v = F_v, \\
C^{-1} \partial_t \sigma - D(\nabla)^T v + B_\sigma \sigma = C^{-1} F_\sigma, \\
v(\mathbf{x}, t) = 0, \quad \sigma(\mathbf{x}, t) = 0, \\
(D_x \sigma)(0, y, z, t) - M_x(0, y, z, t)v(0, y, z, t) = 0, \\
(D_x \sigma)(L_x, y, z, t) + M_x(L_x, y, z, t)v(L_x, y, z, t) = 0, \\
(D_y \sigma)(x, 0, z, t) - M_y(x, 0, z, t)v(x, 0, z, t) = 0, \\
(D_y \sigma)(x, L_y, z, t) + M_y(x, L_y, z, t)v(x, L_y, z, t) = 0, \\
(D_z \sigma)(x, y, 0, t) = 0, \\
(D_z \sigma)(x, y, L_z, t) + M_z(x, y, L_z, t)v(x, y, L_z, t) = 0.
\end{array}
\right.
\qquad (31)$$

The condensed form of the equation will be useful to derive the adjoint state equation in the next section. The developed form will be useful to interpret them and make clear links between incident and adjoint equations.

## 3.2 FWI and adjoint state equations in the elastic approximation

In the elastodynamics approximation, the FWI problem is generally posed as the minimization of the least-squares distance between the observed and calculated particle velocities. The corresponding formulation is thus

$$\min_m J(m), J(m) = \frac{1}{2} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} \left[ (R_i w[m])(\mathbf{x}_r, t) - d_{\text{obs},i}(\mathbf{x}_r, t) \right]^2 dt, \qquad (32)$$

where $m$ gathers the model parameters

$$m = [\rho \ C]^T, \qquad (33)$$

$w[m]$ is solution of the elastodynamics eqs (30) and $R_i$ is a restriction operator such that

$$(R_i w)(\mathbf{x}_r, t) = \int_{\Omega} \delta(\mathbf{x} - \mathbf{x}_r) w_i(\mathbf{x}, t) dx. \qquad (34)$$

The operator $R_i$ extracts the $i$-th component of the wavefield $w$ at the receiver positions $\mathbf{x}_r$.

As such, the FWI problem can be interpreted as a constrained optimization problem, for which the misfit function is based on the least-squares distance between calculated and observed data, and where the calculated data satisfy constraints associated with the solution of the elastodynamics equations. This is the problem we consider in this section. The extension to more general misfit functions is proposed in the discussion. Note also that we consider here a single source position. The extension to multiple sources is straightforward as long as the misfit function is built as the summation of the misfit contribution from each individual source. In this case, the gradient of the misfit
function is simply the sum of the gradient associated with each individual sources. In what follows we consider this is the case and focus on the gradient contribution from a single source using the adjoint state method.

For a general constrained optimization problem with a misfit term $f(x)$ and equality constraints $g(x)$, written as

$$ \min_{x} f(x), \quad \text{s.t.} \quad g(x) = 0, \qquad (35) $$

it is useful operator to consider the Lagrangian operator. This operator is defined as

$$ \mathcal{L}(x, \lambda) = f(x) + \langle g(x), \lambda \rangle, \qquad (36) $$

where $\langle \cdot, \cdot \rangle$ denotes an appropriate scalar product and $\lambda$ a variable adjoint to $x$. This operator is useful because it helps to characterize the solution of such constrained optimization problem. Indeed, a solution can be found as a saddle-point of the Lagrangian operator, that is a couple $x^*, \lambda^*$ which solves the problem

$$ \min_{x} \max_{\lambda} \mathcal{L}(x, \lambda). \qquad (37) $$

The definition of this operator also extends to the case of inequality constraints, and it is a central tool in optimization theory (Bonnans et al. 2006; Nocedal & Wright 2006).

In our case, our interest for the Lagrangian operator is somehow different, as we do not intend to characterize the solution of the FWI problem, but rather to find an alternative (indirect) way to compute the gradient of the misfit function, which turns out to be cheaper in computation cost than a direct method. To this purpose, we first need to introduce the following Euclidean scalar products, depending on if the quantities are space-time dependent, only space dependent, or defined at the boundaries with a time-dependence.

$$ \left\{ \begin{array}{l} \forall (u, v) \in L^2(\Omega \times [0, T])^n \times L^2(\Omega \times [0, T])^n, \\ (u, v)_{\Omega \times [0, T]} = \sum_{k=1}^{n} \int_{0}^{T} \int_{\Omega} u_k(\mathbf{x}, t) v_k(\mathbf{x}, t) \mathrm{d}\mathbf{x}\mathrm{d}t, \\ \\ \forall (u, v) \in L^2(\Omega_y \times \Omega_z \times [0, T])^n \times L^2(\Omega_y \times \Omega_z \times [0, T])^n, \\ (u, v)_{\Omega_y \times \Omega_z \times [0, T]} = \sum_{k=1}^{n} \int_{\Omega_y \times \Omega_z} u_k(y, z, t) v_k(y, z, t) \mathrm{d}y \mathrm{d}z \mathrm{d}t, \\ \\ \forall (u, v) \in L^2(\Omega_x \times \Omega_z \times [0, T])^n \times L^2(\Omega_x \times \Omega_z \times [0, T])^n, \\ (u, v)_{\Omega_x \times \Omega_z \times [0, T]} = \sum_{k=1}^{n} \int_{\Omega_x \times \Omega_z} u_k(x, z, t) v_k(x, z, t) \mathrm{d}x \mathrm{d}z \mathrm{d}t, \\ \\ \forall (u, v) \in L^2(\Omega_x \times \Omega_y \times [0, T])^n \times L^2(\Omega_x \times \Omega_y \times [0, T])^n, \\ (u, v)_{\Omega_x \times \Omega_y \times [0, T]} = \sum_{k=1}^{n} \int_{\Omega_x \times \Omega_y} u_k(x, y, t) v_k(x, y, t) \mathrm{d}x \mathrm{d}y \mathrm{d}t. \end{array} \right. \qquad (38) $$

These scalar products are required to introduce the Lagrangian operator $\mathcal{L}$ associated with our optimization problem. The latter is the sum of: a first term $\mathcal{A}$, associated with the misfit function

$$ \mathcal{A} = \frac{1}{2} \sum_{i=1}^{3} \sum_{r=1}^{N_i} \int_{0}^{T} |(R_i w)(\mathbf{x}_r, t) - d_{\text{obs},i}(\mathbf{x}_r, t)|^2 dt, \qquad (39) $$

a second term $\mathcal{B}$, associated with the constraint that the wavefield $w$ needs to be the solution of the elastodynamics partial differential equations

$$ \mathcal{B} = (S(m)\partial_t w - A_x\partial_x w - A_y\partial_y w - A_z\partial_z w + Bw - F, \lambda)_{\Omega \times [0, T]}, \qquad (40) $$

a third term $\mathcal{C}$, associated with the constraint that the wavefield $w$ needs to satisfy the initial conditions

$$ \mathcal{C} = (w(\mathbf{x}, 0), \mu)_{\Omega}, \qquad (41) $$

a fourth term $\mathcal{D}_x$, associated with the constraints that the wavefield $w$ needs to satisfy the boundary conditions in the $x$ direction

$$ \mathcal{D}_x = ((B_x^- w)(0, y, z, t), v_x^-)_{\Omega_y \times \Omega_z \times [0, T]} + ((B_x^+ w)(L_x, y, z, t), v_x^+)_{\Omega_y \times \Omega_z \times [0, T]}, \qquad (42) $$
a fifth term $\mathcal{D}_y$, associated with the constraints that the wavefield $w$ needs to satisfy the boundary conditions in the $y$ direction

$$ \mathcal{D}_y = ((B_y^- w)(x, 0, z, t), v_y^-)_{\Omega_x \times \Omega_z \times [0, T]} + ((B_y^+ w)(x, L_y, z, t), v_y^+)_{\Omega_x \times \Omega_z \times [0, T]}, \quad (43) $$

and a sixth term $\mathcal{D}_z$, associated with the constraints that the wavefield $w$ needs to satisfy the boundary conditions in the $z$ direction

$$ \mathcal{D}_z = ((B_z^- w)(x, y, 0, t), v_z^-)_{\Omega_x \times \Omega_y \times [0, T]} + ((B_z^+ w)(x, y, L_z, t), v_z^+)_{\Omega_x \times \Omega_y \times [0, T]}. \quad (44) $$

The elastic-FWI Lagrangian operator $\mathcal{L}$ can thus be written

$$ \begin{aligned} \mathcal{L}(m, w, \lambda, \mu, v_x^+, v_x^-, v_y^+, v_y^-, v_z^+, v_z^-) = & \frac{1}{2} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} |(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)|^2 dt + \\ & \underbrace{(S(m)\partial_t w - A_x \partial_x w - A_y \partial_y w - A_z \partial_z w + Bw - F, \lambda)_{\Omega \times [0, T]}}_{\mathcal{B}} + \underbrace{(w(\mathbf{x}, 0), \mu)_{\Omega}}_{\mathcal{C}} + \\ & \underbrace{((B_x^- w)(0, y, z, t), v_x^-)_{\Omega_y \times \Omega_z \times [0, T]} + ((B_x^+ w)(L_x, y, z, t), v_x^+)_{\Omega_y \times \Omega_z \times [0, T]}}_{\mathcal{D}_x} + \\ & \underbrace{((B_y^- w)(x, 0, z, t), v_y^-)_{\Omega_x \times \Omega_z \times [0, T]} + ((B_y^+ w)(x, L_y, z, t), v_y^+)_{\Omega_x \times \Omega_z \times [0, T]}}_{\mathcal{D}_y} + \\ & \underbrace{((B_z^- w)(x, y, 0, t), v_z^-)_{\Omega_x \times \Omega_y \times [0, T]} + ((B_z^+ w)(x, y, L_z, t), v_z^+)_{\Omega_x \times \Omega_y \times [0, T]}}_{\mathcal{D}_z} \end{aligned} \quad (45) $$

The incident field $w$ and adjoint field $\lambda$ live in the time-space domain such that

$$ w \in L^2(\Omega \times [0, T])^9, \quad \lambda \in L^2(\Omega \times [0, T])^9. \quad (46) $$

The adjoint field $\mu$, associated with the initial conditions, lives in the space domain only

$$ \mu \in L^2(\Omega)^9. \quad (47) $$

The adjoint fields $v_x^\pm$, $v_y^\pm$ and $v_z^\pm$, associated with boundary conditions, live on boundaries with a time dependence, such that

$$ v_x^\pm \in L^2(\Omega_y \times \Omega_z \times [0, T])^3, \quad v_y^\pm \in L^2(\Omega_x \times \Omega_z \times [0, T])^3, \quad v_z^\pm \in L^2(\Omega_x \times \Omega_y \times [0, T])^3. \quad (48) $$

We gather the latter in the adjoint field $v = [v_x^- \ v_x^+ \ v_y^- \ v_y^+ \ v_z^- \ v_z^+]$.

Now, we consider $w(m)$ the solution of the incident field eqs (30). We have

$$ \mathcal{L}(m, w(m), \lambda, \mu, v) = J(m), \quad \text{and} \quad \partial_m \mathcal{L}(m, w(m), \lambda, \mu, v) = \nabla J(m)^T, \quad (49) $$

with the convention that partial derivatives $\partial_m$ yield row vectors and partial derivatives $\nabla$ yield column vectors. Developing (49) yields

$$ \int_0^T \langle \partial_m S(m) \partial_t w(m)(\mathbf{x}, t), \lambda(\mathbf{x}, t) \rangle dt + \partial_w \mathcal{L}(m, w(m), \lambda, \mu, v) \partial_m w(m) = \nabla J(m)^T, \quad (50) $$

with $\langle \cdot, \cdot \rangle$ the Euclidean scalar product on $\mathbb{R}^9$. The essence of the adjoint state method is to determine the adjoint fields $\lambda(m)$, $\mu(m)$, $v(m)$ such that

$$ \partial_w \mathcal{L}(m, w(m), \lambda(m), \mu(m), v(m)) = 0. \quad (51) $$

Provided we find such fields, the gradient of the misfit function simply writes

$$ \nabla J(m) = \int_{0}^{T} \langle \partial_{m} S(m) \partial_{t} w(m), \lambda(m) \rangle dt. \quad (52) $$

An equivalent condition to finding $\lambda(m)$, $\mu(m)$ and $v(m)$ as in eq. (51) is the following

$$ \forall q \in L^2(\Omega \times [0, T])^9, \quad (\partial_w \mathcal{L}(m, w(m), \lambda(m), \mu(m), v(m)), q)_{\Omega \times [0, T]} = 0. \quad (53) $$

Basically, we consider the directional derivative in the direction $q$ for any $q$ instead of considering the derivative operator itself. Developing the eq. (53) is the core of the adjoint-state derivation. Consider a first-order Taylor expansion of $\mathcal{L}$ with respect to its second variable $w$. We have,

$$ \forall q \in L^2(\Omega \times [0, T])^9, \quad \mathcal{L}(m, w+q, \lambda, \mu, v) - \mathcal{L}(m, w, \lambda, \mu, v) = (\partial_w \mathcal{L}(m, w, \lambda, \mu, v), q)_{\Omega \times [0, T]} + O(\|q\|^2), \quad (54) $$

The latter formula tells us that the quantity we want to evaluate, namely $(\partial_w \mathcal{L}(m, w, \lambda, \mu, v), q)_{\Omega \times [0, T]}$ is the term which depends linearly on $q$ when considering the difference $\mathcal{L}(m, w+q, \lambda, \mu, v) - \mathcal{L}(m, w, \lambda, \mu, v)$. To compute it, we can thus expand this difference and extract the part which depends linearly on $q$. Interestingly, as we consider a linear forward problem (the elastic system in eq. 30 is linear with respect
to $w$, as well as the boundary conditions and the initial condition), the only nonlinear term in the Lagrangian operator, when considering its dependence to the second variable $w$, comes from the misfit function $J(m)$.

We thus have

$$ \begin{aligned} \mathcal{L}(m, w+q, \lambda, \mu, v) - \mathcal{L}(m, w, \lambda, \mu, v) = & \\ & \underbrace{\frac{1}{2} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} |(R_i(w+q))(x_r, t) - d_{\text{obs},i}(x_r, t)|^2 dt}_{\mathcal{T}} \\ & - \frac{1}{2} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} |(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)|^2 dt + \\ & (S(m)\partial_t q - A_x \partial_x q - A_y \partial_y q - A_z \partial_z q + Bq, \lambda)_{\Omega \times [0,T]} + (q(\mathbf{x}, 0), \mu)_{\Omega} + \\ & ((B_x^- q)(0, y, z, t), v_x^-)_{\Omega_y \times \Omega_z \times [0,T]} + ((B_x^+ q)(L_x, y, z, t), v_x^+)_{\Omega_y \times \Omega_z \times [0,T]} + \\ & ((B_y^- q)(x, 0, z, t), v_y^-)_{\Omega_x \times \Omega_z \times [0,T]} + ((B_y^+ q)(x, L_y, z, t), v_y^+)_{\Omega_x \times \Omega_z \times [0,T]} + \\ & ((B_z^- q)(x, y, 0, t), v_z^-)_{\Omega_x \times \Omega_y \times [0,T]} + ((B_z^+ q)(x, y, L_z, t), v_z^+)_{\Omega_x \times \Omega_y \times [0,T]} \end{aligned} \quad (55) $$

We develop the term noted $\mathcal{T}$ of the right-hand side of eq. (55). Relying on the linearity of the extraction operator $R$, we have

$$ \begin{aligned} \mathcal{T} = & \frac{1}{2} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} |(R_i(w+q))(x_r, t) - d_{\text{obs},i}(x_r, t)|^2 dt = \\ & \frac{1}{2} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} |(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)|^2 dt + \frac{1}{2} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} |(R_i q)(x_r, t)|^2 dt \\ & + \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} [(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)] (R_i q)(x_r, t) dt. \end{aligned} \quad (56) $$

Therefore we have

$$ \begin{aligned} & \frac{1}{2} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} |(R_i(w+q))(x_r, t) - d_{\text{obs},i}(x_r, t)|^2 dt - \frac{1}{2} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} |(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)|^2 dt = \\ & \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} [(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)] (R_i q)(x_r, t) dt + \frac{1}{2} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} |(R_i q)(x_r, t)|^2 dt \end{aligned} \quad (57) $$

The last term $\frac{1}{2} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} |(R_i q)(x_r, t)|^2 dt$ is obviously in $O(\|q\|^2)$, therefore we can write

$$ \begin{aligned} \mathcal{L}(m, w+q, \lambda, \mu, v) - \mathcal{L}(m, w, \lambda, \mu, v) = & \\ & \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} [(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)] (R_i q)(x_r, t) dt + O(\|q\|^2) + \\ & (S(m)\partial_t q - A_x \partial_x q - A_y \partial_y q - A_z \partial_z q + Bq, \lambda)_{\Omega \times [0,T]} + (q(\mathbf{x}, 0), \mu)_{\Omega} + \\ & ((B_x^- q)(0, y, z, t), v_x^-)_{\Omega_y \times \Omega_z \times [0,T]} + ((B_x^+ q)(L_x, y, z, t), v_x^+)_{\Omega_y \times \Omega_z \times [0,T]} + \\ & ((B_y^- q)(x, 0, z, t), v_y^-)_{\Omega_x \times \Omega_z \times [0,T]} + ((B_y^+ q)(x, L_y, z, t), v_y^+)_{\Omega_x \times \Omega_z \times [0,T]} + \\ & ((B_z^- q)(x, y, 0, t), v_z^-)_{\Omega_x \times \Omega_y \times [0,T]} + ((B_z^+ q)(x, y, L_z, t), v_z^+)_{\Omega_x \times \Omega_y \times [0,T]} \end{aligned} \quad (58) $$

Therefore, by identification between the formal expansion in eq. (54) and the actual expansion in eq. (58) we deduce that

$$ \begin{aligned} (\partial_w \mathcal{L}(m, w, \lambda, \mu, v), q)_{\Omega \times [0,T]} = & \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} [(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)] (R_i q)(x_r, t) dt + \\ & (S(m)\partial_t q - A_x \partial_x q - A_y \partial_y q - A_z \partial_z q + Bq, \lambda)_{\Omega \times [0,T]} + (q(\mathbf{x}, 0), \mu)_{\Omega} + \\ & ((B_x^- q)(0, y, z, t), v_x^-)_{\Omega_y \times \Omega_z \times [0,T]} + ((B_x^+ q)(L_x, y, z, t), v_x^+)_{\Omega_y \times \Omega_z \times [0,T]} + \\ & ((B_y^- q)(x, 0, z, t), v_y^-)_{\Omega_x \times \Omega_z \times [0,T]} + ((B_y^+ q)(x, L_y, z, t), v_y^+)_{\Omega_x \times \Omega_z \times [0,T]} + \\ & ((B_z^- q)(x, y, 0, t), v_z^-)_{\Omega_x \times \Omega_y \times [0,T]} + ((B_z^+ q)(x, y, L_z, t), v_z^+)_{\Omega_x \times \Omega_y \times [0,T]} \end{aligned} \quad (59) $$

In what follows, playing with integration by parts and carefully considering boundary terms, we translate the condition (53) into conditions on the adjoint fields $\lambda, \mu$ and $v$. The main rationale behind the mathematical transformations we present hereafter is to isolate the variable $q$, so as to explicitly write the right-hand side of (59) as scalar products between $q$ and quantities yet to determine, depending on $\lambda, \mu$, and $v$. Satisfying the condition (53) will therefore imply that these quantities are equal to 0, which will yield the so-called adjoint-state formula relating $\lambda, \mu$ and $v$. To achieve this, we rely mostly on integration by parts, the symmetry of the operators involved, and their invariance with respect to time or space.
We start with the time-dependent part. Thanks to the symmetry of $S(m)$, the fact that $S(m)$ does not depend on time, and using a time-integration by parts, we have

$$ (S(m)\partial_t q, \lambda)_{\Omega \times [0,T]} + (q(\mathbf{x}, 0), \mu)_{\Omega} = - (q, S(m)\partial_t \lambda)_{\Omega \times [0,T]} + (q(\mathbf{x}, T), S(m)\lambda(\mathbf{x}, T))_{\Omega} + (q(\mathbf{x}, 0), \mu - S(m)\lambda(\mathbf{x}, 0))_{\Omega}. \quad (60) $$

Next, we consider the space-dependent part, starting with the $x$ direction. We use the symmetry of $A_x$, the fact that it does not depend on space, and an integration by parts. We obtain

$$ \begin{aligned} (-A_x \partial_x q, \lambda)_{\Omega \times [0,T]} &= (\partial_x q, -A_x \lambda)_{\Omega \times [0,T]} \\ &= (q, A_x \partial_x \lambda)_{\Omega \times [0,T]} \\ &\quad - (q(L_x, y, z, t), A_x \lambda(L_x, y, z, t))_{\Omega_y \times \Omega_z \times [0,T]} \\ &\quad + (q(0, y, z, t), A_x \lambda(0, y, z, t))_{\Omega_y \times \Omega_z \times [0,T]}. \end{aligned} \quad (61) $$

We thus obtain

$$ \begin{aligned} (-A_x \partial_x q, \lambda)_{\Omega \times [0,T]} &+ ((B_x^- q)(0, y, z, t), v_y^-)_{\Omega_y \times \Omega_z \times [0,T]} + ((B_x^+ q)(L_x, y, z, t), v_y^+)_{\Omega_y \times \Omega_z \times [0,T]} = \\ (q, A_x \partial_x \lambda)_{\Omega \times [0,T]} &+ (q(0, y, z, t), (B_x^-)^T v_y^- + A_x \lambda(0, y, z, t))_{\Omega_y \times \Omega_z \times [0,T]} \\ &+ (q(L_x, y, z, t), (B_x^+)^T v_y^+ - A_x \lambda(L_x, y, z, t))_{\Omega_y \times \Omega_z \times [0,T]} \end{aligned} \quad (62) $$

Similarly, we have for the $y$ direction

$$ \begin{aligned} (-A_y \partial_y q, \lambda)_{\Omega \times [0,T]} &+ ((B_y^- q)(x, 0, z, t), v_y^-)_{\Omega_x \times \Omega_z \times [0,T]} + ((B_y^+ q)(x, L_y, z, t), v_y^+)_{\Omega_x \times \Omega_z \times [0,T]} = \\ (q, A_y \partial_y \lambda)_{\Omega \times [0,T]} &+ (q(x, 0, z, t), (B_y^-)^T v_y^- + A_y \lambda(x, 0, z, t))_{\Omega_x \times \Omega_z \times [0,T]} \\ &+ (q(x, L_y, z, t), (B_y^+)^T v_y^+ - A_y \lambda(x, L_y, z, t))_{\Omega_x \times \Omega_z \times [0,T]} \end{aligned} \quad (63) $$

and for the $z$ direction

$$ \begin{aligned} (-A_z \partial_z q, \lambda)_{\Omega \times [0,T]} &+ ((B_z^- q)(x, y, 0, t), v_z^-)_{\Omega_x \times \Omega_y \times [0,T]} + ((B_z^+ q)(x, y, L_z, t), v_z^+)_{\Omega_x \times \Omega_y \times [0,T]} = \\ (q, A_z \partial_z \lambda)_{\Omega \times [0,T]} &+ (q(x, y, 0, t), (B_z^-)^T v_z^- + A_z \lambda(x, y, 0, t))_{\Omega_x \times \Omega_y \times [0,T]} \\ &+ (q(x, y, L_z, t), (B_z^+)^T v_z^+ - A_z \lambda(x, y, L_z, t))_{\Omega_x \times \Omega_y \times [0,T]}. \end{aligned} \quad (64) $$

The sponge layer term treatment is straightforward as there are no integration by parts nor boundary terms. We have

$$ (Bq, \lambda)_{\Omega \times [0,T]} = (q, B\lambda)_{\Omega \times [0,T]}, \quad (65) $$

from the symmetry of the sponge layer operator $B$. Finally, for the misfit function term, we have

$$ \begin{aligned} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} [(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)] (R_i q)(x_r, t) dt &= \\ \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} [(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)] \int_{\Omega} q_i(\mathbf{x}, t) \delta(\mathbf{x} - \mathbf{x}_r) d\mathbf{x}dt &= \\ \sum_{i=1}^{3} \int_{\Omega} \int_{0}^{T} q_i(\mathbf{x}, t) \sum_{r=1}^{N_r} [(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)] \delta(\mathbf{x} - \mathbf{x}_r) d\mathbf{x}dt &= (q, F_{\lambda})_{\Omega \times [0,T]}, \end{aligned} \quad (66) $$

where $F_{\lambda}(\mathbf{x}, t) = [F_{\lambda,1} \ F_{\lambda,2} \ F_{\lambda,3} \ 0 \ 0 \ 0 \ 0 \ 0 \ 0]^{T}$ with

$$ F_{\lambda,i} = \sum_{r=1}^{N_r} [(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)] \delta(\mathbf{x} - \mathbf{x}_r), \quad i = 1, \dots, 3. \quad (67) $$

Assembling all, we obtain

$$ \begin{aligned} \forall q \in L^2(\Omega \times [0, T])^9, \\ (\partial_t \mathcal{L}(m, w, \lambda, \mu, \nu), q)_{\Omega \times [0,T]} &= (-S(m)\partial_t \lambda + A_x \partial_x \lambda + A_y \partial_y \lambda + A_z \partial_z \lambda + B\lambda + F_\lambda, q)_{\Omega \times [0,T]} + \\ & (q(\mathbf{x}, T), S(m)\lambda(\mathbf{x}, T))_{\Omega} + (q(\mathbf{x}, 0), \mu - S(m)\lambda(\mathbf{x}, 0))_{\Omega} + \\ & (q(0, y, z, t), (B_x^-)^T v_x^- + A_x \lambda(0, y, z, t))_{\Omega_y \times \Omega_z \times [0,T]} + (q(L_x, y, z, t), (B_x^+)^T v_x^+ - A_x \lambda(L_x, y, z, t))_{\Omega_y \times \Omega_z \times [0,T]} + \\ & (q(x, 0, z, t), (B_y^-)^T v_y^- + A_y \lambda(x, 0, z, t))_{\Omega_x \times \Omega_z \times [0,T]} + (q(x, L_y, z, t), (B_y^+)^T v_y^+ - A_y \lambda(x, L_y, z, t))_{\Omega_x \times \Omega_z \times [0,T]} + \\ & (q(x, y, 0, t), (B_z^-)^T v_z^- + A_z \lambda(x, y, 0, t))_{\Omega_x \times \Omega_y \times [0,T]} + (q(x, y, L_z, t), (B_z^+)^T v_z^+ - A_z \lambda(x, y, L_z, t))_{\Omega_x \times \Omega_y \times [0,T]}. \end{aligned} \quad (68) $$
The only possibility to satisfy condition (53) is thus that all the terms in the scalar product in front of q cancel in eq. (68). We thus obtain

$$
\left\{
\begin{array}{l}
S(m)\partial_t\lambda - A_x\partial_x\lambda - A_y\partial_y\lambda - A_z\partial_z\lambda - B\lambda = F_\lambda, \\
S(m)\lambda(\mathbf{x}, T) = 0, \\
A_x\lambda(0, y, z, t) = -(B_x^-)^T(0, y, z)v_x^-, \\
A_x\lambda(L_x, y, z, t) = (B_y^+)^T(L_x, y, z)v_y^+, \\
A_y\lambda(x, 0, z, t) = -(B_y^-)^T(x, 0, z)v_y^-, \\
A_y\lambda(x, L_y, z, t) = (B_y^+)^T(x, L_y, z)v_y^+, \\
A_z\lambda(x, y, 0, t) = -(B_z^-)^T(x, y, 0)v_z^-, \\
A_z\lambda(x, y, L_z, t) = (B_z^+)^T(x, y, L_z)v_z^+,
\end{array}
\right.
\qquad (69)
$$

with the additional condition

$$
\mu = S(m)\lambda(\mathbf{x}, 0). \tag{70}
$$

From (69), we see that the wavefield $\lambda$ is solution of a wavefield equation similar to the one satisfied by the incident wavefield $w$, with a source term $F_\lambda$ which is a force source in the three directions of space. For the direction $x$ (respectively $y$ and $z$), the source force corresponds to the difference between the calculated and observed particle velocity $v_x$ (resp. $v_y$ and $v_z$) localized at the receiver positions through the Dirac function $\delta(x - x_r)$. We still observe the following important differences between this system and the incident field eqs (30).

(i) The sign of the zero-order sponge layer term $B$ is opposite.

(ii) The zero initial condition on $w$ is now a zero final condition on $\lambda$.

(iii) The boundary conditions are not expressed the same.

In what follows, we make explicit the link between the incident and adjoint equations.

We start with re-expressing the boundary conditions. For the x-boundary conditions at $x = 0$, based on the definition of the operators $A_x$ and $B_x^-$ we have

$$
\left\{
\begin{aligned}
D_x \lambda_{\sigma}(0, y, z, t) &= M_x(0, y, z) v_x^-(y, z, t), \\
D_x^T \lambda_{\nu}(0, y, z, t) &= -D_x^T v_x^-(y, z, t),
\end{aligned}
\right.
\qquad (71)
$$

where we decompose the adjoint wavefield into its three first components $\lambda_v$ and its six last components $\lambda_\sigma$ such that $\lambda = [\lambda_v \ \lambda_\sigma]^T$, exactly in the same way the wavefield $w$ is decomposed into particle velocity $v$ and stress $\sigma$. We have $D_x D_x^T = I_3$, therefore the second equation of (71) yields

$$
\lambda_v(0, y, z, t) = -v_x^-(y, z, t), \tag{72}
$$

and the x-boundary condition at $x = 0$ for the adjoint wavefield $\lambda$ is

$$
D_x \lambda_{\sigma}(0, y, z, t) + M_x(0, y, z) \lambda_v(0, y, z, t) = 0. \tag{73}
$$

For the $x$-boundary condition at $x = L_x$, we have

$$
\left\{
\begin{aligned}
D_x \lambda_{\sigma}(L_x, y, z, t) &= M_x(L_x, y, z)v_x^+(y, z, t), \\
D_x^T \lambda_{\nu}(L_x, y, z, t) &= D_x^T v_x^+(y, z, t).
\end{aligned}
\right.
\qquad (74)
$$

therefore, similarly,

$$
\lambda_v(L_x, y, z, t) = v_x^+(y, z, t), \tag{75}
$$

and the $x$-boundary at $x = L_x$ on $\lambda$ is

$$
D_x \lambda_{\sigma}(L_x, y, z, t) - M_x(L_x, y, z) \lambda_v(L_x, y, z, t) = 0. \tag{76}
$$

We observe that there is a sign difference for the boundary conditions of the incident field (30) and the adjoint field. We can apply the same strategy for the y-boundary conditions, and we obtain

$$
\left\{
\begin{aligned}
D_y \lambda_{\sigma}(x, 0, z, t) + M_y(x, 0, z) \lambda_v(x, 0, z, t) &= 0, \\
D_y \lambda_{\sigma}(x, L_y, z, t) - M_y(x, L_y, z) \lambda_v(x, L_y, z, t) &= 0.
\end{aligned}
\right.
\qquad (77)
$$

For the $z$-boundary conditions, the situation is the same, except at $z = 0$, because of the free surface boundary condition applied at this boundary. We have

$$
A_z \lambda(x, y, 0, t) = -(B_z^-)^T(x, y, 0) v_z^-, \qquad (78)
$$
which, based on the definition of $A_z$ and $B_z^-$, is equivalent to

$$ \begin{cases} D_z \lambda_{\sigma}(x, y, 0, t) = 0 \\ D_z^T \lambda_v(x, y, 0, t) = -D_z^T \nu_z^- \end{cases} \qquad (79) $$

Therefore the z-boundary conditions can be expressed as

$$ \begin{cases} D_z \lambda_{\sigma}(x, y, 0, t) = 0 \\ D_z \lambda_{\sigma}(x, y, L_z, t) - M_z(x, y, L_z) \lambda_v(x, y, L_z, t) = 0. \end{cases} \qquad (80) $$

The adjoint system of eqs (69) can thus be rewritten as

$$ \begin{cases} S(m)\partial_t\lambda - A_x\partial_x\lambda - A_y\partial_y\lambda - A_z\partial_z\lambda - B\lambda = F_{\lambda}, \\ S(m)\lambda(\mathbf{x}, T) = 0, \\ D_x\lambda_{\sigma}(0, y, z, t) + M_x(0, y, z)\lambda_v(0, y, z, t) = 0, \\ D_x\lambda_{\sigma}(L_x, y, z, t) - M_x(L_x, y, z)\lambda_v(L_x, y, z, t) = 0, \\ D_y\lambda_{\sigma}(x, 0, z, t) + M_y(x, 0, z)\lambda_v(x, 0, z, t) = 0, \\ D_y\lambda_{\sigma}(x, L_y, z, t) - M_y(x, L_y, z)\lambda_v(x, L_y, z, t) = 0, \\ D_z\lambda_{\sigma}(x, y, 0, t) = 0, \\ D_z\lambda_{\sigma}(x, y, L_z, t) - M_z(x, y, L_z)\lambda_v(x, y, L_z, t) = 0, \end{cases} \qquad (81) $$

with the additional conditions

$$ \mu = S(m)\lambda(\mathbf{x}, 0), \qquad (82) $$

and

$$ \left\{ \begin{aligned} v_x^-(y, z, t) &= -\lambda_v(0, y, z, t), \\ v_x^+(y, z, t) &= \lambda_v(L_x, y, z, t), \\ v_y^-(x, z, t) &= -\lambda_v(x, 0, z, t), \\ v_y^+(x, z, t) &= \lambda_v(x, L_y, z, t), \\ v_z^-(x, y, t) &= -\lambda_v(x, y, 0, t), \\ v_z^+(x, y, t) &= \lambda_v(x, y, L_z, t). \end{aligned} \right. \qquad (83) $$

We now proceed to the change of variable $t' = T - t$, to transform the adjoint system into a prograde system. We introduce

$$ \lambda'(x, y, z, t') = \lambda(x, y, z, t). \qquad (84) $$

We have

$$ \partial_{t'} \lambda'(x, y, z, t') = -\partial_t \lambda(x, y, z, t) \qquad (85) $$

and

$$ \lambda'(x, y, z, 0) = \lambda(x, y, z, T) = 0 \qquad (86) $$

therefore $\lambda'(x, y, z, t')$ is solution of

$$ \begin{cases} S(m)\partial_t'\lambda' + A_x\partial_x\lambda' + A_y\partial_y\lambda' + A_z\partial_z\lambda' + B\lambda' = -F_{\lambda}', \\ S(m)\lambda'(\mathbf{x}, 0) = 0, \\ D_x\lambda_{\sigma}'(0, y, z, t') + M_x(0, y, z)\lambda_v'(0, y, z, t') = 0, \\ D_x\lambda_{\sigma}'(L_x, y, z, t') - M_x(L_x, y, z)\lambda_v'(L_x, y, z, t') = 0, \\ D_y\lambda_{\sigma}'(x, 0, z, t') + M_y(x, 0, z)\lambda_v(x, 0, z, t') = 0, \\ D_y\lambda_{\sigma}'(x, L_y, z, t') - M_y(x, L_y, z)\lambda_v'(x, L_y, z, t') = 0, \\ D_z\lambda_{\sigma}'(x, y, 0, t') = 0, \\ D_z\lambda_{\sigma}'(x, y, L_z, t') - M_z(x, y, L_z)\lambda_v(x, y, L_z, t') = 0, \end{cases} \qquad (87) $$

where $F_{\lambda}'(\mathbf{x}, t') = F_{\lambda}(\mathbf{x}, t)$ and with the additional conditions

$$ \mu = S(m)\lambda'(\mathbf{x}, T), \qquad (88) $$

and

$$ \left\{ \begin{aligned} v_x^-(y, z, t') &= -\lambda_v'(0, y, z, t'), \\ v_x^+(y, z, t') &= \lambda_v'(L_x, y, z, t'), \\ v_y^-(x, z, t') &= -\lambda_v'(x, 0, z, t'), \\ v_y^+(x, z, t') &= \lambda_v'(x, L_y, z, t'), \\ v_z^-(x, y, t') &= -\lambda_v'(x, y, 0, t'), \\ v_z^+(x, y, t') &= \lambda_v'(x, y, L_z, t'). \end{aligned} \right. \qquad (89) $$
Now, the only difference remaining with the incident field eqs (30) is related to the sign in front of the partial derivatives operator $A_x$, $A_y$ and $A_z$, and the sign in the boundary conditions. For the latter, the relationship between the adjoint stress and the adjoint particle velocity at the boundaries are opposed when compared with the incident stress and incident particle velocity

Consider the solution $\lambda^1 = [\lambda_v^1 \ \lambda_\sigma^1]^T$ of the system of eqs (87). Consider the solution $\lambda^2 = [\lambda_v^2 \ \lambda_\sigma^2]^T$ of the incident field eqs (31) with the source term $F_\lambda$. In the following we demonstrate that

$$ \lambda_v^1(\mathbf{x}, t) = \lambda_v^2(\mathbf{x}, t), \quad \lambda_\sigma^1(\mathbf{x}, t) = -\lambda_\sigma^2(\mathbf{x}, t). \tag{90} $$

In other words, solving for the incident field equations with the source term $F_\lambda(x, y, z, T-t)$ yields the adjoint particle velocity field and the opposite of the adjoint stress field.

We have, for $\lambda^1$, dropping the notation $t'$ to $t$ for the sake of simplicity

$$ \left\{ \begin{aligned} & \rho \partial_t \lambda_v^1 + D(\nabla)\lambda_\sigma^1 + B_v \lambda_v^1 = F_{v,\lambda}, \\ & C^{-1} \partial_t \lambda_\sigma^1 + D(\nabla)^T \lambda_v^1 + B_\sigma \lambda_\sigma^1 = 0, \\ & \lambda_v^1(\mathbf{x}, 0) = 0, \quad \lambda_\sigma^1(\mathbf{x}, 0) = 0, \\ & D_x \lambda_\sigma^1(L_x, y, z, t) - M_x(L_x, y, z) \lambda_v^1(L_x, y, z, t) = 0, \\ & D_y \lambda_\sigma^1(x, 0, z, t) + M_y(x, 0, z) \lambda_v^1(x, 0, z, t) = 0, \\ & D_z \lambda_\sigma^1(x, y, z, t) - M_z(x, y, z) \lambda_v^1(x, y, z, t) = 0, \\ & D_z \lambda_\sigma^1(x, y, L_z, t) - M_z(x, y, L_z) \lambda_v^1(x, y, L_z, t) = 0, \end{aligned} \right. \tag{91} $$

On the other hand, for $\lambda^2$, we have

$$ \left\{ \begin{aligned} & \rho \partial_t \lambda_v^2 - D(\nabla)\lambda_\sigma^2 + B_v \lambda_v^2 = F_{v,\lambda}, \\ & C^{-1} \partial_t \lambda_\sigma^2 - D(\nabla)^T \lambda_v^2 + B_\sigma \lambda_\sigma^2 = 0, \\ & \lambda_v^2(\mathbf{x}, 0) = 0, \quad \lambda_\sigma^2(\mathbf{x}, 0) = 0, \\ & D_x \lambda_\sigma^2(0, y, z, t) - M_x(0, y, z) \lambda_v^2(0, y, z, t) = 0, \\ & D_x \lambda_\sigma^2(L_x, y, z, t) + M_x(L_x, y, z) \lambda_v^2(L_x, y, z, t) = 0, \\ & D_y \lambda_\sigma^2(x, 0, z, t) - M_y(x, 0, z) \lambda_v^2(x, 0, z, t) = 0, \\ & D_z \lambda_\sigma^2(x, y, z, t) + M_z(x, y, z) \lambda_v^2(x, y, z, t) = 0, \\ & D_z \lambda_\sigma^2(x, y, L_z, t) + M_z(x, y, L_z) \lambda_v^2(x, y, L_z, t) = 0, \end{aligned} \right. \tag{92} $$

We introduce the field $\tilde{\lambda} = [\tilde{\lambda}_v \ \tilde{\lambda}_\sigma]^T$ such that

$$ \left\{ \begin{aligned} \tilde{\lambda}_v &= \lambda_v^2 - \lambda_v^1, \\ \tilde{\lambda}_\sigma &= \lambda_\sigma^2 + \lambda_\sigma^1. \end{aligned} \right. \tag{93} $$

We demonstrate that $\tilde{\lambda} = 0$. We have, by combination of (91) and (92)

$$ \left\{ \begin{aligned} & \rho \partial_t (\lambda_v^2 - \lambda_v^1) - D(\nabla)(\lambda_\sigma^2 + \lambda_\sigma^1) + B_v (\lambda_v^2 - \lambda_v^1) = 0 \\ & C^{-1} \partial_t (\lambda_\sigma^2 + \lambda_\sigma^1) - D(\nabla)^T (\lambda_v^2 - \lambda_v^1) + B_\sigma (\lambda_\sigma^2 + \lambda_\sigma^1) = 0, \\ & (\lambda_v^2 - \lambda_v^1)(\mathbf{x}, 0) = 0, \quad (\lambda_\sigma^2 + \lambda_\sigma^1)(\mathbf{x}, 0) = 0, \\ & D_x (\lambda_\sigma^2(0, y, z, t) + \lambda_\sigma^1(0, y, z, t)) - M_x(0, y, z) (\lambda_v^2(0, y, z, t) - \lambda_v^1(0, y, z, t)) = 0 \\ & D_x (\lambda_\sigma^2(L_x, y, z, t) + \lambda_\sigma^1(L_x, y, z, t)) + M_x(L_x, y, z) (\lambda_v^2(L_x, y, z, t) - \lambda_v^1(L_x, y, z, t)) = 0 \\ & D_y (\lambda_\sigma^2(x, 0, z, t) + \lambda_\sigma^1(x, 0, z, t)) - M_y(x, 0, z) (\lambda_v^2(x, 0, z, t) - \lambda_v^1(x, 0, z, t)) = 0 \\ & D_y (\lambda_\sigma^2(x, L_y, z, t) + \lambda_\sigma^1(x, L_y, z, t)) + M_y(x, L_y, z) (\lambda_v^2(x, L_y, z, t) - \lambda_v^1(x, L_y, z, t)) = 0 \\ & D_z (\lambda_\sigma^2(x, y, 0, t) + \lambda_\sigma^1(x, y, 0, t)) = 0, \\ & D_z (\lambda_\sigma^2(x, y, L_z, t) + \lambda_\sigma^1(x, y, L_z, t)) + M_z(x, y, L_z) (\lambda_v^2(x, y, L_z, t) - \lambda_v^1(x, y, L_z, t)) = 0 \end{aligned} \right. \tag{94} $$

Therefore, we see that $\tilde{\lambda}$ is solution of

$$ \left\{ \begin{aligned} & \rho \partial_t \tilde{\lambda}_v - D(\nabla) \tilde{\lambda}_\sigma + B_v \tilde{\lambda}_v = 0 \\ & C^{-1} \partial_t \tilde{\lambda}_v - D(\nabla)^T \tilde{\lambda}_v + B_\sigma \tilde{\lambda}_\sigma = 0, \\ & \tilde{\lambda}_v(\mathbf{x}, 0) = 0, \quad \tilde{\lambda}_\sigma(\mathbf{x}, 0) = 0, \\ & D_x \tilde{\lambda}_\sigma(0, y, z, t) - M_x(0, y, z) \tilde{\lambda}_v(0, y, z, t) = 0 \\ & D_x \tilde{\lambda}_\sigma(L_x, y, z, t) + M_x(L_x, y, z) \tilde{\lambda}_v(L_x, y, z, t) = 0 \\ & D_y \tilde{\lambda}_\sigma(x, 0, z, t) - M_y(x, 0, z) \tilde{\lambda}_v(x, 0, z, t) = 0 \\ & D_y \tilde{\lambda}_\sigma(x, L_y, z, t) + M_y(x, L_y, z) \tilde{\lambda}_v(x, L_y, z, t) = 0 \\ & D_z \tilde{\lambda}_\sigma(x, y, 0, t) = 0, \\ & D_z \tilde{\lambda}_\sigma(x, y, L_z, t) + M_z(x, y, L_z) \tilde{\lambda}_v(x, y, L_z, t) = 0 \end{aligned} \right. \tag{95} $$
eq. (95) reveals that $\tilde{\lambda}$ is the solution of the incident field eqs (31) with a zero source term, zero initial condition, and absorbing boundary conditions, that is with no energy injected in the system, neither as source, an initial condition or at the boundaries. For this reason, the solution $\tilde{\lambda}(\mathbf{x}, t)$ is identically null for all $\mathbf{x} \in \Omega$ and $t \in [0, T]$. This concludes the demonstration.

We make few comments about this result. First, it can appear in contradiction with what is usually considered for linear elastic equations. Those are supposed to be self-adjoint, therefore we could think that there exists a transformation of the adjoint field equations showing that the two sets of equations are equivalent. Our findings indicate that this is not the case: there is a change of sign on the adjoint stress variable. However, the incident velocity and adjoint velocity fields can be computed using the same set of equations. The latter observation resolves the contradiction. In practice, second-order elastic equations are usually considered, which involve only the displacement field. Those equations are self-adjoint. We lose this property here because we consider the first-order in time velocity-stress equations for the incident field.

From an implementation point of view, losing this self-adjoint property in this purely elastic case has however no strong implication. Indeed, the same computer code as the one implementing the solution of the incident field equations can be used to compute a solution of the adjoint equations, up to a posterior change of sign to obtain the adjoint stress field. In the following, we provide the gradient formulas for elastic-FWI based on the adjoint-state approach.

### 3.3 Elastic FWI gradient expression

In general, the gradient of the misfit function has to be computed for parameters which are combination of the density $\rho(\mathbf{x})$ and the components of the stiffness tensor $C_{IJ}(\mathbf{x})$, $1 \le I, J \le 6$. Let $\beta(\mathbf{x})$ be such parameter, we denote its dependency on $\rho$ and $C$ as $\beta = \beta(\rho, C)$. The gradient of the misfit function with respect to $\beta$ is obtained by chain rule as

$$ \frac{\partial J}{\partial \beta} = \frac{\partial J}{\partial \rho} \frac{\partial \rho}{\partial \beta} + \sum_{I,J=1}^{6} \frac{\partial J}{\partial C_{IJ}} \frac{\partial C_{IJ}}{\partial \beta}. \qquad (96) $$

The key quantities to compute are thus the gradient of $J$ with respect to $\rho$ and the stiffness tensor coefficients $C_{IJ}$. We have, from eq. (50),

$$ \begin{aligned} \frac{\partial J}{\partial \rho} &= \int_{0}^{T} \left\langle \frac{\partial S}{\partial \rho} \partial_t w(\rho, C), \lambda(\rho, C) \right\rangle dt, \\ \frac{\partial J}{\partial C_{IJ}} &= \int_{0}^{T} \left\langle \frac{\partial S}{\partial C_{IJ}} \partial_t w(\rho, C), \lambda(\rho, C) \right\rangle dt, \end{aligned} \qquad (97) $$

where $w(\rho, C)$ and $\lambda(\rho, C)$ are solutions of the incident and adjoint equations respectively. In eq. (97), we neglect the potential contribution from the source term $C^{-1}F_{\sigma}$. Neglecting this term is not a strong approximation: it can be interpreted as considering that in the case where non-zero moment tensor sources $F_{\sigma}$ are used, they are sufficiently localized in space (which is usual in practice in seismology) and that we compute the gradient with respect to the stiffness tensor coefficients away from these sources. The latter is in practice strongly recommended because updating the subsurface parameters in the vicinity of the source is known to yield instabilities in FWI.

We have

$$ \frac{\partial S}{\partial \rho} = \begin{pmatrix} I_3 & 0 \\ 0 & 0_6 \end{pmatrix}, \qquad (98) $$

therefore from eq. (97) we deduce

$$ \frac{\partial J}{\partial \rho} = \sum_{i \in \{x, y, z\}} \int_{0}^{T} \partial_t v_i(\rho, C)(\lambda_v)_i(\rho, C) dt. \qquad (99) $$

The latter can be rewritten in terms of the displacement field $u$ defined as $\partial_{tt}u = \partial_{t}v$ such that

$$ \frac{\partial J}{\partial \rho} = \sum_{i \in \{x, y, z\}} \int_{0}^{T} \partial_{tt} u_{i}(\rho, C)(\lambda_{v})_{i}(\rho, C) dt. \qquad (100) $$

Regarding the gradient with respect to the components of the stiffness tensor $C$, we first observe that

$$ \frac{\partial (SS^{-1})}{\partial C_{IJ}} = 0, \qquad (101) $$

hence

$$ \frac{\partial S}{\partial C_{IJ}} S^{-1} + S \frac{\partial (S^{-1})}{\partial C_{IJ}} = 0, \qquad (102) $$

and

$$ \frac{\partial S}{\partial C_{IJ}} = -S \frac{\partial (S^{-1})}{\partial C_{IJ}} S. \qquad (103) $$
Therefore

$$ \frac{\partial J}{\partial C_{IJ}} = - \int_{0}^{T} \left\langle S \frac{\partial (S^{-1})}{\partial C_{IJ}} S \partial_{t} w(\rho, C), \lambda(\rho, C) \right\rangle dt. \qquad (104) $$

Using the symmetry of S we further obtain

$$ \frac{\partial J}{\partial C_{IJ}} = - \int_{0}^{T} \left\langle \frac{\partial (S^{-1})}{\partial C_{IJ}} S \partial_{t} w(\rho, C), S \lambda(\rho, C) \right\rangle dt. \qquad (105) $$

Recalling that

$$ S = \begin{pmatrix} \rho I_3 & 0 \\ 0 & C^{-1} \end{pmatrix}, \qquad (106) $$

we have

$$ S^{-1} = \begin{pmatrix} \frac{1}{\rho} I_3 & 0 \\ 0 & C \end{pmatrix}, \qquad (107) $$

and

$$ \frac{\partial S^{-1}}{\partial C_{IJ}} S = \begin{pmatrix} 0 & 0 \\ 0 & \frac{\partial C}{\partial C_{IJ}} C^{-1} \end{pmatrix}. \qquad (108) $$

Hence

$$ \frac{\partial J}{\partial C_{IJ}} = - \int_{0}^{T} \left\langle \frac{\partial C}{\partial C_{IJ}} C^{-1} \partial_{t} \sigma(\rho, C), C^{-1} \lambda_{\sigma}(\rho, C) \right\rangle dt. \qquad (109) $$

Recalling that $\sigma(\rho, C)(\mathbf{x}, 0) = 0$ and $\lambda(\rho, C)(\mathbf{x}, T) = 0$, we have, by integration by parts in time

$$ \frac{\partial J}{\partial C_{IJ}} = \int_{0}^{T} \left\langle \frac{\partial C}{\partial C_{IJ}} C^{-1} \sigma(\rho, C), C^{-1} \partial_{t} \lambda_{\sigma}(\rho, C) \right\rangle dt. \qquad (110) $$

In addition, from the constitutive elastodynamics equations on the incident and adjoint fields, inside the domain of interest where the sponge layer coefficients $b(\mathbf{x})$ are zero, we have

$$ \sigma(\rho, C) = C D^{T} u(\rho, C), \qquad (111) $$

and

$$ \partial_t \lambda_\sigma(\rho, C) = C D^T \lambda_v(\rho, C). \qquad (112) $$

Therefore we have

$$ \frac{\partial J}{\partial C_{IJ}} = \int_{0}^{T} \left\langle \frac{\partial C}{\partial C_{IJ}} D^{T} u(\rho, C), D^{T} \lambda_{v}(\rho, C) \right\rangle dt, \qquad (113) $$

and we thus obtain that

$$ \frac{\partial J}{\partial C_{IJ}} = \int_{0}^{T} \left\langle \frac{\partial C}{\partial C_{IJ}} \epsilon(\rho, C), \lambda_{\epsilon}(\rho, C) \right\rangle dt, \qquad (114) $$

where we introduce the strain field $\epsilon$ and the adjoint strain field $\lambda_\epsilon$ as

$$ \epsilon = D^T u, \quad \lambda_{\epsilon} = D^T \lambda_v. \qquad (115) $$

The formulas (100) and (114) are consistent with conventional formulas derived in Tromp et al. (2005), Yang et al. (2016), and Trinh et al. (2019). In these studies, the gradient with respect to the density is related to the time-domain zero-lag cross-correlation of the incident acceleration $\partial_{tt}u$ and the adjoint displacement field. The gradient with respect to the stiffness tensor coefficients $C_{IJ}$ are related to the time-domain zero-lag cross-correlation of the incident and adjoint strain fields.

# 4 THE ADJOINT STATE METHOD FOR VISCO-ELASTIC FWI

## 4.1 Incident field equations in the visco-elastic approximation

The subsurface is an attenuating medium for the propagation of mechanical waves. Taking into account the attenuation of the subsurface in seismology is a long-standing problem, with various attenuation models proposed, from simple ones, such as the power-law model (Strick 1967; Azimi et al. 1968; Kjartansson 1979), the Kolski-Futterman model (Kolsky 1956; Futterman 1962), the Maxwell body model, the Kelvin-Voigt model (Casula & Carcione 1992), to more evolved ones such as the Zener body model, also known as Standard Linear Solid (SLS) model (Ben-Menahem & Singh 1981; Pipkin 1986), or generalization (i.e. superposition) of the latter (Carcione et al. 1988a, b). Such
generalized SLS models have been considered in seismology to model a quality factor constant over a given frequency band. This corresponds to observations made in seismology over a wide range of frequency bands (Caputo 1967; Liu et al. 1976; Kjartansson 1979). The use of fractional derivatives has also been proposed though less practical as it yields difficulties in terms of numerical solution of associated PDE (see Zhu & Harris 2014 for instance).

For FWI, depending on the environment, taking into account the attenuation in the modelling can be crucial to generate reliable estimate of the subsurface mechanical parameters. In this study, we focus on a formulation derived from the use of a superposition of SLS mechanisms. Using several such mechanisms, we can model frequency-dependent quality factors $Q_{I,J}(\mathbf{x}, \omega)$ (where $\omega$ denotes the frequency) for each of the coefficient of the stiffness tensor $C_{I,J}(\mathbf{x})$. Choosing appropriately the reference frequencies of the SLS, we can fit specific frequency dependence for each of these quality factors. The general assumption in seismology is that the quality factors $Q_{I,J}$ are constant in frequency. Such frequency independence can be fitted by the SLS, but it should be kept in mind that more complex frequency laws can also be fitted with several SLS, which could be found adequate for application of FWI outside the scope of seismology, for instance in medical imaging.

From a mathematical standpoint, each SLS mechanism corresponds to the addition of one ordinary differential equation (ODE) in the original system of elastodynamics equations. The introduction of the corresponding auxiliary variables destroys the self-adjointness property of the elastodynamics system. The adjoint system can however be derived and shown to be stable through the energy analysis of the Cauchy problem. Using L-mechanisms, the visco-elastic system of equations writes

$$ egin{cases} \rho \partial_t v - D(\nabla)\sigma = F_v, \\ \partial_t \sigma - C D(\nabla)^T v - \sum_{l=1}^L C^a \psi_l = F_\sigma, \\ \partial_t \psi_l - D(\nabla)^T v + \omega_l y_l \psi_l = 0, \quad l = 1, \dots, L, \end{cases} \quad (116) $$

where $C^a$ is a 6-by-6 matrix representing the attenuation, such that component-wise we have

$$ C_{I,J}^{a}(\mathbf{x}) = C_{I,J}(\mathbf{x}) Q_{I,J}^{-1}(\mathbf{x}), \quad I, J = 1, \dots, 6, \quad (117) $$

with $Q_{I,J}(\mathbf{x})$ the quality factor associated with the component $I, J$ of the stiffness tensor. The variables $\psi_l \in L^2(\Omega \times [0, T])^6$ are the auxiliary variables associated with the SLS mechanisms, while $y_l \in R$ are scalar calibration parameters which are tuned to ensure the fit of a specific frequency law in the frequency band spanned by the reference frequencies $\omega_l \in R$. This formulation of visco-elastic with SLS mechanisms relies on a specific approximation based on the assumption that $Q \gg 1$ (that might not be fulfilled in very specific cases) and on a simplification proposed by Yang et al. (2016) to obtain constant coefficients $y_l$ independent of the spatial position.

As in the elastic case, the visco-elastic system (116) can be rewritten as a general first-order hyperbolic system

$$ S(m)\partial_t w - A_x \partial_x w - A_y \partial_y w - A_z \partial_z w + Bw = F, \quad (118) $$

with

$$ w = \left[ v \sigma \underbrace{\psi_1 \dots \psi_L}_{L} \right]^T, \quad F = \left[ F_v C^{-1} F_\sigma \underbrace{0 \dots 0}_{L} \right]^T, \quad (119) $$

and

$$ S = \begin{pmatrix} \rho I_3 & 0 & 0 & \dots & 0 \\ 0 & C^{-1} & 0 & \dots & 0 \\ 0 & 0 & I_6 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \dots & I_6 \end{pmatrix}, \quad B = \begin{pmatrix} b I_3 & 0 & 0 & \dots & 0 \\ 0 & b I_6 & -C^{-1}C^a & \dots & -C^{-1}C^a \\ 0 & 0 & \omega_1 y_1 I_6 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & \dots & \omega_L y_L I_6 \end{pmatrix}, \quad (120) $$

and for $i \in \{x, y, z\}$,

$$ A_i = \begin{pmatrix} 0 & D_i & 0 & \dots & 0 \\ D_i^T & 0 & 0 & \dots & 0 \\ D_i^T & 0 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ D_i^T & 0 & 0 & \dots & 0 \end{pmatrix}. \quad (121) $$

Note that we have already included the sponge layers through the definition of the zero-order operator $B$ in eq. (118). We keep the same zero initial conditions $w(\mathbf{x}, 0) = 0$ as for the elastic case. For the boundary conditions, we also keep the same transparent boundary conditions as in eq. (26). The variables $\psi_l$ associated with the SLS mechanisms are defined through ordinary differential equations, therefore they do not require boundary conditions to be specified. As a result, the boundary operators in the visco-elastic case are defined by

$$ B_i^{\pm} = [\pm M_i D_i \underbrace{0_6 \dots 0_6}_{L}], \quad i \in \{x, y\}, \quad B_z^{-} = [0_3 D_z \underbrace{0_6 \dots 0_6}_{L}], \quad B_z^{+} = [M_z D_z \underbrace{0_6 \dots 0_6}_{L}]. \quad (122) $$
The visco-elastic equations thus appear in condensed form as for the elastic eqs (30), with the definitions of the variable w and operators S(m), B(m), F(m), Aᵢ and Bᵢ, i ∈ {x, y, z} overridden by eqs (119), (120), (121) and (122).

In developed form, the visco-elastic equations can be written as

$$
\left\{
\begin{aligned}
& \rho \partial_t v - D(\nabla)\sigma = F_v, \\
& \partial_t \sigma - C D(\nabla)^T v - \sum_{l=1}^{L} C^a \psi_l = F_\sigma, \\
& \partial_t \psi_l - D(\nabla)^T v + \omega_l y_l \psi_l = 0, \quad l = 1, \dots, L, \\
& v(\mathbf{x}, 0) = 0, \quad \sigma(\mathbf{x}, 0) = 0, \quad \psi_l(\mathbf{x}, 0) = 0, \quad l = 1, \dots, L, \\
& (D_x \sigma)(0, y, z, t) - M_x(0, y, z, t)v(0, y, z, t) = 0, \\
& (D_x \sigma)(L_x, y, z, t) + M_x(L_x, y, z, t)v(L_x, y, z, t) = 0, \\
& (D_y \sigma)(x, 0, z, t) - M_y(x, 0, z, t)v(x, 0, z, t) = 0, \\
& (D_y \sigma)(x, L_y, z, t) + M_y(x, L_y, z, t)v(x, L_y, z, t) = 0, \\
& (D_z \sigma)(x, y, 0, t) = 0, \\
& (D_z \sigma)(x, y, L_z, t) + M_z(x, y, L_z, t)v(x, y, L_z, t) = 0.
\end{aligned}
\right.
\qquad (123)
$$

In the following we analyse the impact of modifying the incident field equations for taking into account attenuation in the derivation of the adjoint state equations.

## 4.2 Visco-elastic FWI and adjoint state equations

Using the condensed form of the visco-elastic incident eq. (30), the Lagrangian operator given in eq. (45) remains valid. From the Lagrangian operator, we perform the same derivation as in the elastic-case of the adjoint equations. We obtain

$$
\left\{
\begin{array}{l}
S(m)\partial_t \lambda - A_x^T \partial_x \lambda - A_y^T \partial_y \lambda - A_z^T \partial_z \lambda - B^T \lambda = F_\lambda, \\
S(m)\lambda(\mathbf{x}, T) = 0, \\
A_x^T \lambda(0, y, z, t) = -(B_x^-)^T(0, y, z)v_x^-, \\
A_x^T \lambda(L_x, y, z, t) = (B_x^+)^T(L_x, y, z)v_x^+, \\
A_y^T \lambda(x, 0, z, t) = -(B_y^-)^T(x, 0, z)v_y^-, \\
A_y^T \lambda(x, L_y, z, t) = (B_y^+)^T(x, L_y, z)v_y^+, \\
A_z^T \lambda(x, y, 0, t) = -(B_z^-)^T(x, y, 0)v_z^-, \\
A_z^T \lambda(x, y, L_z, t) = (B_z^+)^T(x, y, L_z)v_z^+
\end{array}
\right.
\qquad (124)
$$

with the additional condition

$$ \mu = S(m)\lambda(\mathbf{x}, 0). \qquad (125) $$

Eq. (124) differ only from eq. (69) (corresponding to the purely elastic case) by the transpose $^T$ on the operators $A_i$, $i \in \{x, y, z\}$ and $B$, as these operators are not symmetric anymore in the viscous case.

In the following, we develop the eqs (124) and (125) to compare the incident and adjoint field equations in the visco-elastic case, as we did previously in the purely elastic case. We start with re-expressing the boundary conditions in the x direction for the adjoint equations.

Based on the definition of $A_x$, and using the notation $\lambda = \left[ \lambda_v \lambda_\sigma \underbrace{\lambda_{\psi,1} \dots \lambda_{\psi,L}}_L \right]^T$, we have for $x = 0$

$$
\left\{
\begin{array}{l}
D_x \lambda_\sigma(0, y, z, t) + \sum_{l=1}^{L} D_x \lambda_{\psi,l}(0, y, z, t) = M_x(0, y, z)v_x^-(y, z, t), \\
D_x^T \lambda_v(0, y, z, t) = -D_x^T v_x^-(y, z, t),
\end{array}
\right.
\qquad (126)
$$

which gives

$$ D_x \left( \lambda_\sigma + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (0, y, z, t) + M_x(0, y, z) \lambda_v(0, y, z, t) = 0. \qquad (127) $$

For $x = L_x$, we obtain

$$
\left\{
\begin{array}{l}
D_x \lambda_\sigma(L_x, y, z, t) + \sum_{l=1}^{L} D_x \lambda_{\psi,l}(L_x, y, z, t) = M_x(L_x, y, z)v_x^+(y, z, t), \\
D_x^T \lambda_v(L_x, y, z, t) = D_x^T v_x^+(y, z, t),
\end{array}
\right.
\qquad (128)
$$
which gives

$$D_x \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (L_x, y, z, t) - M_x(L_x, y, z)\lambda_v(L_x, y, z, t) = 0. \quad (129)$$

Similarly we have for $y = 0$

$$D_y \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, 0, z, t) + M_y(x, 0, z)\lambda_v(x, 0, z, t) = 0, \quad (130)$$

and for $y = L_y$

$$D_y \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, L_y, z, t) - M_y(x, L_y, z)\lambda_v(x, L_y, z, t) = 0. \quad (131)$$

For $z = 0$ we have

$$\left\{ \begin{array}{l} D_z \lambda_{\sigma}(x, y, 0, t) + \sum_{l=1}^{L} D_z \lambda_{\psi,l}(x, y, 0, t) = 0, \\ D_z^T \lambda_v(x, y, 0, t) = -D_x^T \nu_z^{-}(x, y, t), \end{array} \right. \quad (132)$$

and for $z = L_z$ we have

$$D_z \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, y, L_z, t) - M_z(x, y, L_z)\lambda_v(x, y, L_z, t) = 0. \quad (133)$$

The boundary conditions in the visco-elastic case can thus be summarized as

$$\begin{align} D_x \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (0, y, z, t) + M_x(0, y, z)\lambda_v(0, y, z, t) &= 0, \nonumber \\ D_x \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (L_x, y, z, t) - M_x(L_x, y, z)\lambda_v(L_x, y, z, t) &= 0, \nonumber \\ D_y \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, 0, z, t) + M_y(x, 0, z)\lambda_v(x, 0, z, t) &= 0, \tag{134} \\ D_y \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, L_y, z, t) - M_y(x, L_y, z)\lambda_v(x, L_y, z, t) &= 0, \nonumber \\ D_z \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, y, 0, t) &= 0, \nonumber \\ D_z \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, y, L_z, t) - M_z(x, y, L_z)\lambda_v(x, y, L_z, t) &= 0, \nonumber \end{align}$$

which are the same as in the purely elastic case if we replace the adjoint stress $\lambda_{\sigma}$ by the “relaxed” adjoint stress $\lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l}$.

The adjoint system in the visco-elastic case can thus be rewritten as

$$\left\{ \begin{array}{l} S(m)\partial_t \lambda - A_x^T \partial_x \lambda - A_y^T \partial_y \lambda - A_z^T \partial_z \lambda - B^T \lambda = F_{\lambda}, \\ S(m)\lambda(\mathbf{x}, T) = 0, \\ D_x \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (0, y, z, t) + M_x(0, y, z)\lambda_v(0, y, z, t) = 0, \\ D_x \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (L_x, y, z, t) - M_x(L_x, y, z)\lambda_v(L_x, y, z, t) = 0, \\ D_y \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, 0, z, t) + M_y(x, 0, z)\lambda_v(x, 0, z, t) = 0, \\ D_y \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, L_y, z, t) - M_y(x, L_y, z)\lambda_v(x, L_y, z, t) = 0, \\ D_z \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, y, 0, t) = 0, \\ D_z \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, y, L_z, t) - M_z(x, y, L_z)\lambda_v(x, y, L_z, t) = 0, \end{array} \right. \quad (135)$$
with the additional conditions

$$
\left\{
\begin{array}{l}
\mu = S(m)\lambda(\mathbf{x}, 0), \\
\nu_x^-(y, z, t) = -\lambda_v(0, y, z, t), \\
\nu_x^+(y, z, t) = \lambda_v(L_x, y, z, t), \\
\nu_y^-(x, z, t) = -\lambda_v(x, 0, z, t), \\
\nu_y^+(x, z, t) = \lambda_v(x, L_y, z, t), \\
\nu_z^-(x, y, t) = -\lambda_v(x, y, 0, t), \\
\nu_y^+(x, y, t) = \lambda_v(x, y, L_z, t).
\end{array}
\right.
\qquad (136)
$$

As in the elastic case, we can perform a change of variable $t' = T - t$ to change a backward propagating system into a forward propagating system. Dropping the superscript ' for the sake of concision, we obtain the following visco-elastic adjoint system, in condensed form

$$
\left\{
\begin{array}{l}
S(m)\partial_t \lambda + A_x^T \partial_x \lambda + A_y^T \partial_y \lambda + A_z^T \partial_z \lambda + B^T \lambda = -F_{\lambda}, \\
S(m)\lambda(\mathbf{x}, 0) = 0, \\
D_x \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (0, y, z, t) + M_x(0, y, z)\lambda_v(0, y, z, t) = 0, \\
D_x \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (L_x, y, z, t) - M_x(L_x, y, z)\lambda_v(L_x, y, z, t) = 0, \\
D_y \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, 0, z, t) + M_y(x, 0, z)\lambda_v(x, 0, z, t) = 0, \\
D_y \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, L_y, z, t) - M_y(x, L_y, z)\lambda_v(x, L_y, z, t) = 0, \\
D_z \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, y, 0, t) = 0, \\
D_z \left( \lambda_{\sigma} + \sum_{l=1}^{L} \lambda_{\psi,l} \right) (x, y, L_z, t) - M_z(x, y, L_z)\lambda_v(x, y, L_z, t) = 0,
\end{array}
\right.
\qquad (137)
$$

and in developed form

$$
\left\{
\begin{array}{l}
\rho \partial_t \lambda_v + D(\nabla) \left( \lambda_\sigma + \sum_{l=1}^L \lambda_{\psi,l} \right) + B_v \lambda_v = -F_{\lambda,v}, \\
C^{-1} \partial_t \lambda_\sigma + D(\nabla)^T \lambda_v + B_\sigma \lambda_\sigma = 0, \\
\partial_t \lambda_{\psi,l} + C^{-1} C_a \lambda_\sigma + \omega_l y_l \lambda_{\psi,l} = 0, \quad l=1, \dots, L, \\
\lambda_v(\mathbf{x}, 0) = 0, \quad \lambda_\sigma(\mathbf{x}, 0) = 0, \quad \lambda_{\psi,l}(\mathbf{x}, 0) = 0, \quad l=1, \dots, L, \\
D_x \left( \lambda_\sigma + \sum_{l=1}^L \lambda_{\psi,l} \right) (0, y, z, t) + M_x(0, y, z) \lambda_v(0, y, z, t) = 0, \\
D_x \left( \lambda_\sigma + \sum_{l=1}^L \lambda_{\psi,l} \right) (L_x, y, z, t) - M_x(L_x, y, z) \lambda_v(L_x, y, z, t) = 0, \\
D_y \left( \lambda_\sigma + \sum_{l=1}^L \lambda_{\psi,l} \right) (x, 0, z, t) + M_y(x, 0, z) \lambda_v(x, 0, z, t) = 0, \\
D_y \left( \lambda_\sigma + \sum_{l=1}^L \lambda_{\psi,l} \right) (x, L_y, z, t) - M_y(x, L_y, z) \lambda_v(x, L_y, z, t) = 0, \\
D_z \left( \lambda_\sigma + \sum_{l=1}^L \lambda_{\psi,l} \right) (x, y, 0, t) = 0, \\
D_z \left( \lambda_\sigma + \sum_{l=1}^L \lambda_{\psi,l} \right) (x, y, L_z, t) - M_z(x, y, L_z) \lambda_v(x, y, L_z, t) = 0,
\end{array}
\right.
\qquad (138)
$$

with the additional constraints in eq. (136).

The visco-elastic system is not self-adjoint: the equations differ from the incident field eq. (123). The boundary conditions are imposed on the quantity $\lambda_\sigma + \sum_{l=1}^L \lambda_{\psi,l}$ which are a combination of the adjoint stress and the adjoint memory variables. However, this is the quantity which is involved in the spatial derivatives of the equation on the adjoint velocities: therefore it makes sense that the boundary condition is defined on this quantity, as it is required to evaluate the derivatives on the first point for the velocities. Note also that the adjoint system is different from the incident field equation by the localization of the memory variables: they are acting on the equation on the adjoint velocities, instead of the stress. It can be found in the literature however that in terms of attenuation implementation, the two formulations are equivalent.
Regarding FWI code implementation, this has an important implication: a specific adjoint system needs to be implemented. A code dedicated to the solution of the incident field eq. (123) could not be re-used straightforwardly for the adjoint system. This is a strong difference with the purely elastic case.

## 4.3 Gradient expression

The gradient expression for the density is the same as in the purely elastic case. Modifications occur for the gradient with respect to the $C_{IJ}$. The visco-elastic approximation also offers the possibility to invert for the attenuation factors $Q_{IJ}$ and therefore to compute the corresponding gradient.

Let us start with the gradient with respect to the $C_{IJ}$. Based on the Lagrangian definition in the visco-elastic case we have

$$ \frac{\partial J}{\partial C_{IJ}} = \int_{0}^{T} \left\langle \frac{\partial S}{\partial C_{IJ}} \partial_{t} w, \lambda \right\rangle dt + \int_{0}^{T} \left\langle \frac{\partial B}{\partial C_{IJ}} w, \lambda \right\rangle dt. \quad (139) $$

Compared with the pure elastic case, we have a contribution from the zero-order term $B$ coming in $\partial_C B$, which originates from the introduction of the memory variable terms.

From the definition of the matrix $B$ we have

$$ \int_{0}^{T} \left\langle \frac{\partial B}{\partial C_{IJ}} w, \lambda \right\rangle dt = - \int_{0}^{T} \left\langle \sum_{l=1}^{L} \frac{\partial(C^{-1}C^{a})}{C_{IJ}} \psi_{l}, \lambda_{\sigma} \right\rangle dt. \quad (140) $$

Therefore, re-exploiting the derivation from the previous elastic-case we have

$$ \frac{\partial J}{\partial C_{IJ}} = \int_{0}^{T} \left\langle \frac{\partial C^{-1}}{C_{IJ}} \partial_{l} \sigma, \lambda_{\sigma} \right\rangle dt - \int_{0}^{T} \left\langle \sum_{l=1}^{L} \frac{\partial(C^{-1}C^{a})}{C_{IJ}} \psi_{l}, \lambda_{\sigma} \right\rangle dt. \quad (141) $$

Re-arranging the terms, we obtain

$$ \frac{\partial J}{\partial C_{IJ}} = \int_{0}^{T} \left\langle \frac{\partial C^{-1}}{C_{IJ}} \left( \partial_{l} \sigma - \sum_{l=1}^{L} C^{a} \psi_{l} \right), \lambda_{\sigma} \right\rangle dt - \int_{0}^{T} \left\langle \sum_{l=1}^{L} C^{-1} \frac{\partial C^{a}}{C_{IJ}} \psi_{l}, \lambda_{\sigma} \right\rangle dt. \quad (142) $$

Based on the constitutive visco-elastic equations, the first term of the right hand side can be simplified by replacing the stress and the memory variables by the velocity, such that

$$ \frac{\partial J}{\partial C_{IJ}} = \int_{0}^{T} \left\langle \frac{\partial C^{-1}}{C_{IJ}} C D^{T} v, \lambda_{\sigma} \right\rangle dt - \int_{0}^{T} \left\langle \sum_{l=1}^{L} C^{-1} \frac{\partial C^{a}}{C_{IJ}} \psi_{l}, \lambda_{\sigma} \right\rangle dt. \quad (143) $$

We recognize in the first term of the right hand side the gradient with respect to $C_{IJ}$ in the purely elastic case, therefore we obtain

$$ \frac{\partial J}{\partial C_{IJ}} = \int_{0}^{T} \left\langle \frac{\partial C}{\partial C_{IJ}} \epsilon, \lambda_{\epsilon} \right\rangle dt - \int_{0}^{T} \left\langle \sum_{l=1}^{L} C^{-1} \frac{\partial C^{a}}{C_{IJ}} \psi_{l}, \lambda_{\sigma} \right\rangle dt. \quad (144) $$

This expression shows how the visco-elastic approximation brings a correction to the gradient expression with respect to $C_{IJ}$ which depends on the memory variables and the quality factors through the term $\frac{\partial C^a}{C_{IJ}}$.

Finally, we obtain the gradient with respect to the quality factor $Q_{IJ}$ in a similar way. We have

$$ \frac{\partial J}{\partial Q_{IJ}} = \int_{0}^{T} \left\langle \frac{\partial B}{\partial Q_{IJ}} w, \lambda \right\rangle dt. \quad (145) $$

Therefore

$$ \frac{\partial J}{\partial Q_{IJ}} = - \int_{0}^{T} \left\langle \sum_{l=1}^{L} C^{-1} \frac{\partial C^{a}}{Q_{IJ}} \psi_{l}, \lambda_{\sigma} \right\rangle dt. \quad (146) $$

# 5 DISCUSSION

## 5.1 General misfit function

FWI, from a mathematical standpoint, is an ill-posed nonlinear inverse problem. Alternative formulations to FWI have been proposed to mitigate this nonlinearity and make the numerical solution based on local optimization more robust, less dependent on specific workflow and on the initial model. Among these reformulations, numerous attempts consist in modifying the least-squares misfit function into a more generic comparison function (see Pladys et al. 2021 for instance). In this case, we would adopt the general notations

$$ \min_{m} J(m), \quad J(m) = J(Rw[m], d_{\text{obs}}), \quad (147) $$
where R extracts the value of the velocity field at receivers positions

$$ (Rw)(\mathbf{x}_r, t) = \int_{\Omega} \delta(\mathbf{x} - \mathbf{x}_r) v(\mathbf{x}, t) dx = [v_1(\mathbf{x}_r, t) v_2(\mathbf{x}_r, t) v_3(\mathbf{x}_r, t)]^T, \quad (148) $$

and J is a general comparison function between $Rw$ and $d_{\text{obs}}$. Note that we still consider the case of a single source here, the extension to the case of multiple sources remaining straightforward by summation over sources. Replacing the term $\mathcal{A}$ in the Lagrangian definition (45) by the expression in eqs (147) and (148), and applying the same strategy of derivation of the Lagrangian as the one described in eqs (56), (57), (58) and (59), leads to consider the expansion

$$ J(R(w+q), d_{\text{obs}}) = J(Rw + Rq, d_{\text{obs}}) = J(Rw, d_{\text{obs}}) + (\partial_1 J(Rw, d_{\text{obs}}), Rq) + O(\|q\|^2), \quad (149) $$

where $\partial_1$ denotes the derivatives of J with respect to its first argument. Expanding the second term, we obtain

$$ (\partial_1 J(Rw, d_{\text{obs}}), Rq) = \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} \partial_1 J(Rw, d_{\text{obs}})_i(\mathbf{x}_r, t) (Rq)_i(\mathbf{x}_r, t) dt. \quad (150) $$

Thus

$$ (\partial_1 J(Rw, d_{\text{obs}}), Rq) = \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} \partial_1 J(Rw, d_{\text{obs}})_i(\mathbf{x}_r, t) \int_{\Omega} q_i(x, t) \delta(\mathbf{x} - \mathbf{x}_r) dt, \quad (151) $$

and

$$ (\partial_1 J(Rw, d_{\text{obs}}), Rq) = \int_{\Omega} \int_{0}^{T} \sum_{i=1}^{3} \sum_{r=1}^{N_r} \partial_1 J(Rw, d_{\text{obs}})_i(\mathbf{x}_r, t) \delta(\mathbf{x} - \mathbf{x}_r) q_i(\mathbf{x}, t) dt, \quad (152) $$

which leads to

$$ (\partial_1 J(Rw, d_{\text{obs}}), Rq) = (q, F_{\lambda}), \quad (153) $$

where

$$ F_{\lambda,i}(\mathbf{x}, t) = \sum_{r=1}^{N_r} \partial_1 J(Rw, d_{\text{obs}})_i(\mathbf{x}_r, t) \delta(\mathbf{x} - \mathbf{x}_r). \quad (154) $$

We thus observe that in the case of a general misfit function J based on the comparison of particle velocities, the source for the adjoint fields still consists in oriented forces localized at each receiver position $\mathbf{x}_r$. The time signature of these forces depend on the first-order derivative of the misfit function comparison J with respect to calculated data $Rw$.

## 5.2 Gradient with respect to initial and boundary conditions

As mentioned in the introduction, one of the byproduct of the general derivation presented here is the ability to derive formulas for gradient with respect to the initial and boundary conditions, which might be useful in specific contexts. This derivation is straightforward. Let us take for instance the example of the initial condition. Let us denote this initial conditions by

$$ w(\mathbf{x}, 0) = w_0(\mathbf{x}). \quad (155) $$

The Lagrangian of the elastic system now becomes

$$ \begin{aligned} \mathcal{L}(w_0, w, \lambda, \mu, v) = & \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} \left[ |(R_i w)(\mathbf{x}_r, t) - d_{\text{obs},i}(\mathbf{x}_r, t)|^2 dt \right. \\ & \left. + (S(m)\partial_t w - A_x \partial_x w - A_y \partial_y w - A_z \partial_z w + Bw - F, \lambda)_{\Omega \times [0, T]} + \underbrace{(w(\mathbf{x}, 0) - w_0, \mu)_{\Omega}}_{\text{modified term}} \right] \\ & + \left( (B_x w)(0, y, z, t), v_x^{-})_{\Omega_y \times \Omega_z \times [0, T]} + \left( (B_x w)(L_x, y, z, t), v_x^{+})_{\Omega_y \times \Omega_z \times [0, T]} \right. \right. \\ & \left. + \left( (B_y w)(x, 0, z, t), v_y^{-})_{\Omega_x \times \Omega_z \times [0, T]} + \left( (B_y w)(x, L_y, z, t), v_y^{+})_{\Omega_x \times \Omega_z \times [0, T]} \right. \right. \\ & \left. + \left( (B_z w)(x, y, 0, t), v_z^{-})_{\Omega_x \times \Omega_y \times [0, T]} + \left( (B_z w)(x, y, L_z, t), v_z^{+})_{\Omega_x \times \Omega_y \times [0, T]} \right. \right. \end{aligned} \quad (156) $$

In a similar way as for the parameters m, we have

$$ \frac{\partial J}{\partial w_0} = \partial_{w_0} \mathcal{L}(w_0, w(w_0), \lambda(w_0), \mu(w_0), v(w_0)), \quad (157) $$

where $w(w_0)$, $\lambda(w_0)$, $\mu(w_0)$, $v(w_0)$ are solution of the incident field eqs (31) and adjoint field eqs (81). Therefore, we have

$$ \frac{\partial J}{\partial w_0} = -\mu(w_0), \quad (158) $$
where

$$ \mu(w_0) = S(m)\lambda(\mathbf{x}, 0), \qquad (159) $$

with $\lambda$ solution of the adjoint field eqs (81).

Similarly, we could be interested in the gradient with respect to a boundary condition in $x = 0$ such that

$$ B_x w(0, y, z, t) = \phi(y, z, t). \qquad (160) $$

In this case, the Lagrangian operator becomes

$$ \begin{aligned} \mathcal{L}(\phi, w, \lambda, \mu, \nu) = & \sum_{i=1}^{3} \sum_{r=1}^{N_r} \int_{0}^{T} |(R_i w)(x_r, t) - d_{\text{obs},i}(x_r, t)|^2 dt \\ & + (S(m)\partial_t w - A_x \partial_x w - A_y \partial_y w - A_z \partial_z w + Bw - F, \lambda)_{\Omega \times [0, T]} + (w(\mathbf{x}, 0), \mu)_{\Omega} \\ & + \underbrace{((B_x w)(0, y, z, t) - \phi, v_x^-)_{\Omega_y \times \Omega_z \times [0, T]}}_{\text{modified term}} + \underbrace{((B_x w)(L_x, y, z, t), v_x^+)_{\Omega_y \times \Omega_z \times [0, T]}}_{\text{modified term}} \\ & + ((B_y w)(x, 0, z, t), v_y^-)_{\Omega_x \times \Omega_z \times [0, T]} + ((B_y w)(x, L_y, z, t), v_y^+)_{\Omega_x \times \Omega_z \times [0, T]} \\ & + ((B_z w)(x, y, 0, t), v_z^-)_{\Omega_x \times \Omega_y \times [0, T]} + ((B_z w)(x, y, L_z, t), v_z^+)_{\Omega_x \times \Omega_y \times [0, T]}, \end{aligned} \qquad (161) $$

and we have

$$ \frac{\partial J}{\partial \phi} = \partial_{\phi} \mathcal{L}(\phi, w(\phi), \lambda(\phi), \mu(\phi), \nu(\phi)), \qquad (162) $$

where $w(\phi)$, $\lambda(\phi)$, $\mu(\phi)$, $\nu(\phi)$ are solution of the incident field eqs (31) and adjoint field eqs (81). Therefore, we have

$$ \frac{\partial J}{\partial \phi} = -v_x^- (\phi), \qquad (163) $$

where

$$ v_x^- (\phi) = -\lambda_v (0, y, z, t), \qquad (164) $$

with $\lambda$ solution of the adjoint field eqs (81).

More generally, we see that the complementary eqs (82) and (83) give the adjoint quantities which are necessary to express the gradient of the misfit function with respect to the initial or boundary conditions.

## 5.3 Discretization

We shall end this discussion with a comment on a recurring debate among mathematician and geophysicists having worked in this field. This debate is related to the actual implementation of the formulas obtained by the adjoint state technique to compute the FWI gradient. Two possibilities can be considered. The first consist in discretizing directly the formulas for the incident field, adjoint field and gradient formula obtained following the development presented here. This is called the Optimize then Discretize approach. The second consists in discretizing first the incident field equations, then introducing discrete operators and reproduce the same development as the one presented here however using this discrete settings. This is called the Discretize then Optimize approach (see for instance Betts & Campbell 2005).

From a mathematical standpoint, the Discretize then Optimize approach is the one to be chosen. Indeed, this is the only consistent way of obtaining the adjoint field and gradient formulas associated with the discrete incident field which is computed on machine. This way, orders of convergence can be guaranteed for the incident field, adjoint field and the FWI gradient.

However, in practice, the mathematical derivation of the discrete adjoint field and gradient equations is far from being straightforward. An example of such derivation is provided in the 2D acoustic case in Métivier (2009). Already in this relatively simple settings the derivation of the adjoint field equation is far from being trivial. Another consequence of the Discretize then Optimize approach is that the adjoint and gradient equations need to be derived for each discretization scheme of the incident field equations. This constitutes a severe genericity loss for the adjoint approach. Finally, the numerical error introduced by the Optimize then Discretize approach, that is employing a discretization of the adjoint field and gradient equations which might not be consistent with the incident field discretization, is in practice “small”. More specifically, the potential error introduced at this stage is usually below the modelling and instrumental noise levels, which makes this additional numerical error negligible in front of the other error sources in field data applications. In practice, linesearch based optimization schemes require to be fed with descent direction, that are model updates “aligned” with the steepest descent direction: mathematically speaking, directions in the model space for which the scalar product with the gradient is strictly negative (Nocedal & Wright 2006). This is something which is easily guaranteed in practice, despite the various sources of errors. This is the reason why the Optimize then Discretize option is most of the time the one implemented in large 3-D FWI codes actually applied to field data.

A final note might temper this message: the Discretize then Optimize approach can be implemented “automatically” by the usage of automatic differentiation tools (Baydin et al. 2018). These tools directly generate the implementation of the adjoint function of any function
implemented in a computer language which support these tools. These tools rely on an interpretation of any computing code as the composition of elementary algebraic operations (summation, multiplication) for which the adjoint operation is trivial to compute. By chain-rule, it is then possible to derive an adjoint code, which will compute the variables which are the adjoint to the ones computed in a given initial code. The supporting formalism giving the elemental adjoint operations is fully discrete, therefore the use of such kind of automatic tool provide the adjoint quantities the most consistent from the quantities computed in the input code. Such tools are however for now limited to relatively "simple" input code, and might have difficulties in traducing efficiently highly optimized and parallelized codes.

## 6 CONCLUSION

We propose here a complete, rigorous and self-consistent derivation of the FWI gradient formula in the elastodynamics approximation using the adjoint state method. This derivation takes into account the initial and boundary conditions which are usually disregarded in the Lagrangian operator. We show that including them makes it possible to derive consistent expressions for the adjoint state equations. The final conditions for the adjoint state field appears naturally, as well as the boundary conditions. In the purely elastic case, the adjoint field is solution of the same elastic wave equation as the incident field, up to a sign change on the stress variable. In the visco-elastic case, using SLS mechanisms, the system of equations for the incident field is not anymore self adjoint. When the relaxation mechanisms are added to the incident stress equations, the adjoint relaxation mechanisms appear in the adjoint velocity equations. The adjoint boundary conditions are also changed with respect to the incident field boundary equations in the visco-elastic case. We derive the gradient expression for the density, the stiffness tensor coefficients, and the attenuation factor, and show that they are consistent with usual derivations. Finally, owing to the generality of the approach, it is possible to further derive gradient of the misfit function with respect to the initial or boundary conditions in a straightforward way. Given the importance of FWI in seismic imaging at all scales, it appears to us important that this derivation is made available to a broad audience of geophysicists.

## ACKNOWLEDGMENTS

This study was partially funded by the SEISCPE consortium (http://seiscope2.osug.fr), sponsored by AKERBP, DUG, EXXONMOBIL, GEOLINKS, JGI, PETROBRAS, SHELL, SINOPEC, TOTALENERGIES and VIRIDIEN. This study was granted access to the HPC resources provided by the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which was supported by Grenoble research communities, the HPC resources of Cray Marketing Partner Network (https://partners.cray.com), and those of CINES/IDRIS/TGCC under the allocation 046091 made by GENCI.

## DATA AVAILABILITY

No data nor computer codes are associated with this work.

## REFERENCES

* Azimi, S.A., Kalinin, A., Kalinin, V. & Pivovarov, B., 1968. Impulse and transient characteristics of media with linear and quadratic absorption laws, *Izv. Phys. Solid Earth*, **2**, 88–93.

* Baydin, A.G., Pearlmutter, B.A., Radul, A.A. & Siskind, J.M., 2018. Automatic differentiation in machine learning: a survey, *J. Mach. Learn. Res.*, **18**(153), 1–43.

* Becache, E., Fauqueux, S. & Joly, P., 2003. Stability of Perfectly Matched Layers, Group Velocities and Anisotropic Waves, *J. Comput. Phys.*, **188**, 399–433.

* Ben-Menahem, A. & Singh, S., 1981. *Seismic Waves and Sources*, Springer-Verlag New York, New York, NY.

* Bérenger, J.-P., 1994. A perfectly matched layer for absorption of electromagnetic waves, *J. Comput. Phys.*, **114**, 185–200.

* Betts, J.T. & Campbell, S.L., 2005. Discretize then optimize, *Mathematics for Industry: Challenges and Frontiers*, SIAM Society for Industrial and Applied Mathematics Philadelphia, pp. 140–157.

* Bonnans, J.F., Gilbert, J.C., Lemaréchal, C. & Sagastizábal, C.A., 2006. *Numerical Optimization, Theoretical and Practical Aspects*, Springer series, Universitext, https://link.springer.com/book/10.1007/978-3-540-35447-5.

* Bozdağ, E., Peter, D., Lefebvre, M., Komatitsch, D., Tromp, J., Hill, J., Podhorszki, N. & Pugmire, D., 2016. Global Adjoint Tomography: First-generation Model, *Geophys. J. Int.*, **207**(3), 1739–1766.

* Cao, J., Brossier, R., Górszczyk, A., Métivier, L. & Virieux, J., 2022. 3D multi-parameter full-waveform inversion for ocean-bottom seismic data

using an efficient fluid-solid coupled spectral-element solver, *Geophys. J. Int.*, **229**(1), 671–703.

* Caputo, M., 1967. Linear models of dissipation whose Q is almost frequency independent-II, *Geophys. J. Int.*, **13**(5), 529–539.

* Carcione, J., Kosloff, D. & Kosloff, R., 1988a. Wave propagation simulation in a linear viscoacoustic medium, *Geophys. J. Int.*, **93**(2), 393–401.

* Carcione, J., Kosloff, D. & Kosloff, R., 1988b. Wave propagation simulation in a linear viscoelastic medium, *Geophys. J. Int.*, **95**(3), 597–611.

* Casula, G. & Carcione, J., 1992. Generalized mechanical model analogies of linear viscoelastic behaviour, *Bollettino di Geofisica Teorica ed Applicata*, **34**(136), 235–256.

* Cerjan, C., Kosloff, D., Kosloff, R. & Reshef, M., 1985. A nonreflecting boundary condition for discrete acoustic and elastic wave equations, *Geophysics*, **50**(4), 2117–2131.

* Chavent, G., 1974. Identification of function parameters in partial differential equations, in *Identification of Parameter Distributed Systems*, pp. 31–48, eds, Goodson, R. & Polis, M., American Society of Mechanical Engineers, New York.

* Claerbout, J., 1985. *Imaging the Earth's Interior*, Blackwell Scientific Publication.

* Devaney, A., 1984. Geophysical diffraction tomography, *IEEE Trans. Geosci. Remote Sens.*, **GE-22**(1), 3–13.

* Fichtner, A. & van Driel, M., 2014. Models and Fréchet kernels for frequency-(in)dependent Q, *Geophys. J. Int.*, **198**(3), 1878–1889.
Fichtner, A. & Villaseñor, A., 2015. Crust and upper mantle of the western mediterranean-constraints from full-waveform inversion, *Earth planet. Sci. Lett.*, **428**, 52–62.

Fichtner, A., Kennett, B. L.N., Igel, H. & Bunge, H.P., 2009. Full seismic waveform tomography for upper-mantle structure in the Australasian region using adjoint methods, *Geophys. J. Int.*, **179**(3), 1703–1725.

Futterman, W., 1962. Dispersive body waves, *J. Geophys. Res.*, **67**, 5279–5291.

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

Guasch, L., Agudo, O.C., Tang, M.-X., Nachev, P. & Warner, M., 2020. Full-waveform inversion imaging of the human brain, *NPJ Digit. Med.*, **3**(1), 1–12.

Irnaka, T.M., Brossier, R., Metivier, L., Bohlen, T. & Pan, Y., 2022. 3D Multi-component Full Waveform Inversion for Shallow-Seismic Target: Ettlingen Line Case Study, *Geophys. J. Int.*, **229**(2), 1017–1040.

Karaoğlu, H. & Romanowicz, B., 2018. Inferring global upper-mantle shear attenuation structure by waveform tomography using the spectral element method, *Geophys. J. Int.*, **213**(3), 1536–1558.

Kjartansson, E., 1979. Constant Q-wave propagation and attenuation, *J. geophys. Res.*, **84**(B9), 4737–4748.

Koehn, D. *et al.*, 2018. Seismic SH full waveform inversion as new prospection method in archaeogeophysics, in *Conference Proceedings, 24th European Meeting of Environmental and Engineering Geophysics*, pp. 1–5, European Association of Geoscientists and Engineers.

Kolsky, H., 1956. The propagation of stress pulses in viscoelastic solids, *Philos. Mag.*, **1**, 693–710.

Komatitsch, D., 1997. *Méthodes spectrales et éléments spectraux pour l'équation de l'élastodynamique 2D et 3D en milieu hétérogène*, Ph.D. thesis, Institut de Géophysique du Globe de Paris.

Komatitsch, D., Xie, Z., Bozdağ, E., de Andrade, E.S., Peter, D., Liu, Q. & Tromp, J., 2016. Anelastic sensitivity kernels with parsimonious storage for adjoint tomography and full waveform inversion, *Geophys. J. Int.*, **206**(3), 1467–1478.

Kouroussis, G., Verlinden, O. & Conti, C., 2011. Finite-Dynamic Model for Infinite Media: Corrected Solution of Viscous Boundary Efficiency, *J. Eng. Mech.*, **137**(7), 509–511.

Kreiss, H.O. & Lorenz, J., 2004. *Initial-boundary Value Problems and the Navier-Stokes Equations, Classics in Applied Mathematics 47*, SIAM, Philadelphia.

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

Lei, W. *et al.*, 2020. Global adjoint tomography—model glad-m25, *Geophys. J. Int.*, **223**(1), 1–21.

Lions, J.L., 1968. Contrôle optimal de systèmes gouvernés par des équations aux d'érivées partielles, *SIAM Rev.*, **12**, 308–309.

Liu, H.-P., Anderson, D.L. & Kanamori, H., 1976. Velocity dispersion due to anelasticity; implications for seismology and mantle composition, *Geophys. J. Int.*, **47**(1), 41–58.

Lu, Y., Stehly, L., Brossier, R., Paul, A. & Group, A.W., 2020. Imaging Alpine crust using ambient noise wave-equation tomography, *Geophys. J. Int.*, **222**(1), 69–85.

Lysmer, J. & Kuhlemeyer, R., 1969. Finite Dynamic Model for Infinite Media, *J. Eng. Mech. Div.*, **95**(EM4), 859–877.

Marty, P., Boehm, C. & Fichtner, A., 2021. Acoustoelastic full-waveform inversion for transcranial ultrasound computed tomography, in *Medical Imaging 2021: Ultrasonic Imaging and Tomography*, Vol. **11602**, pp. 1160211, International Society for Optics and Photonics.

Métivier, L., 2009. *Une méthode d'inversion non linéaire pour l'imagerie sismique haute résolution*, Theses, Université Paris-Nord - Paris XIII, Fin de rédaction.

Métivier, L. & Brossier, R., 2016. The SEISCOPE optimization toolbox: A large-scale nonlinear optimization library based on reverse communication, *Geophysics*, **81**(2), F11–F25.

Métivier, L., Brossier, R., Labbé, S., Operto, S. & Virieux, J., 2014. A robust absorbing layer for anisotropic seismic wave modeling, *J. Comput. Phys.*, **279**, 218–240.

Nocedal, J., 1980. Updating Quasi-Newton Matrices With Limited Storage, *Math. Comput.*, **35**(151), 773–782.

Nocedal, J. & Wright, S.J., 2006. *Numerical Optimization*, 2nd edn, Springer.

Nouibat, A., Brossier, R., Stehly, L., Cao, J. & Paul, A. & Cifalps Team and AlpArray Working Group, 2023. Ambient-noise wave-equation tomography of the alps and ligurian-provence basin, *J. Geophys. Res.: Solid Earth*, **128**(10), e2023JB026776, doi:10.1029/2023JB026776.

Operto, S., Miniussi, A., Brossier, R., Combe, L., Métivier, L., Monteiller, V., Ribodetti, A. & Virieux, J., 2015. Efficient 3-D frequency-domain mono-parameter full-waveform inversion of ocean-bottom cable data: application to Valhall in the visco-acoustic vertical transverse isotropic approximation, *Geophys. J. Int.*, **202**(2), 1362–1391.

Pipkin, A.C., 1986. *Lectures on Viscoelasticity Theory*, Vol. 7, Springer Science and Business Media.

Pladys, A., Brossier, R., Li, Y. & Métivier, L., 2021. On cycle-skipping and misfit function modification for full-wave inversion: Comparison of five recent approaches, *Geophysics*, **86**(4), R563–R587.

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. & Perkins, C., 2010. Full waveform inversion of a deep water ocean bottom seismometer dataset, *First Break*, **28**, 71–78.

Raknes, E.B., Arntsen, B. & Weibull, W., 2015. Three-dimensional elastic full waveform inversion using seismic data from the sleipner area, *Geophys. J. Int.*, **202**(3), 1877–1894.

Sirgue, L. & Pratt, R.G., 2004. Efficient waveform inversion and imaging: a strategy for selecting temporal frequencies, *Geophysics*, **69**(1), 231–248.

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**, 65–70.

Solano, C.P. & Plessix, R.-É., 2019. Velocity-model building with enhanced shallow resolution using elastic waveform inversion—an example from onshore oman, *Geophysics*, **84**(6), R977–R988.

Stopin, A., Plessix, R.-E. & Al Abri, S., 2014. Multiparameter waveform inversion of a large wide-azimuth low-frequency land data set in Oman, *Geophysics*, **79**(3), WA69–WA77.

Strick, E., 1967. The determination of q, dynamic viscosity and transient creep curves from wave propagation measurements, *Geophys. J. Int.*, **13**(1–3), 197–218.

Tape, C., Liu, Q., Maggi, A. & Tromp, J., 2009. Adjoint tomography of the southern california crust, *Science*, **325**, 988–992.

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

Thrustarson, S., van Herwaarden, D.-P., Krischer, L., Boehm, C., van Driel, M., Afanasiev, M. & Fichtner, A., 2022. Data-adaptive global full-waveform inversion, *Geophys. J. Int.*, **230**(2), 1374–1393.

Trinh, P.T., Brossier, R., Métivier, L., Tavard, L. & Virieux, J., 2019. Efficient 3D time-domain elastic and viscoelastic Full Waveform Inversion using a spectral-element method on flexible Cartesian-based mesh, *Geophysics*, **84**(1), R75–R97.

Tromp, J., 2019. Seismic wavefield imaging of earth's interior across scales, *Nat. Rev. Earth Environ.*, **1**, 1–14.

Tromp, J., Tape, C. & Liu, Q., 2005. Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels, *Geophys. J. Int.*, **160**, 195–216.

Vigh, D., Jiao, K., Watts, D. & Sun, D., 2014. Elastic full-waveform inversion application using multicomponent measurements of seismic data collection, *Geophysics*, **79**(2), R63–R77.

Virieux, J., Asnaashari, A., Brossier, R., Métivier, L., Ribodetti, A. & Zhou, W., 2017. An introduction to Full Waveform Inversion, in *Encyclopedia of Exploration Geophysics*, pp. R1–1-R1-40, eds, Grechka, V. & Wapenaar, K., Society of Exploration Geophysics.

Warner, M. *et al.*, 2013. Anisotropic 3D full-waveform inversion, *Geophysics*, **78**(2), R59–R80.
Wu, R.S. & Toksöz, M.N., 1987. Diffraction tomography and multisource holography applied to seismic imaging, *Geophysics*, **52**, 11–25.

Yang, P., Brossier, R., Métivier, L. & Virieux, J., 2016. A review on the systematic formulation of 3D multiparameter full waveform inversion in viscoelastic medium, *Geophys. J. Int.*, **207**(1), 129–149.

Yuan, H., French, S., Cupillard, P. & Romanowicz, B., 2014. Lithospheric expression of geological units in central and eastern north america from

full waveform tomography, *Earth planet. Sci. Lett.*, **402**, 176–186, Special issue on USArray science.

Zhu, H., Bozdağ, E., Peter, D. & Tromp, J., 2012. Structure of the European upper mantle revealed by adjoint tomography, *Nat. Geosci.*, **5**, 493–498.

Zhu, T. & Harris, J.M., 2014. Modeling acoustic wave propagation in heterogeneous attenuating media using decoupled fractional laplacians, *Geophysics*, **79**(3), T105–T116.