Issue 
J. Space Weather Space Clim.
Volume 4, 2014
Space Weather and Challenges for Modern Society



Article Number  A21  
Number of page(s)  11  
DOI  https://doi.org/10.1051/swsc/2014018  
Published online  09 July 2014 
Regular Article
Methodology for simulation of geomagnetically induced currents in power systems
Geomagnetic Laboratory, Natural Resources, Ottawa, Canada
^{*} Corresponding author: dboteler@nrcan.gc.ca
Received:
23
April
2013
Accepted:
19
May
2014
To assess the geomagnetic hazard to power systems it is useful to be able to simulate the geomagnetically induced currents (GIC) that are produced during major geomagnetic disturbances. This paper examines the methodology used in power system analysis and shows how it can be applied to modelling GIC. Electric fields in the area of the power network are used to determine the voltage sources or equivalent current sources in the transmission lines. The power network can be described by a mesh impedance matrix which is combined with the voltage sources to calculate the GIC in each loop. Alternatively the power network can be described by a nodal admittance matrix which is combined with the sum of current sources into each node to calculate the nodal voltages which are then used to calculate the GIC in the transmission lines and GIC flowing to ground at each substation. Practical calculations can be made by superposition of results calculated separately for northward and eastward electric fields. This can be done using magnetic data from a single observatory to calculate an electric field that is a uniform approximation of the field over the area of the power system. It is also shown how the superposition of results can be extended to use data from two observatories: approximating the electric field by a linear variation between the two observatory locations. These calculations provide an efficient method for simulating the GIC that would be produced by historically significant geomagnetic storm events.
Key words: geomagnetically induced currents (GIC) / electric field / modelling / electromagnetism / methodology
© D. Boteler, Published by EDP Sciences 2014
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
During space weather events the energy transfer from the solar wind leads to electric currents in the magnetosphere and ionosphere. The magnetic fields of these electric currents are superimposed on the Earth’s main field and produce the magnetic field variations seen at the Earth’s surface. These magnetic field variations induce electric currents in the Earth and in technological systems involving long conductors, such as pipelines, phone cables and power systems (Pirjola 2002). In power systems these geomagnetically induced currents (GIC) flow along transmission lines and through transformer windings to ground. The GIC in the transformers creates an extra magnetic field in the transformer core that interferes with the operation of the transformer. This can cause overheating of the transformer, generation of harmonics leading to misoperation of protective relays and increased consumption of reactive power that can cause voltage problems (Molinski 2002; Kappenman 2007). In extreme cases these effects can result in damage to equipment or power blackouts.
Power systems problems resulting from GIC were first reported in 1940 (Davidson 1940). Subsequently, GIC effects have been observed on a number of power systems during major disturbances (Albertson & Thorson 1974; Albertson et al. 1974). In March 1989 a magnetic storm caused widespread line trips, damage to transformers in US and UK and a provincewide blackout in Quebec (Allen et al. 1989; Bolduc 2002). Since then GIC effects have been observed during significant space weather events such as the Halloween storm in 2003 (Pulkkinen et al. 2005) and are attributed as the cause of transformer damage in South Africa (Gaunt & Coetzee 2007). Although the occurrence of significant space weather events is random, the vulnerability of power systems to GIC is increasing due to the use of higher voltage, lower resistance, transmission lines and the increased loading of power systems which leaves little reserve capacity for dealing with unusual conditions. This has led to increased concern about the possible impacts that an extreme space weather event could produce on modern power systems and prompted work to understand the geomagnetic hazard to power systems.
To assess the geomagnetic hazard to power systems requires the ability to simulate the geomagnetically induced currents (GIC) flowing in the power network. The first attempts at modelling GIC (Albertson et al. 1981) started by adapting power system modelling tools such as the Electromagnetic Transients Program (EMTP; Dommel 1993). Later, Lehtinen & Pirjola (1985) developed a standalone method of calculating GIC that has been widely used in the geophysics community (Thomson et al. 2005; Wik et al. 2008; Viljanen et al. 2012). Meanwhile the engineering community continued to apply power system analysis tools for calculating GIC (Kappenman et al. 1981; Radasky et al. 2006; Gilbert et al. 2012). However, the methodology used here for calculating GIC is not clear for anyone not familiar with power system analysis and the power industry software packages have considerable overhead that makes them unsuitable for dedicated GIC simulation applications.
This paper examines the methodology behind the power system analysis software and shows how it can be applied to develop standalone software for simulating GIC. Although the power system modelling software is based on established circuit analysis techniques, their exact application to modelling GIC has never been clearly presented. This paper develops the theory for calculating GIC starting from first principles and shows how the methods normally used for AC modelling can be adapted to model the direct current (DC) geomagnetically induced currents. This clarifies the differences involved in modelling GIC and AC power flow including the placement of the voltage sources in the network. This provides the theoretical foundation for simulation of GIC in a power system.
The practical aspects of GIC simulation are then considered, starting with setting up the model including the resistance values for the transmission lines, transformers and substation grounding. Then it is shown how the GIC simulations can be made using electric field values to calculate the voltage sources or equivalent current sources. It is also shown how superposition of modelling results for northward and eastward electric fields can be combined with time series of magnetic observatory data to provide simulation of GIC for specific events. These techniques are illustrated using a benchmark power system model to show how archived magnetic data can be utilised to simulate the GIC that would be produced during extreme events.
2. Power system modelling
The modelling presented here is based on the Mesh Impedance Matrix method and the Nodal Admittance Matrix method that are the standard approaches used for analysis of alternating current (AC) power flow (Guile & Paterson 1977; Stirling 1978). To model GIC in the power system requires knowledge of the transmission line and transformer resistances and of the grounding resistances of the substations. These resistances are then used to construct a network model which is utilized with the electric fields as input to give the GIC in each branch of the network. Here it is shown how the Mesh Impedance Matrix method and the Nodal Admittance Matrix method can be used for modelling GIC. The Mesh Impedance Matrix method is presented first as there is a more obvious connection between the electric fields and the voltage sources in the lines used in the modelling. It is then shown how the voltage sources can be converted into equivalent current sources and used in the Nodal Admittance Matrix method to calculate GIC.
2.1. Mesh impedance matrix method
In the Mesh Impedance Matrix method the power network is considered as a mesh of loops and a matrix solution is obtained for the current around each loop. The primary element in the mesh impedance network is the loop and the numbering convention is based on the loop number. i_{p} is the current around loop p. r_{p} is the resistance of the transmission line in loop p, and r_{pq} is the resistance common to loops p and q (i.e., the resistance through the transformers at a substation to ground). e_{p} is the induced emf in series with the transmission line impedance in loop p.
2.1.1. Example of a simple network
The derivation of the Mesh Impedance Matrix method can be most clearly illustrated by reference to a simple circuit as shown in Figure 1 where r_{01}, r_{12}, etc. are the resistances to ground, and r_{1}, r_{2}, etc. are the resistances of the transmission lines between substations.
Fig. 1. Singlephase diagram of a power network illustrating mesh resistances, r, induced emfs, e and loop currents, i. 
The starting point for the derivation is Kirchhoff’s voltage law that the algebraic sum of the voltages around any loop is zero. Applying this to the network in Figure 1 gives equations for each loop(1)
This can be written in matrix form:(3)where the diagonal elements of the coefficient matrix are the sum of the resistances round each loop and the offdiagonal elements involve the resistances common to two loops. Matrix inversion can then be used to give the loop currents from which we get the currents in the lines and the currents to ground.
2.1.2. General equations
To develop the matrix equations for a general network, consider a loop p in the middle of a network as shown in Figure 2. This has a line resistance r_{p} and resistances to ground at opposite ends of the loop r_{Ap} and r_{Bp}. e_{p} represents the emf in the line and is obtained by integrating the electric field along the length of the line(4)
Fig. 2. Modelling GIC using the Mesh Impedance Matrix method. 
Current i_{p} circulates in loop p and produces a voltage drop in each resistance. This is the only voltage drop in the line resistance r_{p}, however resistances r_{Ap} and r_{Bp} are also part of other loops and currents in those loops contribute to the voltage drops in r_{Ap} and r_{Bp}. Thus to obtain a solution for loop p requires knowledge of the currents in the other loops and it is necessary to obtain a solution that applies simultaneously to all loops in the network – hence the use of matrix methods.
To develop the matrix equations we start by applying Kirchhoff’s voltage law. Thus for loop p (5)where the first term is the voltage drops produced by loop current i_{p} and the second term is the sum of contributions to voltage drops produced by currents in other loops that share a resistance to ground with loop p with the sign of i_{q} depending on whether the current is in the same or opposite direction to current i_{p} at the shared resistance. In the example shown in Figure 2, r_{1p} and r_{2p} have the value r_{Ap}, r_{sp} and r_{tp} have the value r_{Bp}. Currents in other loops do not contribute to the voltage drops in loop p so the “shared resistances” for these cases are set to zero.
Equation (5) can be generalised to a matrix equation for all loops(6)where(7) (8)with the sign in Eq. (8) depending on whether a loop current i_{q} is in the same or opposite direction to current i_{p} at the shared resistance. The currents are then found by finding the inverse of the impedance matrix [Z] and multiplying by the input voltage sources:(9)
The currents in the lines are given directly by the mesh currents. The currents through a particular resistance to ground are the sum of all mesh currents sharing a path through the resistance to ground.
2.2. Nodal admittance matrix method
In the Nodal Admittance Matrix method the power network is considered as nodes connected together and to ground. Driving voltages (emfs), e, are converted to equivalent current sources, j, as shown in Figure 3. For a voltage source e in series with impedance z, the corresponding equivalent circuit has admittance, y = 1/z, and current source, j = e/z in parallel. A matrix solution is then obtained for the voltage of each node. The node voltages are used to obtain the currents throughout the network.
Fig. 3. Equivalence of series voltage source and parallel current source. 
2.2.1. Example of a simple network
The Nodal Admittance Matrix method can be illustrated by considering Figure 4. Nodes are labelled 1–5 and represent the high voltage connection (bus) at individual substations. The admittances between nodes, y_{12}, y_{23}, etc. represent the admittances of the transmission lines connecting pairs of nodes and y_{1}, y_{2}, etc. represent the admittances of the connection to ground at each substation. Figure 4 represents the same circuit as shown in Figure 1 but with all the components redrawn in terms of admittances and equivalent current sources.
Fig. 4. Singlephase diagram of a power network illustrating nodes and admittances, y, and equivalent current sources, j. 
For the nodal admittance network the starting point for the derivation is Kirchhoff’s current law that the algebraic sum of the currents entering any node is zero, i.e., sum of currents entering on transmission lines equals current flowing to ground. Applying this to the network in Figure 4 gives equations for each node(10)
The current in a line is determined by the current source, the voltage difference between nodes at the ends of the line, and the admittance of the line(11)
Also, using Ohm’s law, the current in the connection to ground can be written in terms of the nodal voltage(12)
Making these substitutions the equations above for each node now become(13)
The current sources in parallel with the line from node k to node n can be represented by “nodal current sources” at each end of the line: one pulling current out of node k and one injecting current into node n. The currents on the lefthand side of Eq. (14) represent the total of the equivalent source currents directed into each node, which we term J_{k}. Thus we can rewrite the equations in the form:(15)
These equations can be written in matrix form(16)
Matrix inversion can be used to solve for the nodal voltages which are then substituted into Eq. (11) to obtain the currents in the lines and also combined with the admittances to ground Eq. (12) to give the currents flowing to ground.
2.2.2. General equations
To develop the general equations for the Nodal Admittance Matrix method consider nodes k and n in the middle of a network as shown in Figure 5. This circuit is identical to that shown in Figure 2, except voltage sources e have been replaced by their equivalent current sources j and impedances have been replaced by the equivalent admittances y. y_{kn} represents the admittance of the transmission line between nodes k and n (note that y_{kn} = y_{nk}), and y_{k} and y_{n} represent the admittances to ground from nodes k and n, respectively.
Fig. 5. Modelling GIC using the Nodal Admittance Matrix method. 
Applying Kirchhoff’s current law we can write an equation for any node k of the form(17)
The current in a line is determined by the current source, the voltage difference between nodes at the ends of the line, and the admittance of the line, given by Eq. (11).
Substituting this into Eq. (17) gives(18)
We make the further substitution(19)where J_{k} is the total of the equivalent source currents directed into node k. Thus we obtain the equation(20)
This equation involves both the nodal voltages, v_{k}, and the currents to ground from each node, i_{k}, as unknowns. The nodal voltage v_{k} is related to the current to ground i_{k} by Ohm’s law so we can substitute for either v_{k} or i_{k} to obtain equations involving only one set of unknowns. In this derivation we make the substitution(21)
Substituting for i_{k} gives equations involving only the node voltages v_{k} as the unknowns:(22)
This can be written in matrix form(24)where the column matrix [V] contains the voltages v_{n} and [Y] is the admittance matrix in which the diagonal elements are the sums of the admittances of all paths connected to node k, and the offdiagonal elements are the negatives of the admittances between nodes k and n, i.e.,(25) (26)
The voltages of the nodes are then found by taking the inverse of the admittance matrix and multiplying by the nodal current sources:(27)
These node voltages can be substituted into Eq. (11) to give the currents in the branches and into Eq. (21) to give the currents to ground from each node.
3. Practical calculations
The above section has shown that both the Mesh Impedance Matrix method and the Nodal Admittance Matrix method can be used for modelling GIC in a power system. The Mesh Impedance Matrix method is conceptually simpler, but the greater computational efficiency of the Nodal Admittance Matrix method makes this the preferred method for simulation of GIC. To perform practical calculations using this theory requires setting up the model of the power system and then calculating the sources to use in the transmission lines for the specified electric fields.
3.1. Setting up the model
Because of their low frequencies (<1 Hz), GIC are usually considered as varying DC currents in the power system and GIC models have only considered the DC resistances of the system. AC power networks utilise threephase power transmission with conductors and transformer windings for each phase connected in parallel. GIC flowing in the three phases meet at the neutral point of the transformer and share a common path through the substation grounding resistance to Earth. The network model can be set up to either calculate the GIC in all three phases or the GIC in each phase individually. In either case it is assumed that the parallel paths of the different phases have identical resistances. Thus the “all phases” model can be constructed by combining the three parallel phase transmission lines resistances into one value (dividing the single line resistance by three) and combining the parallel transformer winding resistances into one value (dividing the single transformer winding resistance by three) and using the substation grounding resistance. The equivalent “single phase” model gives values that are all three times the “all phases” values, i.e., the transmission line and transformer winding resistances per phase and three times the substation grounding resistance.
The actual resistance values vary from system to system, however, there are some general characteristics as shown by Zheng et al. (2014). Higher voltage lines use bundles of conductors for each transmission line resulting in lower resistances. Similarly transformers designed to work at higher voltages handle greater power transfers and are designed with lower resistance windings to reduce the power losses in the transformer. Substation grounding resistance values depend on the soil conditions where the substation is built but are designed to have low values for safety reasons. Neighbouring networks can also influence the GIC in a power system but can be modelled by using a suitable equivalent circuit for the neighbouring network (Boteler et al. 2013). All of these values are used to construct the admittance matrix for the power network.
3.2. Calculating sources in the transmission lines
The other input needed for the Nodal Admittance Matrix calculations is the set of equivalent current sources into each node (Eq. (19)). These are the sum of the equivalent current sources in each line derived from the voltage sources obtained by integrating the electric field along the line (Eq. (4)). When the electric field is uniform or varies linearly with distance along the line this reduces to multiplication of the electric field at the midpoint of the line by the line length. Because the calculations from magnetic field components yield the northward and eastward components of the electric field, it is easier to use these to calculate the voltages produced by these components and then sum these component voltages to give the voltage in the line for the particular direction of the electric field:(28)
Thus calculation of the voltage sources requires calculation of the NorthSouth and EastWest distances spanned by a transmission line.
Calculation of distances on a spherical earth is more complicated than calculation of distances on a flat surface. Consider a transmission line between substations A and B as shown in Figure 6. Assuming a spherical earth, the NS distance is simply calculated from the difference in latitude of substations A and B. However, there is no similar simple relationship for the EW distance because, as shown in Figure 6, lines of longitude converge as they approach the pole. Consequently it is necessary to take into account the latitude of the substations when converting their longitudinal separation into a distance.
Fig. 6. Transmission line between substations A and B superimposed on lines of latitude and longitude. The NorthSouth and EastWest distances between the substations are L _{ N } and L _{ E } respectively (from Horton et al. 2012). 
To get more accurate values and to be consistent with substation latitudes and longitudes obtained from GPS measurements it is necessary to take account of the nonspherical shape of the Earth. The Earth is an ellipsoid with a smaller radius at the pole than at the equator. The precise values depend on the Earth model used. Horton et al. (2012) recommend the WGS84 model (Table 1) which is utilised in the GPS system.
Parameters of the WGS84 Earth model.
For an ellipsoid Earth, the NorthSouth distance is given by:(29)where M is the radius of curvature in the meridian plane and is given by(30)
Substituting in the values from Table 1 gives the expression for the NorthSouth distance in km:(31)where ∆lat is the difference in latitude (degrees) between the two substations A and B, and ϕ is the average of the latitudes of A and B:(32)
Similarly the EastWest distance is given by:(33)where N is the radius of curvature in the plane parallel to the meridian as defined in Eq. (34) and ∆long is the difference in longitude (degrees) between the two substations A and B(34)
Substituting the values from Table 1 gives the following expression for the EastWest distance in km:(35)
The lengths from Eqs. (31) and (35) are then used in Eq. (28) to give the voltage source in the transmission line to be used in the GIC modelling.
4. Simulation of GIC
To illustrate the application of the GIC modelling for simulation of GIC it is convenient to use an example network as shown in Figure 7. This was designed to be used as a benchmark model for testing GIC modelling software, and it includes many features found in real power systems (Horton et al. 2012). Here a distinction is made between “modelling GIC” produced by a specified electric field and “simulation of GIC” produced by a magnetic disturbance: the simulation requiring the use of magnetic field recordings with the GIC modelling to produce GIC values for the duration of the event.
Fig. 7. Map of sample network. Note: Substation 1 and Substation 7 have no connection to ground (from Horton et al. 2012). 
Magnetic observatory data can be used with suitable Earth conductivity models to calculate the electric fields experienced by the power network. There are well established techniques for making these calculations in the time domain or frequency domain using a “uniform field” approximation with onedimensional models of the Earth conductivity structure (Ward & Hohmann 1988; Boteler 1994; Pirjola 2002). There are also more sophisticated methods available for calculating the electric fields produced by more complex source fields (Viljanen et al. 2004). However these require special studies and many power systems do not have sufficient data available to enable their use. Here an examination is made of the GIC simulation that can be performed with magnetic data from one or two magnetic observatories.
4.1. Using magnetic data from one observatory
Magnetic observatory data used with a suitable Earth conductivity model will provide calculations of the northward, E_{X}, and eastward, E_{Y}, components of the electric field. These calculated time series of electric field components E_{X}(t) and E_{Y}(t) can be used in the Nodal Admittance Matrix method calculations for every time increment, t, to calculate GIC. However a more convenient method is to make the calculations of the GIC produced by “reference” electric fields of 1 V/km northward and 1 V/km eastward. Using the principle of superposition, these reference results can be combined to give the GIC for any electric field values (Boteler 2013). First the Nodal Admittance Matrix method calculations are made to obtain the GIC values α_{k} for a northward electric field, i.e., E_{X} = 1 V/km, E_{Y} = 0 and then to obtain the GIC values β_{k} for and eastward electric field, i.e., E_{X} = 0 and E_{Y} = 1 V/km. These values are then scaled by the size of the E_{X}(t) and E_{Y}(t) values relative to the reference values to obtain the current i_{k}(t) to ground at each substation.(36)
Figure 8 shows a time series of E_{X}(t) and E_{Y}(t) values calculated from magnetic data at the Ottawa Magnetic Observatory on November 8, 2004, using the E/B transfer function for Quebec, and the GIC to ground calculated for substations 5 and 6 in the network in Figure 7.
Fig. 8. Calculations made using the Nodal Admittance Matrix method and superposition (Eq. (36)) of results for electric field components E _{ X } and E _{ Y } calculated using data from the Ottawa magnetic observatory to give the GIC at substations 5 and 6 in the sample network of Figure 7. 
4.2. Using magnetic data from two observatories
If data are available from two magnetic observatories on opposite sides of a network then linear interpolation of the electric field calculated for each observatory can be used to estimate the electric fields throughout the network. Consider the case of the benchmark network (Fig. 7) with an observatory A to the North and an observatory B to the South. The magnetic data from each observatory is used to calculate the electric field at the observatory locations. With the time variations of these electric fields there are different values at each time step to use for the interpolation, thus it is not immediately possible to use the scaling of reference results as done for one observatory (Eq. (35)). However, although the actual electric field values from observatory A and observatory B are different their relative contributions to the interpolated electric field at a given location are always the same. The percentage contribution from the observatory A electric field goes from 100% at observatory A to 0% at observatory B. Correspondingly the percentage contribution from the observatory B electric field goes from 100% at observatory B to 0% at observatory A. This is illustrated in Figure 9 where E_{A} is 100 mV/km and E_{B} is 50 mV/km, which shows that the total interpolated electric field can be expressed as the sum of the contribution from Observatory A and the contribution from Observatory B.
Fig. 9. Relative contributions when using linear interpolation between the electric fields calculated for Observatory A and Observatory B. 
Thus, at any location in the network, the electric field can be written as(37)where f_{i} is the fraction of the distance from observatory B to observatory A.
Using the Nodal Admittance Matrix method, partial solutions can be produced for the electric fields calculated from observatory A and for the electric fields calculated from observatory B. Calculations are first made for a northward electric field (E_{XA} = 1 V/km, E_{YA} = 0 V/km) contribution from observatory A, taking account of the factor f_{i}, to give the voltage sources in each line across the network to give the GIC values at each substation α_{Ak}. Repeating the calculations for observatory A, but for an eastward electric field (E_{XA} = 0 V/km, E_{YA} = 1 V/km) gives the GIC values β_{Ak}. Then calculations are made for the contributions from observatory B, again for a northward electric field (E_{XB} = 1 V/km, E_{YB} = 0 V/km) and an eastward electric field (E_{XB} = 0 V/km, E_{YB} = 1 V/km), utilizing the factor (1 − f_{i}) in each case, to give the GIC values α_{Bk} and β_{Bk}. These “reference” values can then be scaled by the actual electric field components E_{XA}(t), E_{YA}(t) and E_{XB}(t), E_{YB}(t) to give the GIC values(38)
Figure 10 shows electric fields for the calculated from magnetic data at the Poste de la Baleine magnetic observatory to the North and the Ottawa Magnetic Observatory to the South on November 8, 2004, using the E/B transfer function for Quebec, and GIC calculated for substations 5 and 6 in the sample network, using the geoelectric fields calculated from the North and South observatories.
Fig. 10. Calculations made using the Nodal Admittance Matrix method and superposition (Eq. (38)) of results for electric field components E _{ X } and E _{ Y } calculated using data from the Poste de la Baleine magnetic observatory (A) to the North and the Ottawa magnetic observatory (B) to the South to give the GIC at substations 5 and 6 in the sample network of Figure 7. 
5. Discussion
This paper has described the methodology for simulation of GIC in a power network utilizing the Mesh Impedance Matrix method or the Nodal Admittance Matrix method. GIC can also be modelled using the LehtinenPirjola (1985) method. The accuracy of results from all these methods depends on how well the input data used represents the characteristics of the magnetic disturbances, the Earth conductivity and the system resistances. Viljanen et al. (2004), Ngwira et al. (2009) and others have used the spherical elementary current system (SECS; Amm & Viljanen 1999) for interpolating the spatial variation of the magnetic field disturbance across the area of a power system. Magnetic field variations can also be fitted using spherical cap harmonic analysis (SCHA; Haines 1985; Haines & Torta 1994) providing an alternative method for interpolation across a power network. These techniques are appropriate when there are a sufficient number of magnetic recording sites to map the spatial structure of the magnetic field variations.
In many cases there may be only one or two magnetic observatories in the vicinity of a power system which does not provide sufficient data to justify use of sophisticated interpolation techniques such as SECS or SCHA. When data is only available from one magnetic observatory these data are used to calculate the electric fields which are then assumed to be uniform across the power network. When data are available from two magnetic observatories they can each be used to calculate electric fields and linear interpolation is then used to calculate the electric field at any location in the power network. In such cases these approximations are forced on the modelling by the limitations of the available data. Whether or not these are reasonable approximations to a real magnetic disturbance depends on the size of the network and the proximity of the magnetic observatory(ies), as well as characteristics of the magnetic disturbance. These approximations are likely to be more reasonable at mid and lower latitudes, away from the more structured magnetic disturbances experienced at high latitudes. However, situations can occur where the geomagnetic monitoring available is insufficient to allow reliable interpolation for GIC modelling. The maximum acceptable spacing between magnetometers will depend on the spatial characteristics (both latitudinal and longitude variations) of the geomagnetic disturbance and will be different for different types of disturbance and in different locations around the world. This is an aspect of GIC simulation that requires further investigation.
Accuracy of GIC simulations can also be affected by spatial variations in the electric fields due to lateral changes in the Earth conductivity structure. The calculations above used a single Earth conductivity model for the whole area of the power system. Where the power system extends over several geological regions then a different 1D Earth conductivity model can be used for each region and has been found to give reasonable results (Viljanen et al. 2004; Marti et al. 2014). However this approach fails to take into account the enhanced electric fields that occur at conductivity boundaries. These are especially important at boundaries where there is a large conductivity contrast such as at a coastline (Pirjola 2013). Determining the electric fields in such situations requires 2D or 3D modelling (Dong et al. 2013). This is a subject that requires further research in connection with GIC applications.
This paper is about the methodology for modelling GIC in a power system. GIC calculations can be made using either the Mesh Impedance Matrix method or the Nodal Admittance Matrix method. Obtaining the solution in the Nodal Admittance Matrix method involves taking the inverse of the nodal admittance matrix and multiplying by the nodal current sources (27) to give the nodal voltages that are then used to calculate GIC. In contrast the Mesh Impedance Matrix method involves taking the inverse of the mesh impedance matrix and multiplying by the voltage sources (9) to give GIC values. The two matrices involved in these calculations should not be confused. The inverse of the nodal admittance matrix is termed the nodal impedance matrix; however, this is not the same as the mesh impedance matrix. This is because the matrices use different numbering schemes: one is based on mesh number and the other is based on node number. The central part of the calculation is taking the inverse of the appropriate matrix. The size of the mesh impedance matrix is determined by the number of loops in the network and is often considerably larger than the nodal admittance matrix which is determined by the number of nodes. For example the simple circuit in Figure 7 has 13 loops but only 8 nodes. For large power networks this size difference becomes important. The Nodal Admittance Matrix method is thus the preferred choice. Sparse matrix techniques can also be used to speed up the calculations.
Some circuit analysis textbooks (Dorf & Svoboda 2006; Hayt & Kemmerly 1993) use terms from topology to classify networks and make a distinction between a circuit made up of a set of meshes that is “planar”, i.e., it can be drawn on a plane, and a circuit of loops that can be nonplanar. The importance of this for power system analysis is that a planar circuit cannot include lines that cross. When analysing multiple voltages levels of a power system there are going to be situations where a power line at a lower voltage levels crosses under a high voltage line. There are also some situations where lines at the same voltage level cross each other. Thus it is important to know if the GIC modelling techniques can handle lines that cross.
A number of tests have been made to check the ability of the various GIC modelling methods to handle crossing lines. Crossing lines were included as a feature in the GIC benchmark model and tests with the different methods all gave the same results. The circuit analysis textbooks state that Nodal Admittance Matrix method based on Kirchhoff’s current law can handle nonplanar circuits (i.e., with crossed lines). For methods based on Kirchhoff’s voltage law, these textbooks distinguish between two methods saying that loop analysis, but not mesh analysis, can handle nonplanar circuits. Examining the derivation of the Mesh Impedance Matrix method there is no reason, in electrical terms, why the planar or nonplanar classification should matter and it gives identical results to the other methods. The term “Mesh Impedance Matrix” method was introduced long before topological definitions started to be applied to circuit analysis. In the topological classification it should perhaps be called a “loop analysis”. The important thing for GIC modelling is that all the methods work for nonplanar networks that include lines that cross.
The GIC simulation examples presented in this paper have been made for the benchmark model which does not allow comparison with measured GIC values. It is obviously desirable to test GIC simulations by making comparisons with GIC recordings at a number of sites on a power network. Such measurements provide a test of the overall simulation, including both the modelling method and the data used in the model. The development from first principles described in this paper provides the theoretical validity of the methodology for modelling GIC. The accuracy of GIC modelling in real situations is more dependent on the accuracy and resolution of the data used as input in the calculations. These include how well the data from a limited number of geomagnetic observatories capture the spatial structure of the geomagnetic disturbance; the accuracy of the Earth conductivity models used to calculate the electric fields, and the accuracy of the resistance values used in the power network model. Therefore comparison of model results with GIC observations is a test of how good are the inputs to the model rather than a test of the GIC modelling methodology itself.
6. Conclusions
GIC in a power network can be calculated using the Mesh Impedance Matrix method or the Nodal Admittance Matrix method can be utilised for modelling GIC.
The Mesh Impedance Matrix method involves construction of the impedance matrix [Z] for the power network, taking the inverse, and multiplying by the input voltage sources [E] obtained by integrating the electric fields along the transmission lines:(39)
This gives the mesh currents [I] that can be used directly to give the GIC in the transmission lines. The GIC through a particular resistance to ground are then obtained by summing all the mesh currents sharing a path through the resistance to ground.
The Nodal Admittance Matrix method involves construction of the admittance matrix [Y] for the power network. The voltage sources in the lines are converted to equivalent current sources which are combined to give the current source into each node [J]. The voltages of the nodes [V] are then found by taking the inverse of the admittance matrix and multiplying by the nodal current sources:(40)
These node voltages are then used to calculate the currents in the transmission lines and the currents to ground from each node.
Although the Mesh Impedance Matrix method is conceptually simpler, the greater computational efficiency of the Nodal Admittance Matrix method makes this the preferred method for simulation of GIC in power systems.
Superposition of GIC model results for “reference” northward and eastward electric fields can be combined with electric fields calculated using date from one or two magnetic observatories to provide time series of GIC values during a space weather event.
A number of factors influence the accuracy of the GIC simulation. The number of magnetic observatories in the vicinity of a power network may not be adequate to describe the spatial variation of the magnetic disturbance across a large power system. Regional 1D Earth conductivity models provide a good first approximation for calculating electric fields; however, these calculations could be improved by taking account of the localised electric field enhancements that occur at Earth conductivity boundaries, such as at a coastline. Finally, the DC resistance values used to construct the system models will influence the accuracy of the GIC results. Thus the accuracy of the GIC simulations is limited, not by the modelling methodology, but by how well the available input data represent the true magnetic variation, the Earth conductivity, and the system characteristics.
Acknowledgments
This work was performed as part of Natural Resources Canada’s Public Safety Geoscience program with support from the Canadian Space Agency and Hydro One. This work benefitted from many useful discussions with Risto Pirjola, Aidan Foss and Larisa Trichtchenko. The referees comments were also very helpful.
References
 Albertson, V.D., and J.M., Thorson, Power system disturbances during a K8 geomagnetic storm: August 4, 1972, IEEE Trans. Power Apparatus and Systems, PAS93, 1025–1030, 1974. [CrossRef] [Google Scholar]
 Albertson, V.D., J.M. Thorson, and A.A. Miske, The effects of geomagnetic storms on electrical power systems, IEEE Trans. Power Apparatus and Systems, PAS93 (4), 1031–1044, 1974. [Google Scholar]
 Albertson, V.D., J.G. Kappenman, N. Mohan, and G.A. Skarbakka, Loadflow studies in the presence of geomagneticallyinduced currents, IEEE Trans. Power Apparatus and Systems, PAS100, 594–606, 1981. [Google Scholar]
 Allen, J., L. Frank, H. Sauer, and P. Reiff, Effects of the March 1989 Solar Activity, Eos, 14, 1479–1488, 1989. [Google Scholar]
 Amm, O., and A. Viljanen, Ionospheric disturbance magnetic field continuation from the ground to the ionosphere using spherical elementary current systems, Earth Planets Space, 51, 431–440, 1999. [Google Scholar]
 Bolduc, L., GIC observations and studies in the HydroQuébec power system, J. Atmos. Sol. Terr. Phys., 64 (16), 1793–1802, 2002. [Google Scholar]
 Boteler, D.H., Geomagnetically induced currents: present knowledge and future research, IEEE Trans. Power Delivery, 9, 50–58, 1994. [Google Scholar]
 Boteler, D.H., The use of linear superposition in modeling geomagnetically induced currents, Proceedings IEEE Power & Energy Society General Meeting, Vancouver, 21–25 July, 2013. [Google Scholar]
 Boteler, D.H., A.J.C. Lackey, L. Marti, and S. Shelemy, Equivalent circuits for modelling geomagnetically induced currents from a neighbouring network, Paper 002195, Proc. IEEE Power & Energy Society General Meeting, Vancouver, 21–25 July, 2013. [Google Scholar]
 Davidson, W.F., The magnetic storm of March 24, 1940 – effects in the power system, Edison Electric Institute Bulletin, 365–366, 374, July, 1940. [Google Scholar]
 Dommel, H., Electromagnetic Transients Program (EMTP) Theory Book, Bonneville Power Administration, Portland, Oregon, 1993. [Google Scholar]
 Dong, B., D. Danskin, R.J. Pirjola, D.H. Boteler, and Z.Z. Wang, Evaluating the applicability of the finite element method for modelling of geoelectric fields, Ann. Geophys., 31, 1689–1698, 2013. [Google Scholar]
 Dorf, R.C., and J.A. Svoboda, Introduction to electric circuits, 7th ed., John Wiley & Sons, 2006. [Google Scholar]
 Gaunt, C.T., and G. Coetzee, Transformer failures in regions incorrectly considered to have low GIC risk, Proceedings, IEEE Power Tech, Lausanne, July, 2007. [Google Scholar]
 Gilbert, J.L., W.A. Radasky, and E.B. Savage, A technique for calculating the currents induced by geomagnetic storms on large high voltage power grids, Proceedings, 2012 IEEE International Symposium on Electromagnetic Compatibility (EMC), 6–10 August, 2012, pp. 323–328, 2012. [Google Scholar]
 Guile, A.E., and W. Paterson, Electrical power systems, vol. 2, 2nd ed., Pergamon Press, Oxford, 1977. [Google Scholar]
 Haines, G.V., Spherical cap harmonic analysis, J. Geophys. Res., 90, 2583–2591, 1985. [CrossRef] [Google Scholar]
 Haines, G.V., and J.M. Torta, Determination of equivalent current sources from spherical cap harmonic models of geomagnetic field variations, Geophys. J. Int., 118, 499–514, 1994. [CrossRef] [Google Scholar]
 Hayt, W.H., and J.E. Kemmerly, Engineering circuit analysis, 5th ed., McGraw Hill, New York, 1993. [Google Scholar]
 Horton, R., D. Boteler, T.J. Overbye, R. Pirjola, and R.C. Dugan, A test case for the calculation of geomagnetically induced currents, IEEE Trans. Power Delivery, 27 (4), 2368–2373, 2012. [Google Scholar]
 Kappenman, J.G., V.D. Albertson, and N. Mohan, Investigation of geomagnetically induced currents in the proposed WinnipegDuluthTwin Cities 500kV transmission line, Final Report for Research Project 12051, prepared for Electric Power Research Institute (EPRI), Palo Alto, July, 1981. [Google Scholar]
 Kappenman, J.G., Geomagnetic Disturbances and Impacts upon Power System Operation. In L.L., Grigsby (Ed.). The Electric Power Engineering Handbook, 2nd Ed., CRC Press/IEEE Press, pp. 161–1622, Chapter 16, 2007. [Google Scholar]
 Lehtinen, M., and R. Pirjola, Currents produced in earthed conductor networks by geomagneticallyinduced electric fields, Ann. Geophys., 3, 479–484, 1985. [Google Scholar]
 Marti, L., C. Yiu, A. RezaeiZare, and D. Boteler, Simulation of geomagnetically induced currents with piecewise layeredearth models, IEEE Trans. Power Delivery, in press, 2014. [Google Scholar]
 Molinski, T.S., Why utilities respect geomagnetically induced currents, J. Atmos. Sol. Terr. Phys., 64 (16), 1765–1778, 2002. [Google Scholar]
 Ngwira, C.M., L.A. McKinnell, P.J. Cilliers, A. Viljanen, and R. Pirjola, Limitations of the modelling of geomagnetically induced currents in the South African power network, Space Weather, 7, S10002, DOI: 10.1029/2009SW000478, 2009. [CrossRef] [Google Scholar]
 Pirjola, R., Review on the calculation of surface electric and magnetic fields and of geomagnetically induced currents in groundbased technological systems, Surv. Geophys., 23, 71–90, 2002. [Google Scholar]
 Pirjola, R., Practical model applicable to investigating the coast effect on the geoelectric field in connection with studies of geomagnetically induced currents, Adv. Appl. Phys., 1 (1), 9–28, 2013. [Google Scholar]
 Pulkkinen, A., S. Lindahl, A. Viljanen, and R. Pirjola, Geomagnetic storm of 2931 October 2003: Geomagnetically induced currents and their relation to problems in the Swedish highvoltage power transmission system, Space Weather, 3, S08C03, DOI: 10.1029/2004SW000123, 2005. [Google Scholar]
 Radasky, W.A., J.G. Kappenman, and J.L. Gilbert, An overview of the development, validation of large power grid models for geomagnetic storms, Proceedings of the 2006 4th AsiaPacific Conference on Environmental Electromagnetics, 1–4 August, pp. 25–28, 2006. [CrossRef] [Google Scholar]
 Stirling, M.J.H., Power System Control, Peter Peregrinus Ltd, London, 1978. [Google Scholar]
 Thomson, A.W.P., A.J. McKay, E. Clarke, and S.J. Reay, Surface electric fields and geomagnetically induced currents in the Scottish Power grid during the 30 October 2003 geomagnetic storm, Space Weather, 3, S11002, DOI: 10.1029/2005SW000156, 2005. [CrossRef] [Google Scholar]
 Trichtchenko, L., A. Zhukov, R. van der Linden, S.M. Stankov, N. Jakowski, I. Stanislawska, G. Juchnikowski, P. Wilkinson, G. Patterson, and A.W.P. Thomson, November 2004 space weather events: realtime observations and forecasts, Space Weather, 5, S06001, DOI: 10.1029/2006SW000281, 2007. [CrossRef] [Google Scholar]
 Viljanen, A., A. Pulkkinen, O. Amm, R. Pirjola, and T. Korja, and BEAR Working Group, Fast computation of the geoelectric field using the method of elementary current systems and planar Earth models, Ann. Geophys., 22, 101–113, 2004. [Google Scholar]
 Viljanen, A., R. Pirjola, M. Wik, A. Ádám, E. Prácser, Y. Sakharov, and J. Katkalov, Continental scale modeling of geomagnetically induced currents, J. Space Weather Space Clim., 2, A17, DOI: 10.1051/swsc/2012017, 2012. [CrossRef] [EDP Sciences] [Google Scholar]
 Ward, S.H., and G.W. Hohmann, Electromagnetic theory for geophysical applications. In : M.N., Nabighian Ed.. Electromagnetic methods in applied geophysics – Theory, vol. 1, Society of Exploration Geophysicists, Tulsa, 131–311, 1988. [Google Scholar]
 Wik, M., A. Viljanen, R. Pirjola, A. Pulkkinen, P. Wintoft, and H. Lundstedt, Calculation of geomagnetically induced currents in the 400 kV power grid in southern Sweden, Space Weather, 6, S07005, DOI: 10.1029/2007SW000343, 2008. [CrossRef] [Google Scholar]
 Zheng, K., D.H. Boteler, R. Pirjola, L.G. Liu, R. Becker, L. Marti, S. Boutilier, and S. Guillon, Effects of system characteristics on geomagnetically induced currents, IEEE Trans. Power Delivery, 29 (2), 890–898, 2014. [CrossRef] [Google Scholar]
Cite this article as: Boteler D: Methodology for simulation of geomagnetically induced currents in power systems. J. Space Weather Space Clim., 2014, 4, A21.
All Tables
All Figures
Fig. 1. Singlephase diagram of a power network illustrating mesh resistances, r, induced emfs, e and loop currents, i. 

In the text 
Fig. 2. Modelling GIC using the Mesh Impedance Matrix method. 

In the text 
Fig. 3. Equivalence of series voltage source and parallel current source. 

In the text 
Fig. 4. Singlephase diagram of a power network illustrating nodes and admittances, y, and equivalent current sources, j. 

In the text 
Fig. 5. Modelling GIC using the Nodal Admittance Matrix method. 

In the text 
Fig. 6. Transmission line between substations A and B superimposed on lines of latitude and longitude. The NorthSouth and EastWest distances between the substations are L _{ N } and L _{ E } respectively (from Horton et al. 2012). 

In the text 
Fig. 7. Map of sample network. Note: Substation 1 and Substation 7 have no connection to ground (from Horton et al. 2012). 

In the text 
Fig. 8. Calculations made using the Nodal Admittance Matrix method and superposition (Eq. (36)) of results for electric field components E _{ X } and E _{ Y } calculated using data from the Ottawa magnetic observatory to give the GIC at substations 5 and 6 in the sample network of Figure 7. 

In the text 
Fig. 9. Relative contributions when using linear interpolation between the electric fields calculated for Observatory A and Observatory B. 

In the text 
Fig. 10. Calculations made using the Nodal Admittance Matrix method and superposition (Eq. (38)) of results for electric field components E _{ X } and E _{ Y } calculated using data from the Poste de la Baleine magnetic observatory (A) to the North and the Ottawa magnetic observatory (B) to the South to give the GIC at substations 5 and 6 in the sample network of Figure 7. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.