Observatory geoelectric fields induced in a two-layer lithosphere during magnetic storms

Abstract

We report on the development and validation of an algorithm for estimating geoelectric fields induced in the lithosphere beneath an observatory during a magnetic storm. To accommodate induction in three-dimensional lithospheric electrical conductivity, we analyze a simple nine-parameter model: two horizontal layers, each with uniform electrical conductivity properties given by independent distortion tensors. With Laplace transformation of the induction equations into the complex frequency domain, we obtain a transfer function describing induction of observatory geoelectric fields having frequency-dependent polarization. Upon inverse transformation back to the time domain, the convolution of the corresponding impulse-response function with a geomagnetic time series yields an estimated geoelectric time series. We obtain an optimized set of conductivity parameters using 1-s resolution geomagnetic and geoelectric field data collected at the Kakioka, Japan, observatory for five different intense magnetic storms, including the October 2003 Halloween storm; our estimated geoelectric field accounts for 93% of that measured during the Halloween storm. This work demonstrates the need for detailed modeling of the Earth’s lithospheric conductivity structure and the utility of co-located geomagnetic and geoelectric monitoring.

Introduction

Magnetic storms are potentially hazardous for the activities and technological infrastructure of modern civilization (e.g., Baker et al. 2008; Cannon et al. 2013; Daglis 2004). Of particular concern is the storm-time induction of geoelectric fields in the Earth’s lithosphere. Corresponding voltage differentials can find their way into electric-power grids through ground connections, driving uncontrolled electric currents that interfere with grid operation, damaging transformers, tripping protective relays, and sometimes causing blackouts (e.g., Boteler et al. 1998; Kappenman 2012; Thomson et al. 2010). The reality of this complicated chain of causes and effects was demonstrated during the great magnetic storm of March 1989 (Allen et al. 1989), which caused the complete collapse of the Hydro-Québec electric-power grid in Canada, temporarily leaving nine million people without electricity (Bolduc 2002). The possible occurrence in the future of an even more intense magnetic storm, one that could cause widespread loss of electric power, has motivated regulatory agencies to issue operational standards intended to mitigate the effects of induction hazards. Concerns across the private sector have motivated reinsurance companies to commission related assessments of risk (e.g., Lloyds 2013; Riswadkar and Dobbins 2010).

In this broad context of practical concerns, the evaluation of induction hazards (e.g., Love et al. 2014; Pulkkinen et al. 2012; Viljanen 1997; Wei et al. 2013) relies on geomagnetic field data collected at ground-based observatories (e.g., Love 2008; Macmillan 2007), facilities operated around the world (e.g., Love and Chulliat 2013; Love and Finn 2011; Minamoto 2013), as parts of national geophysical projects for geomagnetic field mapping and space weather monitoring. Many magnetic observatories have been in continuous operation for decades, and a few have been in operation for over a century. The geomagnetic history recorded in observatory data is useful for analysis of evolving solar-terrestrial interaction and extreme space weather events. Magnetic observatory data are also used as input for the development and testing of algorithms that might be implemented in an operational setting for real-time evaluation of magnetic storm induction hazards.

In contrast, geoelectric field data are usually collected on a campaign basis at multiple stations deployed temporarily as parts of magnetotelluric (e.g., Chave and Jones 2012; Simpson and Bahr 2005) projects undertaken to infer localized lithospheric electrical conductivity (e.g., Schultz 2009; Uyeshima 2007);. More unusual are projects dedicated to the long-term collection of geoelectric data at geographically fixed stations. Such a project has, however, been maintained at the observatories of the Japan Meteorological Agency (JMA) for many years (Minamoto 2013), thus providing an opportunity to study the details of magnetic storm induction. We choose to analyze 1-s resolution magnetic and electric field data that were collected at JMA’s Kakioka observatory and which record intense magnetic storms.

A magnetic storm is seen in observatory data as a temporary enhancement in geomagnetic and geoelectric disturbance, often commencing suddenly, followed by a day or two of growth and evolution, and, finally, a gradual return to a state of relative quiescence. Such nonstationary variation can be described, in terms of mathematics, as a superposition of transient exponential moment functions defined with initial conditions. This is the basis of Laplace transformation from the time domain to the complex frequency domain. Laplace transformation conveniently reduces differential equations in time to algebraic equations in frequency, and time-domain convolution is reduced to simple multiplication. With a parameterized model of the electrically conducting lithosphere, one can obtain a transfer function of complex frequency that describes the physical relationship between transient geomagnetic variation and the resulting induced geoelectric field. Inverse Laplace transformation gives a time-domain impulse-response function that, when convolved with a time series recording geomagnetic variation, gives an estimate of the geoelectric field that is time-causal, a quality of importance for real-time operational algorithms.

Love and Swidinsky (2014) used a simple four-parameter model of lithospheric conductivity: a uniform half-space with a galvanic distortion tensor used to accommodate localized geoelectric polarization arising from regional three-dimensional conductivity. They showed that the corresponding impulse-response function could account for 87% of the geoelectric variation induced at Kakioka during the Halloween magnetic storm of October 2003, but the unmodeled residual 13% includes high-frequency (sub-minute) geoelectric variation that might be of importance for hazard assessment and mitigation. Following on from the work of Love and Swidinsky, and guided by geophysical and geological concepts, we have a practical goal: improved algorithmic estimation of storm-time geoelectric variation. We analyze the predictive utility of a two-layer model of lithospheric conductivity; with each layer parameterized by a separate distortion tensor, the corresponding transfer function can accommodate frequency dependence in geoelectric polarization. We use the nine-parameter time-domain impulse-response function to estimate geoelectric variation recorded in data from Kakioka during five different magnetic storms.

Observatory data

We analyze 1-s magnetic data time series from the Kakioka observatory (geographic coordinates: 36.23 N, 140.19 E). Reported by the JMA since 1983, these data were collected using a fluxgate sensor. Each datum is an instantaneous ‘spot’ measurement of the natural geomagnetic field variation. The magnetic field data have been calibrated for drift in sensor-system response. As a time series, the 1-s data are sufficient to resolve most of the time dependence of rapid and impulsive magnetic variation realized during magnetic storms. Measured vector values can be represented as a discrete sequence, B ^m(t _i ), for time-stamp values t _i ,t _i+1,t _i+2,···, where τ=t _i+1−t _i is the constant 1-s sampling interval. These are reported in terms of geographic-polar components (B _h horizontal intensity; ${B_{d}^{m}}$ declination; ${B_{z}^{m}}$ down). We convert the horizontal components to geographic-Cartesian components, $({B_{x}^{m}} \text {north}; {B_{y}^{m}} \text {east})$.

We also analyze 1-s electric data time series from the Kakioka observatory; 1-s resolution electric data have been measured at Kakioka and reported since 2000. These were collected by measuring the potential difference between pairs of grounded electrodes, oriented along north-south and east-west lines and separated by 190 m. Measured electric vector values are reported in geographic-Cartesian vector components $\textbf {E}^{m} = ({E_{x}^{m}} \text {north}; {E_{y}^{m}} \text {east})$; no measurements were made of the geoelectric vertical component. The electric field data contain a slow drift in baseline that is possibly an artifact of changes in electrode grounding, a common problem (e.g., Ferguson 2012); this part of the time series is not relevant to our analysis, and so, we subtract a trend line from each vector component. Both the magnetic and electric 1-s data are very complete in time for the durations we analyze; we fill a few gaps and replace a few locked values (consecutive values that are evidently not changing in time) by linear interpolation.

Five intense magnetic storms

We analyze data durations recording five intense magnetic storms, each with its own distinctive time-series variation; see Table 1. In Figure 1a, we show the measured $\textbf {B}^{m} = ({B_{x}^{m}} \text {north, blue}; {B_{y}^{m}} \text {east, green})$ geomagnetic time series for the Halloween storm of 29 to 31 October 2003 (Balch et al. 2004, Gopalswamy et al. 2005). This storm caused operational stress in the Swedish (Pulkkinen et al. 2005) and Scottish power grids (Thomson et al. 2005). It commenced with the arrival at the Earth of a coronal mass ejection, which abruptly compressed the magnetopause and generated a positive magnetic impulse in Kakioka ${B_{x}^{m}}$ of about 60 nT at 06:13 UT on 29 October. Complex geomagnetic variation was driven by the arrival of additional coronal mass ejections, solar wind variability, and vigorous magnetospheric convection. This resulted in main-phase intensification of the ring current and a corresponding decrease in Kakioka ${B_{x}^{m}}$, with one minimum coming early in the day of 30 October and another a day later. The ring-current index reached a maximum value of −D s t=383 nT at 23:30 UT on 30 October. The maximum time-rate of change in the north component of the geomagnetic field |Δ B _x | during 29 to 31 October, calculated from first differences of the discrete 1-s data, was 5.58 nT/s; the corresponding average |Δ B _x | was 0.07 nT/s.

Figure 1

Time series recording the Halloween storm of 29 October - 31 October 2003, each with 1-s resolution.(a) Measured Kakioka geomagnetic variation relative to a quiet baseline, measured north ${B_{x}^{m}}$ (blue) and measured east ${B_{y}^{m}}$ (green); (b) east Kakioka geoelectric field variation, measured ${E_{y}^{m}}$ (black) and estimated G E _y (purple); (c) residual difference d _y ; (d) the parts of estimated G E _y corresponding to induction in the top layer of the lithospheric model (red) and the lower half-space (gray); (e) north Kakioka geoelectric field variation, measured ${E_{x}^{m}}$ (black) and estimated G E _x (purple); (f) residual difference d _x ; and (g) the parts of estimated G E _x corresponding to induction in the top layer of the lithospheric model (red) and the lower half-space (gray).

Table 1 Summary of magnetic storm and model properties

In Figure 1b,e we show, respectively, the ${E_{y}^{m}} \text {east}$ and ${E_{x}^{m}} \text {north}$ components of the geoelectric field measured at Kakioka for 29 to 31 October 2003 (black). We note that, in general, when the geomagnetic field is active, so is the geoelectric field, but the two data types are not simply correlated in time. We also note that ${E_{y}^{m}}$ has a much larger variational amplitude than ${E_{x}^{m}}$. This polarization is the results of geomagnetic induction in the three-dimensional conductivity distribution of the lithosphere beneath the observatory site. In the case of Kakioka, this apparently allows geoelectric fields to be more easily induced in the east direction than in the north direction.

Among the other data durations analyzed here, one records the Bastille-Day storm of 15 to 16 July 2000 (Yermolaev et al. 2008), which commenced with a sudden impulse of large amplitude and subsequently exhibited substantial geomagnetic variability; both maximum |Δ B _x | and average |Δ B _x | were comparable to the October 2003 Halloween storm. The storm of 20 to 21 November 2003 had a main-phase maximum −D s t that was greater than that of the Halloween storm but maximum and average |Δ B _x | that were less than those of the Halloween storm. The storm of 05 to 07 November 2001 included two abrupt geomagnetic impulses, a maximum |Δ B _x | greater than that of the Halloween storm, but an average |Δ B _x | less than that of the Halloween storm. In summary, the variety of geomagnetic time-dependencies exhibited by the five storms summarized in Table 1 allows us to test our algorithm for estimation of storm-time geoelectric fields.

Transient induction for a flat Earth

In seeking a practical method for estimating magnetic storm induction of electric fields in the lithosphere beneath a single ground-based observatory, it is natural to consider an idealized model having Cartesian geometry, where a unit positional vector is given by $\boldsymbol {\hat {x}} = (\hat {x}, \hat {y}, \hat {z})$. The Earth’s surface is represented by the plane z=0; the electrically conducting lithosphere, possibly with regional-scale three-dimensional structure, is assumed to occupy the region (z>0); everything above the Earth’s surface (z<0) is assumed to be an insulator.

We consider geomagnetic field variation occurring over timescales of seconds to days and induction within the solid Earth over lithospheric depthscales of hundreds to thousands of kilometers. We can reasonably invoke the quasistatic approximation and neglect displacement currents in Maxwell’s equations. The spatial form and temporal evolution of the geomagnetic field B is then given by the diffusion equation:

$$ \nabla^{2}\textbf{B} - \mu\sigma\partial_{t}\textbf{B} = 0, $$

((1))

where t denotes time, σ is the electrical conductivity (assumed, for now, to be spatially homogeneous), and μ is magnetic permeability. Magnetic time variation is related to an electric field E by Faraday’s law:

$$ \nabla \times \textbf{E} + \partial_{t}\textbf{B} = 0. $$

((2))

Following assumptions commonly made in magnetotelluric analyses, we invoke, as is often standard, the plane-wave approximation (Chave and Weidelt 2012), whereby geomagnetic and geoelectric field variations are assumed to lie in a horizontal plane $(\hat {x}, \hat {y})$ parallel to the surface of the Earth. With this, then, and as a consequence of Equation (1), we can assume that the geomagnetic field has the following separation of variables:

$$ \textbf{B}(t,z) = \textbf{B}_{s}(z) + \textbf{B}_{-}(t)\mathrm{e}^{-kz} + \textbf{B}_{+}(t)\mathrm{e}^{+kz}, $$

((3))

where 1/k is the characteristic spatial depthscale and where B _s (z) is the secular part of the geomagnetic field that changes slowly in time compared to the duration of time considered here; we approximate it as time-steady, and we remark that its existence does not detectably affect our analysis. For convenience of plotting, we subtract a constant B(z) from the magnetic field data. Similarly, the geoelectric field has the form:

$$ \textbf{E}(t,z) = \textbf{E}_{-}(t)\mathrm{e}^{-kz} + \textbf{E}_{+}(t)\mathrm{e}^{+kz}, $$

((4))

and we assume that this does not have a significant secular part.

Transient time signals are conveniently analyzed using the Laplace transform. A hypothetical signal f(t) that varies over time and after some initial condition (t≥0) is conceived as the superposition of exponential moments (e.g., Butkov 1968, Chs. 5 and 8]. The transformation of f(t) into the complex frequency domain is given by:

$$ f(s) = {\mathcal L}\left\{f(t)\right\} = \int_{0}^{\infty}f(t)\mathrm{e}^{-st}dt, $$

((5))

where s is a complex frequency. Inverse transformation back to the time domain is given by

$$ f(t) = {\mathcal L}^{-1}\left\{f(s)\right\} = \frac{1}{2\pi\mathrm{i}}\int_{\gamma - \mathrm{i}\infty}^{\gamma + \mathrm{i}\infty}f(s)\mathrm{e}^{st}ds, $$

((6))

where $\mathrm {i} = \sqrt {-1}$ and where integration is along a contour path in the complex s-plane with the positive real number γ chosen so as to ensure convergence. Calculating an inverse Laplace transform using the Cauchy residue theorem often entails the discovery of a ‘keyhole’ contour path that circumvents singularities, but many inverse transformations, including that needed for this analysis, can be found in reference books. The Laplace transform of the time derivative of the function f(t) is:

$$ {\mathcal L}\left\{\partial_{t}f(t)\right\} = sf(s), $$

((7))

where we have assumed that the pre-response initial condition is assumed to be zero, f(0)=0.

The magnetic diffusion Equation (1) can be reduced to an algebraic equation by using the functional form (3) and the Laplace transform of the time derivative (7), giving:

$$ k^{2}\textbf{B} - s\mu\sigma \textbf{B} = 0. $$

((8))

From this, we obtain the characteristic equation describing the relationship between the spatial wave number k and the complex frequency s:

$$ k = \sqrt{s \mu\sigma}. $$

((9))

Next, recognizing that Faraday’s Law (2) describes an orthogonal relationship between the inducing magnetic field B and the induced electric field E, we use the characteristic Equation (9) to write:

$$ \textbf{R}\textbf{E}(s,z) = C(s, z)\cdot s \textbf{B}(s,z), $$

((10))

where:

$$ \textbf{R} =\left[ \begin{array}{rr} 0 & -1 \\ 1 & 0 \end{array} \right] \hspace{2mm}\text{and}\hspace{2mm} \textbf{R}^{-1} = \left[ \begin{array}{rr} 0 & 1 \\ -1 & 0 \end{array} \right] $$

((11))

are spin matrices that come from the curl operator and are only two-dimensional because we have invoked the plane-wave approximation. C(s,z) is the Schmucker-Weidelt (Schmucker 1973; Weidelt 1972) electromagnetic transfer function with impedance properties determined by the electrical conductivity of the lithosphere. C(s,z) has units of length.

Two-layer model

In developing a model of lithospheric conductivity, we start simple, and then add details and model parameters that might be considered to be, at least, physically plausible. If the model leads to improved predictions of geoelectric fields at Kakioka, then the model can be considered to be ‘useful.’ We begin by recalling Cagniard’s idealized model, which is often presented in textbooks on magnetotellurics: a one-dimensional lithosphere, consisting of a vertical stratification of N horizontal layers, each having uniform conductivity (Cagniard 1953). We denote the conductivities from the top layer to the bottom as σ _n =σ ₁,σ ₂···σ _N and the layer thicknesses as h _n =h ₁,h ₂···h _N , where the bottom Nth layer is assumed to be an infinite half-space, h _N =∞. We denote the horizontal-layer interfaces by the sequence of depths z _n =z ₀,z ₁···z _N−1, where z ₀=0 represents the surface of the Earth and z _N−1 is the top of the Nth layer. For this lithospheric geometry, Wait (1954) showed that, with recursive enforcement of continuity of magnetic and electric fields across each interface, the Schmucker-Weidelt C-functions can be expressed in terms of a recursion formula:

$$ C_{n}(s,z_{n-1}) = \frac{1}{k_{n}}\cdot\frac{k_{n}C_{n+1}(z_{n})+\text{tanh}(k_{n}h_{n})}{1+k_{n}C_{n+1}(z_{n})\text{tanh}(k_{n}h_{n})}, $$

((12))

where:

$$ C_{N}(s,z_{N-1}) = \frac{1}{k_{N}}. $$

((13))

A model with just two horizontal layers is a reasonable incremental extension of the half-space model considered by Love and Swidinsky (2014). For this, the C-function at the Earth’s surface reduces to:

$$ C_{1}(s,z_{0}) = \frac{1}{k_{1}}\cdot\frac{(k_{1}/k_{2})+\text{tanh}(k_{1}h_{1})}{1+(k_{1}/k_{2})\text{tanh}(k_{1}h_{1})}. $$

((14))

This can be simplified even further if the top layer only moderately attenuates magnetic variation - that the electromagnetic skin depth of the top layer is much greater than its thickness (e.g., Simpson and Bahr 2005, p. 25). Indeed, near the Earth’s surface, electrical conductivity can range over many orders of magnitude, but σ ₁≃10⁻³ S/m for the upper, say, h ₁≃30 km, would not be atypical. In this case, for a timescale of, say, 30 s, the electromagnetic skin depth 1/k ₁≃154 km. With this, it is certainly plausible that k ₁ h ₁≃0.19<<1, in which case tanh(k ₁ h ₁)≃k ₁ h ₁, and:

$$ C_{1}(s,z_{0}) \simeq \frac{1+k_{2}h_{1}}{k_{2}+{k_{1}^{2}}h_{1}}. $$

((15))

Still more simplification can be obtained if the electromagnetic skin depth of the lower half-space is much greater than the thickness of the top layer. In the lower lithosphere, olivine has a conductivity of σ ₂≃3×10⁻⁴ S/m (e.g., Yoshino and Katsura 2013). In this case, for a timescale of, say, 30 s, the electromagnetic skin depth 1/k ₂≃280 km. In this case, k ₂ h ₁≃0.11<<1, and:

$$ C_{1}(s,z_{0}) \simeq \frac{1}{k_{2}+{k_{1}^{2}}h_{1}}. $$

((16))

With this, we obtain the expression:

$$ C_{1}(s,z_{0}) \simeq \frac{1}{\mu \sigma_{1}h_{1}}\left(\frac{1}{s + \sqrt{as}}\right) \equiv C_{0}(s), $$

((17))

where we have defined, with simplified notation, the surface transfer function of interest C ₀(s). With these approximations, we identify a new timescale:

$$ \frac{1}{a} = \frac{\mu {\sigma_{1}^{2}} {h_{1}^{2}}}{\sigma_{2}}. $$

((18))

Furthermore, we note that C ₀(s) is functionally dependent on the conductivity of the lower half-space layer σ ₂ and the conductance of the top layer σ ₁ h ₁; there is no separate functional dependence on the conductivity of the top layer σ ₁ nor on its thickness h ₁.

We return to the time domain through inverse Laplace transformation (e.g., Abramowitz and Stegun 1965 eq. 29.3.43),

$$ \mathcal{L}^{-1}\left\{C_{0}(s) \right\} = C_{0}(t) = \frac{1}{\mu\sigma_{1}h_{1}}\text{exp}(at)\text{erfc}\left(\sqrt{at}\right), $$

((19))

which is the impulse-response function (e.g., Kanasewich 1981 Ch. 5) relating the time derivative of the inducing geomagnetic field to the induced geoelectric field. But calculating the time derivative of a discrete time series is numerically unstable to even relatively small amounts of random noise. As a remedy, we choose to consider the impulse-response function for straight-line interpolation between successive discrete data values. For this, we need the Laplace transform of the corresponding ‘ramp’ function. This can be obtained by integrating Equation (19):

$$ {C_{0}^{R}}(t) = {\int_{0}^{t}}C_{0}(t-\tau) d\tau. $$

((20))

With this, we obtain the two-parameter impulse-response function:

$$ {C_{0}^{R}}(t;a,b) = b\cdot\left[\text{exp}(at)\text{erfc}(\sqrt{at})+\frac{2}{\sqrt{\pi}}\sqrt{at}\right], $$

((21))

where:

$$ b=\frac{\sigma_{1}h_{1}}{\sigma_{2}} $$

((22))

is a new depthscale.

We choose, now, to divide the response function into two parts:

$$ {C_{0}^{R}}(t;a,b) = {T_{0}^{R}}(t;a,b)+{H_{0}^{R}}\left(t;ab^{2}\right). $$

((23))

The functional form of the first part:

$$ {T_{0}^{R}}\left(t;a,b\right) = b\cdot\text{exp}(at)\text{erfc}\left(\sqrt{at}\right), $$

((24))

is specified by two independent parameters, the frequency a and the depthscale b, each of which depend, as we have noted, on the conductance of the top layer σ ₁ h ₁ and conductivity of the lower half-space σ ₂. On the other hand, the functional form of the second part of (23) is specified by just one parameter a b ²:

$$ {H_{0}^{R}}\left(t;ab^{2}\right) = \frac{2}{\sqrt{\pi}}\sqrt{ab^{2}t}. $$

((25))

Since a b ²=(μ σ ₂)⁻¹, we note that this can be re-expressed as a function depending only on the conductivity σ ₂:

$$ {H_{0}^{R}}(t;\sigma_{2}) = \frac{2}{\sqrt{\pi\mu\sigma_{2}}}\sqrt{t}, $$

((26))

which we note is the ramp response function for induction in a uniformly conducting lower half-space layer (e.g., Pirjola 1984, eq. 7), previously used by Love and Swidinsky (2014) in their analysis of storm-time induction. From these observations, we recognize that induction within the top layer modifies induction within the lower half-space, but the total inductive effect is not simply an additive superposition: ${T_{0}^{R}}$ depends on both σ ₁ h ₁ and σ ₂, while ${H_{0}^{R}}$ only depends on σ ₂. This asymmetry is a remnant of the derivation of C ₀(s), in which the equations of electromagnetism are solved through recursive matching of boundary conditions, starting from the bottom layer and then working up to the top layer.

Application to discrete time series

To apply these results to linear interpolation between discrete data values, we need the ramp response function corresponding to the sampling interval τ between each data value. We define this as:

$$ C_{\tau}^{R}(t) \equiv {C_{0}^{R}}(t)U(t) - {C_{0}^{R}}(t-\tau)U(t-\tau), $$

((27))

where the Heaviside step function is:

$$ U(t-\tau) = \left\{ \begin{array}{c} 0 \\ 1 \end{array} \right. \text{for } \begin{array}{c} t < \tau \\ t > \tau \end{array} $$

((28))

Next, we recognize that multiplication in the complex frequency domain corresponds to convolution in the time domain, so that for two functions f and g:

$$ {\mathcal L}^{-1}\left\{f(s)\cdot g(s) \right\} = \left(\,f*g\right)(t). $$

((29))

Therefore, discrete values of the induced electric field, denoted E(t _i ), are obtained by discrete linear convolution with the ramp response function (27):

$$ \textbf{R}\textbf{E}(t_{i};a,b) =\left(C_{\tau}^{R}(a,b)* \frac{\Delta \textbf{B}}{\tau}\right)(t_{i}). $$

((30))

where:

$$ \Delta \textbf{B}(t_{j}) = \textbf{B}(t_{j+1}) - \textbf{B}(t_{j}). $$

((31))

In detail, convolution is given by the integral:

$$ \left(\,f*g\right)(t) = {\int_{0}^{t}}f(t-\theta)\cdot g(\theta) d\theta. $$

((32))

With this, Equation (30) can be written as:

$$ \textbf{RE}(t_{i};a,b) = \sum\limits_{j=1}^{i}C_{\tau}^{R}(t_{i}-t_{j};a,b)\cdot \frac{\Delta \textbf{B}(t_{j})}{\tau}, $$

((33))

which can be regarded as an applied form of Equation (10), suitable for algorithmic implementation.

Parameterization for 3D conductivity

Induction in the three-dimensional conductivity distribution of the lithosphere beneath an observatory site results in localized horizontal geoelectric fields having frequency-dependent polarization (e.g., Simpson and Bahr 2005, Ch. 5): Storm-time geoelectric fields will tend to be more easily induced in some azimuthal directions than in others, and the degree to which this is true will depend on the time rate of change of the geomagnetic field (e.g., McKay and Whaler 2006). Conversely, induction in a hypothetical half-space having uniform conductivity or, for that matter, in a Cartesian geometry with conductivity that is only depth dependent, with no lateral conductivity structure, will show no geoelectric polarization. Indeed, in a magnetotelluric analysis of data collected from a geographically distributed set of sensor systems, the frequency-dependent relationship between geomagnetic variation and polarized geoelectric variation is exploited in modeling three-dimensional lithospheric conductivity.

Residual polarization in the measured geoelectric field E ^m not directly accommodated by the lithospheric conductivity model is treated as a ‘galvanic’ distortion of the modeled geoelectric field E:

$$ \textbf{E}^{m}(t_{i})=\textbf{G}\textbf{E}(t_{i}), $$

((34))

where the distortion tensor:

$$ \textbf{G} = \left[ \begin{array}{rr} G_{xx} & G_{xy} \\ G_{yx} & G_{yy} \end{array}\right] $$

((35))

is assumed to be constant in time (frequency independent) (e.g., Groom and Bahr 1992; Jones 2012). Galvanic distortion is often regarded as a nuisance for magnetotelluric analysis. It is sometimes attributed to the quasi-static accumulation of electric charges along the boundaries of spatial heterogeneities in lithospheric conductivity and within conductivity gradients. In the context of induction-hazard science and the need to develop algorithms for estimating storm-time induction of geoelectric fields, lithospheric conductivity is often modeled as being only depth dependent (e.g., Pulkkinen et al. 2007). Recently, however, three-dimensional models have been used in the analysis of storm-time induction (e.g., Püthe and Kuvshinov 2013; Thomson et al. 2009; Viljanen et al. (2014); in their independent analyses of Kakioka geomagnetic and geoelectric data, Love and Swidinsky Love and Swidinsky (2014) and Püthe et al. 2014 each used a time-constant galvanic distortion tensor in their modeling of the polarizing effects of localized three-dimensional lithospheric conductivity.

Given our practical goal of obtaining improved estimates of storm-time geoelectric fields, we choose to extend the parameterizations of the preceding analyses in two different ways. First, we recognize that the development in the ‘Two-layer model’ section is an idealization, the details of which do not need to be slavishly retained. We choose to remove the parametric interdependence of the top ${T_{0}^{R}}$ and half-space ${H_{0}^{R}}$ response functions; instead of two independent parameters {a,b}, we allow for three {a _T ,b _T ,σ _H }, and instead of (23), we have:

$$ {C_{0}^{R}}(t;a_{T},b_{T},\sigma_{H}) = {T_{0}^{R}}(t;a_{T},b_{T})+{H_{0}^{R}}(t;\sigma_{H}). $$

((36))

With this, we can calculate independent top and half-space model electric fields:

$$ \textbf{E}(t_{i};a_{T},b_{T},\sigma_{H})=\textbf{E}_{T}(t_{i};a_{T},b_{T}) + \textbf{E}_{H}(t_{i};\sigma_{H}). $$

((37))

where, from Equation (30):

$$ \textbf{R}\textbf{E}_{T}(t_{i};a_{T},b_{T}) =T_{\tau}^{R}(a_{T},b_{T})*\frac{\Delta \textbf{B}}{\tau} $$

((38))

and:

$$ \textbf{R}\textbf{E}_{H}(t_{i};\sigma_{H}) =H_{\tau}^{R}(\sigma_{H})*\frac{\Delta \textbf{B}}{\tau}. $$

((39))

Second, we choose to assign independent distortion tensors to the model top and half-space electric fields:

$$ \textbf{E}^{m}(t_{i})=\textbf{G}_{T}\textbf{E}_{T}(t_{i}) + \textbf{G}_{H}\textbf{E}_{H}(t_{i}). $$

((40))

To avoid parametric redundancy, we require that each of the distortion matrices satisfy the Frobenius norm:

$$ \text{Tr}(\textbf{G}\textbf{G}^{T}) = 2. $$

((41))

With this, the independent model parameters, denoted {a _T ,b _T ,σ _H ,G}, are nine in number.

Since the response function for the top layer (24) has a time-functional form that is different from that of the half-space (26), by assigning separate distortion tensors to each of the top and half-space transfer functions, we are allowing geoelectric polarization to have frequency dependence, such as can arise from induction in a lithosphere having three-dimensional conductivity. This generalization permits improved fits of the model-estimated geoelectric fields to the measured (and, generally, polarized) geoelectric fields, which is our objective. If this also partly obscures the physical interpretation that might be made of the electrical conductivity properties of the top and half-space model lithospheric layers, this can be accepted so long as the same model parameters provide reasonably accurate estimates of induced geoelectric fields for different magnetic storms.

Parameter estimation

We define the residual discrepancy between the measured geoelectric field and that which we estimate from induction:

$$ \textbf{d}\left(a_{T},b_{T},\sigma_{H},\textbf{G}\right) = \textbf{E}^{m} - \textbf{G}_{T}\textbf{E}_{T} - \textbf{G}_{H}\textbf{E}_{H}, $$

((42))

Using the Kakioka magnetic and electric field data for each storm, summarized in the ‘Observatory data’ section, we obtain a set of model parameters by minimization of the squared misfit:

$$ \varepsilon^{2} = \left.\sum\limits_{i} \textbf{d}(t_{i}) \cdot \textbf{d}(t_{i})\right/\sum\limits_{i} \textbf{E}^{m}(t_{i}) \cdot \textbf{E}^{m}(t_{i}), $$

((43))

which we obtain with a Powell algorithm (Press et al. 1992, Ch. 10.5). We require a=|a|, and for numerical stability in calculating Equation (21), we first calculate $at+\text {log}\left [\text {erfc}(\sqrt {at})\right ]$ and then take its exponential. Given that the number of 1-s data values for a 3-day duration of time is large, 518,400, it is useful to design an efficient computer algorithm for minimizing ε ². We define data-fitting efficiency, or ‘variance reduction,’ as 1−ε ².

Results

In Table 1, we summarize lithospheric model parameters {a _T ,b _T ,σ _H ,G} and data misfits determined, separately, for each of the five storms. In Figure 1b,e, we compare, respectively east and north, components of the measured geoelectric field (E ^m, black) at Kakioka against the estimated geoelectric field (G E, purple) - purple variation approximately equals black variation. Optimal parameters (#4 in Table 1) for the two-layer model leave a misfit residual of only ε ²=0.0686; this means that the model accounts for 93% of the variance in the 2003 Halloween storm 1-s Kakioka E ^m data. These results are an improvement over those of Love and Swidinsky (2014), whose minimal half-space model accounted for 87% of the variance in the Kakioka electric field data for the Halloween storm but where high-frequency variation was not as well estimated. Still, even with the two-layer model presented here, some small residual signal remains, much of which varies over an approximate diurnal timescale. This can be seen in Figure 1c,f, where we plot the residual differences, Equation (42), between the measured and estimated geoelectric fields. The roles played, respectively, by the top and half-space lithospheric layers in affecting the estimated geoelectric field E can be seen in Figure 1d,g, where we plot the top G _T E _T (red) and the half-space G _H E _H (gray) time series. Note that the top layer mostly supports high-frequency electric field variation, while the half-space layer supports variation across a broad range of frequencies.

Geoelectric variation at Kakioka is even more accurately estimated for the other four storms. Still, the optimal model parameters {a _T ,b _T ,σ _H ,G} show what might seem to be substantial variance between estimates for different storms, especially for the model top layer. Since the conductance of the top layer allows modeling of high-frequency electric field variation, the relative amount of such variation among the five different storms ends up determining a relatively different set of parameters for the top layer. We seek to identify a more stable set of model parameters: the average values of {a _T ,b _T ,σ _H ,G} are listed in Table 1 and labeled ‘#A.’ We also list residual misfits ε ² for the average parameters for each of the five storms; in each case, the misfits are only slightly greater than the fits obtained by optimal parameters determined separately for each storm; for example, the #A parameters give a residual misfit of ε ²=0.0816 for the Halloween storm. In more detail, the degree to which the estimated geoelectric fields successfully match those measured at Kakioka can be seen in Figure 2, where we show just 20 min of 1-s E _y variation before and after the commencement of the Halloween storm. Given the complexity of the time series, and the relative simplicity of our model (just nine parameters), it is satisfying that we can estimate geoelectric variation with a high degree of accuracy.

Figure 2

Time series recording 20 min of east geoelectric field variation in the early stage of the 2003 Halloween storm. Shown are measured ${E_{y}^{m}}$ data (black), estimated G E _y variation for optimal model parameters determined from Kakioka data recording the 2003 Halloween storm (#4 in Table 1, purple), and estimated G E _y variation for model parameters determined from five different magnetic storms (#A in Table 1, green).

Interpretation of model lithosphere

Since our analysis here is concerned with the development and testing of an algorithm for estimating the storm-time induction of geoelectric fields, assumptions about the lithosphere and specific model parameterizations can be viewed in purely practical terms - they enable estimation. Detailed inference about the actual nature of the Earth’s lithosphere is not, itself, a goal. At the same time, we recognize that a realistic model of the lithosphere should, in principle, facilitate accurate geoelectric estimation, and so, it is reasonable to compare our simple model with what is either known or independently inferred about the lithosphere.

As noted in the previous section, rapid geomagnetic induction is supported by the conductance of the top layer of the model lithosphere. For the parameters of the average lithospheric model A, the diffusive timescale of the top layer is:

$$ \frac{1}{a_{T}} = 24.08\text{ s}, $$

((44))

which is relatively short, on the order of half the rise time of a magnetic impulse typically seen during an intense storm. The corresponding depthscale and normalized distortion tensor are:

$$ b_{T} = 47.50 \text{km}, \hspace{3mm} \textbf{G}_{T} = \left[ \begin{array}{rr} -0.03 & 0.02 \\ -0.70 & 1.23 \end{array} \right], $$

((45))

while the underlying half-space part of the model has the conductivity parameters:

$$ \sigma_{H} = 3.5\times10^{-4}\text{ S/m}, \hspace{3mm} \textbf{G}_{H} = \left[ \begin{array}{rr} 0.06 & 0.18 \\ -0.28 & 1.37 \end{array} \right]. $$

((46))

Of course, the Japanese archipelago and its regional surroundings are geologically and geophysically complicated. And, indeed, published magnetic sounding and magnetotelluric analyses of the lithosphere beneath and around Honshu (Fuji-ta et al. 1997; Kasaya et al. 2005; Ogawa et al. 1986), the island on which the Kakioka observatory is situated, reveal a complicated three-dimensional conductivity distribution: near the surface, sea water is a very good electrical conductor, 3 S/m; crustal conductivity ranges from about 10⁻² to 10⁻¹ S/m; below this, and across lithosphere depths down to about 30 km, the electrical conductivity is estimated to be relatively low, 10⁻⁴ S/m, but there are possibly isolated parts with dimensions of 10 km or so having conductivity of 10⁻¹ S/m; deeper still, and down to depths greater than 100 km, conductivity ranges from 10⁻² to 10⁻⁴ S/m, with lower conductivity structures possibly corresponding to the subducting Pacific and the overriding Okhotsk plates. The complicated regional distribution of electrical conductivity is almost certainly responsible for the high degree of polarization in the induced geoelectric field measured at Kakioka.

Püthe et al. 2014 have noted a seeming inconsistency between their estimate of the Kakioka galvanic distortion tensor, obtained for a spherical model, and that reported by Love and Swidinsky (2014, eq. 29), obtained for a Cartesian half-space. Püthe et al. speculate that this might be an artifact of the modeling method used by Love and Swidinsky. In fact, the explanation is much more mundane than that; Love and Swidinsky simply made a mistake in their reporting of the tensor elements: unintended rotation of the distortion tensor. Other results reported by Love and Swidinsky, including the average inferred conductivity and fits to the electric field data, are unaffected by this mistake. The elements reported here have been checked carefully for placement and sign; given expected statistical jitter, the half-space elements for G _H are certainly reasonably consistent with those estimated by Püthe et al. (2014).

Conclusions and new questions

The preceding work demonstrates the opportunity potential for algorithmic estimation of geoelectric fields induced in the Earth’s electrically conducting lithosphere during magnetic storms. It also illuminates the possibility that real-time algorithms for geoelectric estimation can be developed using a Laplace-transform formalism, thus ensuring time-causal dependence on (only) past geomagnetic field variation. Still, we can also see some challenges for induction-hazard science. The polarization of geoelectric data recorded at Kakioka requires parameterization for three-dimensional lithospheric conductivity, which we accomplish here, in terms of distortion tensors. From this, we can understand that mapping of hazardous geoelctric fields across the Earth’s surface will require detailed modeling of the three-dimensional lithospheric conductivity. Furthermore, by virtue of localized differences in lithospheric conductivity, storm-time geoelectric variation realized at one observatory site will not always be well correlated with geoelectric variation at another site. It would be of interest to study the correlation between pairs of geoelectric time series and how it falls off as a function of increasing site separation. This could be accomplished by deploying, on a temporary basis, measurement systems separated by different distances and then waiting for an intense magnetic storm to occur. Another issue that remains insufficiently studied is storm-time induction at high auroral-zone latitudes, where geomagnetic activity can have amplitudes much greater than at mid and low latitudes. Unfortunately, very few detailed geoelectric time series are available from auroral latitudes. In light of these limitations, we believe that a modest expansion of long-term geoelectric monitoring would be beneficial for induction-hazard science.