-
Notifications
You must be signed in to change notification settings - Fork 63
Stochastic slip (fakequakes)
STOP. Before you do anything please follow these instructions to install MudPy.
Ok, moving on. In this entry I will explain how to generate stochastic slip distributions using the fakequakes module. All of the files for this example can be found in examples/fakequakes/planar
. I will discuss only a planar geometry. Fakequakes can also generate slip distributions on a 3D geometry, for example from Slab 2.0, I will cover that in another entry. Similarly you can synthesize HR-GPS and strong motion waveforms from the stochastic slip models. This will also be covered elsewhere. For now, the goal is simply to teach you to generate the slip model itself.
The details behind the fakequakes implementation can be found in these two papers:
-
LeVeque, R. J., Waagan, K., González, F. I., Rim, D., & Lin, G. (2016). Generating random earthquake events for probabilistic tsunami hazard assessment. In Global Tsunami Science: Past and Future, Volume I (pp. 3671-3692). Birkhäuser, Cham.
-
Melgar, D., LeVeque, R. J., Dreger, D. S., & Allen, R. M. (2016). Kinematic rupture scenarios and synthetic displacement data: An example application to the Cascadia subduction zone. Journal of Geophysical Research: Solid Earth, 121(9), 6658-6674.
Begin by copying the file subduction.fq.py
into a folder of your choice where you will keep your fakequakes runs. For example on my machine I have a folder /home/dmelgar/fakequakes/run
where I keep all my parameter files.
This is the parameter or control file that defines all the details of your models. Open it in a text editor or IDE (I use Spyder). The first thing we will do is initalize the file structure. The first lines of code look like this:
Set the home
variable to a folder of your choice where you will keep your fakequakes projects. Also set the project_name
variable to some suitable name that describes the project you will be working on.
Next notice the set of action flags in lines 19-26. These are used to define which actions you want the code to take. Set init = 1
and run the parameter file by typing at your terminal $> python /path/to/file/subduction.fq.py
This will simply create the folder structure needed by fakequakes
NOTE MudPy and fakequakes are still predominantly a python 2.7 code. Make sure you use a python 2.x interpreter when running these files. Porting to Python 3.x is underway.
After running the parameter file you can open your file browser and verify that indeed the relevant folders have been created, you should see something like this:
You have now created the folder structure. In order to avoid mistakes and confusion fakequakes and MudPy are very rigid in the sense that certain files have to be in certain folders. We will look at those auxiliary files next.
You have to decide a priori what fault geometry you want your stochastic slip models to be on. An example of such a file is, again, under examples/fakequakes/planar
. Copy this file subduction.fault
to the folder /data/model_info
within your fakequakes project folder.
Open the file in a text editor to see its contents. It should look something like this:
The columns are self explanatory. The fault is broken up into many subfaults. The file contains the geographic coordinates of the centers of each subfault as well as the depth in kilometers. The strike and dip of each subfault are defined as well as the length and width (in meters). Later on if you want to generate slip models on another geometry you need to create a .fault
file with this exact structure. Note also that this file has ~7000 subfaults. This is because we are using relatively small subfaults, 4000 x 4000m over a relatively large area. The number of subfaults and size of the fault model will be entirely dependent on the application.
The second file you will need is a .mod
file that defines the Earth structure. Fakequakes can handle homogenous half-space models and layered Earth models. There are examples of both in examples/fakequakes/planar
. Lets use a layered Earth model, copy the file bbp_norcal.mod
to the folder /structure/
within your fakequakes project folder.
Open the file in a text editor to examine it. This is a very detailed layered Earth model for northern California used by the SCEC BBP project. The file should look like this:
Due to somewhat silly reasons, which I will not explain here, the columns in that file are ordered as follows: layer thickness in kilometers, S-wave speed in km/s, P-wave speed in km/s, density in kg/cm^3, S-wave attenuation, and P-wave attenuation. Do not add headers to this file or unexpected results might occur. Also note that the final layer has thickness of 0. this indicates an infinite halfspace below the previous layers. If you wanted to have just a homogenous halfspace you could have a .mod
file with only one 0 thickness layer.
That's it! It really is that simple you just need those two files, the .fault
file and the .mod
file together with the parameter file to get some nice stochastic slip models.
Before moving on I will just point out that if you are interested in stochastic models with a diversity of magnitudes you don't need many different .fault
files. You can simply have one fault geometry that is big enough to fit the range of magnitudes you are interested in and the code will generate ruptures with random parts of the bigger geometry as the center. This will become clear with the examples below and is also discussed in the Melgar et al. (2016) paper at great length.
Ok, now back to the subduction.fq.py
parameter file. If all you want to do is generate slip models not all of the parameters will be relevant. Many of them are related to waveform synthesis. For simply generating slip models these are the relevant ones:
I will not discuss all of these in detail. For simplicity's sake the relevant ones are:
-
fault_name
is the name of your.fault
file -
model_name
is the name of your.mod
velocity structure file -
distances_name
is the name of a matrix variable which will be created at run-time. This matrix will hold the along-strike and along-dip distances between all subfaults. Simply give it a name that relates to the geometry you are using. -
UTM_zone
, fakequakes projects from geographic to cartesian coordinates using the UTM projection. This requires you define the UTM zone (roughly) of the center of your fault model. -
Nrealizations
defines how many slip models you want per magnitude -
target_Mw
is a vector of desired magnitudes. In the erxample above it will be M7.5 through M8.5 in 0.1 magnitude unit increments. Thus there will be a total of 50 slip models generated (becauseNrealizations=5
) -
max_slip
is some reasonable bound (in meters) of what the maximum allowed slip will be -
time_epi
is the UTC date and time at which your synthetic events will originate -
force_hypocenter
ifFalse
then for every slip model a random subfault will be chosen as the hypocenter. IfTrue
then the hypocenter will be forced to whatever value is inhypocenter
.
If you are interested in what the other parameters they are all explained in Melgar et al. (2016).
Ok, you're ready to go. Now simply change the action flags at the top. Set init = 0
(otherwise you will initialize the folder structure and delete all that you have done). Now set make_ruptures = 1
. Run the file once again by typing at your terminal $> python /path/to/file/subduction.fq.py
you should see output in your terminal:
This first step is a little slow if the fault model has many subfaults. Thankfully you only have to do it once. After the inter-fault distances have been calculated once you will see two .npy
files (binary numpy arrays) appear in the data/distances/
folder. Once these are there simply set load_distances = 1
for future runs and these arrays will be loaded from disk rather than computed again. This speeds up the process greatly if you want to change other parameters in the parameter file and re-run the whole thing.
Next, after the inter-subfault distances are calculated the code will generate the stochastic slip distributions. You should see something like this.
And that's it, when that finishes you have your slip models.
For every rupture two files will be generated. A .log
file with summary information about the rupture and a .rupt
file with the actual slip model. All of these files will be found under output/ruptures
. Note that the value of the variable run_name
from the parameter file is used to label the files. If you change some parameter and want to save the new runs with another name (so as not to overwrite some previous run) simply change run_name
to another value.
The .log
file with the summary information will look something like this, note that the actual values in your run will be different from mine since this is a stochastic method:
The .rupt
file with the slip information looks like this (again values will change from run to run):
The columns are labeled. Some important ones are:
-
dura
, this poorly named column, has the rise time (in seconds) of a particular subfault. Refer to Melgar et al. (2016) to see how the rise times are defined. -
ss-slip
, this column has slip in the "strike-slip" orrake = 0
direction. -
ds-slip
, this column has slip in the "dip-slip" orrake = 90
direction. -
rupt_time
, this column has the time (in seconds) at which a particular subfault starts slipping.
Some observations:
- The total slip at a subfault will be given by (ss-slip^2 + ds-slip^2)^0.5
- Note that even though we set
rake = 90
there are non-zero values in thess-slip
column. Once again this is due to the stochastic nature of the code. On average the rake across the fault model will be 90, however small random deviations from this mean behavior are artificially introduced. - Note some subfaults have zero slip. That's ok. Because this is a stochastic method not all subfaults will always participate in rupture.
- If you are interested only in static or instantaneous slip models you can safely disregard the
rupt_time
column in the.rupt
output file.
I have written a visualization module with some simple functions to view the slip models. If you want to make pretty plots I suggest using GMT. In any case to view one of the slip models in python you can do:
>>> from mudpy import view
>>> rupture='/Users/dmelgarm/FakeQuakes/subduction_test/output/ruptures/subduction.000047.rupt'
>>> view.slip3D(rupture)
This will produce the following 3D plot, again, note that yours will look different:
And that's that! Read the papers at the top of this wiki entry, play around with the parameters in the parameter file and send me an email if you have any questions.