Numerical Calculation of Geoelectric Fields That Affect Critical Infrastructure

David H. Boteler, Risto J. Pirjola

Geomagnetic Laboratory, Natural Resources Canada, Ottawa, Canada.

DOI: 10.4236/ijg.2019.1010053                

Abstract

One of the modern applications of geomagnetism is determining the effect of geomagnetic disturbances on critical infrastructure such as power systems and pipelines. Assessing the geomagnetic hazard to such systems requires calculation of the geoelectric fields produced during geomagnetic disturbances. Such geoelectric fields can then be used as input to system models to calculate the impact on the system. This paper describes what is involved in calculating the geoelectric fields produced during real geomagnetic disturbances. The theory of geomagnetic induction is presented and used to derive the Earth transfer function relating the geoelectric and geomagnetic field variations at the Earth’s surface. It is then shown how this can be used to make practical calculations of the geoelectric fields and how the calculation process can be verified by comparison with analytic solutions obtained with synthetic geomagnetic variation data. The accuracy of the calculated geoelectric fields for geomagnetic risk assessments is limited, not by the accuracy of the calculation methods, but by the availability of geomagnetic field measurements and Earth conductivity information over the whole extent of the affected infrastructure.

Keywords

Geoelectric Fields, Geomagnetic Disturbances, Critical Infrastructure, Hazard Assessment

1. Introduction

Geomagnetism has always been an applied science because of the use of the magnetic field for navigation purposes, e.g. [1] [2]. In modern times a new application of geomagnetism has arisen. This follows the recognition that technological systems using long conductors can be affected by the geomagnetically induced currents (GIC) produced during geomagnetic disturbances [3]. The telegraph was the first system to suffer from natural voltages [4], and in 1859 the “Carrington storm” caused disruption to telegraph systems around the world [5]. Pipelines represent another type of long conductor subjected to geomagnetic induction. In this case, geomagnetic induction gives rise to “telluric currents” in the pipeline that may enhance corrosion and interfere with the corrosion prevention systems for the pipeline [6] [7]. Starting in 1940, geomagnetically induced currents produced in high-voltage power transmission systems during large geomagnetic disturbances have led to equipment damage and power blackouts [8] [9] [10] [11]. Concern about these effects has led to extensive work to understand the geomagnetic effects and to assess the hazard they pose to the operation of these systems, e.g. [12] [13] [14].

Assessing the geomagnetic hazard to grounded networks requires the use of geomagnetic field data and Earth conductivity models to calculate the geoelectric fields experienced by the system. The calculated geoelectric fields can then be used as input to system models to evaluate the geomagnetically induced currents (telluric currents) and the associated voltages and their impacts on the system. The method presented here uses standard theory widely used before, e.g. [8] [15] [16] [17] ; however, the details of the numerical calculations are usually ignored and can lead to errors if not done correctly.

The purpose of this paper is to present the numerical method for calculation of the geoelectric field. However, to show the origin of the relations used, we first present the theory of geomagnetic induction starting from Maxwell’s equations and show how this is used to derive the Earth transfer function relating the geoelectric field to the geomagnetic field variation at the Earth’s surface. Then we show how these calculations can be made in practice using a time series of sampled geomagnetic field data. Further, we show how the numerical calculation can be tested by use of synthetic magnetic field data and comparison of the calculated geoelectric fields with analytic solutions obtained for two test cases. Finally, the source of the Earth conductivity models is described and we discuss how limitations in the availability of such information influence the accuracy of geomagnetic hazard assessments.

2. Theory

The relations between electric and magnetic fields at the Earth’s surface are governed by Maxwell’s equations. For geomagnetic induction in the Earth we are dealing with frequencies ƒ and conductivities σ such that $\sigma \gg \text{2}\pi f{\epsilon }_0$ where ε₀ = 8.854 × 10⁻¹² F/m is the permittivity of free space. Consequently, the displacement current term associated with 2πƒε₀ is negligible and can be dropped from Maxwell’s equations. Thus the equations can be written as

$\nabla \times E=-\frac{\partial B}{\partial t}$ (1)

$\nabla \times B={\mu }_0\sigma E$ (2)

where bold denotes vectors and μ₀ = 4π × 10⁻⁷ H/m is the permeability of free space assumed to be valid in the Earth, too. Any geoelectric field E(t) and geomagnetic field variation B(t) can be expressed as the sum (i.e. an inverse Fourier transform) of their frequency components, ƒ:

$E\left(t\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}E\left(f\right){\text{e}}^{i2\pi ft}\text{d}f}$ (3)

$B\left(t\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}B\left(f\right){\text{e}}^{i2\pi ft}\text{d}f}$ (4)

where $i=\sqrt{-1}$.

The rate of change of the magnetic field denoted by g(t) can then be written as

$g\left(t\right)=\frac{\text{d}B\left(t\right)}{\text{d}t}={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}i2\pi fB\left(f\right){\text{e}}^{i2\pi ft}\text{d}f}$ (5)

Substituting (3), (4) and (5) into (1) and (2) then gives equations

$\nabla \times {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}E\left(f\right){\text{e}}^{i2\pi ft}\text{d}f}=-{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}i2\pi fB\left(f\right){\text{e}}^{i2\pi ft}\text{d}f}$ (6)

$\nabla \times {\displaystyle \underset{-\infty }{\overset{\infty }{\int }}B\left(f\right){\text{e}}^{i2\pi ft}\text{d}f}={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\mu }_0\sigma E\left(f\right){\text{e}}^{i2\pi ft}\text{d}f}$ (7)

Equations (6) and (7) imply that, in the frequency domain, E(ƒ) and B(ƒ) satisfy the following equations

$\nabla \times E\left(f\right)=-i2\pi fB\left(f\right)$ (8)

$\nabla \times B\left(f\right)={\mu }_0\sigma E\left(f\right)$ (9)

Using the standard geomagnetic xyz coordinate system in which the Earth’s surface is the xy plane and x, y and z refer to the northward, eastward and vertically downward directions, respectively, Equation (8) can be written in component form as follows

$\left(\frac{\partial E_z}{\partial y}-\frac{\partial E_y}{\partial z}\right)i+\left(\frac{\partial E_x}{\partial z}-\frac{\partial E_z}{\partial x}\right)j+\left(\frac{\partial E_y}{\partial x}-\frac{\partial E_x}{\partial y}\right)k=-i2\pi f\left(B_{x}i+B_{y}j+B_{z}k\right)$ (10)

where i, j and k are the unit vectors in the x, y and z directions, respectively. For clarity, the dependence of the electric and magnetic components on the frequency ƒ is not explicitly shown in Equation (10). A common assumption made is that the horizontal spatial variations of the fields are much less than the variation with depth (i.e. the skin depth in the Earth is much smaller than the horizontal spatial scales of the fields). This assumption is generally valid when “large-scale” ionospheric-magnetospheric sources are considered or the points of observation are far from the sources, and there are no significant lateral variations in the Earth’s conductivity. In this case, the terms with the variation in the vertical direction (z) dominate on the left-hand side of Equation (10), and so Equation (10) reduces to

$-\frac{\partial E_y}{\partial z}i+\frac{\partial E_x}{\partial z}j=-i2\pi f\left(B_{x}i+B_{y}j+B_{z}k\right)$ (11)

Equating the i, j and k components on both sides of Equation (11) gives

$\frac{\partial E_y}{\partial z}=i2\pi f{B}_x$ (12)

$\frac{\partial E_x}{\partial z}=-i2\pi f{B}_y$ (13)

and B_z = 0.

Thus a “transfer function” K = K(ƒ) can be defined that relates the orthogonal electric E = E(ƒ) and magnetic B = B(ƒ) field components in the frequency domain at the Earth’s surface, i.e.

$E\left(f\right)=K\left(f\right)B\left(f\right)$ (14)

where E(ƒ) and B(ƒ) can represent the relation between either E_x(ƒ) and B_y(ƒ) or −E_y(ƒ) and B_x(ƒ).

This transfer function K(ƒ) is closely related to the “surface impedance” Z(ƒ) commonly used in geoelectromagnetic studies: Z(ƒ) = μ₀K(ƒ). The electric field can also be related to the rate of change of the magnetic field, g(t) = dB(t)/dt. Noting that in the frequency domain g(ƒ) = i2πƒB(ƒ) (see Equation (5)) we obtain from (14)

$E\left(f\right)=K\left(f\right)B\left(f\right)=\frac{K\left(f\right)}{i2\pi f}i2\pi fB\left(f\right)=C\left(f\right)g\left(f\right)$ (15)

where C(ƒ) given by

$C\left(f\right)=\frac{K\left(f\right)}{i2\pi f}$ (16)

is called the “magnetotelluric relation”. It should also be noted that C(ƒ) is equal to the complex skin depth p(ƒ) often used in geoelectromagnetic investigations.

To determine the transfer function, K(ƒ), take the curl of Equation (8). Using Equation (9) and noting that $\nabla \cdot E\left(f\right)=0$, we then obtain at each frequency

(17)

where the propagation constant k = k(ƒ) is defined by

$k=k\left(f\right)=\sqrt{i2\pi f{\mu }_0\sigma }$ (18)

Expanding (17) gives

$\frac{{\partial }^{2}E}{\partial x^2}+\frac{{\partial }^{2}E}{\partial y^2}+\frac{{\partial }^{2}E}{\partial z^2}=k^{2}E$ (19)

Assuming again that the terms involving the variation with z dominate, Equation (19) reduces to

$\frac{{\partial }^{2}E}{\partial z^2}=k^{2}E$ (20)

Equation (20) can be used to obtain an expression for E as a function of z that can be used in Equations (12) and (13) to obtain an expression for the transfer function in the Earth relating E and B. Two types of solution are obtained: the first if the Earth is approximated by a uniform conductivity model, and the second if the Earth is approximated by a layered conductivity model. The second option provides a model that represents the change of conductivity with depth in the Earth. Taking account of lateral changes in conductivity is briefly considered in Section 6.

2.1. Uniform Conductivity Model

A simple first approximation is to model the Earth by a half-space of uniform conductivity σ. In this case, Equation (20) has a solution of the form

$E\left(f\right)=E_0\left(f\right){\text{e}}^{-k\left(f\right)z}$ (21)

where E₀(ƒ) is the (vector) amplitude of the electric field at the Earth’s surface, i.e. at z = 0, and k(ƒ) is given by Equation (18). Substituting the x or y components of E(ƒ) given by (21) into Equations (12) and (13) gives

$\frac{-E_y}{B_x}=\frac{i2\pi f}{k}$ (22)

$\frac{E_x}{B_y}=\frac{i2\pi f}{k}$ (23)

Equations (22) and (23) show that each pair of electric and magnetic field components form an orthogonal right-handed pair where the electric field is rotated 90 degrees counter clockwise from the direction of the magnetic field and the transfer function, K(ƒ) represents the relationship between E_x and B_y and between −E_y and B_x. These equations, utilizing (18), show that in the uniform-Earth case K(ƒ) is given by

$K=K\left(f\right)=\frac{i2\pi f}{k}=\sqrt{\frac{i2\pi f}{{\mu }_0\sigma }}$ (24)

2.2. Layered Conductivity Model

In practice the Earth conductivity varies in all directions, but the greatest variation is with depth and the Earth is often represented by a one-dimensional (1-D) model comprised of N horizontal layers with specified conductivities ( ${\sigma }_1,\cdots ,{\sigma }_N$ ) and thicknesses ( $l_{\text{1}},\cdots ,l_N$ ). Note that the bottom layer is a uniform half-space, so $l_N=\infty$.

As above, we assume that the horizontal variations of the electric and magnetic fields are much less than the variations with depth, i.e. the fields are considered to only depend on the z coordinate. Then the electric field satisfies Equation (20) in each layer with $k_n=\sqrt{i2\pi f{\mu }_0{\sigma }_n}$ ( $n=\text{1},\cdots ,N$ ). As in Section 2.1, we consider an (E, B) pair where E and B denote either E_x and B_y or −E_y and B_x. In the layered-Earth case the solution of Equation (20) needs to allow for upgoing waves resulting from reflections at layer boundaries (not present in the uniform-Earth case). Then the solution for E from (20) is

$E=S_n{\text{e}}^{-k_{n}z}+R_n{\text{e}}^{k_{n}z}$ (25)

where 0 ≤ z ≤ l_n and S_n and R_n are the amplitudes of a downward and an upward propagating wave at the top surface of layer n, respectively (note that, for each layer, the location of z = 0 is set at the top surface of the particular layer, not at the Earth’s surface). The bottom layer N is a uniform half-space, so R_N = 0.

Substituting E given by (25) into (12) or (13) yields

$B=\frac{k_n}{i2\pi f}\left(S_n{\text{e}}^{-k_{n}z}-R_n{\text{e}}^{k_{n}z}\right)$ (26)

As defined by Equation (14), the transfer function K = K(ƒ) relates E and B at the Earth’s surface. Let us denote the ratio of E to B at the top surface of layer n by K_n. The bottom layer (i.e. n = N) is a uniform half-space so the value of K_N is directly obtained from Formula (24) with σ = σ_N. Thus

$K_N=\frac{i2\pi f}{k_N}=\sqrt{\frac{i2\pi f}{{\mu }_0{\sigma }_N}}$ (27)

For each layer above, Equations (25) and (26) are combined and we set $z=0$ to give the ratio of the fields at the top of the layer, i.e. K_n, and set $z=l$ to give the ratio of the fields at the bottom of the layer which, because of continuity of fields across the boundary, is the same as the ratio of the fields at the top of the underlying layer, i.e. $K_{n+1}$,

$K_n={\eta }_n\frac{S_n+R_n}{S_n-R_n}$ (28)

$K_{n+1}={\eta }_n\frac{S_n{\text{e}}^{-k_n{l}_n}+R_n{\text{e}}^{k_n{l}_n}}{S_n{\text{e}}^{-k_n{l}_n}-R_n{\text{e}}^{k_n{l}_n}}$ (29)

where ${\eta }_n=\frac{i2\pi f}{k_n}$ is the “Characteristic Function” of layer n.

Dividing the top and bottom of Equation (29) by ${\text{e}}^{k_n{l}_n}$ gives

$K_{n+1}={\eta }_n\frac{S_n{\text{e}}^{-2k_n{l}_n}+R_n}{S_n{\text{e}}^{-2k_n{l}_n}-R_n}$ (30)

Multiplying by the denominator on the right-hand side gives

$K_{n+1}\left(S_n{\text{e}}^{-2k_n{l}_n}-R_n\right)={\eta }_n\left(S_n{\text{e}}^{-2k_n{l}_n}+R_n\right)$ (31)

Collecting terms in $S_n$ and $R_n$ gives

$\left(K_{n+1}-{\eta }_n\right)S_n{\text{e}}^{-2k_n{l}_n}=\left(K_{n+1}+{\eta }_n\right)R_n$ (32)

Thus we can write

$\frac{\left(K_{n+1}-{\eta }_n\right)}{\left(K_{n+1}+{\eta }_n\right)}{S}_n{\text{e}}^{-2k_n{l}_n}=R_n$ (33)

Using this expression for $R_n$ we can write

$S_n+R_n=\frac{\left(K_{n+1}+{\eta }_n\right)}{\left(K_{n+1}+{\eta }_n\right)}{S}_n+\frac{\left(K_{n+1}-{\eta }_n\right)}{\left(K_{n+1}+{\eta }_n\right)}{S}_n{\text{e}}^{-2k_n{l}_n}$ (34)

and

$S_n-R_n=\frac{\left(K_{n+1}+{\eta }_n\right)}{\left(K_{n+1}+{\eta }_n\right)}{S}_n-\frac{\left(K_{n+1}-{\eta }_n\right)}{\left(K_{n+1}+{\eta }_n\right)}{S}_n{\text{e}}^{-2k_n{l}_n}$ (35)

Substituting these into Equation (28) and cancelling common terms give

$K_n={\eta }_n\frac{\left(K_{n+1}+{\eta }_n\right)+\left(K_{n+1}-{\eta }_n\right){\text{e}}^{-2k_n{l}_n}}{\left(K_{n+1}+{\eta }_n\right)-\left(K_{n+1}-{\eta }_n\right){\text{e}}^{-2k_n{l}_n}}$ (36)

Collecting terms in $K_{n+1}$ and ${\eta }_n$ gives

$K_n={\eta }_n\frac{K_{n+1}\left(1+{\text{e}}^{-2k_n{l}_n}\right)+{\eta }_n\left(1-{\text{e}}^{-2k_n{l}_n}\right)}{K_{n+1}\left(1-{\text{e}}^{-2k_n{l}_n}\right)+{\eta }_n\left(1+{\text{e}}^{-2k_n{l}_n}\right)}$ (37)

Formula (37) is a recursive relation for calculating K_n at the top surface of layer n from $K_{n+1}$ at the top surface of the underlying layer n + 1. The initial value in the application of Formula (37) is K_N given by Equation (27) which is used to calculate K_n at the top surface of the next layer up. This is then used in the calculation for the next layer, and so on, up to the Earth’s surface. The final value is the transfer function K = K₁ relating E and B at the Earth’s surface.

It should be noted that the multi-layer transfer function is exact, as is the uniform-Earth transfer function, but in the multi-layer case the transfer function can only be expressed by the recursive relation (37) with the initial value from (27) whereas the transfer function for a uniform Earth has a simple explicit Formula (24).

The calculations are repeated for each frequency required to build up the full description of the dependence of the transfer function on frequency.

2.3. Expression for the Geoelectric Field

To determine the time variation of the electric field, we use Equation (3) and note that each frequency component of the electric field is given by the transfer function and the magnetic field at that frequency (Equation (14)). Thus we obtain

$E\left(t\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}K\left(f\right)B\left(f\right){\text{e}}^{i2\pi ft}\text{d}f}$ (38)

As shown in Sections 2.1 and 2.2, the Earth transfer function, K(ƒ), can be derived as a function of frequency for a uniform conductivity model or a layered conductivity model of the Earth. To use this to calculate the geoelectric field, note that the B(ƒ) terms in (38) are given by the Fourier transform of the geomagnetic field variation B(t). Thus (38) can be written as

$E\left(t\right)={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}K\left(f\right)\left[{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}B\left(t^{\prime }\right){\text{e}}^{-i2\pi f{t}^{\prime }}\text{d}{t}^{\prime }}\right]{\text{e}}^{i2\pi ft}\text{d}f}$ (39)

This can be more simply written

$E\left(t\right)=F^{-1}\left[K\left(f\right)F\left\{B\left(t\right)\right\}\right]$ (40)

where F and F⁻¹ represent the forward and inverse Fourier transform, respectively.

Thus the calculation of the geoelectric field comprises three steps:

1) Take the Fourier transform of the geomagnetic field variation B(t):

2) Multiply each frequency component by the transfer function of the Earth K(ƒ):

3) Take the inverse Fourier transform to obtain the geoelectric field E(t).

3. Sources of Data

Equation (39) shows that the calculation of geoelectric fields requires a time series of magnetic field measurements B(t) and the Earth transfer function K(ƒ) dependent on the Earth’s conductivity structure.

3.1. Magnetic Field Measurements

Measurements of the Earth’s magnetic field are made at many observatories around the world. Most observatories are members of Intermagnet, which promotes common instrument standards and data collection methods. Data from Intermagnet observatories are available through the Web site http://www.intermagnet.org. In addition, a number of magnetometer chains are operated by universities and research institutes, and these data have been collected together on the web site http://supermag.jhuapl.edu/mag/.

Magnetic field measurements are typically made using Fluxgate magnetometers. A set of three orthogonal magnetometers are used to measure the X, Y and Z magnetic field components in the northward, eastward and vertically downward directions, respectively (Figure 1).

Magnetic field recordings have traditionally been made with a sampling interval of one minute, with higher frequencies filtered out prior to sampling to remove frequencies above the Nyquist frequency to prevent aliasing problems [18]. Recordings have also been made with sampling once every 10 s, or once every 5 s, and the new Intermagnet standard for observatories is to make recordings with a sampling interval of 1 s.

3.2. Earth Transfer Function

The Magnetotelluric (MT) technique is a geophysical method that uses recordings of the geoelectric field and geomagnetic field variations to obtain information about the Earth conductivity structure [19] [20] [21]. Because the fields at

Figure 1. Geomagnetic coordinate system. Magnetic field components, X, Y and Z measured in the northward (x), eastward (y) and vertically downward (z) directions.

different frequencies penetrate to different depths within the Earth, spectral analysis of the geoelectric and geomagnetic field recordings can be used to determine the relationships between E and B at different frequencies and build up a (1-D) profile of the Earth conductivity variation with depth. Magnetotelluric surveys along transects or in arrays can be used to provide 2-D or 3-D models of the Earth conductivity structure.

The results from MT studies can be used in several ways. Published results from MT surveys usually describe the conductivity structure. These can be used to construct a 1-D layered Earth model, which is then used to calculate the Earth transfer function as shown in Sections 2.1 and 2.2. The MT results may show different conductivity values in different zones and a 1-D model can be constructed for each zone with the results being used to obtain “piecewise” calculations of the geoelectric fields affecting critical infrastructure [22]. Alternatively, where they are available, 2-D or 3-D models could be used to determine the transfer functions at the Earth’s surface.

Magnetotelluric studies, in the course of their data processing, generate surface impedance functions for each of their measuring sites. Some studies make these impedance functions available online (e.g. the Lithoprobe studies in Canada available at http://lithoprobe.eos.ubc.ca/ and the Earthscope studies in the US available at http://www.earthscope.org. Using the relation K(ƒ) = Z(ƒ)/μ₀, the surface impedance Z(ƒ) directly gives the Earth transfer function K(ƒ) to be applied to the calculation of the geoelectric field.

4. Numerical Calculations

The first main step in numerical calculation of geoelectric fields, as specified in Section 2.3, is to take the Fourier transform of the magnetic field data for the interval of interest. However, before taking the Fourier transform some preconditioning of the magnetic field data is needed. This follows standard practice for time series analysis, e.g. [18]. To reduce spectral leakage it is desirable to remove the mean and linear trend from the time series and taper the ends by multiplying the time series with a split cosine bell window (or another suitable window). The split cosine bell window is given by

$w_p\left(t\right)=\{\begin{array}{l}\frac{1}{2}\left(1-\cos \left(\frac{2\pi t}{p}\right)\right),0\le t<\frac{p}{2}\\ 1,\frac{p}{2}\le t<1-\frac{p}{2}\\ \frac{1}{2}\left(1-\cos \left(\frac{2\pi \left(1-t\right)}{p}\right)\right),1-\frac{p}{2}\le t<1\end{array}$ (41)

where p is the proportion (usually 10%) of the time series to be tapered [23].

The tapering prevents any spurious frequency components from being introduced by the discontinuity represented by the end of the time series, but has the disadvantage that the calculated geoelectric fields will also be tapered. To overcome this when calculating the geoelectric fields for a particular day, it is recommended that the days on either side be included in the time series for analysis. The geoelectric field is then calculated for the set of 3 days and the results for the first and last day, which are affected by the tapering, are discarded leaving the correct geoelectric field values for the middle day of the set.

The Fourier transform of the geomagnetic data time series can be performed using a Fast Fourier Transform (FFT) that is built into many software packages. There is no standard convention for the definition of a forward transform or an inverse transform and these can differ from one software package to another. Fortunately, in this application where the forward and inverse transforms are always used as a pair, the choice of definition does not matter, although this is not the case when an individual transform is used [24]. Many FFT routines expect the input data to be of a length that is a power of 2 and be comprised of complex numbers. This is easily achieved by padding the geomagnetic data set with zero values at the end to extend the data length to a power of 2 and by assigning the geomagnetic data values to be the real parts and setting the imaginary parts to zero. This complex array is often overwritten by the output array from the FFT routine.

The Discrete Fourier Transforms to go between the samples of the continuous time recording (x) and discrete frequency components (X) can be defined as [24] [25].

$X\left(k\right)=\frac{1}{N}{\displaystyle \underset{n=0}{\overset{N-1}{\sum }}x\left(n\right){\text{e}}^{-i\frac{2\pi kn}{N}}}$ (42)

$x\left(n\right)={\displaystyle \underset{n=0}{\overset{N-1}{\sum }}X\left(k\right){\text{e}}^{i\frac{2\pi kn}{N}}}$ (43)

When performing the forward FFT from the time domain to the frequency domain the output is the complex spectral components at positive and negative frequencies. However, the negative frequency components are output at the end of the array (Figure 2), which is a mathematical property of the Discrete Fourier Transform, not just a property of the FFT routine used. This needs to be taken into account when multiplying by the Earth transfer function values.

Figure 2. Spectral values at negative frequencies are output at end of array.

The output of the FFT is the spectral components of the magnetic field variation at a set of frequencies determined by the time-domain sampling interval, ∆t, and the length of the time series, N, in the time domain that is transformed. These frequencies

$f_i=0,\frac{+1}{N\Delta t},\frac{+2}{N\Delta t},\cdots ,\frac{+\left(\frac{N}{2}-1\right)}{N\Delta t},\frac{+\frac{N}{2}}{N\Delta t},\frac{-\left(\frac{N}{2}-1\right)}{N\Delta t},\cdots ,\frac{-2}{N\Delta t},\frac{-1}{N\Delta t}$ (44)

are those for which the Earth transfer function must be calculated. For each of the positive frequencies multiply the magnetic field spectral value by the Earth transfer function to give the corresponding electric field spectral value. The electric field spectral values for negative frequencies should then be set to the complex conjugates of the values for positive frequencies.

The resulting electric field spectrum (as with the magnetic field spectrum) has an even real part and an odd imaginary part. The odd imaginary part goes through zero at the Nyquist frequency (ƒ_N = 1/2∆t) (and at zero frequency). This requires setting the imaginary part of the E spectrum at ƒ_N equal to zero.

Now perform an inverse FFT on the electric field spectrum. This will give an array of complex values that contains the electric field in the time domain. Just as with the original magnetic field data, the electric field data must be real. Thus the computed array of electric field values should comprise complex values for which the imaginary parts are all zero. In practice, non-zero values are obtained for the imaginary parts but these should be many orders of magnitude (~10) smaller than the real part values. This should be checked as a test that the calculations have been made correctly. Small differences in the above procedure, e.g. not setting the imaginary part at ƒ_N to zero can make the imaginary parts of the computed electric field values much larger than they should be.

5. Verification of the Calculations

Analytic expressions have been derived for synthetic magnetic test data and electric fields at the surface of a uniform Earth and a layered Earth [26] [27].

A synthetic test magnetic field variation has been constructed as the sum of six sine waves:

$B\left(t\right)={\displaystyle {\sum }_{m=1}^6{A}_m\sin \left(2\pi f_{m}t+{\Phi }_m\right)}$ (45)

with amplitudes, A_m approximately proportional to 1/ƒ_m, and phases, Φ_m assigned arbitrary values as shown in Table 1.

Test calculations are performed using two Earth models: 1) a uniform Earth model with a resistivity of 1000 Ω·m; 2) a layered Earth model, used for Québec and consisting of 5 layers with, from the top down, thicknesses and resistivities: 15 km, 20,000 Ω·m; 10 km, 200 Ω·m; 125 km, 1000 Ω·m; 200 km, 100 Ω·m; above a half-space of 3 Ω·m [28].

The formulas for a uniform Earth model (24) or a layered Earth model (37) give the transfer function as a complex number at each frequency. The transfer

Table 1. Parameters of the synthetic test magnetic field variations.

function amplitude is the square root of the sum of the squares of the real and imaginary parts

$\left|K\left(f\right)\right|=\sqrt{{\left[RE\left\{K\left(f\right)\right\}\right]}^2+{\left[IM\left\{K\left(f\right)\right\}\right]}^2}$ (46)

And the phase is given by

$Phase\left\{K\left(f\right)\right\}={\tan }^{-1}\left\{\frac{IM\left\{K\left(f\right)\right\}}{RE\left\{K\left(f\right)\right\}}\right\}$ (47)

In the case of a uniform Earth model the transfer function is given by (24).

Writing the square root of “i” as a complex number $\sqrt{i}=\frac{1+i}{\sqrt{2}}$ we can write Equation (24) as

$K\left(f\right)=\left(1+i\right)\sqrt{\frac{\pi f}{{\mu }_0\sigma }}$ (48)

Thus the real and imaginary parts of K(ƒ) are equal and given by $\sqrt{\frac{\pi f}{{\mu }_0\sigma }}$.

Substituting these into (46) gives

$Amplitude\left\{K\left(f\right)\right\}=\sqrt{\frac{\pi f}{{\mu }_0\sigma }+\frac{\pi f}{{\mu }_0\sigma }}=\sqrt{\frac{2\pi f}{{\mu }_0\sigma }}$ (49)

Substituting into (47), because the real and imaginary parts are equal, gives

$Phase\left\{K\left(f\right)\right\}={\tan }^{-1}\left\{1\right\}={45}^{\circ }$ (50)

This shows that a uniform medium, regardless of the conductivity, σ, gives a transfer function with a phase of 45˚.

For a layered Earth model the recursive Formula (37) gives complex values for K(ƒ) which do not have a simple analytic expression. In this case, the appropriate way to calculate the amplitude and phase is to use the calculated real and imaginary parts of K(ƒ) in Equations (46) and (47).

The transfer function amplitudes |K_m| and phases θ_m for the uniform and layered Earth models, for the six frequencies in the synthetic magnetic data, are shown in Table 2 and Table 3.

Combining the six magnetic field frequency components with the corresponding transfer function values gives the electric field.

Table 2. Transfer function for the uniform Earth model.

Table 3. Transfer function for the layered Earth model.

$E\left(t\right)={\displaystyle {\sum }_{m=1}^6\left|K_m\right|A_m\sin \left(2\pi f_m+{\Phi }_m+{\theta }_m\right)}$ (51)

Multiplying the |K_m| and A_m terms and adding the phases in (51) gives the electric field

$E\left(t\right)={\displaystyle {\sum }_{m=1}^6{E}_m\sin \left(2\pi f_m+{\phi }_m\right)}$ (52)

with the six electric field waves as given in Table 4 for the uniform Earth model, and in Table 5 for the layered Earth model.

The synthetic geomagnetic field variation given by Equation (45) and Table 1 is shown in Figure 3. The exact geoelectric fields obtained from Equation (52) and Table 4 for the uniform Earth model and Table 5 for the layered Earth model are shown in Figure 4(a) and Figure 4(b), respectively. These waveforms provide analytic test cases that can be used to test the accuracy of numerical calculations. It is necessary to emphasize, referring to Section 2, that B expressed by Equation (45) with the parameter values given in Table 1 and shown in Figure 3 and E expressed by Equation (52) with the parameter values given in Table 4 or Table 5 and shown in Figure 4(a) or Figure 4(b) represent either E_x and B_y or −E_y and B_x.

The synthetic magnetic field data, as shown in Figure 3, comprise a set of 3 days of geomagnetic field variations with a sampling interval of 1-minute. This was used as input to the process for numerical computation described in Section 4. The resulting calculated geoelectric fields for Day 2 of the interval were then compared with Day 2 of the exact analytic solutions shown in Figure 4(a) and Figure 4(b). The results of that comparison are shown in Figures 5(a) (uniform Earth) and Figure 5(b) (layered Earth).

Figure 3. Test magnetic field variation.

(a)(b)

Figure 4. Exact electric fields. (a) For the uniform Earth model; (b) For the layered Earth model.

Table 4. Parameters of six electric field waves for the uniform Earth model.

Table 5. Parameters of six electric field waves for the layered Earth model.

(a)(b)

Figure 5. Comparison between numerical and analytic electric field values. (a) For the uniform Earth model; (b) For the layered Earth model.

Making a least-square fit to the data points in Figure 5(a) and Figure 5(b) as follows

$E_{numerical}=a{E}_{analytic}+b$ (53)

gives the values a = 0.999939118 (uniform), a = 1.000000135 (layered), b = 0.3250377 mV/km (uniform) and b = 0.042586 mV/km (layered). The good agreement between the numerical and analytic geoelectric fields shown by these numbers and correlation coefficients of 0.999999343 (uniform) and 0.999999989 (layered) confirm the accuracy of the numerical computation of geoelectric fields.

6. Calculation Using Tensor Transfer Functions

In the situations considered above where the Earth is approximated by a uniform half-space or a layered conductivity model the geoelectric fields produced are always orthogonal to the geomagnetic field variations producing them. However, in reality the Earth has a three-dimensional structure featuring large-scale conductivity changes such as at a coastline or smaller scale geological features with different conductivities. These all create situations where the electric field direction is deflected away from the orthogonal direction as shown in Figure 6.

Figure 6. Deflection of electric fields by conductivity structure in the Earth.

In these cases the geoelectric fields and geomagnetic field variations are related by a tensor transfer function in the frequency domain:

$\left[\begin{array}{c}{E}_x\\ E_y\end{array}\right]=\left[\begin{array}{cc}{K}_{xx}& K_{xy}\\ K_{yx}& K_{yy}\end{array}\right]\left[\begin{array}{c}{B}_x\\ B_y\end{array}\right]$ (54)

Thus the northward component of the geoelectric field, E_x, can be written

$E_x\left(t\right)=E_{xx}\left(t\right)+E_{xy}\left(t\right)$ (55)

where

$E_{xx}\left(t\right)=F^{-1}\left[K_{xx}\left(f\right)F\left\{B_x\left(t\right)\right\}\right]$ (56)

$E_{xy}\left(t\right)=F^{-1}\left[K_{xy}\left(f\right)F\left\{B_y\left(t\right)\right\}\right]$ (57)

Similarly, the eastward component of the geoelectric field, E_y, can be written

$E_y\left(t\right)=E_{yx}\left(t\right)+E_{yy}\left(t\right)$ (58)

where

$E_{yx}\left(t\right)=F^{-1}\left[K_{yx}\left(f\right)F\left\{B_x\left(t\right)\right\}\right]$ (59)

$E_{yy}\left(t\right)=F^{-1}\left[K_{yy}\left(f\right)F\left\{B_y\left(t\right)\right\}\right]$ (60)

It can be seen that these equations are just repeating the process used in the 1-D calculation. Therefore the numerical calculation method described in Section 4 can be used with the components of the tensor transfer function values and the verification process described in Section 5 can be used to test the calculation of the geoelectric field parts given by (56), (57), (59) and (60).

In Equation (54) the minus sign in the relationship between E_y(ƒ) and B_x(ƒ) mentioned in connection with the 1-D transfer function K(ƒ) used in (14) is absorbed into the K_yx(ƒ) term. The K_xx(ƒ) and K_yy(ƒ) terms can be positive or negative depending on the conductivity structure. In an area where the structure is 1-D the tensor transfer function terms K_xx(ƒ) and K_yy(ƒ) are zero, and the diagonal terms are given by K_xy(ƒ) = K(ƒ) and K_yx(ƒ) = −K(ƒ).

7. Discussion

The above verification tests relate to the process for numerical computation of geoelectric fields. The accuracy of electric fields calculated for any particular site also depends on how well the Earth models for that location represent the Earth conductivity structure that is actually there. As mentioned above, the Earth has a complicated 3-dimensional conductivity structure and any Earth model used is inevitably an approximation to the real situation. One-dimensional models only take account of the variation of conductivity with depth, while 2-dimensional and 3-dimensional models attempt to represent some of the lateral variations in conductivity.

The Earth transfer functions used in the above calculations can either be obtained from models or directly from magnetotelluric (MT) measurements. The Earth conductivity models are also derived from MT measurements, so the two approaches should be equivalent. However, 2-D and 3-D models are derived by combining data from an array of MT sites so the Earth transfer function based on the model may not exactly reproduce the transfer functions measured at individual MT sites.

MT measurements are not available everywhere. To provide Earth conductivity information for areas with no MT measurements, zones with similar geologic structure have been identified and the MT results from a survey crossing part of a zone are applied to the rest of the zone. The geoelectric field can then be calculated in each zone. These “piecewise” calculations then provide geoelectric field values across the whole region covered by a power grid and provide a simple approximate way of including Earth conductivity changes in the calculation of geomagnetically induced currents in a power network [13] [22] [29] [30] [31].

8. Conclusions

Calculation of the geoelectric fields produced during geomagnetic disturbances at the Earth’s surface is a key requirement for assessing the geomagnetic hazard to critical infrastructures such as power systems and pipelines.

The geoelectric field can be calculated by taking a Fourier Transform of a time series of magnetic field data, multiplying the magnetic field spectral components by the Earth transfer function to obtain the electric field spectrum, and then taking the inverse Fourier transform to obtain the geoelectric field in the time domain.

Accurate calculations require preprocessing (remove mean and trend and taper with a split cosine bell) of the magnetic field data as well as taking care of the placement of positive and negative frequency terms in the transform array and combining these with the appropriate Earth transfer function terms.

The numerical calculation process has been verified by testing with a synthetic magnetic field time series for which an exact analytic solution of the geoelectric field is available. Comparison of the calculated geoelectric fields with the analytic solution gives effectively perfect agreement (correlation coefficients better than 0.9999).

Uncertainties in calculation of geoelectric fields for critical infrastructure are not a result of the computation process. Instead, they arise from uncertainties in the data used in the calculations. The sparsity of magnetic observing sites means that the magnetic field disturbances experienced by the parts of a power system, pipeline, or other affected infrastructure between observing sites are not necessarily accurately specified. Similarly, Earth transfer functions are based on magnetotelluric surveys that only cover part of a network so they may not exactly represent the Earth conductivity structure throughout the entire system. When more magnetic data and Earth conductivity information become available they can be used with the methods described here to calculate better values of the geoelectric fields affecting critical infrastructure.

Acknowledgements

This work was performed as part of the space weather activities of the Public Safety Geoscience program at Natural Resources Canada. Contribution number: 20180404.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Turner, G. (2011) North Pole, South Pole. The Epic Quest to Solve the Great Mystery of Earth’s Magnetism. The Experiment, New York, 272 p.
[2]	Campbell, W.H. (2003) Introduction to Geomagnetic Fields. 2nd Edition, Cambridge University Press, Cambridge.
[3]	Boteler, D.H., Pirjola, R.J. and Nevanlinna, H. (1998) The Effects of Geomagnetic Disturbances on Electrical Systems at the Earth’s Surface. Advances in Space Research, 22, 17-27. https://doi.org/10.1016/S0273-1177(97)01096-X
[4]	Prescott, G.B. (1866) History, Theory and Practice of the Electric Telegraph. IV Edition, Ticknor and Fields, Boston.
[5]	Boteler, D.H. (2006) The Super Storms of August/September 1859 and Their Effects on the Telegraph System. Advances in Space Research, 38, 159-172. https://doi.org/10.1016/j.asr.2006.01.013
[6]	Gummow, R.A. (2002) GIC Effects on Pipeline Corrosion and Corrosion Control Systems. Journal of Atmospheric and Solar-Terrestrial Physics, 64, 1755-1764. https://doi.org/10.1016/S1364-6826(02)00125-6
[7]	Boteler, D.H. and Trichtchenko, L. (2015) Telluric Influence on Pipelines. In: Revie, R.W., Ed., Oil and Gas Pipelines: Integrity and Safety Handbook, Chapter 21, John Wiley & Sons, Inc., Hoboken, 275-288. https://doi.org/10.1002/9781119019213.ch21
[8]	Boteler, D.H. (1994) Geomagnetically Induced Currents: Present Knowledge and Future Research. IEEE Transactions on Power Delivery, 9, 50-58. https://doi.org/10.1109/61.277679
[9]	Bolduc, L. (2002) GIC Observations and Studies in the Hydro-Québec Power System. Journal of Atmospheric and Solar-Terrestrial Physics, 64, 1793-1802.
[10]	Pulkkinen, A., Lindahl, S., Viljanen, A. and Pirjola, R. (2005) Geomagnetic Storm of 29-31 October 2003; Geomagnetically Induced Currents and Their Relation to Problems in the Swedish High-Voltage Power Transmission System. Space Weather, 3, S08C03. https://doi.org/10.1029/2004SW000123
[11]	Guillon, S., Toner, P., Gibson, L. and Boteler, D. (2016) A Colorful Blackout. IEEE Power & Energy Magazine, 14, 59-71.
[12]	Boteler, D.H. (2001) Assessment of Geomagnetic Hazard to Canadian Power Systems. Natural Hazards, 23, 101-120.
[13]	Viljanen, A., Pirjola, R., Wik, M., ádám, A., Prácser, E., Sakharov, Y. and Katkalov, J. (2012) Continental Scale Modelling of Geomagnetically Induced Currents. Journal of Space Weather and Space Climate, 2, A17. https://doi.org/10.1051/swsc/2012017
[14]	Viljanen, A., Pirjola, R., Prácser, E., Katkalov, J. and Wik, M. (2014) Geomagnetically Induced Currents in Europe. Modelled Occurrence in a Continent-Wide Power Grid. Journal of Space Weather and Space Climate, 4, A09. https://doi.org/10.1051/swsc/2014006
[15]	Wik, M., Viljanen, A., Pirjola, R., Pulkkinen, A., Wintoft, P. and Lundstedt, H. (2008) Calculation of Geomagnetically Induced Currents in the 400 kV Power Grid in Southern Sweden. Space Weather, 6, S07005. https://doi.org/10.1029/2007SW000343
[16]	Boteler, D.H. (2015) The Evolution of Québec Earth Models Used to Model Geomagnetically Induced Currents. IEEE Transactions on Power Delivery, 30, 2171-2178.
[17]	Love, J.J., Lucas, G.M., Kelbert, A. and Bedrosian, P.A. (2018) Geoelectric Hazard Maps for the Mid-Atlantic United States: 100 Year Extreme Values and the 1989 Magnetic Storm. Geophysical Research Letters, 45, 5-14. https://doi.org/10.1002/2017GL076042
[18]	Kavanagh, E.R. (1974) Time Sequence Analysis in Geophysics. Third Edition, Univ. of Alberta Press, Edmonton, 492.
[19]	Kaufman, A.A. and Keller, G.V. (1981) The Magnetotelluric Sounding Method. Methods in Geochemistry and Geophysics Vol. 15, Elsevier Scientific Publishing Company, Amsterdam.
[20]	Simpson, F. and Bahr, K. (2005) Practical Magnetotellurics. Cambridge University Press, Cambridge, 272 p. https://doi.org/10.1017/CBO9780511614095
[21]	Chave, A.D. and Jones, A.G. (2012) The Magnetotelluric Method, Theory and Practice. Cambridge University Press, Cambridge, 552 p. https://doi.org/10.1017/CBO9781139020138
[22]	Marti, L., Yiu, C., Rezaei-Zare, A. and Boteler, D. (2014) Simulation of Geomagnetically Induced Currents with Piecewise Layered-Earth Models. IEEE Transactions on Power Delivery, 29, 1886-1893. https://doi.org/10.1109/TPWRD.2014.2317851
[23]	Bloomfield, P. (2000) Fourier Analysis of Time Series: An Introduction. John Wiley & Sons, New York, 2nd Edition, 269 p.
[24]	Boteler, D.H. (2012) On Choosing Fourier Transforms for Practical Geoscience Applications. International Journal of Geosciences, 3, 952-959. https://doi.org/10.4236/ijg.2012.325096
[25]	Bracewell, R.N. (1978) The Fourier Transform and Its Applications. Second Edition, McGraw-Hill Book Company, New York, 444 p.
[26]	Pirjola, R.J. and Boteler, D.H. (2017) Truncation of the Earth Impulse Responses Relating Geoelectric Fields and Geomagnetic Field Variations. Geosciences Research, 2, 72-92. https://doi.org/10.22606/gr.2017.22002
[27]	Boteler, D.H., Pirjola, R.J. and Marti, L. (2019) Analytic Calculation of Geoelectric Fields Due to Geomagnetic Disturbances: A Test Case. IEEE Access, 7, 147029-147037. https://doi.org/10.1109/ACCESS.2019.2945530
[28]	Boteler, D.H. and Pirjola, R.J. (1998) The Complex-Image Method for Calculating the Magnetic and Electric Fields Produced at the Surface of the Earth by the Auroral Electrojet. Geophysical Journal International, 132, 31-40. https://doi.org/10.1046/j.1365-246x.1998.00388.x
[29]	Viljanen, A. and Pirjola, R. (1994) On the Possibility of Performing Studies on the Geoelectric Field and Ionospheric Currents Using Induction in Power Systems. Journal of Atmospheric and Terrestrial Physics, 56, 1483-1491.
[30]	Viljanen, A., Pulkkinen, A., Amm, O., Pirjola, R., Korja, T. and BEAR Working Group (2004) Fast Computation of the Geoelectric Field Using the Method of Elementary Current Systems and Planar Earth Models. Annales Geophysicae, 22, 101-113. https://doi.org/10.5194/angeo-22-101-2004
[31]	Adám, A., Prácser, E. and Wesztergom, V. (2012) Estimation of the Electric Resistivity Distribution (EURHOM) in the European Lithosphere in the Frame of the EURISGIC WP2 Project. Acta Geodaetica et Geophysica Hungarica, 47, 377-387. https://doi.org/10.1556/AGeod.47.2012.4.1