Content

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

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:

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

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

```
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

```
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 double-couple 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,
magnitude=6.3)
# 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 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 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.)
```

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

```
# 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.

```
# 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.

```
# 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()
```