SciELO - Scientific Electronic Library Online

 
vol.8 número2INTEGRAL MODELLING OF PROPAGATION OF INCIDENT WAVES IN A LATERALLY VARYING MEDIUM: AN EXPLORATION IN THE FREQUENCY DOMAININTERPOLATION AND DENOISING OF SEISMIC SIGNALS USING ORTHOGONAL MATCHING PURSUIT ALGORITHM: AN APLICATION IN VSP AND REFRACTION DATA índice de autoresíndice de assuntospesquisa de artigos
Home Pagelista alfabética de periódicos  

Serviços Personalizados

Artigo

Indicadores

Links relacionados

  • Em processo de indexaçãoCitado por Google
  • Não possue artigos similaresSimilares em SciELO
  • Em processo de indexaçãoSimilares em Google

Compartilhar


CT&F - Ciencia, Tecnología y Futuro

versão impressa ISSN 0122-5383

C.T.F Cienc. Tecnol. Futuro vol.8 no.2 Bucaramanga jul./dez. 2018

http://dx.doi.org/10.29047/01225383.80 

Artículos originales

NUMERICAL CONSIDERATIONS ON THE MODELING OF SOURCE AND BOUNDARY CONDITIONS FOR THE FREQUENCY DOMAIN VISCO-ACOUSTIC WAVE EQUATION SOLUTION

CONSIDERACIONES NUMÉRICAS EN EL MODELAMIENTO DE LA FUENTE Y CONDICIONES DE FRONTERA EN DOMINIO DE FRECUENCIA PARA LA SOLUCIÓN DE LA ECUACIÓN DE ONDA VISCO-ACÚSTICA

Sheryl Avendañoa  b  *  , J.-C Muñoz-Cuartasa  b 

a Group of Computational Physics and Astrophysics.

b Group of Tomography and Seismic, Instituto de Física, Facultad de Ciencias Exactas y Naturales, Universidad de Antioquia, Calle 70 No. 52-21, Medellín, Colombia.

ABSTRACT

Seismic modeling is an important step in the process used for imaging Earth subsurface. Current applications require accurate models associated with solutions of the wave propagation equation in real media. Unfortunately, it is common not to find in the technical literature deep discussions on the impact of specific details associated with the physical modeling of some crucial ingredients of the process, such as seismic source term and boundary conditions. In this paper, we discuss some issues related to the modeling of wave propagation in visco-acoustic media using finite differences. We focus our attention on two major elements of the modeling problem that are associated to the source term and the boundary conditions. We show that the source term can be modeled using a scale parameter that controls the spread of energy and shows that this parameter is a function of frequency and position of the source. As to boundary conditions, we show that Perfectly Matched Layer (PML) parameters are also frequency dependent. For both cases, seismic source scale parameter and PML model parameters we provide values and functions that optimize the performance of the approach for problems where visco-acoustic wave propagation is required. Frequency domain Full Waveform Inversion (FWI), or Reverse Time Migration (RTM) processes that depend fundamentally on the appropriate modeling of the wave-field are potential fields of application of these results.

Key words: Seismic attenuation; Wave propagation modeling; Visco-acoustic medium; seismic source; PML boundary conditions

RESUMEN

El modelamiento sísmico es un ingrediente importante en el proceso de construcción de imágenes del subsuelo. Aplicaciones actuales requieren soluciones precisas asociadas con la solución de la propagación de ondas en medios realistas. Desafortunadamente no es común encontrar en la literatura técnica, discusiones y análisis suficientemente profundas asociadas al impacto detalles específicos asociados con el modelamiento físico de ingredientes claves para este proceso como lo son el término de la fuente y las condiciones de frontera. En este trabajo se discuten algunos de esos detalles relacionados con el modelamiento de la propagación de ondas en un medio visco-acústico usando diferencias finitas. Nuestra atención se enfoca en dos aspectos importantes del modelamiento asociados con el término de la fuente y las condiciones de frontera. Se muestra que el término de la fuente puede ser modelado a través del uso de parámetros de escala que controlan la distribución de energía y se muestra que este parámetro es una función de la frecuencia y de la posición de la fuente. Para las condiciones de frontera se muestra que los parámetros del modelo de fronteras de tipo PML (Perfectly Matched Layer) también son dependientes de la frecuencia. Tanto para el parámetro de escala de la fuente como los parámetros que modelan la PML se provee valores y funciones que optimizan el rendimiento de dichas aproximaciones para problemas donde se requiere el modelamiento de la propagación de ondas en medios visco-acústicos. Inversión de onda completa (FWI) o Migración Reversa en Tiempo (RTM) son escenarios que dependen fuertemente del modelamiento apropiado del campo de onda y son campos potenciales de aplicación para los resultados que se presentan.

Palabras-clave: Atenuación sísmica; Modelamiento de propagación de onda; Medio visco-acústico Fuente sísmica; Condiciones de frontera PML

1. INTRODUCTION

Modeling wave propagation in real sub-surface is a complex task: the detailed propagation physics is complicated and several simplifying approximations are required for computing the solution for many practical situations. For example, it is usual to ignore the effects of attenuation in wave propagation; however in real life, the earth exhibits elastic and viscous behavior, properties that must be considered when one is interested in real and detailed physical processes associated, for example, with wave propagation in the subsurface for seismic exploration.

Ignoring any viscous or elastic behavior in the medium leads to the so-called acoustic wave propagation. Although unrealistic, in general, it is the simplest (but often practical) modeling approximation. Together with the assumption of isotropy, acoustic wave propagation has been widely used in different situations in seismic exploration [1]- [4]. However, to obtain accurate data on the physical properties of the sub-surface in realistic situations, the acoustic wave approximation does not provide any further information and it becomes necessary to do a better modeling of the propagation of waves. In particular, the realistic modeling of wave propagation is an important step in earth subsurface imaging process.

A step forward in the direction of doing realistic modeling of wave propagation is to consider the effects of viscosity. The visco-acoustic media can be defined as a medium without cross propagation, although it exhibits attenuation in the amplitude of the longitudinal wave. This medium presents two phenomena, dissipation and dispersion. The first one is produced by energy absorption such that the amplitude of the wave is reduced especially at high frequencies, while the second is induced by the change in the refractive properties of the media, where the wave velocity depends on the frequency [5].

To describe the attenuation of energy in the seismic wave front [6] proposed a model based on linear solid material rheology and memory variables. Later, [7] followed the same modeling approach, using, however, only one relaxation mechanism. They show how the method works to compensate for attenuation in least-squares reverse time migration. [8] used a visco-acoustic wave equation to compensate for the energy decrease of wave propagation in a realistic media, using an extrapolator based on the propagator of the wave equation in the forward and backward direction. Visco-acoustic modeling has been implemented also in 3D high computing demand simulations [9], [10] and has also been used in the development of specific geometrical configurations (TTI) [11], exhibiting accurate descriptions of the wave propagation in the media.

In general, visco-acoustic phenomena are modeled explicitly as frequency-dependent functions. The solution of the visco-acoustic wave equation in frequency domain must be performed numerically, and because of its simplicity and computational efficiency, the usual approach to solve it is to use finite difference methods. Nevertheless, the propagation problem must take in to account, besides the complex nature of the differential equation, crucial ingredients such as the source term and boundary conditions [12], [13].

In frequency domain, spatial resolution changes every time frequency changes, which complicates the characterization of different parts of the model such as the size of the zones of boundary and the size of the seismic source. All of the foregoing makes solutions of the visco-acoustic wave equation in frequency domain a problem far from trivial. In particular, since at a given frequency the structure of the grid should adapt to the scales associated with the frequency, the energy injection in the media could be radiated at different scales if not considering these effects carefully. The distribution of such energy among grid cells close to the source is a topic that is poorly discussed in the literature, especially the discretization of the source, which is commonly modeled as a Dirac delta function to describe its spatial part, [14]- [16]. This turns to be key in the solution of the numerical forward problem.

The same occurs with boundary conditions: Sponge absorbing boundary condition (ABS) Convolution Perfectly-matched layer (CPML), Perfectly-matched layer (PML), all having parameters that must be tuned, as can be observed in [18], [20] and [21].

When solving the problem of wave propagation in time domain, these are probably easy to fix, but in frequency domain, they are problematic, once again because the size of the region where absorbing boundaries work change with frequency and with the way attenuation is being modeled at the boundary. For example, in [17], these parameters are chosen based on trial and error. How to deal with these parameters? We present in this paper a study of how to deal with PML boundary conditions in frequency domain, showing in particular how to model, implement and control the behavior of the seismic source, and how to estimate the values of the parameters that control the behavior of this model. Furthermore, we show practical functions that allows the determination of the parameters required to model the performance of the PML boundary conditions.

For approaching this analysis, we used the following structure: First, in section 2 we present the basic theoretical considerations relative to our problem, presenting the basic equation used to model visco-acoustic wave propagation. Then, we present some considerations on the source term for seismic wave propagation, as well as considerations associated with the modeling of PML boundary conditions. Then, in section 3, we discuss the results obtained from the modeling of the source term and PML; finally, in section 4 we present our conclusions.

2. THEORETICAL FRAME

THE SOURCE

We consider the propagation of waves in 2D for the case of a visco-acoustic media. For this purpose, the equation for wave propagation in our visco-acoustic medium can be expressed as [16]

where P is the pressure field, K is the acoustic bulk modulus, ξ is the damping function [6] and b is the inverse of the medium density.

There are several choices on how to model attenuation. For example, using the models described in [5], which defines the wavenumber as K(ω)=ω/c(ω)=ω/v(ω) -ic(ω), where c(ω) is the complex velocity, v(ω) the phase velocity and k(ω) the attenuation wavenumber, one can rewrite things in terms of ξ like

Where y is the rate deformation function related to the viscosity of the medium that is often related to the Q attenuation factor, which in the simplest case (Kolsky model, see e.g. 5) is y/ω=1/2Q.

In Eq (1), the source term is

such that F(ω,x) is a body force. The source in these models is generally an impulse at point x s =(x s ,z s). The Ricker wavelet is a possibility to represent a seismic source, such that , F(ω,x)=R(ω)δ(x-x sJ where R(ω) is the Ricker wavelet given by

where δ(x-x s) is a Dirac delta function, but in the discretized domain can be represented with the Kronecker delta. The R 0 is the maximum amplitude of the ricker wavelet and f s is the principal frequency of the source. Now, assuming that the divergence of the force is , F(ω,x)=R(ω)δ(x-x s), the force could be expressed as

Thus, and using Eq. (1), we can write

Here, one can safely discard the second term to the right. When the velocity, density and attenuation are not constant, we propose for the source

Where (i s ,js) is a point in the grid where the source is placed and δi,isδj,js can be approximated as

This approximation comes from the Dirac delta function as the limit (in the sense of distributions) of the sequence of zero-centered normal distributions δσ(x)=1/(√πσ)exp-(x/σ)2 as →0. The factor σ must have a relationship with the grid spacing, i.e. σx= σz~ Δ.

The discretization of the Dirac delta is used because of the changes in grid spacing, which according to the condition Δ=λ/G r =c min /(fG r) depends on the frequency. It is important to make sure that the energy source is distributed evenly within the space to avoid losses by locating the entire amplitude at a point that cannot be located unequivocally in all grids. Therefore, it must meet the property of the Dirac Delta in the discrete case.

On such basis, it may be stated that the solution, is well defined in the discrete domain [14], [15]. Thus, we can define a=ax=az as, using Equation 7:

where s is a parameter that must be tuned to make sure that condition (9) is fulfilled, and at the end the parameter used to control the energy distribution of the source in the model.

PERFECTLY-MATCHED LAYER (PML) ABSORBING BOUNDARY CONDITIONS

Another important ingredient for the solution of Equation 1 is the boundary condition. In this work, we will study the use of perfectly-matched absorbing boundary conditions (PML). The absorbing boundary condition is a virtual boundary, very simple to use, as our media is dissipative in principle.

The PML method consists on using two damping functions to suppress the value of the pressure field in the boundaries at the edges of the square computational domain, yx(x) and y z (z) and a non-physical pressure wavefield P x and P z such that P=P x +P z . This technique is similar to a sponge-like absorbing boundary condition [18], but attenuation occurs at each dimension independently. The effect of these absorbing boundary conditions is that waves propagate in the medium, resembling the behavior of waves that propagate in an infinite medium (no reflections are produced at the borders of the computational domain), such that no need for any particular boundary conditions are required to specify the structure of the solution.

We simply expanded the computational domain ℧ with size characterized by N x xN z to the expanded domain ∂℧ with size N xe xN ze where N Xe =N x +2n xpm l and N ze =N z +2n zpml, where n xpml and n xpml are the number of extra points for the boundary condition (see Figure 1). In the expanded region, the damping functions have the form [17].

Figure 1 PML in Í1 and boundary workspace dil for the numerical solution, in x and z coordinates. 

Where y(ω,x,L z) is the rate deformation function related to the viscosity of the medium introduced a few equations ago, which we will assume as a constant parameter. Note that the functions y zpml (ω,x) and y xpmi (ω,x) are modified deformation functions that apply only in the region of the extended domain. However, the idea is that these deformation functions should be coupled with the deformation functions (and, therefore, to the attenuation) in the original field ℧. Here, m 0 is a parameter that should depend on the frequency and takes a value that makes the amplitude of the wave at the boundary of the domain to fall below a given threshold.

As to the corners, which generate well-known edge problems [17], we use the conventional treatment of averaging the boundary conditions at x and z [22].

Once again, for a given form of the damping function (ω,x,L z) , the parameter m 0 controls the behavior of the PML. We show how this parameter depends on frequency and find a way to fix his value in our implementation of the visco-acoustic modeler. It should be noted that the behavior of the PML also depends on the value of the number of pixels one considers to extend the domain, N xe ,N ze . If these numbers are too small, waves going through the PML region will not attenuate completely, thus producing artificial reflections on the computational domain. Furthermore, if one uses N xe ,N ze that are too large, the attenuation will take place, but there will be a waste of computational resources trying to solve the problem of wave propagation in regions that are of no interest. It would be good to then find a way to identify the smaller values of Nxe ,N ze to make sure that the behavior of the PML is the right one.

3. RESULTS

Now, we discuss the result of our exploration. First, we show how we use Equation 9 to study the effects of our model for implementing the seismic source and controlling its behavior through the use of parameter s. Then, we show how we control the behavior of the PML boundary through parameter m 0 .

Specific details on how we solve the differential equation using a mixed grid are out of the scope of this work, but r the reader may refer to [16] and [19] for more details on its implementation.

THE SOURCE TERM

To get a correct value for s, we calculate Eq. (9) for different frequency values, assuming the position of the source is x s=1Km and z s=1Km in a field of area of 2Km x 2Km, with c=2100 m/s. The frequencies we used are f=1,5,15,30,50Hz. We do not compute it for higher frequencies as the Δ in these cases is smaller, approaching the continuous form of the Dirac delta. With these values, we reach Figure 2.

Figure 2 (a) Integral values for discrete Dirac delta calculate for equation 9 for different frequency values, assuming the position of the source is Xs = 1 Km and Zs = 1 Km; (b) Placed at Xs = 1 Km and Zs = 0.015 Km. 

As we can see in Figure 2a, the values of l(s) remain close to 1 for s>0.9, except for the solution at f=lHz. However, in order to not being too restrictive, we notice that for scales in the range s=[0.65,1.6] l(s) keeps constrained in the range [0.95,1.05], It could be argued that it is enough to ensure the conservation of properties of the delta function around the source. For s (s<0.95) smaller value, the discretization affects the estimation of the spread of the source significantly.

In Figure 2b, we do the same calculation but with the source placed at x s=1Km and z s=15m. We can see in the figure that now the region of the interval of values of s that keep l(s,f)=[0.95,1.05] are reduced to the interval s=[0.6,0.75]. Note how the spread on l(s,f) increases, and how the value of l(s,f) is larger for smaller frequencies.

This experiment shows that the appropriated value of s shall depend on the frequency and the position of the source. Then, looking for a way to find evidence to argue the selection of the value of s, we computed Eq. 9 for different scale factors as a function of frequency.

Figure 3 shows the same situation as that of the previous figure, an area of 2Km x 2Km, with c=2100 m/s. The Figure 3a shows the result at position x s=1Km and z s=1Km. The Figure 3b shows the result at x s=1Km and z s=15m. The scale factors are s=0.65,0.75,0.85,1.0,1.5,2.0 for frequencies in the range f=[1.0,50.0] Hz.

Figure 3 (a) Integral values for discrete Dirac delta calculate for Eq. 9 for different scale factor, assuming the position of the source is x s = 1 Km and z s = 1 Km; (b) Placed at x s = 1 Km and z s = 0.015 Km. 

We can see in Figure 3 that the values of l(s,f) remain close to 1 for the scale factor below 1. We can also note that for low frequencies, the integral has a value close to 1 for scale factors s<l, for scaling factors s=[l,2] the integral has values distant to 1. For scale factor s=[0.65,0.75], l(s,f) fluctuates around 1, while for values in the range s=[0.75,1], we see an asymptotic convergence of l(s,f) to 1, for larger values of the frequencies. This behavior is much more desired as in a multlscale approach, this will ensure better performance of the source for larger frequencies.

It can be, therefore, assumed that the correct value for the scale factor, for our tests with sources located in the center of the domain or at the top, should be s=[0.75,1.0].

THE PML

Now, we focus our attention on the properties of the PML and the selection of the parameters m 0 and N xe ,N ze that control its behavior.

We compute the value of the P-wave amplitude of waves propagating in a 2Kmx2Km field. The media has a constant velocity of 2100m/s, Our Ricker source frequency is 30Hz. The constant quality factor is 0=50, using the damping function of Kolsky (see [5]). The solution uses a cell size Δ=λ/G r, where G r = 7, scale factor s=0.75. The position of the source is x 0= 1 Km and z 0 = 1 Km. We conducted several tests on the dispersion and stability of the solution. It was found that for the scales and frequencies of the experiments presented herein, using G r = 7 offers good solutions of the equation in terms of the Low dispersion. In addition, we tested several frequencies, yielding spectra that behave as expected, amplitude values whose highest value is in the main frequency (see [19]).

Figure 4 shows the real, imaginary parts, and modulus of field P-wave with PML (m 0 = f ) and Figure 5 shows the real, imaginary part, and modulus of field P-wave without PML (m 0 = 0) in a medium with constant velocity, attenuation and density for 15 Hz.

The expected result is a symmetrical and spherical field (in agreement with the result shown in Figure 4). However, in Figure 5 the absence of PML presents a perturbation that is not in agreement with the physical setup of the problem. There are some ripples observed in the figure, which are associated to interference of the incident and reflected waves. These waves are reflected at the borders when no PML is acting in the solution. In this figure, we can see the Importance of the PML so that the results are consistent with the physics of the medium we model (an open boundary problem).

In [17], m 0 is chosen based on trial and error; in this work, we propose defining m 0, n xpml and n zpml as a function of frequency. In Figure 6, we show the modulus of field P-wave (mc =3.0) for 3Hz. The figure at the top shows the full field, the figure at the middle shows an horizontal cut at z=1000m, and the figure at the bottom shows an horizontal cut at z=500m. Figure 7 shows the same, with m 0=3.0 for waves of frequency 10Hz.

In this case, we have used the same value of m 0=3.0 at different frequencies. Initially, for frequencies close to 3Hz, It seems that this choice of m0 was adequate; however, after several experiments, it was noticed that the more we increased the frequency, the PML performance started to deteriorate. This is because the term m0 controls the way in which the functions y pml attenuate the amplitude of the wave; the larger the frequency, the faster the attenuating term falls, so m0 must counteract this behavior at high frequencies. Therefore, m0 is frequency dependent.

We have noticed for low frequencies that m 0 =f Is not enough to achieve good response. Figure 8 shows a horizontal cut on the modulus of field P-wave at z=100m for a numerical solution using m 0 =f (solid line) and m 0= ω =2πf (dashed line) obtained with the 9-point schemes. For f=3,5,10,30Hz, It is found that the best choice for m 0 is m 0= 2 πf, especially at low frequencies. This choice makes sense considering the way attenuation Is modeled in the PML according to Equations 10 and 11, where we see that the explicit frequency dependence of the damping function Is cancelled.

In addition, n xpml and n zpml are also frequency dependent, which for low frequencies (long wavelengths) should have large values compared to the values of N x and N z such that the amplitude of the pressure field is attenuated correctly. On the other hand, for high frequencies (short wavelengths), only a few points n xpml and n zpml . are needed to damp the pressure field in the region of the PML. After a careful study whereby we systematically varied the frequency and the size of the PML, it was found that a convenient approximation for n xpml and n zpml as a function of frequency is

With c1 = -1.24x10-6, c2=3.37x10-4 , c3= -3.07x10-2 and c4=1.07. We found that using this approximation, the behavior of the PML produces the correct attenuation of the wavefield on the extended domain, in a way that no reflections or spurious information appear in the solution coming from the boundaries. What is most interesting in this approximation is that it provides an automatic way to determine the size of the PML region as a function of the frequency (as one would expect) without the need to find a value through trial and error tests.

Figure 4 Field P-wave with m 0 =f 

Figure 5 Field P-wave with m 0 =0. 

Figure 6 Modulus P-wave with m0=3.0 and f=3Hz. The first figure shows the modulus P-wave for all (x,z) and the following figures show the modulus P-wave for a fixed z. 

Figure 7 Modulus P-wave with m0=3.0 and f=10Hz. The first figure shows the modulus P-wave for all (x,z) and the following figures show the modulus P-wave for a fixed z. 

Figure 8 Field P-wave with m 0 =f, ω for f=3,5,10,30 Hz. 

CONCLUSIONS

In this work, we studied the propagation of waves in a visco-acoustic medium through explicit modeling of the attenuation, using damping functions that allow for dispersion depending on the quality factor. We have implemented a finite difference scheme to solve the problem in frequency domain. Special care was taken regarding the numerical dispersion issues of the modeling, for which we used a mixed grid technique and optimal setup of the intercalated grids to minimize numeric dispersion.

We focused our attention on two particular issues that are loosely discussed in the literature when modeling seismic wave propagation. The first one relates to the discreteness of the source function of the seismic source. The second point we discussed is the selection of the parameters that control the behavior of the PML.

Regarding the source, we show that the condition imposed to fix its behavior is given by the normalization condition of the Dirac delta pulse, and it is considered in the implementation presented in Equation 9, where the behavior of the source term is then controlled through a free parameter s that controls the spread of the source around neighboring cells.

We have shown that parameter s depends non-trivially on the frequency and the position of the source in the field. The dependence on the position of the source is mostly due to the effect of the boundary conditions of the problem; concerning the case of seismic exploration, it should represent a field of air or water boundary in one of the borders. It was found, in particular, that if the source is placed deep down from the border, where the upper boundary effects are neglectable, the parameter s is almost non-dependent on the frequency of the waves. A careful study of the behavior of the model source l(s,f) shows that appropriated values of parameter s that lead to an acceptable normalized behavior of the source term in the range of s=[0.75,1].

When we explored the implementation of the PML, it was noticed that this model ingredient depends on three parameters. The attenuation frequency m 0 , and the scale of the PML zone n_xpml and n zpml . We found that, contrary to what is done in other works (see for example [17]), it is not necessary to test by trial and error the values of m 0 that produce the best results for the performance of the PML. We show that m 0 =2 πf is a very much reasonable choice that also provides the best results out of the PML.

For the size of the PML, after a careful set of tests, a set of relations between n xpml and n zpml and the frequency of the waves was found. This relation facilitates the setup of modeling and inversion experiments intended to achieve optimal process performance.

ACKNOWLEDGEMENTS

This work is supported by ECOPETROL and COLCIENCIAS as a part of research project grant No. 0266-2013. The authors want to thank the referees for their observations and suggestions to the paper. JCMC thanks "Estrategia de Sostenibilidad, Universidad de Antioquia".

REFERENCES

[1] Pratt, R. G. (1999). Seismic waveform inversion in the frequency domain, part i: Theory andverification in a physical scale model. Geophysics, 64, 888-901. [ Links ]

[2] dos Santos A.W.G and Pestana, R.C. (2015). "Timedomain multiscale full-waveform inversion using the rapid expansion method and efficient step-length estimation." GEOPHYSICS, 80(4), R203-R216. [ Links ]

[3] Oliveira, A.L.D.; Pestana, R.; Santos, A. (2016). Least-Square Reverse Time Migration (LSRTM) In The Shot Domain. Revista Brasileira de Geofísica, [S.l.], v. 34, n. 3, p. 277-287. [ Links ]

[4] Moradpouri, F., Moradzadeh, A., Pestana, R., Ghaedrahmati, R. and Monfared M.S. (2017). "An improvement in wavefield extrapolation and imaging condition to suppress reverse time migration artifacts." GEOPHYSICS , 82(6), S403-S409. [ Links ]

[5] Wang, Y. (2009). Seismic Inverse Q Filtering. Wiley. [ Links ]

[6] Carcione, J. M., Herman, G. C., and Kroode, A. T. (2002). Seismic modeling. Geophysics, 67:1304-1325. [ Links ]

[7] Dutta, G., Lu, K., Wang, X. and Schuster, G. (2013). Attenuation compensation in least-squares reverse time migration using the visco-acoustic wave equation. In SEG Houston 2013 Annual Meeting. [ Links ]

[8] Huazhong, W., Libin, Z. and Zaitian, M. (2003). Seismic wave imaging in visco-acoustic media. Science in China Ser. A Mathematics, 47:146-154. [ Links ]

[9] Operto, S., Virieux, J., Amestoy, P., L'Excellent, J., Giraud, L. and Ben Hadj, H. (2007). 3D finite- difference frequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver: A feasibility study. Geophysics, 72:SM195-SM211. [ Links ]

[10] Brossier, R., Etienne, V., and Operto, S. and Virieux, J. (2010). Frequency-Domain Numerical Modelling of Visco-Acoustic Waves with Finite-Difference and Finite-Element Discontinuous Galerkin Methods. Acoustic Waves, Eds. Don Dissanayake. [ Links ]

[11] Operto, S., Virieux, J., Amestoy, P., L'Excellent, J., Giraud, L. and Ben Hadj, H. (2007). Finite- difference frequency-domain modeling of viscoacoustic wave propagation in 2d tilted transversely isotropic (tti) media. Geophysics , 74:T75-T95. [ Links ]

[12] Gholamy A. and Kreinovich,V. (2014) Why Ricker Wavelet Are Successful in Processing Seismic Data: Towards a Theoretical Explanation, Departmental Technical Reports, University of Texas at El Paso. [ Links ]

[13] Wang, Y., (2015). Frequencies of the Ricker wavelet, GEOPHYSICS , VOL. 80, NO. 2, P. A31-A37. [ Links ]

[14] Mayo, A. (1984). The fast solution of poisson's and the biharmonic equations on irregular regions. SIAM J. Sci ., 21:285-299. [ Links ]

[15] Tornberg, A. and Engquist, B. (2004). Numerical approximations of singular source terms in differential equations. Journal of Computational Physics. Volume 200, Issue 2, 1 November 2004, Pages 462-488. [ Links ]

[16] Avendaño, S. K. (2017). Full waveform inversion (FWI) in frequency Domain for the wave propagation for Visco-acoustic medium applied in synthetic data. Master thesis in physics, Universidad de Antioquia. [ Links ]

[17] Operto, S. and Virieux, J. (2006). Practical aspects of Frecuecny-domain finite-difference modellling of Seismic wave propagation. Ecole thématique SEISCOPE. [ Links ]

[18] Cerjan, C., Koslo, D., Koslo, R. and Reshef, M. A. (1985). Nonreflecting boundary condition for discrete acoustic and elastic wave equations. Geophysics , 50:2117-2131. [ Links ]

[19] Avendaño, S. K., Ospina M., Muñoz-Cuartas J.C. and Montegranario H. (2018). Modeling visco-acoustic wave equation in frequency domain using Mixed-Grid Finite-Difference method and attenuation-dispersion model defined for Quality factor. Submitted in Brazilian Journal of Geophysics (SBGf). [ Links ]

[20] Berenger, J. P. (1994). A perfectly matched layer for absorption of electromagnetic waves. Journal of Computational Physics, 114:185-200. [ Links ]

[21] Roden, J. A. and Gedney, S. D. (2000). Convolutional PML (CPML). An efficient FDTD implementation of the CFS-PML for arbitrary. Microwave and optical technology letters. 27:334-338. [ Links ]

[22] Johnson, S.G. (2007). Notes on Perfectly Matched Layers (PMLs), online at https://math.mit.edu/~stevenj/18.369/pml.pdf. [ Links ]

Received: April 24, 2018; Revised: July 31, 2018; Accepted: November 01, 2018

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License