**6. Conclusions**

In this study, we have proposed an elastic full-waveform inversion method based on a modified acoustic-elastic coupled equation. This method uses the modified AEC equation to simultaneously compute subsurface pressure and displacement wavefields. It adopts a weighted misfit function to quantify the contributions of pressure and displacement records. With the adjoint operator, it can simultaneously make use of OBS 4C data to reconstruct multiple elastic parameters. The preconditioned TGN algorithm with a multiparameter diagonal Hessian operator is developed to cope with unbalanced illumination and coupling effects. The sensitivity analysis and numerical experiments have validated that the AEC-EFWI method can yield a higher spatial resolution, provide more accurate elastic parameter inversions, and have a higher convergence rate. The discussion reveals the potential of this AEC-EFWI method in inverting elastic parameters using 1C pressure data for OBS and marine towed-streamer cases.

**Author Contributions:** Conceptualization, M.S.; Methodology, M.S.; Software, M.S.; Formal analysis, M.S.; Validation, S.J.; Writing, M.S. and S.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the National Key Research and Development Program of China Project (Grant No. 2018YFC0603502), the National Natural Science Foundation of China (Grant No. 41930105), and the National Key R&D Program of China (Grant No. 2018YFC0310100).

**Acknowledgments:** We thank Pengfei Yu from Hohai University for providing the code for the original acoustic-elastic coupled equation.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **Appendix A. Derivation of the Modified Acoustic-Elastic Coupled Equation**

We start from the standard 2D time-domain displacement–stress elastic wave equation [47], given by

$$\begin{cases} \rho \frac{\partial^2 u\_x}{\partial t^2} = \frac{\partial \tau\_{xx}}{\partial x} + \frac{\partial \tau\_{xz}}{\partial z} \\ \rho \frac{\partial^2 u\_z}{\partial t^2} = \frac{\partial \tau\_{xz}}{\partial x} + \frac{\partial \tau\_{zz}}{\partial z} \\ \qquad \tau\_{xx} = (\lambda + 2\mu) \frac{\partial u\_x}{\partial x} + \lambda \frac{\partial u\_z}{\partial z} + f, \\ \tau\_{xz} = (\lambda + 2\mu) \frac{\partial u\_z}{\partial z} + \lambda \frac{\partial u\_x}{\partial x} + f \\ \tau\_{xz} = \mu \left(\frac{\partial u\_x}{\partial z} + \frac{\partial u\_z}{\partial x}\right) \end{cases} \tag{A1}$$

where *ux* and *uz* denote the horizontal- and vertical-particle displacement components, respectively; *τxx* and *τzz* are the normal stress components; and *τxz* is the shear stress component. *λ* and *μ* are the Lamé constants, and *ρ* is density. *f* indicates the source function that is implemented in *τxx* and *τzz* to generate a pure *P*-wave.

In tensor analysis, the stress tensor **T** can be decomposed into the isotropic pressure −*p***I** and deviatoric *τ<sup>s</sup>* parts,

$$\mathbf{T} = \mathbf{\tau}^s - p\mathbf{I}.\tag{A2}$$

For 2D cases, the pressure wavefield satisfies

$$\begin{split} p &= -\frac{1}{2}tr(\mathbf{T}) \\ &= -\frac{1}{2}(\mathfrak{r}\_{xx} + \mathfrak{r}\_{zz}) = -\left(\lambda + \mu\right)\left(\frac{\partial u\_x}{\partial x} + \frac{\partial u\_z}{\partial z}\right), \end{split} \tag{A3}$$

and the deviatoric stress components become

$$\begin{cases} \tau\_{xx}^s = \tau\_{xx} + p = \mu \frac{\partial u\_x}{\partial x} - \mu \frac{\partial u\_z}{\partial z} \\ \tau\_{zz}^s = \tau\_{zz} + p = \mu \frac{\partial u\_z}{\partial z} - \mu \frac{\partial u\_x}{\partial x} \\ \tau\_{xz}^s = \tau\_{xz} = \mu \frac{\partial u\_x}{\partial z} + \mu \frac{\partial u\_z}{\partial x} \end{cases} \tag{A4}$$

These equations formulate the original acoustic-elastic coupled equation; we refer to the work in [15].

Note that, it is redundant to simultaneously compute the pressure and two normal stress components. Thus, we redefine

$$
\tau\_n^s = \tau\_{xx}^s = -\tau\_{zz}^s, \quad \tau\_s^s = \tau\_{xz}^s. \tag{A5}
$$

where *τsn* and *τss* are the redefined normal and deviatoric stress components. The simplified equation becomes

$$\begin{cases} \rho \frac{\partial^2 u\_x}{\partial t^2} = \frac{\partial \left(\tau\_n^s - p\right)}{\partial x} + \frac{\partial \tau\_s^s}{\partial z} \\ \rho \frac{\partial^2 u\_z}{\partial t^2} = \frac{\partial \tau\_s^s}{\partial x} + \frac{\partial \left(-\tau\_n^s - p\right)}{\partial z} = 0 \\ \qquad p = -\left(\lambda + \mu\right) \left(\frac{\partial u\_x}{\partial x} + \frac{\partial u\_z}{\partial z}\right). \end{cases} \tag{A6}$$

$$\tau\_n^s = \mu \left(\frac{\partial u\_x}{\partial x} - \frac{\partial u\_z}{\partial z}\right)$$

$$\tau\_s^s = \mu \left(\frac{\partial u\_x}{\partial z} + \frac{\partial u\_z}{\partial x}\right)$$

Compared with the original one, this modified equation can provide subsurface 4C elastic wavefields with same accuracy but less computing resources.

For completeness purpose, we also derive the modified AEC equation in 3D. The 3D original AEC equation is given by

$$\begin{cases} \rho \frac{\partial^2 u\_x}{\partial t^2} = \frac{\partial \left(\tau\_{xx}^s - p\right)}{\partial x} + \frac{\partial \tau\_{xy}^s}{\partial y} + \frac{\partial \tau\_{xz}^s}{\partial z} \\ \rho \frac{\partial^2 u\_y}{\partial t^2} = \frac{\partial \tau\_{xy}^s}{\partial x} + \frac{\partial \left(\tau\_{yy}^s - p\right)}{\partial y} + \frac{\partial \tau\_{yz}^s}{\partial z} \\ \rho \frac{\partial^2 u\_z}{\partial t^2} = \frac{\partial \tau\_{xz}^s}{\partial x} + \frac{\partial \tau\_{yz}^s}{\partial y} + \frac{\partial \left(\tau\_{zz}^s - p\right)}{\partial z} \end{cases} \tag{A7}$$

and

$$\begin{cases} \begin{aligned} p &= -\left(\lambda + \frac{2}{3}\mu\right) \left(\frac{\partial u\_x}{\partial x} + \frac{\partial u\_y}{\partial y} + \frac{\partial u\_z}{\partial z}\right) \\ \tau\_{xx}^s &= \frac{2}{3}\mu \left(2\frac{\partial u\_x}{\partial x} - \frac{\partial u\_y}{\partial y} - \frac{\partial u\_z}{\partial z}\right) \\ \tau\_{yy}^s &= \frac{2}{3}\mu \left(2\frac{\partial u\_y}{\partial y} - \frac{\partial u\_x}{\partial x} - \frac{\partial u\_z}{\partial z}\right) \\ \tau\_{zz}^s &= \frac{2}{3}\mu \left(2\frac{\partial u\_z}{\partial z} - \frac{\partial u\_x}{\partial x} - \frac{\partial u\_y}{\partial y}\right) \\ \tau\_{xz}^s &= \mu \left(\frac{\partial u\_x}{\partial z} + \frac{\partial u\_z}{\partial x}\right) \\ \tau\_{xy}^s &= \mu \left(\frac{\partial u\_x}{\partial y} + \frac{\partial u\_y}{\partial x}\right) \\ \tau\_{yz}^s &= \mu \left(\frac{\partial u\_y}{\partial z} + \frac{\partial u\_z}{\partial y}\right) \end{aligned} \tag{A8}$$

Because *τsxx* + *τsyy* + *τszz* = 0, we can similarly provide a 3D modified AEC equation by replacing *τsyy* as −(*τsxx* + *τszz*) to reduce the number of formulas.

#### **Appendix B. Derivation of Gradient Computation for Aec-Efwi Method**

The modified AEC equation can be rewritten as

$$\mathbf{Sw} = \mathbf{f} \tag{A9}$$

where **w** = (*ux*, *uz*, *p*, *<sup>τ</sup>sn*, *τss* )*T* denotes the elastic wavefields, **f** = 0, 0, *fp*, 0, 0*<sup>T</sup>* denotes the source wavelet, and **S** denotes the parameter derivative matrix,

$$\mathbf{S} = \begin{pmatrix} \rho \frac{\partial^2}{\partial t^2} & 0 & \frac{\partial}{\partial x} & -\frac{\partial}{\partial x} & -\frac{\partial}{\partial z} \\ 0 & \rho \frac{\partial^2}{\partial t^2} & \frac{\partial}{\partial z} & \frac{\partial}{\partial z} & -\frac{\partial}{\partial x} \\ (\lambda + \mu) \frac{\partial}{\partial x} & (\lambda + \mu) \frac{\partial}{\partial z} & 1 & 0 & 0 \\ -\mu \frac{\partial}{\partial x} & \mu \frac{\partial}{\partial z} & 0 & 1 & 0 \\ -\mu \frac{\partial}{\partial z} & -\mu \frac{\partial}{\partial x} & 0 & 0 & 1 \end{pmatrix} . \tag{A10}$$

Substituting Equation (2) into Equation (3), the gradient satisfies

$$\frac{\partial E}{\partial \mathbf{m}} = \left[ -\left( \mathbf{S}^{-1} \right) \frac{\partial \mathbf{S}}{\partial \mathbf{m}} \mathbf{w} \right]^T \mathbf{f}'.\tag{A11}$$

With consideration of the self-adjoint assumption **<sup>S</sup>**−<sup>1</sup>*<sup>T</sup>* = **S***<sup>T</sup>*−1, it can be rewritten by

$$\frac{\partial E}{\partial \mathbf{m}} = -\left[\frac{\partial \mathbf{S}}{\partial \mathbf{m}} \mathbf{w}\right]^T \left(\mathbf{S}^T\right)^{-1} \mathbf{f}'.\tag{A12}$$

Here, **S***<sup>T</sup>* is the adjoint operator of the modified AEC-equation, given by

$$\mathbf{S}^{T} = \begin{pmatrix} \rho \frac{\partial^2}{\partial t^2} & 0 & -\frac{\partial}{\partial x} \left(\lambda + \mu\right) & \frac{\partial}{\partial x} \mu & \frac{\partial}{\partial z} \mu \\ 0 & \rho \frac{\partial^2}{\partial t^2} & -\frac{\partial}{\partial z} \left(\lambda + \mu\right) & -\frac{\partial}{\partial x} \mu & \frac{\partial}{\partial x} \mu \\ -\frac{\partial}{\partial x} & -\frac{\partial}{\partial z} & 1 & 0 & 0 \\ \frac{\partial}{\partial x} & -\frac{\partial}{\partial z} & 0 & 1 & 0 \\ \frac{\partial}{\partial z} & \frac{\partial}{\partial x} & 0 & 0 & 1 \end{pmatrix} . \tag{A13}$$

According to the adjoint-state method, we define **wˆ** = **<sup>S</sup>**−<sup>1</sup>*<sup>T</sup>***f** as the solution of the adjoint equation,

$$\mathbf{S}^T \mathbf{\hat{w}} = \mathbf{f}' \tag{A14}$$

where **wˆ** = (*u*<sup>ˆ</sup>*x*, *<sup>u</sup>*<sup>ˆ</sup>*z*, *p*ˆ, *<sup>τ</sup>*<sup>ˆ</sup>*sn*, *τ*ˆ*ss*)*T* is the adjoint variables as used in Equation (8).

For the parameterization of Lamé constants and density **m** = [*λ*, *μ*, *<sup>ρ</sup>*]*<sup>T</sup>*, the partial derivative matrices can be given by

$$\frac{\partial \mathbf{S}}{\partial \lambda} = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial z} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \tag{A15}$$

$$\frac{\partial \mathbf{S}}{\partial \mu} = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial z} & 0 & 0 & 0 \\ -\frac{\partial}{\partial x} & \frac{\partial}{\partial z} & 0 & 0 & 0 \\ -\frac{\partial}{\partial z} & -\frac{\partial}{\partial x} & 0 & 0 & 0 \end{pmatrix}, \tag{A16}$$

and

$$
\frac{\partial \mathbf{S}}{\partial \rho} = \begin{pmatrix}
\frac{\partial^2}{\partial t^2} & 0 & 0 & 0 & 0 \\
0 & \frac{\partial^2}{\partial t^2} & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0
\end{pmatrix} .\tag{A17}
$$

Substituting Equations (A14–A17) into Equation (A12), the gradients of the Lamé constants and density can be given by

$$\begin{split} \frac{\partial E}{\partial \lambda} &= -\frac{p}{\lambda + \mu} \not{p} \\ \frac{\partial E}{\partial \mu} &= -\frac{p}{\lambda + \mu} \not{p} - \frac{\tau\_n^s}{\mu} \mathfrak{k}\_n^s - \frac{\tau\_s^s}{\mu} \mathfrak{k}\_s^s \\ \frac{\partial E}{\partial \rho} &= \frac{\partial^2 u\_x}{\partial t^2} \not{u}\_x + \frac{\partial^2 u\_z}{\partial t^2} \not{u}\_z. \end{split} \tag{A18}$$
