New route in instaseis server to return full Green's tensor

[martinvandriel]

This is needed to couple instaseis to seiscomp which we will use for a new Mars waveform exploration tool. As a nice side effect, this will directly allow to use instaseis with the TDMT (time domain moment tensor inversion) tool in seiscomp, which has a possible broad user base especially if IRIS opens this route as well.

[krischer]

Figure out which definition of the Green's tensor should be returned.

[martinvandriel]

[martinvandriel]

@gempa-jabe @jfclinton We are struggling to understand what components we should be returning to be compatible with TDMT and don't understand the definition above. Should we be returning the ZSS, ZDS, ..., TDS functions? How are these fundamental faults defined in terms of strike, dip, rake or moment tensor?

[gempa-jabe]

How are these fundamental faults defined in terms of strike, dip, rake or moment tensor?

Green's functions are a set of pre computed waveforms for three fundamental faults to represent a deviatoric point source: vertical strike-slip (SS), vertical dip-slip (DS) and a dip-slip with 45° of dip (DD).

Each Green's function consists of 8 components: ZSS, ZDS, ZDD, RSS, RDS, RDD, TSS and TDS where Z,R and T refer to the vertical, radial and tangential components.

Does that help already? I can also show you how the components are combined to form synthetics with a given moment tensor.

[gempa-jabe]

Here is some pseudo code (Python) that creates synthetics for given azimuth and strike, dip and rake.

rdip  = deg2rad(dip)
rstr  = deg2rad(strike)
rrake = deg2rad(rake)

nx = -math.sin(rdip)*math.sin(rstr)
ny =  math.sin(rdip)*math.cos(rstr)
nz = -math.cos(rdip)

dx =  math.cos(rrake)*math.cos(rstr)+math.cos(rdip)*math.sin(rrake)*math.sin(rstr)
dy =  math.cos(rrake)*math.sin(rstr)-math.cos(rdip)*math.sin(rrake)*math.cos(rstr)
dz = -math.sin(rdip) *math.sin(rrake)

Mxx = 2*dx*nx
Mxy = dx*ny+dy*nx
Mxz = dx*nz+dz*nx
Myy = 2*dy*ny
Myz = dy*nz+dz*ny
Mzz = 2*dz*nz

razi = deg2rad(azi)
sinazi = math.sin(razi)
cosazi = math.cos(razi)
sin2azi = math.sin(2*razi)
cos2azi = math.cos(2*razi)

# zss, zds, zdd, rss, rds, rdd, tss, tds are arrays
# len is the number of samples

for i in xrange(len):
    z = 0
    r = 0
    t = 0

    # Mxx
    z += Mxx *  0.5*(cos2azi*zss[i] - zdd[i])
    r += Mxx *  0.5*(cos2azi*rss[i] - rdd[i])
    t += Mxx * -0.5*sin2azi*tss[i]

    # Myy
    z += Myy * -0.5*(zdd[i] + cos2azi*zss[i]);
    r += Myy * -0.5*(rdd[i] + cos2azi*rss[i]);
    t += Myy *  0.5*sin2azi*tss[i];

    # Mxy
    z += Mxy * sin2azi*zss[i];
    r += Mxy * sin2azi*rss[i];
    t += Mxy * cos2azi*tss[i];

    # Mxz
    z += Mxz * -cosazi*zds[i];
    r += Mxz * -cosazi*rds[i];
    t += Mxz *  sinazi*tds[i];

    # Myz
    z += Myz * -sinazi*zds[i];
    r += Myz * -sinazi*rds[i];
    t += Myz * -cosazi*tds[i];

    print z,r,t

[martinvandriel]

Thanks Jan for the quick answer. Let us first clarify, where the interface is: I assume, we pass you a mseed file / stream containing the components ZSS, ZDS, ZDD, RSS, RDS, RDD, TSS and TDS. So what you discribe above would be done on your side and the more important issue for us is how to compute these 8 components.

This raises the question of how the fundamental faults are defined: 'Vertical strike slip' tells the dip and the rake, but not the strike (probably relative to the azimuth towards the receiver to define it independently of the receiver). Similar for the other sources. I could maybe back-engineer from your code above, but would prefer to have a proper definition.

Our idea would be to just call the get_seismograms routine of instaseis for the three fundamental faults and then assemble the data as needed. As the data needed will be buffered on the first call, it should only be marginally slower than the get_seismograms for a single source.

[gempa-dirk]

We can assume Green's functions that are in SAC. Below is an example of their generation using qseis, by Ronjiang Wang, GFZ (http://www.gfz-potsdam.de/sektion/erdbeben-und-vulkanphysik/services/downloads-software/).

The coordinate frame of qseis is:

#                 north(x)
#                  /
#                 /\ strike
#                *----------------------->  east(y)
#                |\                       \
#                |-\                       \
#                |  \     fault plane       \
#                |90 \                       \
#                |-dip\                       \
#                |     \                       \
#                |      \                       \
#           downward(z)  \-----------------------\

Assuming zero azimuth, the x and y components correspond to the radial component pointing away from the source and the transverse components, respectively.

#Mxx   Myy   Mzz   Mxy   Myz   Mzx     file 
 0     0     0     1.0   0     0      'seis_m1'
 1.0  -1.0   0     0     0     0      'seis_m2'
 0     0     0     0     1.0   0      'seis_m3' 
 0     0     0     0     0     1.0    'seis_m4'
 1.0   1.0   1.0   0     0     0      'seis_m6'
-1.0  -1.0   2.0   0     0     0      'seis_cl'

For each moment tensor qseis delivers the files seis_xx.tz, seis_xx.tr and seis_xx.tt.

Association of qseis files with files in the Green's function database:

seis_m1.tt  ->      depth.distance.TSS

seis_m2.tz  ->      depth.distance.ZSS
seis_m2.tr  -> -1 * depth.distance.RSS

seis_m3.tt  ->      depth.distance.TDS

seis_m4.tz  ->      depth.distance.ZDS
seis_m4.tr  -> -1 * depth.distance.RDS

seis_cl.tz  ->      depth.distance.ZDD
seis_cl.tr  -> -1 * depth.distance.RDD

seis_m6.tz  ->      depth.distance.ZEP
seis_m6.tr  -> -1 * depth.distance.REP

The Green's function files are stored in the repository:

$GREENS_FUNCTION_HOME/model/depth/distance/

The depths and distances are given in multiples of 100 m and 1 km, respectively. The formats of the depth and the distance strings are:

depth:    %4.4i
distance: %5.5i

model can be any string describing uniquely the used model, e.g.

prem

A description file containing the distance and depth descriptions must exist in the repository:

$GREENS_FUNCTION_HOME/model.desc

The descriptor file contains the ranges and the intervals of the Green's functions:

# depth [from] [to] [step]
depth 10 30 2
depth 30 700 10
# distance [from] [to] [step]
distance 10 10000 10

[martinvandriel]

Thanks Dirk, I think this should contain all necessary information.

[martinvandriel]

I added a function get_greens_seiscomp according to the definitions above. I noted, that the sign of the radial component seems to be different to Minson & Dreger, but otherwise, the test uses eq 6-8 from that paper and works out smoothly.

Remaining work: add a route to the server that returns these Green's function. @krischer

[martinvandriel]

So we have a http route now as well, mainly based on the seismograms routine with some features removed. Testing is basic, but as tests are extensive in the seismograms route that maybe does not matter to much.

Remaining work: documentation

[martinvandriel]

After some discussion with @gempa-jabe decided to switch the radial sign to be consistent with Minson & Dreger. TODO!

[martinvandriel]

Done in 058beb014bfdf96dee43c46737ea9899854b3eea

Can we close this issue?

[gempa-jabe]

I am fine with the current implementation so far. Thanks a lot.