MODELLING

FIELD OF THE INVENTION

The invention relates finite difference (FD) wave equation modelling of, for example, seismic data, and in particular methods and systems for elimination or accretion of time dispersion in finite difference wave equation modelling of elastic and viscoelastic media (media with absorption).

BACKGROUND

Since approximately 1990, Finite Difference (FD) time-domain techniques have emerged as the primary means to computationally model many scientific and technological problems dealing with acoustic-, elastic-, and electromagnetic-wave propagation. In FD modelling, the solution to the wave equation is approximated using an FD method. The FD method is particularly useful for simulating wave propagation in the heterogeneous earth. It is used routinely to carry out huge computations for a myriad of applications, including exploration seismology where it is the workhorse in numerical modeling and wave-equation-based imaging and inversion algorithms because of its relatively low cost, in particular for low-order FD solution. When higher order, such as fourth order FD solutions are used the amount of temporal dispersion is reduced however the computational cost is greatly increased.

Consequently, there are efforts towards continued refinement in computational techniques for complex wavefield propagation, such as low-symmetry anisotropy, elasticity, and viscoelasticity. For an overview, the reader is referred to, e.g., Etgen, J., S. H. Gray, and Y. Zhang, 2009,“An overview of depth imaging in exploration geophysics”: Geophysics, 74(6), WCA5-WCA17, https://doi.Org/10.1190/1.3223188. Robertsson, J. O. A., J. O. Blanch, K. Nihei, and J. Tromp, eds., 2012,“Numerical modeling of seismic wave propagation: Gridded two-way wave-equation methods”: SEG Geophysics Reprints, Amundsen, L., and 0. Pedersen, 2017, “Time step n-tupling for wave equations”: Geophysics, 82(6), T249-T254. https://doi.Org/10.1 190/geo2017-0377.1 , and Ikelle, L. T., and L. Amundsen, 2018,“Introduction to petroleum seismology”: 2 ^nd edition, SEG, and references therein. The FD time-domain method is accurate and robust. However, one drawback of the FD time-domain method is its numerical temporal dispersion. For low-order FD- accuracy, the scheme has a numerical error which is more pronounced towards higher frequencies or at longer wave propagation distances, and this error is accumulated over time. Tal-Ezer, H., 1986, “Spectral methods in time for hyperbolic equations”: Siam J. Numerical Anal., 23,11-20, introduced the Rapid Expansion Method (REM) to reduce temporal dispersion. He showed that the REM has spectral accuracy in its time integration. For reverse time migration, the REM method is implemented in a time stepping manner as that described in Pestana, R. C. and P. L. Stoffa, 2010, “Time evolution of the wave equation using rapid expansion method”: Geophysics, 75(4), T121- T131 , and Tessmer, E., 2011 , “Using the rapid expansion method for accurate time stepping in modeling and reverse-time migration”: Geophysics, 76(4),S177-S185.

It has been shown that the numerical time dispersion is independent of spatial dispersion - Dablain, M. A., 1986,“The application of high-order differencing to the scalar wave equation”: Geophysics, 51 , 54-66, http://dx.doi.Org/10.1 190/1.1442040.

Stork, C., 2013,“Eliminating nearly all dispersion error from FD modeling and RTM with minimal cost increase”:75th Annual International Conference and Exhibition, EAGE, Extended Abstracts, 10.3997/2214-4609.20130478, pointed out that the numerical time dispersion in FD-modeling can be filtered from seismograms post-modeling using time- variable filter banks. Stork’s approach was further studied by, Dai, N., W. Wu, and L. H., 2014, “Solutions to numerical dispersion error of time FD in RTM”: 84 ^th Annual

International Meeting, SEG, Expanded Abstracts, 4027-4031 , http :/7dx.do i .o rg/ 10.1 190/segam2014-0858.1 and Li, Y. E., M. Wong, and R. Clapp, 2016 “Equivalent accuracy at a fraction of the cost: Overcoming temporal dispersion”, Geophysics, 81 (5), T189-T196.

Anderson, J. E., Brytik, V., and Ayeni, G., 2015,“Numerical Temporal Dispersion Corrections For Broadband Temporal Simulation, RTM and FWI”, in 85th Annual International Meeting, SEG, Expanded Abstracts, proposed a frequency-domain solution that requires frequency resampling before the transform back to the time domain.

US-B-8614930 (Zhang et al. (2013)) and Wang, M., and S. Xu, 2015, “Finite- difference time dispersion transforms for wave propagation”: Geophysics, 80, no. 6, WD19-WD25, http://dx.doi.Org/f 0.119Q/geo2Q15-0059.1 proposed to use a frequency time transform to seismograms post-modeling to eliminate the time-stepping dispersion. In detail, US-B-8614930, discloses a computer implemented method for accounting for temporal dispersion in low-order finite difference seismic propagation comprising: a. transforming a seismic dataset containing artifacts caused by temporal dispersion from time domain to frequency domain to obtain a frequency-domain seismic dataset; b. applying a frequency-domain time varying filter based on an effective phase velocity to the frequency-domain seismic dataset to obtain a filtered frequency-domain seismic dataset; and

c. transforming the filtered frequency-domain seismic dataset from the frequency domain to the time domain to obtain a time-domain filtered seismic dataset, wherein the artifacts caused by temporal dispersion have been reduced in the time-domain filtered seismic dataset.

US-B-8614930 also discloses using the frequency-domain time varying filter to prepare recorded seismic data for migration by a Reverse Time Migration (RTM) algorithm.

The approach of US-B-8614930 was redesigned in Qin, Y., Quiring, S., and Nauta, M., 2017,“Temporal Dispersion Correction and Prediction by Using Spectral Mapping”, 79th Annual International Conference & Exhibition 2017, EAGE, Extended Abstracts, Xu, Z., Jiao, K., Cheng, X., Sun, D., King, R., Nichols, D., and Vigh, D., 2017, “Time- dispersion filter for finite-difference modeling and reverse time migration”, in 87th Annual International Meeting, SEG, Expanded Abstracts, 4448-4452. https://do;.org/10.1190/segam2017-17790820.1 , and Koene, E. F. M., J. O. A Robertsson, F. Broggini, and F. Andersson, 2017, “Eliminating time dispersion from seismic wave modelling”: Geophysical Journal International, ggx563, https://doi.Org/10.1093/gii/ggx563.

Mittet, R., 2017, On the internal interfaces in finite-difference schemes”, Geophysics, 82(4), T159-T182, derived a temporal dispersion correction similar to that of Anderson et al. (2015) and, in addition, he provided a formal derivation.

From the above it is known that the numerical time dispersion is independent of the wavefield’s propagation path, the kind of elastic medium and any spatial modelling errors. Instead of eliminating time dispersion from seismograms, one may add time dispersion to seismograms. Obviously, for seismic exploration, numerical dispersion analysis is of interest for reverse time extrapolation and seismic inversion.

Viscoelastic properties of the medium have a profound influence on wave propagation. Modeling wave propagation in attenuative media with time-domain finite- differences has been described by many authors, e.g., Day, S. M., and J. B. Minster, 1984,“Numerical simulation of attenuated wavefields using a Pade approximant method”: Geophysical Journal International, 78, 105-118; Emmerich, H., and M. Korn, 1987, “Incorporation of attenuation into time-domain computations of seismic wave fields”: Geophysics, 52, 1252-1264; Carcione, J. M., 1993, “Seismic modeling in viscoelastic media”: Geophysics, 58, 1 10-120; Blanch, J. O., J. O. Robertsson, and W. W. Symes, 1995,“Modeling of a constant Q: Methodology and algorithm for an efficient and optimally inexpensive viscoelastic technique”: Geophysics, 60, 176-184; Cheng, N. A. C. H. Cheng, and M. N. Toksoz, 1996,“A time-domain finite-difference method with attenuation by a recursive algorithm”: Earth Resources Laboratory Industry Consortia Annual Report 1996- 12, Accessed January 2, 2018 http://hdl-handle.net/1721.1/75331 ; Bohlen, T., 2002, “Parallel 3-D viscoelastic finite difference seismic modelling”: Computers & Geosciences, 28, 887-899; and Zhu, T., J. M. Carcione, and J. M. Harris, 2013, “Approximating constant-Q seismic propagation in the time domain”: Geophysical Prospecting, 61 , 931— 940.

As the basis for the FD implementation the velocity-stress formulation of the system of differential equations is a popular choice (see Carcione et al.,1988; and Robertsson, J. O., J. O. Blanch, and W.W. Symes, 1994,“Viscoelastic finite-difference modelling”: Geophysics, 59, 1444-1456).

SUMMARY

According to a first aspect of the present invention there is provided a computer implemented method of determining a seismic wavefield in a region of a subsurface of the Earth by eliminating temporal dispersion from a finite difference, FD, wavefield generated from seismic wave propagation data. The method comprises obtaining the FD wavefield containing temporal dispersion by solving a discretized wave equation using the seismic wave propagation data; transforming the FD wavefield from a time domain to a frequency domain using a Fourier transform; and inverse transforming the FD wavefield from the frequency domain using an inverse Fourier transform, wherein the inverse Fourier transform uses a modified frequency based on a temporal discretization of the FD wavefield, to give the seismic wavefield in the time domain.

The modified frequency can be determined by fourier transforming the discretized wave equation in the time domain to obtain a discretized wave equation in the frequency domain; fourier transforming a wave equation in the time domain to obtain a wave equation in the frequency domain; and comparing the discretized wave equation in the frequency domain to the wave equation in the frequency domain. Said step of comparing may comprise, in the wave equation in the frequency domain, identifying a frequency term ; in the discretized wave equation in the frequency domain, identifying an equivalent term; and defining the modified frequency as being equal to said equivalent term.

The wave equation in the frequency domain may be

[w ² + L(x)]u(x, w) =— s(x, w),

where w is the frequency, L(x) is a system operator, ii(c, w) is the seismic wavefield in the frequency domain, and s(x, m) is a source term in the frequency domain, and wherein the discretized wave equation in the frequency domain is

is the equivalent term, u(x, w) is the finite difference wavefield

in the frequency domain, and §(c, w, M) is a finite difference source term.

Said steps of transforming the FD wavefield and inverse transforming the FD wavefield can be implemented by a time-step independent digital filter. The seismic wavefield may be given by

where u(n ) is the seismic wavefield at time step n, F(n, m ) is the digital filter, and u(m ) is the finite difference wavefield at time step m. Said digital filter may comprise a sum of Bessel functions.

Said digital filter may be given by

where

j _m is the m ^,h spherical Bessel function, and J _m is the m ^,h ordinary Bessel function. Alternatively, said digital filter may be given by

where / _m is the m ^,h ordinary Bessel function. According to a second aspect of the present invention there is provided a computer implemented method of determining a seismic wavefield in a region of a subsurface of the Earth, the region comprising anelastic media having a modulus, by eliminating temporal dispersion from a modified finite difference, FD, wavefield generated from seismic wave propagation data. The method comprises obtaining the modified FD wavefield containing temporal dispersion by solving a discretized wave equation having a modified modulus using the seismic wave propagation data; transforming the modified FD wavefield from a time domain to a frequency domain using a Fourier transform; and inverse transforming the modified FD wavefield from the frequency domain to the time domain using an inverse Fourier transform, wherein the inverse Fourier transform uses a modified frequency based on a temporal discretization of the modified FD wavefield, to give the seismic wavefield in the time domain.

The modified modulus in the frequency domain, evaluated at a given frequency, can be defined to equal the modulus in the frequency domain, evaluated at the modified frequency.

The method may comprise determining the modified modulus from an unrelaxed modulus, a relaxed modulus, a set of relaxation frequencies and a corresponding set of relaxation strengths. The modified modulus can have the same form as the modulus of the anelastic media but with a change in the unrelaxed modulus and a shift of the relaxation frequencies and corresponding relaxation strengths.

The modified mo

where M _m (x, t ) is the modified modulus in the time domain,

where M _u is the unrelaxed modulus, AM = M _U - M _r , where M _r is the relaxed modulus, w _; · is the set of (N) relaxation frequencies, a _; · is the corresponding set of relaxation strengths, At is a time step in the discretized wave equation, 5(t) is a Dirac delta function,

and H(t ) is a Heaviside step function. wDί may be much smaller than 1 (wAί « 1), where w is the highest frequency of interest of seismic waves recorded in the seismic wave propagation data. The steps of transforming the modified FD wavefield and inverse transforming the modified FD wavefield can be implemented by a time-step independent digital filter.

According to a third aspect of the present invention there is provided a computer implemented method of adding temporal dispersion to a seismic wavefield in a region of a subsurface of the Earth. The method comprises obtaining the seismic wavefield from seismic recordings; transforming the seismic wavefield from a time domain to a frequency domain using a Fourier transform, wherein the Fourier transform uses a modified frequency based on a temporal discretization of a finite difference, FD, wavefield; and inverse transforming the seismic wavefield back to the time domain using an inverse Fourier transform using an unmodified frequency to give a wavefield, or source signature, containing temporal dispersion.

Said steps of transforming the seismic wavefield and inverse transforming the seismic wavefield can be implemented by a time-step independent digital filter. The wavefield containing temporal dispersion may be given by

where u(n ) is wavefield containing temporal dispersion at time step n, G(n, m) is the digital filter, and u(m) is the seismic wavefield at time step m. Said digital filter may comprise a sum of Bessel functions.

Embodiments of the invention provide a framework for the elimination of temporal numerical dispersion from wave-equation seismograms FD modelled with second-order time integration. The framework can be generalized to FD modelling with higher-order time integration. Numerical dispersion is the separation of different Fourier components of a FD approximation into a train of oscillations that travel with different speeds. It occurs whenever the dispersion relation for the difference approximation is nonlinear.

In accordance with embodiments of the invention, the Fourier transform can be used to relate the solutions of the wave equation and the temporal FD wave equation. Each of the two wavefields is a time-frequency transformation of the other. From this theory, we developed post-modeling filtering transforms for the elimination of the temporal dispersion caused by the inaccuracy of the second-order FD approximation to the time derivatives in the wave equation. The approach is valid for solving the two-way second- order wave equation as well as the system of first-order partial differential equations for stress and particle velocity. Similarly, it is valid for the one-way wave equation. The numerical dispersion correction depends only on the simulation time step, the frequency and the overall simulation time. For elastic media, the dispersion correction is independent of the medium and the propagation path.

For anelastic media, the dispersion correction can be applied when the seismograms are modelled in a frequency-modified model. In the time domain, the modification is applied to the unrelaxed modulus, and the relaxation frequencies and their strengths. The modification depends on the product of relaxation frequency and time step. The schemes for elastic and anelastic media permit the use of low-order FD in time to achieve high numerical accuracy in wavefield propagation. It is also demostrated how numerical time dispersion can be added to seismograms for the wave equations. This transform is the bandlimited inverse of the transform that eliminates dispersion. It is also shown that both transforms can be implemented by using time-step independent digital filters that are pre-computed from sums over special real-valued functions. The digital filter for dispersion accretion is equal to the bandlimited inverse of the filter for elimination of dispersion.

DESCRIPTION OF THE DRAWINGS

Embodiments of the invention will now be described in more detail, and by way of example only, with reference to the accompanying drawings, in which:

Figure 1 a and 1 b are graphs of (a) the time derivative of Ricker signatures and (b) the corresponding normalized amplitude spectra;

Figure 2a and 3b are graphs of the results for acoustic wave equation: (a) the reference, exact solution, and FD modeled trace contaminated with numerical dispersion, and (b) dispersion corrected solution and the error;

Figure 3 is a graph of the Q-factor for the reference model;

Figure 4a and 4b are graphs of the results for viscoacoustic wave equation: (a) the reference, exact solution, and FD modeled trace contaminated with numerical dispersion, and (b) dispersion corrected solution and error;

Figure 5a and 5b are graphs of the results for viscoacoustic wave equation: (a) reference, exact solution, and FD t - i modeled trace contaminated with numerical dispersion, and (b) dispersion corrected solution and error;

Figure 6a is a plot of a digital filter F n, m) that eliminates temporal numerical dispersion from dispersion-contaminated FD wave-equation seismograms. F(n, m), applied to the dispersion-contaminated data u m ) over all time samples m, predicts the dispersion-free data at time sample n;

Figure 6b is a plot of the digital filter G(n, m ) that adds temporal numerical dispersion to dispersion-free wave-equation seismograms or source signatures. G(n, m), applied to the dispersion-free data u m ) over all time samples m, predicts the temporally dispersive wave equation data for second-order FD time integration at time sample n; and Figure 6c is a plot of product of the filters in Figures 6a and Figure 6b, represented as matrices, becomes a symmetric Toeplitz matrix with entries being sinc-function elements.

DETAILED DESCRIPTION

Time integration of wave equations can be carried out with explicit time-stepping using a finite-difference (FD) approximation. The wave equation is the partial differential equation that governs the wavefield that we want to solve for. The FD approximation gives us another partial differential equation with temporal discretization - the one we solve in the computer for the FD wavefield. This approximation in numerical modeling produces a wavefield contaminated with temporal dispersion, particularly at higher frequencies. We show how the Fourier transform can be used to relate the two partial differential equations and their solutions. Each of the two wavefields is then a time-frequency transformation of the other. Firstly, this transformation allows temporal dispersion to be eliminated from the FD wavefield, and secondly, it allows temporal dispersion to be added to the exact wavefield. The two transforms are inverse operations in a bandlimited sense. The transforms can be implemented by using time-step independent digital filters that are pre computed from sums over special real-valued functions (e.g. Bessel functions).

For anelastic materials, the effect of numerical time dispersion in a wavefield propagating in a medium needs special treatment. Dispersion can be removed by using the time-frequency transform when the FD wavefield is modeled in a medium with frequency-modified modulus relative to the physical modulus of interest. In the rheological model of the generalized Maxwell body, the frequency-modified modulus in the time domain can be derived approximately by introducing mathematical approximations to the modulus of the physical wave propagation. The modified modulus then has the same form as the physical modulus, and can be implemented as changes in the unrelaxed modulus and shifts of the relaxation frequencies and their strengths of the physical modulus.

The discussion below starts by defining the wave equation and the corresponding FD wave equation. Then, it is shown that given a solution to one of the equations, the related solution can be derived to the other equation by employing a Fourier transform. This straightforward procedure, in the frequency domain, gives the relationship between the numerically dispersed wavefield and the dispersion-free wavefield. Therefore, the numerically modelled wavefield can be transformed to the frequency domain by using a Fourier transform. Then, it is inverse Fourier transformed to time where it equals the non- dispersive wavefield.

For anelastic media, the material depends on frequency, and the frequency-time transform cannot be applied directly. However, we show that by solving the wave equation for a frequency-modified anelastic medium, producing a modified wavefield, the frequency-time transform can be applied to this numerical wavefield, outputting the sought-after, numerical dispersion-free wavefield for the anelastic medium of interest.

First, Fourier transforms are applied to the wave equation and the FD wave equation so that two wavefields are governed by the same differential operator and the same source field. Then, the wavefields for the two differential equations are equal by definition, such that the procedure allows us to identify and develop transforms for eliminating and adding numerical dispersion to seismograms. Dispersion elimination is the bandlimited inverse process of adding dispersion, or dispersion accretion. In Appendix A, we write the transforms as digital filters that can be pre-computed for any later use. The digital filters are computed as sums over known real-valued functions.

Second, the above approach leads to an understanding of how the dispersion correction can be applied to seismograms in anelastic media when the model for the anelastic modulus is modified. The modification to the modulus is straightforward to define in the frequency domain. However, modeling is done in time domain, and the challenge is to derive the modulus in the time domain. To this end, mathematical approximations are introduced to the frequency-domain modulus in the low-frequency regime which allows inverse Fourier transformation to the time domain.

In the following, firstly the theory for elastic media is developed and secondly the theory for anelastic media. Thirdly the procedures on FD difference modeled data in 1 -D homogeneous acoustic and viscoacoustic media are illustrated.

THEORY FOR WAVE EQUATION IN ELASTIC MEDIA

Let time be denoted by t and the spatial coordinate by x. The wavefield is denoted by u. The fields are vectors for the elastic wave equation and scalars for the acoustic- wave equation. The equations for wavefield propagation can be described by the second- order wave equation, or the system of first-order velocity-stress equations (Virieux, 1986). These equations have identical time dispersion error when the first-order system of equations is solved in a leap-frog manner and the second-order time derivative is approximated by a second-order accurate FD operator derived from a Taylor series (equation 2b). Therefore, we consider the generic formula for the second-order wave equation in time expressed as:

d ² u(x, t )

= L(x)u(x, t ) + s(x, t), #(1)

dt ²

where L(x) is the system operator that contains material parameters and spatial derivatives, and s is the source term. For acoustic wave propagation, L = c ² V ² , where c = c(x) is the velocity and V ² is the Laplacian. The source is often considered to consist of N point sources with signatures a,· at locations x , then s(x, t) = a _j (t)S(x - x _; ·).

The numerical solution of the wave equation (1 ) is discretized in time,

D ² U(X, t)

= L(x)u(x, t) + s(x, t, At) , #(2 a)

Dt ²

where the source s is yet to be determined. We shall see that it can easily be calculated from s.

The classic second-order accurate temporal discretization is:

where At is the time step. Respecting the Nyquist-Shannon sampling theorem, ύ is band- limited in frequency so that W(w) = 0 for all w ³

Here, w= circular frequency (referred to herein simply as“frequency”).

The problem at hand is: We would like to know the solution u x, t) (i.e. the wavefield) to the physical wave equation (1 ), but with the temporal discretization in equation 2 we obtain the numerical wavefield u{x, t). Is there an exact relationship between the physical and numerical wavefield? To answer this question, we analyse the wave equations 1 and 2b in the frequency domain. We define the forward Fourier transform:

with inverse

By transforming equations (1) and (2a) to the frequency domain using 3a, we obtain

[w ² + L(x)]u(x, w) =— s(x, w)#(4)

and

[w ² + L(x)]u(x, w) =—s(x, w, At), #(5) respectively, where

Comparing equations 4 and 5, we see that for the 2nd-order accurate FD scheme, the frequency is modified by the“modification factor” ^ = sinc ^). The modified frequency is based on the temporal discretization of the FD wavefield u(x, t). The modified frequency depends on the time step At of the discretized wave equation. This modification of frequency is well known and leads to numerical time dispersion via a frequency- dependent phase velocity. For the purpose of understanding the time dispersion problem, instead Fourier transform the wave equation 1 by using:

f

F(0) = I dt qcr(ίwί) /(t) ,

which yields

[ώ ² + L(x)]u(x, ) =— s(x, ). #(8)

This non-linear scaling of the frequency axis produces the wave equation 8, which has the same differential operator as the FD wave equation 5, but with differing sources. For the yet-to-be determined wavelet, we make the choice:

s(x, w, At) = s(x, w)#(9a)

with time domain expression:

The time integral maps the original source signature of s to the frequency . At this frequency, the source signature equals the source signature of s at frequency w, such that a conventional inverse Fourier transform produces the FD source signature in the time domain. The effect of the transform in equation 9b is to add temporal dispersion to the FD source signature. In Appendix A, we show that the transform in equation 9b can be implemented by digital filtering.

Now, we have the solution for elastic media to our problem. Since the wavefield u defined by equation 5 has the same system operator and the same source as that for the scaled wavefield u in equation 9, the two wavefields necessarily are identical; whence, u(x, w) = u(x, ) . #(10)

Equation 10 is a key result. We have arrived at the point where, given a solution to one of the wave equations, we can derive the solution to the other wave equation by employing the Fourier transformation. Recall that u is the FD wavefield in time and space. Thus, since u is known for all times, we can Fourier transform u to the frequency domain. Here, this wavefield equals the temporally dispersion-free wavefield u at frequency . Finally, an inverse Fourier transform gives the temporally dispersion-free wavefield in time.

In short, the FD wavefield is not equal to the actual (seismic) wavefield since it includes temporal dispersion. By Fourier transforming the FD wavefield to the frequency domain, it equals the actual wavefield in the frequency domain (equation 10) as a function of the modified frequency (equation 6). By inverse transforming the FD wavefield from the frequency domain using the modified frequency it gives the actual wavefield in the time domain.

The formal equation for elimination of numerical temporal dispersion reads:

In Appendix A, we show that the dispersion correction can be implemented as a time-step independent, digital filter applied to the data. The filter can be pre-computed and applied to any FD data modelled with second-order time integration of the wave equation.

We observe that the dispersion elimination equation 11 is independent of the system operator L(x). Therefore, numerical dispersion related to the temporal discretization of the wave equation can be eliminated by using equation 11 for any modelled wavefield in elastic media.

It is of interest for reverse time extrapolation to add numerical dispersion to the exact wavefield in the case that u is known in time-space. The relation in equation 10 shows that the formal equation for adding numerical temporal dispersion then reads

Equation 12 represents the formal procedure for adding dispersion to data. Its digitalization using filters is given in Appendix A. Observe that the integrals that are involved in equation 12 are exactly those that are present in equation 9b which computes the FD source signature from the original signature. In this respect, the procedure in equation 9b adds temporal dispersion to the original source signature. Adding dispersion can be applied to a) measured or recorded seismic data (i.e. the actual seismic wavefield), or b) source signatures in FD modeling. Importantly, the source signatures can be seismic recordings themselves as well as designed signatures (e.g. Ricker signatures).

An alternative dispersion correction transform can be derived from equation 10 by letting w ® , which also corresponds to applying a shift of variables in

equation 11 , yielding

Equation 13 is similar to equation B-17 in Mittet (2017) and equation 15 in Koene et al. (2017). We notice from equations 1 1 and 13 that the temporal dispersion correction depends only on the simulation time step, the frequency and the overall simulation time. A disadvantage of the approach in equation 13 compared to the one in equation 1 1 regards sampling requirements. The maximum frequency in the argument of the arcsin

function is is the Nyquist frequency. This puts a limit on the frequency content of the seismogram in the simulation. The support of w can be considered as the frequency range up to ^ so that / = f _N , the maximum frequency that can be recorded without aliasing. Furthermore, equation 13 is not necessarily a convenient starting point for deriving digital filters, in contrast to equation 11. The approach disclosed herein to digital filtering is based on utilizing the Jacobi-Anger expansion, which is not applicable for equation 13.

Equation 1 1 has not been derived previously. An advantage of equation 11 is that it lends itself to deriving analytic forms for the digital filter. Since the digital filters are based on equation 11 the filter can be written in terms of an integral which has analytic solution; in an embodiment the digital filter is obtained as a sum of known (Bessel) functions. Previously, starting with equation like equation 13, the integral obtained for the filter has no analytical solution, and therefore must be solved numerically. Therefore, equation 11 is preferred to equation 13 and prior equations of this sort.

THEORY FOR WAVE EQUATIONS IN ANELASTIC MEDIA

In this section is provided a general analysis of the problem of removing numerical time dispersion, again based on considerations in the frequency domain. It is shown that numerical dispersion can be removed if we propagate waves in a frequency-modified medium of the physical medium of interest. To implement this understanding, the challenge is to inverse Fourier transform the modified modulus of the medium from frequency to time. To this end, a series expansion of the modulus is applied which allows inverse Fourier transform of each term in the series. Finally, the procedure can be implemented in an FD simulation for example using an approach such as that presented by Cheng et al. (1996).

General analysis

For anelastic media, the system operator ! = L(x, t) depends on time since some of the material parameters depend on time. For the sake of simplicity, we consider viscoacoustic wave propagation for the pressure field governed by:

where the symbol * denotes temporal convolution. By Fourier transforming equation 14 to the frequency domain, we obtain:

1

w ² + M(x, w)n V u(x, w) =— a(w)d(c). #(15)

p(x)

When we discretize equation 14 using second-order time differences (see equation 2b) we obtain the numerical wavefield u{x, t). In the previous section, for elastic media, we found the relation in the frequency domain between W and u by redefining the frequency in the numerical wave equation. For anelastic media, we anticipate that a modification of frequency also will affect the frequency-dependent material parameters. Therefore, the dispersion elimination or dispersion addition equations 11 and 12, respectively, should not be applied without considerations to the numerical wavefield u{x, t) modelled in anelastic media with the system operator L(x, t). However, consider the wave equation for a modified field in a frequency-modified medium, having a modified modulus:

M _m (x, w ) = M(x, w) , #(16 a)

given in time domain by applying the inverse Fourier transform

Discretizing equation 14 in terms of the modified modulus, gives a wave equation for the modified field:

D ² u _m (x, t ) 1

= M _m (x, t) * V Vu _m (x, t ) + a(t)<5(x)#( 17)

Dt ² p(x) in the frequency domain becomes:

where we require that u _m (a >) = 0 and a(w) = 0 for all The modified frequency w is the same as for elastic media.

We now follow the procedure of redefining the Fourier transform. Applying the transform in equation 7 to the wave equation 14 gives:

2 + M(x, )n u(x, w) =— a(w)5(c) . #(19)

p(x)

We observe that u _m and u now are defined by the same system operator and source ά(w) = a( ) since we have modified the modulus of equation 18 according to equation 16a, yielding:

u _m (x, w ) = u(x, w) . #(20)

Recall that is modelled in time-space, not with the modulus of interest M(x, t), but with the frequency-modified modulus M _m (x, t); has numerical dispersion due to time discretization. The objective is to find the wavefield u in the medium defined by the modulus M(x, t). This wavefield u should see no effect of the numerical dispersion. Equation 20 gives the relationship between the wavefields. Since u _m is known for all times, we Fourier transform u _m to the frequency domain. At these frequencies, u _m equals the temporally dispersion-free wavefield u at frequency w. Finally, an inverse Fourier transform gives the temporally dispersion-free wavefield in time. The formal equation for elimination of numerical temporal dispersion reads:

So far, the theory is good only if we can find an expression for M _m (x, t ) which is implementable in FD schemes. Therefore, we make a closer look at the viscoelastic modulus.

The complex viscoelastic modulus

Finite-difference experts tend to implement realistic attenuation in the time-domain methods by using a generalized Maxwell body or generalized Zener body. Since the rheology of the two models is one and the same, in the present embodiment the complex viscoelastic modulus of the Maxwell body in the form : N

-

M{c, w ) = M _u (x)— DM(c) i Y , #(22 a)

Z_„ '„ w,

w + ϊw

i= i ^J

where M _u is the unrelaxed modulus, AM = M _u - M _r is the relaxation of modulus, M _r being the relaxed modulus, in which the relaxation spectrum consists of N single peaks of strength a,· at discrete relaxation frequencies w _;· . The first part is independent of frequency, whereas the second part has a simple pole at w = -ia _j . In the time domain, the modulus follows by evaluating the inverse Fourier transform:

which yields:

where S{i) is the Dirac delta function and H{i) is the Heaviside unit step function. This well-known result follows by the use of Cauchy’s residue theorem.

By virtue of equation 16, the frequency-modified modulus in the time domain is the solution of:

The integral has no solution as it stands with infinite integration limits. However, to derive the wave equation in the frequency domain (equation 19) we used the Fourier transform in equation 7, which has a limitation on the maximum allowable frequency. For the dimensionless integration variable w in equation 23a, it implies |w| Proceeding with the geophysicist’s customary cavalier attitude towards problems of existence, we make a series expansion of the function (sin( i) + i£) ^_1 in equation 23a around w =

—i arcsinh(&), b = which yields the remarkable appearance of the series:

1 1 — b ² ) arcsinh(&)

+ iO(arcsinh(¾)) ² ,

Vl + b ^{2 w} + i arcsinh(Z ⁾ ) 3

12(1 + b ² )2

#(23 b) where the first term has a simple pole at w = -i arcsinh(6) . Since b is a real number, the inverse hyperbolic sine satisfies arcsinh (b) = log(x + Vx ² + l). In the series, we have disregarded terms that are proportional to ( ίw) ^h ,h = 1,2, ..., which will produce terms in equation 23a proportional to 5 ⁽ⁿ⁾ , the n-th derivative of 5; such terms are not present in equation 22c either. Inserting the series in equation 23b into equation 23a, and assuming that the relaxation frequencies are inside the allowed frequency interval, yields the modified time-domain modulus:

For seismic frequencies, the number (M _u - M _u ’)/M _u is small. However, for propagation over long distances, the effect is notable in the arrival time of the wavefield.

Observe that M _m in equation 23c has the same form as M in equation 22c but with a change in the unrelaxed modulus and shift in relaxation frequencies and their strengths. Therefore, FD experts who implement viscoelastic schemes with the modulus M easily can simulate waves in media with modulus M _m .

The 1-D visco-acoustic wave equation

Now, let’s introduce the viscoacoustic wave equation for numerical simulations. For clarity, we drop the prime-superscript on variables. Further, for simplicity we consider 1 -D wave propagation. We introduce Newton’s law:

where p is the density. For a linear isotropic viscoacoustic material, the time-derivative of the stress-strain relation can be written:

where

ds(x, t) dv(x, t)

. #(24 c)

dt dx

Inserting equation 24a into equation 24b, yields

where

A Fourier transform of equations 24a, 24d and 24e gives, respectively,

Equations 25a-c lead to the wave equation in the frequency domain :

d i d

w ² + M(x, s(c, w) = 0 , #(26)

dx pipe) dx _.

where M(x, w) is given in equation 22a.

The corresponding FD equations are

where equation 27c is derived in Appendix C. Here, 0 < d < 1 ; a parameter of choice dv(x

since equation 27c is derived with the assumption that ^ ' is constant within each time step At. We select 5 such that the FD solution fits the analytical solution to the wave equation in the best possible way. To anticipate the course of events, we reveal that our choice is d = 1. Since time-stepping of v is at t ± y (as seen in equation 27a) this choice implies that v in equation 27c must be interpolated as v(t) = (t + ) + v(t— y)] .

A Fourier transform of the FD equations 27a-c gives:

and allows the construction of the FD wave equation in the frequency domain:

d i d

2 + M(x, w) s(c, w) = 0 , #(29 a)

dx p{x) dx

with modulus:

To obtain the best fit to the modulus of the exact viscoacoustic wave equation, we select 5 = 1. For a>At « 1 the modulus in equation 29b can be approximated by

In this case, with the above assumptions, the modulus of the FD wave equation is approximately equal to the modulus of the exact viscoacoustic wave equation.

The 1-D viscoacoustic wave equation for the frequency-modified medium

To apply the temporal dispersion correction to FD modelled seismograms, equation 16a proposes that we do FD modelling in a frequency-modified medium with change of frequency from w to . When the medium can be described by the complex viscoelastic modulus of the form in equation 22b, it is found that the frequency-modified modulus can be expressed as the time-domain modulus in equation 23c. However, the FD analysis reveals that the modulus seen by FD propagation is that in equation 29b. Ideally, this is the modulus that should be transformed to the time domain when replacing w by w. Under the assumption that wDί « 1 (where w is the highest frequency of interest), as used in equation 29c, the FD related modulus and the true modulus are equal. Therefore, as long as we do not violate wDί « 1 (e.g. wDί < 1/3, as in geophysics 1/3 can be considered to be much less than 1 ) in the numerical modelling, we can expect that the frequency-modified modulus in equation 23c is a good and simple choice for the modulus that enters the numerical modelling. While the model is not fully perfect for high frequencies, being aware of this limitation, will help in making informed decisions to improve modelling with higher frequencies when that challenge surfaces in the future. In other words the discrepancy between the calculated wavefield and the actual (measured) wavefield increases as wDί increases.

EXAMPLES

To demonstrate temporal dispersion elimination by use of equation 1 1 for an acoustic medium and equation 21 for a visco-acoustic medium. The model is 1 -D and homogeneous.

First, to obtain a reference solution, we solve the exact acoustic wave equation B- 1 a (Appendix B) for velocity c = 2000 m/s. We solve for the wavefield in frequency domain by using equation B-1c. The solution in time domain is found by applying the inverse Fourier transform. The source signature shown in Figure 1 is the time-derivative of a Ricker signature (Reference numeral 1 is used for the source signature for the FD equation and reference numeral 2 is used for the source signature for the original differential equation), where the Ricker signature has a central frequency of 30 Hz. Then the recorded seismogram u ^(ex) (t) at distance x = 5 km equals a time-delayed Ricker signature shown by the line numbered 4 in Figure 2a. Second, we solve the FD time domain wave equation for a temporal sampling interval At = 2 ms. The solution W(w) is given analytically by equation B-2d, and an inverse Fourier transform gives the time domain seismogram at distance x = 5 km shown by line number 1 in Figure 1 b. Observe the effect of temporal discretization, which for the second-order time difference in equation 2b leads to numerical dispersion due to phase velocity increasing with frequency. Third, we apply the temporal dispersion correction transform in equation 11 to the numerically dispersed trace obtained in the second step. As shown in Figure 2b, where the dispersion corrected trace u(t) is displayed in solid line (5) and the error u t ) - u ^(ex) (t) displayed in thin line (6), the correction is close to perfect.

The test of the temporal dispersion correction algorithm for viscoacoustic media uses data modelled analytically in the frequency domain and FD time-space domain with a spectral method for the spatial derivative. The medium is homogeneous with velocity

= 2000 m/s and density p = 1000 kgm ^-3 . The wavelet that is fed into the modelling scheme, is the same as in the examples above, as well as the offset, x = 5 km. Three discrete relaxation frequencies are chosen at 1 , 10 and 100 Hz, with peak strengths for attenuation determined so that the Q-factor is close to frequency-independent at Q = 80 as displayed in Figure 3 over the frequency range 0-100 Hz. The reference seismogram, being the analytic solution of the wave equation, given by equation B-3d, is displayed as line 8 in Figure 4a.

Now, to apply our approach to dispersion correction, we need to modify the wave propagation model according to equations 23d-f. First, we apply analytic modelling. We choose a temporal sampling interval, At = 2 ms. The modified relaxation frequencies become, 1 , 9.99 and 94.37 Hz. The change in the unrelaxed modulus is 0.23%. It may seem like a small number, but it has an effect on the travel time over long distances of the wavefield. The solution it _m (a>) is given analytically by equation B-4e, with the proper changes to the modelling parameters as described by equations 23d-f. An inverse Fourier transform gives the time domain seismogram at distance x = 5 km as shown by line 7 in Figure 4a. Then, we apply the temporal dispersion correction transform in equation 22 to the numerically dispersed trace obtained in the previous step. As shown in Figure 4b, where the dispersion corrected trace u(t ) is displayed in solid line 9 and the error u(t ) - u^ ^ex t) displayed in thin line 10, the correction is fairly good. The reason that the result is not perfect, as in the example of the acoustic wave equation, is that the modified model is introduced to the time domain calculation under the assumption that wAί « 1; see equations 29b and 29c.

So far, we have evaluated the dispersion correction methods to FD seismograms modelled with perfect computation of spatial derivatives by Fourier methods. In the final numerical example, we apply the dispersion correction to time-space FD computed data. The FD scheme is that outlined in equations 27a-c, and is implemented with fourth-order spatial derivatives. For numerical stability, we use At = 1.39 ms, which has the effect to produce less numerical temporal dispersion compared to the previous example. We make the same displays as in Figure 4, now plotted in Figure 5. Again, the correction for numerical temporal dispersion is good.

APPENDIX A: DIGITAL FILTERS FOR ELIMINATION AND ACCRETION OF NUMERICAL TIME DISPERSION

Equation 1 1 gives the formal procedure for the dispersion correction. Here, we show how the dispersion correction can be applied as a filtering process. When we interchange the integrals in equation 1 1 , we obtain

where

In the following, we drop the space coordinate x. We discretize the time variables, t = nAt , t' = mAt, and we obtain

The integral in A-2 can be solved analytically in various ways, for example it can be expressed in terms of Beta functions B(x, y) (see Gradshteyn, I. S. and I. M. Ryzhik, 2007, “Table of Integrals, Series, and Products: Seventh edition, Academic Press”, integral 3.892.2) so that the filter F(n, m) can be expressed in at least two forms, one of which is:

with j _m and J _m being spherical and ordinary Bessel functions, respectively, and ^F _c and f _s being coefficients

Special cases are 2 cos

F(0, m ) , #(A - 4d

(1— 4th ² )p

sin(2n)

F(n, 0) #(A - 4e)

Equation A-4a has the alternative solution:

The filter in equations A-4 is independent of the time step interval At, and once computed, it can be reused for any dispersion correction up to time T = nAt. F(n, m ) is displayed in Figure 6a.

Equation 12 represents the formal procedure for adding dispersion to data. It is the bandlimited inverse action of eliminating dispersion. Equation 12 has the same structure as equation 9a which adds dispersion to the original source signatures before they are used in FD modelling. Equation 12 can be written:

where

Again the integral can be solved analytically in several ways, which leads to the digital filter examples:

7 -1

G(n, m ) = J _2n (2m ) + > g _s (n, m, k ) j _n-k (2m) #(A - 7a)

k= 0

where g _s is the coefficient

with special cases are G(0, m) = J ₀ (2m) . #(A - 7 d)

An alternative formula is:

Once computed, the filters in equation A-7 can be reused to add dispersion on any time records up to time T = nAt. G(n, m) is plotted in Figure 6b. The filters F and G are matrices. Their matrix product FG = S becomes a symmetric Toeplitz matrix with entries being sinc-function elements:

2

S(n,m) =— sinc(2(n— m) ) #(A— 9)

p

The product of the two matrices is displayed in Figure 6c. In the limit that At goes to zero, it can be shown that FG = S = /, where / is the identity matrix.

APPENDIX B: ANALYTICAL SOLUTIONS

We derive analytical solutions to various 1 -D wave equations. These solutions are good to use to verify the approaches to elimination of numerical temporal dispersion.

The wave equation for acoustic media

The 1 -D wave equation for the stress field s in a homogeneous medium of velocity c is given by equation 1

In the frequency domain, equation B-1 a reads

with solution

Since the factor— corresponds to temporal integration, we observe that when the wavelet a(w) corresponds to the time derivative of a Ricker signature, then the solution outputs in time domain a time-shifted Ricker wavelet.

The wave equation with temporal numerical dispersion for acoustic media The frequency-domain FD wave equation in 1 -D follows from equation 5 as

d ²

2 + c ² s(c, w) =— a(w, At)S(x) #(B— 2 )

dx ²

It can be written as a wave equation having the frequency-dependent phase velocity c(o At) = - #{B - 2b)

sine

which increases with increasing frequency, and the amplitude-scaled wavelet

The solution of equation B-2a then is

The wave equation for viscoacoustic media

The frequency-domain 1 -D wave equation for the stress field s in a homogeneous viscoacoustic medium of modulus M and density p is given in equation 26 as

where M(w) is given i

By introducing the frequency-dependent phase velocity

equation B-3a obtains the solution

a(w) (ft) !

s(c, w) =— exp . #(b - 3d)

2ίw e(w)

The quality factor Q is given as

The wave equation with temporal numerical dispersion for viscoacoustic media

The frequency-domain FD wave equation in 1 -D is given by equation 29a as

with modulus defined by equation 29b as ft)

When we introduce the modulus

we get the modified phase velocity

Equation B-4a then has the solution

with a given by equation B-2c.

The wave equation with numerical dispersion for frequency-modified viscoacoustic media

The wavefield solution is given by formulas B-4e, where M _u ® M _u ' , a _j ® a _j ' and w _;· ® w) according to equations 23d, 23e, and 23f, respectively.

APPENDIX C: DERIVATION OF EQUATION 27c Here, we derive equation 27c of the description from the definition of X _j (x, t ) in equation 24e. Letting t = kAt in equation 24e, the equation can be evaluated in a discrete manner as:

dn(c,t)

We assume that dx is constant within each time step At and choose for mAt < t < ( m + l)Af.

dv(x, t) 3v(x, (m + 5)Dί)

, #(C - 2)

dx dx

where 0 < d < 1; 5 is chosen so that the FD solution fits the analytical solution in the best possible way. By performing the integral over t, we get

Since X _j holds a summation of exponential functions, it can be updated recursively (Luebbers, R., F. P. Hunsberger, K. S. Kunz, R. B. Standler, and M. Schneider, 1990,“A frequency-dependent finite difference time domain formulation for dispersive materials”: IEEE Trans. Electromagn. Compat., 32, 222-227; Cheng et al. , 1996). With some simple algebra, we obtain the recursive formula:

which equals equation 27c.