A stochastic rupture earthquake code based on the fiber bundle model (TREMOL v0.1): application to Mexican subduction earthquakes

In general terms, earthquakes are the result of brittle failure within the heterogeneous crust of the Earth. However, the rupture process of a heterogeneous material is a complex physical problem that is difficult to model deterministically due to numerous parameters and physical conditions, which are largely unknown. Considering the variability within the parameterization, it is necessary to analyze earthquakes by means of different approaches. Computational physics may offer alternative ways to study brittle rock failure by generating synthetic seismic data based on physical and statistical models and through the use of only few free parameters. The fiber bundle model (FBM) is a stochastic discrete model of material failure, which is able to describe complex rupture processes in heterogeneous materials. In this article, we present a computer code called the stochasTic Rupture Earthquake MOdeL, TREMOL. This code is based on the principle of the FBM to investigate the rupture process of asperities on the earthquake rupture surface. In order to validate TREMOL, we carried out a parametric study to identify the best parameter configuration while minimizing computational efforts. As test cases, we applied the final configuration to 10 Mexican subduction zone earthquakes in order to compare the synthetic results by TREMOL with seismological observations. According to our results, TREMOL is able to model the rupture of an asperity that is essentially defined by two basic dimensions: (1) the size of the fault plane and (2) the size of the maximum asperity within the fault plane. Based on these data and few additional parameters, TREMOL is able to generate numerous earthquakes as well as a maximum magnitude for different scenarios within a reasonable error range. The simulated earthquake magnitudes are of the same order as the real earthquakes. Thus, TREMOL can be used to analyze the behavior of a single asperity or a group of asperities since TREMOL considers the maximum magnitude occurring on a fault plane as a function of the size of the asperity. TREMOL is a simple and flexible model that allows its users to investigate the role of the initial stress configuration and the dimensions and material properties of seismic asperities. Although various assumptions and simplifications are included in the model, we show that TREMOL can be a powerful tool to deliver promising new insights into earthquake rupture processes.

1 Introduction

Rupture models of large earthquakes suggest significant heterogeneity in slip and moment release over the fault plane (e.g., Aochi and Ide, 2011). In order to characterize the seismic source rupture complexity, two main models have been proposed: the asperity model (Kanamori and Stewart, 1978) and the barrier model (Das and Aki, 1977). Asperities are defined as regions on the fault rupture plane that have larger slip and strength in comparison to the average values on the fault plane (Somerville et al., 1999). Asperities also have larger stress drop than the background area (Madariaga, 1979; Das and Kostrov, 1986). Understanding the physical features in the fault zone that produce these high-slip regions is still a challenge.

The most common method for studying seismic asperities is waveform slip inversion. However, information obtained from this method is highly variable due to the inherent nature of the inversion process (see review in Scholz, 2018). The slip inversion results depend on the type of data (such as strong ground motion and geodetic and/or seismic data at different distances) and the inversion technique used. Somerville et al. (1999) used average slip to define asperities. In their criterion, asperities include fault elements for which slip is 1.5 times or more larger than the average slip. By using this criterion, it is possible to estimate the asperity area from a finite-fault slip model. Considering the stress drop for a circular crack model (Δσ) (Eshelby, 1957), the stress drop on an asperity (Δσ_a) can be estimated as $\mathrm{\Delta }{\mathit{\sigma }}_{\mathrm{a}}=(A_{\mathrm{eff}}/{\mathit{A}}_{\mathrm{a}})\mathrm{\Delta }\mathit{\sigma }$, where A_eff and A_a are the rupture effective area and the asperity area, respectively (Madariaga, 1979). The A_eff∕A_a factor (or its reciprocal value) depends on different features with the most relevant one being the type of earthquake. For example, Somerville et al. (1999) found that on average the total area covered by asperities represents 22 % of the total rupture area for inland crustal events. Murotani et al. (2008) showed that A_a∕A_eff is approximately equal to 20 % for plate-boundary events. Similarly, for subduction events, the value of A_a∕A_eff is approximately equal to 25 % (Somerville et al., 2002; Rodríguez-Pérez and Ottemöller, 2013). The previous average values were determined considering values that range from 0.09 to 0.35. This last condition means, for instance, that the reciprocal fraction A_a∕A_eff can deviate from these average values as well (for example, 0.09 to 0.35 for the proportions mentioned above), which leads to great stress contrasts (factors of 2.8 to 11) (Iwata and Asano, 2011; Murotani et al., 2008). Mai et al. (2005) proposed another definition of asperities based on the maximum displacement, D_max. They defined “large-slip” and “very-large-slip” asperities as regions where the slip D lies between $\mathrm{0.33}{D}_{\max }\le D<\mathrm{0.66}{D}_{\max }$ and 0.66D_max≤D, respectively. They found that approximately 28 % of the rupture plane is occupied by large-slip asperities, whereas very-large-slip areas constitute only 7 % of the fault plane. Furthermore, different authors agree that the rupture area of the asperity scales with the seismic magnitude (Somerville et al., 1999; Murotani et al., 2008; Iwata and Asano, 2011; Rodríguez-Pérez and Ottemöller, 2013, among others). The estimation of seismic magnitude is an essential feature for characterizing the energy of an earthquake. In fact, an accurate magnitude estimation is indispensable to conduct both deterministic and probabilistic seismic hazard assessments.

Earthquakes are the most relevant example of self-organized criticality (SOC) (Bak and Tang, 1989; Olami et al., 1992). The concept of SOC can be visualized by imagining a natural system in a marginally stable state, wherein phases of instability may occur that place the system back into a meta-stable state (Barriere and Turcotte, 1994). A popular model representing this process was proposed by Bak and Tang (1989) and is well-known as the “sand pile model”. Some models have been proposed to explain the statistical behavior of earthquake patterns based on the SOC concept: e.g., Caruso et al. (2007), Barriere and Turcotte (1994), Olami et al. (1992), and Bak and Tang (1989). The failure properties of solids have been modeled by simple stochastic discrete models, which are based on the SOC framework. The fiber bundle model, FBM, is one of those models that has been used to reproduce many basic properties of the failure dynamic within solids (Chakrabarti and Benguigui, 1997). Additionally, the FBM has been successfully applied to studies of brittle failure of rocks (Hansen et al., 2015; Monterrubio et al., 2015; Turcotte and Glasscoe, 2004; Moreno et al., 2001).

2 The fiber bundle model

The FBM is a mathematical tool to study the rupture process of heterogeneous materials that was originally introduced by Peirce (1926). Over the years the FBM has been widely used to study failure in a wide range of heterogeneous materials (Hansen et al., 2015; Pradhan and Chakrabarti, 2003). Regardless of the specific FBM type, there are three basic assumptions that all FBMs have in common (Daniels, 1945; Andersen et al., 1997; Kloster et al., 1997; Vázquez-Prada et al., 1999; Phoenix and Beyerlein, 2000; Pradhan et al., 2010; Monterrubio-Velasco et al., 2017).

1. A discrete set of cells (or fibers) is defined on a d-dimensional lattice. In seismology, the bundle can represent a fault system or seismic source wherein each fiber is a section of the fault plane (Moreno et al., 2001) or individual faults (Lee and Sornette, 2000).

2. A probability distribution defines the inner properties of each cell (fiber), such as lifetime or stress distribution.

3. A load-transfer rule determines how the load is distributed from the ruptured cell to its neighbor cells. The most common load-transfer rules are (a) equal load sharing (ELS), in which the distributed load is equally shared to the other cells within the material or bundle, and (b) local load sharing (LLS) whereby the transferred load is only shared with the nearest neighbors.

TREMOL is based on the probabilistic formulation of the FBM, with the failure rate of a set of fibers given by Eq. (1): (Gómez et al., 1998; Moral et al., 2001).

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}U\left(t\right)}{\mathrm{d}t}}=-U\left(t\right)K\left(\mathit{\sigma }\right(t\left)\right),\end{array}$$

where U(t) is the number of fibers that remain unbroken at time t. The hazard rate K(σ(t)) is a function of the fiber stress σ(t). Experimental results show that the hazard rate of materials under constant load can be well-described by the Weibull probability distribution function. This behavior can be summarized in Eq. (2) (Coleman, 1958; Phoenix, 1978; Phoenix and Tierney, 1983; Vázquez-Prada et al., 1999; Moreno et al., 2001; Biswas et al., 2015):

$$\begin{array}{}\text{(2)}& K\left(\mathit{\sigma }\right(t\left)\right)={\mathit{\nu }}_{\mathrm{0}}{\left({\displaystyle \frac{\mathit{\sigma }\left(t\right)}{{\mathit{\sigma }}_{\mathrm{0}}}}\right)}^{\mathit{\rho }},\end{array}$$

where ν₀ is the reference hazard rate, and σ₀ the reference stress. The Weibull exponent, ρ, quantifies the nonlinearity (Yewande et al., 2003). If ${\mathit{\sigma }}_{\mathrm{0}}={\mathit{\nu }}_{\mathrm{0}}=\mathrm{1}$, the expression in Eq. (2) can be simplified to K(σ(t))=σ(t)^ρ. From the probabilistic formulation, two equations arise (Eqs. 3 and 4), which are applied in our algorithm to define the system dynamics. The details of these two equations are described below.

• a)

  Gómez et al. (1998) and Moral et al. (2001) developed a relation to compute the expected rupture time (dimensionless) of the fibers following Eqs. (1) and (2). This expected rupture time interval is defined as δ_k (Eq. 3) and can be applied to any load-transfer rule:

  $$\begin{array}{}\text{(3)}& {\mathit{\delta }}_k={\displaystyle \frac{\mathrm{1}}{{\displaystyle }\sum _{i=\mathrm{1}}^N{\mathit{\sigma }}_i^{\mathit{\rho }}\left(t\right)}},\end{array}$$

  where N is the total number of cells, and σ_i is the load in the ith cell. The dimensionless cumulative time, T, is the sum of δ_k.

• b)

  The failure probability, F_i, which is a function of the load σ_i in each cell, is (Moreno et al., 2001)

  $$\begin{array}{}\text{(4)}& F_i={\mathit{\delta }}_k{\mathit{\sigma }}_i^{\mathit{\rho }}\left(t\right).\end{array}$$

  The dynamic values δ_k and F_i are updated with each time step due to rupture processes and the resulting load transfer.

A suitable FBM algorithm to simulate earthquakes should consider a complex stress field, physical properties of materials, stress transfer between faults (at short and long distances), and dissipative effects. Using the FBM we assume that earthquakes can be considered analogous to characteristic brittle rupture of a heterogeneous material (Kun et al., 2006a, b).

The previous basic concepts about the FBM were considered for the development of the TREMOL code, with the purpose of modeling the behavior of seismic asperities. In the next section, we describe the details of this code.

3 The TREMOL code

Since the main objective of TREMOL is to simulate the rupture process of seismic asperities based on the principles of the FBM, we model two materials with different mechanical properties interacting with each other.

In order to introduce the features of TREMOL we describe three main stages during the application of TREMOL.

• 1.

  Preprocessing
  In this stage we have to assign the following input data:

  • -

    the size of the fault plane,

  • -

    the size of the maximum asperity within the fault plane, and

  • -

    other parameters (load-transfer value π, strength value γ, initial load values σ, and load threshold σ_th).

• 2.

  Processing
  TREMOL uses the data from the preprocessing stage to carry out the FBM algorithm, and by applying Eqs. (3) and (4) the rupture process is computed in the fault plane studied. The asperity size of each earthquake is used by TREMOL to also compute the magnitude of each synthetic earthquake.

• 3.

  Post-processing
  In this stage, TREMOL summarizes the results that are computed in the processing stage and computes the equivalent rupture area (km²). In general, TREMOL output generates a synthetic catalog of earthquakes, which consists of the following:

  • -

    total number of earthquakes that can occur in the fault plane studied,

  • -

    size of the asperity of each earthquake, and

  • -

    magnitude of each earthquake.

In the next sections we describe with more detail each one of the three main stages during the application of TREMOL. An overview of the entire simulation process is shown in Fig. 1.

Figure 1 TREMOL flowchart. At the beginning (preprocess) the algorithm initiates a domain Ω with N_x×N_y cells in which every cell is either part of an asperity or of the background or fault plane. Afterwards (first time step, k=1) a uniform distribution allocates a random stress load and rupture probability to all cells. In addition, asperity cells obtain a random strength value from a uniform distribution. Next (time step k≥2) the failure process starts following the FBM algorithm. After every failure the stress of the broken cell is redistributed via the LLS rule and the number of time steps (k) increases by 1 until the target number of time steps is reached. If the final number of time steps has been reached the simulation stops. At the end, all information about the entire failure process is saved in a database or a synthetic catalog that can be used for statistical analysis. Further details about the algorithm are given in Sect. 3.

Figure 2 Schematic representation of the considered local load rule. The broken cell with load, σ_F, distributes the largest load fraction, σ_O (Eq. 5), to its four orthogonal neighbor cells. The remaining load, σ_D (Eq. 6), is transferred to its four diagonal neighbor cells. Afterwards, the load of the broken cell drops to zero, σ=0. Asperity cells cannot receive any new load.

Figure 3 (a) Spatial distribution of the random initial loads σ(i,j). R_Asp represents a rectangular fault plane of $N_{\mathrm{x}}=N_{\mathrm{y}}=\mathrm{100}$ cells. The color bar indicates the load and the threshold load of σ_th=1. (b) Spatial distribution of the strength γ(i,j). Two main regions can be distinguished in this figure: (1) the asperity region defined as the inner rectangle and (2) a background area or fault plane. While the asperity contains strength values in the range of 3 to 5, the rest of the fault plane has a strength value of 1.

3.1 Preprocessing: input data and initial conditions

In TREMOL, a fault plane is modeled as a rectangle (Ω), which is divided into N_x×N_y cells. Each cell is defined by its position (i, j), where $i\in [\mathrm{1},\mathrm{\dots },N_{\mathrm{x}}]$ and $j\in [\mathrm{1},\mathrm{\dots },N_{\mathrm{y}}]$. In the fault plane Ω earthquakes can occur with different magnitudes. Additionally, it is possible to assign to each fault plane an asperity region (R_Asp).

To define each fault plane (Ω) and its respective asperity region (R_Asp) it is necessary to assign specific properties to their cells. Particularly, it is necessary to define three properties (or values) for each cell of Ω and R_Asp: a load σ(i,j), a strength value γ(i,j), and a load-transfer value π(i,j).

• The load σ(i,j). At the beginning of each realization, TREMOL randomly assigns a value of the load σ(i,j) to each cell of Ω using a uniform distribution function ($\mathrm{0}<\mathit{\sigma }(i,j)<\mathrm{1}$). This assumption simulates a heterogeneous stress field. Moreover, a load threshold σ_th=1 is necessary to create a limit at which a cell must fail (Moreno et al., 2001). At the end of this step any cell within Ω must have a load value between 0 and 1.

• The strength value γ(i,j). This parameter represents an analogy to the concept of hardness or strength. In our model, the algorithm will find it difficult to break a cell if this cell has a value γ>1 since the strength threshold before failure is set as γ_th=1 (see a detailed explanation in Sect. 3.2). As a result, a strength γ>1 may simulate a hard material that needs to be weakened before it can fail. This process can be regarded as similar to material fatigue or creep failure. The strength value for all cells in R_Asp, namely γ_Asp, is chosen in a discrete interval of ${\mathcal{U}}_{\mathrm{D}}=[{\mathit{\gamma }}_{\mathrm{Ref}}-\mathrm{1},{\mathit{\gamma }}_{\mathrm{Ref}}+\mathrm{1}]$, where 𝒰_D is an integer uniformly distributed and γ_Ref is an assigned reference value.

• The load-transfer value π(i,j). This parameter represents the percentage of load that can be distributed from a ruptured cell to its neighbors. In this study, the load in the ruptured cell is called σ_F(i,j). TREMOL uses a local load sharing (LLS) rule considering the eight nearest neighbors. According to previous studies, such as Monterrubio-Velasco et al. (2017), TREMOL redistributes the majority of the load to the four orthogonal neighbors. The load that is transferred to these orthogonal neighbors is called σ_O, and it is defined according to Eq. (5):

  $$\begin{array}{}\text{(5)}& {\mathit{\sigma }}_{\mathrm{O}}(i,j)={\displaystyle \frac{\mathrm{0.98}{\mathit{\sigma }}_{\mathrm{F}}(i,j){\mathit{\pi }}_{\mathrm{F}}(i,j)}{\mathrm{4}}},\end{array}$$

  where π_F is the load-transfer value of the failed cell. Additionally, a small proportion of the load is transferred to the four diagonal neighbors. The value of this load is called σ_D(i,j), and it is defined according to Eq. (6):

  $$\begin{array}{}\text{(6)}& {\mathit{\sigma }}_{\mathrm{D}}(i,j)={\displaystyle \frac{\mathrm{0.02}{\mathit{\sigma }}_{\mathrm{F}}(i,j){\mathit{\pi }}_{\mathrm{F}}(i,j)}{\mathrm{4}}}.\end{array}$$

  The assumption of Eqs. (5) and (6) is in agreement with 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 and Wysession, 2008). Figure 2 is a schematic representation of the load distribution process from the failed cell, σ_F(i,j) (in red), to its nearest neighbors.

In order to differentiate the parameters of the asperity from the rest of the fault plane Ω, we define π_Asp(i,j) and γ_Asp(i,j) that refer only to the cells within R_Asp. For the rest of the fault plane Ω, we are using the same parameters defined previously: π(i,j) and γ(i,j). Figure 3a shows an example of the randomly distributed initial load throughout the fault plane. Figure 3b displays an example of differences between the strength of the asperity and the rest of the fault plane.

Figure 4 Results of one realization by TREMOL. (a) The spatial distribution of avalanches. Patches of the same color indicate one temporal-consecutive Avalanche cluster (synthetic earthquake). (b) Logarithmic representation of the inter-event rate with time. The red dots represent the inter-event rate when the asperity rupture occurs. The blue dots indicate foreshocks.

3.2 Main computational processes

Once the initial information for the entire domain Ω is defined, the core algorithm of TREMOL will realize a transfer, accumulation, and rupture process. While the cells interact with each other, there are two basic failure processes depending on the load of the cell in comparison with the threshold load (Moreno et al., 2001).

• Normal event. If all cells within the system have a load $\mathit{\sigma }(i,j)<{\mathit{\sigma }}_{\mathrm{th}}$, a normal event is generated, and the cell that will fail is randomly chosen considering the individual failure probability of each cell, F(i,j) (Eq. 4).

• Avalanche event. If one or more cells have a load value $\mathit{\sigma }(i,j)\ge {\mathit{\sigma }}_{\mathrm{th}}$, an avalanche event is generated, and the cell that fails is the one with the greatest σ(i,j) value.

Due to the integrated strength property some extra rules for rupture are necessary. The requirement for failure is $\mathit{\gamma }(i,j)=\mathrm{1}$. On the other hand, if a cell with $\mathit{\gamma }(i,j)>\mathrm{1}$ is chosen, its strength is reduced by one unit. This strength condition enables us to simulate a material weakening process during the load-transfer process. Additionally, this condition offers the possibility to produce large load accumulations locally, which are more likely to generate larger ruptures.

When a cell within R_Asp breaks it becomes inactive until the end of the simulation, which means it cannot receive any further load. The large load concentration within the asperity usually produces a very short time interval (Eq. 3), with the result that there is physically not enough time available to reload the stress on an asperity right after its rupture. On the contrary, a cell outside of the asperity region remains active after its failure but its load drops to zero. The simulation ends when all the cells within the asperity have become inactive.

3.3 Output data and post-processing

After every execution TREMOL outputs a catalog detailing where the position (x,y) of the failed cell, the rupture time (Eq. 12), the avalanche event or normal event identification, the mean load, and many other values are saved for each time step. We cluster avalanche events considering the time and space criterion. We assume $a_{i-\mathrm{1}}=(x_{i-\mathrm{1}},y_{i-\mathrm{1}})$ and $a_i=(x_i,y_i)$ are two consecutive avalanche events generated in chronological order. If their Euclidean distance is Δr_i≤r_th (where $r_{\mathrm{th}}=\sqrt[\mathrm{2}]{\mathrm{2}}$), then a_i and a_i−1 will belong to the same cluster. This clustering algorithm is applied to all generated avalanche events. Lastly, we extract a new catalog that shows the size of each cluster, the position of the first element of each cluster related to the nucleation point, and the time when it was initiated. This database is our simulated seismic catalog. Note that the cluster size is given in nondimensional units. However, we use an equivalence between Ω and an effective area A_eff in order to obtain a physical rupture area. Finally, each cell can represent an area in square kilometers. This step is necessary in order to compute an equivalent magnitude, which is comparable with real earthquake magnitudes. For this purpose, we use three magnitude–area relations. In particular, we use Eqs. (7), (8), and (9) obtained by Rodríguez-Pérez and Ottemöller (2013) for Mexican subduction earthquakes:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\log }_{\mathrm{10}}{\mathit{A}}_{\mathrm{a}}=-\mathrm{4.393}+\mathrm{0.991}{\mathit{M}}_{\mathrm{w}},\text{(8)}& {\displaystyle }{\log }_{\mathrm{10}}{\mathit{A}}_{\mathrm{a}}=-\mathrm{5.518}+\mathrm{1.137}{\mathit{M}}_{\mathrm{w}},\text{(9)}& {\displaystyle }{\log }_{\mathrm{10}}{\mathit{A}}_{\mathrm{a}}=-\mathrm{6.013}+\mathrm{1.146}{\mathit{M}}_{\mathrm{w}},\end{array}$$

where A_a is the asperity area (km²). Equation (7) was obtained from asperities defined by the average displacement criterion (Somerville et al., 1999). Equations (8) and (9) were computed from asperities defined by the maximum displacement criterion for a large asperity and a very large asperity, respectively (Mai et al., 2005).

Furthermore, we define the inter-event rate Δν_k as analogous to the rupture velocity:

$$\begin{array}{}\text{(10)}& \mathrm{\Delta }{\mathit{\nu }}_k={\displaystyle \frac{\mathrm{\Delta }{r}_k}{\mathrm{\Delta }{\mathit{\delta }}_k}},\end{array}$$

where Δr_k is the inter-event Euclidean distance between the k event located at (x_k,y_k) and k−1 event in $(x_{k-\mathrm{1}},y_{k-\mathrm{1}})$.

$$\begin{array}{}\text{(11)}& \mathrm{\Delta }{r}_k=\sqrt{(x_k-x_{i-k}{)}^{\mathrm{2}}+(y_k-y_{i-k}{)}^{\mathrm{2}}}\end{array}$$

The inter-event time Δδ_k is computed following

$$\begin{array}{}\text{(12)}& \mathrm{\Delta }{\mathit{\delta }}_k={\mathit{\delta }}_k-{\mathit{\delta }}_{k-\mathrm{1}},\end{array}$$

where δ_k is given by Eq. (3). Figure 4a shows an example of the final spatial distribution of rupture clusters for a particular example. Each cluster is indicated by the same color and represents a simulated earthquake. Figure 4b shows the related inter-event rate. The inter-event rate largely increases when the asperity rupture occurs.

In the post-processing step we additionally computed the rupture duration of the largest simulated earthquake, D_Aval, using the rupture velocity and the effective fault dimensions obtained from finite-fault models (Table 5).

We used Eq. (13) (Geller, 1976) to compute D_Aval:

$$\begin{array}{}\text{(13)}& D_{\mathrm{Aval}}={\displaystyle \frac{L_{\mathrm{Max}}}{V_{\mathrm{r}}}}+{\displaystyle \frac{\mathrm{16}\sqrt{W_{\mathrm{Max}}\times L_{\mathrm{Max}}}}{\mathrm{7}{\mathit{\pi }}^{\mathrm{3}/\mathrm{2}}\mathit{\beta }}},\end{array}$$

where

$$\begin{array}{}\text{(14)}& \mathit{\beta }\approx {\displaystyle \frac{V_{\mathrm{r}}}{\mathrm{0.72}}}.\end{array}$$

Using these considerations, we can assign a physical unit of time (s) to the largest simulated earthquake, $A_{\mathrm{syn}}=L_{\mathrm{Max}}\times W_{\mathrm{Max}}$.

The flowchart in Fig. 1 and the pseudo-codes 1, 2, and 3 summarize the algorithm of TREMOL. A summary of all required parameters to execute the TREMOL code are shown in Table 1.

Table 1 TREMOL preprocessing: input parameters and their definition.

4 Sensitivity analysis

4.1 Methods: parametric study

We performed a sensitivity analysis of the three asperity parameters (γ_asp, S_a−Asp, and π_asp) in order to identify the best combination that produces the best approximation to real data, such as the maximum rupture area, A_syn, and its related magnitude M_syn. In order to investigate the influence of every single parameter, we statistically determined how the results vary with different parameter configurations.

4.1.1 Percentage of transferred load, π_asp – methods

To explore the influence of π_asp, we analyzed 12 values ($\mathrm{0.67}\le {\mathit{\pi }}_{\mathrm{asp}}\le \mathrm{1.0}$, with increments of 0.3). The minimum π_asp=0.67 assigns the same value to an asperity cell and to a background cell. On the other hand, π_asp=1.0 means that the load in a failed asperity cell is fully transmitted to the neighbors (ideal case with no dissipative effects). Note that π_asp=1 does not represent real physical conditions since dissipative effects are ignored completely. On the other hand, if π_asp=0.67 (case 1) the asperity cells would transfer as much load as the cells in the background. The objective is to generate a load concentration within the asperity that corresponds to the largest magnitude. If the asperity cells transfer as much load as the background cells, no such load concentration can be obtained. As a result, we can expect that the mean A_syn for π_asp=0.67 (case 1) is the lowest value in comparison to all other cases.

The input data of this experiment are summarized in Table 2. We assigned a strength to the asperity (R_Asp) ${\mathit{\gamma }}_{\mathrm{asp}}=\mathrm{4}\pm \mathrm{1}$ and a value of γ_Bkg=1 to the rest of the fault plane. These values are chosen after experimental trials, which have shown that the difference is large enough to simulate a significant strength difference with low computational effort. To define the effective area and the asperity size, we chose the values computed for the earthquake of 20 March 2012, M_w=7.4, in Rodríguez-Pérez and Ottemöller (2013): A_eff=2944.2 km² and S_a=0.26. We defined the size of Ω consisting of N_cell=10 000 cells in total. We carried out 50 simulations per π_asp configuration. In addition, we modified the random seed to have different initial load configurations, σ_(i,j), to ensure that the results over π_asp are independent of the initial load conditions σ_(i,j).

Table 2 Input data in order to carry out cases 1 to 12.

4.1.2 Strength parameter, γ_asp – methods

To perform the parametric study of γ_asp, we configured two asperities embedded in Ω. In this experiment, the total size is $\mathrm{\Omega }=\mathrm{200}\times \mathrm{100}$ cells. Afterwards, we located each asperity in the center of the two sub-domains Ω^′ of 100×100 cells. Figure 5 shows a schematic representation of the domains Ω and Ω^′ used in this experiment.

Figure 5 Schematic configuration for the parametric study of γ_asp and S_a−Asp. The size of the domain Ω is $\mathrm{\Omega }=\mathrm{200}\times \mathrm{100}$ cells. Each asperity is located within the center of the two sub-domains Ω^′ of 100×100 cells. The strength parameter γ_asp and degree of heterogeneity for each asperity can be varied according to the material properties.

The separation between the two asperities remains constant. We chose a value of π_asp=0.90 to produce a large contrast between the asperity and the rest of the fault plane (π=0.67) (Monterrubio-Velasco et al., 2017). In order to analyze the influence of γ_asp (and S_a−Asp), the asperity on the right-hand side (Asp. 2) has varying strength values, while the strength of the left asperity (Asp. 1) remains constant. Finally, the maximum ruptured area and magnitude generated in each Ω^′ are computed.

In order to explore how the system behaves when γ_asp changes, we analyzed six different values of ${\mathit{\gamma }}_{\mathrm{asp}}=[\mathrm{2}\pm \mathrm{1},\mathrm{5}\pm \mathrm{1},\mathrm{7}\pm \mathrm{1},\mathrm{9}\pm \mathrm{1},\mathrm{11}\pm \mathrm{1},\mathrm{14}\pm \mathrm{1}]$ (cases 13 to 18). The input data used in this test are summarized in Table 3. We defined the same asperity size for both: $S_{\mathrm{a}\mathrm{1}}=S_{\mathrm{a}\mathrm{2}}=\mathrm{0.22}$. In Fig. 6, we show an example of the spatial configuration of this analysis. The background strength is considered to be ${\mathit{\gamma }}_{\mathrm{bkg}}=\mathrm{1}=$ constant, and the color bar indicates the γ(i,j) values.

Table 3 Main input data in order to carry out cases 13 to 18.

Figure 6 Example of the strength configuration γ(i,j) for the sensitivity analysis of γ_asp. Two asperities with the same size $S_{\mathrm{a}-\mathrm{Asp}}=\mathrm{0.22}$ are defined and embedded in Ω following the schema in Fig. 5. The conservation parameters are π_asp=0.90 and π_asp=0.67. The color bar indicates different γ(i,j) values. The left asperity (Asp. 1) contains constant properties, while the right asperity (Asp. 2) has variable strength values.

4.1.3 Asperity size, S_a−Asp – methods

The modification of the S_a−Asp parameter was based on the same configuration as described in the previous section. We analyzed six different values of the asperity size S_a (cases 19 to 24). In Fig. 7 we show an example of the asperity configuration in which the left asperity (Asp. 1) has a constant size S_a2, while the size of the right one (Asp. 2) increases. In this experiment, we considered ${\mathit{\gamma }}_{\mathrm{asp}}=\mathrm{5}\pm \mathrm{1}$ and π_asp=0.90. The main data related to these six cases are summarized in Table 4.

Table 4 Main input data in order to carry out cases 19 to 24.

Figure 7 An example configuration of different asperity sizes, S_a−Asp. The color bar indicates the strength γ(i,j) values used during the test.

4.2 Model validation – methods

We evaluated the capability of the model to reproduce the characteristics of 10 Mexican subduction earthquakes (eight shallow thrust subduction events, ST, and two intra-slab subduction events, IN). The input data of the effective area, A_eff, and the asperity ratio size, S_a, are given from waveform slip inversions and seismic source studies ($A_{\mathrm{eff}}=L_{\mathrm{eff}}\times W_{\mathrm{eff}}$ and $S_{\mathrm{a}}={\mathit{A}}_{\mathrm{a}}/A_{\mathrm{eff}}$) shown in the database of Mexican earthquake source parameters by Rodríguez-Pérez et al. (2018). This database includes results from two different methodologies: spectral analysis and finite-fault models. From the latter, the database provides estimations of effective fault dimensions, rupture velocity, source duration, number of asperities, stress, and radiated seismic energy on the asperities and background areas. Slip solutions were obtained with teleseismic data for events with $\mathrm{6.4}<{\mathit{M}}_{\mathrm{w}}<\mathrm{8.2}$.

The number of cells was N_cell=10 000 for a domain Ω of 100×100 cells. We modeled the size of Ω proportionally to the size of L_eff and W_eff for each scenario according to the following equations, Eqs. (15) and (16):

$$\begin{array}{}\text{(15)}& {\displaystyle }{N}_{\mathrm{x}}=\sqrt{{\displaystyle \frac{N_{\mathrm{cell}}{L}_{\mathrm{eff}}}{W_{\mathrm{eff}}}}},\text{(16)}& {\displaystyle }{N}_{\mathrm{y}}={\displaystyle \frac{W_{\mathrm{eff}}}{L_{\mathrm{eff}}}}{N}_x,\end{array}$$

where N_x and N_y are the number of cells in the x axis and y axis, respectively. As an example, Fig. 8 presents the size and aspect ratio of Ev. 3 and Ev. 5 (Table 5).

Figure 8 Example of the domain configuration Ω, considering L_eff and W_eff. (a) Example configuration of event 3 and (b) example configuration of event 5. The required data can be found in Table 5.

In some cases, the number of asperities computed in Rodríguez-Pérez et al. (2018) is greater than 1. However, as a first approximation we simplified the problem by modeling only one asperity per earthquake.

In order to study how the asperity size S_a affects the maximum ruptured area, we randomly modified the size as

$$\begin{array}{}\text{(17)}& S_{\mathrm{a}-\mathrm{Asp}}=S_{\mathrm{a}}+(\mathit{\alpha }\cdot (S_{\mathrm{a}}/\mathrm{2}\left)\right),\end{array}$$

where $\mathrm{0}<\mathit{\alpha }<\mathrm{1}$ is a random value. We introduce this assumption because we want to avoid a preconceived final size. In future trials it may be useful to consider the inner uncertainties of finite-fault models. The asperity aspect ratio follows the same proportion as the effective area, $\frac{{\mathrm{\Omega }}_x}{{\mathrm{\Omega }}_y}=\frac{N_{\mathrm{x}\left(\mathrm{Sa}\right)}}{N_{\mathrm{y}\left(\mathrm{Sa}\right)}}$ (Fig. 8).

We carried out 50 realizations per event (Table 5), changing the size S_a−Asp in each one (Eq. 17).

Rodríguez-Pérez and Zúñiga (2016)Mendoza (1989)Rodríguez-Pérez and Ottemöller (2013)Rodríguez-Pérez and Ottemöller (2013)Rodríguez-Pérez and Ottemöller (2013)Rodríguez-Pérez and Ottemöller (2013)Rodríguez-Pérez and Ottemöller (2013)Rodríguez-Pérez and Ottemöller (2013)

Table 5 The finite-fault source parameters used in this work. W_eff and L_eff are the effective fault dimensions (width and length, respectively, according to Mai and Beroza, 2000). A_real is the asperity area, A_eff is the effective rupture area (W_eff×L_eff). Duration is the rupture duration computed from the slip inversion, N_a is the number of asperities, and V_r is the rupture velocity. Ratio is the aspect ratio of the fault area. The type of the event is labeled ST for shallow thrust and IN for intra-slab events.

4.2.1 Modeling the rupture area and magnitude of 10 subduction earthquakes – methods

In this case the number of cells is N_cell=10 000 (100×100). We carried out 50 executions per event and in each execution we randomly changed the size S_a−Asp following Eqs. (15), (16), and (17). The input data of the 10 modeled earthquakes in Table 5 are summarized in Table 6.

Table 6 Main data used for Ev. 1 to Ev. 10.

4.2.2 Case study (Oaxaca, M_w=7.4, 20/03/2012): using different effective areas A_eff for the same event – methods

As reported in Rodríguez-Pérez et al. (2018) for some events, there are several solutions that allow us to analyze the variability in the estimated source parameters (see parameters of events 7, 7a, and 7b in Table 5). In this study, we applied TREMOL to study how the ruptured area and the assessed magnitude change when we use different input data to model the same earthquake. The data related to these three events are summarized in Table 7.

Table 7 Main data for the case study Ev. 7 test, Ev. 7a test, and Ev. 7b test.

4.2.3 Assessing a future earthquake in the Guerrero seismic gap: rupture area and magnitude – methods

We apply our method for the estimation of possible future earthquakes, in particular to compute the expected magnitude, since TREMOL may offer new insights for future hazard assessments. We carried out a statistical test to assess the size of an earthquake that may occur in the Guerrero seismic gap (GG) region.

As input parameters, we used the area found by Singh and Mortera (1991): $L_{\mathrm{eff}}=\mathrm{230}\mathrm{km}\times W_{\mathrm{eff}}=\mathrm{80}$ km. We defined the asperity size ratio S_a as proposed by Somerville et al. (2002) for regular subduction zone events (SB) based on average slip, S_a=0.25. Singh and Mortera (1991), Astiz et al. (1987), and Astiz and Kanamori (1984) proposed a probable maximum magnitude for this region of ${\mathit{M}}_{\mathrm{w}}\approx \mathrm{8.1}-\mathrm{8.4}$. Therefore, using the effective rupture area (L_eff, W_eff, and S_a), we executed the algorithm as in previous sections. The input data related to this analysis are summarized in Table 8. Likewise, we want to estimated the duration D_aval of the event. To compute this value, we used a mean of the V_r from Table 5.

Table 8 Main data for assessing a future earthquake in the Guerrero seismic gap (GG event).

5 Results

5.1 Results: parametric study

5.1.1 Percentage of transferred load, π_asp

Figure 9 shows the mean (black dots) of the maximum ruptured area A_syn, including the upper and lower limits of the standard deviations (blue squares), after the execution of all 12 cases (Table 2) with 50 realizations. The value of A_syn is related to the largest produced cluster in Ω. There are two dominant tendencies identifiable.

1. If π_asp<0.76, the mean of the maximum ruptured area increases continuously more than 1 order of magnitude from 15 to ≈500 km², i.e., an increase of 3333 %. The standard deviation of A_syn for π_asp=0.7 is ≈35 km² (100 % error).

2. If π_asp≥0.76 (cases 4 to 12), the A_syn values remain essentially constant (≈500 km²). Likewise, the upper and lower limit vary around the same order. The standard deviation for this interval is ≈100 km² (20 % error).

Figure 9 Mean of the maximum rupture area (km²), A_syn, for different values of π_asp depicted as black circles. The minimum and maximum limits of the rupture area are represented by blue squares.

Figure 10 Mean and standard deviation of the maximum magnitude over 50 realizations depending on π_asp.

Figure 11 Mean and standard deviation of the ratio $S_{\mathrm{syn}}=A_{\mathrm{syn}}/A_{\mathrm{eff}}$ over 50 realizations for different values of π_asp. The red line indicates the asperity ratio S_a computed for event 7 (Table 5).

Using the mean of A_syn obtained in each case, we computed the corresponding magnitude. The results are given as the mean and standard deviation of the maximum magnitude in Fig. 10 for all 12 cases (see Table 2). Due to the fact that ruptured area and magnitude are correlated (see Eqs. 7, 8, and 9), the pattern in Fig. 10 is very similar to the one in Fig. 9.

Overall, there are three aspects observable.

1. If π_asp≥0.76 (cases 4 to 12), the mean magnitudes show a steady value (≈7.2).

2. If $\mathrm{0.70}\le {\mathit{\pi }}_{\mathrm{asp}}<\mathrm{0.76}$, a transition with an increasing trend with the largest standard deviation is visible.

3. If π_asp=0.67 (case 1), the mean of the maximum magnitude is the lowest.

In this experiment, the initial value of S_a=0.26 remains constant; i.e., the asperity size does not increase randomly (red line in Fig. 11). After executing all configurations, we computed the ratio of $S_{\mathrm{a}-\mathrm{Asp}}=A_{\mathrm{syn}}/A_{\mathrm{eff}}$, relating to the largest ruptured area. We show the mean and standard deviation of this ratio S_a−Asp in Fig. 11. We observed that the ratio of S_a−Asp is always ≈0.10 lower than S_a.

5.1.2 Strength parameter, γ_asp

For each value of γ_asp (Table 3), we performed 50 executions while changing the initial strength parameter of the asperity γ_asp (Fig. 3b). Likewise, we computed the maximum magnitude obtained for each Ω^′. Figure 12 indicates the mean and standard deviation of the computed maximum magnitude with a dependence on γ_asp. The upper subplot (blue markers) shows the results for the left (constant) asperity (Asp. 1). The lower subplot (red markers) shows the results for the right (variable) asperity (Asp. 2).

Figure 12 Statistical results of γ_asp for a configuration similar to Fig. 6. Markers represent the mean value, while the error bars indicate the standard deviation for all 50 executions considering different initial strength configurations. The red markers correspond to the results of the left asperity (Asp. 1) and blue markers to the results of the right asperity (Asp. 2) (Fig. 6). The strength of the left asperity is kept constant, whereas the strength of the right asperity is variable.

We observe in Fig. 12 that the mean magnitude remains essentially independent for ${\mathit{\gamma }}_{\mathrm{asp}}>\mathrm{5}\pm \mathrm{1}$. Additionally, the error bars slightly decrease, while γ_asp increases. Another observation is that when ${\mathit{\gamma }}_{\mathrm{asp}}=\mathrm{2}\pm \mathrm{1}$ the average of the maximum magnitude is the lowest in both asperities. Moreover, there is a transition zone for $\mathrm{2}\pm \mathrm{1}\le {\mathit{\gamma }}_{\mathrm{asp}}\le \mathrm{5}\pm \mathrm{1}$. We observed that ${\mathit{\gamma }}_{\mathrm{asp}}>\mathrm{5}\pm \mathrm{1}$ has a limited influence on the results of the maximum magnitude. The maximum magnitude of ${\mathit{\gamma }}_{\mathrm{asp}}=\mathrm{14}\pm \mathrm{1}$ is approximately 0.3 magnitudes larger than the one of ${\mathit{\gamma }}_{\mathrm{asp}}=\mathrm{5}\pm \mathrm{1}$.

5.1.3 Asperity size, S_a−Asp

Figure 13 shows the mean magnitude and standard deviation as a function of asperity size. The first asperity with the fixed size indicates a relatively constant magnitude of approximately 7.4. Conversely, the second asperity with variable size produces only a slight increase in magnitude. The magnitude of $S_{\mathrm{a}-\mathrm{Asp}}=\mathrm{0.52}$ is approximately 0.5 magnitudes larger than the one of $S_{\mathrm{a}-\mathrm{Asp}}=\mathrm{0.22}$.

Figure 13 Statistical results of S_a−Asp for a configuration similar to Fig. 7. The markers indicate mean and standard deviation for 50 realizations. Red markers correspond to the results of the right and variable asperity and blue markers to the results of the left and stable asperity.

5.2 Results: model validation

5.2.1 Modeling 10 Mexican subduction zone earthquakes

Based on the observations described in the previous section, we used ${\mathit{\gamma }}_{\mathrm{asp}}=\mathrm{5}\pm \mathrm{1}$ and π_asp=0.90 in order to validate the model. We chose ${\mathit{\gamma }}_{\mathrm{asp}}=\mathrm{5}\pm \mathrm{1}$ because it represents the strength interval of $\mathrm{5}\pm \mathrm{1}\le {\mathit{\gamma }}_{\mathrm{asp}}\le \mathrm{14}\pm \mathrm{1}$ with less computational cost. We chose π_asp=0.90 because it represents the relatively constant magnitude for the parameter range $\mathrm{0.76}\le {\mathit{\pi }}_{\mathrm{asp}}\le \mathrm{0.90}$. In addition, π_asp=0.90 enables us to obtain the best approximation to the ratios of S_a−Asp. Both parameter choices ensure an appropriate reproduction of the asperity rupture area, the maximum magnitude, and least computational payload.

Figure 14 depicts a comparison between the (real) asperity area A_real (Table 5) and the area of the largest simulated earthquake, A_syn. We plot the mean (blue dots), the minimum (green triangles), and the maximum (red triangles) of all 50 realizations for each real earthquake event. Black squares represent the real asperity size. The results in Fig. 14 point out that A_syn is almost identical to A_real from Table 5 for the majority of earthquakes. Only three events show significant differences between synthetic and realistic maximum rupture area. Even in these cases, however, A_real is located within the upper and lower limit of A_syn.

Figure 14 Comparison between the real asperity area A_real and the synthetic values (mean and standard deviation) of the largest ruptured event A_syn. Black squares depict the real asperity area from Table 5, whereas blue circles indicate the mean area of 50 executions. Red and green triangles represent the maximum and minimum A_syn.

Figure 15 shows the statistical results of the synthetic maximum magnitude, M_syn, determined for all 10 events. The real magnitudes from Table 5 are given as red markers. Black circles indicate the mean of M_syn of 50 realizations using Eq. (7), whereas blue and green markers indicate the magnitude following the Eqs. (8) and (9), respectively. The error bars represent the standard deviation. We observed that the statistical parameters computed with TREMOL fit the magnitudes shown in Table 5. However, the computed magnitudes depend on the scale relation employed (Eqs. 7, 8, and 9). Figure 16 includes the mean of the three scale relations. Overall, the mean magnitude ${\stackrel{\mathrm{‾}}{M}}_{\mathrm{syn}}$ and the expected magnitude M_w show similar values. Given that the difference between the mean and the expected value (Table 5) is lower than ΔM_w<0.5 for the 10 events, we can confirm that the results of assessing the magnitude by means of TREMOL using a randomly modified asperity size, S_a−Asp (Eq. 17), are reasonable.

Figure 15 Statistical results of the maximum magnitude for the events from Table 5. Red squares depict the real estimated magnitudes from Table 5, while black circles, blue triangles, and green triangles represent the synthetic mean magnitude of 50 executions following Eqs. (7), (8), and (9), respectively. The error bars are the standard deviation of the scale relations.

Figure 16 Statistical results of the maximum magnitude for the events from Table 5. Red squares represent the magnitudes from Table 5, whereas black circles represent the mean magnitude (${\stackrel{\mathrm{‾}}{M}}_{\mathrm{syn}}$) value of all 50 executions following Eqs. (7), (8), and (9). The error bars stand for the standard deviation of the mean for the three scale relations.

Figure 17 shows the real ratio size S_a from Table 5 (black squares) in comparison to the mean of the largest simulated earthquake, S_syn (blue squares). The standard deviation is represented as error bars. The results indicate that in most of the cases the computed S_syn range fits the expected S_a well. Note that for events 3, 7, and 8 the mean values are lower than the reported S_a, while S_syn is overestimated for events 2, 5, and 9. For events 1, 4, 6, and 10 the estimated value of S_syn coincides with the expected one. However, the error bars encompass the expected values in all cases (Fig. 17). Moreover, if we compare Fig. 17 with Fig. 11 we observe that the employed strategy of randomly increasing asperity size (using Eq. 17) generates rupture areas similar to the ones proposed by Rodríguez-Pérez et al. (2018).

Figure 17 Proportion of simulated ruptured area occupied by the largest avalanche, S_syn, in comparison with the real ratio size S_a from Table 5. The real ratio size S_a from Table 5 is represented by black squares and the mean of the largest simulated earthquake, S_syn, by blue squares. The standard deviation is represented as error bars.

We also computed an equivalent rupture duration, D_Aval, using the equation proposed by Geller (1976) to calculate the rise time (Eqs. 13 and 14). Rodríguez-Pérez and Ottemöller (2013) determined the rupture velocity V_r (Ev. 3–8), which is a useful parameter in order to validate our results. Figure 18 shows the results of this analysis. In red we plot the values V_r calculated by Rodríguez-Pérez and Ottemöller (2013) and in blue the D_Aval based on Eq. (13) with V_r provided by Table 5. The equivalent D_Aval using Eqs. (13) and (14) is printed in black. In cases in which we have the reference values, V_r, computed by Rodríguez-Pérez and Ottemöller (2013), we observe that the reference values are always larger than the modeled D_Aval values.

However, it is worth noting that V_r is the mean rupture time that considers the rupture of the whole effective area (A_eff). For the simulated rupture duration, D_Aval, we only consider the rupture length of the largest rupture cluster A_syn. As a result, smaller values than those proposed in Rodríguez-Pérez and Ottemöller (2013) are expected. Nevertheless, the rupture duration shows a clear dependency on the magnitude.

Figure 18 Equivalent rupture duration D_Aval (s) calculated via the rupture velocity by using the size of the largest rupture cluster. Red squares represent the reference values proposed by Rodríguez-Pérez and Ottemöller (2013), while blue squares and black circles depict the synthetic rupture duration computed by means of V_r based on Eqs. (13) and (14), respectively.

5.2.2 Case study (Oaxaca, M_w=7.4, 20/03/2012)

In the cases in which several effective rupture areas were proposed by different studies (see Table 5), it is possible to assess which set of parameters is better in order to simulate an event by means of TREMOL. We tested TREMOL by using three different combinations of L_eff, W_eff, and S_a according to results for Ev. 7 in Table 5. A comparison of these three combinations is visualized in Fig. 19: panel (a) shows the comparison of the ruptured areas, A_real and A_syn; panel (b) shows the mean and standard deviation of the maximum magnitude, M_syn, in comparison to the reference magnitude; and panel (c) shows the ratio S_syn of the simulated events compared to S_a, the real scenarios. Although the three combinations express similar results, the closest approximation between real and synthetic data is generated based on Rodríguez-Pérez and Ottemöller (2013) (Ev. 7).

Figure 19 A comparison between the data from Table 5 and the results by TREMOL for events 7, 7a, and 7b. (a) Maximum ruptured area, A_syn; (b) mean maximum magnitude, M_syn; (c) ratio of maximum event size, S_syn.

5.2.3 Assessing a future earthquake in the Guerrero seismic gap: rupture area and magnitude

Figure 20 Estimation of the characteristics of a future earthquake in the Guerrero seismic gap. (a) Estimated rupture area, (b) rupture duration, (c) average of the mean magnitude considering three scale relations in Eqs. (7), (8), and (9).

In Fig. 20a, we compare the mean of maximum ruptured area, A_syn, including error bars with the reference area, A_a. The rupture area computed in TREMOL shows a possible range from 4000 to 7000 km². This interval is based on a considered size of S_a=0.25. In the subplot of Fig. 20b, we estimated the duration D_aval of the rupture event. The results in Fig. 20b indicate that the duration D_aval is similar to that of the other events of magnitude M_w≈8. The duration may range from 80 to 110 s, while a rupture duration between 90 and 100 s is most likely. Figure 20c shows the mean of the estimated magnitude using Eqs. (7), (8), and (9). TREMOL outputs a possible range of $\mathrm{8.1}\le {\mathit{M}}_{\mathrm{w}}\le \mathrm{8.5}$, which matches the proposed value by Singh and Mortera (1991), Astiz et al. (1987), and Astiz and Kanamori (1984) of ${\mathit{M}}_{\mathrm{w}}\approx \mathrm{8.1}-\mathrm{8.4}$.

6 Discussion

6.1 Discussion: parametric study

6.1.1 Percentage of transferred load, π_asp

In the results, there were two dominant tendencies visible: (1) π_asp<0.76 and (2) 0.76≤π_asp. If π_asp<0.76 the mean of the maximum ruptured area increased continuously more than 1 order of magnitude from 15 to ≈500 km², i.e., an increase of 3333 %. Therefore, the range of π_asp is both crucial and sensitive. A parameter increase of only 15 % affects the size of the biggest earthquake within the system by 3333 %. Considering the large standard deviation of ≈35 km² (100 % error) a parameter configuration based on π_asp<0.76 would be unsuitable for further simulations due to the unstable properties obtained for that range. The second tendency, however, offers the possibility to determine a stable conservation parameter that can be freely chosen in the range of $\mathrm{0.76}\le {\mathit{\pi }}_{\mathrm{asp}}\le \mathrm{1.0}$. The stable state of maximum rupture area is caused by a self-organized critical avalanche size of A_crit≈500km² based on a grid of 100×100 cells with A_eff=2944.2 km² and S_a=0.26. As soon as A_crit is achieved by the system, the largest avalanche will stop increasing in size, whereas other avalanches within the system will be favored to grow. On the other hand, this means that TREMOL breaks the asperity in patches rather than completely during one unique rupture event (see Figs. 4 and 11). This last condition is reasonable considering that the algorithm of FBM used in TREMOL favors clustering the rupture of cells. Therefore, it is reasonable that some cells remain outside of a unique rupture group because they do not satisfy the failure conditions. As a consequence, we think that it is necessary to define an initial area greater than the expected area of the asperity where the asperity rupture can occur. This result also justifies the proposed Eq. (17), wherein the size of the asperity increases randomly up to 50 % larger than the value proposed by Rodríguez-Pérez et al. (2018). Future studies may be useful to better determine the influence of A_crit.

The parametric study indicates that the largest rupture π_asp is produced as long as it is within the range $\mathrm{0.76}<{\mathit{\pi }}_{\mathrm{asp}}<\mathrm{1.0}$. So even though as π_asp increases, large rupture clusters are generated because a large amount of load is transferred to the neighboring cells, thereby producing critical local load concentrations in the system, the particular lower bound is critical. In our simulation short-range interactions convert to long-range processes through the avalanche mechanism in TREMOL v0.1. The explicit interaction range is given by the parameter π and the local load sharing rule, since this produces a load concentration in neighboring cells, promoting ruptures in a local manner (short range). However, the long range is also captured in a more implicit way.

As mentioned in Sect. 3.2, the algorithm searches for a cell to fail that fulfills two different criteria based on the stress and the strength values of the cells. This property results in long-range interactions since the randomness of the initial stress distribution allows cells at large distances to be activated after a sufficient amount of subsequent steps.

6.1.2 Strength parameter, γ_asp

The parameter γ_asp quantifies the “hardness” of the asperity in comparison to the background material. Its value is given as ${\mathit{\gamma }}_{\mathrm{asp}}={\mathit{\gamma }}_{\mathrm{ref}}\pm \mathrm{1}$. The value γ_ref=2 indicates that the strength in the asperity is twice as big as in the background area. The explored range of γ_ref (2, 5, 7, 9, 11, 14) is based on an experimental trial. Figure 12 presents the mean of the maximum magnitude of an event as a function of γ_asp with two tendencies being visible.

1. There is an unstable transition zone of $\mathrm{2}\pm \mathrm{1}\le {\mathit{\gamma }}_{\mathrm{asp}}\le \mathrm{5}\pm \mathrm{1}$ where the maximum rupture has a strong variation. Therefore, a strength value within this range should be avoided.

2. There is a stable zone of $\mathrm{5}\pm \mathrm{1}<{\mathit{\gamma }}_{\mathrm{asp}}<\mathrm{14}\pm \mathrm{1}$, where γ_asp can be freely chosen. However, due to computational costs it is recommended to use the lowest value of ${\mathit{\gamma }}_{\mathrm{asp}}=\mathrm{5}\pm \mathrm{1}$, since the number of necessary time steps to activate the whole asperity increases strongly with the applied asperity strength (see Algorithm 1).

Moreover, as γ_asp increases the simulation requires a larger number of iterations to break a cell in the asperity, thus implying a larger computational cost. Our selection (γ_ref=5) ensures a “stable” maximum magnitude in the lowest computational time. Figure 21 visualizes the magnitude, the number of steps required to activate the whole asperity, and the computational time in seconds for one execution as a function of γ_ref. In this sense, we considered a value of γ_ref=5 to be adequate from a computational point of view and also to ensure relatively constant values of maximum magnitude.

Figure 21 Magnitude, computational time (s), and number of steps to “activate” the whole asperity for one execution as a function of γ_ref.

6.1.3 Asperity size, S_a−Asp

The results of Fig. 13 indicate that asperity size has a significant influence on the maximum magnitude. We emphasize the importance of these results because they show that the parameter S_a−Asp is critical to control the generated magnitude. At the same time, these results provide the appropriate range of values that TREMOL requires to do a reasonable assessment of the maximum rupture area and magnitude of an earthquake.

6.2 Discussion: model validation

The model validation by means of 10 different subduction earthquakes showed that TREMOL is capable of reproducing rupture area and magnitude appropriately – by means of only few input data – in comparison to the results from inversion studies. The computed rupture duration by TREMOL differs from the reference values. The reason may be that the calculation of the rupture duration is based on the largest (critical) rupture area that is not equal to the available asperity area (see Figs. 11 and 6.1.1). Nevertheless, the rupture duration shows a clear dependency on the magnitude. Since TREMOL only requires few input data, it is a powerful tool to simulate future earthquakes, such as those that might take place in the Guerrero gap region. The determination of the magnitude of an earthquake based on the asperity area depends on the scale relation used. We considered it more appropriate that the relation used be related to the tectonic regime to be modeled. For example, other possible relations to be applied for a subduction earthquake regimen could be Strasser et al. (2010) and Blaser et al. (2010). However, if the user wants to include another empirical relation it is possible to add it in the script TREMOL_singlets/postprocessing/calcuMagniSpaceTime Singlets.jl, such as in the example presented by Wells and Coppersmith (1994). It is worth mentioning that the relation proposed by Wells and Coppersmith (1994) was analyzed as part of our tests. However, we found that the magnitude values are lower than the ones reported by Eqs. (7), (8), and (9). Blaser et al. (2010) discussed similarities and differences with Wells and Coppersmith (1994) in detail. Nevertheless, the rupture area is not model sensitive (Fig. 14), so in order to compare real data and simulations it is more appropriate to use the rupture area.

After validating the capability of the model, constraining the input parameters, and analyzing the results, we consider the conceptual basis of TREMOL to be expandable to model other tectonic regimes. For example, the FBM may be applied to study the rupture process in active fault systems and its effect on aftershock production. Likewise, a three-dimensional version can be developed to simulate mainshock rupture and its aftershocks as reported in Scholz (2018), who tested a first prototype of a 3-D version.

6.3 TREMOL: advantages and disadvantages

The algorithm of TREMOL enables the model to store stress history and to simulate static fatigue due to an included strength parameter γ. The vast majority of asperity parameters have already been examined in previous inversion studies and are usually accessible from online databases.

The range of values found in the sensitivity analysis are not unique for the Mexican examples. In fact, the parameters π, γ, and σ are generic for any simulation of similar types of earthquakes. The only information that needs to be defined beforehand for a unique earthquake is the effective area size and the asperity area, which may come from finite-fault models.

Dynamic deterministic modeling of aftershock series is still a challenge due to both the physical complexity and uncertainties related to the current state of the system. In seismicity process simulations the lack of knowledge of some important features, such as the initial stress distribution or the strength and material heterogeneities, generates a wide spectrum of uncertainties. One way to address this issue may be to consider a simple distribution such as a uniform distribution. We think that the validity of this assumption is given by the comparison of the simulated results with real data. It is possible that other distributions might also give similar results. However, the intention of TREMOL v0.1 is to propose a model that can be used to assign values to the unknown properties mentioned before, including different distributions. Therefore, we encourage users to try other distributions and investigate their effects.

The FBM, on the other hand, produces similar statistical and fractal characteristics as real earthquake series, and its parameters can be regarded as analogous to physical variables. Likewise, the FBM is able to simulate failure through static fatigue, creep failure, or delayed rupture (Pradhan and Chakrabarti, 2003; Moreno et al., 2001).

One disadvantage of TREMOL is that its output is highly dependent on the input, which is based on information from kinematic models and therefore contains inherent uncertainties from inversion studies (see Table 5). TREMOL may be able to compensate for some errors, but how far this possibility can be exploited needs to be investigated in the future. Further steps in the advance of the model have just started, which includes the analysis of a machine-learning approach (Monterrubio-Velasco et al., 2018) that will exploit all the possibilities of this technique.

There are still issues that will likely be addressed in future tests, as outlined below.

1. For our validations, we used earthquakes for which a suitable amount of information is available. How can the technique be applied to other events for which little information is available through, for instance, far-field recordings of seismicity?

2. For our validation study, we used a simplified geometry of the real complex asperity geometries. However, other irregular asperity geometries may be introduced in future works.

3. The FBM is a pure statistical model and therefore gives only hints about underlying physical processes. So far, it does not take into consideration physical effects such as pore fluid pressure, soil amplification, stress relaxation of the upper mantle, reactivation of existing faults, volcanic activity, and many more. One strength of the FBM is that an endless number of information layers can be included into the model that would allow us to include physical properties and topography as well.

4. As it currently stands, TREMOL is not able to simulate complete seismic cycles. Rate-and-state friction models such as by Lapusta et al. (2000) and Lapusta and Rice (2003) have the ability to reload stress. TREMOL is still in an early stage of development and thus lacks a reloading feature.

Additional setbacks of TREMOL are that (1) the number of time steps needs to be adjusted manually for every grid resolution and case scenario, and (2) it is based on a sequential algorithm. In order to save the stress history within every cell of the system, a consecutive algorithm is necessary that changes the state of the system with every time step. This limits the integration of a parallel domain, but a parallel distributed memory is a good approach to solve the problem of large domains. As a result high-performance computing facilities are required when very large grid sizes are used (Monterrubio-Velasco et al., 2018).

Overall, the results of TREMOL are promising. However, the results also point out the need for further modifications of the algorithm and more intensive studies. Likewise, many questions are still left to be answered due to the model's early development stage. In the very near future, however, TREMOL may be a true alternative to classical approaches in seismology. The simple integration of layers of information makes TREMOL a simple model that can be easily modified to simulate the most complex scenarios. At the moment, TREMOL cannot compete with state-of-the-art and widely accepted rate-and-state friction-based models, but it is a totally different, complementary, and promising approach that can provide important insights into earthquake physics and hazard assessment from a completely different perspective. The development of TREMOL and similar models should therefore be strongly encouraged and supported.

7 Conclusions

In this study, we present an FBM-based computer code called the stochasTic Rupture Earthquake MOdeL, TREMOL, in order to investigate the rupture process of seismic asperities. We show that the model is capable of reproducing the main characteristics observed in real scenarios by means of few input parameters. We carried out a parametric study in order to determine the optimal values for the three most important initial input parameters.

• π_asp. as long as the fault plane has a conservation parameter of π_bkg=0.67, the conservation parameter of the asperity must be π_asp≥0.76 to ensure a realistic maximum rupture area.

• γ_asp. the best strength interval for the asperity is $\mathrm{5}\pm \mathrm{1}<{\mathit{\gamma }}_{\mathrm{asp}}<\mathrm{14}\pm \mathrm{1}$. However, due to computational costs it is recommended to use the lowest value of ${\mathit{\gamma }}_{\mathrm{asp}}=\mathrm{5}\pm \mathrm{1}$, since the number of necessary time steps to activate the whole asperity increases strongly with the applied asperity strength (see Algorithm 1).

• S_a−Asp. the generated magnitude can be controlled by parameter S_a−Asp. This parameter is dependent on the earthquake of interest and follows results with data from inversion studies.

We also carried out a validation study employing 10 subduction earthquakes that occurred in Mexico. TREMOL proved that it is able to reproduce those scenarios with an appropriate tolerance.

A big advantage of our algorithm is the low number of free parameters (L_eff, W_eff, and S_a) to obtain an appropriate rupture area and magnitude assessment. Our code TREMOL allows its users to investigate the role of the initial stress configuration and the material properties over the seismic asperity rupture. Both characteristics are key factors for understanding earthquake dynamics. The strengths of our FBM model are the simplicity of implementation, the flexibility, and the capability to model different rupture scenarios (i.e., asperity configurations) with varying mechanical properties within the asperities and/or background area or fault plane. Although we simplified the expected complex asperity geometries, irregular asperity geometries may be introduced in future works. Another advantage is the analysis of earthquake dynamics from the point of view of deformable materials that break under critical stress. The results of TREMOL are promising. However, various assumptions and simplifications require further experiments and modifications of the algorithm to cover various tectonic settings. Likewise, the machine-learning application by Monterrubio-Velasco et al. (2018) needs to be incorporated into the model to determine the optimal parameter ranges for different fault types and tectonic regimes. Although many questions are still left to be answered due to the model's early development stage, TREMOL proved to be a powerful tool that can deliver promising new insights into earthquake triggering processes. Our future work will investigate complex asperity configurations, earthquake doublets, and stress transfer in three-dimensional domains.