The Self-Organized Criticality and Periodicity of Temporal Sequences of Earthquakes

A new cellular automaton model has been developed to examine the nature of temporal sequences of earthquakes. The model takes the space dependence of strength on inhomogeneous faults into account and assumes that an earthquake produces a continuously connected area with uniform stress drop. In the model viscous slip is also introduced so that earthquakes may be realizable only when the applied stress accumulates faster than the viscous relaxation. The analysis reveals that the sequences of earthquakes commonly satisfy the power law relation between the intensity and frequency of earthquakes so that earthquakes may be in the state of self-organized criticality. On the other hand, periodicity appears in some sequences that consist of the groups of high seismic activity repeated between calm intervals at an almost constant period. Therefore, self-organized criticality and periodicity coexist in these sequences and the claim that earthquakes are unpredictable because of self-organized criticality may be inadequate.


Introduction
Earthquakes are caused by slips along inhomogeneous faults due to enhanced stress in the surrounding rock.The occurrence of earthquakes involves the fracture of faults and rock, the friction of faults, and the elastic deformation and wave emission induced in rock.A complete formulation that describes all these factors satisfactorily is a complex and di cult problem.So, the processes generating earthquakes have been analyzed in a simpli ed way that focuses on some of these factors.
Using a cellular automaton model, Bak et al. (1987) analyzed sand pile collapse and found a power law relation between the size and frequency of collapses.This analysis led to an important concept of selforganized criticality.Earthquakes are regarded as a representative phenomenon of self-organized criticality in which the Gutenberg-Richter law gives the power law relation between the energy and frequency of earthquakes (Bak and Tang 1989).It is reminded, however, that the concept of selforganized criticality still has something ambiguous (Sornette 1994, Flyvbjerg 1996).
The cellular automaton model for earthquakes have been developed by Ito and Matsuzaki (1990), Nakanishi (1990), and Olami et al. (1992) to have good analogy with the dynamic system consisting of frictional blocks connected with springs.This dynamic system was originally proposed by Burridge and Knopoff (1967).The cellular automaton model by Olami et al. (1992) has been widely applied to various studies of earthquakes, including the conditions of criticality (Carvalho and Prado 2000, Miller and Boulter 2002, Main and Naylor 2010), the occurrences of foreshocks and aftershocks (Hergarten andNeugebauer 2002, Helmstetter et al. 2004), and the comparison with real earthquakes (Peixoto and Prado 2004).
Whether or not earthquakes are predictable is an important problem that is closely related to the selforganized criticality.Bak and Tang (1989) claims that earthquakes are unpredictable because their occurrences are controlled by minor detail far from the origin in the state of criticality.Caruso et al. (2007) agrees with the unpredictability of magnitudes of the next earthquakes using the probability density function.
On the other hand, it has been demonstrated by an experiment for beads that large avalanches of piles can be predicted by detectable variation of internal structure (Ramos et al. 2009).Based on the rst return-time probability of earthquake catalogues, Yang et al. (2004) insists that earthquakes are unlikely phenomena of self-organized criticality.The analysis of time sequence and entropy change of earthquakes before the Tohoku earthquake in 2011 suggested the presence of precursory phenomena of the earthquake (Varotsos 2020).
In connection with the predictability, a cellular automaton model that has a better resemblance with the spring-block system displays a quasi-periodic behavior of earthquakes (Ramos et al. 2006).The periodic occurrence of earthquakes is a basic process of stress relaxation associated with plate motion in the classic concept of elastic rebound theory (Reid 1910).In the present paper the periodicity of earthquakes is examined in more detail using a new cellular automaton model that assumes space-dependent strength of faults and uniform stress drop produced by an earthquake.

Model Of Earthquake Occurrences
For our study of the sequences of earthquakes a new cellular automaton model has been developed.The model has a similar structure to those used by the previous studies such as Olami et al. (1992) and Ramos et al. (2006), but the similarity to the spring-block system is regarded as less important because of a great difference of the system from real faults surrounded by continuous rock.
In our model the stress σ is introduced as a main variable and allocated to each point (x i , y i ) of a twodimensional lattice with a constant distance along the x and y directions.More strictly, σ is the component of stress tensor in the major slip direction.In fact, σ is equivalent to the total force in the previous models (Olami et al. 1992, Ramos et al. 2006) but the nature as a local variable independent of neighbors is more emphasized.It is assumed that σ is increased with time at a uniform and constant rate r by some external actions like plate motion.
To represent the inhomogeneity of faults each lattice point is assumed to have its own strength s as a local property.When σ exceeds s at one of the points a slip takes place and drops σ to zero at this point.Due to the elastic response of surrounding rock the stress drop is assumed to enhance the neighboring stresses at four adjacent points (x i−1 , y j ), (x i+1 , y j ), (x i , y j−1 ) and (x i, y j+1 ) by α times the dropped stress.The factor α should be 1/4 if the total stress is conserved during the slip.As the boundary condition, it is assumed that no stress drop is applied from the outside of the system.When the stress that has been enhanced by an adjacent slip again exceeds the local strength another slip is induced there.Slip may be repeated in this way as a domino.An earthquake is the result of all these sequential slips that form a cluster of m continuous points with the stress of zero.The value of m that is proportional to the fractured area represents the intensity of this earthquake.
It is assumed in our model that every point that has been slipped during the same earthquake should have a free stress of zero.This means that the stress enhancement rule due to adjacent slip is not applied to the points that have been subjected to the stress drop at the same earthquake.This assumption seems to give a better representation of real earthquakes but has not been taken in most previous models.
Viscous slip is also introduced to include the effect of inelastic deformation and the smaller earthquakes than the lattice.If the viscous slip can be formulated by the Maxwell model with relaxation time τ the stress change during a calm period without earthquakes from time t 1 to t 2 is given as 1 Here σ 1 and σ 2 are the stresses at the times t 1 and t 2 , respectively.
It is reminded that the stress enhancement rule due to the adjacent slip is not applied to the points that have been fractured by the same earthquake.This exception prevents the total stress and total energy from being conserved even in the case of α = 1/4.In fact, the stress and energy are lost continuously by viscous slip and discontinuously by the stress drop and elastic wave emission associated with earthquakes.So, the conservation of stress is not rigorously required and regarded as not important in our study.
3 Results Of Simulation

The scheme of simulation
In the numerical simulation the permissible range of strength s is rst prescribed as s l to s u , and s is randomly distributed within the speci ed range over all the lattice points.The initial values of stress σ are randomly given between 0 and the lower limit s l of the strength.Then σ is increased at the rate r with a small step of time.At each step σ is modi ed by the viscous slip following Eq.( 1) and compared with s.If σ exceeds s at one of the lattice points the total area subjected to the stress drop at this event is de ned and the intensity m is evaluated.
Dimensionless variables are introduced in the numerical calculation by scaling time t by τ, stress σ and strength s by s l , and stress increase rate r by s l /τ.The variable parameters are thus the ratio s u /s l and dimensionless rate of r other than the random distributions of strength and initial stress.In most results shown below the stress change is calculated over 101x101 lattice points with the time step of 10 − 6 τ.It has been con rmed that the time step that is as small as this value little affects the results.The factor of stress enhancement due to adjacent slip is xed as α = 1/4.
Some typical results of the simulation are given in the following sections.The rst two examples are the results for relatively inhomogeneous and homogeneous faults with s u /s l = 2 and 1.02, respectively.A common stress increase rate of r = 2s l /τ is assumed there.Then the effect of various stress increase rate is examined.For the convenience of comparison, the random values of strength and initial stress are kept same for the examples presented here even if qualitatively similar results are obtained for other random distribution.

Results for relatively inhomogeneous faults
As the example for relatively inhomogeneous faults, the sequence of earthquakes is analyzed for a larger strength range with s u /s l = 2.The assumed stress increase rate is r = 2 in the unit of s l /τ.The seismic intensity m obtained in this simulation is plotted at the outbreak time t of earthquakes in Fig. 1a.In this sequence earthquakes with various intensities happen almost continuously.It is noted that the overall activity apparently differs before and after the time of about 2. Namely, the periods with high and low seismic activities are more contrasted in the earlier stage.
To examine the singularity of this sequence, the number n of the earthquakes having intensity m is plotted for the earlier (t = 0-2.5)and later (t = 2.5-5) stages in Figs.1b and 1c, respectively.In these plots some data that may have the value of m too large to stay within a continuous range are excluded.The two plots give linear relations between m and n in logarithmic scales so that the following power-law may hold.

2
The index B that is determined with a least-square tting is given in the gures.
Figure 2 gives the snapshots of stress distribution for some selected times.At each point of the snapshots darker colors display higher stress values.At t = 0 stress is randomly distributed among individual lattice points but the distribution changes with repeated earthquakes.Namely, wider areas with the same color arise and spread out.Each of these areas has experienced the same earthquake at the time that is older with darker colors.In this gure wide areas become more abundant during the earlier stage to approach a stable state in the later stage.So, the earlier stage may be a transient or preparatory period before the fault condition is stabilized.
The periodicity of this sequence is examined separately for the earlier and later stages in Figs. 3 and 4, respectively.For this purpose, the sequence of earthquakes that collects the point data of m is rst smoothed out and the average m i of the segment i divided by a common interval of Δt = 0.02 is obtained (Figs.3a and 4a).The amplitude spectrum of this smoothed data is calculated in Figs.3b and 4b with the fast Fourier transform (FFT) using 128 points (125 data points plus the last 3 vacant points of zero).Furthermore, normalized autocorrelation C k is calculated as a function of the time shift kΔt by The autocorrelation is given in Figs.3c and 4c.
If the spectrum has a dominant peak the sequence has a periodic nature with the period of this peak.
Similarly, if the autocorrelation has repeated peaks the sequence has the period equal to the time-shift of the rst peak.The sequence of the earlier stage has a period of about 0.6 consistently with all the time series, spectrum, and autocorrelation in Fig. 3.However, the periodic feature is only transient in the time series.The transiency is supported by rather broad peak of the spectrum and decaying peaks of the autocorrelation.This transient state moves to the later stage that has no clear indication of periodicity in Fig. 4. It is thus considered that singularity is a dominant nature of inhomogeneous faults even if periodicity may appear temporarily.

Results for relatively homogeneous faults
As the example of relatively homogeneous faults, the sequence of earthquakes is given in Fig. 5a for a smaller range of strength of s u /s l = 1.02 with the same stress increase rate of r = 2.In this sequence the group of high seismic activity is repeated many times between calm intervals with an almost constant period of about 0.3.So, the periodicity must be a central feature of this sequence.It is noted, however, that the sequence including the intensity of the greatest earthquake is not same but signi cantly varies among individual groups.In Fig. 5b the number of events meets a linear relation in logarithmic scales and the power law relation ( 2) is still applicable to this case.
Figure 6 gives several snapshots of stress distribution within the same group of high seismic activity around t = 1.2.The rst three snapshots tend to become darker with time as a whole so that the overall stress is getting higher toward the greatest activities containing the widest area in the fourth snapshot.It is noticed there that some traces of big earthquakes are surrounded by several small earthquakes that take place immediately before.After the overall stress is relaxed it again begins to increase.In this way the mean stress level changes cyclically with time The periodicity of the sequence is con rmed more quantitatively in Fig. 7.The spectrum has a sharp peak at the period of 0.28 and the autocorrelation gives repeated peaks that decays slowly.
It is emphasized that criticality and periodicity coexist in this case.The case of the complete homogeneous fault with s u /s l = 1 produces a similar sequence of earthquakes having both criticality and periodicity.

The effect of stress increase rate
When the stress increase rate r is as small as 1 the increased stress is entirely consumed by viscous slips and no earthquakes happen.So, earthquakes are realizable only in the condition that stress is applied su ciently fast.The critical rate of r above which earthquakes appear is numerically estimated to be about 1.
Figure 8a gives the sequence of earthquakes for an intermediate rate of stress increase, r = 1.5, with the condition of relatively inhomogeneous fault, s u /s l = 2.This sequence is again in the state of criticality as is revealed in Fig. 8b by the relation between the intensity and frequency of earthquakes with a little smaller value of B. In Fig. 9 the sequence consists of repeated groups of high seismic activity with the period of about 1.1 consistently with the spectrum and autocorrelation even if the periodic nature may be transient.

Self-organized criticality and periodicity
The sequences of earthquakes that have been obtained in this paper commonly satisfy the power law relation so that earthquakes may be in the state of self-organized criticality.On the other hand, some earthquake sequences contain a periodic component of temporal variation.Therefore, self-organized criticality and periodicity coexist in these examples.
It is told that self-organized criticality is characterized by power law distribution and fast decay of temporal and spatial correlations (Sornette 1994, Flyvbjerg 1996).In fact, a phenomenon that involves power law distribution is usually regarded as the system of self-organized criticality.Our analysis thus suggests that earthquakes could be in the state of self-organized criticality in a usual sense.However, the examples with periodicity may disturb the condition of fast decay of temporal correlation and may be placed outside self-organized criticality.
Apart from the de nition of self-organized criticality, it is practically more important that some earthquakes have periodic nature simultaneously with the power law relation that have been rmly established as the Gutenberg-Richter relation.The periodicity can be directly applied to earthquake prediction and may relate to some other regular properties of earthquakes.
In the examples with periodicity a big earthquake happens within the group of high seismic activity that contains some small precursory events.An analysis of real sequences of earthquakes in Tohoku and Kanto areas, Japan, have revealed that the small earthquakes before a big hazardous earthquake can be identi ed using a self-organized map (Ida and Ishida 2022).Therefore, the claim that earthquakes are not predictable because they are statistical phenomena associated with self-organized criticality (Bak and Tang 1989, Geller et al. 1997, Caruso et al. 2007) may be inadequate.

Implication for real earthquakes
The energy relaxed by an earthquake is proportional to the fractured area (m in our model) multiplied by the mean displacement of slip.Since the mean displacement is approximately proportional to the square-root of the slipped area the energy becomes proportional to m 3/2 .If this approximation is employed the power-law relation ( 2) turns out to be equivalent to the empirical relation between the energy and frequency of earthquakes that is derived from the Gutenberg-Richter law with the b-value equal to index B.
According to the results of our simulation relatively homogeneous faults tends to have smaller values of B compared with inhomogeneous faults.This tendency is consistent with a generally conceived idea of the b-value, but estimated values of B seem to be rather higher than most of the b-values observed for real earthquakes.The discrepancy between calculated and real values is probably caused by the strength that has been distributed at random among individual lattice points in our model.The scales of strength distribution in the model may be too small to represent the inhomogeneity of natural faults.
Our simulation suggests the presence of various types of fault activity, including those with wide ranges of magnitudes, periodic earthquake occurrences, and inelastic slips without seismicity.Some real faults display similar variability as seen in San Andreas Fault of California (U. S. Geological Survey 2004).In Park eld, the middle of the fault, earthquakes of magnitude about 6.0 have taken place regularly since 1857 with a period of 22 years and an earthquake prediction experiment was attempted (Roeloffs and Langbein 1994).In its northwestern part called creep zone the stress accumulation due to plate motion is relaxed by inelastic slip without earthquakes.Further north is occupied by the area of high seismic activities with wide ranges of magnitudes in which the great San Francisco Earthquake of magnitude 7.8 happened in 1906.
It will be a future interesting problem to compare the simulation results with real earthquakes in more details.

Declarations
The author has no competing interest and received no funding to declare for this work.
The sequence of earthquakes for the relatively inhomogeneous fault with s l /s l = 2 and the dimensionless stress increase rate of r = 2.The intensity m is plotted at the outbreak time t of earthquakes (a).The number of earthquakes n having intensity m is plotted in logarithmic scales for the earlier (b) and later (c) stages.
Stress distributions at some selected times for the relatively inhomogeneous fault with s l /s l = 2 and the dimensionless stress increase rate of r = 2. Stress is displayed over 101x101 lattice points by darker colors with higher stress values.

Figure 3 The
Figure 3

Figure 4 The
Figure 4

Figure 5 The
Figure 5

Figure 7 The
Figure 7