An Approach to Model Earth Conductivity Structures with Lateral Changes for Calculating Induced Currents and Geoelectric Fields during Geomagnetic Disturbances

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


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

Mathematical Problems in Engineering
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 threedimensional (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.

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  0 = 4 × 10 −7 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.

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 (H) in the Earth can be written as where Ḣ  is the  component of H. boundary condition is that the magnetic field is assumed to be known and equal to a uniform value  0 .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 Γ 3 .
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 ,  = 1, 2, 3, . . .,   .
Here   denotes the total number of all nodes in the computational domain.Note that the integrations at the boundaries when discretizing ( 5 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 ,  = 1, 2, 3, . . .,   .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.

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 n 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 K is in the  direction and its magnitude equals  0 .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 i can be replaced by j and k, 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.

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 Again the results agree very well with each other.Close to the horizontal interface of the two conductivity structures, the electric field increases at the resistive side whereas decreases at the conductive side.This performance is similar to the coast effect.The physical explanation is that the discontinuous conductivity will lead to an abrupt change of the electric field to guarantee continuity of the induced currents.
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  1 and  2 , the continuity of the induced currents requires that  1  1 =  2  2 .Since  1 ̸ =  2 , the normal component of electric field at the conductivity interface is discontinuous; that is,  1 ̸ =  2 .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.

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  = −1000 km and  = 1000 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  = −100 km to  = 100 km and  = −100 km to  = 100 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  = 25 km to  = 0.The structure for  <  is the Quebec model, and for  >  it consists of two individual models: BC for  <  and "Smaller BC" for  > ."Smaller BC" has the same layer thicknesses as BC but the conductivity value at each layer is five times smaller; that is, the layer conductivities are [0.0004,0.00134, 0.01, 0.002, 0.000666, 0.002, 0.02, 0.2] S/m.The sizes of the computational model in three dimensions are also shown.
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  = 0 to  = −100 km is smaller than that from  = 0 to  = 100 km.It can be explained by the first plot of Figure 8 which shows that the directions of the induced currents around  = 0,  = 0 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  = −100 km and  = 100 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  < 0 side are almost identical for all cases with clear differences occurring from about 20 km to the interface at  = 0 clearly.At the  > 0 side, if the  location of the integration path is close to the interface at  = 0, 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  = 0.
Consider a transmission line parallel to -axis located very close to the conductivity boundary at  = 0.If it is at  = 0+, 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  = 0−, 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  = 0, and extending from  = −100 km to  = 100 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.

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 Γ 1 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.

Figure 2 :
Figure 2: Top parts of the FEM grids in the 2D model.Different colors represent different conductivity layers.The grid sizes can be varied in different layers but should be dense enough taking the frequency of the changing magnetic field and the conductivity in each layer into consideration.

Figure 3 :
Figure 3: Sketch of the computational model used for 3D FEM implementation.All the conductivity values and layer thicknesses are the same with those in Figure 1.For comparison with 2D FEM results, the model size along -axis is chosen small (= 30 km) along which the conductivity structure and geoelectromagnetic fields are constant.As a consequence, all field quantities remain unchanged in this orientation.For real 3D conductivity structures, the model scale should be large enough due to the fact that the electromagnetic fields are three-dimensional.

Figure 4 :
Figure 4: Real parts of the induced currents flowing in different regions.The currents are represented by the curves with arrows.The horizontal conductivity boundaries are indicated by straight lines.These are results from a 2D FEM calculation.The density of the lines represents the magnitude of the currents.The currents will concentrate to the regions with relatively high conductivity values.The currents induced in the first layer of the Quebec structure are very small due to the low conductivity value of this layer.

Figure 5 :
Figure 5: Contour diagrams of induced currents and electric fields.(a) and (c) are results from 2D FEM.(b) and (d) are results from 3D FEM.These results agree very well with each other.The currents in the relatively high conductive layers are larger.The skin effect can be seen from (a) and (b) where the currents attenuate as the depth increases.The electric field is discontinuous at the interface between the two conductivity structures.

Figure 6 :
Figure 6: Electric field at the Earth's surface.Solid lines represent results from 2D FEM and symbols represent results from 3D FEM.(a) Real part and (b) imaginary part.Again the results agree very well with each other.Close to the horizontal interface of the two conductivity structures, the electric field increases at the resistive side whereas decreases at the conductive side.This performance is similar to the coast effect.The physical explanation is that the discontinuous conductivity will lead to an abrupt change of the electric field to guarantee continuity of the induced currents.

Figure 7 :
Figure7: Sketch of the 3D conductivity structure.The structure for  <  is the Quebec model, and for  >  it consists of two individual models: BC for  <  and "Smaller BC" for  > ."Smaller BC" has the same layer thicknesses as BC but the conductivity value at each layer is five times smaller; that is, the layer conductivities are [0.0004,0.00134, 0.01, 0.002, 0.000666, 0.002, 0.02, 0.2] S/m.The sizes of the computational model in three dimensions are also shown.

Figure 8 :
Figure 8: Induced current density at different depths obtained from a 3D FEM calculation.The magnitudes of the currents are different in the different regions.The flowing paths show that the directions of currents change at the interface between the different conductivity structures.Also the skin effect can be seen from the current values.

Figure 9 :
Figure 9: Geoelectric field at the Earth's surface in three individual layered models.The interfaces of the three models are shown as dash lines.The conductivity values in the first layers of each model are also shown in the figure.

Figure 10 :
Figure 10: Geovoltages along different integration paths at the Earth's surface.The paths are all parallel to -axis with different  coordinates.The starting and ending points are at  = −100 km and at  = 100 km.Solid lines represent results from FEM modelling and dash line and dot line represent results from plane wave modelling.

Figure 11 :
Figure 11: Geovoltages along different integration paths parallel to the -axis with different  coordinates.The starting and ending points of these paths are at  = −100 km and at  = 100 km, respectively.All the results are obtained from FEM modelling.
[15,10,of the computational model used for 2D FEM implementation.Quebec has 5 layers.The layer thicknesses are[15,10, 125, 200, ∞] km and the layer conductivities are [0.00005, 0.005, 0.001, 0.01, 0.333] S/m.BC has 8 layers.The layer thicknesses are[4, 6, 5, 20, 65, 300, 200, ∞]km and the layer conductivities are [0.002,0.0067, 0.05, 0.01, 0.00333, 0.01, 0.1, 1] S/m.Only several top layers of a 2D Earth model with the corresponding conductivity values in each layer are shown.The bottom boundary of the model is set at the depth of 1000 km to satisfy the boundary conditions of the magnetic field being sufficiently attenuated.The top boundary Γ 1 , the bottom boundary Γ 2 , and the side boundaries Γ 3 are also shown in the figure.