© D. Boteler, Published by EDP Sciences 2014

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 mis-operation 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 province-wide 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 stand-alone 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 stand-alone 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₀₁, r₁₂, etc. are the resistances to ground, and r₁, r₂, etc. are the resistances of the transmission lines between substations.

Fig. 1. Single-phase 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)

Rearranging gives(2)

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 off-diagonal 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₁₂, y₂₃, etc. represent the admittances of the transmission lines connecting pairs of nodes and y₁, y₂, 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. Single-phase 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)

Rearranging gives(14)

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 left-hand 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)

Regrouping terms gives(23)

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 off-diagonal 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 three-phase 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 mid-point 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 North-South and East-West 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 North-South and East-West 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 non-spherical 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.

For an ellipsoid Earth, the North-South 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 North-South 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 East-West 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 East-West 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 one-dimensional 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 Lehtinen-Pirjola (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 1-D 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 2-D or 3-D 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 non-planar. 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 non-planar 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 non-planar circuits. Examining the derivation of the Mesh Impedance Matrix method there is no reason, in electrical terms, why the planar or non-planar 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 non-planar 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 1-D 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

All Tables

All Figures

	Fig. 1. Single-phase 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. Single-phase 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 North-South and East-West 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