Estimating the seismotelluric current required for observable electromagnetic ground signals

We use a relatively simple model of an underground current source co-located with the earthquake hypocenter to estimate the magnitude of the seismotelluric current required to produce observable ground signatures. The Alum Rock earthquake of 31 October 2007, is used as an archetype of a typical California earthquake, and the effects of varying the ground conductivity and length of the current element are examined. Results show that for an observed 30 nT pulse at 1 Hz, the expected seismotelluric current magnitudes fall in the range ∼10–100 kA. By setting the detectability threshold to 1 pT, we show that even when large values of ground conductivity are assumed, magnetic signals are readily detectable within a range of 30 km from the epicenter. When typical values of ground conductivity are assumed, the minimum current required to produce an observable signal within a 30 km range was found to be ∼1 kA, which is a surprisingly low value. Furthermore, we show that deep nulls in the signal power develop in the non-cardinal directions relative to the orientation of the source current, indicating that a magnetometer station located in those regions may not observe a signal even though it is well within the detectable range. This result underscores the importance of using a network of magnetometers when searching for preseismic electromagnetic signals.


Introduction
Over the past few decades there have been a number of observations of ultra-low frequency (ULF: 10 mHz-10 Hz) magnetic signal anomalies preceding large earthquakes.A few notable examples include: (i) the M s = 7.1 Loma Prieta earthquake that occurred on 17 October 1989, and was preceded by roughly 2 weeks of anomalously high magnetic activity that continued (and increased) until the main shock (Fraser-Smith et al., 1990;Bernardi et al., 1991); (ii) the M s = 6.9 Spitak, Armenia earthquake that occurred on 7 December 1988, and was preceded by several hours of anomalous activity (Molchanov et al., 1992;Kopytenko et al., 1993), and (iii) the M s = 8.0 Guam earthquake that occurred on 8 August 1993, and was preceded by 2 weeks of anomalous activity (Hayakawa et al., 1996).Owing to the great devastation inflicted by large earthquakes (e.g., Bilham, 2010), it would be immensely useful to be able to forecast the approach of a large earthquake on a timescale of hours to days.However, observational repeatability with ULF signals has proven to be a challenging task (Fraser-Smith et al., 1994;Karakelian et al., 2002), and the unique identification of the precursory ULF signals remains a controversial issue due to the plethora of natural signals and any artificial instrumental noise that often confound such observations (Campbell, 2009;Thomas et al., 2009a, b).
In an effort to improve the dependability of individual instrument observations, Bleier et al. (2009) combined data from a number of different instruments in observing the 31 October 2007, M w = 5.6 earthquake that occurred in Alum Rock Park, in San Jose, California (henceforth referred to as the "Alum Rock" earthquake).In that study, Bleier et al. used  ULF range, located ∼ 2 km from the epicenter, to record a series of magnetic pulses with peak magnitudes of a few 10's of nT whose occurrence frequency peaked roughly 2 weeks prior to the main shock.In addition, an air conductivity sensor co-located with the magnetometer showed that air conductivity saturated for ∼ 14 h in the night and morning hours prior to the main shock, and infrared satellite imagery showed anomalously high activity in the vicinity of the epicenter, peaking ∼ 2 weeks prior to the main shock.
Several physical mechanisms have been proposed for the generation of large, underground electrical currents, that may account for the reported observations.These include the electrokinetic and magnetohydrodynamic effects resulting from fluid flowing through rocks, piezomagnetism, stressinduced variations in crustal conductivity, microfracturing, and so on (Draganov et al., 1991;Park et al., 1993;Park, 1996;Fenoglio et al., 1993Fenoglio et al., , 1995;;Johnston, 1997;Merzer and Klemperer, 1997;Molchanov andHayakawa, 1995, 1998;Molchanov et al., 2001;Egbert, 2002;Surkov et al., 2003;Simpson and Taflove, 2005).Recent efforts stemming from rock experiments performed under laboratory conditions, have attempted to unite the assortment of phenomena reported prior to large earthquakes into a single mechanism that originates from the semiconductor-like behavior of pressurized igneous and high-grade metamorphic rocks (e.g., Freund, 2007a, b;Freund et al., 2009).
Regardless of the actual physical mechanism involved in producing the current, the work of Bleier et al. (2009), and other potential ULF precursor observations prompt the question: what is the minimum electrical current necessary to produce an observable magnetic signal on the ground, at a given distance from the epicenter and for an assumed ground conductivity?It is the aim of the present study to answer this question.We use the parameters of the Alum Rock earthquake as an archetypal model of a California earthquake, and assume that the magnetic signal observed on the ground is generated by a small current element co-located with the earthquake hypocenter, and having a dimension which is roughly comparable to the scale size as the rupture zone.We use a set of closed-form solutions to analyze an idealized representation of the problem, which is discussed in Sect. 2 and given in Appendix A for reference.A similar technique has been used in previous work to treat related problems such as the fields generated by electric dipoles, radiating both underground and under the ocean (Dong et al., 2005;Fraser-Smith et al., 1979;Tian and Mata, 1996).The work of Molchanov et al. (1995) studied the closely-related problem of a cylindrically-symmetric, distributed underground current source, which shows excellent agreement with the present work, and is discussed in greater detail in Sect. 4. In Sect. 3 we discuss the choice of our assumed parameter values, present the results for the Alum Rock earthquake case, and show the general distribution of the 3-D fields created by the seismotelluric current element.We provide a discussion and conclusions in Sects.4 and 5, respectively.

Description of the method
In order to study the the electromagnetic wave fields emitted by the seismotelluric current, we use the theoretical approach developed by King et al. (1981) for an antenna lying near a planar interface, and represent the situation using a simple idealized model as shown in Fig. 1a.In this picture, an infinitesimally short, horizontal dipole is located at a depth d in the half-space z > 0, filled by medium 1, which is characterized by its electrical properties µ 1 (magnetic permeability in Henry/m), ε 1 (permittivity in Farad/m), and σ 1 (conductivity in mho/m), all of which are real, and represent the partially-conducting Earth.Medium 2 is the remaining half-space z < 0, which represents air, and is characterized by the electrical properties µ 2 (= µ 0 = 4π × 10 −7 Farad/m), ε 2 (= ε 0 = 8.854 × 10 −12 Farad/m), and σ 2 (=0).If we define the complex permittivity of the medium as ε = ε + iσ/ω (where ω = 2πf is the operating frequency of the dipole), Maxwell's equations in medium j can be written: (1) where j = 1,2 for z > 0 and z < 0, respectively.The volume current density J is assumed to lie in the x-direction, at a depth z = d, and is normalized to have unit electric moment (I l = 1), so that: As a final constraint, the boundary conditions at z = 0 require the continuity of the wave components E x , E y , εE z , µ −1 B x , µ −1 B y , and B z .
In order to solve the system of Eqs.(1-3), the 2dimensional spatial Fourier transform is taken of E as follows: and similarly for B and J .Defining the complex wave number k 2 j = ω 2 µ j εj , and its z-projection γ 2 j = k 2 j − ξ 2 − η 2 where j = 1,2, and imposing the boundary conditions listed above, a set of explicit solutions for the six electromagnetic wave components can be obtained throughout the entire space after a fair amount of non-trivial manipulation.The final set of equations used in the present paper is given in Appendix A, using the same notation as King et al. (1981) (who gives a detailed derivation of the equations).
We have chosen to use this simple, homogeneous model of the earth for a number of reasons: firstly, it enables us to obtain fairly straightforward, closed-form solutions for the fields that only involve a double-integral, and allow us to rapidly explore the effects of a number of geophysical parameters and obtain 'order of magnitude' estimates for the seismotelluric current.Secondly, due to the low frequency of the wave (f = 1 Hz throughout this study), the radiated wave is not affected by small-scale geological features, but rather feels the effects of the average properties of the medium, on the scale size of a wavelength or skin-depth (which, in this case is on the order of the depth of the current element).Finally, assuming a particular geological structure (or equivalently, a conductivity profile) would inevitably constrain our solution to a particular region, and we prefer to keep the solutions as general as possible.

Simulation results
Using the set of equations in Appendix A , we aim to relate the seismotelluric current magnitude at the earthquake hypocenter, to the electromagnetic fields observed at a given point in the space.In order to do so, we must constrain some of the free parameters of our model, those being: (i) the seismotelluric current magnitude, I , (ii) it's depth d, and (iii) scale size l, (iv) the characteristic radiation frequency of the current element, f , (v) the observation location ( x, y, z), and (vi) the Earth's average electrical parameters, (µ 1 ,ε 1 ,σ 1 ).Below, we discuss our approach to constraining these free parameters, and present the ensuing results of our analysis.

Parameter selection
Of the set of free parameters listed above, by far the most elusive is the ground conductivity σ 1 (henceforth referred to as σ , since the air conductivity σ 2 = 0), since it depends not only on the local petrology, but also on its distribution of porosity, temperature, and pressure (e.g., Wait, 1966).Combined with f , the conductivity determines the skin-depth, given by the approximate relation δ = (πf µ 0 σ ) −1/2 , which is the scale length beyond which the radiated electromagnetic wave will begin to be strongly attenuated by the medium.Typical skin depths are shown in Fig. 2 over a 5-order of magnitude variation of σ and 2-order of magnitude variation in f , and vary between ∼ 0.15 km and ∼ 500 km, which is a considerable range.Based on the typical pulse lengths of ∼ 1 s reported by Bleier et al. (2009), we set f = 1 Hz (solid line in Fig. 2) in the analysis performed in this study, with the understanding that our results can be scaled as ∼ e −z/δ for different radiated frequencies and conductivites.The approximate conductivity ranges for some of the typical elements that constitute the crust are shown in Fig. 3, based on Palacky (1993, Fig. 1).In general, the older igneous and metamorphic rocks form the bulk of the crust and are relatively poor conductors, σ ∼ 10 −6 − 10 −3 mho/m.Above this rock stratum is a surface layer of sedimentary rocks that has a wide variation in conductivity, and is also fairly porous and uncompressed allowing for the infusion of mineral-rich ground water (e.g., "Archie's law", Winsauer et al., 1952).This upper layer has typical conductivites in the range ∼ 0.01 − 0.1 mho/m, and extends to a depth of a few hundred meters.Thus, its depth represents only a few percent of the total distance from the earth's surface to the earthquake nucleation depth (∼ 10 km), and its averaged effect on the conductivity is minimal.A survey of representative conductivities in the California region (Lat ∼ 36 • ) has shown that typical conductivities are σ < 10 −3 mho/m down to the depth of the Moho (i.e., within our region of interest, d < 30 km), beyond which conductivities again begin to increase (Hermance, 1995, p. 207).
The other parameter values have a much smaller effect on the attenuation of the radiated fields.We set µ 1 = µ 0 , i.e., assume the earth to be non-magnetic, and ε 1 = 5ε 0 , after typical values of rock permittivitties (e.g., Krauss, 1991, p. 134).Based on these assumed electrical properties of the earth, we analyze the effects of the radiating seismotelluric current element below.

Alum Rock example
We begin by analyzing the situation described by Bleier et al. (2009).We set d = 10 km, and the observation point at ( x = 0, y = 2 km, z = 0) as illustrated schematically in Fig. 1b, where the observing magnetometer station is labeled "Unit 609".Figure 3 shows the observed magnetic field B y for a dipole with unit moment, i.e., I l = 1 A-m as a function of σ , where it is evident that at low σ , the signal passes through the earth virtually unattenuated, whereas for high σ , there is a precipitous drop in power over several orders of magnitude, and virtually no signal can be observed at the surface.The value of σ at which the transition occurs is σ δ = 2.5×10 −3 mho/m, which corresponds to δ = 10 km = d, i.e., the seismotelluric current element is at precisely one skin-depth below the earth's surface.This underscores the importance of the skin-depth as a critical scale length of the system.In fact, we note that the effects of source-depth can be analogously thought of as simply varying the ratio d/δ, which can be inferred from different conductivity profiles using Fig. 2.
Using the results shown in Fig. 3, it is now possible to quantify the approximate seismotelluric current magnitude, by assuming that the observed signal B obs y = 30 nT, following the typical pulse amplitudes observed by Bleier et al. (2009).The moment I l is scaled such that the observed magnetic field (y-axis of Fig. 3) is always equal to B obs y , which results in a σ -dependent value of I l that is relatively flat, but rises sharply when σ > σ δ .Furthermore, since the value of I l is now constrained, we can assume a certain scale size for the seismotelluric current-element, and obtain the required current I .The results of the analysis described above are shown in Fig. 4, where σ is again varied over the entire range of typical crustal conductivities, and the scale length l is varied in the range 0.01 km-10 km, which is consistent with typical rupture lengths of moderate earthquakes scaled by a "safety" factor of ∼ 5−10 (e.g., Bonilla et al., 1984).As expected, the lowest values of I (lower right corner) occur when the rupture length is largest, and σ is smallest, requiring only a few kA to produce the observed 30 nT pulses.On the other hand, for σ > 0.1 mho/m, I becomes unrealistically large, well in excess of ∼ 10 6 A. Focusing on the typical expected range of values, i.e., σ ∼ 10 −3 − 10 −2 mho/m, and l ∼ 0.1 − 1 km (dashed oval), the required current magnitudes fall in the range of a few 10's of kA to a few 100's of kA.We note that the results shown in Fig. 4 scale linearly with I l, so that source-current magnitudes for any other value of B obs y can be simply inferred (e.g., for B obs y = 3 nT, divide Fig. 4 currents by a factor of 10).

Detectability
We now fix the intensity of the radiating element to a constant value, I l = 10 8 A-m (e.g., 100 kA current flowing over a 1 km scale-length), and examine the variation of the observed wave on the earth's surface as a function of lateral distance ( y) from the epicenter.The results are shown in Fig. 5, parameterized by σ , and indicate again that for σ < σ δ , the radiated wave remains largely unaffected by the conducting earth, and decays by roughly a factor of 10 over a 30 km distance.On the other hand, for σ > σ δ , the signal strength is quickly attenuated as a function of distance from the hypocenter.
It is possible to set a lower limit on the seismotelluric current magnitude that can be detected by assuming a certain detectability threshold, say 1 pT in this instance.Using the data in Fig. 5, it is seen that within a 30 km radius about the epicenter, signals will be observed even when the earth is extremely conductive, σ ∼ 0.1 mho/m when I l = 10 8 Am.Since the observed B y scales linearly with I l, and I l = 10 8 A-m produces B y > 10 −9 T within 30 km of the epicenter (for σ ∼ σ δ ), we conclude that for the threshold B thresh y = 10 −12 T, I l = 10 5 A-m.Since rupture length decreases with M, assuming l ∼ 0.1 km results in I ∼ 1 kA, which is a surprisingly low value.

Wave distribution
In order to visualize the distribution of wave power surrounding the seismotelluric current element, we performed the calculation of wave fields for d = 10 km, σ = σ δ , I l = 10 8 Am in the volumetric region | x| < 50 km, | y| < 50 km, −15 < z < 30 km (z is defined positive in the downward direction, as above), and show the magnitude of B y in Fig. 6a.Here we show horizontal slice planes spaced every 3 km, from a depth of 30 km to an altitude of 15 km, and as expected, |B y | peaks near the depth of current element (10 km), and decays as a function of distance from the source.The magnitude of B y is roughly 1 order of magnitude smaller at the surface, than at 10 km, and is strongest at the epicenter.Figures 6b, and  and |E x | (i.e., the strongest field components from an xdirected current element) on the Earth's surface, and indicate that the fields decay by roughly 1 order of magnitude (×10) at ∼ 20 km, and 2 orders of magnitude at ∼ 40 km from the source in the x-direction.However, observations made in non-cardinal directions relative to the orientation of the current-element can fall into deep nulls in the radiated power, where the wave amplitude is attenuated by over 4 orders of magnitude, and thus may not be detected even though they are well within the detectable region.This result of our analysis underscores the importance of using a distributed network of magnetometers when attempting to detect any seismic-related magnetic signal.In a similar vein, we reiterate that the present model assumes a simple homogeneous ground, and any deviations from this assumption (e.g., introduction of small-scale, or structured conductive or magnetic regions) would introduce structure in the magnetic field distribution observed on the ground, in the form of nulls or quiet-regions.This may potentially explain why a few nearby stations did not observe the same signal as station 609 in Bleier et al. (2009).

Discussion
Although the parameter values used in Sect. 3 are modeled after the Alum Rock earthquake, the results are nevertheless more general than they appear and can be used to estimate the expected magnetic signals in other regions.For example, the well-known M w = 6.0 earthquake that occurred on 28 September 2004 in the heavily-instrumented Parkfield, California region produced no observable magnetic signals on  the ground prior to the main shock (Johnston et al., 2006), and possibly no precursory activity of any kind (Bakun et al., 2005).In the context of the present study, one possible reason for the lack of observed precursory magnetic signals is the large value of the ground conductivity in the Parkfield region.Based on magnetotelluric surveys, Unsworth and Bedrosian (2004) estimate that the conductivity may reach values as large as σ ∼ 0.3 − 0.03 mho/m, to a depth of several km near the epicenter, which reduces the skin depth to a range of δ ∼ 0.8 − 3 km at f = 1 Hz (Fig. 2).To crudely estimate the expected magnetic signal at the University of California at Berkeley's PKD electromagnetic observatory (lateral separation of ∼ 20 km from the epicenter), we use Fig. 6a which is calculated for a I l = 10 8 A-m at 1δ, and attenuate by a factor ∼ exp(−d/δ) where d = 7.9 km is the depth of the earthquake, and δ is the skin-depth calculated above, which gives a range of 0.07 − 5 × 10 −5 .The resulting magnetic signal is B y ∼ 0.2 nT −0.15 pT, which straddles the threshold of detectability, and would very likely fall in the noise level and not be detected.The expected signal levels directly above the hypocenter would fall in the range ∼ 7 nT −5 pT, and would very likely be detected by an induction magnetometer (threshold ∼ 1 pT), but not a proton precession magnetometer (threshold ∼ 0.25 nT).This example serves to illustrate the order-of-magnitude magnetic field values might be expected, but of course, the reader is reminded that the current element intensity I l is set at at an arbitrary level in this example, and there is no guarantee that the true source would radiate at this intensity.
The assumption of a point-source radiating as an infinitesimal current element, which has been used throughout this study, deserves some discussion.Although this is certainly a simplified model of the actual situation, it is nevertheless representative of the power reaching the earth's surface, in that any realistic current distribution within the ground can be modeled as a multipole expansion of terms, each of which decays (roughly) inversely proportional to distance, raised to the power of poles.Thus, even if the quadrupole and octupole terms are included, they will decay significantly faster than the lowest order dipole term.Secondly, the real source current is not infinitessimally small but has some finite size which introduces "edge effect" corrections in the vicinity of the dipole.These edge effects have been neglected in our analysis since they significantly complicate the analysis, and are small, so long as the current element lengths are much smaller than the system scale length; in our case the system scale size is on the order of 10 km, whereas the antenna size is ∼ 1 km, or smaller.Finally, the signals reported by Bleier et al. (2009) were time-limited pulses rather than a continuous wave, as we analyzed presently.In order to model an arbitrary time waveform, it can be expanded into its Fourier series, each frequency component propagated from source to observation point as in the present study, and the frequency components summed to produce the resulting signal at the surface.However, the skin depth diminishes rapidly with increasing frequency (e.g., Fig. 2), and thus the dominant frequency components that will propagate to the surface with the least attenuation, will be the lowest frequency component generated by the source, which we have taken to be f = 1 Hz.Thus, performing the analysis on only the lowest frequency component as we have done above, is expected to be representative of the majority of the wave power reaching the observing equipment on the earth's surface.
For illustrative purposes related to the use of a pointsource radiator in the present work, we compare the closelyrelated study of Molchanov et al. (1995) who analyzed the effects of a cylindrically-symmetric, distributed current source flowing respectively in the vertical, azimuthal, and radial directions.Although an identical comparison cannot be made due to a number of differing assumptions, we can roughly scale the current moment of Molchanov et al. (1995), Fig. 3, to ∼ 10 8 A-m, which gives a magnetic field intensity on the ground at a 10 km displacement, of a few nT for the radial and azimuthal current sources (i.e., most similar to our source geometry).This is compared with our value of a few 10's of nT (in the broadside y-direction) or a few nT (in the x-direction) at a 10 km lateral displacement (e.g., Fig. 5), for conductivity values of σ ∼ 10 −4 − 10 −2 mho/m (i.e., deep to shallow σ values in Ibid).This comparison shows very similar results, when considering that our source was located at d = 10 km, as opposed to d = 20 km (Ibid), and is thus expected to be roughly an order of magnitude stronger.Since the aim of the present work was to establish crude "orderof-magnitude" estimates for the field-strength that can be expected due to an underground current source of a prescribed intensity, as well as the expected scaling due to variations in ground conductivity, distance, and current-length, it is comforting to know that the precise details of the source geometry do not significantly alter the field-strength produced by a given current moment.

Conclusions
In this paper, we used a relatively simple model of an underground current source co-located with the earthquake hypocenter to roughly infer the magnitude of the seismotelluric current required to produce observable ground signatures.We used the geometry of the Alum Rock earthquake (Bleier et al., 2009) as an archetype of a typical California earthquake, with a nucleation depth of ∼ 10 km, and varied the ground conductivity and length of the current element.Results show that for a 30 nT pulse at 1 Hz, the expected current magnitudes for the most typical conditions fall in the range ∼ 10 − 100 kA, which is smaller than previous estimates.
Setting the detectability threshold to 1 pT, we showed that even for large values of ground conductivity, the magnetic signal produced by the seismotelluric current would be detectable within a range of 30 km from the epicenter.For typical values of ground conductivity, the minimum current required to produce an observable signal within a 30 km range was found to be ∼ 1 kA, which is a surprisingly low value.
When the magnetic signal was displayed as a 2-D distribution within a 50 km range of the epicenter, our analysis showed that deep nulls in the signal power develop in the non-cardinal directions relative to the orientation of the source current, indicating that a magnetometer station located in those regions may not observe a signal even though it is well within the detectable range.This result underscores the importance of using a network of magnetometers when searching for preseismic electromagnetic signals.
a triaxial induction magnetometer operating in the Published by Copernicus Publications on behalf of the European Geosciences Union.

Fig. 1 .
Fig. 1.Schematic representation of the geometry used in the theoretical model.(a) An x-directed current element is placed at x = 0,y = 0, and depth d, where positive z is defined positive in the downward direciton.(b) Geometry of the Alum Rock earthquake, d = 10 km, and the location of the observing station (Unit 609, after Bleier et al., 2009) displaced by 2 km from the epicenter.

Fig. 2 .
Fig. 2. Skin depth shown as a function of typical crustal conductivities, for radiated waves in the ULF range.

Fig. 3 .
Fig. 3.The calculated value of B w expected at the observation location (Fig. 1b) due to a source of I l = 1 A-m.The conductivity ranges of typical crustal materials are superimposed for reference, as well as the "skin-depth" conductivity, σ δ .

Fig. 4 .
Fig. 4. Current magnitude required to produce a 30 nT field at the observation location, shown as a function of the bulk conductivity of the earth, and estimated length of the current element, l.The dashed region in the center of the figure indicates the most likely values of σ and l, and thus the typical current magnitudes that are most probable.

Fig. 5 .
Fig. 5. Magnetic field amplitude |B y | observed on the earth's surface, due to 10 8 A-m current element, parameterized by bulk conductivity.The curve marked σ δ represents the value of σ = 2.5 × 10 −3 mho/m which gives a skin depth equal to the depth of the current element.

Fig. 6 .
Fig. 6.Spatial distribution of radiated wave fields due to a I l = 10 8 A-m current element at d = 10 km.(a) 3-D distribution of B y , showing horizontal slices at 3 km increments from a depth of 30 km to an altitude of 15 km, (b) 2-D distribution of B y at ground level (z = 0) and (c) similarly for E x .