Natural Hazards and Earth System Sciences The Pulse Azimuth effect as seen in induction coil magnetometers located in California and Peru 2007 – 2010 , and its possible association with earthquakes

The QuakeFinder network of magnetometers has recorded geomagnetic field activity in California since 2000. Established as an effort to follow up observations of ULF activity reported from before and after the M = 7.1 Loma Prieta earthquake in 1989 by Stanford University, the QuakeFinder network has over 50 sites, fifteen of which are high-resolution QF1005 and QF1007 systems. Pairs of highresolution sites have also been installed in Peru and Taiwan. Increases in pulse activity preceding nearby seismic events are followed by decreases in activity afterwards in the three cases that are discussed here. In addition, longer term data is shown, revealing a rich signal structure not previously known in QuakeFinder data, or by many other authors who have reported on pre-seismic ULF phenomena. These pulses occur as separate ensembles, with demonstrable repeatability and uniqueness across a number of properties such as waveform, angle of arrival, amplitude, and duration. Yet they appear to arrive with exponentially distributed inter-arrival times, which indicates a Poisson process rather than a periodic, i.e., stationary process. These pulses were observed using three-axis induction coil magnetometers that are buried 1–2 m under the surface of the Earth. Our sites use a Nyquist frequency of 16 Hertz (25 Hertz for the new QF1007 units), and they record these pulses at amplitudes from 0.1 to 20 nano-Tesla with durations of 0.1 to 12 s. They are predominantly unipolar pulses, which may imply charge migration, and they are stronger in the two horizontal (north-south and east-west) channels than they are in the vertical channels. Pulses have been seen to occur in bursts lasting many hours. The pulses have large Correspondence to: J. C. Dunson (cdunson@quakefinder.com) amplitudes and study of the three-axis data shows that the amplitude ratios of the pulses taken from pairs of orthogonal coils is stable across the bursts, suggesting a similar source. This paper presents three instances of increases in pulse activity in the 30 days prior to an earthquake, which are each followed by steep declines after the event. The pulses are shown, methods of detecting the pulses and calculating their azimuths is developed and discussed, and then the paper is closed with a brief look at future work.


Introduction
Earthquake prediction has long been a controversial topic, and over 100 yr of developing seismological science has yet to produce a practical earthquake forecasting methodology.Occasionally, an observation is made prominent in the scientific community or perhaps in the public News Media, Fraser-Smith et al. (1990) being a notable example.But rebuttals and vigorous debate also remain, Karakelian et al. (2002), Kappler et al. (2010), and Thomas et al. (2010) being notable recent examples.
The consensus stated during a recent EMSEV panel discussion in 2010 was that evidence for the existence of earthquake pre-seismic signals of several sorts is strong, yet hopes for constructing a reliable forecasting network based on these signals in the near future remain dim.Difficulties cited include: (1) the widely variant geological character of the various seismic regions around the world; (2) the sparseness of available measurements; and (3) the wild variability of the pre-seismic signals themselves.Therefore methods that can produce more systematic results out of the sparse pre-seismic signals are appropriate to investigate.QuakeFinder has been funded by government agencies and private donations, and has developed a network of magnetometer sites in California, Peru, and Taiwan, as described by Cutler et al. (2008) and depicted in Fig. 1.
The sensitivity of the Zonge ANT4 coil used in the highresolutions systems is among the best for coils in its class (see Table 1).
The discovery of the pulses in QuakeFinder data was first reported in Bleier et al. (2009).The number of pulses found in our data since this original discovery have steadily increased as we have refined our tools that identify them, but topics such as the statistics of the pulse properties, the character of the bursts, the pulse dispersion characteristics, statistics of their angles of arrival, multi-site observations, etc. have received little comprehensive study as yet.
Our literature search has found scant mention of pulses and even less analysis of them, with typical papers such as Kotsarenko et al. (2004) and Alperovich et al. (2003) mentioning pulses as simply a feature of the data.We find that the majority of work has sought to eliminate pulses (also called spikes) from the analysis, and that many common ULF analysis techniques achieve this result almost by definition.Several of these techniques include: the M.A. indices bank of bandpass filters assays bulk energy changes as in Bernardi et al. (1988), Fraser-Smith et al. (1990), Villante et al. (2009); the Intra-Station Transfer Function, which uses long-term averaging as in Harada et al. (2004) and Izutsu et al. (2007); and the eigenvalue of the Principal Component Analysis method that seeks to find periodic signals in the data as per Gotoh et al. (2002), Serita et al. (2005), Morrison et al. (2004), Haykawa et al. (2007).As far as we could find, there is very little in the way of analysis of pulses in magnetometer data that goes beyond such discussions as found in Geomagnetic Pulsations by J. A. Jacobs, or Whistlers and Related Ionospheric Phenomena by R. A. Helliwell.
Moreover, these well-known "micropulsations" of the "sudden commencement" and "sudden impulse" types have been ruled out as the source of the pulses reported in this paper.Pulses here are only seen in one out of 2-6 local sites, whereas standard "micropulsations" are world-wide and essentially simultaneous.The QuakeFinder network has observed a number of these global events, but pulses reported here are local, not "extra-terrestrial", and often are very powerful.

Detection model
"Why all this sudden interest in pulses?"It is a question we have heard often.Aside from the fact that discovery of a confirmed and repeatable seismic precursor signal would be very important, the shorter duration of the pulses, their appearance in bursts, and in narrow angles of arrival at a three-axis sensor enable additional discriminators to be applied, improving the design of the overall QuakeFinder detection experiment.
The design of the detection experiment can be seen in the application of Bayes' Theorem of Probability.This basic relation can be used to describe the behavior of many detectors, as is shown in Eq. (1): where: -E is the occurrence of an earthquake event.
-D is the occurrence of a detector event.
-P (E) is the probability that an earthquake event will occur that matches our analysis parameters.P (E) varies regionally around the Earth, depending on proximity of the observations to a seismically active area.
-P (D) is the probability that the detector will produce an event, earthquake or not.
-P (D|E) is the conditional probability that a detector event will occur given that an earthquake event occurs.
-P (E|D) is the result, and is the probability of an earthquake occurring given that a detector event has occurred.
Trying to predict rare events such as seismic activity using sensor data requires much caution against false positive results, as does the study of any rare phenomenon such as a medical test for a rare disease.Confidence in the results can take a long time to develop.The process can be thought of as being divided into two steps: (1) first develop confidence in the test (or detector) to prove that it accurately detects the phenomenon; and (2) develop a statistically significant means of using the test (detector) to find more of the phenomenon.
In the case of earthquake pre-cursor signals, little confidence exists that a reliable pre-cursor signal has yet been found.The method here describes an unconfirmed earthquake pre-cursor signal (pulse azimuth clusters), and how a detector might be constructed that reduces the occurrence of false positive results.

Highly discriminative detectors
To see now the importance of minimizing false positive results, suppose for the purpose of argument that we state the results of one run of the analysis are: 1.If an earthquake happens, the probability that the detector will work is 90 % (0.9), i.e., it produces events for nine out of ten earthquakes.
2. But during the control periods of no earthquakes, the detector returns false positives 2 % of the time.
This might lead us to believe that the overall false positive rate is close to 2 %.But that is quite incorrect for the case of rare events such as earthquakes.False positives can easily dominate the answer.Suppose that we assume a probability of 1 big quake every three years in California as stated during a recent oral presentation by Robert Nadeau in his presentation at AGU2009.Let's call that 1/1000 days, or a probability of ∼0.001 per day.We then can write: P (E|D) = 0.9 × 0.001 0.9 × 0.001 + 0.02 × 0.999 This result means that the probability that the detector is reporting a false positive that is 1-0.0429 or almost 96 %! Some have concluded from this insight that making such a detector is intrinsically impossible.But suppose that we add a discriminator such that the detector false alarms drop to 0.2 %.Repeating the calculation, we now get 0.316 and the probability of a false positive is reduced to 69 %.And thus we arrive at a fundamental principal for earthquake forecasting work: A primary focus must be on developing detectors that have strong discrimination characteristics.
Contrasted with a detector that simply seeks anomalous increases in ULF activity with half-hour averages over a period of days, a pulse detector can be much more discriminative, meaning that the quantity of experimental control data compared to the signals of interest can grow more steeply due to the short nature of the signal.The Pulse Azimuth Cluster method also helps here, as it adds another layer of discrimination to the detection.

Method
The pulse azimuth technique uses pairs of time series to measure the peak amplitude ratio (azimuth) of each pulse along orthogonal axes.Classically, the definition of "azimuth" is "the horizontal direction of a compass bearing", but the peak amplitude ratio detector described here is also effective in the two vertical planes: north/vertical; and east/vertical.It is not a measure of the angle from the site to the epicenter or hypocenter, but rather a discriminative way to evaluate the Nat.Hazards Earth Syst.Sci., 11, 2085Sci., 11, -2105Sci., 11, , 2011 www.nat-hazards-earth-syst-sci.net/11/2085/2011/pulses.Consider the amplitudes of the large pulse shown in Fig. 3.These are the raw waveforms of one of the pulses identified by Bleier et al. (2009), and is taken from the vertical and east/west channels of the East Milpitas site.The pulse amplitudes (top panel) are treated as zero-centered by removing the mean from the unsigned counts produced by the 24 bit ADC, and then by using a high pass filter as described in Sect.2.2 below.The amplitudes are large enough here that the nominal geomagnetic signal is scarcely visible in the plot.
An important concept is that the pulse waveforms are unipolar, meaning that they are not like many other magnetometer signals that oscillate and decay around a zero center.The unipolar nature and large amplitude of these pulses allows one to go beyond simple peak amplitude ratios, and instead take a point-by-point measure of the amplitude ratio over the length of the pulse.
The lower panel of Fig. 3 will be described now.Because of the zero centering, the amplitudes of a pair of signals can be treated as having four possible orientations with respect to this zero (x/y): +/+; +/-; -/+; and -/-.Therefore, the top panel is an instance of a negative (-) pulse for the east/west, and a positive (+) pulse for the vertical, hence -/+.There are three such combinations of orthogonal axes, assigned with ordering that follows a right-hand rule sequence: (1) x = north, y = east; (2) x = east, y = down; and (3) x = down, y = north.
Since the ADC samples all channels simultaneously (not multiplexed), we can make use of the quadrant decomposition mechanism of the atan2() function of the ANSI C-library for every pair of samples in the pulse.Comparing this with the standard result, the red line is the atan() bulk estimate of the signless, amplitude squared, energy signal, it is so close here to the blue histogram that the difference is barely visible (the dotted red line at −160 is its 180 degree alias).Examples of analyses using the direction finding approach include Hayakawa et al. (2007) and Izutsu et al. (2007).
Other researchers have treated the vertical channel using the Polarization Ratio concept, e.g., Hayawaka et al. (2007), where the local signal is expected to produce more pronounced increases in the vertical coil than in the horizontal coils.In both of the site histories shown in this paper, however, despite possible association with the steep angle to a hypocenter, more energy appeared in the horizontal coils than in the vertical.
We are reporting occurrences of dozens of consecutive pulse events that arrive in bursts from a given azimuth angle, close to the horizon.
The process of computing the pulse azimuths is now outlined to make the above points clear.There are five steps, also shown as a flow diagram in Fig. 33: 1.Data Exclusion (block out bad data)

Data exclusion
QuakeFinder uses a database approach for masking out known bad data events, though for this paper all data exclusion was done painstakingly by hand.Electrically incorrect conditions were not found in either site's data used in this paper, though there was an outage channel caused by a bulldozer severing our north/south coil at Alum Rock that lasted almost three months.Some known cultural noise, as per Bleier et al. (2009), was also removed from consideration, resulting in the table of data exclusions (Table 2).
Another form of exclusion sometimes used by QuakeFinder is to mask out periods when low frequency S waves from distant earthquakes are shaking the coil.
The Zonge induction coil magnetometer is highly sensitive, and through physical shaking, they can detect S-waves at 0.02 Hz, 0.005 Hz, and ever lower.Fortunately, though these very low frequency signals are often quite strong, their frequency content is well below the high-pass filter cutoff described in the next section.So S-wave exclusion is not warranted in the case of these pulses.

Mean removal and highpass filter
A high-pass filter was applied to reduce the effect that the powerful Pc3-4 pulsations can have on the signal.Originally much attention was paid to a proper mean removal scheme   2006-07-21 08:27:00 2006-07-21 10:08:00  all  CONTAMINATED DATA  2006-08-03 11:10:00 2006-08-03 11:35:00  all  CONTAMINATED DATA  2006-08-04 06:16:00 2006-08-04 07:02:00  all  CONTAMINATED DATA  2006-09-21 20:10:00 2006-09-21 20:15:00  all  CONTAMINATED DATA  2006-11-05 19:20:00 2006-11-05 20:00:00  all  CONTAMINATED DATA  2006-11-25 22:16:00 2006-11-25 22:30: since the 24 bit ADC produces unsigned counts, as shown in Cutler (2008), to ensure that the azimuth directions were realistic.But two realizations now dissuade this approach.The first is that the azimuth angles are not affected very much for the typical mean variations if the signal-to-noise ratio for the pulse exceeds 10:1, which it always does given the thresholds that we have chosen.The second realization is that seeking an absolute azimuth to an epicenter from a set of induction coils is misguided, and that the pulse azimuth method is not in fact direction finding (see below).The focus here is on finding increases of pulse counts from any narrow (∼5 to 30 degrees) range of "angle." So for this step we first convert the unsigned ADC counts into signed, zero-centered values.This is accomplished by selecting a 24 h (86 400 s) period to search for pulses, and then by loading 36 800 s both from the day before and the day after to create a 160 000 s time series.The average is then computed and subtracted, producing signed counts.
Next we use a high-pass filter on this 160 000 second time series that eliminates longer period signals, such as the powerful Pc3-4 geomagnetic pulsations, and also removes the mean.We selected a scaled FIR filter of order 1373 with a cutoff frequency of 0.1 Hz.Its response curve is shown in Fig. 4.

Pulse detection
Once the high-pass filter has been applied, a pair of thresholds is then used to identify pulse events in the zero-centered time series data.A set of parameters has been chosen that allows the pulse detector design to be varied over a significant range of cases.These cases are represented in the standard way as matrices, where each column contains a certain parameter, and each row is another case study.Figure 5 shows a single detected pulse and the parameters.
These parameters are now summarized: Detector type: We have experimented with three types of pulse detectors and used the Threshold method (Thr), which simply applies a straight amplitude threshold to the signal for this paper.
Detector name: This is a single word, used to form all directories and file names unique to this detector.Threshold values: These are pairs of numbers that form a range to be used to control the actions of the pulse detector.There are three pairs, one pair for each axis.The exact interpretation of these values is determined from the selected type of pulse detector to be considered.If a pulse crosses the lower threshold, it is considered a "down" pulse; if it crosses the upper threshold it is considered an "up" pulse.For the detectors discussed in this paper, the two limit numbers are  3, and are in units of counts (24 bit ADC range is 0-16777215) Duration range: These values specify the minimum and maximum length of a pulse that will be considered.The maximum duration is particularly important, as it is used to block continuous wave trains.The pulse detector developed by QuakeFinder is very conservative in that once it detects a pulse excursion, it not only must be back within limits before the maximum pulse length, but there can be no more excursions for another whole maximum pulse length period, before or after the detected pulse, or the entire pulse train is discarded.To get a feel for how much effect this distinction makes, we can note that Bleier et al. (2009) reported anomalous pulses at a rate of 150 per day without any duration or wavetrain discrimination, while the second-generation detector described here peaks at 130 pulses.So it is not a very big difference.
A complete set of these parameters is considered to be the definition of a pulse detector.These have been set to a number of different experimental values, and each time a run is completed, a new listing of pulses is generated.The pulse detectors used for this analysis are listed here in Table 4. Here, an inverse relationship exists in that using a lower threshold identifies more pulses, but it also throws out the very large pulses that slightly overshoot zero and cross the other threshold before returning to zero.If an unrelated noise peak breaks the threshold during the rejection window, an important pulse can be tossed out.Figures 6 and 7 illustrate these issues with the pulse detector: Figure 6 shows the "exp5" detector used on the sensor history of the East Milpitas, CA, site.It shows increases in pulse counts per day against two earthquakes near the site, while Fig. 6 shows the effect of using the tightest detector, "exp6", on the same data set.Note that overall, the pulse    counts are much higher for the tighter threshold, but that the signal looks noisier because the threshold is low enough to cause two effects: (1) more of the regular signal variations exceed the threshold; and (2) big pulses are eliminated more often by a secondary pulse crossing the threshold during the exclusion window, and being excluded.
A recent paper by Bortnik et al. (2010) has proposed modeling these pulse currents as linear flows, and it went on to calculate the magnitude of a current that would produce 2-25 nT pulses, similar to those reported for Alum Rock.Based on this result, we estimate that if pulses reported here are subsurface current flows, the amount of current is on the order of 1-1000 KA.

Azimuth histogram computation
Now we have excluded bad data, applied a high-pass filter, and used the pulse detector to produce a listing of pulses.The detected pulse times can now be used to load the time series for the azimuth computation.

Not direction finding
Restating explicitly, there are two degrees of freedom in the angle of arrival of any signal caused by a local current flow, whether it be lightning, cultural, or seismogenic.They are: (1) the direction from the site to the location of the current flow; and (2) the direction of the current vector at that location.The azimuth method here seeks to identify repeatable angles of arrival whether they point at a hypocenter or not. Figure 8 illustrates the point: Bulk direction finding has been typically limited to a range of 90 degrees due to the fact that squaring the amplitude and computing a ratio of energy in two axes is in effect a signless operation.The traditional formula for such an operation is effectively: 90-atan(|y/x|).It is important to emphasize again that the technique used in this paper is not in fact true direction finding, but rather an algebraic mapping of amplitude ratios onto a circle.
The signals originate from a three-axis set of induction coils buried beneath the surface using the Cartesian right-hand coordinate system on the surface of the Earth: north = X, east = Y , and down = Z.For reasons of familiarity, the azimuth method uses standard compass heading as shown in Fig. 9.
Each pairing of samples within the detected pulse is converted into an angle via the relation: azimuth = atan2 (NS, EW), rounding the angle to the nearest degree, and binning them into 360 one degree-wide bins.These histograms are then treated as circular, where degree 180 wraps around to degree −179.Then each azimuth histogram can be plotted onto a larger plot as a single point to produce a plot such as Fig. 27, where each dot is a separate histogram.Some choices for how to plot each histogram are to plot the peak histogram bin, the median angle value, or the mean angle value.All plots in this paper use the mean value of the histogram.

Azimuth histogram stability
As will be shown below, the SNR of the pulse affects the stability of the resultant azimuth estimate.High SNR pulses, like those at Alum Rock, give a more repeatable azimuth computation than do weaker pulses where the noise dominates more.For the pulses at Alum Rock, each burst of pulses were almost all within an 8 degree window (see Fig. 16 below), whereas in the Tacna case, they have a wider angular variance, on the order of 30 to 45 degrees (Fig. 20, typical).Here are several questions about the properties of these azimuth estimates: 1.How much of the histogram variance is driven by noise? 2. How much is caused by the signal processing? 3. How much is driven by real differences in the waves arriving from the source?
We are not yet far enough along to address these questions very well, but we have simulated by: (1) computing Signalto-Noise-Ratio (SNR), and (2) adding noise to real signals.
A simple model of the relationship of the strength of these pulses vs. the spread of the azimuth histogram can be made for pulses like the one shown in Fig. 3, where we find that the voltage signal to noise ratio is often in excess of 5000 for both axes.That is obviously very good.To get a feel for how SNR affects the stability of the azimuth histogram, we artificially degenerate one of the pulses detected at Alum Rock by whitening (one component is 10× greater than the other, i.e.: atan(0.1)= ∼ 6 degrees), we can produce Fig. 10.We can see that if the weaker signal of the pair has an SNR > 10, the change in angle is very small, i.e., the azimuth method produces a stable estimate.

Azimuth cluster evaluations
The term "Pulse Azimuth Cluster" is used to denote a series of pulses that are detected by a pair of orthogonal sensors, where each pulse projects a fixed ratio of amplitude onto the pair of orthogonal sensors.They are called "clusters" because they appear as bursts over a period of hours or days, and when mapped onto a circle, they are confined in angle between 5 and 45 degrees, (so far).One of the strengths of the method is that the azimuth angle repeats even when successive pulses have widely variant magnitudes.
There are a number of ways to evaluate the properties of these Pulse Azimuth Clusters.We have experimented with summarizing each histogram using mean, median, and peak as stated in Sect.2.4.1 above, the goal of which is to reduce the evaluation of a pulse to specific single angle value, and plot them against time, as is shown in Figs. 19,20,27,28,29,and 31.Our method so far has been mostly to plot and select these angle values using a GUI tool.One quantitative technique we have explored is to bin the angle values into one-degree bins, and then compute a probability density function using the history of the site.This technique can be used to compare seismic and non-seismic time periods.

Results, part one: the pulses
As we discuss results, three earthquakes are also discussed.The earthquake facts are given in Table 5.The relative significance of these earthquakes can be viewed against a background of other earthquakes witnessed by the QuakeFinder network as is shown in Fig. 11.The data sets are largely complete for both earthquakes except for an outage of the north/south coil at the Alum Rock site, which lasted from 28 July 2007 to 14 October 2007.We show results from three sites.Table 6 shows the site parameters.
Now we present the overall pulse counts from the full history at three sites for each channel in Fig. 12.In Fig. 12, you can see the high variability of pulse counts that are produced by using fixed pulse thresholds across multiple sites.The 704 site is obviously quieter than the 702 site.The threshold limits of the exp6 detector (too tight) strip the 702 and 704 pulse counts more than they do the 609 (described above in Fig. 4).And finally, despite the proximity of the BART electric railroad system to the 609 sensor, it has very low overall pulse counts over its five years of operation as compared to the 14 and 6 months of the Peruvian systems.The threshold we use exceeds the nominal BART amplitude discussed in Bleier (2009).
Underscoring what was previously stated, there remains much work to be done in the area of refining the pulse detectors.The "exp5" detector, for example, reveals compelling pre-seismic activity for the 609 East Milpitas site (shown in Fig. 6 above), while the exp8 detector is similarly interesting for the 704 Tacna, Peru, Ch2 History, Fig. 13.
Characterizing these pulses, even with these crude detectors, is an interesting exercise.First we show some typical Other pulse parameters that are of interest are the magnitude, duration, and inter-arrival time.The plots of the pulse magnitude effects take two forms, as is shown in Figs. 17 and  18.
The magnitude for Figs. 17 and 18 is computed using Eq. ( 3 (3) where: S = the amplitude each sample contained within the pulse; and µ = the long-term 160 000 s channel mean described in Sect.2.2 above.
Figure 18 shows the result of averaging the pulse energy per day over the cumulative number of pulses in that day.
More results of a similar nature can be found in the remaining data sets not shown here.

Results, part two: the Azimuth Clusters
Prior to the magnitude = 5.4 earthquake at Alum Rock, a series of pulses was observed that originated from a similar direction at an observation station two kilometers from the epicenter as reported in Bleier et al. (2009).This case shows that the pulses cluster about certain angular locations.See Figs. 27, 28, and 31 to see typical azimuth cluster results.
Figure 19 is a zoom-in of a typical Alum Rock cluster, taken from the period between the time that the north/south coil was repaired and the earthquake occurred.This plot deserves careful study.First note that the entire cluster is confined to a seven degree window.If it is lightning, then it is spread over a period of days, and is originating from some sort of stationary storm system, i.e., it is very likely not lightning.Another interesting feature of the plot is the angular spread.The distance to the hypo-center is approximately 9.5 km, so if these pulses are caused by a series of subsurface current flows, then the signal variance naively implies repeated flows over a ∼1.2 km region.
For the Tacna case, the pulse clusters are substantially different.Figure 20 is one of the nine clusters that were manually identified:     This is the next-to-last cluster before the earthquake in Fig. 28, zoomed-in and trimmed.Given the shorter pulse durations at Tacna, greater distance to the quake, and larger angular variance, we might be tempted to conclude that we are dealing with a different phenomenon altogether.Except that reviews of other QuakeFinder sites reveal similar signal structure, on other commensurate scales (e.g., Fig. 31) In the case of these Tacna pulses, using similar arguments as Alum Rock, we see that all of the clusters (Fig. 28) span multi-day periods, and again we can argue that if it were lightning, it would have originated from a non-moving storm system.If we propose a near-hypocenter source, we see that for ∼46 km to the hypocenter, activity would be occurring over a ∼35 km region.Analysis of the vertical plane azimuths at Tacna reveal a similar clustering near the horizon as shown for the Alum Rock case, meaning larger amplitudes in the horizontal coils than in the vertical.Note now that we have manually divided the data set into two groups.We have pulses that were identified as being part of a Pulse Azimuth Cluster (using the GUI tool), and those pulses that were not.This distinction allows us to contrast the pulse durations at Tacna 704 to see that they are longer than typical, just as Bleier et al. (2009) reported for Alum Rock in Fig. 21.
Going back now to the beginning, we show the first plot we made of the Pulse Azimuth Cluster phenomenon as Fig. 26.This plot was presented by Tom Bleier at AGU 2009, and it is a raw ratio of axis amplitudes made before the atan2() method was employed.The plot in Fig. 27   number of our sensors at present.Little analysis has been done yet on these signals.The reason again that the East Milpitas plots focus on the Vertical and East/West channels is that the North/South coil was severed by a bulldozer and was out for over two months just before the earthquake.Likewise in the 14 days before the M = 6.2 earthquake near the Tacna site, the incidence of Azimuth Clusters increased.This is shown in Fig. 28.an important pitfall to avoid when building the pulse detector).The inset in Fig. 29 shows the large positive pulses that were mistakenly classified as down pulses.But hand detailed analysis has proven that all four of the blue cluster groups in Fig. 29 are detector errors, and are in fact up pulses.
When investigating the properties of these pulses, there are several simple treatments which can offer insight.Typical repeating or man-made signals tend to vary around a mean value, so a plot of a property such as pulse duration would have some sort of central value, and a standard deviation around that value.According to the Central Limit Theorem, that curve would approach the typical bell curve, given enough of the signals.Likewise, if we were to plot the times between a series of periodic events (inter-arrival time), it too would have an almost constant value, with small differences due to noise.
A Poisson Process, on the other hand, is one in which events occur randomly and act as if independent of each other.Given an interval of time, if events occur at random times (aperiodically) during the interval, then the amount of time between events is a random number.If this random number is exponentially distributed, see Montgomery (2007) for more detail, then we are dealing with a Poisson Process.The standard earthquake rate model is an example of a wellknown Poisson Process, as are many other natural phenomena.
Figure 22 shows a comparison of the pulse inter-arrival times taken from the Alum Rock data set.Clearly the plot does not indicate a mean value and the spread around that value as it would for a periodic signal.What it does show looks to be typical data for a Poisson Process.The green line is a preview of what might be a proper fit if further hypothesis testing is carried out in detail.Figure 23 -Pulse Durations at Alum Rock.This time, the green estimate has some qualitatively picked numbers.The reason they are given rather than a '?' question mark as above, is that we are inviting explanations about this result from the community.Fig. 23.Pulse durations at Alum Rock.This time, the green estimate has some qualitatively picked numbers.The reason they are given rather than a "?" question mark as above is that we are inviting explanations about this result from the community.
The significance of the fact that pulse durations behave similarly is somewhat less clear as is shown in Fig. 23.These durations are also not centered around a mean value as they would if they were variations of a repeated phenomena.The pulse cluster amplitude distribution is also interesting.Due to the pronounced asymmetry of the distribution's tails, it is unlikely that the amplitudes will ultimately be shown to have Gaussian properties.One good fit we found was with the Weibull distribution.
The significance of amplitudes that vary in this manner is not known.The Weibull distribution has been used to model dielectric breakdown and conduction in oxides and polymers as in Costa et al. (2006)and Dissado and Fothergill (1992).It has also been widely used for general failure rate modeling.If the pulses are subsurface current flows of a finite length (Bortnik et al., 2010) that occur during higher stress periods approaching mechanical failure (Freund, 2007a), then this might be an appropriate distribution for modeling them.

Discussion
The pulse azimuth detection method described in this paper lacks many desirable features that a fully developed system might have.But it is a first step, and as shown here, it does produce intriguing results -despite the fact that the thresholds have not been finely trimmed separately for the signal strengths at each site.Fortunately, the time required to regenerate pulse sets for the entire ∼500 GB of network data  history has reduced a hundred-fold since the first attempts, and now each run takes only a matter of days using the parallel processing machines of QuakeFinder's Data Center.
The pulse detectors used for this paper identified increases in the counts of uni-polar pulses occurring prior to all three quakes discussed.The significance of the pulses being uni-polar has been the subject of much discussion.Many signals seen in our magnetometers oscillate multiple times before decaying, whereas the uni-polar pulse, by definition, does not.One type of phenomenon that can cause a uni-polar pulses in a coil magnetometer is the migration of charge from one location to another.Lightning is an example of charge moving from cloud to ground or vice-versa.Currents flowing underground, such as those suggested by Freund (2007b), could also produce uni-polar pulses.
These pulses appear to have exponentially distributed inter-arrival times, which suggests that the occurrences may follow the random rate of arrival of a Poisson process, rather than a periodic process.They have similar shapes, and in every case take more than five data points to rise fully (when  sample rate is 32 Hz), and thus they are all band-limited below the Nyquist frequency of our system, meaning that they may be fairly well observed with the QF1005 & QF1007 instruments.
The pulses exhibit dispersion characteristics (not whistler mode), and have the classic look of a sharp impulse that has passed through a low pass filter, where the higher frequencies are delayed more.It will be an interesting exercise to try to determine the "filter" parameters that the media hypothetically applied.Such a result could possibly be compared to the effects that various soil conductivities are known to have from magneto-telluric data as one means to gain insight into whether the pulses are indeed propagating underground.
If the signals are originating from subsurface near a hypocenter, then the idea of evanescent propagation like in Grimalsky et al. (2003), or an analog to the seismic "head wave" discussed in Stein and Wysession (2003) that occurs at interfaces between layers of differing indices of refraction, might be a useful framework for discussing them.
The fact that the pulses do not appear in multiple sites argues against them originating in the ionosphere.A simple comparison of our two closest sites (∼43 km), as seen in   2009) result, as they found that the pulses at Alum Rock did not appear in any of the other California sites.Likewise, inspection of the El Carmen site in Peru (∼800 km away) show no occurrence of the pulses seen in Tacna, Peru, data.Additionally a detailed study was done using two Southern California sites to assess the common occurrence of pulses in the data.Despite the appearance of a number of confirmed uni-polar pulses in our site at Ocotillo Wells (Fig. 31), few were seen in our Julian site 43 km away (Fig. 30).
Similarly, the planetary K p Indices do not show high activity on the dates in question, as is shown in Fig. 32.
We have used lightning data from time to time to try and verify that it is not the source of the signals we are seeing, but  the pulses we are discussing here are very likely not lightning as Bleier et al. (2009) have discussed.In the case of Alum Rock, they argue that the duration of the pulses (mostly exceeding one second) exceed that of the many lightning pulses that we have examined over the years by an order of magnitude.Lightning pulses also rise faster than our sampling rate (32 or 50 samples per second) can depict; pulses here have slow rise times, possibly exhibiting dispersion, and have very little content above 10 Hertz.Unlike Alum Rock with its 2-10 s pulse durations, the case of Tacna shows many of the pulses that are ∼0.2-5s in length, which does more closely match lighting pulse timing.
But residents in Southern Peru tell us that there is no lightning during that time of year and there are other factors that discourage the lightning explanation, among them the fact that that the pulse activity can be shown to significantly drop at the time of the earthquake.
The manmade noise hypothesis is similarly troubled, as the levels of power seen in the pulses is very large, typically 8-20 nT (Bleier et al., 2009), and this level of signal is extremely difficult to produce with man-made equipment at an installed site.Bleier et al. depict  magnetometer using all of the nearby electrical equipment, including a grinder and welder, could not create enough of a signal to exceed the pulse detector's threshold at 0.85 nT by even a single sample.Aside from being very close to the coils, there are few sources that can generate such large energy, and even fewer man-made signals that are not stationary.
For the large signals reported in this paper, the pulse azimuth method produces stable results.But there are limitations, when the SNR is less than 20, that must be taken into account.
As incoherent noise increases in the channels, the estimate varies more.As coherent noise increases in the channels, the system will tend towards giving azimuth estimates at 45 degrees.(X = Y ), (X =-Y ), etc.
Azimuth values near the poles (0, π/2, π, 3π/2 etc.) can be less stable.To understand this idea, think of a pulse that projects almost completely onto one axis.The other axis has almost zero energy in it, thus its contribution to the azimuth estimate is mostly noise.Fortunately all three of the above conditions are easy to detect, though a formal investigation into these limitations has not yet been completed.
Because pulses are clustered in time and angle, repeated occurrences can be treated as a series of Bernoulli trials, where the probability of the next pulse appearing within a window of angle near the preceding pulses can be computed empirically based on the incidence of all pulses in the history of the channel/site.
Each arriving pulse has a probability that it either will arrive in a designated angle-of-arrival window, or it will not.The other outcome is the complementary result: P (in) = Probability next pulse is in window (4) P (out) = Probability not in window And : P (in) + P (out) = 1 The probability that any two successive pulses will arrive in the window is much smaller than for just one (P (in) × P (in)).More generally, the probability that n out of m pulses will arrive in the Azimuth Window can be modeled using the binomial distribution.Experiments with this method so far are promising.One last detail worth mentioning is that the pulses at Alum Rock also exhibit curious spectral characteristics as can be seen in the spectrogram Fig. 25.This is a subselection of pulses detected that day made to highlight the effect.While the small knees outlined appear at near Schumann frequencies, they are narrower in bandwidth and much sharper in peak energy than the Schumann resonance: The list of questions opened up by these observations includes: -Are the pulses related to seismic activity?-Why do pulse azimuth clusters occur in bursts?-What is their observable range?-Does a second site allow location triangulation?-What is their speed of propagation?-Do pulses scale with earthquakes?

Future priorities
The first and foremost priority is to encourage others in Japan, Greece, Italy, et al. who may have similar data sets to try the Pulse Azimuth Cluster method on their data.
A second priority is to pursue the properties of all of the similar pulses found in QuakeFinder data.This includes addressing questions such as: 1. Are their time of day variances?2. Do variances scale with range and magnitude of nearby quakes?
3. Do the various Quakefinder sites report pulses in a manner consistent with the local conditions at the site (conductivity, etc.)?
A third future priority is to use portable systems and trucks to quickly add additional observations to fixed sites that are showing pulse activity before the impending and possibly associated earthquake occurs.In this way, the pulses may be able to be triangulated, allowing the possibilities for a source, and for the speed of propagation to be narrowed down significantly.To do so requires a more portable data and power system, both of which are under intensive development at QuakeFinder.Another initiative is to improve the pulse detectors and also the daily automated monitoring that can alert us to the occurrence of them much more promptly.Pulses can be detected without using the Fourier transform and also without doing system response compensation, and thus are wellsuited for a typical embedded micro-controller to identify and report in near real time.
Since the pulses are very large and have short range, smaller, less expensive coils can allow for denser deployments, as each site costs less.QuakeFinder has developed a new, less expensive coil, the "QFIDO" to evaluate and deploy new, lower cost systems.

Fig. 1 .
Fig. 1.The Calfornia Network of QuakeFinder magnetometer sites (CalMagNet).The blue triangles are the original High School (QF1000) projects sites installed mostly in 2000 and 2001.The green triangles are second generation QF1003 sites, while the gold triangles are the high resolution QF1005 and QF1007 sites, which were installed starting in 2005.The two cases presented in this paper contain data from the East Milpitas site, which can be seen immediately east of the San Francisco Bay Area.

Figure 2 - 3 Fig. 2 .
Figure 2 -The two QuakeFinder sites in Peru.The observations here took place at the Southermost site -Tacna.

Figure 3 -
Figure 3 -Typical Alum Rock Pulse Azimuth Histogram.The top panel is the waveform, and the bottom contrasts typical direction finding results (red) with the histogram result (blue)

Fig. 3 .
Fig. 3. Typical Alum Rock pulse azimuth histogram.The top panel is the waveform, and the bottom contrasts typical direction finding results (red) with the histogram result (blue).

Figure 3 -Fig. 4 .
Figure 3 -Typical Alum Rock Pulse Azimuth Histogram.The top panel is the waveform, and the bottom contrasts typical

Figure 5 -
Figure 5 -Parameters used in pulse detection.The red + and dashed line are the threshold level of pulse detection.T0 is the Start time of the pulse and T1 is the end time.TP is the time of the pulse peak, while P is the peak amplitude.

Fig. 5 .
Fig. 5. Parameters used in pulse detection.The red + and dashed line are the threshold level of pulse detection.T 0 is the Start time of the pulse and T 1 is the end time.T P is the time of the pulse peak, while P is the peak amplitude.

Figure 5 -
Figure 5 -Parameters used in pulse detection.The red + and dashed line are the threshold level of pulse detection.T0 is the Start time of the pulse and T1 is the end time.TP is the time of the pulse peak, while P is the peak amplitude.

Figure 6 -
Figure 6 -Pulse Counts with nominal limits, each point being the pulse count for one day.

Figure 7 -
Figure 7 -Pulse Counts with too-tight limits

Fig. 6 .
Fig. 6.Pulse counts with nominal limits, each point being the pulse count for one day.

Fig. 8 .
Fig. 8. Model of the Pulse Azimuth Effect.The green segments show relative projected amplitude.If φ, the angle to the location of the current flow remains unchanged; but if theta, the angle of the field near that current flow site changes, and then the ratio of energy projected on the coils will change.

Fig. 10 .
Fig. 10.Degeneration of azimuth estimate.As the amount of added noise increases, SNR drops, and the estimate degrades.

Fig. 11 .
Fig. 11.The QuakeFinder experiment, the big picture.All earthquakes witnessed by each and every site are shown here as blue markers, giving range vs. magnitude.The three earthquakes discussed here are circled in red.As more sensors are installed and network density increases, the likelihood of an earthquake happening near a sensor increases.

1Figure 12 - 4 Fig. 12 .
Figure 12 -Overall Pulse Count Results for the five pulse detectors.The det

Figure 12 -
Figure 12 -Overall Pulse Count Results for the five pulse detectors.The detectors were named in the order they were created, so the names are not reflective of the limits.SeeTable 4, above.

Fig. 21 .
Fig. 21.Tacna pulse durations, showing the difference between those identified as being part of a cluster vs. all others.

Figure 21 -Fig. 22 .
Figure 21 -Tacna Pulse Durations, showing the difference between those identified as being part of a cluster vs. all others.

Figure 22 -
Figure 22 -Inter-Arrival Times look exponentially distributed, and may indicate a Poisson Process.(The geometric distribution is simply a discrete analog of the exponential distribution) The green question marks are left in place to indicate that this is an initial result only.
Figure 28 -Pulse History for CMN 704 Tacna Peru Site, clusters identified using GUI tool. 3 4

Fig. 30 . 2 Figure 31 -Fig. 31 .
Fig. 30.Arrival time overlay From Julian, Ca & Ocotillo Wells.All horizontal pulses at both sites were compounded into lists, and these two lists were compared to find timing overlaps.Only 74 out of 14 995 overlap, mostly on specific dates.

Fig. 27 ,
Fig.27, shows that only 74 out of 14 995 pulses analyzed from both horizontal channels have any overlap in time at all.This result echoes theBleier et al. (2009) result, as they found that the pulses at Alum Rock did not appear in any of the other California sites.Likewise, inspection of the El Carmen site in Peru (∼800 km away) show no occurrence of the pulses seen in Tacna, Peru, data.Additionally a detailed study was done using two Southern California sites to assess

Figure 32 -
Figure 32 -Kp Value Plots.The top panel is the Kp history shown for the period discussed in this paper.(Max Kp in each day) The bottom panel is a zoom in to the time of the Tacna earthquake.

Fig. 32 .
Fig. 32.K p value plots.The top panel is the K p history shown for the period discussed in this paper (max K p in each day).The bottom panel is a zoom into the time of the Tacna earthquake.

5 Figure 33 -
Figure 32 -Kp Value Plots.The top panel is the Kp history shown for the period disc

Table 2 .
Table of data exclusions used during this analysis.All were at the East Milpitas site.The Tacna, Peru, site had no outages.
the same.Typical values are shown in Table

Table 5 .
Earthquakes discussed in this paper.

Table 6 .
Detector sites used (the 4th and 5th sites are mentioned only at the end of the paper, but are shown here.)