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
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.
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,
Pyrocko-GF supports the simulation of various dislocation sources, focused on earthquake and volcano studies.
Multiple sources can be combined through the
For convenience, different parameterizations of seismological moment tensor point sources are available.
An isotrope moment tensor for explosions or volume changes.
Double force couple, for pure-shear 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 3-component single force point source.
Excess pore pressure point source.
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 Pyrocko-GF framework with
from pyrocko import gf
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)
A double-couple point source, describing shear ruptures.
dc_source = gf.DCSource(lat=54., lon=7., depth=5e3, strike=33., dip=20., rake=80.)
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)
A compensated linear vector dipole (CLVD) point source.
clvd_source = gf.CLVDSource( lat=48., lon=17., depth=5e3, dip=31., depth=5e3, azimuth=83.)
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.)
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)
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 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.
Boxcar shape source time function.
Triangular shape source time function.
Half sinusoid type source time function.
A simple resonator like source time function.
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.)
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 half-sinusoid source time function.
stf = gf.HalfSinusoidSTF(5., center=0.)
stf = gf.ResonatorSTF(5., frequency=1.0)
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.
In Pyrocko locations are given with five coordinates:
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
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
SatelliteTarget, for the forward modelling of InSAR data, and
GNSSCampaignTargetfor 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')
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')
GNSSCampaignTarget defines station locations and the
three components: east, north and up.
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
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
representing the resulting waveforms can be obtained from the engine’s
# 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() mod_scene.spool()