Moment tensor solutions for earthquakes in North America are routinely determined by Saint Louis University using the codes distributed Computer Programs in Seismology. The basic codes are then invoked using a sequence of shell scripts that acquire the waveforms, remove the instrument response, permit an interactive quality control of ground motions, invert for the source and document the processing and results in detail.
The purpose of this document is to provide the scripts, describe
their function, and provide detail on how they can be modified for
a local use. The result will be a web page such as https://www.eas.slu.edu/eqc/eqc_mt/MECH.NA/20220219003853/index.html.
This tutorial updates the scripts given in Regional Moment Tensor
Inversion.
Because the scripts are designed to support routine moment tensor inversion in a region, the scripts are connected in the sense that a master script invokes supporting scripts after first setting parameters for the supporting scripts.
These bash shell scripts are internally documented with comments and function well. Although one can create GUI's for processing, the value of these scripts is that they define the logical steps for processing. Many of the concepts have been incorporated into other systems.
In order to perform moment tensor inversion, it is necessary that the current distribution of Computer Programs in Seismology (CPS) be compiled and that the PATH variable be modified to point to the executables. In addition two environment parameters are used to point to the CPS directory and to the Green's functions directory. Since I use the bash shell on LINUX and MacOS, I have the following in the .profile file which is read upon login:
LINUX example 1
In ~/.profile
# set PATH so it includes user's PROGRAMS.330/bin if it exists
if [ -d "$HOME/PROGRAMS.330/bin" ] ; then
PATH=":.:$HOME/PROGRAMS.330/bin:$PATH"
fi
In ~/.bashrc
export CPS=${HOME}/PROGRAMS.330
export GREENDIR=/media/sf_rbh/PROGRAMS.310t/GREEN
------------------------
LINUX example 2
In ~/.profile
PATH=:.:$HOME/bin:$PATH:$HOME/PROGRAMS.310t/PROGRAMS.330/bin:
In ~/.bashrc
export GREENDIR='/d/rbh/GREEN'
export CPS='/home/rbh/PROGRAMS.310t/PROGRAMS.330/'
------------------------
OSX example
In ~/.profile
OPATH=$PATH
PATH=:.:$HOME/bin:$HOME/PROGRAMS.310t/PROGRAMS.330/bin:$OPATH
export GREENDIR=${HOME}/PROGRAMS.310t/GREEN
export CPS=${HOME}/PROGRAMS.310t/PROGRAMS.330
If everything is set up correctly, then the shell will be able to find the executables and also be able to access the directory contents:
$ which gsac /home/cps/PROGRAMS.330/bin/gsac $ ls -F $CPS bin/ C.proto* DOC/ readline-8.1/ VOLII/ VOLV/ VOLX/ C* CPSInstall.pdf include/ Setup* VOLIII/ VOLVI/ XVIG/ CALPLOT/ Csrconly* IRIS/ SUBS/ VOLIV/ VOLVII/ C.out dialog-1.1-20080819/ lib/ VOLI/ VOLIX/ VOLVIII/ $ ls -F $GREENDIR CUS.REG/ Models/ WUS.REG/
Here the variable CPS is the location of the CPS software distribution while the variable GREENDIR is the location of the Green's functions required for inversion. These can be anywhere on the disk system, but they must exist.
In the distribution where will be a src directory that contains some utility programs for computing ML and mLg magnitudes from the data. These are provided since these are usful parameters to define the size of an earthquake. They really are part of the package because of a research interest in the relationship between a local magnitude and Mw. It will be necessary to compile these.
In addition the scripts make use of the convert program of the ImageMagick package to convert EPS graphics to PNG for inclusion in the web page created.
Finally, the Generic Mapping Tools (GMT) are use to make some maps
All of the scripts are included in the file MECH.tgz Click on this link to download this file. Then do the following steps:
gunzip -c MECH.tgz | tar xf - cd MECH rm -fr bin mkdir bin cd src make all cd ..
The programs created and their purposes are as follow:
mlmag - pass ground velocities through a Wood-Anderson response and compute ML using the IASPEI formula
mblgmag - pass ground velocities through a WWSSN Short Period response and compute mLg using the IASPEI formula
linreg1 - given an input stream of (x,y) pairs, perform a linear regression of the form y = A + B x. The command line flag -1 uses the residuals, res, as a 1/res weight to implement an L1 inversion_
linreg1w - given an input stream of (x,y,wt) triplets, where wt is a weight, perform a linear regression of the form y = A + B x. The command line flag -1 uses the residuals, res, as a 1/residual weight to implement an L1 inversion
mlwt - given (magnitude,distance) pairs, define a weight that emphasizes observations at distances less than 100 km. the output stream sonsists of (x,y,wt) triplets. The ML scripts will attempt to correct for a difference between regional attenuation and that implicit in the IASPEI ML formula that was defined for Southern California (USA).
stationxml-seed-converter-2.1.0.jar - This is from https://github.com/iris-edu/stationxml-seed-converter. This program will convert the downloaded station xml file from IRIS to datalessSEED than can be used with rdseed and the downloaded miniSEED to create the RESP, pole-zero and Sac files.
FetchData - an IRIS tool (https://service.iris.edu/clients/) that fetches time series data (miniSEED), simple metadata (ASCII) and instrument responses (SEED, RESP and SAC PZs). This can be used to access data from IRIS and cooperating data centers.
fdsnws_fetch - (https://geofon.gfz-potsdam.de/software/fdsnws_fetch/) THIS CODE IS NOT IN THE SOURCE DIRECTORY. IT MUST BE SEPARATELY DOWNLOADED.
To do this use the python pip command which may have to be installed. On LINUX Mint one would fist install pip (tecmint.com/install-pip-in-linux a)s
Install PIP On Debian/Ubuntu/Mint: # apt install python-pip #python 2 # apt install python3-pip #python 3The link gives the commands for other LINUX distributions. Next install the FDSN scripts from pypi.org/project/fdsnwsscripts:
pip oinstall fdsnwsscriptsThe installation will place the scripts into the ${HOME}/.local/bin directory and warn is this is not in the $PATH If not than place the following line in the bash .profile file:
# set PATH to include the .local/bin directory export PATH="$HOME/.local.bin"
Note that Version 2022.17 fdsnws2seed will create a SEED volume which can be used with rdseed but also creates Sac files and Sac pole-zero files.
mseed2sac - an IRIS tool (https://ds.iris.edu/ds/nodes/dmc/software/downloads/mseed2sac/)
Converts miniSEED to SAC format. An optional metadata file may
be supplied that contains values for the SACheaders that are
not available in miniSEED.
query or java -jar CWBQuery.jar $* (https://github.com/usgs/edgecwb/wiki/CWBQuery). This accesses the USGS public server and can provide waveforms in SAC format and responses in SAC pole-zero format. THIS CODE IS NOT IN THE SOURCE DIRECTORY.
The MECH directory provides the required scripts and prototype organization. In my usage, there are similar directories that focus on different geographical reqions. For example I have the organization
MOMENT_TENSOR
|
|----- MECH.NA [Moment tensor inversion for earthquakes in North America]
|
|----- MECH.IT [Moment tensor inversion for earthquakes in Italy]
|
|----- MECH.EU [Moment tensor inversion for earthquakes in other locations in Europe]
Thus the MECH directory is a prototype for any region.
The MECH directory contains the following directories
MECH
├── 0XXXREG
│ ├── DAT.REG
│ │ └── NOUSE
│ ├── GRD.REG
│ ├── HTML.REG
│ ├── MLG.REG
│ ├── ML.REG
│ ├── MTD.REG
│ ├── MTGRD.REG
│ ├── MTGRD.REG.DC
│ ├── MTGRD.REG.DEV
│ ├── MT.OTHER
│ ├── MT.REG
│ ├── NEW2.REG
│ ├── SYN.REG
│ └── TMLG.REG
├── bin
├── PROTO.CWB
├── PROTO.GEOFON
├── PROTO.I
├── PROTO.WS
└── src
└── mseed2sac-2.3
├── doc
├── libmseed
│ ├── doc
│ ├── example
│ ├── matlab
│ └── test
│ └── data
└── src
p>and scripts starting with DODO DOGEOFON DOQUERY DOSOLUTION DOWS5 DOWSSETUP DOCWB DOGEOFON5 DOSETUP DOSOLUTION5 DOWSF DOCWB5 DOGEOFONSETUP DOSETUP5 DOWS DOWSF5
Typical source inversion requires broadband waveform data for a given event. Thus the process starts by specifying the event coordinates. In an automated or GUI based system this information would be provided by the system. For command line inversion, one must edit an initial file, here called DO, which looks like the following:
#!/bin/sh ##### # valid regions # REG Region VELOCITY_MODEL # CUS Central US CUS # WUS Wsstern US WUS ##### # Command syntax: #DO YEAR MO DY HR MN SC MSC LAT LON DEP MAG REG NEIC FELTID STATE/COUNTRY ##### # DO YEAR MO DY HR MN SC MSC LAT LON DEP MAG REG NEIC FELTID STATE/COUNTRY # 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 DOCWB "2022" "05" "17" "08" "07" "34" "000" "48.295" "-121.869" "3.6" " 3.6" "WUS" "NONE" "NONE" "Washington" DOWS "2022" "05" "17" "08" "07" "34" "000" "48.295" "-121.869" "3.6" " 3.6" "WUS" "NONE" "NONE" "Washington" DOWSF "2022" "05" "17" "08" "07" "34" "000" "48.295" "-121.869" "3.6" " 3.6" "WUS" "NONE" "NONE" "Washington" DOGEOFON "2022" "04" "24" "04" "27" "55" "000" "43.05 " "18.15" "19.0" " 4.8" "WUS" "NONE" "NONE" "Bos-Herz" DOCWB5 "2022" "05" "17" "08" "07" "34" "000" "48.295" "-121.869" "3.6" " 3.6" "WUS" "NONE" "NONE" "Washington" DOWS5 "2022" "05" "17" "08" "07" "34" "000" "48.295" "-121.869" "3.6" " 3.6" "WUS" "NONE" "NONE" "Washington" DOWSF5 "2022" "05" "17" "08" "07" "34" "000" "48.295" "-121.869" "3.6" " 3.6" "WUS" "NONE" "NONE" "Washington" DOGEOFON5 "2022" "04" "24" "04" "27" "55" "000" "43.05 " "18.15" "19.0" " 4.8" "WUS" "NONE" "NONE" "Bos-Herz" #NOTE IN AN ACTUAL SCRIPT ONLY 1 WILL NOT HAVE A # COMMENT IN FRONT. # YOU WILL ONLY WANT TO PROCESS ONE EVENT AT A TIME # The use of "" pairs is a habit to insure that the SHELL passes just obne variable for "Noth Dakota" # in stead of two as in North Dakota. The various scripts DCWB check command syntax for exactly 15 arguments.
In this script 15 variables must be specified. These are the year, month, day, hour, minute, second and millisecond of the event origina time. The event latitude, longitude and depth. The event magnitude. The region which is used to select the particular velocity model used to create the Green's functions. Finally there are 3 entries that are used to annotate the final documentation.
Note that the lines are very similar. They differ only in the first command DOCWB, DOWS and DOWSF, etc. These scripts are designed to work with different data sets, e.g., the USGS Continuous Wave buffer, waveform data deposited at IRIS and waveform data deposited at federated data centers. There is also a DOFDSN to obtain waveform data from GFZ. The DOCWB5, DOWS5 and DOWSF5 perform the same functions, but use GMT v5 and greater rather than GMT4 or GMT3.
These scripts do the following:
create a durectory for the event
create a subdirectory structure for the inversion
create a subdirectory structure for downloading waveforms, deconvolving and QC'ing the waveforms.
download the waveform and responses
deconvolve waveforms to make ground velocity in m/s
rotate the wave forms
perform interactive QC
perform the inversion
THE SIMPLEST WAY TO SET UP THE DATA STRUCTURE FOR AN EVENT WITHOUT
ACTUALLY DOWNLOADING WAVEFORMS IS TO RUN
DOWS or DOWS5. If the data are available by
another method, place the sac files and the response files in the EVENT/EVENT/Sac
directory. If the RESP file is available, then run IDOEVT in the
parent directory. If the response is in terms of pole-zero files,
modify the IDOEVT script to identify them and then to use the proper
gsac syntax for deconvolution.
The following links discuss each in detail. Note the differences between DOCWB and DOCWB5 are minor, and the minor differences are indicated.
DOCWB -
This scripts runs to completion with only the interactive QC
oepration.
DOWS
- These scripts create a DOFINISH file that is run
manually to acquire the data and then runs the deconvolution,
rotation, QC and finally the inversion. The reason for the
DOFINISH script is to permit fine tuning. This is especially
the case for acquiring GEOFON data using fdsnws_fetch
which does not support a distance search. Instead of obtaining
all waveform data, this DOFINISH is edited to select the
wavefrms from desired networks.
Each if these scripts invoke DOSETUP, a script to get data using
CWB (DOCWB) or create a script named DOFINISH (DOWS and DOWSF),
and use a script named
DOSOLUTION (or DOSOLUTION5 when using GMT v5 or greater).
DOSETUP -This sets up the directory structure for the event. In addition it associates the REG parameter of the DO script with the model file name and the Green's function directory by defining the shell variables GMODEL and VMODEL which are then placed into the HTML.REG/DOHTML, GRD.REG/DOSTA and EVENT/IDOROT scripts, for example. The case statment in this script has a lot of information that can be used later. For example, the line
CUS) GMODEL="CUS.REG";VMODEL="CUS";FIDREG="cus";REGIONTEXT="Central and Southeastern US" ;;
indicates that the directory containing the Green's functions will be ${GREENDIR}/CUS.REG and that the velocity model will be CUS.mod which will be in ${GREENDIR}/Models.
The script also places the event information in the HTML.REG/DOHTML and EVENTID/IDOEVT scripts after the EVENTID directory is created and the prototypes moved there.
The script may seem complicated, but the requirement is that all moment tensors be determined using the same procedure. Depending on the event location, it may be necessary to use different velocity models.
(data acquisition and preprocessing of waveforms - see below for DOQUERY and DOFINISH)
The data acquisition, deconvolution and rotation is done by the script DOQUERY for DOCWB or DOCWB5, and DOFINISH, which is created by DOWS
DOQUERY -
DOFINISH scripts - created by DOWS, DOWSF or DOFDSN
For routine moment tensor solutions for earthquakes, a grid search over all possible double0couple solutions is performed in the GRD.REG directory. If DIR is the event identifier of the form YYYMMDDHHMMSS, then the DOSOLUTION is just
#####
# Waveform inversion
#####
cd ${MYPWD}/${DIR}/GRD.REG
DOGRD
DODELAY
DOPLTSAC
DOCLEANUP
#####
# in the future do the documentation here
#####
cd ${MYPWD}/${DIR}/HTML.REG
DOHTML
DOSOLUTION5 differs in that DOHTML5 is invoked at the end rather than DOHTML.
DOGRD performs the grid search over a range of source
depths. It used the supporting script DOSTA to find the
Green's functions at the observation distance, and then filters
both. The filtering corner frequencies can be changed by editing
the DOSTA script. In addition a microseism band reject
filter can be applied and the signal time window can be adjusted.
DODELAY uses the timeshifts required to fit the Green's
functions to the observed data to estimate changes in event
location. If there seems to be a large mislocation, perhaps the
Green's function model is not correct or the event location should
be reviewed.
DOPLTSAC (DOPLTSAC5) used the ImageMagick convert
program to make PNG files of the goodness of fit and the waveforms
comparison for the web page.
DOCLEANUP removes temporary files to clean up the
directory to make the archive smaller.
DOHTML (DOHTML5) creates the web page. If there are other moment tensor solutions described in ${MYPWD}/${DIR}/MT.OTHER, these are included in the event web page for completeness. The DOHTMLK can be hand edited to add additional information, such as relocation results.
An example of the HTML created for an event is given in 20201108141007/HTML.REG/index.html and the files in the MT.OTHER directory are in 20201108141007/MT.OTHER to provide an example of how the HTML is written to permit the comparison.
The event directory will also have the subdirectories MTGRD.REG,
MTGRD.REG.DEV and MTGRD.REG.DC. These codes are based
on a paper my Zhu and Ben-Zion
(Zhu, L., and Y. Ben-Zion (2013). Paramaterization of
general seismic potency and moment tensors for source
inversion of seismic wavefrom data, Geophys. J. Int. 194,
839-843. doi: 10.1093/gji/ggt137) that describes how to
perform a grid search for a full moment tensor. The DOMTGRD
in MTGRD.REG searchs for all six moment tensor parameters. The
DOMTGRD in MTGRD.REG.DEV executes wvfmtgrd96 -DEV to
perform a search for the deviatoric moment tensor. Finally the
DOMTGRD in MTGRD.REG.DC searches for the best double couple using
the command wvfmtgrd96 -DC. The wvfmtgrd96 -DC is
equivalent to the CPS wvfgrd96. The wvfmtgrd96 -DEV takes longer
to run since the grid search is over four parameters. The
wvfmtgrd96 searches over 5 parameters and take a very long while
to run.
Within each of these directories, the scripts are run int he following order:
DOMTGRD DODELAY DOPLTSAC or DOPLTSAC5 DOCLEANUP
If these inversions are run, then the DOHTML/DOHTML5 script in HTML.REG will include the results of this proceesing in the HTML.REG/index.html file.
Moment tensor inversion requires good signal-to-noise
ratios and a simple signal. Thus the frequency band and time
window may have to be adjusted.
Other processing scripts may also need to be adjusted. The
following discussion will discuss these scripts. The scripts for
data acquisition, deconvolution and quality control will be
discussed first. Then the scripts for moment tensor inversion will
be discussed. Finally the scripts for creating the HTML
documentation will be discussed last. To make this
discussion easier to follow, the path name consisting of the
origin time will be given. Thus the script location may be in the
directories
20220517080734/20220517080734/Orig 20220517080734/GRD.REG
case ${KCMPNM} in
BHZ) C="BH" ;;
HHZ) C="HH" ;;
HNZ) C="HN" ;;
LHZ) C="LH" ;;
esac
If other channels are used, e.g., EH, add an entry EHZ) C="EH" ;;
fileid list fname dist az concat on format colon
markt on
xlim vel 3.3 -40 vel 3.3 70
#xlim a -10 a 180
r *
cut off
sort up dist
rtr
taper w 0.01
hp c 0.03 n 3
lp c 0.10 n 3
#br c 0.12 0.25 n 4 p 2
qdp 10
ppk q relative perplot 3
wh IHDR20
q
The window and filter settings are the same as those in the DOSTA scripts. I often change the frequency band and decide whether to use a microseism band reject filter. If any of these are changes, then the corresponding values in the DOSTA scripts for the inversion should also be change.
If channels
other than BH, HN or HH are used, then modify the lines
###### define the components
case ${KCMPNM} in
BHZ|HHZ|HNZ) COMP=1 ; NAM=Z;;
BHR|HHR|HNR) COMP=2 ; NAM=R;;
BHT|HHT|HNT) COMP=3 ; NAM=T;;
*) break;;
esac
to have an entry
like BHZ|HHZ|HNZ|EHZ) COMP=1 ; NAM=Z;;. This part of the
code sets up the input for the wvfgrd96 and the other
inversion codes.
An example of running all inversion is given at https://www.eas.slu.edu/eqc/eqc_mt/MECH.NA/20220601150149/index.html
. This was an event in West Texas.
In this case the wvfmtd96 and wvfmtgrd96 -DEV agreed
as did the wvfmt96 and wvfmtgrd96 solutions. As
mentioned the wvfmt96 and wvmtd96 do not have a
robust way to handle time shifts. However this event was well
located and the Green's function model was appropriate and large
time shifts were not required.
The script used iwith GMT5+ after the QC'd waveforms were placed
in the DAT.REG directory is the following:
(cd GRD.REG ;DOGRD ;DODELAY;DOPLTSAC5)
(cd MT.REG ;DOMT ;DODELAY;DOPLTSAC5)
(cd MTD.REG ;DOMTD ;DODLEAY;DOPLTSAC5)
(cd MTGRD.REG.DC ;DOMTGRD;DODELAY;DOPLTSAC5)
(cd MTGRD.REG.DEV;DOMTGRD;DODELAY;DOPLTSAC5)
(cd MTGRD.REG ;DOMTGRD;DODELAY;DOPLTSAC5)
(cd HTML.REG;DOHTML5)
Normally for earthquakes, one would just run the commands
(cd GRD.REG;DOGRD;DODELAY;DOPLTSAC)
(cd HTML.REG;DOHTML5)