Seismic modeling

Abstract

A method of seismic modeling using an elastic model, the elastic model including a grid having a grid spacing sized such that, when synthetic seismic data is generated using the elastic model, synthetic shear wave data exhibits numerical dispersion, the method including: generating generated synthetic seismic data using the elastic model, wherein the generated synthetic seismic data includes synthetic compression wave data and synthetic shear wave data, and modifying the generated synthetic seismic data to produce modified synthetic seismic data by attenuating at least some of the synthetic shear wave data in order to attenuate at least some of the numerically dispersive data.

Claims

1. A method of seismic processing comprising: performing seismic modeling using an elastic model of a geographical structure to provide modified synthetic seismic data, the elastic model comprising a grid having a grid spacing sized such that, when synthetic seismic data is generated using the elastic model, synthetic shear wave data exhibits numerical dispersion, the seismic modeling comprising: generating generated synthetic seismic data using said elastic model, wherein the generated synthetic seismic data comprises synthetic compression wave data and synthetic shear wave data; and modifying the generated synthetic seismic data to produce the modified synthetic seismic data by attenuating at least some of the synthetic shear wave data in order to attenuate at least some of the numerically dispersive data; and the method of seismic processing further comprising using measured seismic data and the modified synthetic seismic data to produce an updated model of the geographical structure; wherein using the seismic model to prospect for hydrocarbons comprises using the seismic model to identify locations for drilling and/or well locations; and further comprising drilling at and/or into the identified locations.

2. A method as claimed in claim 1, wherein the grid spacing is sized such that, when synthetic seismic data is generated using the elastic model, the synthetic shear wave data exhibits numerical dispersion and substantially all of the synthetic compression wave data exhibits no numerical dispersion.

3. A method as claimed in claim 1, wherein the grid spacing is sized to be as large as possible without substantially any of the synthetic compression wave data comprising numerical dispersion.

4. A method as claimed in claim 1, wherein the grid spacing is sized such that the slowest compression wave data is the slowest wave that undergoes substantially no numerical dispersion in the model.

5. A method as claimed in claim 1, wherein the model comprises a time step, wherein the time step is sized such that, when generated synthetic seismic data is generated using the elastic model, the generated synthetic seismic data is substantially stable.

6. A method as claimed in claim 1, wherein the generated synthetic shear wave data comprises a first portion consisting of data that exhibits numerical dispersion and a second portion consisting of data that exhibits no numerical dispersion, and the method comprises attenuating at least some of the first portion.

7. A method as claimed in claim 1, wherein attenuating at least some of the generated synthetic shear wave data comprises applying one or more filters to the generated synthetic seismic data.

8. A method as claimed in claim 1, wherein the elastic model is a 3D elastic model.

9. A method as claimed in claim 1, wherein the elastic model is an isotropic elastic model.

10. A method as claimed in claim 1, wherein the elastic model is an anisotropic elastic model.

11. A method as claimed in claim 1, wherein the model comprises at least one location where at least one reflection in the generated synthetic seismic data occurs, and the method comprises obtaining one or more reflection coefficients for the reflection location and using said one or more reflection coefficients when generating the generated synthetic seismic data.

12. A method as claimed in claim 1, wherein the elastic model is an elastic model representing a real geophysical structure, and wherein the measured seismic data is acquired from the real geophysical structure, the method comprising: comparing the measured seismic data with the modified seismic data; and updating the model based upon the comparison to produce a more accurate model of the geophysical structure.

13. A method as claimed in claim 12, wherein the model that is updated based upon the comparison is a reflectivity model and/or the elastic model on which the generated geophysical data is generated.

14. A method as claimed in claim 12, comprising repeating the steps of generating generated synthetic seismic data, attenuating the generated synthetic seismic data to produce modified synthetic seismic data, comparing the modified synthetic seismic data and updating the model based on the comparison.

15. A method as claimed in claim 12, comprising acquiring the measured seismic data from the geophysical structure.

16. A computer program product comprising computer readable instructions that, when run on a computer, is configured to perform the method of claim 1.

17. A method of prospecting for hydrocarbons, comprising: performing the method of claim 1; and using the seismic model to prospect for hydrocarbons.

18. A method of producing hydrocarbons, comprising: performing the method of claim 1; and producing hydrocarbons.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Certain preferred embodiments will now be described by way of example only and with reference to the accompanying drawings, in which

(2) FIG. 1 shows exemplary data produced using a model according to an embodiment of the present invention.

DETAILED DESCRIPTION

(3) One embodiment of the present invention may involve performing 3D elastic FD simulations of seismic data using a 3D elastic model. These simulations can then be used for further seismic processing. As mentioned below, this model may be isotropic or anisotropic. As should be clear from the above description, the present invention could be performed using any elastic model, e.g. 1D, 2D or 3D with any time-marching algorithm.

(4) The inventors have observed that it is computationally challenging to perform elastic simulations when the shear wave velocities are low, particularly 3D elastic simulations, particularly 3D elastic anisotropic simulations, as the computational requirements become very high, as has been discussed in detail above. The inventors solution to this challenge is to approximate the elastic simulation not by using an equivalent acoustic model (such as done by Chapman et al. (2010, 2014)) but rather by defining the grid based on the lowest acoustic (e.g. compression wave) wave velocity. Basing the grid on the lowest acoustic wave velocity means selecting a grid size that is only just small enough to ensure the acoustic waves display no numerical dispersion. Compared to standard elastic FD modeling, the number of grid points in the model using this approach can be significantly reduced, since in standard elastic modeling the grid is based on the slowest shear wave velocity which is slower than the slowest acoustic wave velocity. Increasing the grid size in this way significantly reduces the memory and compute time required (i.e. the computational efficiency is increased).

(5) However, using such a large grid spacing is not done in the prior art since the modeled shear waves can become dispersive, which is an unwanted effect. In order to address this problem, the inventors have developed a method where filters are applied to the modeled wavefields that attenuates or eliminates the modeled dispersive shear wavefield. Thus, some of the generated modeled wavefield is attenuated or eliminated. In the prior art, there is no attenuation or removal of modeled wavefields. This step may have the disadvantage of actually reducing the amount of data (e.g. the shear wave data may be attenuated or eliminated). However, it does mean that the compression wave data (that is preserved, i.e. not attenuated) is elastic compression wave data and hence is more accurate than, for instance, acoustic modeled compression wave data. Thus, amongst other effects, the invention allows for accurate elastic modeling of compression wave data whilst avoiding the large computational effort.

(6) Regarding a mathematical description of the method, at each location of the grid v.sub.i (i=1, 2, 3) represents the three components of the modeled particle velocity field and σ.sub.ij (i, j=1, 2, 3) represents the components of the modeled stress field. In an embodiment of the invention, spatial filters are designed and are applied to the modeled data, before summing the filtered data, to attenuate shear wave energy while preserving compressional wave energy. These filters depend on the local medium parameters where they are applied to the modeled data as well as frequency.

(7) At each grid location, d.sub.i represents the modeled synthetic data (e.g. particle velocities and stresses, which may total 12 components so that i=1, . . . , 12). Generally, at every location in the grid, say (x,y,z), local filters F.sub.ij are defined (or designed or selected) that are applied to the data
D.sub.i=F.sub.ijd.sub.j  (2)
The filters are chosen such that the new modified data D.sub.i contain little or no dispersive shear waves. The data D.sub.i can be used in further seismic processing, such as FWI and/or RTM and/or AVO analysis.

(8) The filters F.sub.ij can contain (only) two parts. The first part removes the high wavenumber energy in the wavefields, while the second part separates the compression (P-wave) and shear (S-wave) components.

(9) In some embodiments, the elastic model used is an isotropic model. An elastic isotropic model can be described by density ρ, compressional velocity α and shear velocity β, where velocities can depend on frequency ω. The wavefield data d.sub.i that governs the wave propagation is
d.sub.i=(d.sub.1,d.sub.2,d.sub.3,d.sub.4,d.sub.5,d.sub.6)=(v.sub.1,v.sub.2,v.sub.3,σ.sub.13,σ.sub.23,σ.sub.33).

(10) In an isotropic model, the 12 components of the data d.sub.i mentioned above may be reduced to six components due to symmetries and couplings in the model.

(11) The filtered data
D.sub.i=(D.sub.1,D.sub.2,D.sub.3,D.sub.4,D.sub.5,D.sub.6)=({tilde over (v)}.sub.1,{tilde over (v)}.sub.2,{tilde over (v)}.sub.3,{tilde over (σ)}.sub.13,{tilde over (σ)}.sub.23,{tilde over (σ)}.sub.33)
is obtained from equation two. An example of the specific filters F.sub.ij (of which there are 36 components in this example, since i,j=0, . . . , 6) that can be used in this isotropic example are given below.
D.sub.1={tilde over (v)}.sub.1=F.sub.11v.sub.1+F.sub.12v.sub.2++F.sub.13v.sub.3+F.sub.14σ.sub.13+F.sub.15σ.sub.23+F.sub.16σ.sub.33
with filters

(12) $F_{11}=-2{\beta }^2{\omega }^{-2}\frac{d^2}{{\mathrm{dx}}_1^2};F_{12}=-2{\beta }^2{\omega }^{-2}\frac{d^2}{{\mathrm{dx}}_1{\mathrm{dx}}_2};F_{16}=\frac{i}{\mathrm{\rho \omega }}\frac{d}{{\mathrm{dx}}_1}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}$$F_{13}=F_{14}=F_{15}=0$
and
D.sub.2={tilde over (v)}.sub.2=F.sub.21v.sub.1+F.sub.22v.sub.2+F.sub.23v.sub.3+F.sub.24σ.sub.13+F.sub.25σ.sub.23+F.sub.26σ.sub.33
with filters

(13) $F_{21}=F_{12};F_{22}=-2{\beta }^2{\omega }^{-2}\frac{d^2}{{\mathrm{dx}}_2^2};F_{26}=\frac{i}{\mathrm{\rho \omega }}\frac{d}{{\mathrm{dx}}_2}$$F_{23}=F_{24}=F_{25}=0$
and
D.sub.3={tilde over (v)}.sub.3=F.sub.31v.sub.1+F.sub.32v.sub.2+F.sub.33v.sub.3+F.sub.34σ.sub.13+F.sub.35σ.sub.23+F.sub.36σ.sub.33
with filters

(14) $F_{33}=1+2{\beta }^2{\omega }^{-2}\left(\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right);F_{34}=F_{16};F_{35}=F_{26}$$F_{31}=F_{32}=F_{36}=0$
and
D.sub.4={tilde over (σ)}.sub.13=F.sub.41v.sub.1+F.sub.42v.sub.2+F.sub.43v.sub.3+F.sub.44σ.sub.13+F.sub.45σ.sub.23+F.sub.46σ.sub.33
with filters

(15) $F_{43}=\frac{2i{\mathrm{\rho \beta }}^2}{\omega }\frac{d}{{\mathrm{dx}}_1}\left(1+2{\beta }^2{\omega }^{-2}\left[\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right]\right);F_{44}=F_{11};F_{45}=F_{12}$$F_{41}=F_{42}=F_{46}=0$
and
D.sub.5={tilde over (σ)}.sub.23=F.sub.51v.sub.1+F.sub.52v.sub.2+F.sub.53v.sub.3+F.sub.54σ.sub.13+F.sub.55σ.sub.23+F.sub.56σ.sub.33
with filters

(16) $F_{53}=\frac{2i{\mathrm{\rho \beta }}^2}{\omega }\frac{d}{{\mathrm{dx}}_2}\left(1+2{\beta }^2{\omega }^{-2}\left[\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right]\right);F_{54}=F_{12};F_{55}=F_{22}$$F_{51}=F_{52}=F_{56}=0$
and
D.sub.6={tilde over (σ)}.sub.33=F.sub.61v.sub.1+F.sub.62v.sub.2+F.sub.63v.sub.3+F.sub.64σ.sub.13+F.sub.65σ.sub.23+F.sub.66σ.sub.33
with filters
F.sub.61=F.sub.43;F.sub.62=F.sub.53;F.sub.66=F.sub.33;F.sub.63=0;F.sub.64=0;F.sub.65=0.

(17) When these filters are applied to the modeled data d.sub.i at each location, the dispersive shear wave data is attenuated at each location to produce the modified data D.sub.i. The modified data is summed and then can be used for further seismic processing. The modified data D.sub.i, having been filtered, comprises non-dispersive stable compression wave data. In other embodiments, the elastic model used is an anisotropic model, such as a model with transverse isotropy. In geophysics, a common assumption is that the rock formations of the subsurface are locally polar anisotropic or transversely isotropic. Such models can be described by density ρ and five independent stiffness parameters C.sub.12, C.sub.13, C.sub.33, C.sub.44, C.sub.66, or by compressional velocity α, shear velocity β, and Thomsen parameters ε, δ, γ. The Thompson parameters are related to the stiffnesses as

(18) $\mathrm{.Math.}=\frac{C_{11}-C_{33}}{2C_{33}}$$\delta =\frac{{\left(C_{13}+C_{44}\right)}^2-{\left(C_{33}-C_{44}\right)}^2}{2C_{33}\left(C_{33}-C_{44}\right)}$$\gamma =\frac{C_{66}-C_{44}}{2C_{44}}$

(19) The wavefield data d.sub.i that governs the wave propagation is
d.sub.i=(d.sub.1,d.sub.2,d.sub.3,d.sub.4,d.sub.5,d.sub.6)=(v.sub.1,v.sub.2,v.sub.3,σ.sub.13,σ.sub.23,σ.sub.33).
Again, in a transverse isotropic model, the 12 components of the data d.sub.i mentioned above may be reduced to six components due to symmetries and couplings in the model.

(20) The filtered data
D.sub.i=(D.sub.1,D.sub.2,D.sub.3,D.sub.4,D.sub.5,D.sub.6)=({tilde over (v)}.sub.1,{tilde over (v)}.sub.2,{tilde over (v)}.sub.3,{tilde over (σ)}.sub.13,{tilde over (σ)}.sub.23,{tilde over (σ)}.sub.33).
is obtained from equation two. An example of the specific filters F.sub.ij that can be used in this anisoptropic (transverse isoptric) model such that the filtered summed data does not contain significant dispersive shear wave energy are given below.

(21) To simplify notation, we write

(22) $\left[\begin{array}{c}{\tilde{v}}_3\\ {\tilde{\sigma }}_{13}\\ {\tilde{\sigma }}_{23}\end{array}\right]=E\left[\begin{array}{c}{F}_{33}\\ F_{43}\\ F_{53}\end{array}\begin{array}{c}{F}_{34}\\ F_{44}\\ F_{54}\end{array}\begin{array}{c}{F}_{35}\\ F_{45}\\ F_{55}\end{array}\right]\left[\begin{array}{c}{v}_3\\ {\sigma }_{13}\\ {\sigma }_{23}\end{array}\right]$$\mathrm{and}\left[\begin{array}{c}{\tilde{\sigma }}_{33}\\ {\tilde{v}}_1\\ {\tilde{v}}_2\end{array}\right]=E\left[\begin{array}{c}{F}_{66}\\ F_{16}\\ F_{26}\end{array}\begin{array}{c}{F}_{61}\\ F_{11}\\ F_{21}\end{array}\begin{array}{c}{F}_{62}\\ F_{12}\\ F_{22}\end{array}\right]\left[\begin{array}{c}{v}_3\\ {\sigma }_{13}\\ {\sigma }_{23}\end{array}\right]$$\mathrm{where}$$E=\frac{1}{\left(-C_{13}{\omega }^{-2}\left[\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right]+C_{33}\mathrm{qQ}\right)Q-C_{44}\left(q+Q\right){\omega }^{-2}\left[\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right]}$$\mathrm{and}$$Q=\frac{\rho +C_{11}{\omega }^{-2}\left[\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right]-C_{44}{q}^2}{\left(C_{13}+C_{44}\right)q}$$q=\frac{1}{\sqrt{2}}\sqrt{K_1+\sqrt{K_1^2-4K_1}}$$K_1=\frac{p\left(C_{33}+C_{44}\right)-\left(C_{13}^2-C_{11}{C}_{33}+2C_{13}{C}_{44}\right){\omega }^{-2}\left[\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right]}{C_{33}{C}_{44}}$$K_2=-\frac{C_{11}{\omega }^{-2}\left[\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right]+\rho }{C_{33}}$$K_3=-{\omega }^{-2}\left[\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right]-\frac{\rho }{C_{44}}$$C_{11}=C_{12}+2C_{66}$

(23) Now the local filters can be expressed in terms of differential operators and pseudo-differential operators in the following way

(24) $F_{33}=\left(-C_{13}{\omega }^{-2}\left[\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right]+C_{33}\mathrm{qQ}\right)Q$$F_{34}=\frac{\mathrm{iQ}}{\omega }\frac{d}{{\mathrm{dx}}_1}$$F_{35}=\frac{\mathrm{iQ}}{\omega }\frac{d}{{\mathrm{dx}}_2}$$F_{43}=C_{44}\left(q+Q\right)\left(-C_{13}{\omega }^{-2}\left[\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right]+C_{33}\mathrm{qQ}\right)\frac{i}{\omega }\frac{d}{{\mathrm{dx}}_1}$$F_{44}=-C_{44}\left(q+Q\right){\omega }^{-2}\frac{d^2}{{\mathrm{dx}}_1^2}$$F_{45}=-C_{44}\left(q+Q\right){\omega }^{-2}\frac{d^2}{{\mathrm{dx}}_1{\mathrm{dx}}_2}$$F_{53}=C_{44}\left(q+Q\right)\left(-C_{13}{\omega }^{-2}\left[\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right]+C_{33}\mathrm{qQ}\right)\frac{i}{\omega }\frac{d}{{\mathrm{dx}}_2}$$F_{54}=-C_{44}\left(q+Q\right){\omega }^{-2}\frac{d^2}{{\mathrm{dx}}_1{\mathrm{dx}}_2}$$F_{55}=-C_{44}\left(q+Q\right){\omega }^{-2}\frac{d^2}{{\mathrm{dx}}_2^2}$

(25) Finally,

(26) 0$\left[\begin{array}{ccc}{F}_{66}& F_{61}& F_{62}\\ F_{16}& F_{11}& F_{12}\\ F_{26}& F_{21}& F_{22}\end{array}\right]={\left[\begin{array}{ccc}{F}_{33}& F_{34}& F_{35}\\ F_{43}& F_{44}& F_{45}\\ F_{53}& F_{54}& F_{55}\end{array}\right]}^T$
where T denotes transpose. The other terms of F.sub.ij, such as F.sub.13, F.sub.14, F.sub.15, F.sub.23, F.sub.24, F.sub.25, F.sub.31, F.sub.32, F.sub.35, etc. may (or may not) be 0, for instance due to symmetries.

(27) The local filters given above are exact but cumbersome expressions to work with because the differential operators may not be evaluated while under the square root sign. In addition, the filters as they stand are non-local since differential operators appear in the denominators. To fix this, a solution is to expand the filters in terms of Taylor series. As series expansions can be developed in many ways and to varying accuracy, we give one example to demonstrate how the filters F.sub.33, F.sub.34, F.sub.35 may appear in terms of Thomsen parameters after a series expansion:

(28) $F_{33}=1+\frac{2{\beta }^2}{{{\omega }^2\left({\alpha }^2-{\beta }^2\right)}^2}\left[\left(1+\delta \right){\alpha }^2-{\beta }^2\right]\left[\left(1+\frac{\delta }{2}\right){\alpha }^2-{\beta }^2\right]\left(\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right)-\frac{2{\alpha }^2{\beta }^2}{{\omega }^4\left({\alpha }^2-{\beta }^2\right)}\left[3{\alpha }^2\left(\mathrm{.Math.}-\delta \right)-2{\beta }^2\delta \right]{\left(\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right)}^2{F}_{34}=1+\frac{{\alpha }^2}{{\alpha }^2-{\beta }^2}\delta -\frac{2i{\alpha }^2}{{{\omega }^3\left({\alpha }^2-{\beta }^2\right)}^3}\left[\left(1+\delta \right){\alpha }^2-{\beta }^2\right]\hspace{1em}\hspace{1em}\left[\left({\alpha }^2+{\beta }^2\right)\left(\mathrm{.Math.}-\delta \right)-2{\alpha }^2{\beta }^2\mathrm{.Math.}\right]\left(\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right)\frac{d}{{\mathrm{dx}}_1}{F}_{35}=1+\frac{{\alpha }^2}{{\alpha }^2-{\beta }^2}\delta -\frac{2i{\alpha }^2}{{{\omega }^3\left({\alpha }^2-{\beta }^2\right)}^3}\left[\left(1+\delta \right){\alpha }^2-{\beta }^2\right]\hspace{1em}\hspace{1em}\left[\left({\alpha }^2+{\beta }^2\right)\left(\mathrm{.Math.}-\delta \right)-2{\alpha }^2{\beta }^2\mathrm{.Math.}\right]\left(\frac{d^2}{{\mathrm{dx}}_1^2}+\frac{d^2}{{\mathrm{dx}}_2^2}\right)\frac{d}{{\mathrm{dx}}_2}$

(29) FIG. 1 shows an example of the present invention in use. In this example, a simple 3D isotropic model with two layers was used. Of course, more complex anisotropic models can also be used. The uppermost layer is defined by velocities α=1500 m/s (compression wave velocity) and β=0 m/s (shear wave velocity). The uppermost layer may therefore represent sea water or the like. The lowermost layer is defined by velocities α=3000 m/s and β=750 m/s. The density is kept constant in both layers to ρ=1.0 g/cm.sup.3.

(30) As can be seen from FIG. 1, this example demonstrates that the method attenuates the dispersive shear waves appearing within the elastic part of the model (e.g. the lowermost layer) using the present method. A point source is placed at 10 meters depth in the z direction of FIG. 1 and at 0 meters offset in the x direction of FIG. 1. The modeling is done with second order time update and fourth order spatial differentiation operators (Virieux, 1986). The edges of the modeled wavefield are damped to avoid artificial reflections.

(31) FIG. 1 shows the forward modeled v.sub.1-component of the wavefield in a 2D section along the y-location of the source. FIG. 1 shows the v.sub.1-component relative to the horizontal x direction (the offset) and the vertical z direction (the depth). The wavefield component is extracted at 3 Hz. The modeling is carried on for the time needed for the dispersive shear-wave to propagate out of the model into the absorbing layers at the edges of the model. FIG. 1a shows the modeled v.sub.1-component with the dispersive shear-waves from the modeling; and FIG. 1b shows the modified, filtered, reconstructed {tilde over (v)}.sub.1-component where the dispersive shear waves are attenuated according to the present method.

(32) It should be apparent that the foregoing relates only to the preferred embodiments of the present application and the resultant patent. Numerous changes and modification may be made herein by one of ordinary skill in the art without departing from the general spirit and scope of the invention as defined by the following claims and the equivalents thereof.

REFERENCES

(33) Chapman, C. H. and J. O. A. Robertsson, 2013, Correcting an acoustic simulation for elastic effects: U.S. Pat. No. 8,437,219 B2. Chapman, C. H., J. W. D. Hobro and J. O. A. Robertsson, 2010, Elastic corrections to acoustic finite-difference simulations: 80th Annual International Meeting, SEG, Expanded Abstracts, 29, 3013-3017. Chapman, C. H., J. W. D. Hobro and J. O. A. Robertsson, 2014, Correcting an acoustic wavefield for elastic effects: Geophys. J. Int. doi: 10.1093/gji/ggu057 Hobro, J. W. D., Chapman, C. H., and J. O. A. Robertsson, 2011, Elastic corrections to acoustic finite-difference simulations: Plane wave analysis and examples: 73rd EAGE Conference and Exhibition, Expanded Abstracts, 1, 61-65. Moczo, P., Robertsson, J, O. A., and Eisner, L., 2007. The finite-difference time-domain method for modelling of seismic wave propagation. In Wu, R. S., and Maupin, V. (eds.), Advances in wave propagation in heterogeneous Earth, Vol. 48, Advances in Geophysics (ed. R. Dmowska). Oxford: Elsevier-Pergamon, pp. 421-516. Robertsson, J. O. A., and Blanch, J. O., 2011, Numerical methods, finite difference: In Encyclopedia of Solid Earth Geophysics: Springer. Virieux, J., 1986. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method: Geophysics 51, no. 4, pp. 889-901.