Underlying mechanism of precursory activity from analysis of upward earthquake migration

In this paper we analyse the upward earthquake hypocentral migration in the ten known subduction zones and discuss a possible mechanism of such migration. The total time of the migration appears to range from 2.5 to 10 years. It leads to the estimation of the average velocity V z ∼ 60−300 km yr−1. It probably corresponds to the movement of the forcing agent like stress or deformation wave from depths of the upper mantle (600–700 km) to the level of the lithosphere with subsequent initiation of fluid migration inside the crust to trigger shallow earthquakes. Averaged over all zones upward migration travel time is about 5 years (< V z >≈ 120 km yr−1) that coincides approximately with the period of characteristic temperature variation (El Nino) and crustal seismic periodicity in the Pacific region. These findings are helpful for the study of the seismic precursors and analysis of earthquake triggering.


Introduction
Analysis of precursors of the large earthquakes (EQs hereafter) both on the ground surface and inside atmosphere and ionosphere leads to a conclusion on the basic role of the gas/water release before and after EQs (e.g., Molchanov and Hayakawa, 2008). On the other hand, the geodynamic specialists put forward a concept of gas plumes, which originated from the Earth's mantle and are responsible for the transfer of the heat and juvenile gases into the atmosphere (Puscharovsky, 1998;Letnikov, 2002;Osika and Cherkashin, 2008). So the emergence of the EQ precusors can be considered as a result of such gas plume intrusion into the lithosphere and the consequent triggering of EQs Correspondence to: O. A. Molchanov (olmolchanov@mail.ru) there. Main features of the short-term precursors can be explained in this way, i.e., their mosaic distribution in the zone much larger than the size of EQ source (in a so-called preparation zone) and their queue sequence: the first are hydrology/geochemistry precursors (several weeks prior to the main shock), then electromagnetic precursors (several days before) and the last are seismo-acoustic and probably biological precursors (several hours earlier than EQ). At present, a development of preliminary gas/water release is treated in almost every theoretical model of the short-term precursors (see discussion in the book by Molchanov and Hayakawa, 2008).
Tracing of the Earth's interior by seismicity is a conventional method. For example the well-known model of subduction zone, which supposes a curling of the ocean plate down to the depth about 600-700 km, was originated due to such tracing results (Uyeda and Kanamori, 1979;Uyeda, 1982). The main part of the seismic and volcanic activity is concentrated in the subduction zones; therefore, it is natural to analyse the possible formation of gaseous plumes just in these zones. Plume parameters are not well established until now. It can be supposed that their maximum size is of the order of subduction zone and their structure is similar to the "overturned cone" (Osika and Cherkashin, 2008). What is also unclear is the question, are they static or dynamic structures (for example "pulsating plumes"). Probably Japanese scientist Mogi (Mogi, 1973(Mogi, , 2004 was the first who considered the intermittent sequence of several deep and crustal large EQs in the Japan subduction zone, occurring with the time delay of several years. However, Mogi discussed neither a reason for such a relation nor the cause-and-effect direction, i.e., was there evidence of descending or ascending migration of the agent triggering the EQs. Recently Molchanov and Uyeda (2009) have analysed spatial-time distribution of the EQ hypocentres in the Japan, Kurile-Kamchatka and Sunda subduction zones for the depths from 700 km to the surface and found the prevalence of ascending migration with the time delay from 2 to 10 years. Furthermore, by using correlation analyses of EQs and climate in the near-equatorial Pacific area during the last 35 years, Molchanov (2010) demonstrated the evidence that the heat flux from seismicity could induce the longterm climate variations (El Nino effect) and even global warming. We analyse in this paper the EQ hypocentral migration in the ten known subduction zones and discuss a possible mechanism of such migration.

Preparation of data for analysis
We use the data from USGS catalogue covering time period from 1 January 1973 to 1 January 2008 (total duration is T = 35 years) relevant to ten zones with essentially deep seismic activity.
Nine zones are situated in the Pacific area forming the socalled "fire ring" of the seismic and volcanic activity (Fig. 1). The only zone outside it is the Italy-Greece zone.
We divide each zone into 5 layers. By taking into account that usual lower boundary of crust in the Pacific region is near 33 km we change a little the USGS classification as follows: -Layer 1 (crust): d ≤ 38 km.
In each layer, we calculate the seimic energy flux U by averaging it in the time intervals of t = 6 month duration each starting at t k (we have k ∈ [1,70] for the total time T = 35 years): U (t k ,rk,t k ) = Es N k (t k ,t k + t;rk,Es,Es + Es)

≈ ∫Es(∂N/∂Es)dEs
Here rk≡x k , y k , z k are the energy weighted mean values of the flux space coordinates (see paper by Molchanov and Uyeda, 2009), seismic energy is connected with the magnitude from the EQ catalogue by the known relation (Kanamori and Anderson, 1975): and N is the number of EQs in the interval (Es, Es+ Es) or in the corresponding interval (M, M + M). For all layers j = [1,5] we have 350 values U (t k ,r k ,j ) in zone. Then we compute the average energy flux in each layer <U (t k ,rk,j ) >= U a(j ). These values are given in the Table 1 together with additional information. As before and even the global warming . We analyze in this paper the EQ hypocentral migration in the ten known subduction zones and discuss a possible mechanism of such migration.

Preparation of data for analysis
We use the data from USGS catalog covering time period from 1.01.1973 to 1.01.2008 (total duration is Т=35 years) relevant to ten zones with essentially deep seismic activity.  (Molchanov and Uyeda, 2009;Molchanov, 2010), we calculate indices of seismic variability: e(t k ,r k ,j ) = log(U (t k ,r k ,j ))− < log(U (t k ,r k ,j )) > = log(U (t k ,r k ,j )/U n(j )) where U a (j )/U n (j ) ratio is about 1.5. It is rather clear why we use the relative values of seismic energy flux. We are interested in the amplitude variation of triggering agent A(t k ,r k ,j ), which is connected with U (t k ,r k ,j ) through an induction coefficient γ (j ). The induction greatly depends on the medium properties. By supposing that: where α is some exponent, we readily find: Thus, we assume that the change of average triggering amplitude A(j ) in the course of migration from layer to layer is much smaller than change in γ (j ). In this condition, a disclosure of the agent migration from seimic tracing could be a success. Concerning relation (2), we follow the classic tradition of measuring the seismic activity on a logarithmic  scale. Here, unlike our previous paper (Molchanov, 2010), we analyse the unipolar seismic index values e c (t k ,r k ,j ), using simple transform: These unipolar indices are practically coinciding with e(t k ,r k ,j ) for the large values of the seismic energy flux but diminish significantly with the low energy values.

Method of analysis of plume dynamics
Let us have some change of the plume amplitude at the depth d 0 , as А (t, d 0 ). It is known that the amplitude variation in the selected time t = t 0 , is not random fluctuation if the classic criterion is valid: where σ( A(t, d 0 )) is r.m.s. We regard such a perturbation as "packet". If the packet goes up with velocity V, then it appears at the depth d<d 0 in time t = t 0 + t, t = (d 0 -d)/V, i.e the packet propagates along some "trajectory" which is described as follows: Since the nature of the triggering agent is not clear at the moment, we regard it for simplicity as "plume". Some results of seismic tracing for the whole time of analysis, i.e., the values e c (x kj , y kj , z kj ) are presented in Fig. 2. It is reasonable to accept this distribution as the averaged plume pattern. If so, the plume looks as an inverted cone (turned upside), or sometimes as a "flower growing from the mantle". Anyway, the size of the selected zone included onto the analysis exceeds the plume volume in all cases.

Method of analysis of plume dynamics
Let us change some of the plume amplitude at the depth d 0 , as A(t,d 0 ). It is known that the amplitude variation in the is not a random fluctuation if the classic criterion is valid: where σ ( A(t,d 0 )) is r.m.s. We regard such a perturbation as a "packet". If the packet goes up with velocity V , then it appears at the depth , the packet propagates along some "trajectory" which is described as follows: where τ 0 is the time delay of travelling from initial depth d 0 to final depth in the crust d c . If we know the amplitude variation both at d 0 and at d c , then the simplest way to reveal the travelling packet is to calculate the correlation function: which has a maximum at τ = τ 0 if the motion really exists.   Fig.3a.
Surprisingly, these pictures remind of the tracks left by energetic particles in the Wilson bubble chamber.

Results of correlation analysis
We provide below the results of the last, the most informative method of analysis. Some results are presented in Fig.4. An obvious drawback of the method is a necessity to suggest a particular trajectory model. We believe that the most realistic model is the diffusion one (see Discussion section) . However, we computations for all the three trajectory models and show their averaged results together with results of computation in the diffusion model. Moreover, when comparing the data, we apply the similarity criterion for the trajectory model chosen, which is used together with the basic criterion of 2σ level excess. Hence, when plot of normalized correlation for the diffusion/averaged model or for the both exceeds +2 level at = max value and the clear maximum is observed for the both types of the trajectories , we accept it as reliable evidence of the upward/downward migration with the average migration velocity <Vm> = (z 0 -z c )/ max ≈ 600 km/ max. Let us remind that <Vm>, max > for the upward migration and for the downward migration they are negative. at τ = τ 0 + (t 1 − t 0 ) and at τ = τ 0 − (t 1 − t 0 ) in addition to maximum at τ = τ 0 . Note that in a special case of maxima placed at the layer d 0 as t 1 − t 0 = 2τ 0 , the symmetrical correlation humps appear at values τ = ±τ 0 , ±3τ 0 , etc. When the analysed interval is large, i.e., T τ 0 and the maxima in the source layer are placed randomly, the false correlation humps are suppressed. However, the time interval in our case is not so large and we need some criteria showing reliability of correlation humps, as well as improvement of the method itself. Concerning the reliability, we use the relation similar to (5) and accept the correlation hump as a real one when the normalized correlation function exceeds the usual reliability threshold: where the averaging on τ is supposed.
In correspondence with relations (3) and (4) and using for simplicity α = 0.5 we assume: It is clear A(t k ) ≈ (U (t k )/U n ) 1/2 when U (t k ) > U n . Note that relation (7) is an estimate of one-dimension correlation of two parameters. So, d 0 is the averaged depth of the lower layer 5 and d c is the averaged depth in the upper layer 1. In principle we can include several layers into consideration. For example, by accounting for an intermediate layer with the averaged depth d m in addition to upper layer (d c ) and lower layer (d 0 ) we can analyse the one-dimension normalized correlation of three parameters by using equation similar to (7): Here, we meet with the problem of determination, the value τ m being actually the problem of the packet trajectory. By supposing the uniform convection, we obtain τ m = τ * (d 0 − d m )/(d 0 − d c ) and for average values d 0 = 400 km, d c = 30 km and d m = 175 km, we have τ m = τ/2, which we use for an estimate of the 3-parameter correlation. Such an approximation is better than a 2-parameter one, but it does not contain information on the real distribution of the seismic indices in space. The most important is the distribution on the depth. Finally, we analyse normalized 5-layer and two dimension correlation which is described in details in the paper by Molchanov and Uyeda (2009). We use indices of seismic variability in each layer e c (t k ,d k ,j ) that are averaged over the latitude and longitude but not over the depth. Concerning the trajectory, we consider three models: 1. Uniform convection: d(t,τ ) = d 5 − a 1 (t − t k ), a 1 = V c = (d 5 − d 1 )/τ , where d 5 is lower boundary of layer 5 (650 km), d 1 is upper boundary 1 (10 km).

K u r -K a m J a p a n M a r ia n a S u n d a P h ilip p in e B o u g a in v il T o n g a N . T o n g a S . C h ile
It a ly -1 2 -1 0 -8 -6 -4 -2 0 2 4 6 8 1 0 1 2 T im e d e la y , y e a r s -2 0 2 4 C n /rm s (C n ) a n d a v rC n /rm s (a vrC n ),3 m o d e ls ,1 0 z o n e s.

K u r -K a m J a p a n M a r ia n a S u n d a P h ilip p in e B o u g a in v il T o n g a N . T o n g a S . C h ile
It a ly Fig.4a . Results of correlation analysis on the sense of migration and averaged velocity for the following zones using diffusion trajectory model: a) Kurile-Kamchatka(open circles); b) Japan (closed circles); c) Mariana(opened diamonds) ; d) Sunda (closed diamonds); e)Philippines( opened squares); f) Bougainville (closed squares); g)Tonga North (opened triangles); h) Tonga South (closed triangles); i) Chile( opened stars); j)Italy (closed stars). Averaging over all zones is shown by solid line. Estimation of direction and travel time of the migration is shown by arrow using reliability criteria that discussed in the text. Fig.4b. Results of correlation analysis for all zones using averaging over three trajectory models. Legends are the same as in Fig. 4a.
It is evident that upward migration prevails in the seven zones, while both the upward and downward migration is found in two zones (Tonga South and Chile) and the result is not reliable for the Mariana zone. The results are summarized in the Table 2.

K u r -K a m J a p a n M a r ia n a S u n d a P h ilip p in e B o u g a in v il T o n g a N . T o n g a S . C h ile
It a ly -1 2 -1 0 -8 -6 -4 -2 0 2 4 6 8 1 0 1 2 T im e d e la y , y e a r s -2 0 2 4 C n /rm s (C n ) a n d a v rC n /rm s (a vrC n ),3 m o d e ls ,1 0 z o n e s.

K u r -K a m J a p a n M a r ia n a S u n d a P h ilip p in e B o u g a in v il T o n g a N . T o n g a S . C h ile
It a ly Fig.4a . Results of correlation analysis on the sense of migration and averaged velocity for the following zones using diffusion trajectory model: a) Kurile-Kamchatka(open circles); b) Japan (closed circles); c) Mariana(opened diamonds) ; d) Sunda (closed diamonds); e)Philippines( opened squares); f) Bougainville (closed squares); g)Tonga North (opened triangles); h) Tonga South (closed triangles); i) Chile( opened stars); j)Italy (closed stars). Averaging over all zones is shown by solid line. Estimation of direction and travel time of the migration is shown by arrow using reliability criteria that discussed in the text. Fig.4b. Results of correlation analysis for all zones using averaging over three trajectory models. Legends are the same as in Fig. 4a.
It is evident that upward migration prevails in the seven zones, while both the upward and downward migration is found in two zones (Tonga South and Chile) and the result is not reliable for the Mariana zone. The results are summarized in the Table 2.

Logarithmic model: log
3. Diffusion model: we suppose different velocities and trajectories for the upward motion, when a source is placed at the lower boundary d 5 and for the downward boundary, when a source is near the upper boundary d 1 . Upward motion is described by equation: and the downward motion corresponds to: Nat. Hazards Earth Syst. Sci., 11, 135-143, 2011 www.nat-hazards-earth-syst-sci.net/11/135/2011/ Positive values of velocities and delays correspond to the upward motion in all the models while the negative values describe the downward motion. Visual evidence of the packet existance can be found from the depth-time (d, t) diagrams of the variability indices for two zones that are presented in Fig. 3. Surprisingly, these pictures remind us of the tracks left by energetic particles in the Wilson bubble chamber.

Results of correlation analysis
Below, we provide the results of the last, the most informative method of analysis. Some results are presented in Fig. 4. An obvious drawback of the method is a necessity to suggest a particular trajectory model. We believe that the most realistic model is the diffusion one (see Discussion section). However, we compute all the three trajectory models and show their averaged results together with the results of the computation in the diffusion model. Moreover, when comparing the data, we apply the similarity criterion for the trajectory model chosen, which is used together with the basic criterion of 2σ level excess. Hence, when plotting normalized correlation for the diffusion/averaged model or when both exceeds +2 level at τ = τ max value and the clear maximum is observed for both types of trajectories, we accept it as reliable evidence of the upward/downward migration with the average migration velocity < V m >= (z 0 − z c )/τ max ≈ 600 km τ −1 max . Let us be reminded that < V m >, τ max > for the upward migration and for the downward migration they are negative.
It is evident that upward migration prevails in the seven zones, while both the upward and downward migration is found in two zones (Tonga South and Chile) and the result is not reliable for the Mariana zone. The results are summarized in the Table 2.
Our analysis shows that upward migration is a general property of the subduction zones. Characteristic regional velocity is also shown in Fig. 4 after averaging over all zones. Existence of regional (or global) upward migration with the average velocity V ≈ 120 km yr −1 (τ max ≈ 5 years) is clearly seen.

Discussion
The main points of the work are as follows: -Reliable upward migration is found in the 9 zones of total 10 with an exception for Mariana. The travel time varies from 2.5 years to 10 years that corresponds to averaged velocity V ≈ from 65 km yr −1 to 260 km yr −1 . Simultaneous downward and upward migration is discovered in only two zones (Tonga South and Chile) with the average downward migration velocity ≈ 130 km yr −1 and 87 km yr −1 (τ max ≈ −5 years and -10 years).
-Prevailing direction of migration as well as the plume form indicates that the source of migration is situated in the upper mantle.
-Averaged over all zones upward migration travel time is about 5 years (V ≈ 120 km yr −1 ) that coincides approximately with the period of characteristic temperature variation (El Nino) and crustal seismic periodicity in the Pacific region (Molchanov, 2010a). It invokes the possibility of the global heat transfer from the deep layers to the ground surface.
There are several problems in the interpretation of the results. The first is the nature of the gaseous plume variation or a source mechanism. It may be conceivable that in a subduction region an initial pulse could arise near the deep end of the subducting plate due to an explosion-like process related to some phase transition (e.g., Kirby, 1987;Green and Burnley, 1989). In such a case, we may expect that upward migration of EQ hypocentres follows in an analogous process if the mantle behaves visco-elastically. At depths shallower than, say, 200-300 km, free water is released by dehydration of hydrous minerals in subducting slabs, and it may play an important role in triggering the EQs by weakening the fault strength (e.g., Raleigh and Paterson, 1965;Rice, 1992;Byerlee, 1993). The second problem is a mechanism of energy transfer from the depth to the core or the nature of forcing agent that triggers the EQ. Movement of fluids below the astenosphere is hardly possible, and the heat convection is too slow. One possibility may be the deformation wave, which is popular for the explanation of EQ migration in the crust. There are several reports about horizontal migration EQ of hypocentres. The well-known examples are the case of N. Anatolia fault, Turkey, where major EQs migrated westward with a velocity of 60-70 km yr −1 (e.g., Toksoz et al., 1979), and the NW Pacific rim subduction zones where EQs moved with a velocity of 150-270 km yr −1 (Mogi, 1968). Recently, Mukhamediev et al. (2008) reported the westward EQs migration from Mid-Atlantic ridge with a velocity of 200-300 km yr −1 . Their theoretical approach explaining the horizontal deformation wave analyses the diffusion equation for the motion of displacement u x inside the plate induced by initial co-seismic slip u 0 x = u x (t = 0): where the diffusion coefficient D x is as follows: Here h L and h A are the thickness of lithosphere and asthenosphere, respectively, µ is the rigidity, η is the asthenosphere dynamic viscosity, and ν is the Poisson coefficient. The same equation is valid for horizontal stress diffusion. Parameter η here is not well determined but Mukhamediev et al. (2008) reviewed several papers on this subject in which a velocity of the observational horizontal migration V x varies from 60 to 300 km yr −1 and concluded that the most probable value of η ≈ 10 17 Pa s. It leads to the estimate D x ≈ 2700 m 2 s −1 for their accepted values h L = 100 km, h A = 30 km, ν = 0.33, µ = 3 × 10 10 Pa. Note that we could estimate the η value by taking into account the duration of aseismic slip due to the aftershock relaxation time τ 0 ≈ 2η/µ (see e.g., Scholz, 1990). In this case: As a result, we find about the same D x value when assuming the observed value of τ 0 ≈ 3 months. It is rather surprising to discover the same values of the migration velocities V z ≈ 90 − −260 km yr −1 in our case of the upward migration (see Table 2). Moreover, we can estimate the averaged value of the diffusion coefficient directly from observations: < D z >≈ d 2 0 / 4 < τ max > cos 2 (I ) ≈ 3000 m 2 s −1 ∼ D x Here we suggest d 0 = 600 km, < τ max >= 5 years and dip angle of the subduction slab I = 30 • . These estimates encourage us to accept the deformation wave mechanism as a forcing agent for the EQ triggering and to consider the diffusion trajectory model as the most promising. One of the important consequences of the deformation wave concept is a possibility of self-action triggering in the crust layer where the wave induction is especially efficient. It can be found as the appearance of reflection deformation wave. Probably a presence of the downward migration is connected with such a possibility. We could explain also the observed EQ triggering periodicity in this way (Molchanov, 2010) as a result of upward and downward wave interference in the so-called lithosphere -upper mantle resonator. This possibility will be discussed in a separate paper.
Finally, in connection with the problem of precursory activity, the two types of precursors might be expected. First, the long-term deep precursory EQ preceding the crust EQs with a time lead of several years. Indeed it was the basic idea of Mogi (1973Mogi ( , 2004 papers. Secondly, it is evident that the deformation wave stimulates the upward liquid migration in the lithosphere that leads to the appearance of the short-term precursors. It is reasonable to suppose an inhomogeneous structuring of the gas-liquid volumes (so-called "bubbles"). Modelling of the process using the finite automata technique demonstrated a faster elevation of the small bubbles in contrast to the large bubble association (Korovkin et al., 2002). Therefore, the origin of foreshocks, seismo-acoustic, hydrology/chemistry and atmosphere/ionosphere precursory events can be easily explained (see discussion in Molchanov and Hayakawa, 2008).