Abstract

During geomagnetic disturbances, the telluric currents which are driven by the induced electric fields will flow in conductive Earth. An approach to model the Earth conductivity structures with lateral conductivity changes for calculating geoelectric fields is presented in this paper. Numerical results, which are obtained by the Finite Element Method (FEM) with a planar grid in two-dimensional modelling and a solid grid in three-dimensional modelling, are compared, and the flow of induced telluric currents in different conductivity regions is demonstrated. Then a three-dimensional conductivity structure is modelled and the induced currents in different depths and the geoelectric field at the Earth’s surface are shown. The geovoltages by integrating the geoelectric field along specific paths can be obtained, which are very important regarding calculations of geomagnetically induced currents (GIC) in ground-based technical networks, such as power systems.

1. Introduction

During geomagnetic disturbances (GMD), geomagnetically induced currents (GIC) driven by the geoelectric field are flowing in ground-based electrical systems such as electric power transmission networks with neutral grounded transformers. The geoelectric field is the key factor for determining the GIC in power systems, and it can be determined when the sources of the GMD and the Earth conductivity structure are known. For now, it has been usually assumed that the geoelectric field is spatially uniform, which enables a simple calculation of voltages in the transmission lines that drive GIC in the network and through the transformers, for example, [1, 2]. Unfortunately, in the real world, the complexities of ionospheric-magnetospheric source current systems and the inhomogeneities in the Earth’s conductivity structure make the geoelectric field nonuniform [3]. Therefore, the questions arise of how to model and calculate the geoelectric field more accurately and what the exact effects of lateral conductivity variations are on geoelectric fields, on currents flowing within the Earth (called telluric currents), and on GIC.

In geophysical research area, many techniques are used for determining the geoelectric fields and the methods assume different distributions of the source currents and different conductivity structures [4]. As pointed out in [5], here we need to emphasize again that, in electrical engineering research, evaluation of the impacts of GMD on power systems is the primary interest, which is different from the main focus of geophysicists, who mainly want to infer the conductivity profile of the Earth’s interior. Anyway, forward modelling of geomagnetic induction in the Earth introduces many numerical methods and modelling techniques which are good references for research of GMD effects on power systems as well [6].

The methods and techniques to calculate the geoelectric field depend on the model of the Earth conductivity structure and on the GMD source approximation. A widely used model of the Earth is a one-dimensional (1D) conductivity structure, where the Earth is described as a homogeneous or layered semi-infinite conductor with given conductivity values in each layer. The air region is set on the top of the Earth’s surface with the conductivity value equal to zero. The source of GMD is assumed to be a vertically incident plane wave or a line current located at a certain height. In the plane wave case, the geoelectric field occurring at the Earth’s surface can be determined by multiplying the geomagnetic variation by the surface impedance, which is a recursive function including the depths and conductivity values of different layers [7]. This process is known as the “plane wave method” [8]. This technique is invalid if the source differs much from a plane wave [9], which is the situation, for example, for a line current source representing an ionospheric electrojet. In this case, some other techniques, such as the complex image method (CIM) [10] or the Fast Hankel Transform (FHT) method [11], can be used to calculate the electric and magnetic fields at the surface of the Earth. The key point of CIM is that the telluric currents can be approximated by an image of the ionospheric source current located in a complex space, thus enabling closed-form expressions of the electric and magnetic fields at the Earth’s surface. FHT allows fast computations of the integrals involved in the exact expressions of the fields. These approaches bring a lot of results on theoretical understanding of the ground effects of GMD [12].

If the Earth’s conductivity structure is 1D, sophisticated numerical techniques seem to be unnecessary for solving the electromagnetic induction problem because the methods mentioned above can determine geoelectric fields fast and accurately enough. However, the realistic Earth structure is much more complicated than a 1D layered model. One typical structure contains the land-ocean interface and was analyzed in [5, 13, 14]. Although the “coast effect,” which is that the conductivity contrast between land and ocean will enhance the electric field at the landside close to the coastline, has been demonstrated before, simplified assumptions of the ionospheric source make the studies limited. On the other hand, the complicated processes occurring in the near-Earth space during GMD mean that it is impossible to derive detailed functions or formulae for the sources of geomagnetic variations.

In this paper, a typical Earth conductivity structure with lateral variations is modelled for calculating geoelectric fields based on the Finite Element Method (FEM). Due to the features of this structure, two-dimensional (2D) and three-dimensional (3D) implementations are introduced separately. Firstly, the governing equations and boundary conditions in both cases are explained briefly and the numerical results are compared. Then a more complicated Earth conductivity structure which can only be analyzed in 3D is modelled. Some applications and limitations of this approach are discussed at the end.

2. Models and Methods

2.1. Theoretical Background

In electromagnetic calculations related to GMD the displacement currents can be neglected due to the low frequencies involved [15, 16], so calculating the telluric currents and the geoelectric fields becomes a quasistatic field problem. In GMD research, the conductivity structure is generally assumed to be isotropic and the conductivity will have different values in different regions, but in each region the conductivity is assumed to be uniform. The permeability always has its free space value  H/m.

Consider the Earth’s surface as a planar interface between Earth and air. To calculate the geoelectric field, the traditional computational model contains the Earth conductivity region and the air region containing a specific source. Since only the fields at the Earth’s surface and within the Earth are relevant for assessing GIC impacts, the Earth conductivity structure can be seen as a closed region where the Earth’s surface becomes its boundary. Based on the uniqueness of the solution of a boundary value problem [17], only the tangential component of geomagnetic field obtained from geomagnetic observatories in the real world [18] is necessary for determining the electromagnetic fields in the Earth. Under this circumstance, the air region and the external source need not be modelled, which reduces the computation.

To evaluate the approach presented above, a complex conductivity structure with lateral conductivity variations is considered, as shown in Figure 1. A Cartesian coordinate system is used, where points northwards, points eastwards, and points downwards. The Earth’s surface is set at . The conductivity structure under the Earth’s surface is assembled by two different one-layered models where the Quebec layered model is assumed for and the British Columbia (BC) layered model is assumed for [19]. The Quebec and BC models are typical for a resistive and a conducting ground, respectively. In our calculation, the geomagnetic field at the Earth’s surface is assumed to be uniform and to only have the component. These assumptions make all the field quantities have zero derivatives with respect to the coordinate. Therefore, the Galerkin FEM technique with both 2D planar elements and 3D solid elements can be used to solve this boundary value problem.

2.2. 2D Model and FEM Implementation

In the 2D model, the time variation of the geomagnetic field is assumed to be harmonic with the angular frequency . The superscript “” represents symbols in the phasor form. The governing equation of the 2D quasistatic field problem in terms of the magnetic field in the Earth can be written as where is the component of . In the real world, geomagnetic field which is measured at observatories equals .

As explained above, the Earth’s surface is set as the top boundary of the conductivity structure. Then the top boundary condition is that the magnetic field is assumed to be known and equal to a uniform value . The fields at the bottom boundary of the model are assumed to attenuate to zero. To satisfy this boundary condition, the vertical scale of the model is chosen to be 1000 km deep. The horizontal scale of the model is chosen far away from the horizontal interface , minimizing the effects of lateral changes on the geomagnetic fields at the side boundaries. Therefore, the normal derivatives of the magnetic field at the side boundaries equal zero. All the boundary conditions can be rewritten as where is the normal direction at the side boundaries .

Using the Galerkin weighted residual method, (1) can be written as where are the weighting functions and are equal to the basis functions when and . Here denotes the total number of all nodes in the computational domain. Note that the integrations at the boundaries when discretizing (5) are zero because (4) is the second kind homogeneous boundary condition. So only (2) and (3) are needed to be considered and (4) is automatically satisfied. The weighting functions and the basis functions in the element satisfy and where and are the shape functions and . Here denotes the total number of nodes in each element and can be variable due to different types of elements. In a 2D FEM calculation, four-node quadrilateral-shaped elements are applied to mesh the conductivity regions. The FEM grids in a small region containing the conductivity interface near the Earth’s surface are shown in Figure 2. The skin depth is taken into consideration when deciding the densities of the grids, making it dense enough to give the correct results.

2.3. 3D Model and FEM Implementation

For completing the evaluation of the approach, a 3D model with solid elements is considered, and the sketch is shown in Figure 3. In the 3D model, the Maxwell equations are expressed by the magnetic vector potential and the electric scalar potential to decrease the amount of the unknown quantities. Then the governing equations are

The scale of the model along the -axis is set to be 30 km due to the uniformity in this direction. Then the boundary conditions at the front and back planes are where is the normal vector at the boundaries.

The other boundary conditions are similar to those in the 2D model but are expressed by potentials. Equation (2) can be written as where the surface current density is in the direction and its magnitude equals . Equation (3) can be written as

Since the left and right boundaries are far away from the interface between the different conductivity structures, the conductivity can be regarded as horizontally uniform there. The electric scalar potential at one boundary () is forced to be equal to that at the other (). The magnetic fields are parallel to these boundaries; that is, the tangential component of magnetic vector potential is zero. These boundary conditions can be written as

Using the Galerkin weighted residual method and taking the Coulomb gauge into consideration [20], (6) and (7) lead to where can be replaced by and , which represent the unit vectors in the Cartesian coordinate system. It means that (12) consisted of three individual equations. The integrations at the boundaries are nonzero because (9) is the second kind nonhomogeneous boundary condition and should be considered when implementing the discretization.

Eight-node hexahedral-shaped elements are used in the 3D calculation. Results from the same region as in the 2D calculation are obtained for comparison.

3. Numerical Results

The geomagnetic field at the Earth’s surface is assumed to be harmonic with the period of 100 s and the amplitude of 100 nT to simulate a typical GMD. Based on (2), this changing geomagnetic field can be applied as one of the boundary conditions in the 2D FEM calculation, while for (9)   should be converted to the surface current density which has the same magnitude with but the direction of . All the other boundary conditions are applied in the 2D calculation as (3) and (4) and in the 3D calculation as (8), (10), and (11). Geoelectromagnetic fields at the Earth’s surface and induced telluric currents within the Earth can be directly solved from (5) in the 2D FEM, whereas in the 3D FEM they are obtained after the magnetic vector potential and electric scalar potential are solved from (12) and (13).

The 2D FEM results are shown in an area having a horizontal size of 200 km, that is, 100 km away from the conductivity interface to both directions, and a vertical size from the Earth’s surface to the depth of 150 km. The real parts of the induced telluric currents flowing in the two different conductivity regions are shown in Figure 4. The density of the current lines represents the magnitude of the currents. It can be clearly seen from Figure 4 how the flowing paths of the induced currents change due to the lateral conductivity variation. Referring to the conductivity values shown in Figure 1, it can be concluded that the current lines are denser in conductive layers than those in resistive layers. For example, more telluric currents are induced in the third layer of BC structure than in the first layer of Quebec structure.

To compare the results from the 2D and 3D FEM calculations, contour diagrams of the induced currents and electric fields from these two results are shown in Figure 5. These results agree with each other quite well, which is expected if the mesh is refined enough because the 2D and 3D calculations are applied to model the same problem. The contour diagrams show that the induced currents are continuous while the electric fields are discontinuous at the interfaces of different conductivity regions. The skin effect can be seen from Figure 5 where the induced currents decrease as the depth increases. Based on the values of the induced currents and electric fields shown in Figure 5, it can be seen that the induced currents attenuate to a small value at the depth of 150 km in this calculation but will be different if the frequency of the geomagnetic field variation or the conductivity values change. The skin effect and skin depth should be taken into consideration in the modelling process to determine the depth of bottom boundary. This depth is chosen at five to six times of skin depth making the varying electromagnetic field attenuate to zero.

The electric field on the Quebec side increases dramatically near the Earth’s surface due to the lateral conductivity change. For engineering research, the geoelectric fields at the Earth’s surface are significant, and so results from 2D and 3D calculations are shown in Figure 6. It can be seen that both the real and the imaginary part of the electric field close to the conductivity interface increase over 40% at the resistive side and decrease about 60% at the conductive side. The range of this trend of the electric field extends about 60 km at the resistive side and 30 km at the conductive side away from the conductivity interface. This phenomenon is very similar to the well-known “coast effect” which is described, for example, in [21, 22]. The explanation for the abrupt change of the electric field is that, at the interface of two adjacent conductivity regions and , the continuity of the induced currents requires that . Since , the normal component of electric field at the conductivity interface is discontinuous; that is, . It can be seen from (7) that the gradient of electric scalar potential is nonzero at conductivity interfaces to guarantee the continuity conditions of induced currents.

4. Discussion

As mentioned above, the 2D FEM model can be applied in these calculations because the conductivity structures and the geomagnetic field variations are uniform in the direction. For the cases where the conductivity structures vary in three dimensions, 2D and planar models are invalid. Consider a 3D conductivity structure consisting of three different layered models. For the Quebec layered model is still assumed. The structure for is divided into two individual models where for the BC model is unchanged but for the structure has the same thickness of each layer with BC model but the conductivity value in each layer becomes five times smaller. The scales of the model along the -axis and the -axis are unchanged. The scale along -axis needs to extend far away from horizontal interface and is set between  km and  km. A sketch of the 3D model is shown in Figure 7.

The geomagnetic field is still assumed to be harmonic with the period of 100 s and with the amplitude of 100 nT and to only have the component. All the governing equations and boundary conditions are the same as in Section 2.3. Figure 8 shows the induced currents obtained from four different depths (1 km, 5 km, 10 km, and 20 km) in the range from  km to  km and  km to  km. It can be seen from Figure 8 that the magnitudes of the induced currents are different in different conductivity regions. The skin effect can also be seen in the figure.

The contour diagram of geoelectric field at the Earth’s surface is shown in Figure 9. The interfaces and the conductivity values in the first layer of the three layered models are also shown in the figure. It can be concluded based on the computational results that geoelectric fields increase at the low conductive region close to interfaces of the horizontal changes, for example, the area from  km to . The maximum amplitude of the geoelectric field in Figure 6 is around 885 mV/km, while in Figure 9 it is 1085 mV/km. In other words, the existence of “Smaller BC” model increases the geoelectric field over 20% at the Quebec side. Another interesting phenomenon is that although the conductivity difference between Quebec and BC (0.00005 S/m and 0.002 S/m) is larger than that between Quebec and Smaller BC (0.00005 S/m and 0.0004 S/m), the geoelectric field at the Quebec side in the range from to  km is smaller than that from to  km. It can be explained by the first plot of Figure 8 which shows that the directions of the induced currents around , deviate towards the -axis a little more; that is, the concentration of induced currents increases the geoelectric fields.

The geoelectric fields will generate geovoltages between different locations at the Earth’s surface, which drive GIC in power systems. The geovoltages are obtained by integrating the geoelectric field along specific paths. Results from FEM modelling presented in this paper and from plane wave modelling, which is used, for example, in [23], are compared. Since the geomagnetic field only has component, the geoelectric field calculated by plane wave method only has component. As shown in Figure 10, for results from FEM modelling, five integration paths are chosen where all of them are parallel to -axis but located at different coordinates (100 km, 50 km, 0, −50 km, and −100 km). The coordinates of the starting points and the ending points of the integration paths are set at  km and  km, respectively. While for results from the plane wave modelling, there are only two different results where one is for the area Quebec + Smaller BC and the other is for Quebec + BC.

The geovoltages at side are almost identical for all cases with clear differences occurring from about 20 km to the interface at clearly. At the side, if the location of the integration path is close to the interface at , the geovoltages from FEM modelling are much different from those obtained by plane wave modelling (compare the black solid line, the green solid line, and the green symbol line with the dash line and the dot line). The red solid line and the red symbol line coincide with the plane wave modelling better because they differ to locations distant from the interface at .

Consider a transmission line parallel to -axis located very close to the conductivity boundary at . If it is at , the geovoltage can have a 20% difference by using FEM modelling compared to plane wave modelling, as concluded from the black line and the dash line in Figure 10, while if it is at , the geovoltage difference is 15% by comparing results shown by the black line and the dot line. Generally speaking, horizontal conductivity changes have dramatic influences on geoelectric fields and geovoltages. This should be taken into consideration when evaluating the impacts of GIC on power systems which are constructed close to conductivity boundaries.

The geoelectric field calculated by the plane wave method only has the component, making the geovoltage zero if the integration path is parallel to -axis. However, as shown in Figure 11, the geovoltages are nonzero due to the influences of conductivity variations. The results are obtained by choosing the integration paths parallel to the -axis with different coordinates (−100 km, −50 km, 0, 50 km, and 100 km). It can be clearly seen from Figure 11 that the geovoltages in the direction can reach over 10 V. In other words, transmission lines located close to the boundary at , and extending from  km to  km, can experience voltages of 10 V possibly leading to significant GIC.

The 2D and 3D FEM calculations in Section 2 are implemented at a typical laptop computer. For 2D modelling, the amount of nodes is 16865 and the computational time is around ten seconds. For 3D modelling, the amount of nodes is 23668 and the computational time is several tens of seconds, while for the 3D modelling in this section, because the scale of the model in direction is much larger than that in the 3D modelling in Section 2, the need of computational resource increases dramatically preventing the use of the same laptop computer. The calculation is implemented at a supercomputer and the amount of nodes is 346053. The computational time is over two hours. This is expected due to the number of unknowns being over a million.

5. Conclusions

An approach to model earth conductivity structures for calculating the geoelectric field during GMD is presented in this paper. The main difference of this approach with other techniques is that the air region and the source region need not be modelled. The key step in engineering research of GIC is determination of the geoelectric field which causes voltages at the Earth’s surface. The electromagnetic fields within the Earth can be assured unique and correct due to the uniqueness theorem. To evaluate the boundary conditions, an Earth structure with two-dimensional conductivity changes is modelled and computed in 2D and 3D FEM with different kinds of grids. Then a 3D conductivity structure is modelled and induced currents at different depths, geoelectric field, and geovoltages at the Earth’s surface are obtained. Based on the comparisons of the results from 2D and 3D calculations, several conclusions can be drawn.(1)For a 2D conductivity structure discussed in this paper, when the geomagnetic field is assumed to be uniform along the -axis, both a 2D model expressing Maxwell’s equations in terms of the magnetic field and a 3D model using potentials can be used. The results from these two techniques agree very well. Under these circumstances, a 2D model is more appropriate due to the simplicity and smaller computational requirements.(2)For a 3D conductivity structure, which is more realistic in the real world, only 3D elements can be used. This demonstrates a limitation in the 2D model. Compared to those in 2D conductivity structure cases, the computational time increases significantly in the 3D conductivity structure calculations even though the simulations are run on a supercomputer. This is a limitation in the 3D model.(3)The approach presented in this paper can be a useful tool for modelling conductivity structures with lateral variations. It is clearly shown how the induced currents flow in different conductivity regions and the impact on the electromagnetic fields near the interfaces between the regions is also shown. Taking the highly computational demand into consideration, this technique is more applicable to handle the problems with conductivity changes, while for horizontally uniform structures other methods and techniques are fast and accurate enough.

In the numerical cases of this paper, the geomagnetic field is assumed to be uniform and only have one component. In the real world, the geomagnetic field variations during GMD can be obtained from geomagnetic observatories located at different places all over the world, instead of the simplifying assumption in this paper. The geomagnetic fields from observatories are the total fields which are partly produced by the source currents in the space and partly by induced currents in the Earth. These real-time data are applied as the boundary conditions on the top boundary in (2) and (9). The observed geomagnetic fields are more likely nonuniform making the 2D implementation invalid even for a 2D conductivity structure. In this case 3D models and calculations are more needed.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors thank David Boteler, Donald Danskin, and Larisa Trichtchenko for reviewing and discussing this paper. This work was supported by the National Natural Science Foundation of China (51177045) and the Program of the Co-Construction with Beijing Municipal Commission of Education of China (GJ20120006). The editor thanks two anonymous referees for assisting in evaluating this paper.