Toy model of accumulation and radiation of strain energy driven by big data of earthquakes

We propose the pointwise toy-model, which demonstrates the accumulation and radiation of the (seismic) strain-energy. The calculated absorbed-seismic-energy is mapped to the accumulated strain-energy to balance the radiated strain-energy, because the accumulation of the strain-energy is caused by the subduction of one plate below another plate rather than the absorption of the seismic energy radiated from the hypocenter, which are located at other points. The distribution function of the accumulated strain-energy can be approximated with the Γ-distribution, where the high accumulated-strain-energy-tail is approximated by the inverse-power-law like tail. Meanwhile, the distributions of the radiated and accumulated strain-energies deviate from the Γ-distribution in the both low and high radiated strain-energy-domain. In particular, they have the inverse-power-law (Pareto) type tails in both high radiated and accumulated strain-energy-tails. The accumulated and radiated strain-energies calculated in the long period indicate that the present linear-map must be improved in order to balance the accumulated strain-energy with the radiated strain-energy.


Introduction
In the long history of the seismology [1], the studies of earthquakes have been developed by investigations of the dynamics of the plate, exclusively, because the earthquakes occur at the boundary of the plate as a result of the plate tectonics [2]. Then, physical quantities such as the fault-slip-distance and seismic-wave (P-wave, S-wave, coda wave, Rayleigh wave etc.) [3,4] are significant to characterize the dynamics of earthquakes. Meanwhile, one conclusion we have reached to is that the earthquake is categorized as complex system, which is beyond mathematical modelings and axioms. In particular, the platewise dynamics [5], that causes the earthquakes, make its mathematical modeling more difficult than the pointwise dynamics. Provided that such platewise dynamics attribute to the chaotic dynamics of earthquakes [6], we cannot predict the space-time dynamics of earthquakes, accurately. Meanwhile, recent developments of analyses of the big-data represented by the artificial intelligence (AI) are expected to give us new insights on such complex dynamics of the earthquakes, whose mathematical modelings involve difficulties. Then, we can expect that AI-learning obtained using big data-set of earthquakes might be also useful in order to predict the spacetime dynamics of earthquakes. Such our expectations, however, are readily blown up, when we remind that earthquakes are chaotic aging-system [7]. Indeed, the AI-learning postulates that the recurrences of events are learned, repeatedly, so that the AI-learning of the earthquakes is impossible, unless the recurrences of events (i.e, a series of earthquakes from the main shock to after shocks) are observed as chronological datum, repeatedly. In other words, the similar aging system must be observed in earthquakes with the different main shock, repeatedly, in order to enable AI-learning to predict spacetime evolutions of earthquakes. Such recurrences of the aging process in the earthquakes, however, have not been reported, yet. Consequently, the pointwise toy-model is still significant to avoid mathematical modeling of the complex dynamics of the platewise earthquakes. Provided that the pointwise toy-model with minimum set of parameters gives us useful insights on the dynamics of the earthquakes, the consideration of the pointwise toy-model is significant for understanding of the characteristics of the earthquakes in itself. Based on such an attitude, the epidemic-type aftershock sequence (ETAS) model by Ogata [8], Helmster and Sornette [9], spring-mass model by Abe-Kato [10] or temporally correlated network model by Abe and Suzuki [7,11] were proposed. In this paper, we discuss the new ponitwise toy-model, which demonstrates the statistical characteristics of the accumulated and radiated strain-energy [12,13].The time evolution of the strain-energy, which changes in accordance with the absorption [14] and radiation of the seismic energy [15] derived from occurrence of the earthquake, is formulated by mapping the absorbed seismic-energy to the accumulated strainenergy. Here, we assume that the absorption and radiation are evaluated by the function which depends on the inverse-power-law (IPL) of the distance between the location of hypocenter and focused location. Of course, the accumulation of the strain-energy by the only occurrences of earthquakes at other points is physically implausible, because the accumulation of the strain-energy is caused by the subduction of one plate below another plate. Then, the absorbed seismic-energy tends to be underestimated using our pointwise toy-model, so that the mapping of the absorbed seismic-energy to the accumulated strain-energy is essential to balance the accumulated strain-energy with the radiated strain-energy. Hereafter, we call our proposing pointwise toymodel as the accumulated and radiation of the strain-energy (ARE) model for convenience. The physical space is expressed with set of spheres with the radius Î  r 3 , whose center is located on the surface of the Earth. Here, physical domain of the Earth ( Î Ì   x 3 3 ) is bounded by Ω 2 (Ω 2 : surface of the Earth). Then, the neighborhood of the point ≔ ( x y , , 0 o and radius r), as shown in left half of figure 1. Here, the time evolution of the strainenergy, whose location satisfies [km] (z c : focal depth), is considered, then, those of the deep-focusearthquakes satisfying z c > 100 km are neglected, as shown in figure 2, whereas absorptions of the seismic energy derived from all the earthquakes during specific time-interval are considered. The absorbed seismic-energy is calculated, when the location of the hypocenter at t (Î +  : time) is outside ( ) x S r, o , whereas the radiated strainenergy is calculated, when the location of the hypocenter The absorbed seismic-energy is accumulated until the next earthquake occurs inside ( ) x S r, o . Now, such an accumulated absorbed-energy is called as the accumulated seismic-energy. Finally, the accumulated seismicenergy is mapped to the accumulated strain-energy.
Since map of the accumulated seismic-energy to the accumulated strain-energy has not been known yet, then, it is determined from the balance between the accumulated and radiated strain-energies using observed data-set of hypocenter. As a result, the ARE model is data-driven, partially. The seismic energy radiated from the hypocenter (x c ) is transferred to other spheres ( , 0, 40 (unit: • E, • N, km) and r = 100 [km], as shown in the left half of figure 1. In our study, the absorption of the seismic energy is, however, calculated using all the earthquakes, which occurred in the whole earth and are listed in the catalog by Japan Meteorological Agency. The radiation or absorption of the seismic energy in the upper domain ( ) x S r, o satisfying z > 0 seems to be implausible at glance, because there is no ground in z > 0. We, however, consider the re-radiation or reabsorption via the reflection of the radiated and absorbed seismic-energies occur on the ground. In short, the calculation of the radiation and absorption of the seismic energy in the upper domain of ( ) x S r, o (z > 0) are performed in a same manner with that in the lower domain of ( ) x S r, o (z 0). The undetermined map and parameters in the ARE model will be determined by big data-set of hypocenter, namely, data-set of locations of

ARE model
The seismic energy is radiated from the hypocenter with the energy that depends on its magnitude. Gutenberg [17] defined the relationship between the magnitude and radiated seismic-energy as 10 where E and M are the radiated seismic(strain)-energy from the hypocenter and magnitude of the earthquake, respectively. Then, E is the observable quantity, because M is observed. Other formulation of the radiated seismic-energy from the hypocenter on the basis of the slip of the plate is proposed by Fukuyama [18]. Here, the seismic energy, which is absorbed or radiated by the sphere ( ) x S r, o , depends on the positional relationship between x c (location of the hypocenter) and x o (center of the sphere ( ) x S r, o ). The right half frame of figure 1 shows that the seismic energy is absorbed inside˜( ) x S r, o (˜Î  S 2 : surface of S) when r l and [ ] j j Î 0, max and radiated inside S when l < r, because we assume that the seismic-wave [19,20] propagates toward the isotropic direction [3], whereas the heterogeneous ground of the Earth yields the anisotropic propagation of the seismic wave [21]. Additionally, we neglect the seismic energy-flux shaded by neighboring spheres and its spatial damping, which tend to decrease the seismic energy-flux.
Here, we assume that the seismic energy radiated from the hypocenter decays in accordance with the IPL of the distance ( ≔ | | x x ) measured from the hypocenter such as: is the IPL number related to the spatially decaying rate of E, ℓ ∞ is the scale-factor for non- . In later numerical analyses, ℓ ∞ = r = 100 [km] is used. From equation (2), the seismic-energy-flux, which is integrated in unit time, is calculated as where e ≔ x/|x|. The radiated seismic-energy from the hypocenter behaves as the point-wise energy-source such as the light source. Therefore, the radiated seismic-energy is absorbed in the irradiated part on S ( [ ] j j Î 0, max ), as shown in the right half frame of of figure 1. Then, the ARE model defines the absorbed and radiated seismic-energies (ε > 0: absorption, ε < 0: radiation) as: satisfiesˆ·ˆ= n x 0 (ˆÎ  n 2 ), as shown in the right half frame of figure 1. Quantities with superscriptˆare normalized by ℓ ∞ . Figure 2 shows ε versus l/r in cases of q = 1, 2 and 3. The absolute value of the radiated seismic-energy decreases in accordance with the increase of l/r in the case of q = 1 (i.e., q < 2). The absolute value of the radiated seismic-energy is constant in regardless of l/r in the case of q = 2. Then, the inverse power law with q = 2 satisfies the conservation of the seismic energy radiated from the hypocenter (E) inside S, when l/r < 1 is satisfied. The absolute value of the radiated seismic-energy increases, as the distance from the hypocenter and x o increases. Then, the radiated seismic-energy tends to be overestimated in the range of 2 < q. On the other hand, the absorbed seismic-energy increases from zero, as l/r increases from l/r = 1, in the case of q = 1, whereas the absorbed seismic-energy decreases in the range of 1 < l/r when 1 < q, as shown in cases of q = 2 and 3 in figure 2. Hereafter, the radiated seismic-energy is interpreted as the radiated strain-energy, whereas the absorbed seismic-energy is stored as a part of the accumulated strain-energy. The effects of spacial decaying-rate (q) on the distribution function of the accumulated and radiated strain-energies are discussed in section 3.
Reminding that the strain-energy is accumulated by the subduction of the plate, spontaneously, rather than successive absorptions of the radiated strain-energies, the accumulated seismic-energy (Ẽ a ) in S, which is the summation of the absorbed seismic-energy by the next radiation of the seismic energy in S, is presumably smaller than the accumulated strain-energy (E a ) in S, because the subduction of the plate is primary contribution to the accumulated strain-energy. Provided that there is the correlation betweenẼ a and E a , we can consider the map˜ E E a a . In section 3.1, we consider the numerical procedure to calculate E a via mappingẼ a .
3. Statistical analyses of data-set of accumulated and radiated seismic energies 3.1. Preparation of data-set of accumulated and radiated strain-energies In section 2, the ARE model formulated the absorbed and radiated seismic-energies, as shown in equation (4).
The catalog of the earthquakes collected by Japan Meteorological Agency stores big-data-set of earthquakes in Japan including information of hypocenter such as Î  x c 3 (location of hypocenter), t k (origin-time of k-th earthquake), M (magnitude). Using these data-set, we can calculate ( )  x t, o at each origin-time (t ä U k t k ). In order to obtain data-set of accumulated and radiated strain-energies, namely, ( ( , , a k o r k o , the following procedures 1-3 are executed:  (1) and (4) 2.Ẽ a and E r are calculated using the following rule: where dt is Lebesgue's measure, ( ) d x is the Kronecker's delta function, T is set as T = τ 0 = 1 day, when On the other hand, T is set as Here, τ m is the waiting time until the next earthquake occurs inside ( ) x S o , when the earthquake occurs inside ( )

The accumulated seismic-energy, namely,˜(
o is shown in figure 3. The concrete form of map is discussed later using data-set of hypocenter. E a (Ẽ a ), E r and E s ≔ E a − E r (strain-energy) are set as zero, when t k = T m .

Statistical result of (E a ,E r )
The statistical characteristics of ( (     (8)). The parameters defining f Γ are shown in Tabs I and II. Dashed lines in f (E a ) versus E a and f (E r ) versus E r correspond to the Pareto distributions (˜( , whose parameters ϖ(q) and ñ(q) are defined in table 4.
, whose parameters ς(q) are defined in table 4.

The distribution function of
where N is the total number of data-set of ( ) E E , f E a versusẼ a and ( ) f E r versus E r are fitted by Γ-function ( f Γ in equation (8)) and Pareto type distribution, whose parameters are defined in tables 1, 2 and 4. Figure 5 shows that E a = E r is satisfied in most part of (˜) E E , a r , as predicted by the fact that the accumulated strain-energy is usually larger than the accumulated seismic-energy. (˜) f E a is well fitted by the Γ-function in the lowẼ a domain (i.e.,˜< E 10 a 11 ), whereas (˜) f E a follows the Pareto type distribution in the highẼ a -tail (i.e., 10 11 < E a ). Meanwhile, ( ) f E r is not fitted by the Γ-function in the low E r (i.e., E r < 10 14 ), whereas ( ) f E r follows the Pareto type distribution at high E r tail (i.e., 10 14 < E r ). The highẼ a and E r tails of (˜) f E a and ( ) f E r become fatter, as q increases, because e a and |e r | increases around l/r = 1, as q increases, as shown in figure 2. Finally, we can confirm that linear map ofẼ a with inclination 10 3 fits E r with the best accuracy amongẼ 10 m a (m = − 6, −3, 0, 3, 6). Then, we use map of Table 1. Parameters u, v and x 1 which define Γ-function in equation (8) fitting to (˜) f E a in figure 5, in cases of q = 1, 1.5, 2, 2.5 and 3.  Table 2. Parameters u, v and x 1 which define Γ-function in equation (8) fitting to ( ) f E r in figure 5, in cases of q = 1, 1.5, 2, 2.5 and 3.  Table 3. Parameters u, v and x 1 which define Γ-function in equation (8) fitting to ( ) f E a in figure 6, in cases of q = 1, 1.5, 2, 2.5 and 3.  and ( ) f E r obtained using data-set before the Great East Japan Earthquake on 11th Mar. in 2011 and those obtained using data-set before the Great East Japan Earthquake on 11th Mar. in 2011are demonstrated in order to investigate the effects of the big earthquake in appendix B.  From results of ( ) f E a and ( ) f E r in figure 6, we conjecture that the dynamics of E a and E r might follow superdiffusion such asĹe vy-flight [22] in ( ) E E , a r -plane. Such an investigation of motion of point-( a r -plane is set as our future study.

Temporal variations of E s in prolonged interval
In previous discussions,˜( )  (5). x E r o on x o . In short, regimes with high E a and E r concentrate around the Japan-subduction-zone, Additionally, it is interesting that both E a and E r decrease, concentrically, around Japan-subduction-zone. In order to evaluate of differences between E r and E a , we define new parameter ρ such as becomes narrower in comparison of that when T m =8 years. Figure 9 showsẼ a * and E r * versus t/day (left frame) and E a * and E r * versus t/day (right frame) at three locations of x o , namely, (36°N, 140°E) (point A: around  .Ẽ a * and E r * versus t/day (left frame) and E a * and E r * versus t/day (right frame) in the case of q = 2.

Concluding Remarks
In this paper, we proposed the ARE model to understand the statistical characteristics of the accumulated and radiated strain-energies. Here, the absorbed and radiated seismic energies are modeled by the radiated seismicenergy which follows the IPL law of the distance from the hypocenter. The accumulated strain-energy is calculated by mapping the accumulated seismic-energy to balance with the radiated strain-energy. The distribution function of the accumulated strain-energy (E a ) can be approximated by not Γ-distribution but Pareto type distribution at the high E a tail, when the IPL parameter defining the energy-flux, q, is changed. On the other hand, the distribution function of the radiated strain-energy (E r ) deviates from Γ-distribution in small domain of E r , where the distribution function of the radiated strain-energy has the Pareto type distribution at high E r tail similarly to that of E a . The accumulated E a and E r during longer period indicated that the linear mapping ofẼ a to E a is too rough approximation to satisfy that E a ; E r . In order to increase the accuracy of mapping, additional parameters related to the locality (ground characteristics etc) might be useful.
We readily confirm that ( ) f E r versus E r approach Γ-distribution with the fat tail, as ò decreases. In short, deviations of ( ) f E r from Γ-distribution in the regime of low E r decrease, as ò increases. As a result, we can conclude that data-set ( ) E E , a r , whose total number at ( ) x y , o o is small, attributes deviation of ( ) f E r from Γdistribution. On the other hand, the high E a -tail of ( ) f E a decreases markedly, when ò = 2 × 10 6 .