This repository hosts MCQsim, a multi-cycle earthquake rupture simulator for generation the of earthquake sequences along predefined faults or fault systems of arbitrary geometric complexity. MCQsim uses graphical user interfaces (GUIs) for most of data input and output. These GUIs were developed in MATLAB. Stand-alone versions of these GUIs, along with the MATLAB GUIs can be found here, e.g., under "01_Mstuff". MCQsim itself is written in C and parallelized via MPI. It works on large-scale computational infrastructure as well as personal computers. The code can be found under "00_Cstuff".
The most recent version of MCQsim, provided here, employs QuadTree and H-matrix to a) decrease the size of the stiffness kernel, and b) decrease computation time. This update dramatically increased code efficiency and scalability, especially for larger-scale simulations with 100k+ fault elements. Future releases are planned to include pore-pressure changes as well as surface topography, further guided by the specific needs by our fellow users (i.e., you).
Please check the following publication for details on MCQsim model formulation and cite when using the code
Zielke, O., and P.M. Mai (2023). MCQsim: A Multicycle Earthquake Simulator. Bull. Seismol. Soc. Am. 113(3), 889-908, doi:10.1785/0120220248. *.pdf (open access)
Zielke, O., and P.M. Mai (2025). Does Subsurface Fault Geometry Affect Aleatory Variability in Modeled Strike-Slip Fault Behavior?. Bull. Seismol. Soc. Am. XX(Y), 1-17, doi:10.1785/0120240152. *.pdf
Contact: olaf.zielke@kaust.edu.sa
Following are a few example outputs, generated with MCQsim for a "toy example" listric normal fault that features geometric roughness and heterogeneous strength distribution. We present animated GIFs of the single-event quasi-dynamic evolution of a M7.54 event (slip and slip-rate). We further present the earthquake catalog and catalog derivates for the same fault.
| M7.54 normal faulting EQ (slip) | M7.54 normal faulting EQ (rate) |
|---|---|
 |
 |
Earthquake catalogs, created with MCQsim, can serve as an input for PSHA, providing an earthquake rupture forecast for a modeled fault (system). As such, they provide internally consistent, physics-based branches of the PSHA logic tree. We recently streamlined the post-processing workflow to feed MCQsim earthquake catalogs into the OpenQuake PSHA engine. Below is an example seismic hazard map for toy model listric fault, using QGIS and OpenQuake's IRM toolkit.
Following below are the results of strong scaling tests that we performed on Shaheen3@KAUST. For these tests, we created 5-kyr-long catalogs along two intersecting faults with average slip rates of ~3mm/yr and parameterized by 10k, 50k, 100k, and 200k elements. We show these scaling tests here a) to showcase MCQsim performance on our HPC facilities, and b) to provide users with first-order guidence on simulation preparation (e.g., CPU number as a function of fault element number i.e., fault system size). Below, we show scaling in terms of speed up factor and efficiency (relative to the respective simulation with the least amount of CPUs). Speed up closely follows the slope of ideal scaling, before beginning to deviate from it as CPU number exceeds some threshold. The inset table indicates the "serial" part of code that is currently not parallelized (e.g., catalog output to screen and file). The "odd" shapes in effciency curves for the 10k and 50k element simulations reflect the internal CPU organization within a single node of Shaheen3 (each node consisting of 192 CPUs, single-node indicated by gray area).
IMPORTANT: Keep in mind that the results of this scaling test are specific to the model setup and Shaheen3. It is best-practice to perform these tests yourself on the HPC facilities you are going to use.
| Element number | Resolved Elements | Stiffness matrix (H-mat) | Stiffness matrix (classical, N^2) |
|---|---|---|---|
| 10,000 | 2,770 | 1.1 Gb | 3.6 Gb |
| 50,000 | 3,656 | 14.1 Gb | 90 Gb |
| 100,000 | 3,248 | 45.5 Gb | 360 Gb |
| 200,000 | 3,775 | 182.1 Gb | 1,440 Gb |
Our H-matrix implementation drastically reduces the number of "resolved elements" and therefore the number of individual fault interactions that need to be computed.
Here is a simple way to download specific folders from GitHub, as opposed to downloading the entire project:
open a browser tab and go to
https://download-directory.github.io
then copy and paste the GitHub folder you want to download, e.g.,
https://github.com/OlafZielke-EQ/MCQsim/tree/main/TREAD/
into the corresponding field and the download will start automatically.
The MCQsim worklflow can be broken down into three components, namely pre-simulation, simulation, and post-simulation. The first and third component rely heavily on MATLAB-based GUIs that we developed for this purpose. The actual simulation is "terminal" based.
MCQsim simulations require a number of files to be located in the same folder as the executable. They are created with the MATLAB input GUIs.
| File name (input) | File content |
|---|---|
| FaultModel.txt | Parameter file. Contains details about how to run the simulation e.g., catalog length. ASCII format. |
| FaultModel_FLT.txt | Summary file. Contains information about fault geometry. ASCII format. |
| FaultModel_FLTtrig.dat | Data file. Contains the 3-D fault geometry, including the list of all triangular fault elements. Binary format. |
| FaultModel_Summary_RoughStrength.txt | Summary file. Contains information about used fault roughness, strength, and friction configuration. ASCII format. |
| FaultModel_1_Rough.dat | Data file. Contains the 3-D fault geometry. This file is equivalent to _FLTtrig.dat except that fault roughness might be added. Binary format. |
| FaultModel_1_Strgth.dat | Data file. Contains friction and strength parameter for each fault element. Binary format. |
| FaultModel_BND.txt (optional) | Summary file. Contains information about of boundary surface geometry. ASCII format. |
| FaultModel_BNDtrig.dat (optional) | Data file. Contains the 3-D geometry of boundary surface triangular elements. Binary format. |
The parameter file has a specific structure that must not be changed. Do not include or remove lines and do not enter spaces in the descriptions. Following is an example parameter file "FaultModel.txt".
Run_Parameter_Information
--------------------------
InputName: FaultModel
Realization_Number: 1
CatalogTypeNr(1/2/3): 3
PlotCatalog2Screen(0/1): 1
RandSeedValue: 1001
MinElementNum4Catalog: 1
ContinueCatalog(0/1): 0
LoadPrevious_Kh_mat(0/1/2): 1
Kh_mat_file_2_load: FaultModel_Khmat.dat
--------------------------
StoreSTF4LargeEQs(0/1): 1
MinMagnitude4STF: 7.0f
UseRuptPropag(0/1): 0
MinMagnitude4RupProp: 6.5f
--------------------------
IntSeisLoadStep(days): 1.0f
pow2forInterSeisSteps: 8
--------------------------
Visco_AfterSlip(years): 3.0f
Visco_DeepRelax(years): 20.0f
--------------------------
HealingDuringCoseisFraction: 0.05f
OvershootFraction: 1.2f
PreStressFraction: 0.90f
MinCoSeisSlipRate(m/s): 5.0E-3
--------------------------
EQrecordLength(years): 5000.0f
--------------------------
A brief explanation of the parameter file entries
| Parameter information | Description |
|---|---|
| InputName: | This is the "root" part of file name e.g., "FaultModel" of file "FaultModel_FLT.txt". All input files of the simulation need to share this input (see input files example above). |
| Realization_Number: | MCQsim was built to facilitate parameter space investigations, where multiple realizations can be generated with input GUI "B_MCQsim_RoughStrgth" (while only changing selected parameter). The realization number refers to an index in the file name of the "*Roughn.dat" and "*Strgth.dat" file. For example. For example, for realization 1 would, these files would be named "FaultModel_1_Roughn.dat" and "FaultModel_1_Strgth.dat"; for realization 42, this would be "FaultModel_45_Roughn.dat" and "FaultModel_45_Strgth.dat" |
| CatalogTypeNr(1/2/3): | Catalog type refers to the amount of information that is provided in the catalogs. Type 1 contains only general event characteristics such as time, magnitude, hypo-center, moment-centroid, rupture area, mean slip, and mean stress drop. Type 2 further contains the event's moment rate function. Type 3 further contains information for each fault element that was active (slipping) during the event. This includes the element's rupture start time, slip in strike- and dip-direction, coseismic (static) stress change, and stability type. MCQsim first generates a "RAW" catalog that corrresponds to Type 3. This "RAW"catalog can be processed as part of the post-simulation workflow to generate other catalog types that are smaller in size (not every investigation may have the need to consider element-based information). |
| PlotCatalog2Screen(0/1): | This switch (0 = no, 1 = yes) indicates if catalog information should be printed to screen during the simulation. It is typically left at 1. |
| RandSeedValue: | If desired, some parameters can be randomized during simulation. For example, friction parameters of activated fault elements may be permitted to vary (within specified range) at the end of the coseismic phase. Reusing the same seed value allows to create the exact earthquake sequence. IMPORTANT each process/CPU gets its own seed value (the process' MPIrank is added to "RandSeedValue"). Hence, in order to create the exact same catalog, it is important to not only use the same seed value, but also the same number of processes. |
| MinElementNum4Catalog: | Minimum number of fault elements that need to be activated in an event to add the event to the catalog. If set to 1, all events are included in the catalog; if set to 10, only events with 10+ active elements are included in the catalog. IMPORTANT The element number refers only to the writing part (the smaller events still exist in the simulation, they are just not written to file). |
| ContinueCatalog(0/1): | This switch (0 = no, 1 = yes) indicates if an existing catalog should be extended. |
| LoadPrevious_Kh_mat(0/1/2): | This switch indicates if a new stiffness matrix "*Khmat.dat" should be created or not, and whether it should be saved to file (e.g., to be used when continuing an existing catalog). switch = 0 means to create a new stiffness matrix but NOT storing it to file. switch = 2 means to create a new stiffness matrix and storing it to file. switch = 2 means to load an existing stiffness matrix. A typical setup for this and the previous entry would be ContinueCatalog 0 and LoadPrevious_Khmat 1_ or ContinueCatalog 1 and LoadPrevious_Khmat 2_ |
| Kh_mat_file_2_load: | The file name used to either read (if LoadPrevious_Kh_mat 2) or write (if LoadPrevious_Kh_mat 1) the stiffness matrix "*Khmat.dat" |
| StoreSTF4LargeEQs(0/1): | This switch (0 = no, 1 = yes) indicates whether the quasi-dynamic rupture evolution of individual events should be saved to file. |
| MinMagnitude4STF: | The minimum magnitude an event must have in order to write its quasi-dynamic rupture evolution to file. |
| UseRuptPropag(0/1): | Not used in current code version. |
| MinMagnitude4RupProp: | Not used in current code version. |
| IntSeisLoadStep(days): | The minimum time step (in days) used for post-seismic i.e., interseismic loading and stress redistribution. |
| pow2forInterSeisSteps: | MCQsim uses the minimum time step (previous row) and a binary search to determine the loading step to the next event. The value of pow2forInterSeisSteps refers to just that. For example with a value of 8, MCQsim first jumps 2^8 days ahead in time to identify if the combination of inter-seismic loading and post-seismic stress redistribuiton caused at least one element to fail. If not, loading proceeds. If yes, then the loading step is halfed, etc. |
| Visco_AfterSlip(years): | Postseismic response ALONG THE FAULT, employing a simple Maxwell spring-dashpot model. This is the time it takes to decrease a postseismic signal to 1/e (approx. 35%) of its original value. Providing this measure in years (rather than using viscosity) seems more intuitive. For example the postseismic signal/effect (here after slip) decreased by ~65% within the first 3 years. Higher value corresponds to higher viscosity. |
| Visco_DeepRelax(years): | Postseismic response ALONG THE BOUNDARY (if used), employing a simple Maxwell spring-dashpot model. This is the time it takes to decrease a postseismic signal to 1/e (approx. 35%) of its original value. Providing this measure in years (rather than using viscosity) seems more intuitive. For example the postseismic signal/effect (here visco-elastic relaxation) decreased by ~65% within the first 20 years. Higher value corresponds to higher viscosity. |
| HealingDuringCoseisFraction: | Fraction by which the friction coefficient of activated (unstable or conditionally stable) fault elements increases coseismically (during each coseismic iteration step) once these elements stopped to slip in the event (while the rupture itself is still ongoing). |
| OvershootFraction: | Describes the ratio of dyanamic and arrest friction coefficient. |
| PreStressFraction: | Faults are pre-stressed at the beginning of a simulation to this fraction of their frictional strength (unless an existing catalog is continued). This is done to speed up the intial "burn in" phase i.e., decreasing the computational time for the modeled fault system to reach its dynamic equilibrium. |
| MinCoSeisSlipRate(m/s): | Minimum slip rate to be considered coseismic slip. |
| EQrecordLength(years): | Record length in years. If a catalog is continued, this value refers to the added length and not to the total length. |
Assuming all files are present, a simulation can be started via:
mpirun -n 6 ./MCQsim FaultModel.txt
This command will start MCQsim as an MPI run, using n = 6 CPUs. The parameter file "FaultModel.txt" contains relevant information to start/control the simulation. MCQsim then runs and produces a number of output files. Most importantly the " *Catalog.dat" file and a number of " *.srfb" files.
MCQsim creates a number of output files, as listed below. They can be visualized using our MATLAB-based GUIs. Catalog files can be opened with GUIs starting with a "Z" such as "Z_Catalog_MCQsim_TREAD" or "Z_MCQsim_LoadCatalog". They also allow you to export the catalog to either *.mat or *.txt file. The PreRun file as well as the quasi-dynamic rupture models of individual events can be opened with GUIs starting with "Y" such as _"Y"PreRunAndSRF_MCQsim_Tread.
| File name (output) | File content |
|---|---|
| FaultModel_BrchInfo.dat | Data file. Contains information about the QuadTree structure. Useful for OpenQuake pipeline. Binary format. |
| FaultModel_Khmat.dat (optional) | Data file. Contains stiffness matrix for elastic interaction. Allows continuing an existing simulation/catalog without re-calculating Kh. Binary format. !large file size! |
| FaultModel_PreRunData.dat | Data file. Friction properties, fault strength, long-term slip-rate of fault elements at start of simulation. Binary format. |
| FaultModel_PostRunState.dat | Data file. Stress state after simulation is over. Allows continuing/extending an existing catalog. Binary format. |
| FaultModel_RAWCatalog.dat | Data file. Earthquake catalog. Binary format. |
| FaultModel_M7.38472_t3845.8374.srfb | Data file. Quasi-dynamic earthquake rupture for selected large event. Name provides event magnitude and time. Tthis is a binary version of a *.srf file (standard ruputre format; Graves, 2002). Binary format. |