An analytical expression for early electromagnetic signals generated by impulsive line-currents in conductive Earth crust, with numerical examples

Changes in the electromagnetic (EM) field after an earthquake rupture but before the arrival of seismic waves (“early EM signals”) have sometimes been reported. Quantitative evaluations are necessary to clarify whether the observed phenomena are accounted for by known theories and to assess whether the phenomenon can be applied to earthquake early warning. Therefore, analytical expressions for the magnetic field generated by an impulsive line-current are derived for a conductive half-space model, and for a two-layer model; the somewhat simpler situation of a conductive whole-space is also considered. By analyzing the expressions obtained for the generated EM field, some expected features of the early EM signals are discussed. First, I verify that an early EM signal arrives before the seismic waves unless conductivity is relatively high. Second, I show that early EM signals are well approximated by the whole-space model when the source is near the ground surface, but not when it is at depth. Third, I show that the expected amplitudes of early EM signals are within the detection limits of commonly used EM sensors, provided that ground conductivity is not very high and that the source current is sufficiently intense. However, this does not mean that the EM signals are easily distinguishable, because detector sensitivity does not account for additive noise or false positive detections.


Introduction
Seismic fault slip mainly generates ground motion, but electromagnetic (EM) variations have been reported during and prior to some earthquakes [e.g., Eleman 1966, Iyemori et al. 1996, Honkura et al. 2002, Abdul Azeez et al. 2009, Honkura et al. 2009, Okubo et al. 2011, Gao et al. 2014]. The term "early EM signal" means any measurable seismogenic EM variation that precedes the first seismic phase arrival. A few reports include early EM signals that appear related to earthquakes. For example, Iyemori et al. [1996] observed changes in the geomagnetic field at the time of the 1995 Mw 7.1 Hyogo-ken Nanbu (Kobe) earthquake at two stations nearly 100 km from the source region. The EM variations had maximum amplitudes of 0.6-1.0 nT, and began ~10 s before the Pwave arrival (~10 s after the estimated seismic origin time). For the 2007 Mw 6.9 Iwate-Miyagi Nairiku earthquake, Okubo et al. [2011] observed changes in the geomagnetic field ~3 seconds before the arrival of seismic waves at a hypocentral distance of 27 km; the maximum signal amplitude was ~0.2 nT. The reliability of other previously reported early EM signals is disputed, as these studies did not sufficiently discuss discrimination between signal and background geomagnetic fluctuations (i.e., false positive detections). However, if early EM signals can be reliably observed, then they could be incorporated into early earthquake warning systems, which currently rely solely on seismic observations [e.g., Kanamori 2005].
Various mechanisms have been suggested to explain the conversion of mechanical energy to EM energy. Some examples include the motional induction effect [e.g., Eleman 1966, Gershenzon et al. 1993, the electrokinetic effect [e.g., Ishido andMizutani 1981, Pride 1994], the piezoelectric effect [e.g., Utada 2000a,b, Huang 2002], and the piezomagnetic effect [e.g., Yamazaki 2011a,b]. To assess the detectability of early EM signals, it is essential to estimate expected amplitudes and to determine the distinguishing signal characteristics (if any) associated with each mechanism. The most straightforward way to study early EM signals is, of course, to study the measured EM fluctuations themselves; however, strictly observational approaches are difficult because large earthquakes and sensitive geomagnetic recording stations are both rare, meaning the undisputed data set is small. Most observations of early or coseismic EM signals have come from temporary sites designed for EM prospecting, including magnetotelluric (MT) instruments [e.g., Honkura et al. 2002, Abdul Azeez et al. 2009, Gao et al. 2014. Successful detections of early EM signals are rather rare at continuous stations; among the only clear examples are Iyemori et al. [1996] and Okubo et al. [2011]. One potential way to capture more early EM signals is to deploy instruments near the epicenter of a large, shallow earthquake, because many aftershocks are expected [e.g., Ujihara et al. 2004, Kuriki et al. 2011]. However, because aftershocks are necessarily smaller than the main shock, and early EM signals related to the main shock are already near the limits of conventional detectors, aftershocks are difficult to detect in practice. Consequently, it is difficult to compile a suitably large data set of early EM signals. It is therefore expedient to use theoretical models of EM variations to guide observational efforts.
In the present work, subsurface electric currents induced by ground motions are examined as possible early EM source mechanisms. The situation is schematically shown in Figure 1. During an earthquake rupture, electric currents can be induced by the motion of conductive crustal material in the ambient geomagnetic field [e.g., Gershenzon et al. 1993, Yamazaki 2012, Gao et al. 2014]. The objective is to determine the expected temporal variations in the generated EM field. An electric current in a vacuum generates a magnetic field that can be easily calculated using the Biot-Savart law. Some earlier numerical calculations were performed by assuming that the medium was a vacuum or a perfect resister [e.g., Taira et al. 2009, Okubo et al. 2011]. However, an electrical current in a conductive medium should be dealt with by solving the time-dependent Maxwell's equations, to determine, for example, if the conductive medium acts as an insulator against EM propagation.
The reminder of this paper is organized as follows. In Section 2, the governing equations are described for specific situations. In Section 3, a set of analytical expressions is derived for the EM field generated by an impulsive line-current in a conductive half space. In Section 4, some characteristics of early EM signals in a conductive medium are discussed, along with some numerical examples. The validity of using a simpler model for evaluating early EM signals is also considered. Finally, the main conclusions are presented in Section 5.

Definition of the problem
In this paper, I derive analytical expressions for the early EM signals arising in the two model geometries: a two-layer model comprising a half-space with uniform conductivity, representing the Earth's crust, overlain by a perfect resistor, representing the air ( Figure 2a) and a whole-space model consisting of a uniformly conductive infinite medium (Figure 2b), which is the simplest system that can model the effect of electrical conductivity. A line current is considered an approximation of the electric current induced by the coseismic motion of fault-zone material, which is often more conductive than the surrounding crust, as illustrated in Figure 1. Considering the situation whereby the fault length is orders of magnitude greater than the spatial scale of conductivity variation across the fault, it is reasonable to assume that the electric current occurs as a line of infinite length, thus reducing the problem to two dimensions.
The half-space model (Figure 2a) is the simplest realistic approximation of the Earth. A realistic model must involve at least two layers: a conductive ground and resistive air. In contrast, the whole-space model ( Figure 2b) is somewhat unrealistic; however, solving Maxwell's equations in the whole-space model is considerably easier than for the half-space model. For this reason, I also consider the whole-space model throughout the following, so that any situations in which the simplified model is valid can be understood.
When the electric current density I is given as a function of location (x, z) and time t, the electric and magnetic fields, E and B respectively, are determined via Maxwell's equations: YAMAZAKI Figure 1. A schematic illustration of the subsurface-current generating mechanism analyzed in this paper. Seismic motion of a fault block containing a zone with high conductivity subject to the ambient geomagnetic field induces an electric current according to Faraday's law. where f and n are, respectively, the electric and magnetic permeabilities.
It is assumed that n is constant such that n =n 0 = 4r×10 -7 H/m. It is also assumed that the electric current density is the sum of an induced current density and a conductive current density; the conductive current density obeys the isotropic Ohm's law (3) where v is the electrical conductivity, and I ext is the external current density. Additionally, the displacement current term is ignored; that is, the second term of Equation (2). This approximation yields (4) Introducing a vector potential that satisfies B =∇×A, Equations (1) and (4) are combined to form a single expression The two-dimensionality of the problem allows a coordinate system to be chosen such that the y-axis is oriented parallel to the external electric current, and the positive z-direction is vertically downward (see Figure 2). In this coordinate system, only the y-component of I ext is non-zero; the y-component of I ext is denoted as I y ext . The electric current density of an impulsive linecurrent at depth d is therefore given by (6) where d is the Dirac delta function. In this situation, we can chose a solution of Equation (5) of the form A = (0, A,0), with the function A satisfying the scalar equation (7) The EM field is then determined by solving Equation (7) under appropriate boundary conditions. The first boundary condition comes from the EM field being zero at infinite distance from the source current: i.e., A(x,z → ±∞, t)=const. When considering the two-layer model, conditions on the continuity of A and ∂A/∂z are also imposed at the boundary z = 0 to ensure the continuity of E and B at any instance. The non-zero components of the electric and magnetic fields are then given by E y = ∂A/∂t, B x = −∂A/∂z, and B z = ∂A/∂x.

Derivation of solutions
The governing equations can be simplified by the substitution of regularized variables t * ≡n 0 vt, x * ≡n 0 vx, z * ≡ n 0 vz , and d * ≡ n 0 vd. Using these variables, the source term is written as (8) The differential operators become ∂/∂t=n 0 v∂/∂t * and ∇= n 0 v∇ * . The governing equation including the source term can then be rewritten as (9) in a conductive medium (i.e., v >0), where I 0 * ≡n 0 2 vI 0 , and in a vacuum (i.e., v =0).
To find the solution, consider the integral expression and the corresponding Fourier expression of A (12) In these identities, g * and ~* represent the wavenumber and angular frequency, respectively. Substituting these expressions into Equations (9) and (10), the governing equations for A (g * , z * ,~*) are obtained as (13) where Solutions to Equations (13) and (14) are given as .
) ) ) ) ) ) ) ) ) ) ) ) ) ) ) and (2) follows. In the case of the whole-space model ( Figure  2b), which is described solely by Equation (13), the solution is given by In the case of the two-layer model (Figure 2a), which is described by Equation (13) for z>0 and by Equation (14) for z<0, upward-and downward-propagating terms should be added to Equation (16) [e.g., Stoyer 1977]. These additional terms are determined such that they satisfy the boundary conditions, and the following solution is obtained: Observations of variations in the EM field are usually obtained near to the ground surface (i.e., z ≅ 0). Because of the continuity of the EM field at ground level (z=0), it is sufficient to consider A either for z >0 or for z < 0. For this reason, only the time-domain solution for A at z >0 is considered.
To evaluate the integral in Equation (12) over ~, I refer to the formulas for the inverse Fourier transforms and where u is defined by Equation (15), and g is a positive real number. The function erfc represents the complementary error function, The result shown in Equation (18) can be derived by changing the integration variable from ~* to u and choosing an appropriate integral path over which to apply the Cauchy integral theorem (see Figure 3). Equation (19) can then be derived from Equation (18) after some algebraic manipulation. In the following, I consider phenomena only for t >0 because, from Equations (18) and (19), it is clear that the solution is trivial at negative times.
To evaluate the integration of Equation (12) over the variable g * , I use the easily derived results formulas: Here, erfi is the imaginary error function defined by Equation (21) is derived by integration by parts using the identity: Using Equations (18), (19), (21) and (22), the integration of Equation (12) using the function given by Equation (17) can be evaluated to yield an analytical expression for the vector potential in the time domain, specifically where z * + = z * ± d * . In cases where the source current persists for a finite amount of time, the corresponding vector potential can be determined by integrating A over time t. For example, in the situation where a source electric current exists over the time period ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) ) Q V ) and (16) intensity I 0 , the vector potential, denoted by A STEP , is calculated as and the corresponding EM field is given by derivatives of A STEP . Note that in this calculation evaluated in the next section, we must ensure that the differences between regularized and non-regularized variables (i.e., t and t * ) are correctly accounted for.

Discussion and numerical examples
This section considers three points based on the analytical expressions derived for early EM signals. First, the onset times of early EM signals are considered, in order to impose principal limitations on the detectability of early EM signals. Second, the validity of using vacuum and conductive whole-space models is considered, in order to assess the viability of further simplifying the calculations. Third, the detectability of early EM signals is discussed for a few specific cases.
The numerical examples below use the following parameters. The current intensity I 0 and hypocentral distance d are set to 50 A and 10 km, respectively. For the proposed current generation mechanism (i.e., induction due to coseismic fault motion; Figure 1), I 0 = |SvV × B amb |, where S, V, and B amb represent the cross-sectional area of the conductor, the ground motion velocity, and the ambient geomagnetic field, respectively. An example situation for which I 0 = 50 A is |V| = 0.1 ms −1 , S = 10 km 2 , v = 0.1 Sm −1 , and |B| = 50000 nT, with V perpendicular to B. Given the numerical considerations of Iyemori et al. [1996], these values might be typical for a Mw 7.1 earthquake. However, for simplicity the following calculations are per-formed by approximating that the electric currents appear on a line, as in Equation (6); the discussions below still hold under this approximation.
The value of I 0 does not become considerably larger, even for earthquakes with much higher moment magnitudes. As long as we assume that induction is due to ground motion (i.e., Figure 1), the generated current intensity depends on the ground velocity, not the seismic moment. In addition, early EM signals do not scale with seismic moment, but instead scale with seismic moment released at the first stage of the rupture. Yet initial seismic moment release is never particularly large, even at Mw ≥ 9 [e.g., Ammon et al. 2005, Koketsu et al. 2011. Therefore, the ampli-  (18). When t < 0, we can use the path colored green; integration along Cconverges to zero for this case. When t >0, we adopt the black integral contour. Integration over C +1 and C +3 converges to zero, and C +2 becomes Gauss's integral in this case. tudes of early EM signals are not expected to increase, even for great earthquakes.

Characteristics of the solution
An important feature of temporal variations in early EM signals is the "delay" of the signal. First, let us analyze the form of the analytical expression for EM potential given in Equation (23). The arguments of the exponentials and error function terms in Equation (23) are proportional to r *2 /t * = nvr 2 /t or its square root, where r represents a value with the same order of magnitude as x and z. Therefore, it is expected that Equation (23) will have its extremal value at around t * = r *2 , which is equivalent to t = nvr 2 ≅ 10 −6 × vr 2 . The time at which the vector potential takes this extremal value is regarded as the delay time of the signal arrival.
For the case of r = 10 km = 10 4 m, for example, the delay time assuming v =10 -1 , 10 -2 , and 10 -3 Sm −1 is of the order of 10, 1, and 0.1 s, respectively. A typical value for the arrival time of a seismic wave at a hypocentral distance of 10 km is about 2 seconds. These results imply that arrivals of early EM signals will precede those of seismic waves unless conductivity is rather high (i.e., 10 -2 Sm −1 ).

Effect of the ground surface
Let us consider the temporal variations in the early EM signals for the two source geometries shown in Figure 4 using the analytical solutions developed above.
I consider the case whereby the temporal variation of the source current is given by the Heaviside step function H(t), which is zero for t < 0 and unity for t > 0. The intensity of the current I 0 and hypocentral distance d are assumed as 50 A and 10 km, respectively. This intensity value is realistic for the proposed current-generation mechanism (i.e., induction due to coseismic fault motion generation; see Figure 1), and is approximately equivalent to that expected for a ground motion velocity, fault width, thickness of the conductor around the fault, conductivity of the conductor, and ambient magnetic field intensity of 0.1 ms −1 , 10 km, 1 km, 0.1 Sm −1 , and 50000 nT, respectively.
We first consider the case whereby the source current is located directly beneath the observation point ("case 1", shown in Figure 4a). In this situation, the xcomponent of the magnetic field B x can be written as and the z-component is zero. Plots of B x for v = 1.0 × 10 −1 , 10 -2 , 10 -3 , and 10 -4 Sm −1 are shown in Figure 5.
Note that, regardless of the conductivity, this function converges to 1/2rz as expected from the Biot-Savart law for a steady current. For small values of v (i.e., 10 -3 or 10 -4 Sm −1 ), B x rapidly converges to its final value B x (0,0,t=∞). In contrast, for large values of v, the convergence is rather slow. These results demonstrate that predictions made using the Biot-Savart law are good approximations of the time-varying magnetic field for small v, but are poor for large v.
It should be noted that magnetic field variations calculated for the whole-space model shown in Figure  2b are also considerably different from those calculated for the two-layer model (Figure 2a). When we consider the whole-space model, the temporal variation of B x is given by Expanding the exponential, we can confirm that its convergence speed is of the order of t -1 . In contrast, from Equation (25) we see that the convergence speed of B x for the two-layer model is of the order of t -1/2 . These differing convergence rates are apparent in Figure 6, which compares B x for the two models assuming v = 10 −2 Sm −1 .
For comparison, we also consider the situation where both the source current and observation site are located on the boundary of the layer ("case 2", shown in Figure 4b). In this case, the z-component of the magnetic field B z for the two-layer model can be expressed as Numerical examples of the temporal evolution of B z are shown in Figure 7 for various conductivity values.
In contrast to the results obtained for the buried line-source, the whole-space solution provides a reasonable approximation to that of the two-layer model for the situation shown in Figure 4b when the source and receiver are located at the same depth. In this geometry, the convergence speed of B z is proportional to T -1 for both the two-layer and whole-space models. This result implies that the temporal variation in B z obtained using the whole-space model can be considered a good approximation to that predicted by the twolayer model, although there are some small differences. The level of agreement between the two models is apparent in Figure 8, which compares the temporal evolution of B z computed for half-space, conductive whole-space, and resistive whole-space models assuming v = 10 −2 Sm −1 .

Detectability of early EM signals
To assess the detectability of early EM signals, we examine expected amplitudes and shapes for a few specific cases. For earthquake early warning, the situation demonstrated in case 1 (Figure 4a) is most important, because seismic P-waves cannot be detected at observation sites before reaching at the ground surface in this situation. For this reason, we focus primarily on the situation in Figure 4a, and specifically on the generated magnetic field B x ( Figure 5). Note, however, that many of the following characteristics are unchanged for the geometry in Figure 4b and the resultant B z field in Figure 7.
The numerical examples show that a high resistivity is required for the signal to be detected. In resistive cases (i.e., v =10 -4 or 10 -3 S/m), an early EM signal appears almost at the seismic origin time, and the temporal pattern is not considerably different from a step function. Recall that this result was obtained by assuming that temporal variations in the source current were also given by a step function. The results for the resistive case therefore imply that the temporal variations in B largely reflect temporal variations in the source current. It was assumed that the source current is generated by EARLY EM SIGNALS BY LINE CURRENTS Figure 6. A comparison of the temporal variation in the horizontal component of the magnetic field for the source-receiver geometry shown in Figure 4a for a conductive whole-space (solid black line) and two-layer model (dotted line). For reference, the solution for a resistive whole-space (i.e., a vacuum) is also shown.   Figure 4b for whole-space (solid black line) and two-layer (dotted line) models. For reference, the solution for a resistive whole-space (i.e., a vacuum) is also shown. motional induction, so that the intensity of the source current is proportional to the slip velocity. Therefore, it can be concluded that temporal variations in B x represent temporal variations in the slip velocity (i.e., the earthquake source-time function) when v ≤ 10 -3 S/m. In contrast, in conductive cases, temporal variations in B x are quite emergent, and appear considerably different from the source-time function. Smoothed variations in the magnetic field are rather difficult to distinguish from background geomagnetic variations; hence, detection of early EM signals is expected to be difficult in a conductive medium.
Concerning the amplitudes of early EM signals, the present technique yields estimates of ~0.5-1.0 nT. Whether or not signals of these amplitudes can be resolved is primarily dependent on the type of instrument used. Most devices that have observed a coseismic or pre-seismic EM signal are either flux-gate magnetometers [e.g., Iyemori et al. 1996] or induction coils [e.g., Ujihara et al. 2004, Abdul Azeez et al. 2009]. Their minimum resolutions are 0.1 nT or finer; however, considering fluctuations in (or stability of ) the data, their practical resolution limit is ~1 nT. Finer-scale instruments, such as SQUID magnetometers [e.g., Katori et al. 2013], have better resolution; however, these sensors are costly, and dense installations are thus infeasible at present. Taking background noise into account, we can estimate the general detection limit of early EM signals at ~0.5-1.0 nT.
We can therefore tentatively conclude that signals are detectable in principle. Numerical examples demonstrate that 0.5-1.0 nT is the expected range of early EM signal amplitudes, assuming that the current intensity (I 0 = 50 A) is realistic. As mentioned at the beginning of this section, the current intensity (I 0 = 50 A) corresponds to a Mw 7 earthquake. Therefore, a Mw 7 or larger earthquake could be detectable before the P arrival if EM signals could be observed at distances of ~10 km.
It must be noted, however, that the above discussion only provides one possible scenario in which early EM arrivals can be detected. This is by no means comprehensive. A wide variety of V, S, and v values are possible in the scenario depicted in Figure 1. Exact values depend on the electrical structure of the fault and individual earthquake source processes. In this sense, the detection limit is not rigid. The ability to distinguish signal from noise is another limiting factor, though addressing this challenge falls outside the scope of this work.

Conclusions
In this paper, I have derived analytical expressions for the time-varying EM field generated by an impulsive line-current in both a two-layer and whole-space model. Some quantitative features of the predicted early EM signals are demonstrated through analysis of the expressions obtained. Concerning the arrival times of early EM signals, an important result is that the arrivals of early EM signals precede the arrivals of the first seismic waves by a few seconds unless conductivity is high (i.e., 10 −2 Sm −1 ). Consequently, with a low conductivity, early EM signals could be detectable in principle. Concerning the viability of the simplified model used in the calculations, the temporal evolution of early EM signals is well approximated by a wholespace model when the source is near the surface; for a source at depth, this approximation breaks down and it is necessary to consider a two-layer model instead. Concerning the detectability of early EM signals, we expect amplitudes of >0.5 nT if ground conductivity is not very high and the source current is sufficiently intense. Given that this exceeds the noise floor of most common types of magnetometers, early EM signals could be detected in principle. However, accurate realtime detection, and discrimination between signal and noise, are outstanding challenges that will be explored in future work.