# 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

Auxiliary functions

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]

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]

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.

: 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.

: 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.

: 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.