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 schemeSources
MexHatSource
: Mexican hat wavelet sourceSinSqrSource
: Sine squared sourceGaussSource
: Gauss derivative sourceblast_source
: A source blasting in all
directionsAuxiliary functions
lame_lamb
: Calculate the lambda Lame
parameterlame_mu
: Calculate the mu Lame parameterxz2ps
: Convert x and z displacements to
representations of P and S wavesmaxdt
: Calculate the maximum time step for
elastic wave simulationsscalar_maxdt
: Calculate the maximum time
step for a scalar wave simulationTheory
A good place to start is the equation of motion for elastic isotropic materials
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
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
Substituting the stress formula in the above equation yields
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
And the corresponding stress components are
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.
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.
Parameters:
The x, z coordinates of the source
The area bounding the finite difference simulation
The number of nodes in the finite difference grid
The amplitude of the source (\(A\))
The approximate frequency of the source
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:
The x, z coordinates
indexes
()¶Get the i,j coordinates of the source in the finite difference grid.
Returns:
The i,j coordinates
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.
Parameters:
The x, z coordinates of the source
The area bounding the finite difference simulation
The number of nodes in the finite difference grid
The amplitude of the source (\(A\))
The peak frequency of the wavelet
The delay before the source starts
Note
If you want the source to start with amplitude close to 0,
use delay = 3.5/frequency
.
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).
Parameters:
The x, z coordinates of the source
The area bounding the finite difference simulation
The number of nodes in the finite difference grid
The amplitude of the source (\(A\))
The frequency of the source
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:
The x, z coordinates
indexes
()¶Get the i,j coordinates of the source in the finite difference grid.
Returns:
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:
The x, z coordinates of the source
The area bounding the finite difference simulation
The number of nodes in the finite difference grid
The amplitude of the source
The frequency of the source
The delay before the source starts
The type of source to use, like
MexHatSource
.
Returns:
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:
The \(\mu\) Lame parameter at all the grid nodes
The \(\lambda\) Lame parameter at all the grid nodes
The value of the density at all the grid nodes
The x, z limits of the simulation area, e.g., the shallowest point is at zmin, the deepest at zmax.
The time interval between iterations
Number of time steps to take
A lists of the sources of waves for the particle movement in the x and
z directions
(see MexHatSource
for an example
source)
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!
If not None, than yield a snapshot of the displacements at every snapshot iterations.
Number of grid nodes to use for the absorbing boundary region
The intensity of the Gaussian taper function used for the absorbing boundary conditions
If True, will yield P and S wave panels instead of ux, uz. See
xz2ps
.
Yields:
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:
The \(\mu\) Lame parameter at all the grid nodes
The value of the density at all the grid nodes
The x, z limits of the simulation area, e.g., the shallowest point is at zmin, the deepest at zmax.
The time interval between iterations
Number of time steps to take
A list of the sources of waves
(see MexHatSource
for an example
source)
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
If not None, than yield a snapshot of the displacement at every snapshot iterations.
Number of grid nodes to use for the absorbing boundary region
The intensity of the Gaussian taper function used for the absorbing boundary conditions
Yields:
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\)).
Parameters:
The P wave velocity
The S wave velocity
The density
Returns:
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\)).
Parameters:
The S wave velocity
The density
Returns:
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:
The x, z limits of the simulation area, e.g., the shallowest point is at zmin, the deepest at zmax.
The number of nodes in the finite difference grid
The maximum velocity in the medium
Returns:
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:
The wave velocity at all the grid nodes
The x, z limits of the simulation area, e.g., the shallowest point is at zmin, the deepest at zmax.
The time interval between iterations
Number of time steps to take
A list of the sources of waves
(see MexHatSource
for an example
source)
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
If not None, than yield a snapshot of the displacement at every snapshot iterations.
Number of grid nodes to use for the absorbing boundary region
The intensity of the Gaussian taper function used for the absorbing boundary conditions
Yields:
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
Where w_a are the centered differences weights
Parameters:
The x, z limits of the simulation area, e.g., the shallowest point is at zmin, the deepest at zmax.
The number of nodes in the finite difference grid
The maximum velocity in the medium
Returns:
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)
Parameters:
The x and z displacement panels
The x, z limits of the simulation area, e.g., the shallowest point is at zmin, the deepest at zmax.
Returns:
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.