Fault Scarp Dating Tool - a MATLAB code for fault scarp dating using in-situ chlorine-36 supplemented with datasets of Yavansu and Kalafat faults

We publish a MATLAB code used to analyze concentration profile of cosmogenic 36Cl accumulated in-situ through a rupture history of the fault scarps in western Turkey (Mozafari et al., 2019). The code is a version of the forward modeling Matlab code -Fault Scarp Dating Tool- (Tikhomirov, 2014). The code models a 36Cl profile accumulated in the fault scarp surface through a guessed rupture history, and compares the modeled and measured 36Cl profiles with statistical tests. Rupture histories are randomly generated in bounded solution space using Monte-Carlo method or optimized using Random Walk algorithm to achieve the best fit of the modeled and measured 36Cl profiles. The code has a user-friendly interface, a build-in help and an example of input data.


Data
The following data are presented:

Experimental design, materials, and methods
The full description of the methods used to acquire the datasets of the Kalafat and Yavansu fault scarps is published in the main publication [1]. The description includes the sampling in the field, mechanical and chemical preparation and details on sample measurement with accelerator mass-spectrometry. In this publication we would like to present the computation code used to analyze the datasets.
Fault Scarp Dating Tool (FSDT) is the MATLAB® code for fault scarp dating and earthquake history reconstruction. The code applies forward modeling and Monte-Carlo algorithms to determine the most probable displacement history: the age of displacement events, the displacement values (slips) and the uncertainty of ages and slips. Displacement histories are randomly generated, transformed to 36 Cl concentration profile and compared to measured concentrations with goodness-of-fit tests. The current code is an improved and updated version of the previous FSDT code published in [2].
The program has build-in help files and a user-friendly graphic interface. Output of the program is presented graphically and numerically and can be saved in a variety of file formats. The program was tested on MATLAB® versions 8.5 and higher. To run the code, a MATLAB® assembly should include Parallel Computing Toolbox™. The source code can be distributed and modified under Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License.

Production of cosmogenic 36 Cl in the fault scarp
A mathematical model of cosmogenic isotope production in a fault scarp is essential to transform exposure history into 36 Cl distribution. Classic surface exposure dating approach for calculation of accumulated isotope concentration cannot be applied to a fault scarp, due to complex and changing geometry of the scarps. In our previous publication, we formulated a mathematical model for shielding factor calculation in the fault scarp geometry and gave example with high-energy neutrons [3]. In our new code, this model is broadened to all particles which share 36 Cl production: high-energy neutrons, thermal and epithermal neutrons, fast and negative muons. In this new model, high-energy neutrons flux is calculated in the form of the shielding factors over the sample stripe volume. Relying on highenergy neutron attenuation length perpendicular to the scarp surface, thermal and epithermal neutron fluxes are calculated over the samples above colluvium using two-group diffusion approach [4]. To derive these fluxes for samples under colluvium, high-energy neutron flux over the colluvium surface and attenuation length of high-energy neutrons perpendicular to colluvium surface are additionally precalculated. Fast muon flux is calculated over the sample stripe with a code modified after [5] to fit fault scarp geometry. This modified code follows the muon production model by [6,7] and the muon scaling by [5]. The rate of stopped negative muons over the sample stripe is calculated as a derivative of one exponential approximation of fast muon flux with respect to scarp geometry. Using fast and negative muon fluxes, the production of thermal neutrons by muons is calculated. Production of thermal neutrons from decay of 238 U and (a,n) reactions by U and Th series is calculated following [8]. In total, 36 Cl production rate on Ca, K, Ti, Fe and 35 Cl by all significant reaction channels are derived over the sample stripe volume. These calculations of the production rate are done for the different relative positions of footwall and colluvial wedge, which correspond to periods of quiescence between ruptures.

Two-step computation
Direct modeling of one exposure history of a fault scarp takes about a half of an hour with our code on a modern office PC. The computation is done at a huge amount of points of the scarp: several points are required to average concentration over the sample volume, the sample stripe usually contains dozens of samples, and few fault scarp configurations interchange each other through scarp history. To reduce computation time and to apply the Monte-Carlo method, we divided the computation process into two stages.
At the first stage, the particle fluxes and the shielding factors are calculated at nodes of 3-dimentional grid, which covers all possible positions of the sample stripe and all possible fault scarp configurations in the past. The dimensions of the grid are the depth perpendicular to the scarp surface, the position along the fault scarp height and the fault scarp height along its surface. Additionally, sample coefficients for all production channels and scaling factors for high-energy neutrons [5] are calculated at this stage. The result of the first stage is saved in a separate file and serves as a database for the second calculation stage.
At the second stage, the code randomly generates an exposure history with a predefined number of displacement events. The particle fluxes and the shielding factors are interpolated to each particular scarp configuration using the database. 36 Cl distributions accumulated during each scarp configuration are calculated and merged to represent the exposure history. In comparison to direct modeling, the modeling of an exposure history with the database support takes about a tenth of a second. Creation of one database requires several hours, but in the light of numerous simulations with the Monte-Carlo method, two-stage computation gives tremendous time profit.

Monte-Carlo and Random Walk algorithms
FSDT code applies simple Monte-Carlo or Random Walk algorithms to search the most probable exposure history. The ages of displacement events, the colluvium positions, the beginning of exposure and the erosion rate are the desired parameters which uniquely describe the exposure history in our model. The ages of displacement events are equal to ages of rupture events for normal fault scarp which configuration was affected by ruptures only. Each rupture event is represented by one earthquake or by multiple earthquakes occurred in the short time span which cannot be resolved by the method. Due to aspects of the modeling, the position of colluvium level prior to displacement event is more convenient to operate than the normal component of slip. The same reference point and the position axis which were defined for sample position on the scarp are used for the position of colluvium level. The colluvium position prior to the event can be also calculated as the total cumulative slip less the cumulative slip prior to the event. The beginning of exposure is an age when the production of 36 Cl in the sample stripe becomes appreciable by the dating method. This modeling parameter cannot be certainly referred to any dating or geological characteristic of the fault scarp; the closest characteristic would be the exposure age of the fault scarp top.
The Random Walk algorithm of the code generates new exposure history by variation of each desired parameter according to a normal distribution. Parameter values of previous exposure history are used as means of normal distributions, while deviations are fixed through the run of the algorithm. The reduced chi-square value is calculated for each generated exposure history to indicate the goodness of fit of the modeled and measured 36 Cl profile: where Ne i is the measured 36 Cl concentration, atoms g À1 , Nm i is the modeled 36 Cl concentration, atoms g À1 , se i is the relative standard error of AMS measurement of 36 Cl concentration, sp i is the relative standard error of parent element concentrations, n is the number of samples, and k is the number of calculated parameters. Our code offers two types of conditions which can be set on chi-square value in the Random Walk mode: threshold of chi-square value and minimization of chi-square value. At the threshold condition, new exposure history is saved and accepted as the next point of algorithm if its chi-square value is less than or equal to the threshold. At the minimization condition, all generated exposure histories are saved, while new exposure history is accepted as the next point if its chi-square value satisfies following condition: where c 2 min is the minimal found chi-square value;c 2 new is the chi-square value of new exposure history; U is uniformly distributed value on the open interval (0,1); and a is the decision parameter.
The Monte-Carlo algorithm of the code randomly samples solution space bounded by uniform and/ or normal distributions of desired parameters. Uniform distribution is preferred when no independent information on a parameter value is available, while normal distribution is used to fix a parameter relying on independent data, for instance the age of the historical earthquake or weathering boundary.
Only chi-square threshold condition is available in Monte-Carlo mode. Exposure histories with chisquare value less than or equal to the given threshold are saved, while the rest are rejected. The Monte-Carlo algorithm of the code is less advanced than the Random Walk algorithm and requires longer time to find the best-fit solution.
Both Monte-Carlo or Random Walk algorithms cannot perform direct optimization of the number of events. Therefore, scenarios with different number of events should be considered to find the most probable exposure history. However, indirect optimization of the number of events is an inherent property of the code; extra events tend to converge through the simulation. The chi-square indicator could be insufficient to distinguish between scenarios with different number of events. To support it, the corrected Akaike information criterion (AICc) is calculated as follows, Such formulation of AICc differs from the standard definition of AICc by a constant, which can be omitted within analysis of one data set. As it seen from the formula above, AICc is more strict indicator to additional modeling parameters than chi-square. In our model the minimal value of AICc indicates the exposure history with the optimal number of events.