Quantitative inﬂuence of coast effect on geomagnetically induced currents in power grids: a case study

– In recent years, several magnetic storms have disrupted the normal operation of power grids in the mid-low latitudes. Data obtained from the monitoring of geomagnetically induced currents (GIC) indicate that GIC tend to be elevated at nodes near the ocean-land interface. This paper discusses the inﬂuence of the geomagnetic coast effect on GIC in power grids based on geomagnetic data from a coastal power station on November 9, 2004. We used a three-dimensional (3D) Earth conductivity model to calculate the induced electric ﬁeld using the ﬁnite element method (FEM), and compared it to a one-dimensional (1D) layered model, which could not incorporate a coastal effect. In this manner, the GIC in the Ling’ao power plant was predicted while taking the coast effect into consideration in one case and ignoring it in the other. We found that the GIC predicted by the 3D model, which took the coastal effect into consideration, showed only a 2.9% discrepancy with the recorded value, while the 1D model underestimated the GIC by 23%. Our results demonstrate that the abrupt lateral variations of Earth conductivity structures signif-icantly inﬂuence GIC in the power grid. We can infer that high GIC may appear even at mid-low latitude areas that are subjected to the coast effect. Therefore, this effect should be taken into consideration while assessing GIC risk when power networks are located in areas with lateral shifts in Earth conductivity structures, such as the shoreline and the interfaces of different geological structures.


Introduction
Geomagnetic disturbances (GMD) resulting from solar activities give rise to geomagnetically induced currents (GIC), which disrupt the normal operation of the power grid, highspeed railway, and oil and gas pipelines (Pirjola, 1985(Pirjola, , 2000Kappenman & Albertson, 1990;Boteler et al., 1998;Molinski, 2002;Girgis & Vedante, 2012;Liu et al., 2016). Several geomagnetic storms with adverse impacts on power grids have occurred in high-latitude regions such as North America and northern Europe (Kappenman, 1996;Bolduc, 2002;Wik et al., 2008). Notably, in March 1989, a GIC event led to a blackout across the entire Canadian province of Quebec for 12 h, resulting in huge socioeconomic costs (Bolduc, 2002). In recent years, some mid-low latitude areas like Brazil, Australia, Japan and China have also witnessed disruptive GIC incidents (Kappenman, 2004;Gaunt & Coetzee, 2007;Liu et al., 2009bLiu et al., , 2013Pulkkinen et al., 2010) with the expansion of long-distance and large-capacity transmission lines as well as the growing scale of the power grids in these countries. The Ling'ao nuclear power station is located in Shenzhen city, Guangdong, China. On November 9, 2004, the station recorded a maximum GIC amplitude of 75.5 A through the neutral point of its transformer No. 1, which is even higher than that seen in some high-latitude areas. A preliminary analysis of the GMD event and its induction of an unusually large current indicates that in addition to the usual factors -transmission line direction, length, and corner effect of the power grid -coast effects may have increased the amplitude of the GIC (Liu et al., 2011).
The conductivity of seawater is much greater than that of adjacent land. During a geomagnetic storm, the current densities on either side of the land-sea interface differ greatly, which markedly enhances the induced electric field perpendicular to the coast and increases the GIC in power grids situated near coastal areas (Fischer, 1979;Gilbert, 2015). This ''coast effect'' phenomenon has been discussed in the literature. Gilbert (2005) modeled the electric field distribution near the oceanland interface by using a two-region model and demonstrated that the coast effect had a significant impact on GIC in the power grid. Olsen & Kuvshinov (2004) built a thin-shell model and quantitatively demonstrated that the geoelectric field within 50 km of the coast is significantly enhanced during a geomagnetic storm. This model was simplified by Pirjola (2000Pirjola ( , 2013, who confirmed this impact of the coast on GIC. Beggan (2015) produced synthetic models of the auroral electrojet in different locations over UK and investigated the effects of varying the 2D thin-sheet model. Divett et al. (2017) explored the impact of GIC on New Zealand's South Island electrical transmission network by developing a thin-sheet conductance (TSC) model for the region, a geo-electric field model, and a GIC network model. Wang et al. (2016) modeled geomagnetic induction hazards by using a 3D conductivity model of Australia. Fujita et al. (2018) presented a nationwide intensity map of the geomagnetically induced electric field in Japan by using a realistic geomagnetic storm as the induction source, a fine mesh (0.125°· 0.125°), and a realistic 3D heterogeneous ground resistivity model. The results are significant for the GIC calculation in Japan in future. Another study showed that the local enhancement of geoelectric fields resulting from the abrupt lateral change in conductivity between sea water and land can be more than 10 V km À1 , which is significantly greater than the electric field strength of 1 V km À1 that is typically associated with an induced electric field (Horton et al., 2012). Although the coast effect is already a known phenomenon and previous research has hinted that large GICs may occur in coastal areas owing to the amplification of the geoelectric field, the coast effect influencing GIC in a real power grid has yet to be quantified in mid-low latitude areas, such as China. Because the GIE influenced by the coast effect is local, it is still difficult to determine the contribution of the coast effect on GIC.
To determine the contribution of the coast effect on GIC, we have conducted a mathematical study of a geoelectric field induced during a magnetic storm. A one-dimensional (1D) Earth conductivity model proved insufficient to accurately describe the coast effect. We then adopted a three-dimensional (3D) Earth conductivity model to predict the electric field induced during the magnetic storm that occurred at Ling'ao power station on November 9, 2004. The results of a similar study conducted in the 1D layered model were used for comparison. Next, we modeled the GIC through the neutral point of the transformer at the power plant using a GIC-equivalent model of the power grid. From a comparison of the two sets of quantitative results derived from the 1D and 3D models with the recorded data, we were able to show that the coast effect does indeed contribute to amplifying GIC in power grids.
Our results will help operation staff make a more holistic assessment of GIC-associated risk and recognize the danger point in power grids. Furthermore, they will provide a reference that could be used to employ proper defensive measures against storm disaster at various geographical locations.
2 Experimental methods

Basic modeling approaches
The modeling of GIC in power grids can be divided into two separate steps (Pirjola, 2000): (1) ''geophysical part'': the determination of the geoelectric field and (2) ''engineering part'': the computation of the GIC produced by the electric field. GIC modeling approaches have matured over the years (Horton et al., 2012;Butala et al., 2017). The effect of a geomagnetically-induced electric field on a power grid is taken to be equivalent to a set of voltage sources imposed on its transmission lines between various grounded points. The value of the voltage equals the integral of the geoelectric field along the line, which converts the GIC calculation into a circuit problem for step 2.
In earlier studies, the plane wave method was commonly utilized to calculate GIC without considering lateral variation of Earth conductivity structures (Pirjola, 2002;Wik et al., 2008;Liu et al., 2009b;Watari et al., 2009). As a result of the low frequency of magnetic storms, the electromagnetic field can achieve a penetration depth of several hundred kilometers. At this depth, lateral variations in conductivity on a small scale are so small that the 1D model can be considered a feasible approximation of Earth conductivity structures. Ngwira et al. (2008) verified that the modeling of GICs is effectively improved by the use of a multilayered ground conductivity model. Results show that the newly developed GIC model is only appropriate for the specific geomagnetic station and GIC site pair (Ngwira et al., 2009). However, as the scale of the average power grid continues to expand, Earth conductivity structures underneath the power grid will increasingly exhibit diversity and complexity. In this case, the 1D model becomes progressively incapable of incorporating the large electric field variations found at interfaces such as coastlines and can no longer accurately represent the GIC distribution. Gannon et al. (2017) indicated that the use of 3-D conductivity in electric field estimates can result in large differences in field magnitude and direction. For the purpose of accurately studying the coast effect on GIC in coastal power grids, we developed a 3D conductivity model based on the previously described finite element method (FEM) . Only in a few geometric models can we solve the ground electromagnetic field using analytic method. While numerical methods, especially FEM that is widely used at present, are very suitable for dealing with the problems of complex physical distribution and model boundaries. It has been demonstrated that the results from FEM and the results of analytical expressions or Fast Hankel Transform (FHT) calculations are identical to within 1%, when a sufficiently large model was chosen for FEM with a sufficiently good resolution of elements (Dong et al., 2013). In addition, the program of FEM is strongly universal. So using FEM is convenient to solve the electromagnetic field in various media.
In practical calculations, it is assumed that the source of a disturbance is the surface current at a certain height in the air. Because the conduction current in the Earth is far greater than the displacement current, we neglected displacement current density and considered only conduction current density for the simulation in the eddy current field. As the weighted residual principle is aimed at discretizing the differential equation for algebraic equations, we used the Galerkin method to perform the problem. The sequence of weight functions and of primary functions were made identical. When the potential function is taken as an unknown quantity, the basic equations in the time domain can be expressed as : where the Earth conductor is taken as an isotropic structure; l 0 is the vacuum permeability; V is the whole solving domain; S is the field boundary; w is the scalar weight function and is identical to the primary function; e x , e y , and e z are the unit vectors of the rectangular coordinate system; W represents three-vector weight functions, which are obtained by multiplying the unit vectors by w (Biro & Preis, 1989); e n is the unit normal vector at the interface; A and u, respectively, express vector potential and scalar potential; r is the Earth's conductivity; and K is the surface current density.
Using the interpolation method, we built the approximate functions of A and u. These were combined with the basic equations to solve A and u for every node. Given these values, the electric field strength can be easily derived.
Suppose that the shape function on each unit is N i (i = 1, 2,. . .n; n is the total number of nodes in each unit), the interpolation functions of A and u for each unit are expressed in equations (3) and (4): In the calculation of the three-dimensional field, the hexahedral eight-node element is selected to create the domain mesh. The shape function is where i = 1, 2, 3,. . .,8. The detailed discretization process is presented in Wang (2011). n, g, and f are the local coordinates of the standard cell. When solving the field distribution, we paid greater attention to the changes in electric field at interfaces between adjacent conductivity structures. The interfaces needed to meet the following conditions:

3D Earth conductivity model
Since boundary conditions are more complicated when the boundary is located at points of lateral variation in Earth's conductivity, the 3D ''center to periphery'' model, assumes the boundary is far enough from the solution domain to ensure accuracy and speed in calculation. The central region of the 3D model represents the actual conductivity structure, whereas at the periphery, the model is simplified to a 1D conductivity layered structure that does not take lateral differences into consideration.
From the 3D model, the mathematical calculation of the geoelectric field was performed using FEM. Model boundary conditions were determined as follows.
The ground was set as the upper boundary. The tangential component of the intensity of the magnetic field acting on it was obtained from the geomagnetic observatory. Surface current density K was applied to the model. The boundary conditions are represented in equations (10) and (11). Equation (10) represents the second boundary condition: where 1 l 0 r Â A ¼ H is the strength of the magnetic field on the ground and e n is the outer normal vector of the ground with e n = Àe z . We assume that B x and B y are the northern and eastern components of the geomagnetic field, respectively, where B x = L cos D and B y = L sin D; L and D are the horizontal components of the geomagnetic field and magnetic declination, respectively.
Thus, based on (10), the southern and eastern components of the surface current density, ÀK x and K y , can be obtained and can together be taken as the boundary condition for the upper boundary. Equation (11) indicates that only a tangential component of the surface current density of the conductive medium exists.
The lower boundary is determined by the penetration depth of the electromagnetic field, where the electric field strength is assumed to decrease to zero. The boundary condition is shown in (12): The boundary conditions around the sides of the model depend on the upper boundary conditions. For example, ground current density K x , corresponds to the eastern component B y .
There are no lateral variations in conductivity around the boundaries. For this reason, telluric current is induced only in the direction of Àx, parallel to the eastern and western borders. The boundary condition for the east-to-west boundary, to which the magnetic field is perpendicular, is expressed in (13): The induced telluric current flows vertically in or out through the northern and southern boundaries, meaning the northern and southern boundary conditions are parallel to the magnetic field. This is represented in (14).
If the ground current density K y is applied to the upper boundary, the method to determine boundary conditions is similar to that described above.

Power grid model to study GIC at Ling'ao nuclear power station
Guangdong Province is located in the developed coastal region of China. This region is dominated by a 500-kV power grid that includes long transmission lines with small resistances; therefore, it is prone to large GIC during magnetic storms. Figure 1 depicts the power grid model that we used to study the GIC at Ling'ao nuclear power station. The model covers 21 nodes and 22 transmission lines including 5 singlecircuit and 17 double-circuit lines. The equivalent parameters of the 500-kV power-grid node and transmission line in the part of Guangdong Province used for GIC calculation are listed in Tables 1 and 2, respectively. The frequencies of the GIC are in the range 0.001-0.01 Hz. In the context of a power grid, GIC can be treated as quasi-direct current (DC). The total resistance in the circuit depends on the resistances and number of transformers, the resistances of the shunt reactors, and the grounding resistances of substations. For the same voltage class, we used reference values of typical DC parameters in power grids. GIC can be blocked by a series-compensated capacitor because of its quasi-DC characteristics. Therefore, the substation with a blocker can be regarded as the boundary of the model. Since, in the Guangdong and Guangxi power grid, the transmission lines from Hezhou to Luodong and from Maoming to Jiangmen are equipped with series-compensated capacitors, Hezhou and Maoming were set as boundaries.

GIC calculation from a 1D model of Earth conductivity structures
The 1D layered Earth conductivity model for Guangdong Province based on the geological structure and data relating to magnetotelluric (MT) sounding (Wenlu et al., 1987;Xuecheng, 1996) is summarized in Table 3.
We obtained geomagnetic data from the Zhaoqing Geomagnetic Observatory (112.30°E, 23.10°N), which is not very far from Ling'ao nuclear power station (114.58°E, 22.64°N). The records were made between 12:49:12 and  22:52:22 on November 9, 2004. Each geomagnetic data point was sampled over a period of 1 s and was chosen so as to include the moment at which the current reached its peak amplitude. Data were collected over one hour and used to calculate the geomagnetic field. Both the geomagnetic data and the layered Earth conductivity model were input into the calculation program written using the MATLAB platform. The results are shown in Figure 2.

GIC calculation from a 3D model of Earth conductivity structures using FEM
The GIC of a node is mainly affected by its adjacent node (Pirjola, 2005); consequently, appropriately shrinking the size of a model can improve its computational efficiency. In order to achieve this, we selected seven nodes around the Ling'ao nuclear power station. Neither the local plane wave method (Marti et al., 2014) nor the GIC calculation method for a uniform geoelectric field can accurately reflect the enhancement of geoelectric field at the interface between different conductivity structures. Therefore, we used a 3D FEM to explore GIC distribution near the coastline. Li et al. (1986) recorded observation data of 14 MT points along three profiles including Jiaoling-Raoping, Zijin-Haifeng, and Zengcheng-Shenzhen in Guangdong Province. Liu et al. (2013) obtained MT data from Yunfu and Jieyang in Guangdong Province, respectively, and established a one-dimensional geoelectric structure model by using the inversion method. Based on the MT sounding data and geological structures in the above references, we established a 3D conductivity model for coastal areas that takes lateral shifts in conductivity structure into consideration, as shown in Figures 3 and 4. The tangential component of the geomagnetic field was obtained from geomagnetic observatories. Since the impacts of GIC are associated with only the Earth and its surface, the 3D conductivity structure can be regarded as a closed region with the Earth's surface as its boundary. In this case, the air region and the external source need not be modeled, which greatly reduces computational cost . The modeling region is shown in Figure 3 and extends for about 684 km in the north-south direction and 881 km in the eastwest direction. The line with the double arrowhead represents the shoreline. The coastline approximately moves along the 45°direction from the southwest to northeast. We selected Raoping, Jiesheng, Kuiyong, and Zhanjiang as the demarcation points and connected them to the shoreline. The land, which consists of three different geoelectric structure parts, is in the north, and the dotted line represents the boundary between different structures. The sea is in the south and its conductivity is uniform, whereas the geoelectric structure under the sea is the same as that of the land adjacent to it. The depth of sea is taken as 1.2 km. Based on measured data and a calculated value for penetration depth, we set the lower border at 400 km underground (Keding, 1994). The yellow markers indicate the location of the MT point, and the red marker indicates the Ling'ao substation.
In Figure 4, we built the standard conventional Cartesian geomagnetic coordinate system, in which the x-axis, y-axis, and z-axis represent the northward, eastward, and downward directions, respectively. Different conductivities can exist in a plane parallel to xoy, representing the variation in the lateral direction. Different conductivities in the z direction of the planes parallel to yoz represent the variation with depth. The slash denotes the coastline. The shaded areas represent the sea. The grid indicates different geoelectric structures. The ground forms the upper boundary. The boundaries of the model are far enough from the Ling'ao nuclear power station. The induced currents are parallel or perpendicular to the boundary surface. Note that Figure 4 is a schematic diagram of the model, and the grid represents different values of conductivity. It is not a practical mesh generation.
The uniqueness theorem of electromagnetic fields (Sheng, 1991;Yang, 2002) shows that field distribution is uniquely determined by the basic equations in the region and boundary conditions. Equations (1) and (2) show that the equations in the field are eddy-current equations, in which conductivity is obtained from measurement. For boundary conditions, the upper boundary condition is derived from the actual variation of geomagnetic field; the lower boundary is the farthest boundary that a changing electromagnetic field can reach, and the electric field strength decreases to zero considering the skin effect. The boundary conditions are consistent with the behavior of current and field at the boundary. The computational domains were distributed into eight-node hexahedral elements by FEM, which is a mature numerical method used for complicated current systems and/or spatially varying conductivity models (Dong et al., 2013). To guarantee computational precision, the meshes were denser near the interface and at the C. Liu et al.: J. Space Weather Space Clim. 2018, 8, A60 surface, where large changes in conductivity occur, and sparser in the regions where conductivity is nearly uniform. All of these factors can ensure the reliability of the solution.
3 Results 3.1 GIC calculated from the 1D layered Earth conductivity structures Figure 2 shows the calculated results obtained from the 1D model and the measured values of the GIC through the neutral point of the transformer in the Ling'ao substation between 22:00 and 23:00 on November 9, 2004. The curves formed by calculated values measured values have approximately similar trends. The maxima for both the recorded data and model appear at 22:50:35. However, the calculated values are smaller than the measured values, particularly at peak values. The recorded total peak value of GIC flowing through the neutral point of Ling'ao Station was 114.33 A, whereas the value calculated using the layered-earth conductivity model was 88.1 A over three phases. This calculated value was 23% lower than the measured value. Based on this preliminary analysis, we postulated that the coast effect may be the source of this discrepancy.
3.2 Electric field distribution at the moment of maximum GIC in Guangdong Province from the 3D Earth conductivity structures Figure 5 shows the partial electric field distribution in Guangdong Province at the moment of maximum GIC at Ling'ao nuclear power station between the 9th and 10th of November, 2004. The electric field in the high-resistance area (the land side) appears distinctly distorted. There is a significant increase in magnitude of the electric field strength and its direction is perpendicular to the coastline. The results confirm that the lateral differences in Earth conductivity structures greatly influence geoelectric field distribution, and that this is particularly evident near the coastline. This coast effect is the likely cause of the significantly elevated GIC at Ling'ao nuclear power plant.
3.3 Electric potential in stations near Ling'ao power station Figure 6 depicts the electric potential of each station relative to the Ling'ao power station. Although the curves vary similarly over time, their respective peak values at 90 s are obviously different. The maximum value of 132 V appears at the Dongguan substation. Due to the same recorded geomagnetic data, the electric potential at different substations follows a similar trend. However, changes in conductivity structures, length of transmission lines and angle between transmission lines, and the direction of the electric field affect their electric potential values and cause strong variations in peak times. It can be seen from Figure 5 that the direction of Dongguan-Ling'ao is similar to that of the electric field, which explains the large electric potential difference.
We then coupled the potentials above with the model of the power grid to obtain GIC at Ling'ao nuclear power station. To ensure accuracy, we set the same boundary conditions of the power grid that were applied to the simplified model including the GIC for transmission lines Zengcheng-Shipai, Huizhouer-Shipai, Shantou-Huizh, Shajiao-Zhongshan, and Shajiao-Guang'nan at 180 s. The line GIC were calculated by the method described in Section 3.1.

Overall comparison of measured and calculated GIC
To highlight the influence of the coast effect, the three results are plotted together in Figure 7. The curves express the GIC calculated from the 3D model of Earth conductivity structures, the GIC calculated from the 1D model of Earth conductivity structures and the measured GIC. Figure 7 shows that the peak times for all predictions are near 90 s, but the values vary. The GIC calculated from the 3D model is closer to the measured data than the GIC calculated from the 1D model. The maximum error of À3.32 A for the former comparison is only À2.9% of the measured peak GIC of 114.33 A, in contrast to the 23% error in the latter comparison. This result indicates that coast effect does indeed make an important contribution to the GIC in power grids.

Discussion and conclusion
We built 1D and 3D models of Earth conductivity structures to calculate the distribution of the electric field near ocean-land interfaces and GIC at Ling'ao nuclear power station during the magnetic storm of November 9, 2004. The 3D simulation results closely coincide with the measured data and the accuracy of the 3D model is 20.1% greater than that of the 1D model. This validates the proposition that the coast effect increases the geoelectric field near the seashore and the value of GIC in the power grids. When an objective power grid   Fig. 6. Electric potential of each station relative to that at the Ling'ao substation.
C. Liu et al.: J. Space Weather Space Clim. 2018, 8, A60 is situated near the shoreline without considering the coast effect, the predicted GIC in the power grid will be underestimated and the potential danger point will be ignored. Substations located in these areas should use transformers with a strong anti-GIC capability.
In China, because of the complex tectonics and greatly varying geological structures, the lateral shifts in Earth conductivity structures are an important factor that cannot be neglected during GIC risk evaluation. The 1D conductivity model is no longer sufficient for estimating GIC and assessing related risks in a power grid. Instead, 3D Earth conductivity models should be adopted. The present results will deepen our understanding of GMD-related disasters occurring in power grids located in mid-low latitudes, and provide a strong analytical approach for the assessment of risk as well as design of prevention measures. Our results are also of great significance to the investigation of factors that affect the ground currents of high-voltage direct-current transmission-line grounding electrodes, a study that would require more elaborate MT sounding data.