Seismic Forecasting Using a Brownian Passage Time Distribution

Abstract

Seismic forecasting using a Brownian Passage Time distribution is presented in this chapter. Seismic forecasting is concerned with the probabilistic assessment of general seismic hazard, including the frequency and magnitude of earthquakes in a given area over a given period of time. Seismic forecasting generally look for trends that lead to an earthquake. The estimation of the time that a strong earthquake will occur requires the determination of the distribution that the earthquake recurrence time follows. Brownian Passage Time distribution describes reliably the physical processes related with earthquakes’ occurrence. The model assumes that the evolution of the stress loading between two earthquakes depends on the constant loading rate, and a random component, which follows the Brownian Relaxation Oscillator. Its hazard function is in good agreement with the temporal evolution of earthquake occurrence as the hazard rate is very low after an earthquake, and then increases as time passes, and takes a maximum value at the mean recurrence time and since then, it decreases asymptotically exhibiting a pure quasi-periodic temporal occurrence.

1. Introduction to seismic forecasting

Seismic forecasting is concerned with the probabilistic estimation of the frequency and magnitude of seismic events in a given area over a given period of time using advanced statistical and scientific methods. This can be distinguished from seismic prediction, which is the specification of the time, location, and magnitude of a future seismic event with sufficient precision that a warning can be issued. The two can be further distinguished from seismic warning systems, which is the detection of a possible seismic event in real time to regions that might be affected. Seismic warning systems focus on a very short time outlook ranging from few minutes to few days, seismic prediction looks at specific future time with acceptable time ranging from few minutes to few hours while seismic forecasting makes use of probability estimations with time scales ranging from few days up to several decades.

In this chapter on seismic forecasting, literature body on seismic formation process will be reviewed to develop a theoretical framework that validates and upholds ideas that will be further used to develop a stochastic model that can be used to forecast future seismic events. The elastic rebound theory will be validated, and an analysis of recurrence time models will be done so as to select the best model that can be used in seismic forecasting. Given a homogeneous, consistent, and complete past seismic data of a region, the unbiased maximum-likelihood estimates of model parameters can be estimated and used as input parameter to seismic forecasting.

1.1 History of seismic forecasting

Since the ancient times attempts to predict seismic events were made, with people associating such events with the spiritual world. In some societies, such events were considered a sign of bad luck or punishment for disobedience from the supernatural beings. The scientific revolution was a game changer in this mistrial subject with scientist being optimistic that a practical method of seismic prediction would soon be found. By the end of the nineteenth century continued failure leads to many people questioning whether it was even possible to predict a seismic event. Scientific evidence of few predictions of large seismic events has not occurred, and a few claims of success prediction remain controversial. As a result, emphasis has been shifted from seismic prediction to seismic forecasting. Due to the high level of destruction and loss of lives, after larger seismic events, a lot of scientific and national government resources have been pooled and allocated into seismic forecasting rather than prediction as it has proved to be useful in seismic risk mitigation in areas such as establishment of building codes, insurance rate structures, awareness, preparedness programs, and public policy related to seismic events. Statistical methods used for seismic forecasting look for trends or patterns that lead to a seismic event. The trends involve many complex variables, and the advanced statistical techniques are needed to understand them. These approaches tend to have relatively long time periods, making them useful for seismic forecasting.

1.2 Elastic-rebound model

Previously, it was thought that ruptures of the surface were the result of strong ground shaking rather than the converse suggested by Harry Fielding Reid [1], the first scientist to explain the seismic formation process after the great 1906 San Francisco earthquake. The theory postulated that steady tectonic force causes strain to accumulate slowly in a rock and eventually become large up to a threshold constant value called the elastic limit. The elastic limit is the maximum strain that a rock can withstand without breaking. When the limit is exceeded, an earthquake will occur. At that time, a sudden movement occurs along the fault line, releasing the accumulated energy, and the rocks snap back to their original undeformed shape. After an event another cycle starts. When the accumulated strain is great enough to overcome the strength of the rocks, an earthquake occurs again. The duration of an “earthquake cycle” is the ratio of event-strain release to tectonic strain rate. Rocks in a fault plane are subjected to tectonic force caused by tectonic plate movement. Fault plane subjected to plate tectonic moves at the rate of a few centimeters per year, over a time period of decades. The stored energy is released during the rupture partly as heat, partly in damaging the rock, and partly as seismic waves.

The elastic-rebound model plays a central acceptable role in our understanding of earthquake mechanics. This approach, which has some observational grounding, has been the basis for long-term forecasting models. A number of statistical models have been proposed for seismic forecasting. Discrete probability models are proposed to forecast a number of events in a given time interval and continuous probability models have been proposed to forecast the time until the next seismic event.

It is not a good idea to model seismic recurrence using a normal distribution as it gives positive probability to negative intervals. Better models are the Weibull gamma and log-normal distributions, and can be used as alternatives. In this chapter, we are going to consider the Weibull distribution as it is practically convenience as proven by its wide application in statistical quality control. Nishenko and Buland [2] identified several theoretical distributions after normalizing 15 characteristic seismic sequences data by looking at the shape of a generic distribution for recurrence intervals.

Nishenko and Buland [2], agreed with Hagiwara’s [3] model preferences of the exponential, Weibull, gamma, and lognormal distributions by analyzing the reliability of the distributions. They further pointed that the log-normal provides the best fit to the distribution of normalized intervals, and there was no difference in the estimated parameters of the log-normal normal distribution in seismic data from different regions with different time scales. They all postulated that log-normal distribution was a good model for seismic risk assessment.

The elastic rebound model provides a framework seismic recurrence modeling. Extending advanced probabilistic modeling to this theory provides a basic model for seismic forecasting. The idea is the inclusion of the seismic random perturbations in the elastic rebound model.

The same approach was used by Kagan and Knopoff [4] in modeling time-dependent model in statistical seismic analysis. Empirical analysis does not use a Poisson process in many cases. Matthews et al. [5] disagreed with Kagan and Jackson [6] on the idea that seismic recurrence time intervals are shorter than the mean interval. The use of exponential distribution in modeling seismic recurrence is debatable, and currently, there is no much literature and it is still not clear whether seismic recurrence has a strong central tendency or not.

1.3 Recurrence time models

The recurrence time models that can be considered are the exponential, Weibull, gamma, and log-normal distributions. Use of these distributions has been motivated primarily on the grounds of familiarity, simplicity, and convenience.

1.3.1 Exponential distribution

The exponential distribution is the probability distribution of the time between seismic events in a Poisson point process. A number of events in a given time interval follow a poison distribution as it assume that events occur continuously and independently at a constant average rate. This means the time between events follows an exponential distribution. One-parameter exponential distribution has a property of memoryless implying the distribution of a waiting time until a certain event does not depend on how much time has elapsed already. This property disqualifies the exponential distribution as a possible model for seismic recurrence time model as it opposes the elastic rebound model.

1.3.2 Gamma distribution

Gamma distribution is a generalized exponential distribution. This distribution contains the densities of the sum of $n$ independent exponentials. The shape parameter $\mathrm{\theta }>1$ has zero hazard rate at a time zero and increases to a finite asymptotic level that is always smaller than the mean recurrence rate, which is not the case with seismic events. This makes the Gamma distribution not a potential model for seismic recurrence time models.

1.3.3 Weibull distribution

Weibull distribution contains the densities of the minimum of $n$ independent exponential distribution with rate $t$ that have independent occurrence times, and then, the distribution has a Weibull (t, n). It is a distribution for which seismic occurrence rate is proportional to a power of time. The shape parameter, $k$, is that power plus one. A value of $k<1$ indicates that the seismic occurrence rate decreases over time that is there is significant “infant mortality,” or seismic event occurrence is high early and decreasing over time as the seismic events are “weeded” out of the population. This may be because the seismic building force stabilizes over time. A value of $k=1$ indicates that the seismic occurrence rate is constant over time. This is consistent with the elastic rebound model but the Weibull distribution reduces to an exponential distribution. A value of $k>1$ indicates that the seismic occurrence rate increases with time. This happens if there is an increase in seismic building process activities. The hazard-rate functions either start at zero and increase to $\infty$ or vice versa, depending on the parameter k, this disqualifies the weibull distribution as a potential seismic recurrence time model as the elastic rebound model do not support infinite or zero hazard rate.

1.3.4 Log-normal distribution

The log-normal distribution is obtained by taking the exponential of a normal distribution. The asymptotic hazard rate is always zero and hazard-rate functions that increase from zero at time $t=0$ and then eventually decrease to zero. The probability density function puts least weight in the left tail. The mean residual life increases without bound as $\mathrm{t}\to \infty$. Its asymptotic properties disqualify the log-normal family as a reliable seismic recurrence model as this means the longer the time since the last seismic event, the longer it will be expected until the next seismic event. As suggested by Davis et al. [7] that the log-normal model provided only a slightly better fit than the gamma or Weibull models.

1.3.5 The Brownian relaxation oscillator

The Brownian Relaxation Oscillator can be used for seismic forecasting of time of next event. If seismic building process are fixed and tectonic forces load at a constant rate, seismic events will occur after a fixed time interval. The point process will form identical events. The only significant variable in such a deterministic model will be the “strain state,” and after a long time of events, it forms a cycle of loading and instant relaxation oscillating over time forming a deterministic relaxation oscillator model. The strain state will go to zero soon after a seismic event and increase upward at a constant rate up to a fixed elastic limit value. Immediately after exceeding that value an event occurs that relaxes the strain level to zero. The cycle will continue over time making the time of next seismic event predictable.

This is in agreement with the elastic rebound model proposed by Reid [1] with strain level meaning the same as cumulative elastic strain. The strain level could also mean cumulative moment deficit or total stress level. The strain level can be described as the absolute rapture potential. Figure 1 shows a diagrammatic representation of a deterministic relaxation oscillator.

Figure 1.

Deterministic relaxation oscillator.

Let S(t) be the strain in the fault plane at time $t$ measured from a value S₀ after an event. The strain in the fault is the sum of initial strain after an event S₀ and an increase in strain φ due to tectonic loading. The elastic rebound model assumes the strain in the rock load constantly at a rate λ. An event will occur when stain level reaches a constant elastic limits value S_E.

$$\phi =\mathit{\lambda t}$$E1

$$S\left(t\right)=s_0+\mathit{\lambda t}$$E2

Since φ will be set to 0 after an event and starts to increase constantly at a rate λ until the next event. If events occur at constant intervals t_E and if $t$ is the clock time, then

$$S\left(t\right)=s_0+\lambda \left(t-t_E\right)$$E3

$$S\left(t\right)=s_0+\mathit{\lambda t}-\lambda t_E$$E4

The deterministic process can be expressed as

$$S\left(t\right)=s_0+\mathit{\phi t}-{\mathit{\phi t}}_E$$E5

The graphical representation of such a process is shown in Figure 1 aforesaid.

The seismic formation process has other complex variables, which can be represented in the model with a random error term or white noise εt a random perturbation process. The stochastic relaxation oscillator can be expressed as

$$S\left(t\right)=s_0+\mathit{\phi t}-{\mathit{\phi t}}_E+{\epsilon }_t$$E6

The above stochastic relaxation oscillator equations show that the stain level is made up of a sum of three components.

s_o is the initial constant stain pre-existing in a fault plane. The value of s_o is a constant but differs from fault plane to faults depending on the type of rock.

$\mathit{\phi t}-{\mathit{\phi t}}_E$ is the difference between the strain level at any given time and the strain level after the previous seismic event. The value is as a result of stain accumulation due to seismic building process and depends on time since the last seismic event. Seismic forecasting involves being able to find the best probability model that best fits this component. The component is assumed to follow a Brownian Passage-Time Distributions as it depends on time past since the last event. In this model $s_0\le {\mathit{\phi t}}_E$, this explains why some events may result in aftershock or event after of aftershock because in some events the relaxation may not release all the stains in the fault giving a possibility of an event after another event.

The component ${\epsilon }_t$ is a random term that takes into the model the effect of other random variables that affect seismic stain building as a result of seismic perturbation. The term is also called the white noise and follows a Brownian motion or Gaussian or normal distribution with mean zero and constant variance. This term results in the formation of a graph with a path up or above the smooth theoretical graph of S_t.

The graphical representation of such a process is shown in red color in Figure 2. Such a path is like the drunken man’s path. The value of constant variance determines how far the process moves from the line representing the deterministic model.

Figure 2.

Brownian relaxation oscillator with random term.

The seismic strain builds between the time interval $t_0\text{and}{t}_E$ and builds at a rate. This means the seismic recurrence intervals will have average length $\mathrm{\mu }=\frac{\mathrm{\delta }}{\mathrm{\lambda }}$. Indeed, the deterministic oscillator with identical recurrence intervals has a variance, $\mathrm{\sigma }=0.$

1.4 Brownian passage-time distributions

Let T denote the first passage time to level $s>0$ by Brownian motion with drift rate $\mathrm{\lambda }>0$ and diffusion rate ${\mathrm{\sigma }}^2$. The probability distribution of T has a well-known closed form with probability density function given by

$$f\left(t;s;\mathrm{\lambda };\mathrm{\sigma }\right)=\frac{s}{\sqrt{2\mathrm{\pi }{\mathrm{\sigma }}^2{\mathrm{t}}^3}}exp\left(-\frac{{\left(s-\mathrm{\lambda }t\right)}^2}{2{\mathrm{\sigma }}^{2}t}\right),t\ge 0$$E7

The distribution is also called the inverse Gaussian or the Brownian passage time.

The cumulative Gaussian probability function is given by

$$F\left(t\right)=P\left(T\le t\right)$$E8

$$\text{Let}\alpha =\frac{\sigma }{\mu }$$

$$=\mathrm{\Phi }\left\{\frac{1}{\alpha }\left(\sqrt{\frac{t}{\mu }}-\sqrt{\frac{\mu }{t}}\right)\right\}+{\mathit{exp}}^{\left(\frac{2}{{\alpha }^2}\right)}\mathrm{\Phi }\left\{\frac{1}{\alpha }\left(\sqrt{\frac{t}{\mu }}+\sqrt{\frac{\mu }{t}}\right)\right\}$$E9

Using the log-likelihood function is given by

$$L\left(s;\mathrm{\lambda };\mathrm{\sigma }/t_i\right)=\prod _{i=1}^n\left\{\frac{s}{\sqrt{2\mathrm{\pi }{\mathrm{\sigma }}^2{t_i}^3}}exp\left(-\frac{{\left(s-\mathrm{\lambda }{t}_i\right)}^2}{2{\mathrm{\sigma }}^2{t}_i}\right)\right\}$$E10

we find that expectation of $T$, $E\left(T\right)=\mu$ and variance of $T$, $\mathit{var}\left(T\right)={\sigma }^2={\left(\mathit{\mu \alpha }\right)}^2$ where $\mathrm{\alpha }$ is the ratio of the standard deviation to the mean, that is, the coefficient of variation or the aperiodicity of the failure-time distribution. This is also called the ratio of the sample standard deviation to the sample mean. $\mathrm{\alpha }$ is a measure of irregularity in the seismic event sequence and determines the shape of the Brownian Passage-Time Distributions.

1.4.1 Time dependence

According to Matthews et al. [5], the Brownian passage-time distributions quantify occurrence-time probabilities for a steadily loaded system subject to random perturbations. The distribution is used to answer questions like “what is the effect of past time since the last event on conditional probabilities for the next event?”

1.4.2 Hazard rate

The hazard rate of Brownian passage-time distributions is given by

$$h\left(t\right)=\frac{f\left(t\right)}{1-F\left(t\right)}$$E11

1.4.3 Residual life distribution

To examine the effect of elapsed time on occurrence probabilities, we may consider the residual life conditioned on $>t$. The conditional probability of a seismic event at any given time, $t$given that $\text{time}P$ has passed after an event is given by

$$R\left(P\right)=P\left(P\le t+P/P>t\right)$$E12

$$=\frac{F\left(t+P\right)-F\left(t\right)}{1-F\left(t\right)}$$E13

residual life density

$$=r\left(p\right)=\frac{\partial R\left(P\right)}{\partial P}$$E14

$$=\frac{f\left(t+P\right)}{1-F\left(t\right)}$$E15

mean residual life is

$$m\left(p\right)={\int }_0^{\infty }pr\left(p\right)\partial p$$E16

The properties of the Brownian passage-time distributions that makes it appropriate for seismic forecasting is the shape of the hazard-rate function and behavior of the residual life as time increases, asymptotic mean residual life, and hazard rate. The value of $h\left(0\right)$ is also important, as it governs the likelihood of immediate re-rupture after an event. All BPT hazard-rate functions share a general shape, The hazard-rate function, $h\left(t\right)$, of all models in the Brownian passage-time distribution family always starts at 0 at $t=0$, then goes upward a highest value at some time after the probability distribution’s mode, and then goes down approaching a asymptotically a fixed value $h_{\infty }$ as time approaches infinity. The value can be found by taking the limit as $t\to \infty$ as,

$$h_{\infty }=\frac{1}{2\mu {\alpha }^2}$$E17

Brownian failure process eventually attains a quasi-stationary state in which residual time to failure becomes independent of passed time an applicable property in seismic forecasting. The Brownian passage-time distribution describes the failure probability of a Brownian relaxation oscillator as a function of elapsed time and the statistical properties of the failure time series.

Brownian relaxation oscillator can be used to estimate the probability of a seismic event in the next time interval in a given region. Seismic data from the past events can be utilized to estimate the parameters $s;\mathrm{\lambda };\mathrm{\sigma }$, which differs from region to region. Once the parameter has been estimated a good approximation of $t$, the time until the next event can then be calculated. This is the proposed model for seismic forecasting.

1.5 Model parameter estimation

1.5.1 Maximum-likelihood estimates

Suppose we have a consistent, complete, and homogeneous past seismic data of time between seismic events, $t_1;t_2;t_3;\dots ..t_n$ of a given region. The estimates of model parameters s; λ and σ can be found by looking for values that maximize the parameters using the sample data. Such values can be found by equating to zero the first derivatives of f(t;s;λ;σ) with respect to each parameter.

$$\frac{\partial f\left(t;s;\mathrm{\lambda };\mathrm{\sigma }\right)}{\partial s}=\prod _{i=1}^n\frac{\partial \frac{s}{\sqrt{2\mathrm{\pi }{\mathrm{\sigma }}^2{t_i}^3}}exp\left(-\frac{{\left(s-\mathrm{\lambda }{t}_i\right)}^2}{2{\mathrm{\sigma }}^2{t}_i}\right)}{\mathrm{\partial s}}=0$$E18

$$\frac{\partial f\left(t;s;\mathrm{\lambda };\mathrm{\sigma }\right)}{\mathrm{\partial \lambda }}=\prod _{i=1}^n\frac{\partial \frac{s}{\sqrt{2\mathrm{\pi }{\mathrm{\sigma }}^2{t_i}^3}}exp\left(-\frac{{\left(s-\mathrm{\lambda }{t}_i\right)}^2}{2{\mathrm{\sigma }}^2{t}_i}\right)}{\mathrm{\partial \lambda }}=0$$E19

$$\frac{\partial f\left(t;s;\mathrm{\lambda };\mathrm{\sigma }\right)}{\mathrm{\partial \sigma }}=\prod _{i=1}^n\frac{\partial \frac{s}{\sqrt{2\mathrm{\pi }{\mathrm{\sigma }}^2{t_i}^3}}exp\left(-\frac{{\left(s-\mathrm{\lambda }{t}_i\right)}^2}{2{\mathrm{\sigma }}^2{t}_i}\right)}{\mathrm{\partial \sigma }}=0$$E20

The aforesaid equations will give maximum-likelihood estimates of s;λ and σ that can be used in the probability density function $f\left(t;s;\mathrm{\lambda };\mathrm{\sigma }\right)$ to estimate the variable time $t$ of future events.

1.5.2 Magnitude forecast

The time of a seismic event is estimated from the above model $f\left(t;s;\mathrm{\lambda };\mathrm{\sigma }\right)$, but every seismic event must have a forecasted value of the expected magnitude. Suppose we have consistent, complete and homogeneous past seismic data of magnitudes $m_1;m_2;m_3;\dots ..m_n$ of a given region.

The unbiased estimate of the future expected magnitude is the mean of past magnitude

$$m_{\mu }=\frac{{\sum }_{i=1}^n{m}_i}{n}$$E21

From 1.5.1 and 1.5.2 above, the time and magnitude of future events can be forecasted.

2. Conclusion

The Brownian relaxation oscillator with random term for seismic random perturbation modeled above represents a model that can be used for seismic forecasting. The inclusion of the error term in the model gives allowance of a deviation of the forecast from the actual event. Since the error term is known to be normally distributed with mean zero and constant variance, the expectation of the deviation is zero. The error term also represents other seismic formation variables and errors in model parameter estimation.

For one to accurately use the model, a consistent, complete, and homogeneous with unified magnitude, past seismic data are needed. Such data are always unavailable, because a complete catalog of the full population of events begins from the start of the earth planet formation. Statistical method utilizes a very small sample beginning “yesterday” when seismic monitoring started to infer population parameters. If such model parameters can be accurately estimated, seismic events can be accurately forecasted using the Brownian relaxation oscillator with random perturbation.

Conflict of Interest

The author declares no conflict of interest.