FIELD OF INVENTION 
     The present invention relates to methods for determining the earth&#39;s subsurface resistance to the flow of electrical current, which is useful for identifying subsurface geological features such as porosity, pore fluid, and matrix properties, as well as man-made features. 
     BACKGROUND 
     Subsurface resistivity studies are often used to obtain information about subsurface geology, being considered inexpensive, fast, and non-invasive. Such studies are based on Ohm&#39;s law relating voltage (V), current (I), and electrical resistance (R) as V=IR. That is, for a given volume of material, a voltage applied across that material will result in a current flow through that material in proportion to its electrical resistance (or more generally, its electrical impedance). 
     Resistivity, typically denoted by the Greek letter “ρ,” is a material constant that is related to electrical resistance. It is more particularly a normalized form of electrical resistance that takes account of the cross-sectional area (A) of the volume—the current flowing perpendicular to this area, and the resistance decreasing in direct proportion to this area—, and the length (L) of the volume—the current flowing parallel to this length, and the resistance increasing in direct proportion to this length (ρ=RA/L). 
     Subsurface resistivity is typically measured by staking four electrodes into the earth, spaced apart from each other, using a current source to force a known current to flow through the earth between two of the electrodes, and measuring the voltage difference across the other two electrodes. The value of V/I is proportional to the resistivity (V/I=ρ L/A), and the constant of proportionality (L/A) is determined by a complex computerized calculation. 
     While such resistivity studies are relatively inexpensive and fast compared to known alternatives, they still require visiting the site and using specialized equipment to make measurements. And while thought of as being non-invasive, such resistivity studies do require staking electrodes into the ground. Accordingly, there is a need for a method for determining geological subsurface resistivity that further reduces cost and increases speed by eliminating the need to visit the site or disturb the surface. 
     SUMMARY 
     A method for determining geological subsurface resistivity. The method is based on determining the electrical resistance of a geological volume of material that has received a first lightning strike, including obtaining a first set of lightning parameters, the first set of lightning parameters including an indicium of the current of the first lightning strike at a first initial time and an indicium of the current of the first lightning strike at a first decay time subsequent to the first initial time, and making a first inference of the resistance, at least in part, from the first set of lightning parameters. 
     The first inference may be normalized so as to make it proportional to the current of the first lightning strike at the first initial time. 
     Where the geological volume has received a second lightning strike, the method may further include obtaining a second set of lightning parameters, the second set of lightning parameters including an indicium of the current of the second lightning strike at a second initial time and an indicium of the current of the second lightning strike at a second decay time subsequent to the second initial time, and making a second inference of the resistance, in part, from the second set of lightning parameters. 
     The first inference may be normalized so as to make it proportional to the current of the first lightning strike at the first initial time, and the second inference may be normalized so as to make it proportional to the current of the second lightning strike at the second initial time. 
     The first and second inferences may be averaged. The average may be determined after normalizing the first inference so as to make the first inference proportional to the current of the first lightning strike at the first initial time, and normalizing the second inference so as to make the second inference proportional to the current of the second lightning strike at the second initial time. 
     It is to be understood that this summary is provided as a means of generally determining what follows in the drawings and detailed description and is not intended to limit the scope of the invention. Objects, features and advantages of the invention will be readily understood upon consideration of the following detailed description taken in conjunction with the accompanying drawings. 
    
    
     
       BRIEF DESCRIPTION OF DRAWINGS 
         FIG. 1  is a pictorial representation of the conditions precedent to a lightning strike. 
         FIG. 2  is a schematic of an electrical circuit modeling the conditions of  FIG. 1 . 
         FIG. 3  is a schematic of the electrical circuit of  FIG. 2  at the time of initiation of a lightning strike. 
         FIG. 4  is a schematic of the electrical circuit of  FIG. 3  at times subsequent to the time of initiation of the lightning strike. 
         FIG. 5  is a histogram of reported values of peak lightning current for a number of lightning strikes in a test area. 
         FIG. 6  is a histogram of resistance values computed according to the present invention for lightning strikes represented in  FIG. 5 . 
         FIG. 7  is a histogram of resistivity values computed according to the present invention for the lightning strikes represented in  FIG. 6 . 
     
    
    
     DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS 
     The present invention resulted from the inventor&#39;s recognition that geological subsurface resistivity ought to be related to lightning data, and determination of that relationship. Heretofore, lightning had been considered a purely meteorological phenomenon. 
       FIG. 1  illustrates a natural system  10  that forms a lightning strike. A cloud  12  develops a negative charge relative to the (earth) ground  14  underneath, which correspondingly develops a positive charge, the difference in charge defining a voltage “V” between the cloud and the ground. No current flows between the earth and the cloud as the voltage V increases, until the voltage becomes high enough to ionize the atmosphere  16  between the cloud and the ground. 
       FIG. 2  illustrates a circuit approximation of the lightning system  10  according to the invention. The capacitor C (having capacitance C) stores the charges in the cloud and the ground shown in  FIG. 1 , and the resistor R (having resistance R) represents geological subsurface resistance. 
       FIG. 3  shows the circuit of  FIG. 2  as soon as the voltage V has reached the atmospheric ionization voltage, referred to herein as V o . At that instant (t o ), a lightning strike  18  discharges through the atmosphere  16  into the ground, corresponding to an initial current I o  flowing out of the capacitor, jumping the gap between point “P 1 ” and point “P 2 ,” and passing through the resistance R. 
       FIG. 4  shows the circuit of  FIG. 3  at later times (t), with voltage V t  and current I t . 
     This circuit approximation can be used to relate the geological subsurface resistance R, which can be related to resistivity ρ, to the current in a lightning strike. 
     In an electrical circuit in which the voltage across a capacitor discharges through a series resistance R, the voltage decays exponentially with a time constant RC:
 
 V   t   =V   o   e   −t/RC   Equation (1)
 
Since V t =RI t , and assuming R is a constant:
 
 I   t   =I   o   e   −t/RC   Equation (2)
 
     I o  is the peak current of the lightning strike, which is a measured and recorded parameter for lightning strikes around the world, made available to the public. For example, Vaisala owns and operates the “National Lightning Detection Network” in the U.S., and operates the Canadian Lightning Detection Network, which together provide a comprehensive historic lightning database for lightning strikes in the United States and Canada. Another reported lightning parameter is known as the “peak-to-zero” time. 
     The circuit approximation provides a means for relating the geological subsurface resistance R, which can in turn be related to resistivity ρ, to lightning parameters as follows. 
     Taking the natural logarithm of both sides of Equation (2) provides a first expression for R, i.e., 
     
       
         
           
             
               
                 
                   
                     
                       
                         ln 
                         ⁡ 
                         
                           ( 
                           
                             I 
                             t 
                           
                           ) 
                         
                       
                       = 
                       
                         
                           ln 
                           ⁡ 
                           
                             ( 
                             
                               I 
                               0 
                             
                             ) 
                           
                         
                         - 
                         
                           t 
                           / 
                           RC 
                         
                       
                     
                     ; 
                     
                       and 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       so 
                     
                   
                   ⁢ 
                   
                     
 
                   
                   ⁢ 
                   
                     R 
                     = 
                     
                       t 
                       
                         C 
                         ⁡ 
                         
                           [ 
                           
                             
                               ln 
                               ⁡ 
                               
                                 ( 
                                 
                                   I 
                                   0 
                                 
                                 ) 
                               
                             
                             - 
                             
                               ln 
                               ⁡ 
                               
                                 ( 
                                 
                                   I 
                                   t 
                                 
                                 ) 
                               
                             
                           
                           ] 
                         
                       
                     
                   
                 
               
               
                 
                   Equation 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     ( 
                     3 
                     ) 
                   
                 
               
             
           
         
       
     
     Also,
 
 C=∈A/d   Equation (4)
 
where “∈” is the dielectric constant of air (≈8.854×10 −12  farads/m); “d” is the distance between the plates of the capacitor C, corresponding to the distance between the cloud and the ground; and “A” is the effective area of the plates of the capacitor that are used to model the cloud and the ground. It is recognized that the area “A” can be approximated by a circle of diameter “d,” so that A=πd 2 . Then,
 
 C=∈πd,   Equation (5)
 
and Equation (3) becomes:
 
     
       
         
           
             
               
                 
                   R 
                   = 
                   
                     t 
                     
                       ε 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       π 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         d 
                         ⁡ 
                         
                           [ 
                           
                             
                               ln 
                               ⁡ 
                               
                                 ( 
                                 
                                   I 
                                   0 
                                 
                                 ) 
                               
                             
                             - 
                             
                               ln 
                               ⁡ 
                               
                                 ( 
                                 
                                   I 
                                   t 
                                 
                                 ) 
                               
                             
                           
                           ] 
                         
                       
                     
                   
                 
               
               
                 
                   Equation 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     ( 
                     6 
                     ) 
                   
                 
               
             
           
         
       
     
     Now note that V o =V BD  d, where V BD  is the breakdown voltage of the atmosphere  16 , which is approximately 3.0 MV/m. Since V o =I 0  R, then 
     
       
         
           
             
               
                 
                   R 
                   = 
                   
                     
                       
                         V 
                         BD 
                       
                       ⁢ 
                       d 
                     
                     
                       I 
                       0 
                     
                   
                 
               
               
                 
                   Equation 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     ( 
                     7 
                     ) 
                   
                 
               
             
           
         
       
     
     Equating the two expressions for “R” (Equations (6) and (7)) yields: 
                           V   BD     ⁢   d       I   0       =     t     ε   ⁢           ⁢   π   ⁢           ⁢     d   ⁡     [       ln   ⁡     (     I   0     )       -     ln   ⁡     (     I   t     )         ]             ⁢     
     ⁢     Or   ,     
     ⁢       d   2     =       tI   0         V   BD     ⁢   ε   ⁢           ⁢     π   ⁡     [       ln   ⁡     (     I   0     )       -     ln   ⁡     (     I   t     )         ]                       Equation   ⁢           ⁢     (   8   )                 
And from Equation (7):
 
     
       
         
           
             
               
                 
                   
                     R 
                     2 
                   
                   = 
                   
                     
                       tV 
                       BD 
                     
                     
                       
                         I 
                         0 
                       
                       ⁢ 
                       ε 
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         π 
                         ⁡ 
                         
                           [ 
                           
                             
                               ln 
                               ⁡ 
                               
                                 ( 
                                 
                                   I 
                                   0 
                                 
                                 ) 
                               
                             
                             - 
                             
                               ln 
                               ⁡ 
                               
                                 ( 
                                 
                                   I 
                                   t 
                                 
                                 ) 
                               
                             
                           
                           ] 
                         
                       
                     
                   
                 
               
               
                 
                   Equation 
                   ⁢ 
                   
                       
                   
                   ⁢ 
                   
                     ( 
                     9 
                     ) 
                   
                 
               
             
           
         
       
     
     Equation (9) solves for the resistance R in terms of how the lightning current in a particular lightning strike  18  changes over time. In particular, recall from Equation (3) that: 
     
       
         
           
             
               t 
               
                 [ 
                 
                   
                     ln 
                     ⁡ 
                     
                       ( 
                       
                         I 
                         0 
                       
                       ) 
                     
                   
                   - 
                   
                     ln 
                     ⁡ 
                     
                       ( 
                       
                         I 
                         t 
                       
                       ) 
                     
                   
                 
                 ] 
               
             
             = 
             RC 
           
         
       
     
     The quantity {t/[ln(I o )−ln(I t )]} is constant; therefore “t” can be a time corresponding to any desired current I t . It is beneficial to choose a time “t” for which measured lightning current parameters exist, particularly that time known as the “peak-to-zero” time, because this is a reported lightning parameter. 
     Corresponding to the “peak-to-zero” time is the current at this time, which may be referred to as the “zero” current with the understanding that the current is never identically zero but only approaches zero asymptotically according to Equation (3). 
     The “zero” current is not a reported lightning parameter, but the “peak-to-zero” time is necessarily the time at which the current reaches “zero” current. So in order to measure and specify a “peak-to-zero” time, it is necessary to measure and specify “zero” current. Hence, for any given peak current, the corresponding “zero” current is a known or determinable parameter. 
     In the lightning detection networks, lightning strikes are detected at multiple detection stations that triangulate the detected signals to determine the strike locations, and thereby to estimate the magnitudes of the strikes based on the magnitudes of the detected signals, falling off as the reciprocal of the square of the distance. At any particular lightning detection station, the strike may be far enough away that the signal never rises above the background noise and is therefore not detectable. At the other extreme, the strike is close enough to the lightning detection station that the signal is above the background noise and is therefore detectable all the way from its peak value to its “zero” value. Between the extremes, a strike may be close enough to the detection station for the peak current to be detectable, and for some time thereafter as the current decays, but not after it has decayed so much that it reaches “zero.” 
     It is assumed herein that if a “peak-to-zero” time is recorded for a lightning strike, the corresponding “zero” current was either actually detected or was estimated from the detected rate of decay. It is also assumed that the “zero” current is a fixed percentage of the peak current (equivalently a certain number of time constants of the exponential rate of decay), so that meaningful comparisons of the “peak-to-zero” time for one lightning strike with the same time for another lightning strike can be made. That is, it is assumed that “zero” current=KI o , where K is a fixed value which is (significantly) less than 1.0. 
     One conservative way to establish the level of “zero” current reported by a lightning detection network would be, in view of all the lightning strikes for which “peak-to-zero” data are reported, to define the “zero” current as being the smallest value of current that remains detectable above the background noise at all the detection stations used for reporting the data. That would ensure that an actual “zero” current was actually detected and measured to establish the “peak-to-zero” time for every reported lightning strike. An estimate of the “zero” current using this rationale is to take the “zero” current as being half the smallest of the reported peak currents. 
       FIG. 5  shows a histogram of reported values of peak current for 204 lightning strikes in a test area of about 32,000 m 2  (measurements taken from the year 2000 through 2011). The largest peak current is 140 kiloamps and the smallest peak current is about 1 kiloamp. Taking half of this latter value as “zero” current, Equation (9) can be solved for each of these lighting strikes. 
       FIG. 6  shows the computed resistance values, which, discarding obvious outliers, can be seen to fall broadly within a range 1.0 and 2.5 kilohms, with the significant majority of the values falling more narrowly between about 1.3+/−0.3 kilohms. Thus it is apparent that a usefully precise measure of ground resistance R in the test area can be obtained by use of Equation (9). 
     Resistivity ρ is related to resistance R as ρ=RA/L. Here, A/L defines the geometry of the volume of earth (“VOE”) into which the lightning current flows, and specifies a size of the VOE in linear dimensions (e.g., meters). As an approximation, the charged area of the earth is a circle with a radius proportional to the elevation of the cloud  12  ( FIG. 1 ), and the depth of penetration of the current is proportional to the radius of this circle. Thus the area of the volume of earth through which the current flows is proportional to the square of the elevation of the cloud, and the length of the volume of earth through which the current flows is proportional to the elevation of the cloud. Taking the elevation of the cloud as an unknown constant,  FIG. 7  shows a histogram of resistivity for the R values shown in  FIG. 6  based on this approximation. 
     It should be recognized that, while Equation (9) relates R to peak current I o , it is not essential to use peak current in the computation. It can be observed from Equations (2) and (3) more generally that: 
     
       
         
           
             
               R 
               = 
               
                 t 
                 
                   C 
                   [ 
                   
                     ln 
                     ( 
                     
                       
                         I 
                         ⁡ 
                         
                           ( 
                           
                             t 
                             = 
                             
                               t 
                               1 
                             
                           
                           ) 
                         
                       
                       - 
                       
                         ln 
                         ( 
                         
                           I 
                           ⁡ 
                           
                             ( 
                             
                               t 
                               = 
                               
                                 t 
                                 2 
                               
                             
                             ) 
                           
                         
                         ] 
                       
                     
                   
                 
               
             
             , 
             
               
                 for 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 any 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 values 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 of 
                 ⁢ 
                 
                     
                 
                 ⁢ 
                 
                   t 
                   2 
                 
               
               &gt; 
               
                 
                   t 
                   1 
                 
                 . 
               
             
           
         
       
     
     It follows that an equivalent to Equation (9) can be obtained and used to solve for R for any two values of current, together with the corresponding times, or equivalently, a first current at an initial time t 1 , and the rate of decay of the current which establishes the current at any subsequent time t 2 . 
     Computations according to the principles described herein are preferably performed on a computer; either a general purpose digital computer or computer system suitably programmed to perform such computations, such as a PC or Macintosh, or a special purpose computer suitably hardwired, or a combination of the two. Data for the computations may be downloaded from the lightning databases over a local area network, or over a wide area network such as the Internet, or they may be provided to the computer or computer system in computer readable media, such as hard disks and CD-ROMs. 
     It is to be understood that, while a specific method for determining geological subsurface resistivity has been shown and described as preferred, other methods in accordance with the principles described herein could be utilized without departing from the principles of the invention. 
     The terms and expressions which have been employed in the foregoing specification are used therein as terms of description and not of limitation, and there is no intention in the use of such terms and expressions to exclude equivalents of the features shown and described or portions thereof, it being recognized that the scope of the invention is defined and limited only by the claims which follow.
A method for determining geological subsurface resistivity. The method includes obtaining a set of lightning parameters associated with a lighting strike received by a geological volume of material, the set of lightning parameters including an indicium of the current of the lightning strike at a first initial time and an indicium of the current of the lightning strike at a first decay time subsequent to the first initial time, and inferring the resistance of the volume of geological material, at least in part, from the set of lightning parameters.