Practical guide for developing
1D seismic velocity models from observed 3-component seismograms of first-arriving
teleseismic P waves using receiver function analysis and synthetic
The guide below consists of three independent sections:
a code rfgen_asc previously distributed here has been
superceeded by the code rfsyn which uses the method of spectral smoothing
identical to that in recfunk code family.
Primary codes were written by Jeffrey Park, with modifications and
service routines by Vadim Levin and Nikolai Shapiro.
ALL CODES ARE FREE , except for SAC analysis
package which is copyrighted.
However, only SAC format is used in the computations. As long as you can manipulate SAC headers yourself - no proprietary software would be needed.
All operations are performed on a UNIX workstation, either SUN or Linux PC.
Compilers for C and FORTRAN will be needed.
I/O is in the form of ASCII timeseries, and also in binary SAC format.
Images are in PostScript (thus you'll
need a viewer...)
recfunk.ascii: a FORTRAN code that will transform a collection of 3 component records
(vertical, radial, transverse) into ASCII tables of 2 component RFs (radial and transverse),
binned with backazimuth or epicentral distance.
A rudimentary manual for recfunk code is here . The code recfunk.ascii differs only in that it doesn't generate graphical screen output,
and thus may be compiled without OpenLook libraries and PLOTIT package.
recfunk_average.ascii: a FORTRAN code that will transform a collection of 3 component records
(vertical, radial, transverse) into averaged 2 component RFs (radial and transverse).
The code takes same input as the recfunk above.
It will average all
observations, and will output two ASCII files, for radial and transverse receiver functions. Files contain
two-column tables (time, RF_value).
rfsyn: a FORTRAN code computing a receiver-function response of
a stack of anisotropic layers to a plane wave incident from below via a
instructions for use are here
jp2tab C-language model format converter.
mk_macro C-language program that helps prepare data for >recfunk, courtesy of Nikolai Shapiro (CU Boulder).
a GNU ZIP-compressed TAR file with source codes and
SUN Solaris executables is here .
On LDEO or Yale SUN networks include ~vadim/bin in your path and the codes will
GMT graphics package. See GMT website for instructions
on setting it up.
A procedure described below does the trick or getting data ready.
It will be easy to follow on a SUN platform. Not tested under Linux...
Data for the analysis described below are 3 components of motion associated with
teleseimic compressional (P) waves.
A "receiver function" is defined here as a timeseries of shear motion within the coda
of the teleseismic P wave. The main assumption is that the
vertical component contains primarily compressional motion, and thus deconvolving it
from the horizontal components removes the signature of the earthquake source, plus all
compressional reverberations, and leaves only shear-polarized motion.
- Step 1. Obtaining data.
Most likely source of data is IRIS DMC archive. The tool developed
by DMC staff, called WEED, may be uploaded from their website . The tool allows the user to specify time windows of data (e.g., 200-sec window centered on a time of the first P wave arrival) for specific station-earthquake combinations, to format the
request to the DMC, and to unpack the SEED volume the archive
will send in responce.
- Step 2. Getting data ready.
To use the SAC macro generator included in this package (assuming you have SAC package installed)
unpack the SEED into SAC format. It is essential to keep data for a specific
station in a separate directory. Once data is unpacked, execute the code mk_macro. As presently compiled, it anticipates file naming convention
of the IRIS DMC (e.g., 1995.339.06.57.37.0000.CN.DRLN..BHE.SAC).
It will make two files: rewr_sac.m - a macro for SAC that will
rotate, demean, detrend and rename the data;
and in_recfuink input file for recfunk codes.
Once you have made rewr_sac.m - fire up SAC and execute
See example of in_recfunk.
- Step 3. Data quality control.
While it is not necessary given the procedure described above, looking at seismograms is advisable, to guard against unanticipated features.
In particular, things to check are:
presence of meaningfull signal on all 3 components;Remove file
names of suspect traces from in_recfunk.
how the time windows are cut. Start time and duration of the "signal" window may be adjusted within in_recfunk for each record;
NOTE No assumptions are made at this stage with respect to the origin of the
For the purposes of developing vertical 1D profiles of seismic properties via forward modeling
with synthetic seismograms (described in the following section)
a single-trace RFs will be required. These could be extracted from the tables generated by recfunk
with any text manipulation application (awk, perl, sed), by pulling out values for a chosen backazimuthal
or epicentral-distance bin.
- Data requirements.
SAC formatted records of 3 components of motion (Z,R,T) are needed for the computation of RFs.
Headers of SAC-formatted files
need to contain receiver and source positions, and fields "gcarc" and "baz" need
to be filled in.
The time window for analysis needs to have equal amounts of
pre-event noise and signal, e.g., 100 sec of each. Analysis windows
are defined explicitly for each triplet of records in the input file in_recfunk
NAMING CONVENTION for data files is as follows:
3 files are expected, with vertical, radial and transverse components of motion.
File names have to end with either r,t or z.
- Step 1.Set up a directory with data. Create a file called in_recfunk (name mandatory!)
with file names and time windows for analysis. It is practical to use full path names in in_recfunk
so that the code could be run in a separate directory.
See example of in_recfunk.
- Step 2. Run the code. recfunk is interactive, and promts the user for input.
The code writes out 4 ASCII files with backazimuthal and epicentral sweeps of radial and transverse RFs,
File names are generic: out[rt]_baz.grid out[rt]_epi.grid.
Each contains RFs for each bin, separated by > (for use with GMT plotting commands).
See example of the format. Columns in the table are time, baz, RF_value.
- Step 3. Visualize receiver functions. This could be done in many ways since the output of recfunk
is a set of ASCII tables. I like GMT, and use following scripts.
To plot backazimuthal sweeps, use RFwig-baz.plot shellscript.
A sample command looks like this:csh RFwig-baz.plot outr_baz.grid outt_baz.grid 0.3 (radial first, transverse second!), and will generate a file called RFwig-baz.ps
To plot backazimuthal sweeps, use RFwig-baz.plot shellscript.
Resulting plot will look something like this epicentral sweep of RFs band-limited at ~0.5Hz, computed for Lake Vostok in Antarctica.
- Step 4. Examine the plots and decide whether:
data is OK. Sometimes an odd record yields a wild oulier that swamps bins around it;
bin sizes in backazimuth and epicentral distance have been selected well;
time windows for plotting are selected right;
gain is approriate for the strength of the signal;
- Step 5. Re-run recfunk and/or plotting scripts untill happy.
Alternatively, an averaged receiver function may be created for the chosen subset of records
with the help of recfunk_average.ascii code. For this purpose, user is advised to
create a separate in_recfunk file that would have only records for a chosen direction and/or distance
ranges. This is especially important when T component RFs are to be averaged, as they
vary with direction in many cases. Files created by recfunk_average.ascii are called outr.trace,
outt.trace, and will be overwritten every time the code is invoked.
"Modeling" in this context means finding a vertical distribution
of seismic properties that would yield a synthetic seismogram that
looks close to the observed one. It is thus a classic case of
Step 1. Set up a model. Look here for model format description.
Step 2. Set up the data (radial and transverse RFs) in the form
of 2-column (time, value) ascii tables. These will be used for comparison plots only.
Step 3. Decide on modeling parameters:
Backazimuth (direction from the station to the source, clockwise from north) (e.g., 45 deg). This is important if models are anisotropic
phase velocity of the incident P wave
Typical phase velocities for a source at depth of 33 km:
30 deg - 12.6 km/s
60 deg - 16.3 km/s
90 deg - 23.88 km/s
120 deg - 58.56 km/s (PKP)
Step 4. Compute synthetics with datasynmod-plot script, by issuing a command like this:
csh datasynmod-plot 1 45 essmod 18 data_radial data_transverse
The output is written to a file called "datasynmod.ps"
NB: Modify plotting parameters if desired.
Step 5. Look at the plot like this one (in file "datasynmod.ps", generated by the above command from the model ) to
see if the waveforms match.
Red lines show synthetics, blue lines show data.
Model values plotted in the box may be changed by modifying the awk commands in the script .
Step 6. Modify the model and repeat 4, 5 until satisfied.