An example with the 2015 M7.8 Nepal earthquake
Here I will walk you step by step through running the Nepal earthquake. This is one of the pre-canned examples which you can find in the MUD/examples folder. Note this is using only the high-rate GPS stations. If you want to reproduce the full results that we calculated in the Galetzka et al., 2015 paper in Science, that's more complicated because there we used InSAR. That's for another time. The model you will obtain here is different (this is inversion after all!) from the one in the paper but this is a useful end-to-end exercise.
Parameter files used for inversion have a .inv extension. This is not required I just do it so I have an easy way of identifying what is what. I have a folder separate from everything else where I keep my parameter (.inv) files. So create a folder where you will save all your parameter files and then copy the file nepal_version_2.inv.py from $MUD/examples/nepal_inverse/ into that folder. In my case all my inversion files are in /Users/dmelgar/Slip_inv/run/ this is what that folder looks like for me after copying the relevant file:
Open the nepal_version_2.inv.py file and change the variables that define where the inversion results and auxiliary files will be. This is a different folder than where you save your parameter files, like so:
home='/Users/dmelgar/Slip_inv/'
project_name='Nepal_example'
run_name='example'
home is where all the different inversions will live. project_name will be the folder name for this inversion project, and run_name will be the name of the inversion output files. After this is done set the flags to only initalize the folder structure by doing:
init=1 #Initalize project
make_green=0 #Compute GFs
make_synthetics=0 #Compute synthetics for a given model at given stations
G_from_file=0# =0 read GFs and create a new G, =1 load G from file
invert=0 # =1 runs inversion, =0 does nothing
Nothing else matters for now. Next step is to run it. Go to the terminal and run the script by typing:
python /Users/dmelgar/code/MudPy/examples/nepal_inverse/nepal_version2.inv.py
You should now see the entire inversion project folder structure. On my system it looks like this:
You need 3 auxiliary files to run the inversion (i) a fault geometry file (ii) a layered Earth structure file and (iii) a Green's function list file (GFlist) that specifies the station locations and types of data. MudPy requires the auxiliary files to be in specific locations so from $MUD/examples/nepal_inverse/ copy the following files:
(i) nepal_10.fault is the fault geometry and should be placed in home/project_name/data/model_info/
(ii) avouac.mod is the Earth structure model and should be placed in home/project_name/structure/
(iii) nepal_simple_v2.gflist is the GFlist file and should go in home/project_name/data/station_info/
Now you need data to invert! I've provided 5Hz GPS waveforms in $MUD/examples/nepal_inverse/data/. There is data for 5 sites and three components of motion. These are provided as SAC files. Copy all the data into home/project_name/data/waveforms
You need to edit the file by hand. The GFlist file should reflect the type of data that each waveform corresponds to (displacement, velocity, tsunami, etc.) and where the data is. The GFlist file for my system looks like:
#station lon lat "static,disp,vel,tsun,insar" "Static file, displacement file, velocity file,tsunami file,insar file" "static sigmas(n,e,u),displacement sigmas(n,e,u),velocity sigmas(n,e,u),tsunami sigma,tsunami sigma"
KKN4 85.2788 27.8008 0 1 0 0 0 /foo/bar /Users/dmelgar/Slip_inv/Nepal_example/data/waveforms/KKN4.disp /foo/bar /foo/bar /foo/bar 1 1 1 1.0 1.0 1.0 1 1 1 0 0 0 0 1 0
NAST 85.3277 27.6567 0 1 0 0 0 /foo/bar /Users/dmelgar/Slip_inv/Nepal_example/data/waveforms/NAST.disp /foo/bar /foo/bar /foo/bar 1 1 1 1.0 1.0 1.0 1 1 1 0 0 0 0 1 0
RMTE 86.5971 26.991 0 1 0 0 0 /foo/bar /Users/dmelgar/Slip_inv/Nepal_example/data/waveforms/RMTE.disp /foo/bar /foo/bar /foo/bar 1 1 1 1.0 1.0 1.0 1 1 1 0 0 0 0 1 0
SNDL 85.7989 27.3849 0 1 0 0 0 /foo/bar /Users/dmelgar/Slip_inv/Nepal_example/data/waveforms/SNDL.disp /foo/bar /foo/bar /foo/bar 1 1 1 1.0 1.0 1.0 1 1 1 0 0 0 0 1 0
SYBC 86.7124 27.8143 0 1 0 0 0 /foo/bar /Users/dmelgar/Slip_inv/Nepal_example/data/waveforms/SYBC.disp /foo/bar /foo/bar /foo/bar 1 1 1 1.0 1.0 1.0 1 1 1 0 0 0 0 1 0
CHLM 85.3141 28.2072 0 1 0 0 0 /foo/bar /Users/dmelgar/Slip_inv/Nepal_example/data/waveforms/CHLM.disp /foo/bar /foo/bar /foo/bar 1 1 1 1.0 1.0 1.0 1 1 1 0 0 0 0 1 0
The five numbers 0 1 0 0 0 indicate the data type. These stations are all labeled as displacement waveforms. One station can have multiple data types. There are 5 columns that specify where the waveform (or static data) for that station lives. You'll note 4 of those columns are just empty place holders labeled /foo/bar/ this is just to indicate there's no data for that data type. Change the displacement column (as above) to point to where you are keeping the waveform data.
The parameter file should work as provided. It will calculate 5Hz Green functions and synthetics and lowpass filter them to 0.5Hz. It will use Laplacian spatial smoothing (regularization) and it will try 20 values of regularization strength (the regularization parameter) logarithmically spaced between 10^-2 and 10^1. The epicenter are already defined and it will assume a rupture velocity of 3.2km/s and use only a single time window. These are all parameters to be changed and explored to find better fits to the data but for now do not change them.
In preparation to run the inversion the path of least resistance is to set ncpus=1 this will use only 1 core on your machine and make the inversion a bit slow but should work fine. If you get tired of the inversion being slow take a look at this wiki entry on how to enable parallel processing. Most modern laptops, for example, have 4-12 cores that you can use all at once.
In the parameter file set:
init=0 #Initalize project
make_green=1 #Compute GFs
make_synthetics=1 #Compute synthetics for a given model at given stations
G_from_file=0# =0 read GFs and create a new G, =1 load G from file
invert=0 # =1 runs inversion, =0 does nothing
Then run from the terminal by typing:
python /Users/dmelgar/code/MudPy/examples/nepal_inverse/nepal_version2.inv.py
This will will calculate the Gfs and synthetics. You can se the output in home/project_name/GFs/dynamic/
In the parameter file set:
init=0 #Initalize project
make_green=0 #Compute GFs
make_synthetics=0 #Compute synthetics for a given model at given stations
G_from_file=0# =0 read GFs and create a new G, =1 load G from file
invert=1 # =1 runs inversion, =0 does nothing
And once again run from the terminal:
python /Users/dmelgar/code/MudPy/examples/nepal_inverse/nepal_version2.inv.py
With some luck you will get a lot of output like this:
Running inversion 16 of 20 at regularization levels: ls =2.3357214690901213 , lt = 0.0
... calcualting variance reduction...
... Static ABIC requested
... computing and saving synthetics...
... writing model results to file /Users/dmelgar/Slip_inv/Nepal_example/output/inverse_models/models/example.0015.inv
... inversion wall time was 0:00:00.699570, total wall time elapsed is 0:00:11.449828
The output files will be in home/project_name/output/inverse_models/models here's what it looks like for me
Because we asked for 20 different values of the regularization parameter there will be 20 sets of output files labeled 0000 through 0019. Look at the .log files, we know the Nepal earthquake had a mgnitude of M7.8 so we would like an inversion close tot hat value. The file example.0012.log should look like:
Project: Nepal_example
Run name: example
Run number: 0012
Velocity model: avouac.mod
Fault model: nepal_10.fault
G name: nepal
GF list: nepal_simple_v2.gflist
Solver: nnls
lambda_spatial = 0.78475997035146106
lambda_temporal = 0.0
Beta(degs) = 45
Mean rupture velocity (km/s) = 3.2
Number of rupture windows = 1
L2 = 1.0776498818937843
VR static(%) = nan
VR displacement(%) = 77.469276358268516
VR velocity(%) = nan
VR tsunami(%) = nan
VR InSAR LOS(%) = nan
Lm = 24.947717300397361
ABIC = 32900.955540993862
M0(N-m) = 1.0163035547271009e+21
Mw = 7.9380156297427131
This is close to the target magnitude and has a variance reduction of 77.4% for the displacement waveforms. 100% is a perfect fit so that value is quite high and a good sign. Now let's plot the waveform fits and the slip model.
At your terminal start an ipython shell by typing
ipython
Then within this ipython shell load the plotting functions by doing
from mudpy import view
Plot the waveform fits by doing
view.synthetics(home,project_name,'example','0012','nepal_simple_v2.gflist','d',None,[0.5],[0,60],'lon',1,'gps')
This should make the following plot:
The black lines are the data used in the inversion and the red lines are what is predicted by the model. The peak values are indicated by the numbers next to each waveform.
You can also make a quick and dirty plot of the slip distribution by doing:
view.slip3D('/Users/dmelgar/Slip_inv/Nepal_example/output/inverse_models/models/example.0012.inv')
which should make this plot:
If you want to make nicer plots (with GMT for example) the slip data is in the file example.0012.inv.
Now Is this model good enough? That's up to each researcher. This is a good start, the static offsets are well modeled and some of the dynamic features too. However you can see that the synthetics (red waveforms) arrive "before" the data meaning the rupture speed is probably too fast. You should now try running the model with 4 time windows. That will probably improve the fits by allowing slip at later "slower" times.