Numerical Calculation of Geoelectric Fields That Affect Critical Infrastructure

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.


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 in-International Journal of Geosciences duced 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.

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 where ε 0 = 8.854 × 10 −12 F/m is the permittivity of free space. Consequently, the displacement current term associated with 2πƒε 0 is negligible and can be dropped from Maxwell's equations. Thus the equations can be written as (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 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

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 ( 1 , , N σ σ and thicknesses ( 1 , , N l l  ). Note that the bottom layer is a uniform half-space, so N l = ∞ .
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 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) 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 For each layer above, Equations (25) and (26) Collecting terms in n S and n R gives Substituting these into Equation (28) and cancelling common terms give n n k l n n n n n n k l n n n n Collecting terms in Formula (37) is a recursive relation for calculating K n at the top surface of layer n from 1 n K + 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 1 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.

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 ob- This can be more simply written where F and F −1 represent the forward and inverse Fourier transform, respectively.

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.

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.

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  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(ƒ)/μ 0 , the surface impedance Z(ƒ ) directly gives the Earth transfer function K(ƒ ) to be applied to the calculation of the geoelectric field.

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 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 rec- 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].
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.

Verification of the Calculations
And the phase is given by 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 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.   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

Calculation Using Tensor Transfer Functions
In the situations considered above where the Earth is approximated by a uni-  In these cases the geoelectric fields and geomagnetic field variations are related by a tensor transfer function in the frequency domain: Thus the northward component of the geoelectric field, E x , can be written Similarly, the eastward component of the geoelectric field, E y , can be written 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(ƒ ) .

Discussion
The above verification tests relate to the process for numerical computation of 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

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 netotelluric 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.