Sign convention for ZRT seismograms

[liamtoney]

Hi Lion, I'm using Syngine to compute force source Green's functions. I'm requesting displacement seismograms in ZRT. I was assuming the convention to be:

• Z positive up
• R positive away from source
• T positive 90° clockwise from the line pointing from source to receiver

I built a sandbox to test this by applying unit forces in the up, south, and east directions at a source due west of a receiver:

Given my sign convention above, the first two waveform panels make sense. But in the last one, there's a negative first motion for a force in the direction from source to receiver. This suggests that radial is defined to be positive towards the source. Am I missing something, is this in fact the convention, or is this perhaps a bug?

Please click here for the code I used to produce the plots.

import matplotlib.pyplot as plt
import obspy
from obspy import UTCDateTime

BASE_URL = 'https://service.iris.edu/irisws/syngine/1/query?'
MODEL = 'iasp91_2s'
T0 = -20  # [s]
GF_DURATION = 40  # [s]


def get_syngine_waveforms(
    station_lat, station_lon, source_lat, source_lon, source_depth, origin_time, forces,
):
    url = BASE_URL + '&'.join(
        [
            'format=miniseed',
            'components=ZRT',
            'units=displacement',
            'model=' + MODEL,
            'origintime=' + origin_time.format_iris_web_service(),
            'starttime=' + (origin_time + T0).format_iris_web_service(),
            'endtime=' + (origin_time + GF_DURATION).format_iris_web_service(),
            'receiverlatitude=' + str(station_lat),
            'receiverlongitude=' + str(station_lon),
            'sourcelatitude=' + str(source_lat),
            'sourcelongitude=' + str(source_lon),
            'sourcedepthinmeters=' + str(source_depth),
            'sourceforce=' + ','.join([str(force) for force in forces]),
        ]
    )

    return obspy.read(url)


STATION_LAT = 0  # [deg.]
STATION_LON = 0.01  # [deg.]

SOURCE_LAT = 0  # [deg.]
SOURCE_LON = 0  # [deg.]

SOURCE_DEPTH = 0  # [m]

# [Newtons] f_r, f_theta, f_phi (r is positive upwards, theta is positive to the south,
# and phi is positive to the east)
UP = (1, 0, 0)
SOUTH = (0, 1, 0)
EAST = (0, 0, 1)

# Plot "station map"
fig, ax = plt.subplots()
ax.scatter(
    SOURCE_LON, SOURCE_LAT, color='green', label='Source',
)
ax.scatter(
    STATION_LON, STATION_LAT, color='red', label='Station',
)
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')
ax.legend()
ax.set_title(f'Source depth = {SOURCE_DEPTH:g} m')
fig.tight_layout()
fig.show()

# Grab and plot waveforms
for forces in UP, SOUTH, EAST:
    st = get_syngine_waveforms(
        station_lat=STATION_LAT,
        station_lon=STATION_LON,
        source_lat=SOURCE_LAT,
        source_lon=SOURCE_LON,
        source_depth=SOURCE_DEPTH,
        origin_time=UTCDateTime(2016, 5, 22, 7, 57, 34),
        forces=forces,
    )
    fig = plt.figure()
    st.plot(fig=fig, equal_scale=True)
    fig.suptitle(f'$f_r$, $f_\\theta$, $f_\\phi$ = {forces}')
    fig.tight_layout()
    fig.show()

[krischer]

I'll defer this @martinvandriel because he wrote most of the rotation logic.

[martinvandriel]

Thanks for reporting, I guess you are right that there is some issue here.

T positive 90° clockwise from the line pointing from source to receiver

Is this a general convention? I guess I would have assumed ZRT to be right handed. In that case the theta trace in your figure would also be of wrong sign. Just to clarify.

[liamtoney]

You are right that the convention I wrote down above is indeed not right-handed. But I also may be wrong... for example, the IRIS rotation docs here seem to show a right-handed coordinate system with:

• Z positive down
• R positive away from source
• T positive 90° clockwise from the line pointing from source to receiver

As currently implemented, the Instaseis sign convention is also right-handed, which is good — it just might not be the right-handed convention that people are expecting. I'd imagine they'd be expecting what's shown in the image above, i.e. what's on IRIS. So perhaps this could just be resolved with some documentation clarifications?

[qnishida]

I was also confused about the definition. On the same page, the below figure shows Z positive upward.

In the case of the moment tensor, the polarization of the Z component was consistent with the upward positive. I compared the moment tensor response with the Mineos (explosion source, gcarc = 80 degrees). The first arrival shows the teleseismic P wave.

On the other hand, in the case of a vertical single force (gcarc = 80 degrees), the radial component of P wave was flipped.

The coordination should be consistent for these two cases.

If these differences are specifications, clear documentation would be appreciated.

[claudiodsf]

Hi, I also stumbled upon this problem while working with single forces.

My conclusion is that Syngine produces correct ZNE traces but wrong ZRT results, with radial component flipped.

Click here to show a modified version of @liamtoney code, for the pure vertical force (1, 0, 0) at a station to the east.

#!/usr/bin/env python
import matplotlib.pyplot as plt
from obspy import UTCDateTime
from obspy.geodetics import gps2dist_azimuth
from obspy.clients.syngine import Client
client = Client()

SOURCE_LATITUDE = 0.
SOURCE_LONGITUDE = 0.
SOURCE_DEPTH_IN_KM = 0.
STATION_LATITUDE = 0.
STATION_LONGITUDE = 0.01

FORCE_R_IN_N = 1e10
FORCE_THETA_IN_N = 0
FORCE_PHI_IN_N = 0

MODEL = 'iasp91_2s'
COMPONENTS = 'ZRT'

ORIGIN_TIME = UTCDateTime(2016, 5, 22, 7, 57, 34)
T0 = -20  # seconds
GF_DURATION = 40  # seconds

COMPONENTS = 'ZNE'
st_zne = client.get_waveforms(
    model=MODEL,
    sourcelatitude=SOURCE_LATITUDE, sourcelongitude=SOURCE_LONGITUDE,
    sourcedepthinmeters=SOURCE_DEPTH_IN_KM*1e3,
    sourceforce=[FORCE_R_IN_N, FORCE_THETA_IN_N, FORCE_PHI_IN_N],
    receiverlatitude=STATION_LATITUDE,
    receiverlongitude=STATION_LONGITUDE,
    origintime=ORIGIN_TIME,
    starttime=ORIGIN_TIME+T0,
    endtime=ORIGIN_TIME+T0+GF_DURATION,
    units='displacement',
    components=COMPONENTS
)
COMPONENTS = 'ZRT'
st_zrt = client.get_waveforms(
    model=MODEL,
    sourcelatitude=SOURCE_LATITUDE, sourcelongitude=SOURCE_LONGITUDE,
    sourcedepthinmeters=SOURCE_DEPTH_IN_KM*1e3,
    sourceforce=[FORCE_R_IN_N, FORCE_THETA_IN_N, FORCE_PHI_IN_N],
    receiverlatitude=STATION_LATITUDE,
    receiverlongitude=STATION_LONGITUDE,
    origintime=ORIGIN_TIME,
    starttime=ORIGIN_TIME+T0,
    endtime=ORIGIN_TIME+T0+GF_DURATION,
    units='displacement',
    components=COMPONENTS
)

# rotate to ZRT using obspy
_dist, az, baz = gps2dist_azimuth(
    SOURCE_LATITUDE, SOURCE_LONGITUDE, STATION_LATITUDE, STATION_LONGITUDE
)
print(f'azimuth: {az}, back-azimuth: {baz}')
st_zrt_obspy = st_zne.copy().rotate(method='NE->RT', back_azimuth=baz)

for st, title in zip(
    [st_zne, st_zrt, st_zrt_obspy], ['ZNE', 'ZRT_syngine', 'ZRT_obspy']
):
    fig = plt.figure()
    st.plot(fig=fig)
    fig.suptitle(f'{title}, back azimuth: {baz}')
    fig.tight_layout()
plt.show()

The ground motion recorded at this station should be upwards (positive Z) and eastward (positive E), which is the case, looking at the "ZNE" traces produced by Syngine.

In the "ZRT" frame, the same motion should be away from the source, so: positive Z and positive R.

Now, Syngine gives negative R (wrong), while rotating the traces with ObsPy, one gets positive R (correct).

Note that if you rotate in ObsPy using azimuth instead of back azimuth, you get the same result of Syngine.
Could that be the problem?

[liamtoney]

I just updated the minimal working example at the start of this issue, given here:

▶ Please click here for the code I used to produce the plots.

So that it still runs. The plots produced are the same, so I didn't update them.

[liamtoney]

Thurin et al. (2023) is an erratum which cites this issue:

The Instaseis GitHub repository (https://github.com/krischer/instaseis/) was established on 21 August 2014. The error in the radial‐component force‐generated synthetic seismograms was identified in the Instaseis GitHub repository on 7 July 2020 (#77). The erroneous seismograms were used within part of our study in 2022, and the issue has not been fixed in Instaseis as of today (6 April 2023). We made the correction in our moment tensor (and force) estimation code (MTUQ: https://github.com/uafgeotools/mtuq) on 25 January 2023, such that the MTUQ will flip the radial component of force‐generated Instaseis synthetic seismograms, thereby correcting the Instaseis issue. (A similar correction is applied within the landslide‐modeling code — lsforce (Toney and Allstadt, 2021), which also access force generated synthetic seismograms via Syngine and Instaseis.) Any previous or ongoing study that uses AxiSEM force (not moment tensor) synthetics, via Instaseis or Syngine, should re‐examine their analysis to see if it was impacted by flipped‐sign radial component synthetic seismograms. This includes the downloadable databases available at http://ds.iris.edu/ds/products/syngine/.

References

Thurin, J., Tape, C., & Modrak, R. (2023). Erratum to Multi‐Event Explosive Seismic Source for the 2022 $M_\mathrm{w}$ 6.3 Hunga Tonga Submarine Volcanic Eruption. The Seismic Record, 3(2), 168–170. https://doi.org/10.1785/0320230014