Modeling active fault systems and seismic events by using a fiber bundle model – example case: the Northridge aftershock sequence

Earthquake aftershocks display spatiotemporal correlations arising from their self-organized critical behavior. Dynamic deterministic modeling of aftershock series is challenging to carry out due to both the physical complexity and uncertainties related to the different parameters which govern the system. Nevertheless, numerical simulations with the help of stochastic models such as the fiber bundle model (FBM) allow the use of an analog of the physical model that produces a statistical behavior with many similarities to real series. FBMs are simple discrete element models that can be characterized by using few parameters. In this work, the aim is to present a new model based on FBM that includes geometrical characteristics of fault systems. In our model, the faults are not described with typical geometric measures such as dip, strike, and slip, but they are incorporated as weak regions in the model domain that could increase the likelihood to generate earthquakes. In order to analyze the sensitivity of the model to input parameters, a parametric study is carried out. Our analysis focuses on aftershock statistics in space, time, and magnitude domains. Moreover, we analyzed the synthetic aftershock sequences properties assuming initial load configurations and suitable conditions to propagate the rupture. As an example case, we have modeled a set of real active faults related to the Northridge, California, earthquake sequence. We compare the simulation results to statistical characteristics from the Northridge sequence determining which range of parameters in our FBM version reproduces the main features observed in real aftershock series. From the results obtained, we observe that two parameters related to the initial load configuration are determinant in obtaining realistic seismicity characteristics: (1) parameter P , which represents the initial probability order, and (2) parameter π , which is the percentage of load distributed to the neighboring cells. The results show that in order to reproduce statistical characteristics of the real sequence, larger πfrac values (0.85< πfrac < 0.95) and very low values of P (0.0< P ≤ 0.08) are needed. This implies the important corollary that a very small departure from an initial random load configuration (computed by P ), and also a large difference between the load transfer from on-fault segments than by off-faults (computed by πfrac), is required to initiate a rupture sequence which conforms to observed statistical properties such as the Gutenberg–Richter law, Omori law, and fractal dimension.


Introduction
Most earthquakes occur when adjacent blocks move along fractures in the Earth's crust, as a consequence of stress build-up arising from the regional strain and the stress change caused by a preceding earthquake or by the tectonic stress accumulation (Stein et al., 1994).These fractures, or faults, are discontinuous geological features consisting of a number of discrete segments (Segall and Pollard, 1980), which can be up to hundreds of kilometers in total length.Faults are the weakest parts of the crust and thus are more likely to release accumulated stresses, by means of slipping, than non-fractured crust.Earthquakes may occur many times on the same fault over millions of years.Known active faults have ruptured several times in the last 120 000 years and are considered likely to move again (Wallace, 1981).It has been observed that earthquakes are strongly correlated in space around the active fault systems (Kroll, 2012).
Fault systems have a statistical self-similar structure over a wide range of scales (Kagan and Knopoff, 1980;Sadovskiy et al., 1984;Hirata and Imoto, 1991) which can be described by means of fractal geometry, as introduced by Mandelbrot and Pignoni (1983).Geometrical fractal structures such as faults arise from self-organized criticality (SOC) phenomena over large temporal periods (Bak and Creutz, 1994).
SOC systems have been studied as a means to explain seismicity by several authors (Barriere and Turcotte, 1994;Scholz, 1991).In particular, models based on cellular automata have been used to describe SOC behavior in earthquake series (Aki, 1965;Barriere and Turcotte, 1994;Castellaro and Mulargia, 2001;Georgoudas et al., 2007;Adamatzky and Martínez, 2016).The fiber bundle model (FBM), a model based on cellular automata (Peirce, 1926;Daniels, 1945;Coleman, 1956), provides a conceptual and numerical description of the rupture process in heterogeneous media (Kun et al., 2006b).In this article, we present a model that improves over the base model presented in Monterrubio-Velasco et al. (2017) by projecting the geometrical fault systems in a FBM algorithm.We developed a novel and powerful simulation model with simple and straightforward input parameters that is capable of providing simulation scenarios that are statistically consistent with real cases.The model circumvents the high complexity related to the derivation of deterministic models of by including faults as an abstraction that try to explain weaker regions that are more likely to generate earthquakes.In this version of the model, the faults are represented as a projection on the surface as a 2-D approximation.The extension to 3-D is being considered for a later phase of the study.We test the FBM's ability to capture the statistics coming from empirical laws (i.e., the GR law, the fractal capacity dimension, the MO law, and the Hurst exponent) by means of a parametric study.As a case study, we consider the Northridge aftershock sequence (17 January 1994, M W = 6.7) along with the geometry of its fault system.The statistical characteristics from the Northridge aftershock sequence are compared with the statistics obtained in synthetic catalogs generated with our FBM.Even though we focus on one aftershock case, the applicability of our model can be extended to analyze other aftershock sequences using particular fault system as input configuration.For example, in Monterrubio-Velasco et al. (2018), this model was applied to classify, via machine learning algorithms, three different aftershock sequences.Lastly, we study the event productivity as function of on-fault (or weakling areas) and off-fault parameters.

Background
This section gives a general description of the fiber bundle model and the statistical relations used in this work.

Fiber bundle model: general description
The basic components necessary to construct an FBM are (Andersen et al., 1997;Phoenix and Beyerlein, 2000;Pradhan and Chakrabarti, 2003;Sornette, 1989;Kloster et al., 1997): 1. Defining a discrete set of cells located in a regular (ndimensional) lattice.For our purpose, this lattice will be two-dimensional, representing the geographical study area in an azimuthal view.
2. Assigning a probability distribution function for the inner properties of each cell.This failure law will define the probability distribution function of the stress (in the static case) or the probability distribution function of the rupture time (in the dynamic version) (Hansen et al., 2015).
3. Establishing the load-sharing rule.This component is crucial in a FBM since the model shows a fundamental change depending on the manner of the load transfer after a cell fail (Pradhan et al., 2010).The most used sharing rules are the equal or global load sharing (Turcotte et al., 2003) and the local sharing rule (LLS).This last sharing rule favors the stress concentrations and promotes nearest neighbors to reach a critical rupture state.
FBM was developed in two versions that simulate material rupture by different effects: (1) the static version in which the fiber strength is time independent (Vázquez-Prada et al., 1999;Kun et al., 2006a;Pradhan et al., 2010), and (2) the dynamic version where the study of the material rupture is a time-dependent process, such as stress rupture, creep rupture, static fatigue, or delayed rupture (Coleman, 1956;Moral et al., 2001b).In this research work, we will use a dynamic FBM with LLS in the probabilistic formulation developed in Moreno et al. (2001).

Probabilistic FBM
According to laboratory studies, the Weibull distribution describes the hazard rate (κ) on materials subjected to a constant load or stress, σ (Coleman, 1956;Moreno et al., 2001): where ν 0 is the hazard rate under stress σ 0 , and the ρ exponent is the Weibull index defined in the range 2 ≤ ρ ≤ 50 (Yewande et al., 2003;Kun et al., 2006a;Nanjo and Turcotte, 2005).In the present work, we use a dimensionless representation of quantities, so that we can use normalized stresses and Eq. ( 1) can be written as κ(σ ) = σ ρ , following Moreno et al. (2001) and Monterrubio-Velasco et al. (2017).Gómez et al. (1998) introduced a probabilistic approach as an alternative formulation to the dynamic FBM which we follow here.
The FBM model simulation starts by discretizing a hypothetical surface in a bi-dimensional array (N x × N y ).At the first step, the load of the cells in the bundle is assigned following a uniform distribution σ (x,y) = U [0, 1), x = 1, . .., N x , and y = 1, . .., N y .The notation U[0,1) is a mathematical form to represents a uniform distribution with values in the range equal to or greater than 0 and lower than 1.This assumption simulates an initial heterogeneity in the load cell properties.The hazard rate assigned at each cell is computed using Eq. ( 1) (Moreno et al., 1999;Pradhan et al., 2010).Furthermore, a rupture probability, F (x,y) , is computed for each individual element or cell, and at each step.This value is load dependent and is defined by where δ is the time interval (or inter-event time) until the next rupture occurs, valid for any load-sharing rule (Moral et al., 2001a), computed as (Moreno et al., 2001;Moral et al., 2001a) . (3) From Eq. (3), δ is a dimensionless quantity since σ (x,y) and ρ are dimensionless.
The time of occurrence or cumulative time T (k) is defined as where k = 1, . .., k max , being k max the maximum number of steps.Our model uses Eqs. ( 2), (3), and (4) to compute the rupture probability and the inter-event time.

Previous aftershock models
In our previous work, we developed a FBM version to simulate spatial and magnitude aftershock patterns.Following the general assumptions proposed in Correig et al. (1997), we modified it in order to define a computational domain that describes a particular geographical area.Moreover, the initial load values are ordered according to a probability value P .P = 0 represents a random spatial distribution of initial loads (heterogeneous), and P = 1 implies a homogeneous distribution in agreement with a proxy of Coulomb stress changes produced by a main event in the center of the computational domain (Monterrubio-Velasco et al., 2017).A local load-sharing rule including the eight nearest neighbors, and a threshold load (σ th ≡ 1) are established a priori.When the load in a cell exceeds this threshold load, a critical point is reached and a imminent rupture occurs (called avalanches).If more than one cell exceeds σ th , the cell to fail is that which exhibits the maximum load.The rupture algorithm is sequential, since at each time step one cell has to fail.The inter-event time (δ) must be updated at each discrete step as described in Eq. ( 3).After a cell fails, they distribute its respective percentage of load (σ T = σ F • π) to its eight neighbors.Perpendicular neighbors will receive the largest amount of load ((σ T • 0.98)/4), while diagonal neighbors get a load of (σ T •0.02)/4.The different weights were chosen considering what is expected for the maximum shear stress directions with respect to the main stress orientation that gives rise to both synthetic and antithetic faulting (e.g., Stein et al., 1994).In Monterrubio-Velasco et al. (2017), a large number of numerical experiments were carried out to define the appropriate values of these quantities.Varying the weights results in slight changes but as long as the main directional properties in weight are preserved the output does not change significantly.

Statistical and fractal relations
In order to quantify the resemblance between synthetic catalogs and real seismic catalogs, we use statistical measures which are relevant for evaluating the SOC behavior.These measures are represented by power laws in magnitude (GR law), time (MO law), and space (e.g., fractal dimension).In Appendix A1, we introduce these relations, describing their applicability, as well as the interpretation and the methods of quantification.Table 1 summarizes the characteristics, abbreviations, and usefulness of the empirical relations used in this work.

Extending the model
In our previous work (Monterrubio-Velasco et al., 2017), we did not incorporate the information of the fault geometry which is a fundamental property to describe a particular tectonic region.Therefore, the main contribution of this research work is extending the model by the addition of the fault system geometry by prescribing parameters that quantify "weakness" properties, i.e., the capacity to produce load concentrations that generate a rupture.In contrast to the classical definition that uses measures such as dip, strike, and slip, here, the faults are considered the weakest parts of the Earth's crust.The parameter π quantifies the load transfer and it controls the amount of load distributed from a failed cell to its neighbors.Since the seismic rupture is not conservative, the parameter π(x, y) defines the percentage of load lost at each discrete step.The output of the model is a synthetic catalog with statistical properties that change depending on the input parameters.

Algorithm
In the Appendix, three pseudo-codes are included to describe the model algorithm.Our model is coded in Julia language (Bezanson et al., 2017) for the sequential version and in Python language for the paralleled version of the loops.
-The initial conditions are as follows: In this work, and as a first attempt, we simplified the real 3-D domain by choosing a bi-dimensional surface to represent the epicentral distribution.Moreover, we assume that considering that the seismicity of southern California is shallow and mostly restricted to the planar strike-slip faults, the two-dimensional approach can be used as a simplification.Therefore, the 2-D Cartesian grid is a rectangular domain of N x × N y = N T square cells.The domain is a planar representation of the study area for which we define three values at each cell in x ∈ [1, . .., N x ] and y ∈ [1, . .., N y ] that assign their properties: a. σ (x,y) : a discrete value of load where the initial load distribution is taken from a uniform distribution function with values in the range [0,1).
b. π (x,y) : a load-transfer value that defines the percentage of the load distributed to the neighbors after a cell fails.Since the fault geometry is an abstraction that simulates weaker regions in the model domain, the values of π (x,y) that define the percentage of load distributed after a cell fails have either a background value (π back ) or a fault value (π frac ).
It was found that an appropriate value of this parameter must fall within the range 30 ≤ ρ ≤ 50.
e. P , the heterogeneity of the initial load distribution: this parameter has a global influence on the final scenario but we have used the best configuration as tested in Monterrubio-Velasco et al. (2017).
No external load is received after the initial load assignation, so that our model describes the relaxation process after a mainshock.Therefore, we do not discuss or simulate the mainshock or foreshocks.The load increase in a cell is due to internal load transfer processes.In a companion study (Monterrubio-Velasco et al., 2019), we present the TREMOL v0.1 code which studies the case of asperities.That model was developed as a mainshock simulator based on the FBM to analyze the case of Mexican subduction earthquakes.A version of our model that includes the effect of tectonic loading is still in progress.
-The rupture conditions are as follows: In this work, we choose the values proposed in Monterrubio-Velasco et al. (2017, 2019) in order to reduce the degrees of freedom.The only two parameters that change during the execution of the model are σ (x,y) and F (x,y) .Throughout the iterations of the simulation, we identify two possible outcomes: normal events, i.e., minor or background ruptures, and avalanches, i.e., a collection of spatiotemporally clustered events that result in large rupture (Monterrubio-Velasco et al., 2017).We remark that, in the present work, avalanches are actual secondary ruptures, whereas normal events are minor events of low magnitude that produce the rupture of a single cell.The consecutive rupture of avalanches events will produce a rupture cluster with a size determined by their area in (cell) units S(N A ).
-The completion of a simulation is as follows: A FBM simulation is terminated when any cell in the system is unable to exceed the threshold σ th .We have empirically determined that the total number of steps k max when this situation occurs is typically k max ≈ 3N T /4.Beyond this value, the system no longer generates loads that overpass σ th .Hence, we take this value as a terminal condition in our simulations.After the simulation is completed, we obtain a synthetic seismic catalog which includes the number of simulated earthquakes (N A ) with their corresponding area S(N A ) (number of events that produces a single rupture), occurrence time t A (Eq. 3), and their spatial location (x, y).In Sect.3.2, we discuss how the avalanche size S(N A ) can be converted into magnitude.It should be noted that two realizations with identical parameters result in different seismic catalogs due to the random component in the initial load.A detail procedure is explained in the pseudo-codes provided in Appendix A2.

Experimental procedure
The discrete planar faults of a particular region are modeled by using an image of the real fault system.This image is mapped in the domain (see the example in Fig. 1).A parametric study is employed to determine the best range of values to produce synthetic catalogs with appropriate statistical characteristics.In this work, we use the following parameters: ρ = 30, σ th = 1, σ D = 0.02/4, and σ N = 0.98/4.We also fix a background π back = 0.65 value for all non-faulted cells and assume a square grid, with same lateral size on the x axis and y axis, i.e., N x = N y = √ N T .These values are considered from the results obtained in Monterrubio-Velasco et al. (2017).The reasons to choose these values are summarized as follows: -30 < ρ ≤ 50 produces features observed in real aftershock time series (Monterrubio-Velasco, 2013;Monterrubio et al., 2015).
σ th = 1 is the upper bound of the uniform distribution and appears as a natural threshold.
-In Monterrubio-Velasco ( 2013), Monterrubio et al. (2015), and Monterrubio-Velasco et al. (2017), a range of 0.63 < π < 0.70 was determined experimentally to produce ruptures that mimic features of real catalogs without considering any difference between regions (i.e., no difference between faulting and non-faulting regions).We considered π back = 0.65 as a mid-range value to be assigned to background cells.
The map of faults has a real physical size in km 2 .So, after executing the Algorithm 1 in Sect.A2.1 (Appendix A2), in the post-processing analysis, we assign an area at each cell in km 2 , namely A cell .To compute the avalanche area S(N A ) in km 2 , we use the relation for j = 1, . .., N A .We computed an equivalent magnitude using the scale magnitude-area relation proposed by Hanks and Bakun (2008) in Eq. ( 6).
where A j , for j = 1, . .., N A is the rupture area expressed in km 2 and M W is the moment magnitude.This relation is specific for events in a crustal-plate-boundary tectonic regime (Stirling and Goded, 2012).
Careful attention has been given to minimum magnitudes which depend on size of the cell A cell , i.e., are proportional to  2018).In order to make results comparable for different grid sizes, we filter out events rupturing less than a minimum amount of cells for the finer grids.This is done because we want to show the influence of different grid sizes over the statistical results, and we want to compare our results with real cases where the smallest magnitudes are not resolved by the seismograph networks.
The epicentral location of the simulated aftershocks is the position of the first avalanche event ((E x (j ), E y (j ))) in the cluster.
We define (r) as series of the Euclidean distance between two consecutive epicenters and τ (t) the inter-event time series corresponding to two consecutive avalanches (see Eq. 4).
Table 2 defines the model parameters and provides optimal search ranges.

Test case: Northridge aftershock sequence
In order to validate and compare our synthetic seismic catalogs with real seismicity, we modeled as an example case the fault system geometry and the seismic properties of the Northridge aftershock sequence (Fig. 2).The Northridge   (Thio and Kanamori, 1996).The mainshock occurred due to the rupture of a previously unrecognized blind reverse fault with a moderate southward dip (Savage and Svarc, 2010; Scientists of the U.S. Geological Survey, 1994).The Northridge mainshock occurred at a fault belonging partly to a large fault system in the transverse ranges.This fault system is under compression in the NNW direction related to the "big bend" of the San Andreas Fault (Norris and Webb, 1990;Hauksson et al., 1995).The Northridge earthquake was followed by a sequence of aftershocks, between 17 January and 30 September 1999, including eight aftershocks of magnitude M W ≥ 5 and 48 of 4 ≤ M W ≤ 5. We computed the statistical parameters of the aftershocks using the data recorded by the Southern California Seismic Network.We consider the aftershock time period from the mainshock (17 January 1994 at 04:31 UTC) until 1 year later (17 January 1995).In space, we consider events that occur in a square area of 0.6 • × 0.6 • taking as center the mainshock epicenter (Turcotte, 1997) In Fig. 2, we show the spatial representation of the faults analyzed in this work, and the seismicity of magnitude larger than 2.0 during 1981-2006.We also assumed that the faulting area is the same, independently of the model domain size .However, the number of cells modify the size of each cell, resulting in 0.077, 0.043, and 0.027 km 2 for N x = [180, 240, 300] cells, respectively.
In Table 3, we show the statistical parameters for the Northridge sequence computed for different threshold magnitudes from M min = 1.5 to M min = 3.5.These values will be used as a reference to determine the sets of FBM param-  eters that best reproduce the Northridge statistics.Figures 3  and 4 show the fitting of the GR law, MO law, and Hurst exponents (H , H τ , and H Mag ) for different minimum magnitudes (M min > 1.5, > 2.0, > 2.5, > 3.0, and > 3.5).We note that Wiemer and Wyss (2000) calculated the minimum magnitude of completeness in the Los Angeles area as M c ≈ 1.5.

Results
We divide the results and their analysis in three domains: space: fractal capacity dimension, D 0 , for the epicentral distribution of synthetic earthquakes, and Hurst exponent for epicentral distance between consecutive synthetic earthquakes, H ( ); time: inter-event times H (τ ) and MO parameters (p, c, K).

Parametric analysis over the synthetic series
For each parameter, we list its observed properties to facilitate the reading.We are analyzing simultaneously three parameters (the size of the domain N measured in cells units, the initial order configuration P , and the percentage of transfer load in a fault-on cell π frac ).However, it is worth noting that they are coupled between them.Then, we analyze how these three parameters are modifying the statistical and the fractal properties of synthetic catalogs.

Fractal dimensions of synthetic catalogs
The first analysis is related to the fractal capacity dimension, D 0 (Sect.A1.1), shown in Fig. 5.
-As N increases, D 0 becomes more sensitive to P (where P is the initial load values ordered according to a probability value), because the larger the area, the more abundant and scattered the events, and the effect of P over the simulated events increases.The effect of an increase in P is a reduction of D 0 .On the other hand, as P tends to zero, the events are more scattered because randomness in the initial load causes sparsity in the event distribution.
π frac does not seem to have a large impact in D 0 across our experimental range.
-For P ≤ 0.16 and π frac > 0.9 (for all N values), our FBM simulations yield a D 0 compatible with that computed for the Northridge series with a M min ≈ 2 (Table 3).

The b value, mean, and maximum magnitude of synthetic catalogs
The b value is clearly influenced by all three parameters (N , P , and π frac ), as is shown in Fig. 6.
- The difference in the b values as a function of N might be due to two possible causes: 1. the number of events included in the statistical fit and 2. the size of the earthquake-simulated avalanches, S(N A ).
In fact, related to this observation in Fig. 7, we observe that as P increases, the maximum magnitude also increases, also modifying the frequency-magnitude distribution and the b value.
As an example, Fig. 8 shows the frequency-magnitude distribution for P = [0, 0.08, 0.16, 0.24, 0.32, 0.38], π frac = 0.9, π frac = 0.65, and N = 300.We observe that as P increases, the productivity of intermediate size events decreases, and the maximum magnitude increases.Large P values imply that the probability to find cells with large loads clustered in the middle increases (Monterrubio-Velasco et al., 2017).So, in this condition, it is more likely to generate larger earthquakes.This behavior is similar to that observed in the characteristic earthquake distribution (Wesnousky, 1994).
The results of the mean and maximum magnitudes are depicted in Fig. 7.  -From Fig. 7a, we observe that the mean magnitude is independent of π frac and to a lesser extent of P .It is worth pointing out that the b values are similar to those shown in Table 3 when M min ≈ 2 or 2.5.
-Fig.7b shows that, in our model, an aftershock with a magnitude such as that of the largest Northridge aftershock magnitude (M W = 5.9) is obtained for a nonunique combination of parameters.When P > 0.08 and π frac > 0.7, this magnitude is overestimated.Given that the largest aftershock has a magnitude M W = 5.9, larger-magnitude values are not described in the series.
From Figs. 6 and 7, we observe that the best range of b values, similar to that obtained for Northridge sequence, is for P < 0.16, π frac > 0.90, and N ≥ 240.-The re-scaled range analysis of the (r) series reveals their independence on N and π frac but shows a slightly higher dependence with P .In general, as P increases, H ( ) also increases.As P decreases, H ( ) → 0.5, implying that the system tends to a random behavior of the inter-event epicentral distribution.However, H ( ) always remains similar to the values of Table 3 for P ≈ 0. Gkarlaouni et al. (2017) showed that the seismicity in the Corinth rift (Greece) corresponded to H ( ) → 0.5 as the threshold magnitude decreased.

Hurst exponent of synthetic catalogs
-The analysis of H (τ ) reveals in general values > 0.5, which implies a persistence in the dynamic system of inter-event times; i.e., the behavior of future inter-event time can be extrapolated from previous behavior.As P increases, the persistence of H (τ ) also increases, which may be related to the MO trend (discussed below).The influence of π frac over H (τ ) is not clear; however, the number of events decreases when we take a larger cutoff magnitude, and this fact could affect the re-scaled range statistics.

MO parameters in synthetic catalogs
The MO empirical law requires careful analysis in our FBM implementation.First of all, we must take into account that we are using dimensionless time (see Eq. 3).Our cumulative time, T j , is computed using Eq.(4).For example, considering the input values of P = 0.08, π frac = 0.90, and N x = 180, if we include all simulated events (minor events and synthetic aftershocks) together (e.g., Fig. 10a), we get a satisfactory MO fit, with p = 1.1, c = 1.5, K = 11 991.3, and rms = 675.0.But if we use only the time occurrence of the simulated aftershocks (S(N A )), the parameters depart from the expected trend (e.g., Fig. 10b), obtaining p = 1.9, c = 4.9, K = 1571.9,and rms = 10.8.From Fig. 10b, two regions are distinguished by the density of events, since for the first interval of time the density is larger (blue) than for the subsequent region (green).As a consequence, the MO fit is deviated adjusting for events in the blue region.Real aftershocks do not always follow a single MO decay trend (Utsu and Ogata, 1995).Moreno et al. (2001) developed an alternative model to understand this phenomenon called leading and cascade (LA-CAS) events.This model proposes a separation of earthquakes in two groups: one that strictly follows the MO hypothesis (called leading aftershocks; LAs) and those called cascades, which are the events that occur between two consecutive LAs.Note that to obey the MO relation, the inter-event times must increase monotonically.In our previous work (Monterrubio et al., 2015), we tested the LA-CAS algorithm to study the temporal behavior of three real aftershock sequences.We applied the LA-CAS algorithm to the synthetic earthquakes series, as illustrated in Fig. 10b.After segregating the events in LA and CAS, we obtain a better fit to the MO relation for the LA series (Fig. 11, LA: p = 0.9, c = 0.1, K = 5.1, and rms = 1.2).The full parametric results computed from the MO relations are shown in Fig. 12a and  b.
The MO parametric results provide information of the temporal behavior in the simulated series.We observe that P < 0.08 implies p and c values close to the expected Northridge values (see Table 3).However, after segregating LA and CAS, the number of events decreases, so the value of K is much lower than expected.These results indicate that series with an initial load configuration organized with a probability P > 0.08 depart from a typical MO behavior.This occurs because the aftershock series produced with larger initial organization probabilities (P > 0.08) tend to generate temporal series with very short elapsed times (Eq.3), and larger avalanche clusters, which does not follow a typical MO distribution.The closest behavior to the observed MO parameters of Northridge (Table 3) occurs for P ≤ 0.08 and π frac > 0.7.

Trigger and shadow regions
The load-transfer value π(x, y) is highly relevant to reproduce temporal, magnitude, and spatial patterns of real series (Monterrubio et al., 2015).Considering this fact, we are also interested in studying the implication of this value in the aftershock productivity, in particular for off-fault regions.To test the productivity as a function of π bkg , two extreme values are considered for π bkg = 0.25 and π frac = 0.65.In Fig. 13, we observed that activity in background (non-faulting) cells largely decreases for small π bkg values.This occurs because for low π values, the probability to produce an event with larger ruptured area decreases (Fig. 14).The results suggests that variations in π bkg for different regions of the domain can lead to producing shadow and triggered regions, giving a scenario closer to a real case (King et al., 1994;Stein, 1999;Hainzl et al., 2014).The lack of the depth (3-D) in our model could produce bias in the shadow region interpretation and is a limitation to a closer description of the phenomenon.However, in this study, our first attempt is more focused on the parametric implications of the fault regions included in the model as "weak" areas.We expect that the integration of triggered and shadow regions will be plausible in future implementations to improve the results.

Results summary: synthetic catalogs
Lastly, we estimated the error between the real and synthetic statistical values using a measure similar to the Euclidean distance r E−M min .However, in our case, r E−M min is referred to a normalized parametric space because of the different www.solid-earth.net/10/1519/2019/units in the parameters, and it is defined as where ps = [D 0 , M W , b value, M max , M min , H ( ), H (τ ), p, c] is the vector that contains the values of the series generated with a given set of input parameters P , π frac , and N. Similarly, the vector ps M min contains the values for the Northridge series considering four different minimum mag-nitudes M min = [1.5, 2.0, 2.5, 3.0] (Table 3).We computed r E−M min for the 162 combinations of P , π frac , and N, each one with three realizations.In Table 4, we show the minimum of r E−M min for M min = [1.5, 2.0, 2.5, 3.0].The minimum Euclidean distance occurs when we consider the NOR series with M min = 2.0.It is worth mentioning that the minimum magnitude of the synthetic aftershocks considered in this work is also ≈ 2.0.The results show that the most appropriate set of parameters to model this data series is P = 0, π frac = 0.90 using N = 300 cells.Figure 15 shows the Eu- Table 4. Minimum Euclidean distance r E−M min (Eq.7) using four different M min NOR series (Table 3).
M min r E−M min N P π frac 1.5 0.78 240 0.16 0.6 2.0 0.63 300 0 0.95 2.5 0.88 300 0.08 0.9 3.0 1.53 300 0.16 0.7 clidean distance between the sequence generated by NOR using the set of parameters obtained with M min = 2.0 and its real values.In Fig. 16, we show an example of a spatial distribution of events and its related GR relation, using the set P = 0, π frac = 0.95, and N = 300.As shown in Fig. 16, the largest aftershocks have its epicenter on the fault's cells (M > 3.5).The epicenters are depicted with a blue star.This scatter plot also shows that the smallest events usually occur spread out.The relation of the magnitude and the cumulative number of events, generated in this example, shows a GR fit with a similar b value to that computed in Table 3.

Discussion
The main goal of this study is to integrate prior knowledge of the spatial geometry of faults in the implementation of the FBM algorithm, improving the model previously proposed in Monterrubio-Velasco et al. (2017).As it was pointed out, to introduce the fault system geometry we assume some cells to be weaker than the rest representing faults in the bi-dimensional array.This "weakness" is assigned by a single parameter called π frac .The lack of the depth dimension may leave out information about the full phenomenon.However, it is a first attempt to test the FBM as an aftershock simulator including complex features.The main advantage of the present model is the reduced number of equations to be solved in comparison with deterministic models for similar purposes and the low number of parameters used to describe the model dynamics (π frac , π bkg , ρ, P , and N T ).
To validate our model, we used as an example the geome-  try of the Northridge fault system and the statistics of the aftershocks.Note that this model version describes the relaxation process after a mainshock.Therefore, we do not discuss or simulate the mainshock or foreshocks.In particular, we explored the power laws' exponents (D 0 , b, H, p parameters) in relation to the model parameters (Sect.A1).
Other models have been proposed to describe, with a simplified mechanism, the statistics of earthquakes, such as the two-fractal overlap model (Bhattacharya et al., 2009) or the Olami-Feder-Christensen (OFC) model (Olami et al., 1992).
In particular, the OFC model has a very similar algorithmic to our proposed model (Kawamura et al., 2012), but our model yields similar results with fewer input parameters and it is simpler to implement.As a statistical modeling tool, we need a parametric analysis to properly fit observational data.In our study, we have searched the range of values that generate synthetic series capable of reproducing the statistical relations of real aftershock series.In particular, we explored three (π frac , P and N) of the five free parameters to quantify their leading role in the model.We point out that π bkg and ρ are assumed as constants following results in Monterrubio et al. (2015) and Monterrubio-Velasco et al. (2017).In agreement with Monterrubio et al. (2015), we also confirm that the transferred load value π is the most critical parameter in order to reproduce temporal, magnitude, and spatial patterns of real series.Our results also suggest that variations in π for different regions of the domain might generate shadow regions (King et al., 1994;Stein, 1999;Hainzl et al., 2014).The initial load configuration, controlled by P , results as a determinant to describe the final statistical features in the model.In particular, the results indicate that P and π frac are inversely proportional.As we increase π frac , a small value of P is required to reproduce aftershock statistics.If the fault geometry is not considered in the model (π bkg = π frac ), the particular range of 0.60 < π < 0.70 found in Monterrubio-Velasco et al. (2017) is required to capture statistical patterns.
The results are sensitive to the size of the domain.An exhaustive parametric analysis using machine learning techniques to classify the synthetic series as function of the input parameters (the size N, P , and π frac ) was carried out in Monterrubio-Velasco et al. (2018).In Fig. 17 from Monterrubio-Velasco et al. (2018), we show the mean error of three different ML classification algorithms (random forests, supported vector machine, and flexible discriminant analysis) as a function of the domain (grid) size.The figure shows that as the grid size is increased, the classification error decreases, meaning that large grid sizes allow us to distinguish among the different properties.In other words, for small grid size, the difference is indistinguishable, while larger grid sizes are able to capture the differences.We observe the results using as classification two input parameters: P (in red) and π frac (in blue).When we use the P parameter, we observe that the size domain has to increase in order to reduce the mean classification error, and it becomes the minimum for N ≥ 300.On the other hand, if we want to classify the synthetic catalogs considering π frac , the figure shows that the error classification reaches a minimum value for lower grid sizes N ≥ 200.So, if we consider the case of P = 0, and the classification is based on π frac , then a proper grid size used to model aftershocks, including faults, is N ≥ 200.We can confirm that an optimization of the parametric search using classification machine learning techniques can be very useful in this stochastic model.
Considering the example of Northridge, our results suggest that the best combination of parameters to approximate to real cases depends on the minimum magnitude of the real catalogs, as shown in Table 4. Related to the completeness magnitude, Davidsen and Baiesi (2016) define the shortterm aftershock incompleteness (STAI) as a phenomenon that arises from overlapping waveforms and/or detector saturation, such as events that are missed in the coda of preceding ones.One important consequence of STAI is an increase in the local magnitude of completeness, since small events are not well recorded.It is worth noting that in this work we are not analyzing the STAI phenomena because we are not explicitly modeling this process.We use the Northridge catalog obtained by the Southern California Seismic Network (SCSN), and we analyze it as a "final" catalog.In our statistics and analysis applied to the real catalog, we consider different magnitude cutoffs, as shown in Table 3.The cutoff magnitude is not related to the time.On the other hand, it is noteworthy that our model is not affected by the STAI, because this phenomenon arises from overlapping waveforms, www.solid-earth.net/10/1519/2019/Solid Earth, 10, 1519-1540, 2019     and in our approach we are not considering this physical process.To modify the minimum magnitude in the synthetic catalogs, we only filter the events with small rupture areas.
The usefulness of this stochastic model is its capability to generate a large number of scenarios with statistical properties similar to real cases, with low computational cost and a low number of free parameters.

Conclusions
We present a novel model simulation of aftershock sequences that incorporates a 2-D spatial distribution of faults.The representation of faults is carried out by assigning weak cells embedded in a background of "normal" cells.However, this model fulfills statistical properties of aftershock when it is well tuned.We choose statistical relations which describe the aftershocks' behavior in space, magnitude, and time.By means of a parametric study, we have found the range of values that generates synthetic series capable of reproducing the statistical relations of real aftershock events.In particular, we have used the Northridge fault system geometry projected on the surface and its aftershock sequence as a study case.We conclude that the initial load configuration (quantified by parameter P ), which specifies the randomness in the background load distribution and the ratio of transferred load for a faulting cell π frac are the key parameters that control the earthquake's statistical patterns in FBM-simulated events.Moreover, these parameters are complementary; i.e., in absence of fault geometry information (π frac = π back ), values in the range 0.08 < P < 0.32 ensure statistical compatibility with real aftershocks.In particular, for π frac = π and N = 180, we recover the results obtained previously without fault information (Monterrubio-Velasco et al., 2017).On the other hand, when fault geometry is available, as in the case of the Northridge fault system, the results obtained in this work show that, in order to reproduce statistical characteristics of the real sequence, larger π frac values (0.85 < π frac < 0.95) and very low values of P (0.0 < P ≤ 0.08) are needed.This implies the important corollary that a very small departure from an initial random load configuration is required to initiate a rupture sequence which conforms to observed statistical properties such as the Gutenberg-Richter law, Omori law, and fractal dimension.In summary, the proposed model is a useful tool to model aftershock scenarios by means of its inherent statistical patterns in time, space, and magnitude.Moreover, the model circumvents the high complexity related to the derivation of deterministic models of earthquake rupture phenomena.We demonstrated the ability of the model to simulate statistically consistent data with those of real scenarios using only a few input parameters.Our model can be an alternative to the study of the complex behavior of earthquakes.Future work will focus on optimization of the parametric search using machine learning techniques and extensions towards a 3-D FBM version that incorporates the depth dimension into the model.

Figure 1 .
Figure 1.Map (azimuthal view) of the Northridge fault system.It is digitalized to include in our computing domain where each pixel of the map corresponds to a cell (x, y) in our array.

Figure 2 .
Figure 2. Map that includes the seismicity of magnitude larger than 2.0 during 1981-2006.The yellow star indicates the Northridge epicenter (M W 6.7, 1994).Red lines depict the faults of this region considering an approximated area of 0.6 • × 0.6 • .Blue circles indicate aftershock locations and the size of their magnitude.
The smallest array of N = 180 cells tends to produce similar b values, independently of P and π frac .On the other hand, as N increases, the b value is more sensitive to P and π frac .-For values of N = 240 and 300 cells, the b value shows a clear dependency on π frac .In these finer grids (N = 240), we consider an initial random load distribution P = 0, we observe that the synthetic b approaches b values computed for the Northridge sequence, if π frac > 0.90.This last observation is important to justify the influence of include fractures regions in our model.Because π frac > 0.90, it indicates a clear "nonconservative" properties in the fault-on cells, and under this assumptions b value is closer to the expected value computed for Northridge.

Figure a and
Figure a and b show the results of the Hurst exponent for inter-event distance H ( ) and inter-event time H (τ ) (Eq.A6).

Figure 6 .
Figure 6.The (a) b value (Eq.A8) and (b) the b-value error (Eq.A9) computed for different synthetic series considering three input parameters (N , P , π frac ).From top to bottom, N = 180, N = 240, and N = 300 (cells per lateral size in the domain).

Figure 10 .
Figure 10.Modified Omori relation (Eq.A11) (dashed red line) fitting for one simulation with P = 0.08, π frac = 0.90, N = 180.(a)Synthetic aftershocks series for M min = 2.0 (black dots).Two regions are distinguished, where the density of events in the blue region is much larger than that in the green.(b) All events generated (minor normal and avalanches) (in black dots).A high and homogeneous data density along the plot is observed.

Figure 11 .
Figure 11.MO fitting for the leading aftershock (LA) series computed for the same synthetic series of Fig. 10.

Figure 12 .
Figure 12.Modified Omori parameters computed in the LA synthetic series: (a) p, (b) c, (c) K, and (d) rms.Each marker is obtained by different synthetic series considering three input parameters (N , P , π frac ).From top to bottom, N = 180, N = 240, and N = 300 (cells per lateral size in the domain).

Figure 13 .
Figure 13.Epicentral spatial distribution for two examples: (a) P = 0, π frac = 0.95, N = 180, and π bkg = 0.25; and (b) P = 0, π frac = 0.95, N = 180, and π bkg = 0.65.The different circles sizes represent the equivalent magnitude according to the legend.The center of each circle is the epicentral or nucleation point of each synthetic event.The major events occur over the faults.

Figure 17 .
Figure 17.Mean error of three different ML classification algorithms (random forests, supported vector machine, and flexible discriminant analysis) as a function of the domain size (figure from Monterrubio-Velasco et al., 2018).

Table 1 .
Empirical relations, interpretation, and main parameters.

Table 2 .
Model parameters.The earthquake shook the San Fernando Valley, which is 31 km northwest of Los Angeles, near the community of Northridge.This earthquake is the largest recorded in the Los Angeles metropolitan area in the last century.The depth of the hypocenter was 18 ± 1 km.The seismic moment, M o , was estimated at 1.58 × 10 23 N m, with a stress drop of 27 MPa

Table 3 .
Statistical parameters of the real catalog of Northridge aftershocks using different threshold magnitudes (M min ).