Cellular Automaton Equivalence of Megathrust Earthquake Process

This paper links the seismic process of a megathrust earthquake with a cellular automaton (CA) and paves the way for understanding the seismic process by investigating CA. During the seismic process of the Great East Japan Earthquake (GEJE) of magnitude 9, anomalous noise with large amplitude in the frequency range of 3 Hz to 10Hz is observed in the Fourier spectrum. The anomalous noise and the normal ground vibrations in the GEJE process are generalized by introducing (cid:0)-tremor, which is dened as the curvature of the Fourier amplitude spectrum of the velocity. It is found that the binarization of the vibration velocity conserves the (cid:0)-tremor, and the binarized velocity and velocity are equivalent. The binarization allows the recorded data to be directly compared with CA since CA is inherently binary, and since both the binarized velocity deviation and the CA are considered to follow the same master equation which represents the thermodynamics of a stochastically uctuating non-equilibrium system. It is found that the recorded data in the GEJE process is a subset of the CA and the former can be reproduced by the latter.


Introduction
Strong earthquakes are a major concern in disaster management, and various measures are being taken for strong earthquakes. Earthquake Early Warning system in Japan warns people when an earthquake of 5 or greater is expected on the Japan seismic scale. When an earthquake is detected, the system analyzes the data captured by seismographs near the epicenter to estimate the epicenter, the magnitude of the earthquake and the seismic intensity. The estimated information is quickly released so that people can move to safe places or evacuate from dangerous places before strong surface waves arrive. Regarding building regulations, the seismic standards of the Building Standard Law in Japan require minor damage in medium-scale earthquakes with a seismic intensity of 5 or greater, and no collapses in large-scale earthquakes with a seismic intensity of 6 to 7.
On the other hand, earthquakes generally last less than a minute, and the dominant state of ground motion is seismically silent. Therefore, in order to understand the seismic process, it is necessary to investigate the silent state. Nonvolcanic tremor is one of the notable discoveries regarding the silent state. Obara investigated the seismically silent period in southwest Japan and identi ed the nonvolcanic tremor, the weak but noticeable signal with typical frequency range from 1 Hz to 10 Hz (Obara, 2002). In 2003, Rogers and Dragert related tremors to ground slip events. Tremor activity accompanied by a slip event was observed approximately every 12 months for 6 consecutive years at Cascadia subduction zone interface (Rogers & Dragert, 2003). Regarding the mechanism of the long duration tremor, Peng and Chao observed the tremor induced by an earthquake and discussed that tremor occurred as a simple frictional response to the driving force (Peng & Chao, 2008). This paper focuses on the weak ground vibrations of micron/second scale in the seismic process of the Great East Japan Earthquake (GEJE) of magnitude 9 occurred at the epicenter of latitude 38.06N, longitude 142.51E, and depth 24 km at 14:46 on March 11, 2011. The duration of the seismic process under investigation is from 2006 to 2016. The rest of the paper is organized as follows: First, the GEJE process is characterized. The dynamic parameter which describes the unique feature of GEJE is de ned and extended to include the stable and normal ground vibration. Second, the de ned dynamic parameter is shown to be conserved in the binarization of ground vibrations. The ground vibration and the binarized ground vibration are shown to be equivalent. Third, the binarized ground vibration is considered as a stochastically uctuating non-equilibrium thermodynamic system that satis es the master equation. A cellular automaton (CA) is also considered as a system described by the same master equation. Finally, the ground vibration is reproduced by CA, and has been shown to be equivalent to CA.

Characterization Of Geje Process
The ground vibration of the GEJE process is characterized, by observing the anomalous noise with large amplitude in the frequency range of 3 Hz to 10Hz in Fourier spectrum.

Observation of ground vibration signals
Ground vibration velocity data acquired every 0.05 seconds at the seismic station KSN is downloaded in chronological order from the web site of F-net, broadband seismograph network of National Research Institute for Earth Science and Disaster Resilience (NIED, 2019). KSN is located at latitude 38.98N, longitude 141.53E, and altitude 260 m. The downloaded data is converted to piecewise deviation which is the difference between the velocity and the velocity averaged over the subsequent 10 data points. The piecewise deviation uctuates around zero, and its squared average is the dispersion in statistics. The piecewise deviation data is divided into blocks of 1024 data, which corresponds to the data acquisition time of 51.2 seconds, and the Fourier amplitude of each block is calculated. The Fast Fourier Transform (FFT) algorithm is applied with no overlap, and no ltering. The upper bound of the frequency domain is 10Hz, which is half the data acquisition frequency. The lower bound is 0.02 Hz which is determined by the block size 1024. Therefore, the FFT with the sampling frequency of 20Hz and the block size of 1024 is equivalent to an FFT with a 0.02-10 Hz bandpass lter. Fig. 1 shows a comparison of the velocities and spectrograms in the up-down (UD), north-south (NS), and east-west (EW) direction. The velocity data was acquired at KSN every 0.05 seconds from Mar. 3, 2011to March 11, 2011. The period includes the magnitude 9 Great East Japan Earthquake (GEJE) occurred at 14:46 on March 11, 2011. In the spectrogram range from 1 Hz to 10 Hz, there are noticeable signals shown as the vertical brown lines. In the quiet period before the earthquake of magnitude 7.3, the timing of the vertical brown lines in the spectrograms ( Fig. 1 (a4), (b2), and (c2)) respectively matches the timing of the wave clusters which have larger amplitude than surroundings ( Fig. 1 (a3), (b1), and (c1)). Since the UD component contains greater number of vertical brown lines than the other components, we focus on the UD component in the later sections.
The third brown line in Fig. 1 (a4), which corresponds to the velocity deviation in Zone_A in Fig. 1 (a3), constructs a ner spectrogram structure. Fig. 2 shows the details of the Zone_A of 12500 second duration. The velocity deviation and its spectrogram are shown in Fig. 2 (a) and 2 (b), respectively. Fourier amplitude spectrum and its 10 moving averages are respectively indicated by the black and red lines in the Log10-Log10 plots of Fig. 2 (c1) to 2 (c6). The velocity deviations in Fig. 2 (d1) to 2 (d6) are the source data for the Fourier spectrum. The number below the velocity deviation graph indicates the time interval (seconds x 20). The velocity deviations are extracted from the beginning, center, end, and their intermediates of the period shown in Fig. 2 (a), and are chronologically exhibited from left to right. The rst and last Fourier spectra show small negative curvatures in the range 1 Hz to 10 Hz (Fig. 2 (c1) and 2 (c6)). The rest of the spectra show large values and negative curvatures in the range from 1 Hz to 10 Hz (Fig. 2 (c2) to 2 (c5)). The curvature widens in the center of the zone and narrows in the rest of the zone.
The amplitude of the velocity deviation is small at the beginning and end, and large in the central zone ( Fig. 2 (d1)-(d6)).

De nition of α-tremor
The curvature of the Fourier amplitude spectrum is de ned as the ratio of P ni − P i to | P 2 − P 1 | in Fig. 3 (a), where P i is the point of the (frequency, spectrum) coordinate system. The frequency of P 1 and P 2 are 2.97 Hz and 9.8 Hz, which correspond to the 152th and 502th point on the frequency axis, respectively. The average of the spectrum value of the nearest 5 points are assigned as the spectrum value for the P 1 and P 2 . P i is a point in the 2.97-9.8 Hz range. P ni is determined so that the line from P i to P ni is perpendicular to the line connecting P 1 and P 2 . If the spectrum value of P i is greater than that of P ni , the curvature is negative. Otherwise, the curvature is non-negative. The curvature is independent of the scale change since the line length in log10 plot is invariant to the scalar multiplication of the coordinate values.
We de ne α-tremor as the product of "-1" and the curvature of which absolute value is greater than the absolute value of other curvatures in the frequency range of 2.97 to 9.8 Hz ( Fig. 3 (a)). An arbitrary ground velocity signal is classi ed as either positive α-tremor or non-positive α-tremor.
The α-tremor for the velocity data acquired at KSN during March 03, 2011 to March 11, 2011 is exhibited in Fig. 3 (b). As expected, the timing of the positive peaks of α-tremor is similar to the timing of the brown line in Fig. 1 (a4).
It should be noted that the piecewise velocity deviation is equivalent to the raw velocity data in evaluating the α-tremor. Fig. 4 (a) compares the Fourier amplitude spectrum of the velocity deviation data to the Fourier amplitude spectrum of the raw velocity data. In the range of 2.97 Hz to 9.8 Hz, the Fourier spectrum of the deviation velocity (black line) matches the Fourier spectrum of the raw data (green line) by 80%. Therefore, we may select either piecewise deviation velocity data or raw velocity data to obtain the characteristic Fourier spectrum in the range 2.97 Hz to 9.8 Hz. The orange and red lines in Fig. 4 (a) are the 10-moving averages of the black and green lines, respectively. The source data of the spectrum, which are the velocity deviation and velocity acquired at KSN during the period from March 1, 2012 to March 10, 2012, are shown in Fig. 4 (b) and 4 (c), respectively.

Seismological signi cance of α-tremor
In the previous sections, we have discussed the uctuation of ground vibrations for the data recorded by F-net's broadband seismograph (NIED, 2019) installed at KSN. In order to understand the relation between α-tremor and seismic events, α-tremor is re-examined for the velocity data recorded by the high-sensitivity seismographs of Hi-net (NIED-2, 2019) consisting of nearly 800 stations with an average spacing of 20 km. The Hi-net seismographs are installed at the bottom of boreholes at a depth of 100-3500 m to reduce the noise generated by winds, ocean waves, and human activity. The natural frequency of the seismograph is 1 Hz. The locations of the three Hi-net seismic stations near the KSN, namely IWTH27, MYGH03, and IWTH18, are shown in Fig. 5.
The Fourier amplitude spectrum of the velocities recorded from 21:18 to 11 minutes on 2011/03/06 are compared between KSN and IWTH27 in Fig. 6. Since the two seismic stations are 6 km apart and close to each other, the signals arriving at the stations are similar. The negative curvature of the spectrum in the range 2.97 Hz to 9.8 Hz, positive α-tremor, is observed from 21:18 to 21:24 (blue rectangle in Fig. 6 (a1)). The velocity plots corresponding to α-tremor are relatively dense, as shown in the blue rectangle in Fig. 6 (a2). In this gure, we do not see seismic p-waves, seismic s-waves, and the nonvolcanic tremor that is detectable by plotting the envelope which is the root mean square trace of the 2 Hz-16 Hz bandpassltered velocity data (Obara & Hirose, 2006). Therefore, it is appropriate to think of positive α-tremor as a kind of noise signal rather than an earthquake or nonvolcanic tremor.
The positive α-tremor is excluded from Hi-net, and not found in the spectrum of IWTH27 ( Fig. 6 (b1)). However, instead, Hi-net detects an earthquake from 21:26 to 21:27, one minute after the end of the series of positive α-tremor ( Fig. 6 (b2)), suggesting a relation between positive α-tremor and the earthquake.
In order to understand the positive α-tremor in terms of seismic events, microearthquakes occurred in the neighborhood of the positive α-tremor are investigated. The top plot in Fig. 7 (a) shows two microearthquakes E1 and E2 in the seismogram recorded at Hi-net IWTH27 for 12540 seconds starting from 2011/03/06 20:18. The time range is similar to the range discussed in Fig. 2 (a) and (b), during which the clusters of positive α-tremor are observed at F-net KSN. The microearthquake E1 is the earthquake occurred one minute after the end of the positive α-tremor cluster shown in Fig. 6 (a1). At the timing of E1 and E2, no earthquake is recorded at F-net KSN, as shown in the bottom plot in Fig. 7 (a).
By applying a bandpass lter in the range of 4 Hz to 9.5 Hz, relatively high-amplitude and low-amplitude clusters emerge in the F-net KSN seismogram ( Fig. 7 (b) top). The rst high-amplitude cluster shown at the left side of a red vertical line in the middle plot of Fig. 7 (b) corresponds to the rst cluster of positive α-tremors. One minute after the end of the cluster, the microearthquake E1 occurs as shown in the bottom plot of Fig. 7 (b). It should be noted that the microearthquake E2 disappeared after applying the bandpass lter (bottom plot in Fig. 7 (b)). As shown in Fig. 7 (d1), E2 has no signi cant frequency component from 4 Hz to 9.5 Hz in the Fourier amplitude spectrum. E2 has the frequency components primarily around 1 Hz, which is dissimilar to the positive α-tremor spectrum in Fig. 6 (a1).
The Fourier amplitude spectrum of the microearthquake E1 shows that it has wide range of frequency components from 1 Hz to 10 Hz ( Fig. 7 (c1)), which is qualitatively similar to the Fourier spectrum of the positive α-tremor in Fig. 6 (a1). The similarity in frequency components implies similarities in material composition, system size and dynamics between E1 and the positive α-tremor.
Since the source velocity data of the E1 spectrum clearly shows p-wave and s-wave ( Fig. 7 (c2)), the epicenter of E1 is estimated. The epicenter of E1 is graphically estimated as the intersection of three circles, which are de ned by the radius D calculated by Omori's formula: D = kT, with the locations of three seismic stations as the origin. Where T is the difference between the arrival times of p-wave and swave at the station, and k is an empirical factor equal to 8 (Kato & Okamoto, 2016). Fig. 8 (a) shows that the estimated E1 epicenter is close to the GEJE epicenter, suggesting that E1 belongs to the events in the seismic process of GEJE.
In summary, the positive α-tremor is an anomalous noise, but it occurs one minute before the microearthquake E1 and has the Fourier frequency component similar to E1 whose epicenter is close to the GEJE epicenter. Consequently, α-tremor is an important and characteristic constituent factor of the seismic process of GEJE.

Binarization Of Velocity Deviation Data
The velocity deviation data is binarizable without losing the α-tremor property. In Fig. 9, the Fourier amplitude spectrum and spectrogram calculated from the velocity deviation are compared to those calculated from the binarized velocity deviation. Fig. 9 (a1) and 9 (b1) shows the binarization procedure. If each velocity deviation data in Fig. 9 (a1) is greater than the mean of the data set under consideration, the deviation data is converted to 1, otherwise the deviation data is converted to 0. The binarized data can be expressed as the time sequence of 0 and 1 as shown in Fig. 9 (b1). The clear negative curvature in the frequency range 1 Hz to 10 Hz, shown in both the Fourier spectrum of the velocity deviation and the binarized data, implies that the α-tremor is preserved in the binarization ( Fig. 9 (a2) and (b2)). Fig. 9 (a3) is the Fourier amplitude spectrogram duplicated from Fig. 2 (b), of which source data is the velocity deviation shown in Fig. 2 (a). The source data is binarized and its spectrogram is calculated as shown in Fig. 9 (b3). The qualitative similarity between the spectrogram of the binarized data and the spectrogram of the source data suggests that the α-tremor is preserved in the binarization ( Fig. 9 (a3) and (b3)). Therefore, the binarized velocity and velocity deviation are equivalent as long as the ground vibration is considered as a α-tremor uctuation.

Non-equilibrium Thermodynamics Of Ground Vibrations
Vibrational states are de ned by the binarized velocity deviation, and the stochastic dynamics of transition of the state are considered to be described by the master equation which represents the thermodynamics of a stochastically uctuating non-equilibrium system. The thermodynamics of ground vibrations in the GEJE process is investigated.

De nition of ground vibration state
Since the α-tremor is conserved in the binarization of the velocity signal, the essential of the ground motion is the distribution of the signal rather than the shape of the signal. Therefore, it is reasonable to de ne the ground vibration state at the speci ed time interval by counting the "1" clusters of the binarized velocity deviation in the time interval. In de ning the vibration state, the binarized velocity sequence ( Fig.  10 (a)) is divided into blocks with 10 data points, and the number of clusters of "1" is counted in each block. In order to preserve the total number of the cluster, the rule shown in Fig. 10 is applied. In the 10data block, we scan the cell from left to right and count one if the sequence of "10" is found. At the end of the scan, at the 10th data point, we count one only if the 11th data point is "0" (Fig. 10 (b)). The counting rule restricts the maximum number of clusters in a block to ve, and de nes ve vibrational states s 1 , s 2 , s 3 , s 4 and s 5 , each containing 1, 2, 3, 4, and 5 clusters ( Fig. 10 (c1)-(c5)).
The number of data points 10 per block is determined by examining samples of binarized velocity data. If a block of a particular size is completely occupied by 1s, then the number of clusters of 1s in the block is 1. If this is the case for all blocks, no uctuation in state can be detected. The block size needs to be increased to detect the characteristics of the state. If the blocks are very large and each block contains all possible patterns of 0s and 1s, then all the blocks will look similar and no uctuation of state will be detected. In this case, the block size must be reduced. After examining a few cases, it is found that the 10 data points per block is adequate to preserve the characteristic of the state change.
Thermodynamics of the ground vibration Fig. 11 (b) shows the rst 100 data of the binarized velocity deviation data of Fig. 9 (b1). The 10 data in each row of Fig. 11 (b) are the binarized velocity deviation data splitted from Fig. 11 (a), and constitute the vibration state of a time interval of 0.5 seconds. In general, the vibration state shows a different pattern of binary sequence for each row and contains a different number of clusters for each row (Fig. 11  (c)). Since each row corresponds to a different time, the state of ground vibration uctuates over time.
The time series of the number of clusters in Fig. 11 (c) shows the history of the state transitions. Since the number of clusters in a state is de ned as a state index, the square brackets that make up the pair of two numbers indicate that the state of the number in the lower row has transitioned to the state of the number in the upper row. The transition rate matrix W ij de nes the total number of transitions from istate to j-state so that the W ij count is incremented by 1 when a transition from i-state to j-state occurs ( Fig. 11 (d)). The result of the W ij counting for the 100 data is shown in Fig. 11 (e). Fig. 11 (f) shows the unnormalized probability density vector, of which component p i is the total number of i-state.
The state of ground vibration, which uctuate over time, implies that the state is non-equilibrium. It is known that the thermodynamics of a uctuating nonequilibrium system are described by the master equation (Eq. (1)), and the entropy production rate (Eq. (2)) which is similar to the entropy of the second law of thermodynamics of equilibrium systems (Haitao, Y. & Jiulin, D., 2014).
where W ij and p i are coherent with those in Fig. 11. J ij and F ij are called the ow from i-state to j-state and thermodynamic force, respectively. Fig. 12 shows EPR, the Fourier amplitude spectrum, W ij contour plot, and t − t diagram. These are calculated from the binarized data of the velocity deviation in Fig. 2 (d1) -2 (d6).
The W ij contour of Fig. 12 (b3) with the yellow area in the lower left is different from the W ij contours of Fig. 12 (b2), (b4) and (b5) which have the yellow area in the upper right, but the EPRs in Fig. 12 (b2)-(b5) are all greater than 0.08. Therefore, the W ij contour does not correspond to the EPR.
On the other hand, the large EPR, and the gentle negative curvature of the Fourier spectrum, which implies a small positive α-tremor correspond to the W ij contour with the yellow area in the lower left ( Fig. 12 (a3) and (b3)). The large EPR, and the sharp negative curvature, which implies a large positive α-tremor correspond to the W ij contour with the yellow area in the upper right ( Fig. 12 (a2), (b2), (a4), (b4), (a5) and (b5)). Consequently, the W ij corresponds to the pair of (EPR, α-tremor).
Macroscopic thermodynamics in the GEJE process Fig. 13 (a), which shows the thermodynamics of the UD ground vibration recorded at KSN from 2006-01-01 00:00 to 2016-12-31 23:59, compares the macroscopic EPR calculated every 10 days (EPR(10days)) with the macroscopic trend of positive α-tremor calculated every 51.2 seconds. The logarithm is applied to α-tremor and the positive component is extracted. A sharp decrease in EPR(10days) is seen at the timing of the strong peak of the positive α-tremor ( Fig. 13 (a)). Since EPR(10days) = 0 corresponds to a thermodynamic equilibrium or time independent system, the sharp decrease in EPR(10days) implies a abrupt change in the system from non-equilibrium toward equilibrium. The sharp decrease in EPR(10days) is observed at the timing of GEJE and in the period of approximately 3 years before and after GEJE.
EPR(10days) takes high, medium, and low values at the time intervals (b), (c), and (d) in Fig. 13(a). The characteristics of these time intervals are shown in Fig. 13 (b), (c), and (d). In all cases, the density of states uctuates around a value of 0.5. Fig. 13 (b) and (c), which are the cases of large and medium EPR(10days), both show that the frequency of occurrence is concentrated around α-tremor= -1 and EPR(10days)= 0.3 (left graph). In the positive α-tremor region, the frequency of occurrence is very low (center graph). On the other hand, for a small EPR(10days), the center graph in Fig. 13 (d) shows a recognizable concentration of occurrence frequency in the positive α-tremor region (red down arrow).

Equivalence Of Ground Vibrations And Ca
The binary sequence generated by periodic (d, p)-CA184 has a form similar to the form of the binarized velocity deviation and is considered to follow the thermodynamic equations Eq. (1) and (2). The thermodynamic parameters are calculated for the CA and compared with the parameters of the recorded data.

Cellular automaton
The CA-184 or Rule 184 cellular automaton is brie y reviewed, and (d, p)-CA184 is de ned as a stochastic extension of CA-184. Then, the data generation protocol using (d, p)-CA184 with periodic boundary conditions is discussed.
CA-184 (Nishinari & Takahashi, 1999) CA-184 which is frequently introduced for modeling tra c jams is the one-dimensional array of cells, each containing "1" or "0". The time evolution of the cells, which depends only on the cells on both side of it, follows the rule-184: where t and j are natural numbers, indicating time and space, respectively. The denominator U t + 1 j , which is the cell value at (time, space) = (t + 1,\,x j ), is determined by the numerator U t j − 1 U t j U t j + 1 , which is a sequence of three binary numbers at time t. The rule-184 states that when there is the sequence "10" (one followed by zero) at time t, it will become "01" (zero followed by one) at time t + 1. In other words, assuming "0" and "1" denote empty cell and the cell occupied by "1", respectively, if there is an empty cell on the right side of the occupied cell at time t, the "1" in the occupied cell moves right and occupies the empty cell at time t + 1. Therefore, the occupied cells move to the right as time passes. There are eight possible time evolution patterns as shown on the right side of Eq. (3). The binary number "10111000" obtained by arranging the denominators in order is 184 in decimal notation, which is the origin of the name rule-184.
(d, p)-CA184 : CA with probabilistic uctuations Probabilistic uctuation needs to be introduced in CA for expressing the complex phenomena such as ground vibration on the earth. Introducing the moving or hop probability p, the rule-184 (Eq. (3)) can be rephrased that if there is a sequence "10" in the numerator at time t, then with probability p = 1, the sequence become "01" at time t + 1. Extending the probability p to an arbitrary real number, p-CA184 is de ned as the following: if there is a sequence "10" in the numerator at time t, then with probability p, the sequence become "01" at time t + 1. In other words, when p-CA184 advances, it advances with probability p. The CA-184 with the moving probability p is called Asymmetric Simple Exclusion Process (ASEP). However, ASEP is called p-CA184 in this paper since p represents the state characteristics in the phase transition process described in later sections. In addition, since the initial density of states d plays an essential role in the phase transition, we de ne p-CA184 calculated with d as (d, p)-CA184. Here, the density of states is the ratio of the number of "1" to the total number of cells in the cell set under consideration.
Periodic (d, p)-CA184 and data generation The (d, p)-CA184 introduced in this paper consists of 10 cells in space, (x i , 1 ≤ i ≤ 10)) and has periodic boundary conditions (Fig. 14 (a)-(b)). At the boundaries, x 10 and x 1 , the cell value at time t + 1 is determined by referring the cells on both sides at time t. The boundary conditions are determined by applying the rule-184 to U t 10 U t 1 U t 2 ) /U t + 1 1 and U t 9 U t 10 U t 1 ) /U t + 1 10 ( Fig. 14 (b)). The CA data is generated by reading the cell value of x 10 at each time step t j (Fig. 14 (c)).

Comparison of recorded data with CA
The binary sequence generated by the periodic (d, p)-CA184 is considered to follow the thermodynamic equations Eq. (1) and Eq. (2). As in the case of the binarized velocity deviation, the thermodynamic parameters are calculated for the CA and compared directly with the thermodynamic parameters of the recorded data.
Recorded data and CA data are compared for the three 10-day intervals with large, medium and low EPR(10days) values shown in Fig. 13 (a). In order to compare dynamic parameters directly, the scale for calculating those parameters is uni ed. In this section, for both recorded and CA data, the EPR is calculated for the rst 1000 data points of the data block for which α-tremor is calculated. Here, one data block for either the recorded data or the CA data is de ned as 1024 data points equidistant in chronological order. For the recorded data, the one data block corresponds to 51.2 seconds of a data acquisition time. This section compares the 16875 data blocks of the recorded data (corresponding to the data length of 10 days) with the CA data blocks and evaluate the degree of agreement between them.
Eq. (2), the de nition of EPR, shows that EPR is calculated from transition rate matrix W ij and state vector p i . In other words, we need to know the EPR and W ij to determine p i . On the other hand, referring to the discussion for Fig. 12, W ij corresponds to a pair of EPR and α-tremor. So, we replace W ij with αtremor and EPR in evaluating the p i matching. Adopting α-tremor is advantageous since α-tremor is observable, is directly related to ground vibrations, and provides insights on earthquakes. Consequently, the EPR, α-tremor, p i , and additionally the density of states are evaluated in the comparison of the recorded with the CA data.
CA data characteristics CA data of periodic (d, p)-CA184 is generated by the method shown in Fig. 14. The initial density of states "d" is selected to be 0.5 to be coherent with the density of states of the recorded ground vibration data which uctuates around 0.5 as shown in Fig. 13 (b)-(d). 2000 CA data is generated for each hop probability "p" of 0.2, 0.3, 0.5, and 0.87, and the characteristics of the generated data are shown in Fig.  15. Fig. 15 shows frequency of (α-tremor, EPR) pairs, frequency of α-tremor, and density of states. In all gures in Fig. 15, the frequency of occurrence of α-tremor is concentrated around -1, with small values in the positive region, similar to the characteristics of the recorded ground vibration in Fig. 13.
If "p" is greater than or equal to 0.5, a larger positive α-tremor is observed (red down arrow in the α-tremor frequency plot in Fig. 15 (c) and (d)) than if "p" is less than 0.5. Therefore, the large "p" which is the nearequilibrium condition corresponds to more noticeable positive α-tremor distribution. On the other hand, for the recorded data, the noticeable positive α-tremor distribution is also seen in Fig. 13 (d) for the other near-equilibrium condition, the small EPR(10days). Consequently, the α-tremor frequency distribution in near-equilibrium systems is similar in both the 10-day scale recorded data and the 2000-data-point scale (d = 0.5, p)-CA184 data.
10-day comparison of recorded data and CA data Fig. 16 compares periodic (d = 0.5, p < 1)-CA184 with ground vibration data recorded for 10 days. The comparison is performed on 16875 recorded data blocks, each data block corresponding to 51.2 seconds. An evaluation vector containing 8 components: EPR, α-tremor, density of states, p 1 , p 2 , p 3 , p 4 , and p 5 is calculated for each data block of the CA and the recorded data. Then, the Pearson correlation between the evaluation vector of the recorded data block and the evaluation vector of the CA data block is calculated. For each evaluation vector of the recorded data block, Pearson correlations are calculated for all CA data blocks. The CA data block with maximum absolute value of the Pearson correlation is selected as the "optimal match", and compared with the recorded data block as shown in Fig. 16. In Fig.  16, the p i components are differentiated by colors so that the values of p i can roughly be read. Fig. 16 (b) shows the comparison regarding the 10-day time interval with the medium EPR(10days) value of 0.054. The CA data consists of 10000 data blocks generated for (d, p) of (0.5, 0.2), (0.5, 0.22), (0.5, 0.3), (0.5, 0.5), and (0.5, 0.87). For every (d, p), 2000 data blocks are generated. EPR, α-tremor, density of states, p 1 , p 2 , p 3 , p 4 , and p 5 simultaneously match between the recorded and CA data: 99.9% of the 16875 pairs of CA and recorded data blocks show a Pearson correlation greater than 0.9. Fig. 16 (a) shows the comparison regarding the 10-day time interval with the high EPR(10days) value of 0.06. To handle the high value in EPR, 4000 data blocks are generated for each of p = 0.2, 0.22, 0.3, 0.5 and 0.87, and 2000 data blocks are additionally generated for each of p = 0.15 and 0.17. Total 24000 data blocks are included in the CA data. 99.7% of the 16875 pairs of CA and recorded data blocks show a Pearson correlation greater than 0.9. In the all cases of the 10-day comparison between the recorded and the CA data, EPR, α-tremor, density of states, p 1 , p 2 , p 3 , p 4 , and p 5 are simultaneously in good agreement. For a given recorded data block, we can nd a CA data block that closely matches the recorded data block. It is implied that the recorded data block is a subset of the periodic (d = 0.5, p < 1)-CA184 data block. It is also implied that by appropriately selecting the pop probability "p", the ground vibration states can be reproduced with the periodic (d = 0.5, p < 1)-CA184, and the ground vibration in the GEJE process is equivalent to the periodic (d = 0.5, p < 1)-CA184.
It is worth noting that the match between the CA and the recorded data for the small EPR(10days) in Fig.   16 (c) is achieved by introducing the additional CA data generated for (d, p)=(0.5, 0.7). This implies that the small EPR(10days) system is closer to the equilibrium system of p=1 than the medium or large EPR(10days) system. This is consistent with the discussion for Fig. 13 that the smaller the EPR(10days), the closer to equilibrium. Therefore, both the CA and the recorded data provide similar thermodynamic information in a 10-day scale.
It should also be noted that if the evaluation vector of the recorded data block and the CA data block match, good matches are obtained for detailed dynamic properties such as the Fourier amplitude spectrum, the contour plot of W ij and the t − t diagram. With respect to the case of large EPR(10days), With respect to the case of low EPR(10days), Fig. 17 (b1) and (b2) compare the detailed characteristics of the best matched CA data block and the corresponding recorded data block, respectively. Both show large frequency components near 10 Hz in the Fourier spectrums, a large yellow distribution in the upper right of the transition rate matrix map, and a sparse distribution of blue cells in the t − t diagram.
The above comparisons show that there is coherences between the best matched CA and the corresponding recorded data with respect to the Fourier spectrum, the transition rate matrix, and the t − t diagrams of the binary sequence.

Conclusions
The weak ground vibrations in the GEJE process are characterized by α-tremor. Then, the invariance of αtremor is shown in the transformation from velocity to binarized velocity deviation. The binarized velocity deviation and the raw velocity signal are considered as equivalent as long as the ground vibration is considered as a α-tremor uctuation. Subsequently, vibrational states are de ned by the binarized velocity deviation, and the stochastic dynamics of transition of the state are considered to be described by the master equation. The binary sequence generated by the periodic (d, p)-CA184 has a form similar to the form of the binarized velocity deviation and is considered to follow Eq. (1) and (2). So, the thermodynamic parameters can be calculated for the CA and compared with the parameters of the recorded data. The binarized data recorded in the GEJE process is directly compared with the CA of periodic (d = 0.5, p < 1) -CA184.
The EPR, α-tremor, density of states, and state vector containing 5 components, calculated for each of the 16875 recorded data blocks are simultaneously in good agreement with those calculated for the (d = 0.5, p < 1)-CA184 data blocks. Over 99.7% of the 16875 pairs of recorded data blocks and CA data blocks show a Pearson correlation of 0.9 or greater. Therefore, the ground vibration in the GEJE process can be reproduced by periodic (d = 0.5, p < 1)-CA184 by carefully selecting the density d and the state condensation factor p, and the ground vibration is equivalent to the periodic (d = 0.5, p < 1)-CA184.
It is con rmed that if the CA data and the recorded data match for EPR, α-tremor, density of states, and state vector containing 5 components, then the CA data and the recorded data match also for the Fourier spectrum, transition rate matrix, and t − t diagram.
A direct link between the seismic process of a megathrust earthquake and a cellular automaton is established, which will provide CA as an alternative method to understand the seismic process.

Figure 11
Transition rate matrix and probability density of state. (a) The rst 100 data of the binarized velocity deviation data of Fig. 9 (b1). (b) t-t diagram, or pile of blocks containing 10 data points divided from (a).
Page 21/25 (c) The number of clusters in the block, or the index "i" of the vibration state s i . Chronological transition sequence from the lower state s i to the upper state s j . (d) The procedure for calculating W ij , which is a component of the transition rate matrix. (e) W ij calculated for (c). (f) Probability density distribution of states in the 100 data in (c). The numbers are not normalized. p i is the total number of the states s i in the 100 data.

Figure 17
Comparison of detailed parameters between the recorded data and its best-matched CA data. The left is the Fourier amplitude spectrum, the center is the lled contour of transition rate matrix, and the right is the t-t diagram, or the rst 100 data points of the binary sequences. (a1) Detailed parameters for the recorded data block of 51.2 seconds starting at 00:34:33 on 2009-04-30. The elapsed time is included in the time interval with large value of EPR(10days). (a2) Detailed parameters of the CA data block best-matched with the data block of (a1) that is generated by (d=0.5, p=0.5)-CA184. (b1) Detailed parameters for the recorded data block of 51.2 seconds starting at 03:38:01 on 2011-12-24. The elapsed time is included in the time interval with low value of EPR(10days). (b2) Detailed parameters of the CA data block bestmatched with the data block of (b1) that is generated by (d=0.5, p=0.87)-CA184.