PyrockoGF  Geophysical forward modelling with precalculated 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.
Pyrocko contains a flexible framework to store and work with precalculated
GFs. It is implemented in the pyrocko.gf
subpackage. 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:
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 precalculation¶
Calculating and managing PyrockoGF stores is accomplished by Pyrocko’s fomosto application. More details how GF stores are setup and calculated can be found in the fomosto tutorial.
Downloading GF stores¶
Calculating and quality checking GF stores is a time consuming task. Many precalculated stores can be downloaded from our online repository, the Green’s mill (greensmill.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 steplike surface displacements.
fomosto download kinherd global_2s_25km
Using GF Stores¶
GF stores are accessed for forward modelling by the PyrockoGF Engine
. Here is how we can startup the engine for modelling:
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¶
PyrockoGF 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 

An isotrope moment tensor for explosions or volume changes. 

Double force couple, for pureshear earthquake ruptures. 

Full moment tensor representation of force excitation. 

A pure compensated linear vector dipole source. 

Volumetric linear vector dipole, a rotational symmetric volume source. 

A 3component single force point source. 

Excess pore pressure point source. 
Finite sources¶
Source 
Short description 

Rectangular fault plane. 

Ring fault for volcanic processes, e.g. caldera collapses. 

Relative parameterization of a twin double couple source. 

Excess pore pressure line source 
First import the PyrockoGF framework with
from pyrocko import gf
Explosion source¶
An isotropic explosion point source, which can also be used for dislocations due to volume changes.
explosion = gf.ExplosionSource(lat=42., lon=22., depth=8e3, volume_change=5e8)
Double couple¶
A doublecouple point source, describing shear ruptures.
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.
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.
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.
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.
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.
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 postprocessing. If no specific STF is defined a unit pulse response is assumed.
STF 
Short description 

Boxcar shape source time function. 

Triangular shape source time function. 

Half sinusoid type source time function. 

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 surrounding it ( is the sampling interval).
stf = gf.BoxcarSTF(5., center=0.)
Triangular STF¶
A triangular shaped source time function. It can be made asymmetric.
stf = gf.TriangularSTF(5., peak_ratio=0.5, center=0.)
Half sinusoid STF¶
A halfsinusoid source time function.
stf = gf.HalfSinusoidSTF(5., center=0.)
Resonator STF¶
stf = gf.ResonatorSTF(5., frequency=1.0)
Modelling targets¶
PyrockoGF 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 threecomponent 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 steplike 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, andthe
GNSSCampaignTarget
for e.g. steplike GPS displacements.
# 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 socalled lineofsight of the radar. See the forward modelling examples for detailed instructions of usage.
# east/north shifts as numpy.ndarrays in [m]
# lineofsight 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 PyrockoGF¶
Forward modelling, given a source and target description, is handled in the
socalled 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 steplike 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.
# 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:
# 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:
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()