We need to know the spatial and temporal behaviors of the earthquake’s rupture in order to study the physics of earthquake sources.
This is one of the fundamental methods proposed by Professor M. Kikuchi and Professor Hiroo Kanamori [Link] to calculate teleseismic body-wave inversion for modeling the earthquake source.
Nowadays, there are many methods and approaches to computing more information about the detailed rupture process. An excellent paper review has been done by Professor Satoshi Ide [Link 1, Link 2] entitled “Slip Inversion” as a book chapter in Treatise on Geophysics [Link]. This paper provides us detail explanation and development history of the methods/theories about slip inversion of the earthquake.
There is also a helpful book written by Prof. Agustin Udias, Prof. Raul Madariaga, and Prof. Elisa Buforn, entitled “Source Mechanisms of Earthquakes, Theory, and Practices” [Link]. This book also provides the detail theories, methodology, and some technical issues to provide the source mechanism of the earthquake.
You can also see some finite fault model database that I have compiled before. [Click Here].
Please prepare these programs first:
- Install Rdseed [Link]
- Install the Seismic Analysis Code (SAC) [Link]
- Install Numerical Recipes [Link]. Numerical Recipes is available by open access to users in developing countries, including Indonesia. In this case, we need Numerical Recipes in Fortran.
- Kikuchi’s and Kanamori program from the website maintained by ERI, Univ. of Tokyo [Link].
Step by step:
- Download package of Numerical Recipes
- Extract (easily using tar -xzvf)
- Go to the directory of Numerical Recipes
- Enter directory “f”
- Delete a file named libnrf.a
- Open file named jacobi.f and ludcmp.f
- In those files, change value in line 4 (NMAX = 500) to be NMAX = 10000 (NMAX must be set at a larger value than 1000)
- Close those files
- Open file Makefile
- In line 86 change $(FLIB) to “flib”
- Close those files
- Open terminal, type: make flib
- A new libnrf.a is available now
- Copy this file (libnrf.a)
Next, download the inversion program from the website. You can also download the sample files. Extract both of “Program” and “Sample” to directories.
- Go to directory “Program”
- Make a new directory named Recipes (easily using the command: mkdir Recipes)
- Paste the libnrf.a but rename it to be recipes.a
- Next, go to directory “Subroutine”
- Open a file named “Readme”
- This is a file described how to update our library file
- Copy all of the commands in the “Readme” file
- Open the terminal, paste, and it will be run
Next, go to directory “Sample”, cut the folder Input and Output and then paste in directory “Program”. Now, the directory “Program” has 5 sub-directories, i.e. Main, Recipes, Subroutine, Input, and Output.
Now, we will compile all of our Fortran (main programs). Easy, there is a Makefile in the directory Program. But, we need to change the Makefile to be executable. Open terminal (Ctrl+Alt+T), type sudo su, enter your Linux password, chmod 777 Makefile), the close terminal.
Then you can easily compile the Fortran program that you need. For example type: make green to get an executable program to calculate Green’s function using a given velocity model.
In this stage, you can also modify the .bashrc file to write a path to this directory, thus you can use the program anywhere in your Linux environment.
Preparing Our Seismic Data
- Select our case study (an earthquake).
- We can download seismic data from some institutions, e.g. IRIS, GFZ, Japan’s network, BMKG network, etc. I have compiled some of the data resources on this website [Link].
- In this tutorial, we use IRIS’s Wilber 3 facility.
- Easy, select your earthquake.
- Select stations with epicentral distance 30-90 degree.
- Why? Because for data from stations with a distance less than 30 degrees, the waves suffer strong perturbations by the structure of the upper mantle. For data from stations larger than 90 degrees, the waves are diffracted because they reach the nucleus of the Earth.
- Download the data with some length for example 5 minute before P arrival and 60 minutes after P arrival or 3 minutes after S arrival. Then, select waveform with sampling every 0.05 seconds (20 Hz) or usually named as BH* channel. Seismic data usually configured to 5 levels of samples: HH a 100Hz; SH a 50Hz; BH a 20Hz; LH a 1Hz & VH a 0.1Hz [more information, click here: Link].
- Note: waveforms from IRIS wilber3 have event information inside its header (evla, evlo and evdp). If you use waveforms from other agencies (e.g. GFZ or BMKG), please check whether the event information has been put in the data header or not. If not, easy, you can use the command: ch (change header) to fill the event information. IRIS’s waveforms also may have arrival time data on its header.
- Download the seismic data in a format of seed (because we also need the RESP file) [Link]
- Make a new directory named Data, put the seed file inside it.
- Extract this rdseed data using the command: rdseed -f filename.seed -R -d -o 1 -p [see this tutorial: Link]
- Then the folder contains a lot of files including SAC files, RESP files, and pole-zeros files.
- Select waveforms with same (uniform) data, for example only select 00.BH* data. Remove all other data except data with 00.BH* (00.BHE, 00.BHN and 00.BHZ)
- Make a script named “mk_conv.farm.pl” according to the manual from the code’s website.
- Change the value of high_pass_co, low_pass_co, Duration, pre_ev, and sample (depend on your purposes).
- Don’t forget to provide files named “hypo” and “sacmacro“
- In line 61, change the coda of the SAC file, according to your data, for example, if we download the data from IRIS, the coda will be “M.SAC” instead of “Q.SAC“.
- In the end of “sacmacro” file, add command “quit”
- Or you can also add some additional commands such as rmean, rtr and taper in sacmacro file (read the function of this command in manual of SAC).
- Run the Perl program (make sure “Perl” have been installed in your Linux)
- Use command: perl mk_conv.farm.pl
- Now two new files are available i.e. response and i_conv.farm
- Make sure there is no “error” in every step.
- Go to the main directory
- Copy some input data e.g. jb.ptime, jb.stime, jb.table to directory “Data”.
- In directory Data, run ./conv.sac.farm < i_conv.farm
- This will result in a new file named fort.22
- Then we will convert fort.22 to fort.1 (main input of inversion program)
- Use command: ./rotSH to produce fort.1
- Again, make sure there is no error in every step.
This is the main purpose of this tutorial, but if you have prepared your input data very well, you can continue easily to this step by reading the manual. There are there different sets of inversion program i.e.:
Set 1: Inversion allowing mechanism changes.
Set 2: Fine tuning of the source time function and slip distribution.
Set 3: Inversion with a fixed fault mechanism.
For each step, you can see the input, stdin, and output in the manual which is available on the website. Use command: ./program name < stdin [general format]. [more information click here: Link]
To run the program:
- ./green < i_green
- ./inversion < i_inversion
- ./graphics < i_graphics
The parameters that you need to set in i_green, i_inversion and i_graphics depend on your cases and situations. You can read more descriptions in the manual.
For more description of theory and methodology, read these papers:
Kikuchi, M., & Kanamori, H., Bull. Seism. Soc. Am., 72, 491-506, 1982. [Link]
Kikuchi, M., & Kanamori, H., Phys. of the Earth and Plan. Int., 43, 205-222, 1986. [Link]
Kikuchi, M., & Kanamori, H., Bull. Seism. Soc. Am., 81, 2335-2350, 1991. [Link]
Kikuchi, M., Kanamori, H. & Satake, K., J. Geophys. Res., 98, 15797-15808, 1993. [Link]
If you have further questions, don’t hesitate to send me a message (email@example.com).