Pyrocko-GF - Geophysical forward modelling with pre-calculated Green’s functions¶

Introduction ¶

Many seismological methods make use of numerically calculated Green’s functions (GFs). Often these are required for a vast number of combinations of source and receiver coordinates, e.g. when computing synthetic seismograms in seismic source inversions. Calculation of the GFs is a computationally expensive operation and it can be of advantage to calculate them in advance. The same GF traces can then be reused many times as required in a typical application.

To efficiently reuse GFs across different applications, they must be stored in a common format. In such a form, they can also be passed to fellow researchers.

Furthermore, it is useful to store associated meta information, like e.g. travel time tables for seismic phases and the earth model used, together with the GF in order to have a complete and consistent framework to play with.

Overview of the components and interaction for Pyrocko-GF

Pyrocko contains a flexible framework to store and work with pre-calculated GFs. It is implemented in the pyrocko.gf sub-package. Also included, is a powerful front end tool to create, inspect, and manipulate GF stores: the fomosto tool (“forward model storage tool”).

Media models ¶

Where Pyrocko uses layered 1D earth models in the calculation of GFs it uses the TauP .nd (named discontinuity) format. These files are ascii tables with the following columns:

Structure of a named discontinuities file (.nd).¶

depth [km]    Vp[km/s]    Vs[km/s]    density[g/cm^3]    qp    qs

The fomosto application description holds an exemplary earth model configuration. Users can define own input models before calculation. Also a number of predefined common variations of AK135 and PREM models are available. They can be listed and inspected using the cake command line tool.

Green’s function stores ¶

GF pre-calculation¶

Calculating and managing Pyrocko-GF stores is accomplished by Pyrocko’s fomosto application. More details how GF stores are set-up and calculated can be found in the fomosto tutorial.

Downloading GF stores¶

Calculating and quality checking GF stores is a time consuming task. Many pre-calculated stores can be downloaded from our online repository, the Green’s mill (greens-mill.pyrocko.org).

The available stores include dynamic stores for simulating waveforms at global and regional extents, as well as static stores for the modelling of step-like surface displacements.

Downloading a store with fomosto¶

fomosto download kinherd global_2s_25km

Using GF Stores¶

GF stores are accessed for forward modelling by the Pyrocko-GF Engine. Here is how we can start-up the engine for modelling:

Import and initialise the forward modelling engine.¶

from pyrocko.gf import LocalEngine

engine = LocalEngine(store_dirs=['gf_stores/global_2s/'])

A complete list of arguments can be found in the library reference, LocalEngine.

Source models ¶

Pyrocko-GF supports the simulation of various dislocation sources, focused on earthquake and volcano studies.

Note

Multiple sources can be combined through the CombiSource object.

Point sources¶

For convenience, different parameterizations of seismological moment tensor point sources are available.

Source	Short description
`ExplosionSource`	An isotrope moment tensor for explosions or volume changes.
`DCSource`	Double force couple, for pure-shear earthquake ruptures.
`MTSource`	Full moment tensor representation of force excitation.
`CLVDSource`	A pure compensated linear vector dipole source.
`VLVDSource`	Volumetric linear vector dipole, a rotational symmetric volume source.
`SFSource`	A 3-component single force point source.
`PorePressurePointSource`	Excess pore pressure point source.

Finite sources¶

Source	Short description
`RectangularSource`	Rectangular fault plane.
`RingfaultSource`	Ring fault for volcanic processes, e.g. caldera collapses.
`DoubleDCSource`	Relative parameterization of a twin double couple source.
`PorePressureLineSource`	Excess pore pressure line source

First import the Pyrocko-GF framework with

Import all object from pyrocko.gf.¶

from pyrocko import gf

Explosion source¶

An isotropic explosion point source, which can also be used for dislocations due to volume changes.

Initialise a simple explosion source with a volume¶

explosion = gf.ExplosionSource(lat=42., lon=22., depth=8e3, volume_change=5e8)

Double couple¶

A double-couple point source, describing shear ruptures.

Initialise a double-couple source.¶

dc_source = gf.DCSource(lat=54., lon=7., depth=5e3, strike=33., dip=20., rake=80.)

Moment tensor¶

A moment tensor point source. This is the most complete form of describing an ensemble of buried forces to first order.

Initialise a full moment tensor.¶

mt_source = gf.MTSource(
   lat=20., lon=58., depth=8.3e3,
   mnn=.5, mee=.1, mdd=.7,
   mne=.6, mnd=.2, med=.1,)

# Or use an event
mt_source = MTSource.from_pyrocko_event(event)

CLVD source¶

A compensated linear vector dipole (CLVD) point source.

Initialise a CLVD source.¶

clvd_source = gf.CLVDSource(
    lat=48., lon=17., depth=5e3, dip=31., depth=5e3, azimuth=83.)

VLVD source¶

A volumetric linear vector dipole, a uniaxial rotational symmetric moment tensor source. This source can be used to constrain sill or dyke like volume dislocation.

Initialise a VLVD source.¶

vlvd_source = gf.VLVDSource(
   lat=-30., lon=184., depth=5e3,
   volume_change=1e9, clvd_moment=20e9, dip=10., azimuth=110.)

Rectangular fault¶

Classical Haskell finite source model, modified for bilateral rupture.

Initialise a rectangular fault with a width of 3 km, a length of 8 km and slip of 2.3 m.¶

km = 1e3

rectangular_source = gf.RectangularSource(
    lat=20., lon=44., depth=5*km,
    dip=30., strike=120., rake=50.,
    width=3*km, length=8*km, slip=2.3)

Ring fault¶

A ring fault with vertical double couples. Ring faults can describe volcanic processes, e.g. caldera collapses.

Initialise a dipping ring fault.¶

ring_fault = gf.RingFault(
    lat=31., lon=12., depth=2e3,
    diameter=5e3, sign=1.,
    dip=10., strike=30.,
    npointsources=50)

Source Time Functions ¶

Source time functions describe the normalized moment rate of a source point as a function of time. A number of source time functions (STF) are available and can be applied in pre- or post-processing. If no specific STF is defined a unit pulse response is assumed.

STF	Short description
`BoxcarSTF`	Boxcar shape source time function.
`TriangularSTF`	Triangular shape source time function.
`HalfSinusoidSTF`	Half sinusoid type source time function.
`ResonatorSTF`	A simple resonator like source time function.

Boxcar STF¶

A boxcar source time function. In the plot, each point is representative of the STF’s integral in the time interval $[-\Delta t/2, +\Delta t/2]$ surrounding it ( $\Delta t$ is the sampling interval).

Initialise an boxcar STF with duration of 5 s and centred at the centroid time.¶

stf = gf.BoxcarSTF(5., center=0.)

Triangular STF¶

A triangular shaped source time function. It can be made asymmetric.

Initialise a symmetric triangular STF with duration 5 s, which reaches its maximum amplitude after half the duration and centred at the centroid time.¶

stf = gf.TriangularSTF(5., peak_ratio=0.5, center=0.)

Half sinusoid STF¶

A half-sinusoid source time function.

Initialise a half sinusoid type STF with a duration of 5 s and centred around the centroid time.¶

stf = gf.HalfSinusoidSTF(5., center=0.)

Resonator STF¶

Initialise a resonator STF with duration of 5 s and a resonance frequency of 1 Hz.¶

stf = gf.ResonatorSTF(5., frequency=1.0)

Modelling targets ¶

Pyrocko-GF Targets are data structures holding observer properties to tell the framework what we want to model, e.g. whether we want to model a waveform or spectrum at a specific receiver site or displacement values at a set of locations. Each target has properties (location, depth, physical quantity) and essentially is associated to a GF store, used for modelling. The target also defines the method used to interpolate the discrete, gridded GF components. Please also see the Pyrocko GF modelling example.

Note

In Pyrocko locations are given with five coordinates: lat, lon, east_shift, north_shift and depth.

Latitude and longitude are the origin of an optional local Cartesian coordinate system for which an east_shift and a north_shift [m] can be defined. A target has a depth below the surface. However, the surface can have topography and the target can also have an elevation.

Waveforms¶

Objects of the class Target are used to calculate seismic waveforms. They define the geographical location (e.g. the station), component orientation (e.g. vertical or radial), physical quantity, and optionally a time interval

# Define a list of pyrocko.gf.Target objects, representing the recording
# devices. In this case one three-component seismometer is represented with
# three distinct target objects. The channel orientations are guessed from
# the channel codes here.
waveform_targets = [
    gf.Target(
       quantity='displacement',
       lat=10., lon=10.,
       store_id='global_2s_25km',
       codes=('NET', 'STA', 'LOC', channel_code))
    for channel_code in ['E', 'N', 'Z']

See the forward modelling example for a complete Python script and further explanation.

Static surface displacements¶

Modelling of step-like surface displacements is configured with StaticTarget objects. The resulting displacements have no time dependence, but can hold many locations. Special forms derive from the StaticTarget class:

the SatelliteTarget, for the forward modelling of InSAR data, and
the GNSSCampaignTarget for e.g. step-like GPS displacements.

Initialising a StaticTarget.¶

# east and north are numpy.ndarrays in meters
import numpy as num

km = 1.0e3
norths = num.linspace(-20*km, 20*km, 100)
easts = num.linspace(-20*km, 20*km, 100)
north_shifts, east_shifts = num.meshgrid(norths, easts)

static_target = gf.StaticTarget(
    lats=43., lons=20.,
    north_shifts=north_shifts,
    east_shifts=east_shifts,
    interpolation='nearest_neighbor',
    store_id='ak135_static')

The SatelliteTarget defines the locations of displacement measurements and the direction of the measurement, which is the so-called line-of-sight of the radar. See the forward modelling examples for detailed instructions of usage.

Initialising a SatelliteTarget.¶

# east/north shifts as numpy.ndarrays in [m]
# line-of-sight angles are NumPy arrays,
# - phi is _towards_ the satellite clockwise from east in [rad]
# - theta is the elevation angle from the horizon

satellite_target = gf.SatelliteTarget(
    lats=43., lons=20.,
    north_shifts=north_shifts,
    east_shifts=east_shifts,
    interpolation='nearest_neighbor',
    phi=phi,
    theta=theta,
    store_id='ak135_static')

The GNSSCampaignTarget defines station locations and the three components: east, north and up.

Forward modelling with Pyrocko-GF ¶

Forward modelling, given a source and target description, is handled in the so-called Engine using the process() method.

Initialisation of the engine requires setting the folder, where it should look for GF stores. This can be configured globally by setting the store_superdirs entry in file ~/.pyrocko/config.pf or locally using the initialization arguments of the LocalEngine.

Note, that modelling of dynamic targets (displacement waveforms) requires GFs that have many samples in time and modelling of static targets (for step-like displacements) usually only one. It is therefore meaningful to use dynamic GF stores for dynamic targets and static stores for static targets.

Forward modelling dynamic waveforms¶

For waveform targets, Pyrocko Trace objects representing the resulting waveforms can be obtained from the engine’s response.

forward model wave forms of a DoubleCouple point.¶

# Setup the LocalEngine and point it to the GF store you want to use.
# `store_superdirs` is a list of directories where to look for GF Stores.
engine = gf.LocalEngine(store_superdirs=['/data/gf_stores'])

# The computation is performed by calling process on the engine
response = engine.process(dc_source, waveform_targets)

# convert results in response to Pyrocko traces
synthetic_traces = response.pyrocko_traces()

# visualise the response with the snuffler
synthetic_traces.snuffle()

Forward modelling static surface displacements¶

For static targets, the results are retrieved in the following way:

forward model static surface displacements of a rectangular fault¶

# Get a default engine (will look into directories configured in
# ~/.pyrocko/config.pf to find GF stores)
engine = gf.get_engine()

response = engine.process(rectangular_source, satellite_target)

# Retrieve a list of static results:
synth_disp = response.static_results()

For regularly gridded satellite targets, the engine’s response can be converted to a synthetic Kite scene:

forward modelling from an existing kite scene.¶

from pyrocko import gf
from kite import Scene

km = 1e3
engine = gf.LocalEngine(use_config=True)

scene = Scene.load('sentinel_scene.npz')

src_lat = 37.08194 + .045
src_lon = 28.45194 + .2

source = gf.RectangularSource(
    lat=src_lat,
    lon=src_lon,
    depth=2*km,
    length=4*km, width=2*km,
    strike=45., dip=60.,
    slip=.5, rake=0.,
    anchor='top')

target = gf.KiteSceneTarget(scene, store_id='ak135_static')

result = engine.process(source, target, nthreads=0)

mod_scene = result.kite_scenes()[0]
mod_scene.spool()

Navigation

GF pre-calculation¶

Downloading GF stores¶

Using GF Stores¶

Point sources¶

Finite sources¶

Explosion source¶

Double couple¶

Moment tensor¶

CLVD source¶

VLVD source¶

Rectangular fault¶

Ring fault¶

Boxcar STF¶

Triangular STF¶

Half sinusoid STF¶

Resonator STF¶

Waveforms¶

Static surface displacements¶

Forward modelling dynamic waveforms¶

Forward modelling static surface displacements¶

Navigation