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

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 strain-energy. 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 toy-model 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\in {{\mathbb{R}}}^{3}$, whose center is located on the surface of the Earth. Here, physical domain of the Earth (${\boldsymbol{x}}\in {{\mathbb{X}}}^{3}\subset {{\mathbb{R}}}^{3}$) is bounded by Ω² (Ω²: surface of the Earth). Then, the neighborhood of the point ${{\boldsymbol{x}}}_{o}:= \left({x}_{o},{y}_{o},{z}_{o}=0\right)$ is expressed with the spherical domain $S\left(r,{{\boldsymbol{x}}}_{o}\right)$ ($S\left(r,{{\boldsymbol{x}}}_{o}\right)$ has a center at ${{\boldsymbol{x}}}_{o}:= \left(x,y,0\right)$ and radius r), as shown in left half of figure 1. Here, the time evolution of the strain-energy, whose location satisfies ${z}_{c}\in \left[0,100\right]$ [km] (z_c : focal depth), is considered, then, those of the deep-focus-earthquakes 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 ($\in {{\mathbb{R}}}_{+}$: time) is outside $S\left(r,{{\boldsymbol{x}}}_{o}\right)$, whereas the radiated strain-energy is calculated, when the location of the hypocenter (${{\boldsymbol{x}}}_{c}\in {{\mathbb{X}}}^{3}$) at t is inside $S\left(r,{{\boldsymbol{x}}}_{o}\right)$ (i.e., ${{\boldsymbol{x}}}_{c}(t)\in S\left(r,{{\boldsymbol{x}}}_{o}\right)$. The absorbed seismic-energy is accumulated until the next earthquake occurs inside $S\left(r,{{\boldsymbol{x}}}_{o}\right)$. Now, such an accumulated absorbed-energy is called as the accumulated seismic-energy. Finally, the accumulated seismic-energy is mapped to the accumulated strain-energy.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 (${{\boldsymbol{x}}}_{c}\in {{\mathbb{X}}}^{3}-S\left(r,{{\boldsymbol{x}}}_{o}\right)$). x _o are equally spaced in X − Y plane by setting ${{\boldsymbol{x}}}_{o}=\left(100+i,10+j,0\right)$ $(i\in [0,50]\cap {{\mathbb{Z}}}_{+},\,j\in [0,40]\cap {{\mathbb{Z}}}_{+})$ (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 $S\left(r,{{\boldsymbol{x}}}_{o}\right)$ satisfying z > 0 seems to be implausible at glance, because there is no ground in z > 0. We, however, consider the re-radiation or re-absorption 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 $S\left(r,{{\boldsymbol{x}}}_{o}\right)$ (z > 0) are performed in a same manner with that in the lower domain of $S\left(r,{{\boldsymbol{x}}}_{o}\right)$ (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 hypocenter, origin-time, magnitudes of the earthquakes, which occurred in Japan from 1st Jan. 2000 to 31th Dec. 2016 (17 years period) [16]. The number of the total events, which are calculated in our study, is set as 2 639 486. Finally, the ARE model certainly gives us some insights on the statistical characteristics of the dynamics of the earthquakes in terms of the seismic energy.

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

$$\begin{eqnarray}&&{\mathrm{log}}_{10}E=4.8+1.5M,\end{eqnarray} \tag{ 1 }$$

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 $S\left(r,{{\boldsymbol{x}}}_{o}\right)$, depends on the positional relationship between x _c (location of the hypocenter) and x _o (center of the sphere $S\left(r,{{\boldsymbol{x}}}_{o}\right)$). The right half frame of figure 1 shows that the seismic energy is absorbed inside $\tilde{S}\left(r,{{\boldsymbol{x}}}_{o}\right)$ ($\tilde{S}\in {{\mathbb{R}}}^{2}$: surface of S) when r ≤ l and $\varphi \in \left[0,{\varphi }_{\max }\right]$ 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:= \left|{\boldsymbol{x}}\right|$) measured from the hypocenter such as:

$$\begin{eqnarray}&&f\left(x\right):= \displaystyle \frac{E}{4\pi {\left(x/{{\ell }}_{\infty }\right)}^{q}}=\displaystyle \frac{E}{4\pi {\hat{x}}^{q}},\end{eqnarray} \tag{ 2 }$$

where $q\in {{\mathbb{R}}}_{+}$ is the IPL number related to the spatially decaying rate of E, ℓ_∞ is the scale-factor for non-dimensionalization of x and $\hat{x}=| \hat{{\boldsymbol{x}}}| := | {\boldsymbol{x}}/{{\ell }}_{\infty }|$. 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

$$\begin{eqnarray}&&{\boldsymbol{f}}\left(\hat{{\boldsymbol{x}}}\right):= \displaystyle \frac{E}{4\pi | {\boldsymbol{x}}/{{\ell }}_{\infty }{| }^{q}}\displaystyle \frac{{\boldsymbol{x}}}{\left|{\boldsymbol{x}}\right|}=\displaystyle \frac{E}{4\pi {\hat{x}}^{q}}{\boldsymbol{e}}.\end{eqnarray} \tag{ 3 }$$

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 ($\varphi \in [0,{\varphi }_{\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:

$$\begin{eqnarray}\varepsilon =\left\{\begin{array}{ll}{\rm{(Absorption)}} & \\ {e}_{a}:= \left(2\pi \right){\displaystyle \int }_{0}^{{\varphi }_{\max }}\left(-{\boldsymbol{f}}\left(\hat{{\boldsymbol{x}}}\right)\cdot \hat{{\boldsymbol{n}}}\right){\hat{r}}^{2}\sin \varphi d\varphi =\left(2\pi \right){\displaystyle \int }_{0}^{{\varphi }_{\max }}-\left|{\boldsymbol{f}}\left({\boldsymbol{x}}\right)\right|\displaystyle \frac{\hat{l}\sin \varphi \cos \varphi -\hat{r}\sin \varphi }{\left|\hat{{\boldsymbol{x}}}\right|}{\hat{r}}^{2}d\varphi . & \left(1\lt l/r\right)\\ {\rm{(Radiation)}} & \\ {e}_{r}:= \left(2\pi \right){\displaystyle \int }_{0}^{\pi }\left(-{\boldsymbol{f}}\left(\hat{{\boldsymbol{x}}}\right)\cdot \hat{{\boldsymbol{n}}}\right){\hat{r}}^{2}\sin \varphi d\varphi =\left(2\pi \right){\displaystyle \int }_{0}^{\pi }-\left|{\boldsymbol{f}}\left(\hat{{\boldsymbol{x}}}\right)\right|\displaystyle \frac{\hat{l}\sin \varphi \cos \varphi -\hat{r}\sin \varphi }{\left|\hat{{\boldsymbol{x}}}\right|}{\hat{r}}^{2}d\varphi & \left(l/r\leqslant 1\right)\end{array}\right.\end{eqnarray} \tag{ 4 }$$

where $\hat{{\boldsymbol{n}}}:= {\boldsymbol{n}}/{{\ell }}_{\infty }$ is the unit normal vector on $\tilde{S}\left(\hat{r},{\hat{{\boldsymbol{x}}}}_{o}\right)$, and ${\varphi }_{\max }:= {\cos }^{-1}\left(r/l\right)$ satisfies $\hat{{\boldsymbol{n}}}\cdot \hat{{\boldsymbol{x}}}=0$ ($\hat{{\boldsymbol{n}}}\in {{\mathbb{R}}}^{2}$), as shown in the right half frame of figure 1. Quantities with superscript $\hat{}$ 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 (${\tilde{E}}_{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 ${\tilde{E}}_{a}$ and E_a , we can consider the map ${\tilde{E}}_{a}\to {E}_{a}$. In section 3.1, we consider the numerical procedure to calculate E_a via mapping ${\tilde{E}}_{a}$.

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 ${{\boldsymbol{x}}}_{c}\in {{\mathbb{X}}}^{3}$ (location of hypocenter), t_k (origin-time of k-th earthquake), M (magnitude). Using these data-set, we can calculate $\epsilon \left(t,{{\boldsymbol{x}}}_{o}\right)$ at each origin-time (t ∈ ⋃ _k t_k ).

In order to obtain data-set of accumulated and radiated strain-energies, namely, $\left({E}_{a}\left({t}_{k},{{\boldsymbol{x}}}_{o}\right),{E}_{r}\left({t}_{k},{{\boldsymbol{x}}}_{o}\right)\right)$, the following procedures 1-3 are executed:

• 1.  
  ${e}_{a}\left({t}_{k},{{\boldsymbol{x}}}_{o}\right)$ and ${e}_{r}\left({t}_{k},{{\boldsymbol{x}}}_{o}\right)$ are calculated using data set of x _c and M from equations (1) and (4)
• 2.  
  ${\tilde{E}}_{a}$ and E_r are calculated using the following rule:

  $$\begin{eqnarray}\begin{array}{rcl}{\tilde{E}}_{a}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right) & = & {\int }_{{t}_{m}({{\boldsymbol{x}}}_{o})}^{{t}_{m}({{\boldsymbol{x}}}_{o})+T}{dt}\delta \left(t-{t}_{k}\right){e}_{a}\left({t}_{k},{{\boldsymbol{x}}}_{o}\right),\\ {E}_{r}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right) & = & {\int }_{{t}_{m}({{\boldsymbol{x}}}_{o})}^{{t}_{m}({{\boldsymbol{x}}}_{o})+T}{dt}\delta \left(t-{t}_{k}\right)\left|{e}_{r}\left({t}_{k},{{\boldsymbol{x}}}_{o}\right)\right|,\quad m\in {{\mathbb{Z}}}_{+},\end{array}\end{eqnarray} \tag{ 5 }$$

  where dt is Lebesgue's measure, $\delta \left(x\right)$ is the Kronecker's delta function, T is set as T = τ₀ = 1 day, when $0\lt \left|{e}_{r}\left({t}_{k}\right)\right|$ during $t\in \left[{t}_{m}\left({{\boldsymbol{x}}}_{o}\right),{t}_{m}\left({{\boldsymbol{x}}}_{o}\right)+{\tau }_{0}\right]$. On the other hand, T is set as $T={\tau }_{m}\left({{\boldsymbol{x}}}_{o}\right)$, when ${e}_{r}\left({t}_{k}\right)=0$ during $t\in \left[{t}_{m}\left({{\boldsymbol{x}}}_{o}\right),{t}_{m}\left({{\boldsymbol{x}}}_{o}\right)+{\tau }_{m}\left({{\boldsymbol{x}}}_{o}\right)\right)$ and $0\lt \left|{e}_{r}\left({t}_{k}={t}_{m}\left({{\boldsymbol{x}}}_{o}\right)+{\tau }_{m}\left({{\boldsymbol{x}}}_{o}\right)\right)\right|$. Finally, ${t}_{m+1}\left({{\boldsymbol{x}}}_{o}\right)$ is set as ${\min }_{k}{t}_{k}| {t}_{k}\gt {\tau }_{m}+T$. Here, τ_m is the waiting time until the next earthquake occurs inside $S\left({{\boldsymbol{x}}}_{o}\right)$, when the earthquake occurs inside $S\left({{\boldsymbol{x}}}_{o}\right)$ at $t={t}_{m}\left({{\boldsymbol{x}}}_{o}\right)$. As a result, ${T}_{m}={t}_{m}\left({{\boldsymbol{x}}}_{o}\right)+T={t}_{m+1}\left({{\boldsymbol{x}}}_{o}\right)$.
• 3.  
  The accumulated seismic-energy, namely, ${\tilde{E}}_{a}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right)$ is mapped to the accumulated strain-energy ${E}_{a}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right)$ to balance ${E}_{a}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right)$ with the radiated strain-energy ${E}_{r}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right)$. The image of map of ${\tilde{E}}_{a}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right)\to {E}_{a}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right)$ is shown in figure 3. The concrete form of map is discussed later using data-set of hypocenter. E_a (${\tilde{E}}_{a}$), E_r and E_s := E_a − E_r (strain-energy) are set as zero, when t_k = T_m .

Zoom In Zoom Out Reset image size

3.2. Statistical result of (E_a ,E_r )

The statistical characteristics of $\left({E}_{a}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right),{E}_{r}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right)\right)$ are investigated on the basis of their data-set created in accordance with procedures 1-3 in section 3.1. The map ${\tilde{E}}_{a}\to {E}_{a}$ is determined from data-set of $\left({\tilde{E}}_{a}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right),{E}_{r}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right)\right)$. As discussed in Introduction, $\left({E}_{a}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right),{E}_{r}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right)\right)$ in the hypocenter within z = 100 [km] are considered. For reference, the locations of hypocenter in physical domain including Japan-subduction-zone are shown in figure 4. We can confirm that a large number of hypocenter also exist in the range of z > 100 [km], whereas we do not calculate $\left({\tilde{E}}_{a}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right),{E}_{r}\left({T}_{m},{{\boldsymbol{x}}}_{o}\right)\right)$ in the range of z > 100 [km], because the ratio of the event-number satisfying 100 < z [km] to the total number of events is small, adequately.

Zoom In Zoom Out Reset image size

The distribution function of $({\tilde{E}}_{a},{E}_{r})$, namely, $f\left({\tilde{E}}_{a},{E}_{r}\right)$ is defined by:

$$\begin{eqnarray}&&f\left({\tilde{E}}_{a},{E}_{r}\right):= \displaystyle \frac{1}{N}\sum _{{{\boldsymbol{x}}}_{o}\in {\mathbb{X}},{T}_{m}\in \left[0,{T}_{{\rm{fin}}}\right]}\delta \left({\tilde{E}}_{a}-{\tilde{E}}_{a}\left({{\boldsymbol{x}}}_{o},{T}_{m}\right)\right)\delta \left({E}_{r}-{\tilde{E}}_{r}\left({{\boldsymbol{x}}}_{o},{T}_{m}\right)\right),\end{eqnarray} \tag{ 6 }$$

where N is the total number of data-set of $\left({E}_{a},{E}_{r}\right)$, T_fin is the final time of t. Similarly, $f\left(\tilde{{E}_{a}}\right)$, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ are calculated by

$$\begin{eqnarray}\displaystyle \begin{array}{rcl}f\left({\tilde{E}}_{a}\right): & = & \displaystyle \frac{1}{N}\sum _{{{\boldsymbol{x}}}_{o}\in {\mathbb{X}},{T}_{m}\in \left[0,{T}_{{\rm{fin}}}\right]}\delta \left({\tilde{E}}_{a}-{\tilde{E}}_{a}({{\boldsymbol{x}}}_{o},{T}_{m})\right),\\ f\left({E}_{a}\right): & = & \displaystyle \frac{1}{N}\sum _{{{\boldsymbol{x}}}_{o}\in {\mathbb{X}},{T}_{m}\in \left[0,{T}_{{\rm{fin}}}\right]}\delta \left({E}_{a}-{E}_{a}({{\boldsymbol{x}}}_{o},{T}_{m})\right),\\ f\left({E}_{r}\right): & = & \displaystyle \frac{1}{N}\sum _{{{\boldsymbol{x}}}_{o}\in {\mathbb{X}},{T}_{m}\in \left[0,{T}_{{\rm{fin}}}\right]}\delta \left({E}_{r}-{E}_{r}({{\boldsymbol{x}}}_{o},{T}_{m})\right).\end{array}\end{eqnarray} \tag{ 7 }$$

From equations (6) and (7) $f\left({\tilde{E}}_{a},{E}_{r}\right)$, $f\left({\tilde{E}}_{a}\right)$, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ do not have temporal resolutions.

Next, $f\left(\tilde{{E}_{a}}\right)$, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ are fitted by Γ-distribution (f_Γ), which is defined by

$$\begin{eqnarray}&&{f}_{{\rm{\Gamma }}}\left(x;u,v\right):= \displaystyle \frac{{v}^{u}}{{\rm{\Gamma }}(u)}{\left(x-{x}_{1}\right)}^{u-1}\exp \left(-\displaystyle \frac{x-{x}_{1}}{v}\right),\end{eqnarray} \tag{ 8 }$$

Additionally, $f\left(\tilde{{E}_{a}}\right)$, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ are fitted by the Pareto type distribution such as $f\left(\tilde{{E}_{a}}\right)\propto {\tilde{E}}_{a}^{\varpi (q)}$, $f\left({E}_{a}\right)\propto {E}_{a}^{-\varsigma (q)}$ and $f\left({E}_{r}\right)\propto {E}_{r}^{-\varrho (q)}$.

Firstly, figure 5 shows $f\left({\tilde{E}}_{a},{E}_{r}\right)$ versus ${\tilde{E}}_{a}$ and E_r together with $f\left({\tilde{E}}_{a}\right)$ versus ${\tilde{E}}_{a}$ and $f\left({E}_{r}\right)$ versus E_r in cases of q = 1 (upper-left frame), 1.5 (upper-middle frame), 2 (upper-right frame), 2.5 (lower-left frame) and 3 (lower-right frame). Here, lines corresponding to linear recurrence lines ${E}_{r}={10}^{m}{\tilde{E}}_{a}$ (m = − 6, − 3, 0, 3 and 6) are added. Additionally, $f\left({\tilde{E}}_{a}\right)$ versus ${\tilde{E}}_{a}$ and $f\left({E}_{r}\right)$ 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 $\left({\tilde{E}}_{a},{E}_{r}\right)$, as predicted by the fact that the accumulated strain-energy is usually larger than the accumulated seismic-energy. $f\left({\tilde{E}}_{a}\right)$ is well fitted by the Γ-function in the low ${\tilde{E}}_{a}$ domain (i.e., ${\tilde{E}}_{a}\lt {10}^{11}$), whereas $f\left({\tilde{E}}_{a}\right)$ follows the Pareto type distribution in the high ${\tilde{E}}_{a}$-tail (i.e., 10¹¹ < E_a ). Meanwhile, $f\left({E}_{r}\right)$ is not fitted by the Γ-function in the low E_r (i.e., E_r < 10¹⁴), whereas $f\left({E}_{r}\right)$ follows the Pareto type distribution at high E_r tail (i.e., 10¹⁴ < E_r ). The high ${\tilde{E}}_{a}$ and E_r tails of $f\left({\tilde{E}}_{a}\right)$ and $f\left({E}_{r}\right)$ 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 ${\tilde{E}}_{a}$ with inclination 10³ fits E_r with the best accuracy among ${10}^{m}{\tilde{E}}_{a}$ (m = − 6, −3, 0, 3, 6). Then, we use map of ${E}_{a}={10}^{m(q)}{\tilde{E}}_{a}$ ($m\left(q=1\right)=2.2$, $m\left(q=1.5\right)=2.56$, $m\left(q=2\right)=2.85$, $m\left(q=2.5\right)=3.05$ and $m\left(q=3\right)=3.15$), where $m\left(q\right)$ is used as the best-fitting parameter, and plot $f\left({E}_{a},{E}_{r}\right)$ in figure 6. Figure 6 shows $f\left({E}_{a},{E}_{r}\right)$ versus E_a and E_r together with $f\left({E}_{a}\right)$ versus E_a and $f\left({E}_{r}\right)$ versus E_r in cases of q = 1 (upper-left frame), 1.5 (upper-middle frame), 2 (upper-right frame), 2.5 (lower-left frame) and 3 (lower-right frame). Linear recurrence-lines (E_r = 10^−m(q) E_a , m = −6, −3, 0, 3 and 6) are added to contours of $f\left({E}_{a},{E}_{r}\right)$. $f\left({E}_{a}\right)$ versus E_a are well fitted by Γ-function (f_Γ in equation (8)). The parameters defining f_Γ for $f\left({E}_{a}\right)$ are shown in table 3. Additionally, parameters ς, which defines the Pareto type distribution, are shown in table 4. The bulk tendency of $f\left({E}_{a}\right)$ versus E_a is similar to that of $f\left({\tilde{E}}_{a}\right)$ versus ${\tilde{E}}_{a}$. Meanwhile, the effects of cut-of-event-number at the specific domain ($S\left({x}_{o},{y}_{o}\right)$) on f(E_a , E_r ), $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ are investigated in appendix A. Finally, we can confirm ${\mathbb{E}}\left({E}_{a}\right)\sim {\mathbb{E}}\left({E}_{r}\right)$ (${\mathbb{E}}\left(x\right)$: expectation-value of x). Differences between $f\left({E}_{a},{E}_{r}\right)$, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1. Parameters u, v and x₁ which define Γ-function in equation (8) fitting to $f\left({\tilde{E}}_{a}\right)$ in figure 5, in cases of q = 1, 1.5, 2, 2.5 and 3.

q	1.00	1.50	2.00	2.50	3.00
u	32.1	14.7	10.6	10.4	10.8
v	0.159	0.230	0.277	0.292	0.309
x1	6.41	7.81	8.08	7.92	7.75

Table 2. Parameters u, v and x₁ which define Γ-function in equation (8) fitting to $f\left({E}_{r}\right)$ in figure 5, in cases of q = 1, 1.5, 2, 2.5 and 3.

q	1	1.5	2	2.5	3
u	116	107	84.9	56.3	44.9
v	0.110	0.114	0.127	0.153	0.169
x1	0.989	1.60	3.11	5.38	6.65

Table 3. Parameters u, v and x₁ which define Γ-function in equation (8) fitting to $f\left({E}_{a}\right)$ in figure 6, in cases of q = 1, 1.5, 2, 2.5 and 3.

q	1	1.5	2	2.5	3
u	31.9	14.6	10.5	10.7	10.7
v	0.159	0.232	0.277	0.288	0.309
x1	8.63	10.4	10.9	10.9	10.9

Table 4. Parameters ϖ, ϱ and ς which define the Pareto distribution in figures 5 and 6, in cases of q = 1, 1.5, 2, 2.5 and 3.

q	1	1.5	2	2.5	3
ϖ	0.64	0.638	0.59	0.552	0.504
ϱ	0.575	0.598	0.546	0.56	0.512
ς	0.673	0.614	0.587	0.552	0.493

From results of $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ in figure 6, we conjecture that the dynamics of E_a and E_r might follow super-diffusion such as ${\rm{L}}\acute{{\rm{e}}}$ vy-flight [22] in $\left({E}_{a},{E}_{r}\right)$-plane. Such an investigation of motion of point-$\left({E}_{a},{E}_{r}\right)$ in $\left({E}_{a},{E}_{r}\right)$-plane is set as our future study.

3.3. Temporal variations of E_s in prolonged interval

In previous discussions, ${\tilde{E}}_{a}\left({{\boldsymbol{x}}}_{o}\right),{E}_{a}\left({{\boldsymbol{x}}}_{o}\right),{E}_{r}\left({{\boldsymbol{x}}}_{o}\right),{E}_{s}\left({{\boldsymbol{x}}}_{o}\right)$ are set as zero, when t = T_m , because of equation (5). Now, we define ${\tilde{E}}_{a}^{* }\left({{\boldsymbol{x}}}_{o},{T}_{m}\right),{E}_{a}^{* }\left({{\boldsymbol{x}}}_{o},{T}_{m}\right),{E}_{r}^{* }\left({{\boldsymbol{x}}}_{o},{T}_{m}\right)$, and ${E}_{s}^{* }\left({{\boldsymbol{x}}}_{o},{T}_{m}\right)$ as ${\tilde{E}}_{a}^{* }\left({{\boldsymbol{x}}}_{o},{T}_{m}\right)={\sum }_{k\,=\,1}^{m}{\tilde{E}}_{a}\left({{\boldsymbol{x}}}_{o},{T}_{k}\right)$, ${E}_{a}^{* }\left({T}_{m}\right)={\sum }_{k\,=\,1}^{m}{E}_{a}\left({{\boldsymbol{x}}}_{o},{T}_{k}\right)$, ${E}_{r}^{* }\left({{\boldsymbol{x}}}_{o},{T}_{m}\right)={\sum }_{k}{E}_{r}\left({{\boldsymbol{x}}}_{o},{T}_{k}\right)$ and ${E}_{s}^{* }\left({{\boldsymbol{x}}}_{o},{T}_{m}\right)={\sum }_{k\,=\,1}^{m}{E}_{r}\left({{\boldsymbol{x}}}_{o},{T}_{k}\right)+{E}_{a}\left({{\boldsymbol{x}}}_{o},{T}_{k}\right)$ without setting them as zero at t = T_m . Here, E_a is mapped by ${E}_{a}={10}^{m(q)}{\tilde{E}}_{a}$, as discussed in section 3.2. Figure 7 shows contours of ${\tilde{E}}_{a}^{* }\left({T}_{m}\right)$ and ${E}_{r}^{* }\left({T}_{m}\right)$ when T_m = 8 years (left-half frames) and 16 years (right-half frames), in the case of q = 2. We can confirm the dependencies of both ${E}_{a}\left({{\boldsymbol{x}}}_{o}\right)$ and ${E}_{r}\left({{\boldsymbol{x}}}_{o}\right)$ 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

$$\begin{eqnarray}&&\rho \left({T}_{m}\right):= \displaystyle \frac{{E}_{a}^{* }\left({T}_{m}\right)-{E}_{r}^{* }\left({T}_{m}\right)}{{E}_{r}^{* }\left({T}_{m}\right)}.\end{eqnarray} \tag{ 9 }$$

figure 8 shows contours of ρ, when T_m = 8 years (left frame) and 16 years (right frame). ρ is smaller than unity around Japan-subduction-zone when T_m = 8 years, whereas regime satisfying $\rho \in {{\mathbb{R}}}_{-}$ when T_m = 16 years, becomes narrower in comparison of that when T_m =8 years. Figure 9 shows ${\tilde{E}}_{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 Tokyo), (35°N, 136°E) (point B: around Osaka city), (34°N, 130°E) (point C: around Fukuoka city). ${E}_{r}^{* }$ are markedly larger than ${\tilde{E}}_{a}^{* }$ in cases of points A-C, whereas ${E}_{r}^{* }\gt {E}_{a}^{* }$ at point A, ${E}_{r}^{* }\gt {E}_{a}^{* }$ at point B and E_a > E_r (4000 < t/day) and E_a < E_r (t ≤ 4000/day) at point C are obtained. Therefore, linear mapping of ${\tilde{E}}_{a}\to {E}_{a}$ is too rough to satisfy E_a ≃ E_r or E_r ≤ E_a , because the radiation occurs after E_a increases by the critical value. Another map $\psi ({\tilde{E}}_{a},{{\boldsymbol{x}}}_{o})\to {E}_{a}({{\boldsymbol{x}}}_{o})$ must be necessary by adding the variable of recurrence such as x _o together with ${\tilde{E}}_{a}$ in order to decrease the distance $\left|{E}_{a}-{E}_{r}\right|$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 seismic-energy 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 ${\tilde{E}}_{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 discussed $f\left({E}_{a},{E}_{r}\right)$, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ in section 3.2. Data-set of $\left({E}_{a},{E}_{r}\right)$ is composed by following procedures 1-3 in section 3.1. Here, we investigate the effects of event-number on $f\left({E}_{a},{E}_{r}\right)$, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$. In order to accomplish our aim, the lower-limit of number of data-set is considered. Setting ${N}_{t}\left({x}_{o},{y}_{o}\right):= \left|{\bigcup }_{m}\left({E}_{a}\left({T}_{m},{x}_{o},{y}_{o}\right),{E}_{r}\left({T}_{m},{x}_{o},{y}_{o}\right)\right)\right|$ as the total number of data-set at $S\left({x}_{o},{y}_{o}\right)$, we calculate $f\left({E}_{a},{E}_{r}\right)$, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ using data-set at $S\left({x}_{o},{y}_{o}\right)$, which satisfies ≤ N_t ($\epsilon \in {{\mathbb{R}}}_{+}$:cut-off-event-number). Figure 10 shows $f\left({E}_{a},{E}_{r}\right)$, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ when = 2 × 10⁴ (upper-left frame), 2 × 10⁵(upper-middle frame), 5 × 10⁵ (upper-right frame), 10⁶ (lower-left frame) and 2 × 10⁶ (lower-right frame) in the case of q = 2. We readily confirm that $f\left({E}_{r}\right)$ versus E_r approach Γ-distribution with the fat tail, as decreases. In short, deviations of $f\left({E}_{r}\right)$ from Γ-distribution in the regime of low E_r decrease, as increases. As a result, we can conclude that data-set $\left({E}_{a},{E}_{r}\right)$, whose total number at $\left({x}_{o},{y}_{o}\right)$ is small, attributes deviation of $f\left({E}_{r}\right)$ from Γ-distribution. On the other hand, the high E_a -tail of $f\left({E}_{a}\right)$ decreases markedly, when = 2 × 10⁶.

Zoom In Zoom Out Reset image size

We are interested in the change in the statistical characteristics of earthquakes owing to the Great East Japan Earthquake on 11th Mar. in 2011. Hence, differences between $f\left({E}_{a},{E}_{r}\right)$, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ 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 2011. Figures 11 and 12 show $f\left({E}_{a},{E}_{r}\right)$, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ 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 2011, respectively, together with $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ in figure 6. We can readily confirm that the variance of $f\left({E}_{a},{E}_{r}\right)$ obtained using data-set before the Great East Japan Earthquake on 11th Mar. in 2011 is larger than $f\left({E}_{a},{E}_{r}\right)$ obtained using data-set after the Great East Japan Earthquake on 11th Mar. in 2011 in all the cases of q, as shown in figures 11 and 12. Actually, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ obtained using data-set before the Great East Japan Earthquake on 11th Mar. in 2011are lower than those in figure 6 in their high energy-tails in all the cases of q. Meanwhile, $f\left({E}_{a}\right)$ and $f\left({E}_{r}\right)$ obtained using data-set after the Great East Japan Earthquake on 11th Mar. in 2011 is lower than those in figure 6 in the low-energy-domains of E_a and E_r . Consequently, we can conclude that the earthquakes with large E_a and E_r occur after the Great East Japan Earthquake on 11th Mar. in 2011, exclusively, whereas the earthquakes with small E_a and E_r occur before the Great East Japan Earthquake on 11th Mar. in 2011.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size