Geophysical Curiosities

Geophysical Curiosity 1
When you look at the formula for the P-wave velocity, namely \( \sqrt[]{ \frac{K+4/3 \mu }{\rho}} \), with K the elastic modulus, which I rather call the stiffness or resistance to pure compression (only change in volume) and μ the resistance to deformation (distortion without a change in volume), you realize that the P-wave is not a pure compressional wave. Pure compressional waves occur in media which have a shear modulus μ=0, such as water and air.

Geophysical Curiosity 2
We know that seismic energy can be propagated in the earth in two ways: P-waves and S-waves. Now in an isotropic elastic medium the P-wave particle displacement is in-line with the propagation direction and in case of an S-wave the particle displacement is perpendicular to the propagation direction. This particle displacement can be in any direction as long as it is perpendicular to the propagation direction.

In anisotropic media, the situation is quite different.
For P-waves the particle displacement is deviating slightly from the propagation direction (ray direction), as the displacement is normal to the wave-front. The ray vector is namely no longer perpendicular to the wave-front.
For S-waves it is even more complicated. It was noticed in field experiments that the vertical S-wave reflection time from a horizontal interface differed depending on the polarization direction. Hence, the velocity did depend on the polarization. Therefore, it was concluded that two S-wave modes exist, one with the particle displacement parallel (S-fast), the other with the particle displacement perpendicular (S-slow) to fracture orientation (or fine layering). In addition, the S-wave can, and that is quite odd, only propagate with these two polarizations. So, if a source generates an S-wave with polarization in between, this S-wave will not propagate at all, but will split into the two modes discussed before and each mode will propagate with its own velocity. This is called “shear-wave splitting” or “birefringence”. And again, the polarizations are not exactly perpendicular to the propagation direction. Instead, they are perpendicular to the P-wave polarization, which as mentioned before, is neither in the propagation direction.

Geophysical Curiosity 3
Keep it simple, you stupid” This slogan is of course based on the well-known political slogan “It’s the economy, you stupid” and here refers to “why introducing anisotropy parameters” in the description of wave propagation, given that a full description of wave propagation is already complicated enough. The reason is that in general we don’t honour the real Earth when we try to model geophysical observations. In 1992, I gave an “oration” when I became Professor at the University of Utrecht, called “Do we do justice to the Earth” with a wink to my wife’s occupation as a judge. It must be clear, that the title refers to all the approximation we make to make geophysical problems tractable. If we want to model exactly the seismic response of real rocks, that means including small-scale inhomogeneities, we need complex mathematical algorithms and still can only do the modelling in detail for specific, well-described cases. And even then, you might wonder whether in practice we know or even can use such detail.

It is exactly for that reason that effective media are used. In these media, the detail is “summarised/replaced” by “simple” anisotropy parameters. It might appear “making it more complicated”, but after having familiarised yourself with it, you realize that it makes the description of wave propagation more realistic and especially more applicable. Nice examples are shown in the AVA exercise in the “Quantitative Reservoir Characterization” course.

Geophysical Curiosity 4
How confusing can you make it?” Take the description of anisotropy. The confusion is mainly due to using misleading symbols and different notations for the same quantities. Two examples: 1) using the symbol C for incompressibility and 2) the use of \( \delta \)(2) in one set of publications and \( \delta \)1 in other publications, both describing the same anisotropy parameter for propagation in the (x1,x3) plane. So, when reading publications make sure you are aware of the definitions.

In describing wave propagation, the elastic and anisotropy parameters occur in the so-called stress-strain relationship which needs to be combined with Newton’s first law of motion to obtain the wave equation. In principle (“worst case”: triclinic symmetry) there are 21 different (Cij) parameters needed to describe the medium. Luckily in practice we can get away with fewer parameters and in the “best case” with only 2. For almost all practical applications we have the following cases: orthorhombic (7 Cij parameters), Transverse Isotropic (5 Cij parameters) and Isotropic (2 Cij parameters). Orthorhombic describes a finely layered medium with vertical fractures, Transverse Isotropic describes either fine layering or a set of vertical fractures only, whereas Isotropic describes a medium without any sub-seismic scale layering nor fractures.

For weak anisotropy, we use the following parameters:
Orthorhombic: Κ, \( \mu \), \( \epsilon \)(1), \( \epsilon \)(2), \( \delta \)(1), \( \delta \)(2) , \( \delta \)(3), \( \gamma \)(1), \( \gamma \)(2),
Transverse Isotropy: VTI: (Κ, \( \mu \), \( \epsilon \), \( \delta \), \( \gamma \)), HTI: (K, \( \mu \), \( \epsilon \)(v), \( \delta \)(v), \( \gamma \)(v))
Isotropy: K, \( \mu \)

Unhappily, all rocks should be considered triclinic, but in reality, we can “ignore” that. We wouldn’t be able to derive the needed 21 parameters from seismic observations anyway. So, let’s forget it and restrict ourselves to at most orthorhombic or tilted orthorhombic in case of dipping structures.

One quantity we can observe quite well with modern wide- or multi-azimuth acquisition and that is the orientation of fracture sets. Namely, the NMO velocity of waves travelling perpendicular to fractures is lower than waves travelling parallel to fractures. In addition, if S-waves can be observed (multi-component land or OBC acquisition) shear-wave splitting also indicates the orientation of fractures.

Geophysical Curiosity 5
We are all aware of gravity, especially when you fall of your bike after too many drinks. Especially then it looks that gravity is a powerful force, although strictly speaking it is an acceleration. Interesting enough, gravity is the weakest force in nature. How can that be? You would agree when you realize that it needs the whole earth, that is a mass of 6 1024 kg, to make me feel it. Also, according to Newton, if I meet another person, I should feel the “attraction” between us, but I don’t. This shows how weak gravity is.

Geophysical Curiosity 6
In Seismic processing often, the word “Stacking” is mentioned. Traditionally, stacking referred to applying an NMO correction to offset traces and adding them to form a stack trace, with improved signal-to-noise ratio.
But do we still do stacking? Yes, we still do stacking, but in a different way. The stacking we do takes place in the imaging part of processing. Namely, in a modern processing scheme the data is prepared for imaging and after each offset is imaged, then the results are stacked. This can be most easily seen in case image gathers are built and displayed. After inspection and possibly re-imaging, these are the traces that are stacked.

Geophysical Curiosity 7
A P-wave reflected by an interface will change in amplitude according to the reflection coefficient. Hence, when the acoustic impedance (ρvp) increases across the boundary, a positive amplitude (meaning the wavefront starts with a compression) will remain a positive amplitude upon reflection.
But how about a S-wave? What will happen to the amplitude when reflected by an interface where the shear impedance (ρvs) increases?
The first thought would be that the reflection coefficient will be positive, so the amplitude will keep the same polarity.
But, no, no, that is not the case. How confusing. Why can that be?
The answer lies in the fact that we use in the industry a moving reference frame.
So, for P waves, the positive direction is downwards on the way down and upwards on the way up.
What does this mean for S waves, say SH? It means on the way down the positive shear direction is to the right and on the way up the positive shear is of course also to the right. But if the, say SH wave front start with a movement to the right and the shear impedance increases across the interface, on the way up the movement will be in the same direction. However, now that direction is negative due to the choice of the frame we are using. So, the reflection coefficient is negative for an increase in shear impedance.
I hope it is no longer confusing.

Geophysical Curiosity 8
The local angle domain. (LAD, Koren & Ravve, Geophysics, 2011, V76, N01).

When we consider fractures, we look at the source-receiver azimuth as an important parameter. Using either the NMO velocity as a function of azimuth (low resolution) or the AVA as a function of azimuth (high resolution) we determine the orientation of the fracture system(s). But is this correct?

It might be, but not always. Namely, what is important is the local azimuth of the “incident ray” and the “reflected ray” at the reflection point, because they are directly influenced by the local fracture orientation at both sides of the interface. So, in case the “local azimuth” is significantly different from the surface source-receiver azimuth, it is important to consider it.

But how often does this occur? I think seldom, but yes in cases of complicated geology and stress fields. Now those occur often in places where significant hydrocarbons can be found and where drilling is very expensive: deep water with reservoirs below salt or basalt.

So, there are cases when considering the local angle domain (LAD) is a must for obtaining a correct image of the subsurface.

Geophysical curiosity 9
I thought I knew a reasonable number of statistical distributions, until I went to Wikipedia ( and realised I only knew a few. (So, it is true that you don’t know what you don’t know).

More interesting, however, is that I recently came across a distribution called logistic distribution. What the name means I am lost, as logistics refers to moving goods from one place to another. However, I found it a very useful distribution for geophysics, as it can handle limited interval distributions of geophysical parameters of different magnitudes. Again, a description can be found in Wikipedia (

An interesting application can be found in an article by J. Eidsvik, D. Bhattacharjya and T. Mukerji: Value of information of seismic amplitude and CSEM resistivity. (Geophysics, 2008, Vol73, N4)

Geophysical curiosity 10
Systematically organised fractures or cracks play an important role in seismic wave propagation and in fluid flow of reservoirs. If of a size significantly below seismic wavelength their effect on seismic wave propagation can be described by anisotropy.

For a VTI medium we estimate the Thomsen weak anisotropy parameters \( \epsilon, \delta, \gamma, \eta \) from travel time information of P and PS or S wave.

For a HTI medium, caused by the presence of a single set of cracks, we estimate the Thomsen-type parameters \( \epsilon \)(v), \( \delta \)(v), \( \gamma \)(v), \( \eta \)(v) from either the azimuthal NMO velocities or AVA gradients of P and PS waves or from shear wave splitting (birefringence), in which case the fracture infill (brine, oil, gas, dry) has no influence. Note that for P or PS waves the infill has an influence.

There are several crack models, based on isolated infinite fractures with linear slip boundary conditions (Schoenberg), isolated penny-shaped cracks (Hudson) and hydraulically connected cracks and matrix pore space (Thomsen). Fortunately, seismic data are not sensitive to the shape of the fractures.

However, seismic is sensitive to pore fill & saturation and whether the pore spaces are hydraulically connected, which depends on the pore throat micro-geometry and the fluid properties. At seismic (low) frequencies, connection (pressure equilibrium) between fracture and medium pore space can be assumed. “Squirt flow” occurs at frequencies between sonic and ultrasonic bands. In that case Gassmann theory should be replaced by Biot theory. At high frequencies (ultrasonic) no fluid interaction occurs between fracture and matrix pore space. At very high frequencies Rayleigh scattering becomes important. Note that at intermediate frequencies, dispersion will occur due to the fluid flow between fracture and matrix pore space.

An empty crack is compliant because of its 2D shape, equant (non-flat) pores are stiffer. If a crack is isolated or connected to other cracks of similar shape and orientation than a pore fluid will stiffen the crack(s). If connected with equant matrix pore space, the fluid will “squirt” into the stiffer matrix pore space. Hence, the pores effect the compliance, that is the anisotropy.

There several ways of describing the presence of cracks or fractures. One way is to describe the influence of fractures on the wave propagation in terms of the parameters \( \epsilon \)(v), \( \delta \)(v), \( \gamma \)(v), \( \eta \)(v) (penny shaped cracks). The other way is to consider the fracture properties or excess fracture compliance parameters, thus the additional compliances due to fractures. Two sets are used. One uses KN (excess normal fracture compliance) and KT (excess tangential compliance). The other set uses similar, but dimensionless and normalised parameters, \( \Delta \)N (normal weakness) and \( \Delta \)T (tangential weakness) together with the Lamé parameters \( \lambda \) and \( \mu \) (linear slip theory). \( \Delta \)T gives estimate of the crack density and \( \Delta \)T /\( \Delta \)N is sensitive to fluid saturation. Obviously, the sets are related.

Crack density can be derived from the seismic observations of P, PS and/or AVA gradients (azimuths parallel/strike or normal to the fractures) or from shear wave splitting.

The main difference, except for shear waves, is between pore space filled with fluids (brine, oil) or filled with gas (dry). The dry-crack model is close to elliptical (\( \epsilon \)(v)= \( \delta \)(v)), for fluid-filled cracks \( \delta \)(v) is close to the crack density for typical background Vs/Vp ratios (around 0.5). In the (\( \epsilon \)(v), \( \delta \)(v)) plane all possible pairs of \( \epsilon \)(v), and \( \delta \)(v) are confined to a narrow area.

An important aspect is the frequency of the seismic wave. Low frequent (seismic) means that there is enough time for the pressure to equalize in fractures and matrix (equant) porosity, so the entire pore volume is involved. For moderately high frequencies (ultra-sonic) only the fracture volume is involved.

A complete characterization of fractures can be based on wide-azimuth 3D P-wave data or a combination of P and PS data or S data.

Geophysical curiosity 11
What happens with reflections beyond the critical angle?

Let’s look at AVA, that means the reflection coefficient with angle of incidence (offset). This can for a flat interface between two homogeneous layers be calculated exactly using the Zoeppritz equation, but also approximated quite accurately to first order in the velocity and density contrasts using an expression formulated in terms of horizontal slowness p (Aki & Richard, Quantitative Seismology, eq. 5.44):

\( R(p) = \frac{1}{2}(1 - 4 \beta^2p^2) \frac{ \Delta \rho }{ \rho } + \frac{1}{2(1 - \alpha^2p^2)} \frac{ \Delta \alpha }{ \alpha } - 4 \beta^2p^2 \frac{\Delta \beta}{\beta } \)

Aki & Richards mention that the linearized expression is only valid when none of the angles are near 90°. When the P-wave velocity below the interface is larger than above, the angle of transmission \( \Theta \)t can be 90° or “larger”.

However, the expression can also be given in terms of the average angle \( \Theta=( \Theta_1+ \Theta_t)/2 \)

\( R( \Theta)= \frac{1}{2} (1-4 \frac{ \beta^2 }{ \alpha^2 } sin^2 \Theta) \frac{ \Delta \rho }{ \rho } + \frac{1}{2cos^2 \Theta } \frac{ \Delta \alpha }{ \alpha } -4 \frac{ \beta^2 }{ \alpha^2 } sin^2 \Theta \frac{ \Delta \beta }{ \beta } \) ,

but now the transmitted angle can become complex:
\( \Theta_t= \frac{ \pi }{2} -i\ cosh^{-1} ( \frac{ \alpha_2 }{ \alpha_1 } sin \ \Theta_r) \) , \( \frac{ \alpha_2 }{ \alpha_1 } sin\ \Theta_r \geq1 \)

As a consequence, the reflection coefficient becomes complex. This means that the amplitudes (trace samples) become complex.

This makes us realize that for AVA, we calculate at an angle of incidence \( ( \Theta_i) \) the reflection amplitude for each sample. That means looking at a sequence of samples, forming a wavelet, the same wavelet, but scaled with the reflection coefficient, will be returned. But only in case the reflection coefficient is real. If complex, all samples will have the same phase change and the reflected wavelet will have a different appearance. As the phase change is angle (offset) dependent, the wavelet change will be angle (offset) dependent, complicating the analysis somewhat.

Ref: Downton & Ursenbach, Geophysics, 2006, V71, N4.

Geophysical curiosity 12
What is the difference between updating a velocity model in the data or in the image domain? Are we using the same data?

Updating of velocities means that we have an estimated velocity-depth model and look for diagnostics and ways which will allow improving the velocities.

Data domain
We have observed the arrival times at all receivers for all shots. We know the shot and receiver locations at the surface and a rough/first estimate of the velocity model, which will be represented by a grid of cells, with for each cell \( (i) \) a velocity \( (v_i) \). Then we apply raytracing from each shot to the receivers. That means that the ray will go through a subset of all cells. For each ray we can calculate the arrival time by adding all the \( \Delta t_i \) calculated from the pathlength \( (l_i) \) in each cell divided by the cell velocity \( v_i \). We will find a difference between the calculated/modelled and the observed travel time. This difference could come from anywhere along the ray. So, we simply divide the difference over the total ray path, leading to a small velocity update in each cell traversed. Having many shot-receiver combinations, means that (hopefully) each cell will be traversed by many rays. Each traverse will give a small contribution to the update. If they agree the sum will be a significant update, if they disagree the sum will be “zero”, that means no update. That’s all there is to it in the data domain. So we don’t do any migration/imaging.

Image domain
We migrate the recorded data using the first estimate of the velocity-depth model. The results are Common Image Gathers (CIG). They will show residual move-out (diagnostic), which will be used to update the velocity model.

Although in both cases we use the same data, the advantage of using the image domain is that the signal-to-noise ratio will usually be much better. The disadvantage is that it needs more computer time.

Geophysical curiosity 13

Question: Is seismic a hologram?

An optical hologram is not an image; it is a recording of a wavefield. Remember how a hologram is made: A single frequency laser light is split into two parts. One part is mirrored directly onto the recording plate; the second one is reflected / diffracted by the object. What is recorded on the plate is the interference pattern between the two.
When we beam the laser onto the “developed” plate and look through we see the object in 3D. Looking from different directions we see the object from these different directions. Even more fascinating is that when we break off a part of the plate, we can still see the total object from all directions, although the smaller the remaining piece the poorer the quality of the image.


Remember that in seismic we can consider the subsurface to consist of a collection of diffraction or scattering points. This would imply that a single shot recorded by a single geophone would record diffraction energy from any point of the subsurface. Hence, the single record contains information about the entire subsurface and it would be possible, in principle to reconstruct/image the complete subsurface. This is like the concept that a single point on a hologram is enough to observe the object in 3D. But again, we would prefer a larger piece of the plate, i.e. more observations, to obtain a better-quality image.

Lesson learned
If we would apply holographic seismic imaging, we could do with a much sparser and therefore cheaper survey.

“Extended resolution: Neidell is right”, E.A. Robinson, TLE 2018V37, N1
‘holistic migration”, E.A. Robinson, TLE 1998, V17, N3

Added statements by Öz Yilmaz
“Spatial under sampling can be compensated by temporal oversampling”
“Temporal oversampling enables the extension of the spectral bandwidth beyond the Nyquist frequency”.

Geophysical curiosity 14

Migration or Wavefield reconstruction versus Imaging

It is important to make a distinction between wavefield reconstruction and imaging.

Migration is a two-step process: 1) wavefield reconstruction and 2) imaging.

Wavefield reconstruction is the determination of the wavefield over a volume of interest. Imaging involves the making of a picture of the geometrical distribution of the reflecting surfaces within the medium.

This is best illustrated by the fact that we perform a downward continuation, forward in time, of the wavefield generated by a source and a downward continuation, but now backward in time, of the recorded wavefield, followed by the application of an imaging condition to obtain the reflection point in space and its reflection strength.

Again, we have various methods for wavefield continuation, the more sophisticated the more accurate (acoustic, elastic, visco-elastic) and various so-called imaging conditions (cross-correlation, deconvolution, least-squares deconvolution).

Geophysical curiosity 15

Up to now I didn’t know of the existence of quaternions and that they can also be used in geophysics. What is a quaternion?
A complex number is defined as a+ib, a quaternion as q=a+bi+cj+dk, hence it is a hypercomplex number.
It allows a Fourier domain representation of multi-component seismic data. For example, so-called “Projection Onto Convex Sets” (POCS) reconstructs all three seismic components at once.
Applications are in image disparity methods (edge detection), interpolation, time-lapse analysis, separation of up- and down-going waves (VSP) and calculation of attributes.
However, so far I have encountered very few “real” applications.

Please feel free to comment or add points of view and indicate whether you would like it to be mentioned on the website or not

Back to main page

© Breakaway

Last modified: Friday, 17 August 2018, 12:34 PM