A non-stationary earthquake probability assessment with the Mohr – Coulomb failure criterion

From theory to experience, earthquake probability associated with an active fault should be gradually increasing with time since the last event. In this paper, a new nonstationary earthquake assessment motivated/derived from the Mohr–Coulomb failure criterion is introduced. Different from other non-stationary earthquake analyses, the new model can more clearly define and calculate the stress states between two characteristic earthquakes. In addition to the model development and the algorithms, this paper also presents an example calculation to help explain and validate the new model. On the condition of best-estimate model parameters, the example calculation shows a 7.6 % probability for the Meishan fault in central Taiwan to induce a major earthquake in years 2015–2025, and if the earthquake does not occur by 2025, the earthquake probability will increase to 8 % in 2025–2035, which validates the new model that can calculate non-stationary earthquake probability as it should vary with time.


Introduction
Owing to our imperfect understanding and the natural randomness of earthquake, several models have been proposed for estimating earthquake probability in a given period of time.Among them, the Poisson model might be the one that is mostly used in many applications (e.g., Weichert, 1980;Ang and Tang, 2007;Ashtari Jafari, 2010).However, it must be noted that the Poisson model is "memoryless" (Devore, 2008), meaning that the calculation (or the probability) is only a function of the length of time, but irrelevant to when the last earthquake occurred.However, it seems that the recurrence of a characteristic earthquake associated with the same fault should not be memoryless.That is, the earthquake probability should gradually increase with time after a major earthquake has occurred, which is a non-stationary process.
The scope of this study is to develop a new non-stationary model for earthquake probability assessment, mainly motivated/derived from the Mohr-Coulomb failure criterion.Meanwhile, a comprehensive review of other non-stationary earthquake models is also given in the paper (Sect.2), followed by the derivations of our new non-stationary analysis (Sect.3).Moreover, an example calculation is also presented to help demonstrate the new model, and to validate the model's robustness (Sect.4).

An overview of non-stationary earthquake models
In this section, we would like to provide a comprehensive review of non-stationary earthquake analyses and models.Specifically, we referred to the analyses as "statistical models" and "physical models".

Statistical models
Basically, those statistical models developed are more or less a derivative of the Poissonian calculation.For example, Vere-Jones and Ozaki (1982) proposed the use of a time-variant model parameter for the Poissonian calculation, making their model become non-stationary although the calculation is still Poissonian in essence.Similarly, another work suggested the use of adjusted return period (related to current time and original return period) in a Poissonian calculation, in order to modify the Poisson model from being stationary to nonstationary (Wang et al., 2013).Another type of the modification is to use different models than the exponential distribution to model the earthquake's inter-occurrence interval as a random variable.(Note that for an event modeled by a Poisson process, the number of events in a given period of time is a discrete random variable following the Poisson distribution; meanwhile the time when the next event would recur is a continuous variable following the exponential distribution.)For example, the log-normal distribution (Ferráes, 2005), Weibull distribution (Yakovlev et al., 2006), and Gamma distribution (Gómez and Pacheco, 2004) all have been suggested for the replacement of the exponential distribution, with them all featuring a non-stationary calculation after such modifications.
Based on given earthquake data, it must be noted that the statistical models are all empirical in a sense.In other words, the models do not consider earthquake mechanisms, such as tectonic stress accumulation and the resistance of fault planes.

Physical models
In consideration of earthquake mechanics, several nonstationary earthquake analyses have also been proposed from a different perspective.It must be noted that the models are not entirely a "product" of physics, but are somehow based on the concepts/theories working together with empirical relationships.Specifically, we would like to introduce three of them in the following that are more related to our nonstationary earthquake model.
, where X(t) is the stress at time t, λ is the longterm stress increment rate, and σ is the magnitude of a Brownian motion W (t).
The first one we would like to introduce is the timepredictable model (Shimazaki and Nakata, 1980).Figure 1 is a schematic diagram illustrating the basics of the model.Essentially, the model relies on a best-estimate relationship between coseismic fault slip (or displacement) and time.For instance, given the last event with fault slips as shown in Points A and B of Fig. 1, the next event should recur at the time of Points C and D. As a result, the concept of the method is that the recent event with a smaller fault slip should accompany a smaller stress drop, and with a constant stress increment with time, it should lead to a shorter time for the stress to re-reach a stress level (or failure stress state) that could induce earthquakes.In other words, if there is a larger slip rate from the recent earthquake, a longer time should be needed for the next event to recur.
The next of this "physical-model" group is the Brownian model (Ellsworth et al., 1999;Matthews et al., 2002).In contrast to the previous model, the Brownian model is not based on a constant stress increment, while considering the stress increments between two consecutive events a stochastic process.As shown in Fig. 2, the model considers the stress-time series to be a combination of a long-term stress increment and a Brownian motion that simulates transient stress randomness, and with such a function we can estimate the date of the next earthquake by examining if the stress reaches the failure state within a given time.However, it must be noted that the so-called failure state of the model is a concept or a hypothesis, which cannot be formulated or calculated with a pressure or force even though the strength properties of fault planes are given.
The third model we would like to introduce is the negative binomial model (Tejedor et al., 2015).As with the previous analyses, the model is also based on two imaginary stress states.As shown in Fig. 3, the essence of the model is that the stress change in unit time could be modeled by two scenarios: "stress increases" and "stress does not increase".As a result, there are many possible "stress routes" (as shown in Fig. 3) between two consecutive events, and the time needed for each route and its probability can be calculated with a given earthquake return period.Finally, the inter-occurrence interval calculated can be proved as a negative binomial distribution for such a non-stationary probability assessment.
To sum up, the three physical models are all facilitated with stress states, although they cannot be calculated with a force/pressure from earthquake theories/hypotheses.Although we do share this perspective for our model development, the biggest difference is that our model defines and calculates the two stress states more clearly, on the basis of the Mohr-Coulomb failure criterion that are well established, and used in rock mechanics, structural geology, etc.
3 The new non-stationary earthquake probability assessment

Overviews of Mohr-Coulomb failure criterion and elastic rebound theory
The Mohr-Coulomb failure criterion is a model describing the response of materials subject to external stresses (Pariseau, 2007), and it is commonly applied to rock mechanics as well as other applications.Figure 4 is a schematic diagram illustrating the essentials of the model.Basically, as the Mohr circle is below the failure envelope, a shear failure is not expected.By contrast, as long as the Mohr circle is in contact with the failure envelope, a shear failure could occur in rock.
In this paragraph, we would like to elaborate on the Mohr-Coulomb failure criterion.Using Fig. 4 as an example, Mohr circle A presents the major principal stress σ 1 (the bigger one) from the vertical direction, and the minor principal stress σ 3 (the smaller one) from the horizontal direction.But with the external compressional stress increasing in the horizontal direction, later on the horizontal stress will be larger than the vertical stress, making the horizontal and vertical stresses become σ 1 and σ 3 instead at the time.Specifically, this is the case for the thrust-fault earthquake owing to tectonic compression in the horizontal direction.More explanation will be given in the following.
On the other hand, it is generally accepted that tectonic activities are the main reason causing rock failures under the ground, resulting in an earthquake with the release of accumulated strain energy.Afterwards, the energy reaccumulates and re-releases until the next earthquake, becoming the foundation of the so-called elastic rebound theory (Keller, 1996) that was proposed by Reid in the early twentieth century (Reid, 1910).More importantly, those nonstationary earthquake models mentioned earlier were somehow motivated by the elastic rebound theory.between the two stress states, many "stress routes" can be present, and the probability of each route can be estimated with the model, then the probability distribution for the interval between two consecutive events can be developed.

The basics and algorithms of the new non-stationary model
The two earthquake theories above were the main motivation and the foundation of the new non-stationary model: (1) based on the Mohr-Coulomb failure criterion, the rock subject to the stress state as Mohr circle C (see Fig. 4) should fail and cause an earthquake, at which we refer to it as failure state; (2) from the elastic rebound theory, the stress state right after a characteristic earthquake should be restored to Mohr circle A, which is called the initial state at time t 0 , in the process towards the next (characteristic) earthquake.As a result, the problem of evaluating the earthquake probability within a given time t * after the last event (or after t 0 ) is becoming a problem as follows: what is the chance of the major principal stress at time t * (denoted as σ 1_t * ) being greater than the stress at the failure state (denoted as σ 1_failure )?In that sense, the question can be formulated as follows: Clearly, the problem is governed by two variables, σ 1_t * and σ 1_failure , and their relationships with other parameters will be detailed later.Note that the notations used in the following derivations are summarized at the end of the paper (Table A1).
-The major principal stress at failure state, σ 1_failure : based on the Mohr-Coulomb failure criterion, the major principal stress σ 1 at failure state (Point C in Fig. 4) can be expressed as a function of the minor principal stress σ 3 , and two strength parameters of the shearing plane, i.e., cohesion c and friction angle ϕ (Pariseau, 2007): -The minor principal stress at failure state, σ 3_failure : the minor principal stress at failure is attributed to the overburden earth pressure (in vertical direction) above the seismogenic depth d, which can be estimated with the following formula based on rock mechanics: where γ is rock unit weight, and σ v denotes the stress in the vertical direction.It must be noted that the case shown in Fig. 4 and Eq. ( 3) is for a thrust-fault earthquake specifically; that is, at the initial state at t 0 , the vertical stress is larger than the horizontal stress because the rock's lateral earth coefficient is smaller than 1.0 (Pariseau, 2007).However, near the failure state, the vertical stress will be smaller than the horizontal stress because of the increasing external force in the horizontal direction from tectonic compression.
On the other hand, the Mohr circles of the initial and failure states for the case of strike-slip faults are shown in Fig. 5, indicating σ 3_failure is equal to γ × d × K for this case, where K is the rock's lateral earth coefficient.
As to the case of normal-fault earthquakes subject to tension, the Mohr circles of the initial and failure states are shown in Fig. 6.
With the three cases shown in Figs.4-6, it should be understood that the directions of major and minor principal stresses can vary with time; that is, the major (bigger) and minor (smaller) principal stresses can be from vertical and horizontal directions at the initial state, while they change to horizontal and vertical directions near the state of failure, which is exactly the case of thrustfault earthquakes as shown in Fig. 4. By contrast, for the normal-fault case subject to horizontal tension as shown in Fig. 6, the stresses in the vertical and horizontal di- rections will always be the major and minor principal stresses, either at the initial state or at failure state, given the stress in the horizontal direction is always smaller than that in the vertical direction.
-The major principal stress at time t * , or σ 1_t * : for the thrust-fault earthquake due to tectonic compression in the horizontal direction, the key task of the new analysis is to estimate the major principal stress at time t * after the last event.For this case, the major principal stress at time t * can be formulated as follows: where ASI represents the annual stress increment.Note that both σ 1_t * and σ 3_initial of this thrust-fault case are the forces in the horizontal direction.As explained previously, this is because at the initial state the horizontal force is smaller than the vertical force due to the rock's lateral earth coefficient less than 1.0.But with the increment of external stress in the horizontal direction, the horizontal stress at time t * will exceed the force in the vertical direction, making the horizontal force become the major principal stress at the time, thus denoted as Moreover, given the rock's lateral earth coefficient

The return period t˜and its relationship with σ 1_t *
In addition to γ , d, K, etc., the return period t of characteristic earthquakes is another input of the non-stationary analysis.Moreover, the mean value and standard deviation of σ 1_t * can be expressed as a function of t, and is used for developing its probability density function for the non-stationary probability assessment.
From the meaning of return period, it is understood that the event will recur when return period t is due.As a result, the major principal stress at return period t (denoted as σ 1_ t ) should be equal to σ 1_failure .Note that for the thrust-fault case as explained previously, σ 1_ t , σ 3_initial , and σ 1_failure in the following equation all denote the forces in the horizontal direction: Therefore, the mean value of ASI (denoted as µ ASI ) can be derived as follows: where E[ ] denotes the mean value of a variable in probability and statistics.
On the other hand, as the variability of annual stress increment is equal to n in terms of coefficient of variation or COV (i.e., standard deviation/mean value), its standard deviation (denoted as s ASI ) can be derived as follows with its mean value from Eq. ( 6): With the mean (Eq.6) of ASI, we can continue deriving the mean value of the major principal stress at time t * : Similarly, the standard deviation of the major principal stress at time t* (denoted as s σ 1_t * ) can be derived as follows with s ASI in Eq. ( 7): V [ ] denotes variance in probability and statistics, which is the square of standard deviation.
In order to establish the probability density function of σ 1_t * , the information about what type of probability distribution the variable is following is as essential as its mean value and standard deviation.But to the best of our knowledge, there is no study providing tangible evidence to the distribution model of σ 1_t * , so we suggest the use of the normal distribution for this random variable, as the normal distribution is usually recommended for a probability analysis when the variables' distribution is unknown (Abramson et al., 2002).

Summary
Figure 7 is a schematic diagram illustrating the essentials of the non-stationary assessments.The fundamental idea of the analysis is to estimate the probability distribution of the major principal stress at time t * after the last event (or after t 0 ), and compare it to the stress that could cause earthquakes.Simply speaking, the governing equation of the analysis can be formulated as Eq. ( 1), and the mean value and standard deviation of the major principal stress at time t * can be derived as Eqs.( 8) and ( 9), respectively.
Most importantly, different from other non-stationary models using the initial and failure states conceptually, the new model can formulate and calculate the forces of the two stress states with the well-established Mohr-Coulomb failure criterion.
To sum up, the new non-stationary model is governed by a total of six parameters as follows: return period ( t), two faultplane strength parameters (c and ϕ), rock unit weight (γ ), earthquake seismogenic depth (d), and COV of annual stress increment (n).

Presumption and limitations
The elastic rebound theory is a plausible explanation for earthquakes, but specifically speaking, it is more of a theory for main shocks.As a result, the new non-stationary analysis of the study motivated by such a theory is more applicable to main shocks, a situation similar to other non-stationary models which are also more applicable to main shocks rather than dependent shocks (Shimazaki and Nakata, 1980;Ellsworth et al., 1999;Matthews et al., 2002;Tejedor et al., 2015).
On the other hand, like any other temporal earthquake probability analyses, our model cannot predict the magnitude of the recurring event, either.In other words, the (earthquake) temporal model usually governed by a given return period to predict earthquake probability in a given time cannot predict earthquake magnitudes or energy release.Again, such a framework is similar to other earthquake temporal analyses only focusing on the earthquake probability in a given period of time, but not focusing on the probability distribution of earthquake magnitude or energy release when the event recurs.

The purpose of the model demonstration
In this section, we would like to present an example calculation using the new non-stationary model, even though the six model parameters (such as COV of annual stress increment) cannot be obtained without a much more extensive investigation.This situation is similar to other non-stationary models.For example, although the Brownian model introduced earlier is acceptable with its concept/methodology, it is difficult to use it in some case studies because the model parameters, such as the long-term stress increment rate and the magni- tude of a Brownian motion, are hard to determine in a totally objective manner.This might be the reason that no real case study using the Brownian model has been reported yet so far, although the model is generally acceptable in terms of methodology.
Nevertheless, we still consider it is beneficial to provide a model demonstration to accompany the methodology for a better presentation/understanding of the new model proposed.Moreover, given the earthquake probabilities did vary with time from the example calculation, the demonstration helps validate the non-stationary analysis, the main purpose of providing the example calculation in the study.

The Meishan fault in central Taiwan
Specifically, we used the Meishan fault in central Taiwan as the example calculation.One of the reasons is that the return period of the characteristic earthquake induced by the active fault has been proposed as 162 years (Lin et al., 2008), and it has been used in a couple of applications (Wang et al., 2012(Wang et al., , 2013)).For example, Wang et al. (2013) used the data in their seismic hazard assessment with a focus on the characteristic Meishan earthquake, evaluating the annual rate of PGA > 0.23 g in nearly cities of the active fault in central Taiwan.
Figure 8 shows the location of the Meishan fault in central Taiwan.In 1906, this strike-slip fault caused the socalled Meishan earthquake at M w = 6.4.At that time, around 1200 people were killed because of the earthquake.
However, to the best of our knowledge, there is no consensus before the 1906 earthquake on how many characteristic earthquakes, and when, had been induced by the active fault.As a result, it is impossible to estimate the standard deviation of the earthquake return period from conventional statistical approaches with only one "confirmed" data point available.
Similarly, without extensive investigations and studies in the future, we believe no research team can provide a totally transparent estimate on COV of annual stress increment near the Meishan fault in central Taiwan, which is another key model parameter for the new non-stationary analysis.As a result, the set of the best-estimate model parameters used in the following example calculation is inevitably related to our judgment.Nevertheless, we would like to emphasize that the purpose of providing the example calculation is to help understand and validate the new non-stationary analysis.Admittedly, the model parameters used will not be perfectly agreeable, but at the same time we believe the problem cannot be addressed by any research team at this point without extensive investigations and studies on annual stress increment variability (n), the rock's lateral earth coefficient (K), the fault plane's cohesion (c) and friction angle (ϕ), and the earthquake return period ( t) and its range or standard deviation.

The best-estimate data for this model demonstration
Table 1 summarizes our best-estimate model parameters for the example calculation.First, we used a typical range (see Table 1) from rock mechanics as our best-estimate friction angle and cohesion (Pariseau, 2007), given the strength parameters of the fault plane few miles under the ground are not clear.Similarly, a probable range of 0.2-0.5 was used as our best estimate for the rock's lateral earth coefficient (K), given no site-specific studies and data have been reported.As for the earthquake focal seismogenic depth, we considered the depth should be close to 6 km, like the last Meishan earthquake (Ng et al., 2009).However, in order to account .The earthquake probability in three 10-year periods for the example calculation, with the earthquake return period of 162 ± 50 years and other input data in Table 1 Figure 9.The earthquake probability in three 10-year periods for the example calculation, with the earthquake return period of 162 ± 50 years and other input data in Table 1.
for this uncertainty, we used a best-estimate range as 4-8 km in the analysis.
A similar situation was encountered in the determination of the best-estimate return period as mentioned previously.As a result, on the basis of 162 years as the mean value from the literature (e.g., Wang et al., 2013), 12 best-estimate ranges from 162 ± 10 years and 162 ± 120 years were used in the model demonstration.Another purpose of carrying out the 12 analyses is to examine how sensitively this model parameter could affect the non-stationary analysis.
As for COV of annual stress increment (n), to the best of our knowledge, there is no research so far that can really answer the question.As a result, the range of 0.25-1 was used as our best estimate characterizing the variability of annual stress increment near the Meishan fault in central Taiwan.More discussion about the input-data characterizations is given in Sect.5.1.

Monte Carlo simulation
Because our input data were characterized by a range rather than a single value, it is difficult to solve the non-stationary probability (i.e., Eqs. 1, 8, and 9) with analytical approaches.Therefore, we used Monte Carlo simulation (MCS) to solve the problem as it has been used in many different applications.For more details about Monte Carlo simulation, readers can refer to the textbooks of Ang and Tang (2007), Abramson et al. (2002), among many others.

The result
With our best-estimate model parameters summarized in Table 1, Fig. 9 shows the average probabilities for three 10year periods from Monte Carlo Simulation with a sample size of 5000.For example, the non-stationary model shows a 7.6 % probability for the Meishan fault to induce a major earthquake in years 2015-2025, under the return period of .The earthquake probability in three 10-year periods for the example calculation, with the earthquake return period of 162 ± 100 years and other input data in Table 1 Figure 10.The earthquake probability in three 10-year periods for the example calculation, with the earthquake return period of 162 ± 100 years and other input data in Table 1.
162 ± 50 years, along with other model parameters given in Table 1.Then, if the event does not occur by 2025, the earthquake probability in 2025-2035 will increase to 8 %; similarly, if the event does not occur by 2035, the probability will further increase to 8.4 %.
Figure 10 shows the result for another scenario under a return period of 162 ± 100 years.Similar to the previous calculation shown in Fig. 9, the earthquake probabilities in the three 10-year periods with different starting dates increase with time.Therefore, from the two example calculations, the new analysis indeed calculated a non-stationary probability, which provided some validation to the new non-stationery model derived from the Mohr-Coulomb failure criterion.
But interestingly, although the two calculations (Figs. 9 and 10) show the average probabilities (the average value of 5000 randomizations from Monte Carlo simulation) are very close to each other, they also demonstrate that the average probabilities are smaller with a larger range in earthquake return period (i.e., 162 ± 100 years) as input.Although such a result seems counterintuitive, the cause and explanation are given in the following section with the support of more calculations.
Figure 11a shows the results of the 12 calculations with return periods from 162 ± 10 to 162 ± 120 years.We can see that the output probability range indeed increases with the input range of earthquake return periods.However, the average probability (output) is almost identical, or the average probability is insensitive to the input data from 162 ± 10 to 162 ± 120 years.
Our speculation of the results shown in Fig. 11a is as follows: although the input-data range is symmetrical (e.g., 162 ± 50 years), the range of the output can be highly asymmetrical owing to the non-linear analysis.To prove this, we processed the 5000 data from our MCS calculation as Fig. 11b, and we can see that the output distribution is not symmetrical as input, causing the average value and the middle value (i.e., (maximum + minimum)/2) to be quite differ- ent.Therefore, the cause of a slightly larger average probability (output) with a smaller range of earthquake return period (input) is due to the non-linear relationship of the calculation, making the output range become asymmetrical although the input distribution is symmetrical.An additional assessment of the new model is to compare the example calculation with the commonly-used Poisson model (although such a stationary model is not realistic).Surprisingly, as seen in Figs. 9 and 10, the two models are in good agreement, calculating the earthquake probabilities in 6 % (Poisson) and 8 % (new model) for the next Meishan earthquake to occur within 2015-2025 in central Taiwan.This additional information may somehow add extra support to the model's robustness.

Input-data characterizations
As many analyses, input-data characterizations are equally or even more challenging than developing the model itself. .Schematic diagram illustrating the stationary process after combining many non-stationary processes; taking T = t 0 and T = t 1 for example, the sum of that many non-stationary probabilities will be close to each other, although the probability is very low for fault D at T = t 0 , and it is very low for fault A at T = t 1 .
As to this new model, we hope to see more studies focusing on the model-parameter calibrations, such as COV of annual stress increment (n), the rock's lateral earth coefficient (K) and the fault plane's strength parameters (c and ϕ) close to the focal point of the earthquake.However, this is a task beyond the scope of this study, and it cannot be completed without extensive investigations and studies.On the other hand, one could argue why a geologically well-investigated fault (e.g., the Chelungpu fault) in Taiwan is not chosen as the model application to reduce uncertainty, and here is our response: the engineering parameters of the model (i.e., lateral pressure coefficient K, strength parameters c and ϕ, etc.) are not clear either for those so-called well-investigated faults.As a result, no matter which fault is selected as the model application, engineering judgment must be involved in the determinations of those model parameters, more or less creating the same level of uncertainty when it comes to site characterizations on the model parameters K, n, c, ϕ, etc.

Should earthquake occurrence follow a stationary or non-stationary process?
Although characteristic earthquakes related to a given active fault should be non-stationary, in the 1970s a study has provided statistical evidence to the opposite: earthquake is stationary (Gardner and Knopoff, 1974).However, it must be noted that the study was not focusing on characteristic earthquakes, but it was based on the regional seismicity in California.
Figure 12 is a schematic diagram that helps explain the difference between the two problems.For each fault, the recurring earthquake should be a non-stationary process, and the non-stationary earthquake probability would be reset at the last event and gradually increase with time.By contrast, the seismicity in a region would become stationary with so many non-stationary processes present.For example, at T = t 0 (see Fig. 12), the sum of that many stationary probabilities should be close to that at T = t 1 (or at any moment), although the earthquake probability induced by fault D should be very low at T = t 0 , while others are higher.
The relationship can be simply explained with the "patronand-bank" analogy.Each patron (analogy to each fault) going to the bank is obviously a non-stationary process, with the probability increasing with time since the very last visit.But for the banks (analogy to the seismicity), it is a stationary process irrelevant to time, as so many patrons or so many non-stationary processes are being dealt with at one time.

Summary and conclusion
Given that the characteristic earthquake associated with an active fault should be a non-stationary process, this paper introduces a new non-stationary analysis to evaluate earthquake probability within a given period of time.Different from previous models, the new analysis more clearly defines and calculates two earthquake stress states, on the basis of the well-established Mohr-Coulomb failure criterion.
In addition, this paper also presents a model demonstration accompanying the new non-stationary model.On the condition of the best-estimate model parameters, the active fault has a 7.6 % probability of inducing the next Meishan earthquake in 2015-2025, and if the earthquake does not occur by 2025, the non-stationary probability will increase to 8 % in 2025-2035, which validates the new non-stationary model that can calculate non-stationary earthquake probability varying with time.Major principal stress at initial state σ 3_initial Minor principal stress at initial state σ 3_t * Minor principal stress at t *

Figure 1 .
Figure 1.Schematic diagram for the time-predictable model: (a) best-estimate relationship between cumulative coseismic slips and time, and (b) the earthquake-time prediction facilitated with a failure state and a constant stress increment.

Figure 2 .
Figure2.Schematic diagram showing the essentials of the Brownian model; within the two imaginary stress states, the model considers the stress-time series should be random and can be modeled by a long-term stress increment and a Brownian motion as X(t) = λt + σ W (t), where X(t) is the stress at time t, λ is the longterm stress increment rate, and σ is the magnitude of a Brownian motion W (t).

Figure 3 .
Figure3.Schematic diagram illustrating the negative binomial model; between the two stress states, many "stress routes" can be present, and the probability of each route can be estimated with the model, then the probability distribution for the interval between two consecutive events can be developed.

Figure 4 .Figure 5 .
Figure 4. Schematic diagram illustrating the Mohr-Coulomb failure criterion; circle A represents the initial state after a thrust-fault earthquake or at t 0 , circle B denotes stress states at t * after t 0 , and circle C is the stress state corresponding to the failure state that causes rock failure and earthquakes.

Figure 6 .
Figure 6.The Mohr circles for the non-stationary earthquake analysis for normal-fault earthquakes.

Fig. 7 .Figure 7 .
Fig. 7.The essentials of the new non-stationary model: Developing the probability distribution of the major principal stress at time t* (i.e., 1_ * t  Figure 7.The essentials of the new non-stationary model: developing the probability distribution of the major principal stress at time t * (i.e., σ 1_t * ) after the last event or after t 0 .

Fig. 8 .Figure 8 .
Fig. 8.The location of the Meishan fault in central Taiwan Figure 8.The location of the Meishan fault in central Taiwan.
Fig.9.The earthquake probability in three 10-year periods for the example calculation, with the earthquake return period of 162 ± 50 years and other input data in Table1 Fig.10.The earthquake probability in three 10-year periods for the example calculation, with the earthquake return period of 162 ± 100 years and other input data in Table1

Fig. 11 .Figure 11 .
Fig. 11.The result of the 12 calculations using return periods from 162 ± 10 to 162 ± 120 years: a) the average probability and the range from the 12 MCS calculations, and b) the asymmetrical distribution in the output data because of the non-linear relationship between input and output in this analysis

Fig. 12 .
Fig. 12. Schematic diagram illustrating the stationary process after combining many non-stationary processes; taking T = t0 and T = t1 for example, the sum of that many non-stationary probabilities will be close to each other, although the probability is very low for Fault D at T = t0, and it is very low for Fault A at T = t1 Figure 12.Schematic diagram illustrating the stationary process after combining many non-stationary processes; taking T = t 0 and T = t 1 for example, the sum of that many non-stationary probabilities will be close to each other, although the probability is very low for fault D at T = t 0 , and it is very low for fault A at T = t 1 .

Table 1 .
Summary of the model parameters used in the analyses.
* K is the rock's lateral earth coefficient; * * n is COV of annual stress increment.

Table A1 .
Notations.t0Thetime when the last event occured t *The time interval after the last event σ 1_t *Major principal stress at time t * σ 1_failure Major principal stress at failure state σ 3_failure Minor principal stress at failure state * Standard deviation of σ 1_t * σ 1_initial