Calculating seismic observables

pyprop8 can calculate three kinds of seismic observables: seismic spectra, displacement seismograms, and static offsets (e.g. for simulating GPS or InSAR data). This document describes how to set up and perform such simulations. For additional details, see the API reference.

Problem setup

In order to calculate seismograms or seismic spectra, we need to specify three things:

  1. The properties of the earth model;

  2. The seismic source; and

  3. The receivers at which the wavefield is to be recorded.

Earth model

In pyprop8, the earth model is assumed to be an infinite, layered halfspace. Each layer is characterised by four parameters: a layer thickness (‘thickness’, km), a P-wave velocity (‘vp’, km/s), an S-wave velocity (‘vs’, km/s), and a density (‘rho’, g/cm^3).

To represent a model, we need to create an instance of LayeredStructureModel. The easiest way to achieve this is to initialise the class, passing a list (or list-like object) containing the properties of each layer. Two formats are supported.

  1. Each entry in the list is a tuple (thickness, vp, vs, rho). The first entry in the list is assumed to correspond to the uppermost (surface) layer, with subsequent entries representing successively deeper layers. The final entry in the list specifies the properties of the underlying half-space, and its thickness must be specified as np.inf. This is the default format for LayeredStructureModel.

  2. Alternatively, each entry in the list may be a tuple (ztop, vp, vs, rho), where ztop represents the depth (in km) to the interface defining the top of the layer. In this case, list entries need not be ordered. One entry must have ztop=0, defining the properties of the surface layer; the ordering and thickness of all other layers will be automatically inferred. If there is duplication of interface depths within the list, only the first will be retained. The properties associated with the lowermost interface specified are deemed to represent the underlying halfspace. To use this format, pass the optional argument interface_depth_form=True when initialising LayeredStructureModel.

The choice of format is simply a matter of user preference and convenience; it has no impact on the manner in which the simulations proceed. Thus:

layers = [(3.0,    1.8,   0, 1.02),
          (5.0,    4.5, 2.4, 2.5),
          (np.inf, 8.0, 4.5, 3.0)]
model = pp.LayeredStructureModel(layers)

will behave identically to:

layers = [(0.0, 1.8,   0, 1.02),
          (3.0, 4.5, 2.4, 2.5),
          (8.0,   8, 4.5, 3.0)]
model = pp.LayeredStructureModel(layers, interface_depth_form=True)

It will be noted that the surface layer in this case has a shear-wave velocity of 0 km/s. This signifies that it is a fluid – this particular model might correspond to a location within an ocean. Only the uppermost layer is allowed to be a fluid.

Once you have created your model, you can call print(model) to obtain a graphical representation:

>>> print(model)
------------------------------------------------------- z = 0.00 km
  vp = 1.80 km/s       FLUID        rho = 1.02 g/cm^3
------------------------------------------------------- z = 3.00 km
  vp = 4.50 km/s   vs = 2.40 km/s   rho = 2.50 g/cm^3
------------------------------------------------------- z = 8.00 km
  vp = 8.00 km/s   vs = 4.50 km/s   rho = 3.00 g/cm^3

Receivers

Next we need to specify receiver locations. pyprop8 supports two formats, which result in different computation paths internally and hence different costs.

The first is designed to be used when the wavefield is sought at isolated, pre-determined locations - perhaps locations chosen to correspond to a deployment of instruments. In this case, we need to create an instance of ListOfReceivers. This requires the x- and y-coordinates of each station, and the depth at which instruments are buried (which should be given as ‘0’ if receivers lie on the surface):

stations = np.array([[ 5.3,  6.2],
                     [-2.1,  0.3],
                     [ 1.5, -3.1]]) # 3 pairs of (x,y) coordinates
depth = 4 # Seafloor instruments
receivers = pp.ListOfReceivers(stations[:,0],stations[:,1],depth)

The computation algorithm currently requires all receivers to lie at a common depth. If multiple depths are required, it will be necessary to divide receivers between multiple ListOfReceivers objects and perform simulations for each.

ListOfReceivers also accepts an optional argument, geometry. By default, this has the value geometry='cartesian', implying that x- and y-coordinates are expressed relative to a Cartesian kilometre grid (i.e. a location (1,2) lies 1km east and 2km north of some arbitrarily-chosen origin point). The source location (see below) is assumed to be specified within the same coordinate system, and this is used to determine the source-receiver distances and azimuths used internally within the code. Alternatively, geometry='spherical' may be specified. In this case, x- and y-coordinates are treated as degrees longitude and degrees latitude. This is convenient for approximating real-world scenarios, and source-receiver distances and azimuths are calculated assuming great-circle propagation on the surface of a sphere of radius 6371km. However, pyprop8 remains based on a flat-Earth approximation, and its validity degrades as the source-receiver distance increases.

Alternatively, we may wish to achieve a general characterisation of the wavefield throughout a region. For this case, pyprop8 provides RegularlyDistributedReceivers, which assumes that receivers lie on a regular polar grid, centred upon the event location. This creates geometric simplifications that can be exploited by the computation algorithm, substantially accelerating calculations. For this case, we must specify ranges for radius and azimuth, and the number of grid points for each:

# 5 equally spaced radii between 10km and 50km (i.e. [10 20 30 40 50])
rmin, rmax = (10, 50)
nr = 5
# 8 equally spaced azimuths between 0 and 315 degrees (i.e. every 45 degrees)
phimin, phimax = (0, 360)
nphi = 8
# Receivers may still be buried
depth = 0
receivers = pp.RegularlyDistributedReceivers(rmin,rmax,nr,phimin,phimax,nphi,depth)

This will result in a regular polar grid of 40 stations. By default, it is assumed that the minimum and maximum azimuths are measured in degrees (counter-clockwise from the x/East axis, when viewed from above); pass an optional argument degrees=False to use radians.

As already mentioned, the use of RegularlyDistributedReceivers requires the event to lie at the centre of the polar grid. The location of this central point can be specified with respect to a global Cartesian coordinate system by passing optional arguments x0= and y0=. If these are omitted, the grid is assumed to be centred on the origin of the Cartesian system, and the source location should also be specified (see below) as (x,y) = (0,0).

Seismic source

pyprop8 assumes that the seismic source acts at a single point in space (and time; but see the discussion of source time-functions, below). To represent a seismic source, we need to create an instance of PointSource. This requires us to specify the source location (an x-coordinate, a y-coordinate, and a depth); the source mechanism (expressed as a co-located moment tensor and force vector); and an event time:

event_x, event_y, event_dep = ( 5, -3, 10) # Spatial location
M = np.array([[ 1, 0, 0],
              [ 0,-1, 0],
              [ 0, 0, 0]]) # Moment tensor, expressed in Cartesian system
F = np.array([[ 0],
              [ 0],
              [ 0]]) # Force vector, expressed in Cartesian system
event_time = datetime.datetime.fromisoformat("2021-08-17T03:45:37") # Date/time
source = pp.PointSource(event_x, event_y, event_dep, M, F, event_time)

The lateral (x/y) source coordinates should be specified in the same coordinate system as is used for the receivers (see above). This may be either a Cartesian kilometre grid, with arbitrary origin; or longitude and latitude coordinates specified in degrees.

Moment tensors and force vectors are expressed in a Cartesian, z-up system. The function utils.rtf2xyz() is available to convert moment tensors from the spherical-polar definition that is ubiquitous in global seismology (e.g. if catalogue source mechanisms are to be used). Additionally, the function utils.make_moment_tensor() is available to construct a moment tensor (in a spherical-polar system) given strike, dip and rake angles and magnitude information. Note that both moment tensor and force vector must be supplied, although all entries may be zero.

The moment tensor is represented by a 3x3 array, and the force vector by a 1x3 array. It is also possible to specify multiple moment tensor/force vector combinations within a single instance of PointSource:

M = np.array([[[ 1, 0, 0],
               [ 0,-1, 0],
               [ 0, 0, 0]],
              [[ 0, 0, 0],
               [ 0, 0, 0],
               [ 0, 0, 0]]]) # Two 3x3 moment tensors
F = np.array([[[ 0],
               [ 0],
               [ 0]],
              [[ 1],
               [ 0],
               [ 0]]]) # Two 3x1 force vectors
source = pp.PointSource(event_x, event_y, event_dep, M, F, event_time)

In this case, the shape of the moment tensor array becomes Nx3x3, and that of the force vector Nx1x3. As described below, pyprop8 will perform separate simulations for each of the N moment tensor/force vector pairs, returning N sets of simulated data. This capability is exploited internally to assist in calculation of derivatives, but may have limited value for general use: any computational efficiencies are modest, and potentially offset by the growth in memory demands.

Obtaining seismic spectra

Once the model, receivers and source have been created, obtaining seismic spectra (i.e., the spectrum of the wavefield at each receiver) can be as simple as:

om_min, om_max = (0,5*2*np.pi)
nfreq = 1000
omegas = np.linspace(om_min,om_max,nfreq) # Frequencies at which spectrum is sought
spectra = compute_spectra(model, source, receivers, omegas, squeeze_outputs=False) # Do the calculation

This creates the (complex) array spectra, containing a spectrum of the wavefield at each receiver. The shape of the array depends on the choice made when specifying the receivers.

  1. If receivers is an instance of ListOfReceivers, spectra will have shape (source.nsources, receivers.nstations, 3, nfreq), where source.nsources is the number of moment tensor/source vector pairs specified within the source object, receivers.nstations is the total number of receivers, and nfreq is the number of frequency points at which evaluation was requested. The third dimension indexes the three components of motion: radial, transverse, and vertical.

  2. If receivers is an instance of RegularlyDistributedReceivers, spectra will have shape (source.nsources, receivers.nr, receivers.nphi, 3, nfreq), where receivers.nr and receivers.nphi are the number of grid points in the radial and azimuthal directions, respectively.

However, if squeeze_outputs=True (which is the default if this argument is omitted), numpy.squeeze() will be applied to the output of compute_spectra. This discards any dimensions that have size 1.

Derivatives

To obtain derivatives of spectra with respect to source parameters, we need to first create an instance of DerivativeSwitches. This is used to specify the derivatives that are sought, for example:

derivs = DerivativeSwitches(moment_tensor=True, force=False, x=True, y=True, z=True)

We can then call compute_spectra, passing this object as an optional derivatives argument:

spectra, derivatives = compute_spectra(model, source, receivers, omegas, derivatives=derivs, squeeze_outputs=False)

Notice that when derivatives is set (or more precisely, when it receives any value other than None), compute_spectra now returns two arrays. The spectra array is organised precisely as described above; the derivatives array has an additional dimension. Again, there are two possibilities:

  1. If receivers is an instance of ListOfReceivers, derivatives will have shape (source.nsources, receivers.nstations, derivs.nderivs, 3, nfreq), where source.nderivs is the total number of derivatives requested.

  2. If receivers is an instance of RegularlyDistributedReceivers, derivatives will have shape (source.nsources, receivers.nr, receivers.nphi, derivs.nderivs, 3, nfreq).

If squeeze_outputs=True, numpy.squeeze() will also be applied to this array, discarding any dimensions of size ‘1’.

The index of specific derivatives (e.g. the ‘depth’ derivative) within the array will vary depending on the components that are requested. Therefore, DerivativeSwitches provides a mechanism for automatically determining the appropriate index. For example:

# Derivatives wrt the six independent moment tensor components
dmt = derivatives[:,:,derivs.i_mt:derivs.i_mt+6,:,:]
# Derivatives wrt the event depth
ddepth = derivatives[:,:,derivs.i_dep,:,:]

Obtaining seismograms

pyprop8 provides a separate routine to compute seismograms (i.e., time series). Fundamentally, this is simply a matter of obtaining spectra and then taking the Fourier transform; however, compute_seismograms() handles various additional book-keeping and processing tasks. At the simplest:

nt = 200 # Number of time-series points
dt = 0.5   # Sampling interval (s)
tt, seismograms = compute_seismograms(model, source, receivers, nt, dt, squeeze_outputs=False):

This computes an nt-point displacement time series for each receiver, sampled every dt seconds. Two arrays are returned. The first (tt) has dimension (nt, ) and contains the time-points at which the seismogram is obtained (i.e. the sequence [0, dt, 2*dt, ..., (nt-1)*dt]). The second is again dependent on the manner in which receivers are specified:

  1. If receivers is an instance of ListOfReceivers, seismograms will have shape (source.nsources, receivers.nstations, 3, nt), where source.nsources is the number of moment tensor/source vector pairs specified within the source object, receivers.nstations is the total number of receivers, and nfreq is the number of frequency points at which evaluation was requested. The third dimension indexes the three components of motion: x, y, and z.

  2. If receivers is an instance of RegularlyDistributedReceivers, seismograms will have shape (source.nsources, receivers.nr, receivers.nphi, 3, nt), where receivers.nr and receivers.nphi are the number of grid points in the radial and azimuthal directions, respectively.

Again, using the default option of squeeze_outputs=True will apply numpy.squeeze() to the the output array.

By default, seismograms represent ground displacment relative to a Cartesian coordinate system. By passing the optional argument xyz=False it is possible to instead obtain seismograms relative to a polar basis, i.e. radial/transverse/vertical motion.

Source time-functions

By default, pyprop8 is based on the assumption that the source behaves like a step-function in time: all energy is released instantaneously. This is common in seismological simulations, but is not a realistic representation of most seismic sources. It is therefore important to convolve the seismograms with a source time-function that provides a more reasonable representation of energy release. It is efficient to implement this convolution as a multiplication in the frequency domain, and compute_seismograms() has the facility to apply such a transformation to the seismic spectra prior to taking the Fourier transform.

To use this, it is necessary to determine the frequency spectrum of the desired source time-function. This should be implemented as a Python function (or other callable) that takes an angular frequency as its sole argument, and returns the (complex) value of the source spectrum at that point, e.g.:

def stf(om):
    if om==0:
        f = 1
    else:
        f = np.sin(om)/om
   return f

This callable should then be passed via the source_time_function keyword argument, e.g. source_time_function=stf. For convenience and illustration, some standard functions are provided in the utils module.

Derivatives

Computing derivatives of seismograms with respect to moment tensor components follows the pattern already described in the context of spectra. As before, it is necessary to create an instance of DerivativeSwitches, and pass it to compute_seismograms() via the keyword argument derivatives:

derivs = DerivativeSwitches(moment_tensor=True, force=False, x=True, y=True, depth=True)
tt, seismograms, derivatives = compute_seismograms(model, source, receivers, nt, dt, derivatives=derivs, squeeze_outputs=False)

Again, this will cause compute_seismograms() to return a third array, derivatives.

  1. If receivers is an instance of ListOfReceivers, derivatives will have shape (source.nsources, receivers.nstations, derivs.nderivs, 3, nt), where source.nderivs is the total number of derivatives requested.

  2. If receivers is an instance of RegularlyDistributedReceivers, derivatives will have shape (source.nsources, receivers.nr, receivers.nphi, derivs.nderivs, 3, nt).

If squeeze_outputs=True, numpy.squeeze() will be applied to this array. Again, it is recommended to index the array using the functionality provided in DerivativeSwitches.

Obtaining static offset measurements

Finally, it is also possible to simulate only the static offset (i.e. the part of seismic displacement that remains once transient motion has ceased). This follows the same pattern as the functions already described; usage can be as simple as:

static = compute_static(model, source, receivers)

This will return an an array of shape (source.nsources, receivers.nstations, 3) or (source.nsources, receivers.nr, receivers.nphi, 3), depending on the manner in which receivers is specified. This contains displacements expressed relative to a Cartesian x/y/z basis, indexed by the final dimension. Alternatively, it is possible to specify one or more ‘line of sight’ vectors, which determine the direction(s) in which displacement is to be measured. Thus, for example:

los = np.array([[1,0],
                [1,0],
                [0,1]])
static = compute_static(model, source, receivers,los_vector=los)

will return an array of shape (source.nsources, receivers.nstations, 2) or or (source.nsources, receivers.nr, receivers.nphi, 2), containing displacements relative to the two specified directions. Again, passing squeeze_outputs=True will eliminate unnecessary dimensions from the array.

As before, derivatives may be obtained by passing an additional derivatives=derivs argument, where derivs is an instance of DerivativeSwitches. Thus:

static, derivatives = compute_static(model, source, receivers, los_vector=los, derivatives=derivs)

will result in derivatives having shape (source.nsources, receivers.nstations, derivs.nderivs, nlos) or (source.nsources, receivers.nr, receivers.nphi, derivs.nderivs, nlos), where nlos is the number of line-of-sight vectors provided (or ‘3’ in the default case).