Abstract:
A method of predicting an earthquake in a seismically active region below an inosphere, including the steps of measuring ( 10, 12 ) a relative fluctuation of plasma density ( 26 ) in the inosphere ( 20 ), inferring a relative amplitude of an acousto-gravity wave in the inosphere ( 20 ), and inferring an earthquake magnitude from the relative amplitude of the acousto-gravity wave.

Description:
FIELD AND BACKGROUND OF THE INVENTION 
     The present invention relates to a method for predicting earthquakes. 
     Large earthquakes are among the most destructive of natural disasters. The toll that these disasters have taken in both lives and property is well known and need not be recited here. Many methods have been proposed in the past for predicting earthquakes, to provide enough warning for the people affected to prepare. One geophysical method, for example, based on the gravitational field turbulence associated with a seismic event, was proposed by Chiba (Jiro Chiba,  Proceedings of The Institute of Electrical And Electronic Engineers,  1993  International Carnahan Conference on Security Technology , pp. 215-218). Another geophysical method, based on measurements of natural low frequency radio signals, was proposed by Hayakawa (M. Hayakawa,  Phys. Earth. and Plan. Int. , vol. 77 pp. 127-135 (1993)). These methods provide only about 20 seconds advance notice of an earthquake. 
     There is thus a widely recognized need for, and it would be highly advantageous to have, a method for predicting earthquakes that provides more advance warning than known methods. 
     SUMMARY OF THE INVENTION 
     According to the present invention there is provided a method of predicting an earthquake in a seismically active region below an ionosphere, comprising the steps of: (a) measuring a relative fluctuation of plasma density in the ionosphere; (b) inferring a relative amplitude of an acoustic-gravity wave in the ionosphere; and (c) inferring an earthquake magnitude from the relative amplitude of the acoustic-gravity wave. 
     The present invention is based on the discovery of an earthquake precursor that appears in the ionosphere above the epicenter several hours before the earthquake. FIGS. 1A through 1F show two way travel times of reflections from the ionosphere, as a function of frequency, measured by a vertical ionosonde, in the 24 hours before an earthquake that occurred in the Vrancha region of Romania on May 30-31, 1990. FIG. 1A shows the structure of the normal ionosphere, 24 hours before the earthquake. Only one reflecting layer is present, and it reflects only at frequencies up to about 15 MHz. As the time of the earthquake approaches, more reflecting layers develop, and they reflect over a wider range of frequencies. Because the reflection frequency is proportional to the square root of the ion density, this wider range of frequencies indicates that the ionosphere becomes more inhomogeneous above the epicenter in the 12 to 16 hours before the earthquake. These inhomogeneities are believed to be produced by acoustic-gravity waves (AGW) associated with the buildup of strain in the Earth&#39;s crust immediately before the earthquake. 
     This phenomenon provides up to 20 hours warning in advance of an earthquake, but it does not provide a measure of the magnitude of the impending earthquake. According to the present invention, a measure of the magnitude of the earthquake is provided by measurements of oblique scattering of radio waves from the ionospheric inhomogeneities. These measurements are interpreted using a theory of the coupling of AGW and ionospheric plasma density that is presented herein. 
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS 
     The invention is herein described, by way of example only, with reference to the accompanying drawings, wherein: 
     FIGS. 1A through 1F show reflections from the ionosphere, as a function of frequency, at 2 hour decrements in the 24 hours before an earthquake; 
     FIG. 2 is a schematic illustration of the measurement geometry of the present invention; 
     FIG. 3 is a flow chart of the present invention. 
    
    
     MODULATION OF PLASMA DENSITY BY ACOUSTIC-GRAVITY WAVES 
     During its movement through the ionospheric cold dense plasma, an AGW of frequency ω and wavenumber vector {right arrow over (k)} causes the neutral particle density N m  and the rate of ionization per unit volume q({right arrow over (r)}) to vary as follows: 
     
       
           N   m   =N   m0   +δN   m  exp[ i ( {right arrow over (k)}·{right arrow over (r)}−ωt )]  (1)  
       
     
     
       
           q ( {right arrow over (r)} )= q   0   +δq  exp[ i ( {right arrow over (k)}·{right arrow over (r)}−ωt )]  (2)  
       
     
     (A. Gossard and Y. Huk,  Waves in the Atmosphere , 1975). Here, δN m  and δq are the disturbances of neutral particle density and of the rate of ionization per unit volume due to the AGW, {right arrow over (r)} is the position vector, and t is time. Let us now estimate the changes in plasma concentration caused by the AGW. Here we must consider two cases, the electrostatic case and the electrodynamic case. 
     In the D-layer and lower E-layer, the problem is electrostatic. At these altitudes, the ionization rate is as in equation (2), but the recombination rate is of the form α ef N e   2 , where α ef  is the recombination coefficient, on the order of 0.2×10 −6  cm 3 s −1  to 10 −6  cm 3 s −1 . At these altitudes (50-100 km) ions move with the speed of neutral particles, V i ≈V m ≈50-100 m/s, and the speed of charged particles is less than the speed of the AGW. In other words, the charged particle motion has no influence on the change of plasma concentration, i.e., V i &lt;&lt;V A , where V A  is the speed of sound. In this case, from the equation of conservation of particles,                  ∂   N       ∂   t       =     q   -       α   ef          N   2       -     div        (     N          V   →     α       )                 (   3   )                                
     it follows that the characteristic time of change of particle density due to recombination, τ R , is much less than the characteristic time of particle transport, τ T . In fact, for N e0 =N i0 =N 0 =10 10 m −3 , V m ≈50 m/s and the characteristic scale of quasiregular plasma inhomogeneity Λ≈200 km at these altitudes, we obtain τ R =(α ef N) −1 ≈100 s-500 s and τ T =Λ/V m ≈4000 s. With τ T &lt;&lt;τ R , we can ignore the movement of charged particles in (3), i.e.,                  ∂   N       ∂   t       =     q   -       α   ef          N   2                 (3a)                                
     Introducing (2) to (3a) and linearizing, we find that the amplitude of plasma inhomogeneities due to AGW can be estimated as:                δ                   N     e   ,   i         =         N     e   ,   i       -     N       0      e     ,   i         ≈     δ                 q      2        α   ef            N       0      e     ,   i           ω   2     +       (     2        α   ef          N       0      e     ,   i         )     2                     (   4   )                                
     Here, q 0 =α ef N 0e   2  (N 0e ≈N 0i , plasma is quasi-neutral), δq=q 0 δN m /N m0 =q 0 δ A , where δ A  is the relative amplitude of the AGW. Using this notation, we finally have:                δ                   N     e   ,   i         =       δ   A        2        α   ef            N       0      e     ,   i     3         ω   2     +       (     2        α   ef          N       0      e     ,   i         )     2                   (   5   )                                
     for AGW with frequency ω&lt;2α ef N 0 , δN/N 0 =δ A /2. Hence, for lower frequencies, the plasma disturbances are proportional to their source, i.e., to the amplitude of the AGW. For higher frequencies, i.e., ω&gt;2αefN 0 , δN/N 0 =2δ A F(ω), where F(ω)=(α ef N 0 /ω) 2 . Because in the D-layer and the lower E-layer, α ef N 0 =0.002 s −1 -0.01 s&lt;&lt;ω, F(ω)&lt;&lt;1, and δN/N 0 &lt;&lt;δ A . Hence, for high frequencies, plasma disturbances are weak, i.e., much less than the amplitude of the AGW. 
     In the F-layer, the problem is electrodynamic. For these altitudes, the amplitude of the moving plasma inhomogeneities depends on the orientation of the wave vector {right arrow over (k)} of the AGW relative to the geomagnetic field {right arrow over (B)} 0  and relative to the plasma drift velocity {right arrow over (V)} d . As was shown in C. O. Hines,  Can. J. Phys . Vol. 38 pp. 1441-1481 (1960), the amplitude of moving plasma disturbances (MPD) is maximal when the condition of spatial resonance is “working” (i.e., the plasma drift velocity is equal to the AGW phase velocity): 
     
       
         ω= {right arrow over (k)}·{right arrow over (V)}   d   (6)  
       
     
     The MPD are created by the transport processes during interactions of charged particles with neutrals. The latter are modulated according to equation (1) by the AGW, so the charged particles also are modulated by the AGW. Because at these altitudes (&gt;150 km) magnetic and electric fields are significant, the interaction between the AGW and the plasma is electrodynamic. In this case, because the frequency of ion-neutral interactions, ν im , is on the order of N m (T m +T i ) ½ , where T m  is the temperature of the neutral species and T i  is the temperature of the ions, the modulation of N m  and T m  by AGW produces a corresponding change in ν im  and hence the corresponding change in the ion velocity {right arrow over (V)} i . This modulation of {right arrow over (V)} i  causes the redistribution of plasma density and the creation of δN. 
     Let us consider that the plasma is isothermal, i.e., T e =T i =T, which is the case in the lower ionosphere, up to about 200 km, and choose a coordinate system in which the magnetic field {right arrow over (B)} 0  is parallel to the z-axis. The particle movement is considered in a coordinate system (x,y,z) which is at rest with respect to average neutral flow. The background plasma is quasi-regular and homogeneous. If {right arrow over (V)} 1   m  is the velocity of neutral species in the field of the AGW, {right arrow over (k)} is the wavevector of the AGW, δ A =δN m /N m0 =V 1   m /V ph , where V ph =ω/k is the phase velocity of the AGW, and δ T =δT m /T m0  is the relative amplitude of the modulation of the temperature by the AGW, then we can accurately represent the frequency of interaction as a periodic function due to modulation cause by the AGW (D. F. Martin,  Proc. Roy. Soc. Canada , Vol. A209, p. 216 (1950)): 
     
       
         ν im =ν im   0 [1+(δ A +δ T /2)exp{ i ( {right arrow over (k)}·{right arrow over (r)}−ωt )}]=ν im   0 +δν im   (7a)  
       
     
     
       
         ν ei =ν ei   0 [1+(δ i +3δ T /2)exp{ i ( {right arrow over (k)}·{right arrow over (r)}−ωt )}]=ν ei   0 +δν ei   (7a)  
       
     
     Because of the difference between {right arrow over (V)} e  and {right arrow over (V)} i  in the magnetic ({right arrow over (B)} 0 ) and electric ({right arrow over (E)} 0 ) fields (E 0 ≈10 mV/m-50 mV/m) between particles (electrons and ions), the ambipolar field is created, i.e., δ{right arrow over (E)}=−∇φ (φ is the potential of the ambipolar field), and also the changes of the tensors of particle mobility: 
     
       
         
           {circumflex over (M)} 
           α 
           ={circumflex over (M)} 
           0 
           α 
           +δ{circumflex over (M)} 
           α 
           ,α=e,i  
         
       
     
     
       
           {right arrow over (E)}={right arrow over (E)}   0   +δ{right arrow over (E)}   (8)  
       
     
     where                δ                     M   →     α       =         δ                   ν   α         ν   α   0            (           M   p   α         0       0           0         M   p   α         0           0       0         -     M   ∥   α             )               (   9   )                                
     Here, p denotes the Pedersen component of electrons (α=e) and ions (α=i) mobility, i.e., in the direction of the electric field, and M ∥   α  is the mobility of charged particles along the magnetic field. Finally, in the electrodynamic case, we can present (taking into account that all disturbances are weaker than the background values) in the linear case, the disturbances of charged particle (electrons and ions) velocities (M. G. Gel&#39;berg,  Geom. And Aeron . Vol. 20, pp. 271-274 (1980)).                δ                     V   →     α       =       q        (       δ                     M   ^     α            E   →     0       +         M   ^     0   α        δ                   E   →         )       +         B   0       Q   α              M   ^     0   α        δ                     V   →     m       -           D   α          B   0         Q   α              M   ^     0   α          ∇   δ                     N     N   0                   (   10   )                                
     Here, D α  is the coefficient of diffusion of electrons (α=e) or ions (α=i); Q α =ω Hα /ν α , ω Hα  is the gyrofrequency if plasma particles in the magnetic field; ν e =ν em +ν ei ; ν i =ν im +ν ei . In the case when disturbances of the plasma particles due to their modulation by the AGW are weak, i.e., δN&lt;&lt;N 0 , we can in equation (3) exclude the ambipolar field. Linearizing equation (3) and taking into account equation (10), an equation can be obtained which, for longitudinal AGW, gives                  δ                 N       N   0       =       (     ω   -       k   _     ·       V   _     d         )                δ   A          kV   ph          cos   2        χ     +       (       δ   A     +       δ   T     /   2       )          (         k   →     ·       E   →     0           Q   i          B   o         )               (     ω   -       k   →     ·       V   →     d         )     2     +     (       D   α          k   ∥   2       )                   (   11   )                                
     Here, {right arrow over (V)} d  is the drift velocity of charged particles, and χ is the angle between {right arrow over (k)} and {right arrow over (B)} 0 . From (11) it follows that if Q i V ph /V d cos 2 χ&lt;&lt;1, then the influence of modulation exceeds the effects of interaction between the plasma and the neutral species. In this case, the angle ψ={fraction (π/2+L )}−χ is less than ψ k =sin −1 {square root over (V d +L /Q i +L V ph +L )}. For an altitude of 200 km, Q i =10 −1  to 10 −2 , V ph ˜600 ms −1 , E 0 =30 mV/m, and we obtain ψ k ˜7° to 20°. Expression (11) has a maximum for AGW satisfying the condition ω−{right arrow over (k)}·{right arrow over (V)} d ≈D α k ∥   2 , where k ∥  is the component of {right arrow over (k)} parallel to the magnetic field:                  (       δ                 N       N   0       )     max     =           δ   A        k                   V     p                 h            cos   2        χ     +       (       δ   A     +       δ   T     /   2       )          (         k   →     ·         E   →     0     /     Q   i              B   0       )             D   α          k   ∥   2                 (   12   )                                
     The above analyses illustrate the principal possibility of electrodynamic and electrostatic mechanisms of the redistribution of plasma by longitudinal AGW and, finally, of the creation of plasma inhomogeneities. 
     DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     The present invention is of a method of predicting earthquakes as much as 12-16 hours in advance. 
     The principles and operation of earthquake prediction according to the present invention may be better understood with reference to the drawings and the accompanying description. 
     The present invention is based on N. D. Philipp and N. Sh. Blaunshtein, Power of H E -Scatter Signals,  Proc. XII Union Conference on Radiowave Propagation, Kharkov , USSR, 1978, pp. 166-169, which is incorporated by reference for all purposes as if fully set forth herein. 
     Referring now to the drawings, FIG. 2 illustrates schematically the geometric setup of the apparati used in the present invention. A radio frequency transmitter  10  and a radio frequency receiver  12  are separated by up to 1000 km or more on the surface of the Earth. Together, transmitter  10  and receiver  12  constitute an oblique ionosonde. Transmitter  10  transmits radio frequency energy in a solid angle  22 . Receiver  12  receives radio frequency energy from a solid angle  24 . In particular, receiver  12  receives radio frequency energy from transmitter  10  that is reflected from anomalous plasma inhomogeneities in the portion of an ionospheric anomaly  20  that lies within the volume of overlap  26  of solid angles  22  and  24 . The separation of transmitter  10  and receiver  12  is chosen to be consistent with the elevations of solid angles  22  and  24 , the power of transmitter  10  and the sensitivity of receiver  12 ; the preferred separations are between 500 km and 1000 km. The radio frequency energy should be transmitted at a frequency greater than the ionospheric plasma frequencies, so frequencies greater than about 20 MHz are used. 
     Between transmitter  10  and receiver  12  are deployed one or more vertical ionosondes  14 . When an anomalous reflection pattern such as that illustrated in FIG. 1D is detected by one of vertical ionosondes  14 , transmitter  10  and receiver  12  are activated to measure the relative fluctuation of plasma density, δN/N 0 , of equation (11), using equation (6) of Philipp and Blaunshtein:                P   R     =       P   T                αλ   3          (     δ                   N   /     N   0         )       2       32      e          2      π            r   4                ∫   V                G   T          G   R          sin   2           sin   2          (     θ   /   2     )              exp        (       -   8          π   2          α   2          ψ   2          sin   2          θ   2       )               V                   (   13   )                                
     In this equation, P T  is the transmitter power, P R  is the receiver power, r is the distance from transmitter  10  and receiver  12  to overlap volume  26  (if transmitter  10  and receiver  20  are not equidistant from overlap volume  26 , the term r 4  should be replaced with r T   2 r R   2 , where r T  is the distance from transmitter  10  to overlap volume  26  and r R  is the distance from receiver  12  to overlap volume  26 ), λ is the wavelength of the radio frequency radiation, G T  is the transmitter gain, G R  is the receiver gain, the angles θ, ψ and χ are defined in Philipp and Blaunshtein, and the integration is performed numerically over overlap volume  26 . (The relative fluctuation of plasma density is denoted as δN/N 0  herein, instead of the notation “ΔN/N” used in Philipp and Blaunshtein.) Note that the latest arrival time measured by vertical ionosonde  14  provides an upper bound on the altitude of anomaly  20 , and the integration in equation (13) should be limited to altitudes below this upper bound. α is the ratio L/λ, where L is the longitudinal extent (along the earth&#39;s magnetic field) of the inhomogeneities in anomaly  20 . This ratio may be estimated using equation (11) of Philipp and Blaunshtein, as described therein. 
     Solid angles  22  and  24  are scanned to provide a map of δN/N 0  in the ionosphere between transmitter  10  and receiver  12 . The epicenter of the impending earthquake is predicted to be directly below the point of maximum δN/N 0 , to within about 5 km. δ A  is estimated using equation (12), on the assumption that δ T  is much smaller than δ A . The ionospheric parameters of equation (12), such as wave vector {right arrow over (k)}, are well known to those ordinarily skilled in the art, and may be found, for example, in M. G. Gel&#39;berg,  The Inhomogeneities of the High - Latitude Ionosphere , Nauka, Novosibirsk, 1986. To ensure the validity of an analysis based on equation (12), the integration in equation (13) is limited to altitudes at which the electrodynamic model is valid, i.e., above about 80 km. It has been found empirically, for low-magnitude earthquakes, that for earthquakes of magnitude greater than about 3 on the Richter scale, the earthquake magnitude is approximately 0.8 times δ A . If the predicted magnitude is greater than about 5, the populace near the predicted epicenter is warned to evacuate or take other precautionary measures. 
     Earthquakes tend to occur in fault zones, which are linear geological features. Therefore, transmitter  10  and receiver  12  preferably are deployed at opposite ends of a fault zone, and are used to probe the ionosphere above the fault zone. For example, to monitor the San Andreas fault of California, transmitter  10  may be deployed in San Diego and receiver  12  may be deployed in San Francisco (or vice versa). In the case of a relatively restricted area of seismic activity, such as the Vrancha Region of Romania, transmitter  10  and receiver  12  may be deployed on opposite sides of the area of seismic activity. In the latter case, it is not necessary to scan solid angles  22  and  24 . Instead, solid angles  22  and  24  are kept fixed, and are made wide enough for overlap volume  26  to cover the entire area of seismic activity. 
     As an alternative to deploying an oblique ionosonde, a sufficiently dense array of vertical ionosondes  14  may be deployed, with each vertical ionosonde  14 , in addition to monitoring the reflectivity structure of the ionosphere thereabove, also serving as a combined transmitter and receiver to measure δN/N 0  in a solid angle directly thereabove. The upper limit of integration for equation (13) is provided by the highest altitude inferred from the maximum two way travel time of the reflection patterns received by vertical ionosondes  14 . The lower limit of integration for equation (13) is provided by either the lowest altitude inferred from the minimum two way travel time of the reflection patterns received by vertical ionosondes  14 , or by the lowest altitude of validity of the electrodynamic model (˜80 km), which ever is higher. 
     FIG. 3 is a flow chart of the present invention. In box  100 , vertical ionosondes  14  monitor the reflectivity structure of the ionosphere for anomalies. Upon detection of an anomalous reflection pattern (box  102 ), transmitter  10  and Receiver  12  are activated to measure δN/N 0  (box  104 ). δ A  is estimated from δN/N 0  in box  106 , and the magnitude of the earthquake is estimated in box  108 . 
     While the invention has been described with respect to a limited number of embodiments, it will be appreciated that many variations, modifications and other applications of the invention may be made.