Viscoelastic seismic modeling and Q estimation for an

attenuating media

Nimisha Vedanti1,2 , Ajay Malkoti2

1
CSIR-National Geophysical Research Institute, Hyderabad, India-500007
2
Academy of Scientific and Innovative Research-NGRI, Hyderabad, India-500007
[email protected]

Abstract

In seismic modelling experiments, the propagating seismic wave experiences attenuation

when real media is simulated. However, for computational ease we often ignore this effect

and model the Earth as either an acoustic or an elastic medium. In this paper we briefly

discuss the importance of considering attenuation in seismic modelling by using a

viscoelastic wave equation. Further we also briefly discuss the effects of attenuation on

seismic data and its incorporation in seismic modelling. We see that attenuation has a large

effect on synthetic seismogram which should be properly addressed in advanced seismic

processing and imaging methods.

1 | Page
Introduction

Attenuation is a property of a medium due to which the energy of a seismic wave can be

dissipated in the form of heat and thus leading to a reduction in the wave amplitude. Factors and

phenomena responsible for the attenuation can be broadly classified as either extrinsic, or

intrinsic. In extrinsic attenuation, redistribution of energy leads to a reduction of amplitudes,

while in intrinsic attenuation the waves suffer energy loss due to conversion of mechanical

energy into heat. Commonly observed wave phenomena like geometric spreading, scattering,

leaky modes, etc. can be cited as examples of extrinsic attenuation and they do not contribute

towards intrinsic attenuation. In general, seismic attenuation is quantified in terms of quality

factor ‘Q’, which can be defined as a ratio of stored energy to dispersed energy, as it measures

relative energy loss per oscillation cycle of seismic waves. In ideal scenario, ‘Q’ is related to the

physical state of the rock and it increases with increase in density and velocity.

Since attenuation results in a loss of seismic energy, the recorded seismic traces are

required to be compensated for an exact of amount of “Q”. To achieve this, determination of a

precise ‘Q model’ is required, especially in the high attenuating regions like Deccan Volcanic

Province (DVP), where thick basalt sequences cause high attenuation to the propagating seismic

waves (Vedanti et al., 2018, 2015). In these areas, data processing using conventional methods

may fail because of low signal to noise ratio especially in the sub-basalt formations. In a

synthetic study, carried out by Malkoti et al., (2015), it is shown that due to attenuation,

amplitudes of the late arrivals or amplitudes from deeper reflectors are highly diminished in the

generated seismogram, which makes processing and interpretation very challenging in sub-basalt

layers. In recent past, attempts have been made to apply advanced techniques like FWI to
improve the seismic imaging. This technique needs ‘complete wave field information’ with

precise amplitudes to obtain the accurate gradients and thus, it needs a ‘Q’ structure of the

domain. There are several techniques mentioned in the literature for estimation of ‘Q’ from the

acquired seismic data. Most of these techniques are based upon certain data attributes, such as

spectral ratio technique depends mainly on the amplitude, however, the seismic attenuation has a

high influence on several other attributes known as the first order effect. Hence, considering the

medium as an attenuating media has its own consequences on ‘Q’ estimation techniques and it is

quite likely to obtain different values of ‘Q’ from the same data using different methods (Tonn,

1991). To understand this problem of the estimation of “Q”, we first need to understand the

theory of viscoelasticity and the aspects we should consider while incorporating attenuation in a

seismic modelling experiment. Thus, in this paper, we briefly discuss these aspects of

incorporating attenuation in seismic modelling and demonstrate its need by using a simple Earth

model.

Viscoelastic Models

Viscoelasticity is a property of a medium that exhibits both viscous and elastic

characteristics when it undergoes deformation. In a pure elastic medium, the process of

‘application and removal of the load’ follows energy conservation; however in a viscoelastic

material, it involves energy loss. The energy is dissipated during the loading and unloading cycle

and thus obtained hysteresis curve area can be used to estimate the attenuation or the quality

factor ‘Q’. There are several approaches available in literature to model the quality factor, which

includes simple damping, frictional models, complex moduli, time dependent moduli etc.

(Carcione, 2007; Kjartansson, 1979; Liu et al., 1976; Tal‐Ezer et al., 1990). However, in this

3 | Page
paper, we preferred to follow more reasonable ‘time dependent moduli’ approach for which it is

required to define the basic characteristics of a viscous material.

Ideal viscoelastic material exhibits the characteristics functions (creep and stress-relaxation

functions) as shown in Fig. 1. These characteristic functions are obtained under the two kinds of

tests, where the first test is, “Creep Recovery Test”, and the second is “Stress Relaxation Test”.

In former, a constant stress is applied on the sample for a fixed duration of time and created

stress changes are observed over the body. In the latter test, a constant strain/deformation is

applied on the body and then the stress variation with time is observed.

Figure 1: Characteristic functions for an ideal viscoelastic material

Different mechanical models can be defined to describe such a viscoelastic material. In these

models, number of parameters may vary from 1 to 𝑁 depending upon different configurations,

which can be used to approximate the rheology of the material by satisfying the above defined

tests functions. While doing so, we try to mimic the attenuation by using minimum number of

elements to incorporate the correct amount of attenuation, while preserving the characteristics.

Thus, we start with a simplest family of models i.e., the 1-element/parameter family that includes
single springs or dash pots. In a pure elastic medium there is no dissipation of energy, so it can

be represented with a spring. However, a viscous material can lose the energy and hence it is

modeled as a dashpot, which acts like a damper. As mentioned earlier, springs represent the

elastic and dash-pots represent the viscous components. Stress-strain relation for these two can

be written as

For spring, σ = kϵ
dϵ
For Dashpot, σ =η
dt

To represent the Earth more accurately, these two models shall be used in different combinations

to form models with higher number of parameters. 1 spring and 1 dashpot model is called as 2

parameters family. The spring and dashpot components can be arranged either in parallel or in

series which are called as Kelvin-Voigt model and Maxwell model, respectively. These

arrangements should properly model the characteristics shown in Figure 1. However, the Kelvin-

Voigt model cannot model either creep function or stress-relaxation function correctly, which is

a disadvantage while Maxwell model only has acceptable stress relaxation function. This

limitation led to the inclusion of higher number of parameters such as 3-parameter model e.g.,

Zener model (Fig. 2), also known as Standard Linear Solid (SLS) model, and 4-parameters

model e.g., Burger model.

Figure 2: The Zener model with three elements (wiki).

5 | Page
The most useful model to compute the time dependent moduli among these models is a

Generalized Standard Linear Solid (GSLS) model for which the general stress-strain relation can

be written as:

L L
∂i σ(t) ∂i ϵ(t)
a0 σ(t) + ∑ ai = b 0 ϵ(t) + ∑ b i (1)
∂t i ∂t i
i=1 i=1

Where, ai , bi are constants for a linear material. The complex modulus (𝑀∗ ) for a general

viscoelastic material can be calculated by applying the sinusoidal varying oscillating stress of

given frequency, e.g., σ(ω) = σ0 exp(iωt). After certain amount of time the initial effect is

negligible and the strain will also be in the form of ϵ(ω) = ϵ0 exp(iωt). The complex modulus of

the material is given by (Bland, 1960; Christensen, 2012)

σ0 (ω) a0 + ∑L=1 ai ωi
M(ω) = 0 = (2)
ϵ (ω) b0 + ∑Li=1 bi ωi

Thus, the better approach to model time dependent moduli is to assemble many Zener models in

parallel, which gives us a “Generalized Zener model” (Fig. 3).

Figure 3: Standard Linear Solid model with Zener elements

Viscoelastic wave formulation

There are many formulations available in literature to define the above mentioned viscoelastic

materials in time domain. In this paper, we only discuss the memory variable approach which is

based upon a SLS model (Carcione et al., 1988; Liu et al., 1976). As we know that the stress (σij )

strain (ϵkl ) relationship for a viscoelastic medium can be written in the form of convolution

integral as

σij = Λ̇ (t) ⋆ εij (t)δij + Ṁ ⋆ εij (3)

Where Λ and M are the relaxation functions. A general relaxation function (𝛹) can be modeled

with the help of SLS model comprised of L relaxation mechanisms (Carcione, 2007).

L
1 τlε t
ψ(t) = MR [1 − ∑ (1 − l exp (− l ))] H(t − t′) (4)
L τσ τσ
l=1

MR is the relaxed modulus, H(t) is the heavy side function, τlσ and τlε stress and strain relaxation

times, respectively. The stress-strain and memory variable equations can be obtained by

substituting the relaxation function (Eq. 4) into Eq. 3. Further, separate relaxation functions

should be considered for the P-wave and S-wave. Thus, it will yield following relations:

L
∂σij 1
= Λ ε̇k,k δij + 2M ε̇i,j − ∑ rijl
∂t L (5)
l=1

∂rijl 1 l 1 1
=− l
rij + l Λl ε̇k,k δij + l M l (ε̇i,j + ε̇j,i ) (6)
∂t τσ τσ τσ

7 | Page
Where, we have Λ = Π − 2M; Λl = Πl − 2M l , and Λ R = ΠR − 2MR . To simplify the equations

we have assumed that, Π = ΠR (1 − ∑Ll=1 Tpl ); M = MR (1 − ∑Ll=1 Tsl ); Πl = ΠR Tpl ; and M l =

MR Tsl . Here ΠR and 2MR are relaxed modulus for respective waves functions, rijl , is known as

1 τlεp 1 τlεs
the memory variable, and Tpl and Tsl stands for L (1 − ) and L (1 − ), respectively.
τlσ τlσ

Eq. (5) and (6), along with the continuity equation, forms a complete set of equation for

viscoelastic modelling. Thus the complete set of wave equation for 2D viscoelastic wave can be

written in expanded form as follows.

∂vx ∂σxx ∂σxz

ρ = + + ρfx (7)
∂t ∂x ∂z

∂vz ∂σzx ∂σzz

ρ = + + ρfz (8)
∂t ∂x ∂z
L
∂σxx 1 l
= (Λ + 2M) ε̇x,x + Λ ε̇z,z − ∑ rxx (9)
∂t L
l=1

L
∂σzz 1 l
= Λ ε̇x,x + Λ ε̇y,y + (Λ + 2M) ε̇z,z − ∑ rzz (10)
∂t L
l=1

L
∂σxz 1 l
= M ε̇x,z + M ε̇z,x − ∑ rxz (11)
∂t L
l=1

l
∂rxx 1 l 1 1
= − l rxx + l (Λl + 2M l ) ε̇x,x + l Λl ε̇z,z (12)
∂t τσ τσ τσ

l
∂rzz 1 l 1 1
= − l rxx + l Λl ε̇x,x + l (Λl + 2M l ) ε̇z,z (13)
∂t τσ τσ τσ

l
∂rxz 1 l 1
= − l rxz + l M(ε̇x,z + ε̇z,x ) (14)
∂t τσ τσ
The relaxation times can be determined by minimizing the error between Q(ω) and given Q0

(Blanch et al., 1995):

ω2
2
ϕ = ∫ [Q−1 (ω, τlσ , τlε ) − Q−1
0 (ω)] dω (15)
ω1

ℜe[Mc (ω)] 1
where, Q(ω) = ℑm[Mc(ω)] ; M C (ω) =ℱ { ∂t [ψ(t)] } and τlσ = ω is the stress relaxation time.
l

These equations can be solved using the finite difference method. Here we use a synthetic data

set to demonstrate the importance of the method.

Measurement of attenuation

In general attenuation measurements are carried out using the amplitude information from the

seismic data. Thus, in this paper we only discuss the amplitude based spectral ratio technique, for

‘Q’ estimation. The data used in this study was generated by designing a synthetic Vertical

Seismic Profiling (VSP) survey. The VSP geometry used is shown in Fig (4), where we have laid

the receivers vertically and placed the source at the top for the seismic wave simulation.

Spectral ratio technique

This method is very common among the geophysicists and has many variants. The

fundamental principle of this method is to compare the spectral characteristics of the

seismogram while assuming that the amplitude of a wave can be described by the following

relationship:

𝜔𝑡
𝐴(𝜔) = 𝐺(𝑡)𝑅(𝑡)𝐴0 (𝜔) exp (− ) (16)
2𝑄

where, 𝑨𝟎 is the initial amplitude of the wave which has reduced to 𝑨 after travelling for

time t in the given medium of attenuation characterized by 𝑸. 𝑮 and 𝑹 represent the

9 | Page
 reduction in amplitude due to geometric spreading and reflectivity, respectively. Assuming

two receivers placed some distance apart records the amplitude 𝑨𝟏 (𝝎) and 𝑨𝟐 (𝝎)

respectively at time 𝒕𝟏 and 𝒕𝟐 respectively. We obtain

𝐴2(𝜔) 𝜔Δ𝑡
ln ( )= − (17)
𝐴1 (𝜔) 2𝑄

Where, 𝚫𝒕 = 𝒕𝟐 − 𝒕𝟏 . When we plot the logarithm of amplitude ratio with the frequency, it

𝚫𝐭 𝝅
represents the equation of straight line with a slope as equal to , which can be in turn
𝑸

utilized to estimate the attenuation 𝑸.

Figure 4: Arrangement of source (in red asterisk) and receiver (blue triangles) along with

absorbing boundaries (along edges) as used in simulation.

First order effect of attenuations

When a medium offers attenuation to seismic waves, certain phenomena come into play which

can be understood as the first order effect of attenuation. Here we describe the most important

ones and their mitigation.

1) Frequency dependence of Q: The attenuation experienced by a seismic wave is dependent

on the frequency and attenuation mechanism. In a seismic modelling experiment, it’s

advised to assume constant attenuation over the seismic frequency range (McDonal et al.,

1958).

2) Frequency dependence of the reflectivity: Attenuation causes the reflectivity of the

interfaces to become frequency dependent and thus can affect the attenuation estimation.

Hence one must constrain the experiments to include only/nearly normal incident rays. It

can be achieved using VSP geometry.

3) Velocity dispersion causing the travel time difference (drift): As mentioned earlier,

attenuation causes different frequencies to be attenuated differently. The higher

frequencies attenuate faster in comparison to lower frequencies. This causes changes in

seismic wavelet and it gets broader with depth. It also causes shifting of the peak

amplitude and thus the events experiences a drift. This effect can be corrected or

compensated by providing the appropriate time shift or by matching with well log.

Numerical Simulation of Seismic Wave Propagation in Viscoelastic media

Following the above mention concepts, we carried out a synthetic VSP modelling for the elastic

as well as viscoelastic medium. The outputs of these simulations are compared in Fig. 5. Some

11 | Page
of the important parameters for the model and for the simulation are shown in the Table 1, along

with their corresponding values. It can be seen that the receivers are arranged vertically as in

VSP geometry (Fig. 4). The simulation parameters, i.e. time step (𝑑𝑡) was taken according to

stability condition and to minimize the grid dispersion, more than 6 grids/wavelength were used.

To suppress the edge reflections, damping type absorbing boundaries of (Cerjan et al., 1985)

were applied on all the sides. The synthetic VSP seismogram generated by using viscoelastic

formulation, considering attenuation in the media, is used for ‘Q’ estimation. The results

generated after taking care of all the above mentioned effects are shown in Fig. 5.

We have also used a two layered model and generated a VSP record for this model. Source and

simulation parameters were same as for homogeneous. We have chosen two nearby locations to

picked events (down going and upcoming) and computed the attenuation for down going as well

as upcoming wave. Fig. 6 demonstrates the key steps involved in attenuation measurement using

the Spectral Ratio method for a two layer case. The steps are namely- selection of traces,

clipping the waveform part, Fourier Transform of the clipped part, and estimation of attenuation

using the slope of linear fit. For the seismic wave simulation, we have used the FDwave package

developed by (Malkoti et al., 2018a, 2018b).

Table 1: List of parameters used in seismic simulation carried out for a viscoelastic media.

Parameter Value

Model physical parameters

Velocity, 𝑽𝒑 2000

Velocity, 𝑽𝒔 1700

Density, 𝝆 1900

Quality factor, Q 70

Source Parmeters

Source signature Ricker

Central frequency, 𝒇𝟎 15Hz

Zero time offset/shift, 𝒇𝟎 0.07sec

Total time length, 𝑻 1sec

Simulation parameters

Model size (x,z) 3km x 3km

Grid spacing, 𝚫𝒉 5m

Time step, 𝚫𝒕 0.1msec

Absorbing boundary nodes 40

13 | Page
Figure 5: A zoomed section of the VSP gather generated using elastic wave formulation (in red

color) and viscoelastic formulations (in blue color).

Results

The Fig. 5 shows a very simple experiment using homogeneous model to demonstrate the

difference between the elastic and viscoelastic seismogram. The difference due to the first order

effects is very prominent for late phases. This type of mismatch, if not taken care of, can lead to

the erroneous results. In Fig. 6, we have shown a successful application of this method for Q

estimation. As we can see in this figure that considering the first order effects of attenuation and

following the prescribed solutions, we can estimate the value of attenuation quite precisely.
Figure 6: An example of attenuation estimation using spectral ratio. The complete traces selected

for the estimation (a); The clipped part used for the estimation at two places (b); Fourier

spectrum of the given clipped traces (c), The least square fitting for the determination of slope to

estimate Q (d). Here ‘Trace 1 refers’ to elastic VSP trace and Trace 2 refers to viscoelastic VSP

trace.

15 | Page
Conclusions

In this paper we have discussed the behavior of an attenuating media and how to model it, using

the theory of viscoelasticity, we have shown that the seismic attenuation has a large effect on

synthetic seismogram, which was generated for a synthetic VSP survey. Further, we

demonstrate an application of Spectral ratio technique to estimate precise value of ‘Q’ using the

synthetic seismic data. We have also discussed the first order effects of attenuation, which should

be considered while incorporating seismic attenuation in the seismic modelling experiment.

Further, details on the theory and seismic wave simulation algorithm being used are available in

(Carcione et al., 1988) and (Malkoti et al., 2018a, 2018b) respectively.

References

Blanch, J., Robertsson, J., Symes, W., 1995. Modeling of a constant Q: Methodology and
algorithm for an efficient and optimally inexpensive viscoelastic technique. Geophysics
60, 176–184. https://doi.org/10.1190/1.1443744

Carcione, J.M., 2007. Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic,
Porous and Electromagnetic Media. Elsevier.

Carcione, J.M., Kosloff, D., Kosloff, R., 1988. Wave propagation simulation in a linear
viscoelastic medium. Geophysical Journal International 95, 597–611.

Cerjan, C., Kosloff, D., Kosloff, R., Reshef, M., 1985. A nonreflecting boundary condition for
discrete acoustic and elastic wave equations. Geophysics 50, 705–708.
https://doi.org/10.1190/1.1441945

Kjartansson, E., 1979. Constant Q-wave propagation and attenuation. Journal of Geophysical
Research: Solid Earth 84, 4737–4748. https://doi.org/10.1029/JB084iB09p04737

Liu, H.-P., Anderson, D.L., Kanamori, H., 1976. Velocity dispersion due to anelasticity;
implications for seismology and mantle composition. Geophysical Journal International
47, 41–58.

Malkoti, A., Vedanti, N., Kunagu, P., Tiwari, R., 2015. Modeling viscoelastic seismic wave
propagation in Deccan flood basalt, western India, in: SEG Technical Program Expanded
Abstracts 2015, SEG Technical Program Expanded Abstracts. Society of Exploration
Geophysicists, pp. 3764–3768. https://doi.org/10.1190/segam2015-5898555.1

Malkoti, A., Vedanti, N., Tiwari, R.K., 2018a. An algorithm for fast elastic wave simulation
using a vectorized finite difference operator. Computers & Geosciences 116, 23–31.

Malkoti, A., Vedanti, N., Tiwari, R.K., 2018b. Viscoelastic modeling with a vectorized finite
difference operator. Computers & Geosciences (Revised submitted).

McDonal, F., Angona, F., Mills, R., Sengbush, R., Van Nostrand, R., White, J., 1958.
Attenuation of shear and compressional waves in pierre shale. GEOPHYSICS 23, 421–
439. https://doi.org/10.1190/1.1438489

Tal‐Ezer, H., Carcione, J., Kosloff, D., 1990. An accurate and efficient scheme for wave
propagation in linear viscoelastic media. Geophysics 55, 1366–1379.
https://doi.org/10.1190/1.1442784

Tonn, R., 1991. The determination of the seismic quality factor Q from VSP data: a comparison
of different computational methods. Geophysical Prospecting 39, 1–27.

17 | Page
Vedanti, N., Lakshmi, K., Dutta, S., Malkoti, A., Pandey, O.P., 2015. Investigation of
Petrophysical Properties and Ultrasonic P-and S-Wave attenuation in Deccan Flood
Basalts, India, in: 2015 SEG Annual Meeting. Society of Exploration Geophysicists.

Vedanti, N., Malkoti, A., Pandey, O.P., Shrivastava, J.P., 2018. Ultrasonic P- and S-Wave
Attenuation and Petrophysical Properties of Deccan Flood Basalts, India, as Revealed by
Borehole Studies. Pure and Applied Geophysics 175, 2905–2930.
https://doi.org/10.1007/s00024-018-1817-x