Finite difference solution of the 2D wave equation (`fatiando.seismic.wavefd`)¶

Finite difference solution of the 2D wave equation for isotropic media.

Warning

This code is experimental and poorly tested! We cannot guarantee that the results are accurate. We strongly discourage use of this code for research purposes. Still, the general behaviour of the waves seems consistent with theory and can be useful for illustrating how seismic waves behave (reflection and refraction, for example).

elastic_psv: Simulates the coupled P and SV elastic waves using the Parsimonious Staggered Grid method of Luo and Schuster (1990)
elastic_sh: Simulates SH elastic waves using the Equivalent Staggered Grid method of Di Bartolo et al. (2012)
scalar: Simulates scalar waves using simple explicit finite differences scheme

Sources

MexHatSource: Mexican hat wavelet source
SinSqrSource: Sine squared source
GaussSource: Gauss derivative source
blast_source: A source blasting in all directions

Auxiliary functions

lame_lamb: Calculate the lambda Lame parameter
lame_mu: Calculate the mu Lame parameter
xz2ps: Convert x and z displacements to representations of P and S waves
maxdt: Calculate the maximum time step for elastic wave simulations
scalar_maxdt: Calculate the maximum time step for a scalar wave simulation

Theory

A good place to start is the equation of motion for elastic isotropic materials

\[\partial_j \tau_{ij} - \rho \partial_t^2 u_i = -f_i\]

where $\tau_{ij}$ is the stress tensor, $\rho$ the density, $u_i$ the displacement (particle motion) and $f_i$ is the source term. But what I’m interested in modeling are the displacements in x, y and z. They are what is recorded by the seismometers. Luckily, I can use Hooke’s law to write the stress tensor as a function of the displacements

\[\tau_{ij} = \lambda\delta_{ij}\partial_k u_k + \mu(\partial_i u_j + \partial_j u_i)\]

where $\lambda$ and $\mu$ are the Lame parameters and $\delta_{ij}$ is the Kronecker delta. Just as a reminder, in tensor notation, $\partial_k u_k$ is the divergent of $\mathbf{u}$. Free indices (not $i,j$) represent a summation.

In a 2D medium, there is no variation in one of the directions. So I’ll choose the y direction for this purpose, so all y-derivatives cancel out. Looking at the second component of the equation of motion

\[\partial_x\tau_{xy} + \underbrace{\partial_y\tau_{yy}}_{0} + \partial_z\tau_{yz} - \rho \partial_t^2 u_y = -f_y\]

Substituting the stress formula in the above equation yields

\[\partial_x\mu\partial_x u_y + \partial_z\mu\partial_z u_y - \rho\partial_t^2 u_y = -f_y\]

which is the wave equation for horizontally polarized S waves, i.e. SH waves. This equation is solved here using the Equivalent Staggered Grid (ESG) method of Di Bartolo et al. (2012). This method was developed for acoustic (pressure) waves but can be applied without modification to SH waves.

Canceling the y derivatives in the first and third components of the equation of motion

\[\partial_x\tau_{xx} + \partial_z\tau_{xz} - \rho \partial_t^2 u_x = -f_x\]

\[\partial_x\tau_{xz} + \partial_z\tau_{zz} - \rho \partial_t^2 u_z = -f_z\]

And the corresponding stress components are

\[\tau_{xx} = (\lambda + 2\mu)\partial_x u_x + \lambda\partial_z u_z\]

\[\tau_{zz} = (\lambda + 2\mu)\partial_z u_z + \lambda\partial_x u_x\]

\[\tau_{xz} = \mu( \partial_x u_z + \partial_z u_x)\]

This means that the displacements in x and z are coupled and must be solved together. This equation describes the motion of pressure (P) and vertically polarized S waves (SV). The method used here to solve these equations is the Parsimonious Staggered Grid (PSG) of Luo and Schuster (1990).

References

Di Bartolo, L., C. Dors, and W. J. Mansur (2012), A new family of finite-difference schemes to solve the heterogeneous acoustic wave equation, Geophysics, 77(5), T187-T199, doi:10.1190/geo2011-0345.1.

Luo, Y., and G. Schuster (1990), Parsimonious staggered grid finite-differencing of the wave equation, Geophysical Research Letters, 17(2), 155-158, doi:10.1029/GL017i002p00155.

class fatiando.seismic.wavefd.GaussSource(x, z, area, shape, amp, frequency, delay=None)[source]¶

Bases: fatiando.seismic.wavefd.MexHatSource

A wave source that vibrates as a Gaussian derivative wavelet.

\[\psi(t) = A 2 \sqrt{e}\ f\ t\ e^\left(-2t^2f^2\right)\]

Parameters:

x, z
: float
The x, z coordinates of the source
area
: [xmin, xmax, zmin, zmax]
The area bounding the finite difference simulation
shape
: (nz, nx)
The number of nodes in the finite difference grid
amp
: float
The amplitude of the source ($A$)
frequency
: float
The approximate frequency of the source
delay
: float
The delay before the source starts

Note

If you want the source to start with amplitude close to 0, use delay = 3.0/frequency.

coords()¶

Get the x, z coordinates of the source.

Returns:

(x, z)
: tuple
The x, z coordinates

indexes()¶

Get the i,j coordinates of the source in the finite difference grid.

Returns:

(i,j)
: tuple
The i,j coordinates

class fatiando.seismic.wavefd.MexHatSource(x, z, area, shape, amp, frequency, delay=0)[source]¶

Bases: object

A wave source that vibrates as a Mexican hat (Ricker) wavelet.

\[\psi(t) = A(1 - 2 \pi^2 f^2 t^2)exp(-\pi^2 f^2 t^2)\]

Parameters:

x, z
: float
The x, z coordinates of the source
area
: [xmin, xmax, zmin, zmax]
The area bounding the finite difference simulation
shape
: (nz, nx)
The number of nodes in the finite difference grid
amp
: float
The amplitude of the source ($A$)
frequency
: float
The peak frequency of the wavelet
delay
: float
The delay before the source starts

Note

If you want the source to start with amplitude close to 0, use delay = 3.5/frequency.

coords()[source]¶

Get the x, z coordinates of the source.

Returns:

(x, z)
: tuple
The x, z coordinates

indexes()[source]¶

Get the i,j coordinates of the source in the finite difference grid.

Returns:

(i,j)
: tuple
The i,j coordinates

class fatiando.seismic.wavefd.SinSqrSource(x, z, area, shape, amp, frequency, delay=0)[source]¶

Bases: fatiando.seismic.wavefd.MexHatSource

A wave source that vibrates as a sine squared function.

This source vibrates for a time equal to one period (T).

\[\psi(t) = A\sin\left(t\frac{2\pi}{T}\right)^2\]

Parameters:

x, z
: float
The x, z coordinates of the source
area
: [xmin, xmax, zmin, zmax]
The area bounding the finite difference simulation
shape
: (nz, nx)
The number of nodes in the finite difference grid
amp
: float
The amplitude of the source ($A$)
frequency
: float
The frequency of the source
delay
: float
The delay before the source starts

Note

If you want the source to start with amplitude close to 0, use delay = 3.5/frequency.

coords()¶

Get the x, z coordinates of the source.

Returns:

(x, z)
: tuple
The x, z coordinates

indexes()¶

Get the i,j coordinates of the source in the finite difference grid.

Returns:

(i,j)
: tuple
The i,j coordinates

fatiando.seismic.wavefd.blast_source(x, z, area, shape, amp, frequency, delay=0, sourcetype=<class 'fatiando.seismic.wavefd.MexHatSource'>)[source]¶

Uses several MexHatSources to create a blast source that pushes in all directions.

Parameters:

x, z
: float
The x, z coordinates of the source
area
: [xmin, xmax, zmin, zmax]
The area bounding the finite difference simulation
shape
: (nz, nx)
The number of nodes in the finite difference grid
amp
: float
The amplitude of the source
frequency
: float
The frequency of the source
delay
: float
The delay before the source starts
sourcetype
: source class
The type of source to use, like MexHatSource.

Returns:

[xsources, zsources]

Lists of sources for x- and z-displacements

fatiando.seismic.wavefd.elastic_psv(mu, lamb, density, area, dt, iterations, sources, stations=None, snapshot=None, padding=50, taper=0.002, xz2ps=False)[source]¶

Simulate P and SV waves using the Parsimonious Staggered Grid (PSG) finite differences scheme of Luo and Schuster (1990).

This is an iterator. It yields panels of $u_x$ and $u_z$ displacements and a list of arrays with recorded displacements in a time series. Parameter snapshot controls how often the iterator yields. The default is only at the end, so only the final panel and full time series are yielded.

Uses absorbing boundary conditions (Gaussian taper) in the lower, left and right boundaries. The top implements the free-surface boundary condition of Vidale and Clayton (1986).

Parameters:

mu
: 2D-array (shape = shape)
The $\mu$ Lame parameter at all the grid nodes
lamb
: 2D-array (shape = shape)
The $\lambda$ Lame parameter at all the grid nodes
density
: 2D-array (shape = shape)
The value of the density at all the grid nodes
area
: [xmin, xmax, zmin, zmax]
The x, z limits of the simulation area, e.g., the shallowest point is at zmin, the deepest at zmax.
dt
: float
The time interval between iterations
iterations
: int
Number of time steps to take
sources
: [xsources, zsources] : lists
A lists of the sources of waves for the particle movement in the x and z directions (see MexHatSource for an example source)
stations
: None or list
If not None, then a list of [x, z] pairs with the x and z coordinates of the recording stations. These are physical coordinates, not the indexes!
snapshot
: None or int
If not None, than yield a snapshot of the displacements at every snapshot iterations.
padding
: int
Number of grid nodes to use for the absorbing boundary region
taper
: float
The intensity of the Gaussian taper function used for the absorbing boundary conditions
xz2ps
: True or False
If True, will yield P and S wave panels instead of ux, uz. See xz2ps.

Yields:

[t, ux, uz, xseismograms, zseismograms]

The current iteration, the particle displacements in the x and z directions, lists of arrays containing the displacements recorded at each station until the current iteration.

References:

Vidale, J. E., and R. W. Clayton (1986), A stable free-surface boundary condition for two-dimensional elastic finite-difference wave simulation, Geophysics, 51(12), 2247-2249.

fatiando.seismic.wavefd.elastic_sh(mu, density, area, dt, iterations, sources, stations=None, snapshot=None, padding=50, taper=0.005)[source]¶

Simulate SH waves using the Equivalent Staggered Grid (ESG) finite differences scheme of Di Bartolo et al. (2012).

This is an iterator. It yields a panel of $u_y$ displacements and a list of arrays with recorded displacements in a time series. Parameter snapshot controls how often the iterator yields. The default is only at the end, so only the final panel and full time series are yielded.

Uses absorbing boundary conditions (Gaussian taper) in the lower, left and right boundaries. The top implements a free-surface boundary condition.

Parameters:

mu
: 2D-array (shape = shape)
The $\mu$ Lame parameter at all the grid nodes
density
: 2D-array (shape = shape)
The value of the density at all the grid nodes
area
: [xmin, xmax, zmin, zmax]
The x, z limits of the simulation area, e.g., the shallowest point is at zmin, the deepest at zmax.
dt
: float
The time interval between iterations
iterations
: int
Number of time steps to take
sources
: list
A list of the sources of waves (see MexHatSource for an example source)
stations
: None or list
If not None, then a list of [x, z] pairs with the x and z coordinates of the recording stations. These are physical coordinates, not the indexes
snapshot
: None or int
If not None, than yield a snapshot of the displacement at every snapshot iterations.
padding
: int
Number of grid nodes to use for the absorbing boundary region
taper
: float
The intensity of the Gaussian taper function used for the absorbing boundary conditions

Yields:

t, uy, seismograms
: int, 2D-array and list of 1D-arrays
The current iteration, the particle displacement in the y direction and a list of the displacements recorded at each station until the current iteration.

fatiando.seismic.wavefd.lame_lamb(pvel, svel, dens)[source]¶

Calculate the Lame parameter $\lambda$ P and S wave velocities ($\alpha$ and $\beta$) and the density ($\rho$).

\[\lambda = \alpha^2 \rho - 2\beta^2 \rho\]

Parameters:

pvel
: float or array
The P wave velocity
svel
: float or array
The S wave velocity
dens
: float or array
The density

Returns:

lambda
: float or array
The Lame parameter

Examples:

>>> print lame_lamb(2000, 1000, 2700)
5400000000
>>> import numpy as np
>>> pv = np.array([2000, 3000], dtype='float64')
>>> sv = np.array([1000, 1700], dtype='float64')
>>> dens = np.array([2700, 3100], dtype='float64')
>>> lamb = lame_lamb(pv, sv, dens)
>>> print "[ {:g}  {:g} ]".format(lamb[0], lamb[1])
[ 5.4e+09  9.982e+09 ]

fatiando.seismic.wavefd.lame_mu(svel, dens)[source]¶

Calculate the Lame parameter $\mu$ from S wave velocity ($\beta$) and the density ($\rho$).

\[\mu = \beta^2 \rho\]

Parameters:

svel
: float or array
The S wave velocity
dens
: float or array
The density

Returns:

mu
: float or array
The Lame parameter

Examples:

>>> print lame_mu(1000, 2700)
2700000000
>>> import numpy as np
>>> sv = np.array([1000, 1700], dtype='float64')
>>> dens = np.array([2700, 3100], dtype='float64')
>>> mu = lame_mu(sv, dens)
>>> print "[ {:g}  {:g} ]".format(mu[0], mu[1])
[ 2.7e+09  8.959e+09 ]

fatiando.seismic.wavefd.maxdt(area, shape, maxvel)[source]¶

Calculate the maximum time step that can be used in the simulation.

Uses the result of the Von Neumann type analysis of Di Bartolo et al. (2012).

Parameters:

area
: [xmin, xmax, zmin, zmax]
The x, z limits of the simulation area, e.g., the shallowest point is at zmin, the deepest at zmax.
shape
: (nz, nx)
The number of nodes in the finite difference grid
maxvel
: float
The maximum velocity in the medium

Returns:

maxdt
: float
The maximum time step

fatiando.seismic.wavefd.scalar(vel, area, dt, iterations, sources, stations=None, snapshot=None, padding=50, taper=0.005)[source]¶

Simulate scalar waves using an explicit finite differences scheme 4th order space. Space increment must be equal in x and z.

Parameters:

vel
: 2D-array (defines shape simulation)
The wave velocity at all the grid nodes
area
: [xmin, xmax, zmin, zmax]
The x, z limits of the simulation area, e.g., the shallowest point is at zmin, the deepest at zmax.
dt
: float
The time interval between iterations
iterations
: int
Number of time steps to take
sources
: list
A list of the sources of waves (see MexHatSource for an example source)
stations
: None or list
If not None, then a list of [x, z] pairs with the x and z coordinates of the recording stations. These are physical coordinates, not the indexes
snapshot
: None or int
If not None, than yield a snapshot of the displacement at every snapshot iterations.
padding
: int
Number of grid nodes to use for the absorbing boundary region
taper
: float
The intensity of the Gaussian taper function used for the absorbing boundary conditions

Yields:

i, u, seismograms
: int, 2D-array and list of 1D-arrays
The current iteration, the scalar quantity disturbed and a list of the scalar quantity disturbed recorded at each station until the current iteration.

fatiando.seismic.wavefd.scalar_maxdt(area, shape, maxvel)[source]¶

Calculate the maximum time step that can be used in the FD scalar simulation with 4th order space 1st time backward.

References

Alford R.M., Kelly K.R., Boore D.M. (1974) Accuracy of finite-difference modeling of the acoustic wave equation Geophysics, 39 (6), P. 834-842

Chen, Jing-Bo (2011) A stability formula for Lax-Wendroff methods with fourth-order in time and general-order in space for the scalar wave equation Geophysics, v. 76, p. T37-T42

Convergence

\[\Delta t \leq \frac{2 \Delta s}{ V \sqrt{\sum_{a=-N}^{N} (|w_a^1| + |w_a^2|)}} = \frac{ \Delta s \sqrt{3}}{ V_{max} \sqrt{8}}\]

Where w_a are the centered differences weights

Parameters:

area
: [xmin, xmax, zmin, zmax]
The x, z limits of the simulation area, e.g., the shallowest point is at zmin, the deepest at zmax.
shape
: (nz, nx)
The number of nodes in the finite difference grid
maxvel
: float
The maximum velocity in the medium

Returns:

maxdt
: float
The maximum time step

fatiando.seismic.wavefd.xz2ps(ux, uz, area)[source]¶

Convert the x and z displacements into representations of P and S waves using the divergence and curl, respectively, after Kennett (2002, pp 57)

\[P: \frac{\partial u_x}{\partial x} + \frac{\partial u_z}{\partial z}\]

\[S: \frac{\partial u_x}{\partial z} - \frac{\partial u_z}{\partial x}\]

Parameters:

ux, uz
: 2D-arrays
The x and z displacement panels
area
: [xmin, xmax, zmin, zmax]
The x, z limits of the simulation area, e.g., the shallowest point is at zmin, the deepest at zmax.

Returns:

p, s
: 2D-arrays
Panels corresponding to P and S wave components

References:

Kennett, B. L. N. (2002), The Seismic Wavefield: Volume 2, Interpretation of Seismograms on Regional and Global Scales, Cambridge University Press.

Finite difference solution of the 2D wave equation (fatiando.seismic.wavefd)¶

Finite difference solution of the 2D wave equation (`fatiando.seismic.wavefd`)¶