Visco-pseudo-elastic TTI FWI/RTM formulation and implementation

Abstract

A method, including: obtaining, with a computer, an initial geophysical model; modeling, with a computer, a forward wavefield based on the initial geophysical model with wave equations including a second order z-derivative in a rotated coordinate system that accounts for a tilted transverse isotropic (TTI) medium; modeling, with a computer, an adjoint wavefield with adjoint wave equations including a second order z-derivative in a rotated coordinate system that accounts for a tilted transverse isotropic (TTI) medium, wherein the wave equations and the adjoint wave equations include relaxation terms accounting for anelasticity of earth in an update of a primary variable and an evolution relationship for the relaxation terms; and obtaining, with a computer, a gradient of a cost function based on a combination of a model of the forward wavefield and a model of the adjoint wavefield.

Claims

1. A method, comprising: obtaining, with a computer, an initial geophysical model; determining a stiffness ratio (a.sub.l); modeling, with a computer, a forward wavefield based on the initial geophysical model with wave equations including a second order z-derivative in a rotated coordinate system that accounts for a tilted transverse isotropic (TTI) medium: modeling, with a computer, an adjoint wavefield with adjoint wave equations including a second order z-derivative in a rotated coordinate system that accounts for a tilted transverse isotropic (TTI) medium, wherein the wave equations and the adjoint wave equations include relaxation terms accounting for anelasticity of earth in an update of a primary variable and an evolution relationship for the relaxation terms, and further wherein the evolution relationship is $\frac{\partial R_l}{\partial t}=a_l\stackrel{.}{\sigma }-{\omega }_l{R}_l,$ in which R.sub.l is a relaxation term, {dot over (σ)} is a first derivative of stress, a.sub.l is the stiffness ratio, and ω.sub.l is a characteristic relaxation frequency, and wherein: the wave equations are $\frac{\partial \stackrel{->}{\sigma }}{\partial t}=\stackrel{.}{\stackrel{->}{\sigma }}+\stackrel{.fwdarw.}{S_{\sigma }}-\underset{l=1}{\overset{\mathrm{NMECH}}{.Math.}}{\omega }_l\stackrel{.fwdarw.}{m_l},\frac{\partial \stackrel{.}{\stackrel{->}{\sigma }}}{\partial t}=C\stackrel{->}{\sigma },\mathrm{and}$ $\frac{\partial \stackrel{.fwdarw.}{m_l}}{\partial t}=-{\omega }_l\stackrel{.fwdarw.}{m_l}+a_l\stackrel{.}{\stackrel{->}{\sigma }},\mathrm{and}$ the adjoint wave equations are $-\frac{\partial \overleftarrow{\sigma }}{\partial t}=\stackrel{.}{\overleftarrow{\sigma }}+\underset{l=1}{\overset{\mathrm{NMECH}}{.Math.}}{a}_l\overleftarrow{m_l},-\frac{\partial \stackrel{.}{\overleftarrow{\sigma }}}{\partial t}=C\overleftarrow{\sigma }+\overleftarrow{S_{\sigma }},\mathrm{and}$ $\frac{\partial \overleftarrow{m_l}}{\partial t}=-{\omega }_l\overleftarrow{m_l}+a_l\stackrel{.}{\overleftarrow{\sigma }},$ with t being time, {dot over (σ)} being a first derivative of stress, a.sub.l being the stiffness ratio, and ω.sub.l being a characteristic relaxation frequency, σ being stress, m.sub.l being a memory response, C being material properties, and S.sub.σ being application of a source to stress; obtaining, with a computer, a gradient of a cost function based on a combination of a model of the forward wavefield and a model of the adjoint wavefield; and updating the initial geophysical model, with the computer, with an adjustment determined from the gradient of the cost function to obtain an updated geophysical model.

2. The method of claim 1, further comprising generating a subsurface image of the updated geophysical model that includes subsurface structures.

3. The method of claim 2, further comprising causing a well to be drilled at a location derived from the subsurface image.

4. The method of claim 1, wherein the primary variable is based on stress.

5. The method of claim 1, further comprising, based at least in part on the updated geophysical model, estimating a subsurface property that indicates hydrocarbon deposits in a subterranean geologic formation.

6. The method of claim 1, wherein the subsurface property is at least one of velocity, density, or impedance.

7. A method, comprising: obtaining, with a computer, an initial geophysical model: determining a stiffness ratio (a.sub.l); modeling, with a computer, a forward wavefield based on the initial geophysical model with wave equations including a second order z-derivative in a rotated coordinate system that accounts for a tilted transverse isotropic (TTI) medium; modeling. with a computer, an adjoint wavefield with adjoint wave equations including a second order z-derivative in a rotated coordinate system that accounts for a tilted transverse isotropic (TTI) medium; wherein the wave equations and the adjoint wave equations include relaxation terms accounting for anelasticity of earth in an update of a primary variable and an evolution relationship for the relaxation terms, and further wherein the evolution relationship is $\frac{\partial R_l}{\partial t}=a_l\stackrel{.}{\sigma }-{\omega }_l{R}_l,$ in which R.sub.l is a relaxation term, {dot over (σ)} is a first derivative of stress, a.sub.l is the stiffness ratio, and ω.sub.l is a characteristic relaxation frequency, and wherein: the wave equations are $\frac{\partial \stackrel{->}{\sigma }}{\partial t}=\stackrel{.}{\stackrel{->}{\sigma }}+\stackrel{.fwdarw.}{S_{\sigma }}-\underset{l=1}{\overset{\mathrm{NMECH}}{.Math.}}{\omega }_l\stackrel{.fwdarw.}{m_l},\frac{\partial \stackrel{.}{\stackrel{->}{\sigma }}}{\partial t}=C\stackrel{->}{\sigma },\mathrm{and}$ $\frac{\partial \stackrel{.fwdarw.}{m_l}}{\partial t}=-{\omega }_l\stackrel{.fwdarw.}{m_l}+a_l\stackrel{.}{\stackrel{->}{\sigma }},\mathrm{and}$ the adjoint wave equations are $-\frac{\partial \overleftarrow{\sigma }}{\partial t}=\stackrel{.}{\overleftarrow{\sigma }}+\underset{l=1}{\overset{\mathrm{NMECH}}{.Math.}}{a}_l\overleftarrow{m_l},-\frac{\partial \stackrel{.}{\overleftarrow{\sigma }}}{\partial t}=C\overleftarrow{\sigma }+\overleftarrow{S_{\sigma }},\mathrm{and}$ $\frac{\partial \overleftarrow{m_l}}{\partial t}=-{\omega }_l\overleftarrow{m_l}+a_l\stackrel{.}{\overleftarrow{\sigma }},$ with t being time, {dot over (σ)} being a first derivative of stress, a.sub.l being the stiffness ratio, and ω.sub.l being a characteristic relaxation frequency, a being stress, m, being a memory response, C being material properties, and S.sub.σ being application of a source to stress; obtaining, with a computer, an imaging condition based on a combination of a model of the forward wavefield and a model of the adjoint wavefield; and generating a subsurface image from the imaging condition.

8. The method of claim 7, wherein the subsurface property is at least one of velocity, density, or impedance.

9. The method of claim 7, wherein the primary variable is based on stress.

10. The method of claim 9, further comprising causing a well to be drilled at a location derived from the subsurface image.

11. The method of claim 7, further comprising, based at least in part on the subsurface image, estimating a subsurface property that indicates hydrocarbon deposits in a subterranean geologic formation.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) While the present disclosure is susceptible to various modifications and alternative forms, specific example embodiments thereof have been shown in the drawings and are herein described in detail. It should be understood, however, that the description herein of specific example embodiments is not intended to limit the disclosure to the particular forms disclosed herein, but on the contrary, this disclosure is to cover all modifications and equivalents as defined by the appended claims. It should also be understood that the drawings are not necessarily to scale, emphasis instead being placed upon clearly illustrating principles of exemplary embodiments of the present invention. Moreover, certain dimensions may be exaggerated to help visually convey such principles.

(2) FIG. 1 illustrates an exemplary method embodying the present technological advancement.

(3) FIG. 2 illustrates an exemplary three-layer VTI true model.

(4) FIG. 3 illustrates a comparison of a far-offset modeled trace.

(5) FIG. 4A illustrates a comparison of stress plots generated via forward simulation with and without attenuation.

(6) FIG. 4B is a comparison of traces of the reflected pressure data generated via forward simulation with and without attenuation.

DETAILED DESCRIPTION

(7) Exemplary embodiments are described herein. However, to the extent that the following description is specific to a particular embodiment, this is intended to be for exemplary purposes only and simply provides a description of the exemplary embodiments. Accordingly, the invention is not limited to the specific embodiments described below, but rather, it includes all alternatives, modifications, and equivalents falling within the true spirit and scope of the appended claims.

(8) A first order implementation is efficient for the models without rotations, such as vertical transversely isotropic (VTI) models. There are several ways to address propagation in the TTI media. For the implementation of the stress-velocity in a TTI media, a rotated staggered or Lebedev grid can be used. It is significantly more expensive than a staggered grid implementation used for the VTI media. It has also been suggested to implement Q related algorithms in frequency domain by pseudo-spectral method, see for example (J. Ramos-Martinez* et al. 2015, Valenciano et al. 2011) or using Low Rank modification of the pseudo-spectral method (Sun et al.).

(9) Modeling of waves in a TTI media is computationally cheaper using a collocated grid (second-order formulation) as compared to a staggered grid (first-order formulation). A time domain second order isotropic system has been discussed in literature before (see, Askan et al. 2007).

(10) An alternative way to look at the second order TTI system has been suggested in (Yi Xie et al. 2015). However, this is different from the present technological advancement. The starting point of the derivation in Yi Xie et al. is the fractional derivative model, wherein the present technological advancement can start with a formulation that uses SLS physical mechanisms. As a result, the final form of the equations in Yi Xie et al. have a dependence on the exponent fractional derivative, while the formulation of the present technological advancement does not.

(11) The general form of a second-order stress formulation for a VTI/TTI system without attenuation can be summarized as

(12) $\begin{array}{cc}\frac{{\partial }^2}{\partial t^2}\left(\begin{array}{c}{\sigma }_H\\ {\sigma }_V\end{array}\right)=A\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2}\right)\left(\begin{array}{c}{\sigma }_H\\ {\sigma }_V\end{array}\right)+B\frac{{\partial }^2}{\partial z^2}\left(\begin{array}{c}{\sigma }_H\\ {\sigma }_V\end{array}\right)+\frac{\partial \stackrel{.fwdarw.}{S_{\sigma }}}{\partial t},& \left(1\right)\end{array}$
where constant density is assumed and a symmetry condition is used, and σ.sub.11≡σ.sub.22. In the TTI media, the z-axis is aligned along an axis of symmetry n, then the governing equation becomes

(13) $\begin{array}{cc}\frac{{\partial }^2}{\partial t^2}\left(\begin{array}{c}{\sigma }_H\\ {\sigma }_V\end{array}\right)=A\left(\frac{{\partial }^2}{\partial x^2}+\frac{{\partial }^2}{\partial y^2}+\frac{{\partial }^2}{\partial z^2}\right)\left(\begin{array}{c}{\sigma }_H\\ {\sigma }_V\end{array}\right)+B\frac{{\partial }^2}{\partial {\tilde{z}}^2}\left(\begin{array}{c}{\sigma }_H\\ {\sigma }_V\end{array}\right)+\frac{\partial \stackrel{.fwdarw.}{S_{\sigma }}}{\partial t},& \left(2\right)\end{array}$
where the current configuration ({tilde over (x)},{tilde over (y)},{tilde over (z)}) is related to the reference configuration (x,y,z) via

(14) 0$\begin{array}{cc}\frac{{\partial }^2}{\partial {\tilde{z}}^2}(.Math.)=[n_x^2\frac{{\partial }^2}{\partial x^2}+n_y^2\frac{{\partial }^2}{\partial x^2}+n_z^2\frac{{\partial }^2}{\partial z^2}+2n_x{n}_y\frac{{\partial }^2}{\partial \mathrm{xy}}+2n_y{n}_z\frac{{\partial }^2}{\partial \mathrm{yz}}+2n_x{n}_z\frac{{\partial }^2}{\partial \mathrm{xz}}]& \left(3\right)\end{array}$

(15) Here, A and B are matrices of material properties.

(16) Since the VTI and TTI systems are related in a straightforward way via a rotation tensor, the following will focus on the VTI system for simplicity. Accordingly, the above governing equations in second order form can be recast into first order form as follows:

(17) $\begin{array}{cc}\frac{\partial \stackrel{.fwdarw.}{\sigma }}{\partial t}=\stackrel{.}{\stackrel{.fwdarw.}{\sigma }}+\stackrel{.fwdarw.}{S_{\sigma }},& \left(4\right)\\ \frac{\partial \stackrel{.}{\stackrel{.fwdarw.}{\sigma }}}{\partial t}=C\stackrel{.fwdarw.}{\sigma },& \left(5\right)\end{array}$
where A and B from (2) are merged into C to shorten the notation, and the corresponding adjoint equations are as follows:

(18) $\begin{array}{cc}-\frac{\partial \overleftarrow{\sigma }}{\partial t}=\stackrel{.}{\overleftarrow{\sigma }}& \left(6\right)\\ -\frac{\partial \stackrel{.}{\overleftarrow{\sigma }}}{\partial t}=C\overleftarrow{\sigma }+\overleftarrow{S_{\sigma }}& \left(7\right)\end{array}$

(19) S.sub.σ is the source term. When the data is recorded by a receiver, the actual signal that the ship (in the context of a marine application) can be applied to either the stress term or the velocity term. Here S.sub.σ means that the seismic source is applied to stress.

(20) With the introduction of attenuation, we can derive the governing forward wave equations as follows:

(21) $\begin{array}{cc}\frac{\partial \stackrel{.fwdarw.}{\sigma }}{\partial t}=\stackrel{.}{\stackrel{.fwdarw.}{\sigma }}+\stackrel{.fwdarw.}{S_{\sigma }}-\underset{l=1}{\overset{\mathrm{NMECH}}{.Math.}}{\omega }_l\stackrel{.fwdarw.}{m_l},& \left(8\right)\\ \frac{\partial \stackrel{.}{\stackrel{.fwdarw.}{\sigma }}}{\partial t}=C\stackrel{.fwdarw.}{\sigma }& \left(9\right)\\ \frac{\partial \stackrel{.fwdarw.}{m_l}}{\partial t}=-{\omega }_l\stackrel{.fwdarw.}{m_l}+a_l\stackrel{.}{\stackrel{.fwdarw.}{\sigma }},& \left(10\right)\end{array}$

(22) and the corresponding adjoint equations are

(23) $\begin{array}{cc}-\frac{\partial \overleftarrow{\sigma }}{\partial t}=\stackrel{.}{\overleftarrow{\sigma }}+\underset{l=1}{\overset{\mathrm{NMECH}}{.Math.}}{a}_l\overleftarrow{m_l}& \left(11\right)\\ -\frac{\partial \stackrel{.}{\overleftarrow{\sigma }}}{\partial t}=C\overleftarrow{\sigma }+\overleftarrow{S_{\sigma }}& \left(12\right)\\ \frac{\partial \overleftarrow{m_l}}{\partial t}=-{\omega }_l\overleftarrow{m_l}+{\alpha }_l\stackrel{.}{\overleftarrow{\sigma }}.& \left(13\right)\end{array}$

(24) With respect to the equations (8)-(13), Q is accounted for in the form of {right arrow over (m.sub.l)} (which are the memory variables). Modeling of viscoelastic behavior can be done in many ways. One such realization is a phenomenological model based on a series and parallel configuration of springs and dashpots, also referred to as the Standard Linear Solid (SLS) approach.

(25) As an example of the second order we summarize strain-related formulation by substituting

(26) $\begin{array}{cc}A=\left(\begin{array}{cc}{c}_{11}& 0\\ \widetilde{c_{13}}& c_{44}\end{array}\right)\mathrm{and}& \left(14\right)\\ B=\left(\begin{array}{cc}{c}_{11}-c_{44}& \widetilde{c_{13}}\\ -\widetilde{c_{13}}& c_{33}-c_{44}\end{array}\right),\mathrm{where}{c}_{13}=\sqrt{\left(v_{\mathrm{pz}}^2-v_{\mathrm{sz}}^2\right)*\left(v_{\mathrm{pn}}^2-v_{\mathrm{sz}}^2\right)}-v_{\mathrm{sz}}^2,c_{44}=v_{\mathrm{sz}}^2,c_{33}=v_{\mathrm{pz}}^2,c_{11}=v_{\mathrm{px}}^2,\mathrm{with}{v}_{\mathrm{px}}=v_{\mathrm{pz}}\sqrt{1+2\mathrm{.Math.}},v_{\mathrm{pn}}=v_{\mathrm{pz}}\sqrt{1+2\delta },\widetilde{c_{13}}=c_{13}+v_{\mathrm{sz}}^2,& \left(15\right)\end{array}$
into (2) for A and B correspondently.

EXAMPLE

(27) FIG. 1 illustrates an exemplary method for a visco-pseudo-elastic TTI FWI using Equations (8) to (13). In step 101, a geophysical model is assumed. A geophysical model gives one or more subsurface properties as a function of location in a region.

(28) In step 103, amplification factors a.sub.l are computed.

(29) In steps 105A and 105B, the forward wavefield model Equations (8)-(10) and the adjoint model Equations (11)-(13) are solved.

(30) In step 107, the gradient of the cost function is obtained from a convolution of the forward and adjoint equations in order to arrive at the gradients of the objective function with respect to the inversion parameter(s).

(31) In step 109, the gradient of the cost function (which provides the rate of the change of the cost function in a given direction) is then used to update the geophysical model in order to minimize the cost function. Step 109 can include searching for an updated geophysical property model that is a perturbation of the initial geophysical property model in the gradient direction that better explains the observed data. The iterative process of FIG. 1 can be repeated until predetermined convergence between measured data and the updated model is obtained.

(32) When the updated assumed model converges, the process proceeds to step 111. In step 111, an updated subsurface model is used to manage hydrocarbons. As used herein, hydrocarbon management includes hydrocarbon extraction, hydrocarbon production, hydrocarbon exploration, identifying potential hydrocarbon resources, identifying well locations, determining well injection and/or extraction rates, identifying reservoir connectivity, acquiring, disposing of and/or abandoning hydrocarbon resources, reviewing prior hydrocarbon management decisions, and any other hydrocarbon-related acts or activities. Based at least in part on the updated geophysical model, step 111 can include estimating a subsurface property that indicates hydrocarbon deposits in a subterranean geologic formation.

(33) The subsurface property can include velocity, density, or impedance. The following describes a modeling experiment to illustrate a comparison of simulated reflected data via a conventional first order formulation and simulated reflected data via the present technological advancement. FIG. 2 illustrates a 3 layer true VTI model (a VTI model is used for this example because there is no suitable and available conventional techniques to compare to the present technological advancement; for the present technological advancement, the title is set to zero in order to be applicable to a VTI model). The model extents in (x,y,z) directions are 2047.5 m×2060 m×2200 m, with grid spacing of Δx=19.5 m, Δy=Δz=20 m. There are two horizontal reflectors located at depths of 660 m and 1540 m. The properties in each layer of the three layers are listed in the order (V.sub.p0,ρ,ε,δ) and listed in the order of increasing depth are (1500.0, 1.0, 0.0, 0.0), (1700.0, 1.0, 0.2, 0.10) and (1500.0, 1.0, 0.2, 0.10). V.sub.s0 is generated using
V.sub.s0=√{square root over (4/3∥ε−δ∥)}V.sub.p0.
A uniform background Q=100 is assumed.

(34) FIG. 3 illustrates a comparison of far-offset modeled traces of the reflected pressure data generated via forward simulation using the true model of FIG. 2 by the conventional first order technique and the present technological advancement. 301 is a trace generated by the present technological advancement and 302 is the difference between the present technological advancement and the conventional first order technique. As expected, the difference is miniscule.

(35) The present technological advancement is also applicable to RTM. The RTM image can be based on the same principles as the first gradient from FWI, with a difference between the two being the imaging condition and is not directly related to Q.

(36) For RTM, a Q model is built, forward modeling and back propagation with the Q model is performed using the relaxation mechanisms as described above for the present technological advancement, and then the RTM imaging condition is performed to form the image.

(37) In another example, a modeling experiment illustrates the comparison of the reflected data simulated via first order formulation with the proposed second order formulation. This example utilized the three layer VTI model of FIG. 2. This VTI model is described above and its characteristics are not repeated here.

(38) Seismic attenuation attenuates the amplitude as well as distorts the phase of the propagating wave as clearly seen from the above example. The proposed second-order formulation accounts for these phase and amplitude corrections in a form consistent with the usual standard first-order formulation.

(39) FIG. 4A illustrates a comparison of stress plots generated via forward simulation with and without attenuation. 401 is the response where the medium is anelastic, 402 is the response without anelasticity, and 403 is the difference between 401 and 402.

(40) FIG. 4B is a comparison of traces of the reflected pressure data generated via forward simulation with and without attenuation. The plot compares a far offset trace simulated by the proposed second-order formulation. 404 is the response where the medium is anelastic, 405 is the response without anelasticity and the difference is 406.

(41) Derivation

(42) The following discussion explains a derivation of the second order visco-acoustic stress-only equations with a desired dispersion relation.

(43) In the absence of viscous effect, the wave propagation can be formulated in the following second order stress-only equation:

(44) $\frac{{\partial }^2\sigma }{\partial t^2}=C\sigma$

(45) For a plane-wave mode σ=exp [i(ωt−q.Math.x)]ξ, where q is the space wave number, ω is temporal frequency, and {right arrow over (ξ)} is the polarization vector, the following dispersion relation holds

(46) $\frac{{\omega }^2}{q^2}=C\left(q\right)$

(47) With N visco-elastic mechanisms (NMECH) to mimic the physical visco-elastic effect, the desired dispersion relation is

(48) ${\omega }^2=C\left(q\right)\left(1-\underset{l=1}{\overset{\mathrm{NMECH}}{.Math.}}\frac{a_l{\omega }_l}{{\omega }_l+i\omega }\right)$

(49) To achieve this goal, the original governing equation is modified to

(50) $\frac{\partial \sigma }{\partial t}=\stackrel{.}{\sigma }-\underset{l=1}{\overset{\mathrm{NMECH}}{.Math.}}{\omega }_l{R}_l$$\frac{\partial \stackrel{.}{\sigma }}{\partial t}=C\sigma$$\frac{\partial R_l}{\partial t}=a_l\stackrel{.}{\sigma }-{\omega }_l{R}_l.$
The first two equations are the governing equations. The third equation shows the evolution of the memory variables R with time.

(51) The normal equation for the above system is

(52) 0$i\omega \hat{\sigma }=\widehat{\stackrel{.}{\sigma }}-\underset{l=1}{\overset{\mathrm{NMECH}}{.Math.}}{\omega }_l{\hat{R}}_l$$i\omega \widehat{\stackrel{.}{\sigma }}=C\left(q\right)\hat{\sigma }$$i\omega {\hat{R}}_l=a_l\widehat{\stackrel{.}{\sigma }}-{\omega }_l{\hat{R}}_l.$

(53) Solving these equations and eliminating {dot over ({circumflex over (σ)})} and {circumflex over (R)}.sub.l yields

(54) $i\omega \hat{\sigma }=\frac{C\left(q\right)}{i\omega }\hat{\sigma }-\underset{l=1}{\overset{\mathrm{NMECH}}{.Math.}}\frac{a_l{\omega }_l}{{\omega }_l+i\omega }\frac{C\left(q\right)}{i\omega }\hat{\sigma },$
which leads to the desired dispersion relation.
Computer Implementation

(55) In all practical applications, the present technological advancement must be used in conjunction with a computer, programmed in accordance with the disclosures herein. Preferably, in order to efficiently perform FWI, the computer is a high performance computer (HPC), known as to those skilled in the art. Such high performance computers typically involve clusters of nodes, each node having multiple CPU's and computer memory that allow parallel computation. The models may be visualized and edited using any interactive visualization programs and associated hardware, such as monitors and projectors. The architecture of system may vary and may be composed of any number of suitable hardware structures capable of executing logical operations and displaying the output according to the present technological advancement. Those of ordinary skill in the art are aware of suitable supercomputers available from Cray or IBM.

(56) Conclusion

(57) The present techniques may be susceptible to various modifications and alternative forms, and the examples discussed above have been shown only by way of example. However, the present techniques are not intended to be limited to the particular examples disclosed herein. Indeed, the present techniques include all alternatives, modifications, and equivalents falling within the spirit and scope of the appended claims.

REFERENCES

(58) The following documents are each incorporated by reference in their entirety: 1. Tarantola, A., (1984) “Inversion of seismic reflection data in the acoustic approximation,” Geophysics, v. 49, no. 8, pp. 1259-1266. 2. Hak, B., and Mulder W. A. (2011) “Seismic attenuation imaging with causality,” Geophysical Journal International, v. 184, no. 1, pp. 439-451. 3. Day, S. M. and Minster J. B. (1984) “Numerical simulation of attenuated wavefields using a Padé approximant method,” Geophysical Journal of International, v. 78, pp. 105-118. 4. Ursin, B. and Toverud, T. (2002) “Comparison of seismic dispersion and attenuation models,” Studia Geophysica et Geodaetica, v. 46, pp. 293-320. 5. Carcione, J. M., Kosloff D. and Kosloff R. (1988) “Viscoacoustic wave propagation simulation in the Earth,” Geophysics, v. 53, no. 6, pp. 769-777. 6. Carcione J. M. (2007) “Wave fields in real media: Wave propagation in anisotropic anelastic, porous and electromagnetic media” Handbook of Geophysical Exploration: Seismic Exploration, v. 38, Elsevier, 2.sup.nd Edition. 7. Muller, T. M., Gurevich, B. and Lebedev, M. (2010) “Seismic wave attenuation and dispersion resulting from wave-induced flow in porous rocks—A review,” Geophysics, v. 75, no. 5, pp. 147-164. 8. Bai, J., Yingst D., Bloor R. and Leveille J. (2012) “Waveform inversion with attenuation,” Society of Exploration Geophysicists Extended Technical Abstract, 5 pgs. 9. Royle, G. T. (2011) “Viscoelastic Orthorhombic Full Wavefield Inversion:

(59) Development of Multiparameter Inversion Methods,” Society of Geophysicist Extended Abstract, pp. 4329-4333. 10. Robertsson, J. O. A., Blanch J. O. and Symes W. W. (1994) “Viscoelastic finite-difference modeling,” Geophysics, v. 59, no. 9, pp. 1444-1456. 11. Hestholm, S., Ketcham, S., Greenfield, R., Moran, M. and McMechan G., (2006) “Quick and accurate Q parameterization in viscoelastic wave modeling,” Geophysics, v. 71, no. 5, pp. 147-150. 12. Charara, M., Barnes, C. and Tarantola, A. (2000) “Full waveform inversion of seismic data for a viscoelastic medium” Methods and Applications of Inversion: Lecture Notes in Earth Sciences, v. 92, pp. 68-81. 13. Käser, M., Dumbser, M., Puente, J. and Igel, H. (2007) “An arbitrary high-order Discontinuous Galerkin method for elastic waves on unstructured meshes—III. Viscoelastic attenuation,” Geophysics Journal International, v. 168, no. 1, pp. 224-242. 14. JafarGandomi, A. and Takenaka, H., (2007) “Efficient FDTD algorithm for plane-wave simulation for vertically heterogeneous attenuative media,” Geophysics, v. 72, no. 4, pp. 43-53. 15. Quintal, B., (2012) “Frequency-dependent attenuation as a potential indicator of oil saturation,” Journal of Applied Geophysics, v. 82, pp. 119-128. 16. Betts, J. T. and Campbell, S. L., (2005) “Discretize then optimize,” Mathematics in Industry: Challenges and Frontiers A Process View: Practice and Theory, Ferguson, D. R. and Peters, T. J., eds., SIAM Publications. 17. Thevenin D. and Janiga G., (2008) “Optimization and computational fluid dynamics,” Springer-Verlag. 18. J. Ramos-Martinez, A. A. Valenciano, Sean Crawley & M. Orlovich (PGS) (2015)

(60) “Viscoacoustic compensation in RTM using the pseudo-analytical extrapolator,” SEG Technical Program Expanded Abstracts, pp. 3954-3958. 19. Valenciano, A. A., Chemingui, N., Whitmore, D and Brandsberg-Dahl, S. (2011) “Wave equation migration with attenuation and anisotropic compensation,” 81st Annual International Meeting, SEG, Expanded Abstracts, pp. 232-236. 20. Aysegul Askan, Volkan Akcelik, Jacobo Bielak, and Omar Ghattas, (2007), “Full Waveform Inversion for Seismic Velocity and Anelastic Losses in Heterogeneous Structures,” Bulletin of the Seismological Society of America, v. 97, no. 6, pp. 1990-2008, December 2007, doi: 10.1785/0120070079. 21. Junzhe Sun, Sergey Fomel, and Tieyuan Zhu, (2014), “Viscoacoustic modeling and imaging using low-rank approximation,” SEG Technical Program Expanded Abstracts, v. 80, no. 5, pp. A-103-A108. 22. Yi Xie, James Sun, Yu Zhang, Joe Zhou (2015) “Compensating for visco-acoustic effects in TTI reverse time migration,” SEG New Orleans Annual Meeting, pp. 3996-3998. 23. Same as #26 US2016/291,178 (below) 24. Denli, H., Akcelik, V., Kanevsky, A., Trenev, D., White, L. and Lacasse, M. (2013) “Full-Wavefield Inversion for Acoustic Wave Velocity and Attenuation,” SEG Annual Meeting, Society of Exploration Geophysicists, pp. 980-985. 25. U.S. Pat. No. 8,194,498 (2012) Du, X. et al., 22 pgs. 26. U.S. Patent Publication 2016/291,178 (2016) Xie, Y. et al., 16 pgs. 27. U.S. Patent Publication 2015/362,622 (2015) Denli, H. et al., 20 pgs.