The Lehtinen-Pirjola Method Modified for Efficient Modelling of Geomagnetically Induced Currents in Multiple Voltage Levels of a Power Network

. The need for accurate assessment of the geomagnetic hazard to power systems is driving a requirement to model geomagnetically induced currents (GIC) in multiple voltage levels of a power network. The Lehtinen-Pirjola method for modelling GIC is widely used but was developed when the main aim was to model GIC in only the highest voltage level of a power network. Here we present a modification to the Lehtinen-Pirjola (LP) method designed to provide an efficient method for modelling GIC in multiple voltage levels. The LP method calculates the GIC flow to ground from each node. However, 15 with a network involving multiple voltage levels many of the nodes are ungrounded, i.e. have infinite resistance to ground which is numerically inconvenient. The new modified Lehtinen-Pirjola (LPm) method replaces the earthing impedance matrix [Ze] with the corresponding earthing admittance matrix [Ye] in which the ungrounded nodes have zero admittance to ground. This is combined with the network admittance matrix [Yn] to give a combined matrix ([Yn]+[Ye]), which is a sparse symmetric positive definite matrix allowing efficient techniques, such as Cholesky decomposition, to be used to provide the nodal 20 voltages. The nodal voltages are then used to calculate the GIC in the transformer windings and the transmission lines of the power network. The LPm method with Cholesky decomposition also provides an efficient method for calculating GIC at multiple time steps. Finally, the paper shows how software for the LP method can be easily converted to the LPm and provides examples of calculations using the LPm method.

for understanding the impact of geomagnetic disturbances on power networks is the ability to model the GIC produced in a network by specified geoelectric fields. In 1985, Lehtinen and Pirjola published a landmark paper that provides the first description of a stand-alone method for modelling GIC. The  method (hereafter referred to as the 'LP method') has been widely used in the geophysics and space weather community and provided the basis for GIC studies in many countries (e.g. Pirjola and Lehtinen, 1985;Mäkinen, 1993;Mäkinen et al., 1993;Thomson et al., 2005;Wik et al., 2008;35 Viljanen et al., 2012;Torta et al., 2014;2017;Divett et al., 2018).
The LP method was designed at a time when mostly only the highest voltage levels of a power network were considered in GIC calculations. This was because the transmission lines at the lower voltage levels have higher resistance so will experience smaller GIC values. However, in a desire to provide more comprehensive modelling of GIC in a power network, many modern studies are now looking to model GIC in multiple voltage levels in a power network. The LP method has been effectively 40 used for such studies (e.g., Mäkinen, 1993;Mäkinen et al., 1993;Viljanen et al., 2012;Divett et al., 2018); however, using the LP method for multiple voltage levels involves many ungrounded nodes, thus having infinite resistance to ground, which is numerically inconvenient. Also, the main focus of the LP method was the GIC flow to ground through the transformer primary windings, which was the desired output when modelling a single voltage level of a power network. However, models for multiple voltage levels require calculation of the nodal voltages which are then used to calculate the GIC in the transformer 45 windings .
In this paper we show how the LP method can easily be modified to efficiently model GIC in multiple voltage levels of a power network by converting the LP method to calculate the nodal voltages directly. First we summarise the steps in the LP method and then show how these can be modified to give the modified Lehtinen-Pirjola method (hereafter referred to as the 'LPm method'). We also show that the LPm method involves inversion of a matrix that is symmetric positive definite allowing 50 the use of efficient methods including sparse matrix techniques. Then we show how software for GIC calculations using the LP method can easily be converted to the LPm method and provide example calculations for the benchmark model introduced  The GIC modelling method derived by Lehtinen and Pirjola (1985), the 'LP' method, is produced by starting with Kirchhoff's current law that the net current flowing into a node, k, on branches from other nodes, n, is equal to the current flowing to ground from node k (LP equation 8): and relates the current in a branch ikn to the driving electromotive force (emf) ekn (if there is one), the voltage difference between the nodes at the ends of the line vk and vn, and the admittance of the branch (LP equation 7) ykn: Substituting (2) into (1) gives (LP equation 9): Note that, when considering multiple voltage levels, branches in the network consist of not just transmission lines but also transformer windings. The transmission lines experience the driving emf produced by the magnetic field variations, whereas the transformer windings do not.
The driving emf in each transmission line is represented by an equivalent current source = (4) 80 The equivalent current sources are then summed to give the current source directed into each node (LP equation 13).
Making this substitution in (3) gives 85 The first summation represents the dependence of current ik on voltage vk so gives diagonal elements of a network admittance 90 The second summation represents the dependence of current ik on all the other nodal voltages vn so gives the off-diagonal elements of the (symmetric) network admittance matrix 95 100 (note that the superscript n in the left hand side is not an index). Combining the above equations gives (LP equation 11): The values of [I e ] were the desired output when modelling a single voltage level of a power network. However, if there was more than one transformer at a substation (as usually occurs), it was necessary to split the current in proportion to the admittances of the transformer windings to determine the fraction of the current that flowed in each transformer winding. 140 Now, when modelling the GIC in multiple voltage levels of a power network, many of the nodes are ungrounded. However, the LP method needs to specify an earthing impedance for each node. This is done by adding 'virtual' connections to ground from each ungrounded node (Mäkinen, 1993;Pirjola, 2005 voltages are then used to calculate the GIC flow in the branches using equation (2). This is the equation to use for the GIC in the transmission lines. For branches of the network that are transformer windings, there is no driving emf so equation (2) reduces: 150 to give the GIC flow in the transformer windings.

Lehtinen-Pirjola Modified Method
When modelling GIC in multiple voltage levels of a power network, it is necessary to calculate the nodal voltages before 155 calculating the GIC in the transmission lines and transformer windings. In the Lehtinen-Pirjola modified (LPm) method the matrix equations are modified to provide a solution in terms of the nodal voltages. This also has the advantage that there is no need to add virtual earthing connections to ground from the ungrounded nodes.
To convert the currents flowing to ground [I e ] provided by the LP method to nodal voltages [V n ] start with equation (15) and 160 make the substitution from equation (12) (Pirjola, 2007) The LP method allows for the [Z e ] matrix to have off-diagonal elements representing the voltage produced at node i by currents 170 flowing to ground from other nodes (Pirjola, 2008). However, if the circuit is set up with a node at the neutral point of each substation, this does not happen (see . In this case, [Z e ] becomes diagonal with elements equal to the earthing resistances ri of the nodes and the inverse of [Z e ] is simply the earthing admittance matrix [Y e ] given by Then equation (19) can be rewritten as: The voltages of the nodes are then found by taking the inverse of the sum of the admittance matrices and multiplying by the 185 nodal current sources These node voltages can then be substituted into (2) and (17) to give the GIC in the transmission lines and the transformer 190 windings. The GIC flow to ground is simply given from Ohm's law using the neutral point node voltage and the admittance to ground (equation (18) The LPm method involves inversion of a matrix ([Y e ] + [Y n ]) which is symmetric (i.e., Hermitian as the elements are real) and positive definite and can thus be solved using a particularly efficient case of lower-upper (LU) decomposition, the Cholesky 195 decomposition (Press et al., 2007). Note that most of the nodes within the network are unconnected, meaning that [Y n ] has many zeros. This is also the case with [Y e ], so the Cholesky decomposition enables the use of sparse matrix methods (Stott and Alsaç, 1987;Press et al, 2007), thus providing an efficient way to model GIC in multiple voltage levels of a power network.

Calculation of GIC Time Series
GIC modelling is now being used, not just to assess the GIC for specified electric field values, but also to determine the 200 variation of GIC throughout a geomagnetic disturbance. If the network configuration does not change during that time (not always the case), then the matrix inversion does not need to be recalculated at every time step.
If the electric field is assumed to be uniform across the network then linear superposition can be used to calculate the GIC. (A uniform electric field would be produced, e.g., if calculations are made using data from a single magnetic observatory and a one-dimensional (1-D) earth conductivity model.) The GIC modelling can be made for two cases: i) a northward electric field 205 of 1 V/km, and ii) an eastward electric field of 1 V/km. For each location k in the network, this modelling gives reference GIC values, αk and βk for the northward and eastward electric fields that can be scaled by the actual electric field values at each time step and then combined to give the time series of GIC values at that location.
This concept can be extended for using two magnetic observatories , but this still requires use of a single 210 1-D earth conductivity model for the whole network.
In practice there is considerable variability in the earth conductivity structure across a power network. There are many modelling techniques for calculating the electric fields in such cases, ranging from use of multiple 1-D earth models (Marti et al., 2014) to use of magnetotelluric transfer functions and 3-D earth conductivity models (Weigel, 2017). In these cases, the 215 electric fields across the network can change from place to place and from one time step to the next. This will result in a  The advantage is that the solution of a triangular set of equations is quite trivial, as equation (27) can be solved by forward substitution and equation (28) can be solved by back substitution. Also, once the Cholesky decomposition has been done, it is possible to solve with as many right-hand sides as required, one at a time. Thus the LPm method provides a much faster and versatile way of calculating a time series of GIC values.

Comparison between the LP and LPm Methods
To illustrate the differences between the LP method and the LPm method consider the circuit for a substation with a twowinding transformer and three autotransformers as shown in Figure 1. The LP method requires the addition of virtual connections to ground from nodes 1 and 2, as explained above. However, in the LPm method the connection to ground is expressed as an admittance value. For the ungrounded nodes the admittance to ground is zero, which can easily be included 245 in the earthing admittance matrix without having to add virtual connections to the circuit. previous section). An advantage of the LPm method is that the matrix calculation yields the nodal voltages directly without the need to obtain them from the earthing currents as required in the LP method. Finally, the same step is performed in the LP and LPm methods to use the nodal voltages to calculate the GIC in the transmission lines and transformer windings.

6 Example Calculation using the LPm method
To illustrate the use of the LPm technique we present the calculation of GIC in the benchmark model of Horton et al. (2012) shown in Figure 3. The following tables will also provide values for testing when converting software from the LP method to the LPm method.
To construct the network admittance matrix [Y n ] and the earthing admittance matrix [Y e ] it is first necessary to assign node 265 numbers to the buses and neutral points at each substation. The modelling does not depend on any particular choice of node numbers. Here we use the node number assignment shown in Table 1. Note that some buses are not included in the model for calculating the GIC because they are connected to transformer windings in a delta configuration so there is no path for the GIC to flow. Also, there is no neutral point node at substation 7 because there are no transformers at this site as it is just a switching station connecting transmission lines. 270 The network admittance matrix [Y n ] is constructed using the transmission line information and transformer information presented in Tables II and III of Horton et al. (2012). In the network admittance matrix [Y n ] the values are given by equations 8 (for the diagonal elements) and 9 (for the off-diagonal elements). Note that the values are presented for a single phase of the power network and it is assumed that the other two phases are identical. The network admittance matrix values for the benchmark model are given in Table 2. 275 The earthing admittance matrix [Y e ] is constructed using the substation grounding resistance, rg, presented in Table I Table 3.
Note that the general theory is expressed in terms of impedance and admittance which can have reactive components, but, in practice, at the frequencies applicable to GIC the reactive components are negligible, and the network characteristics can be 285 described as purely resistive or conductive.
The resulting inversion of the matrix gives ([Y n ] + [Y e ]) -1 shown in Table 4. The voltage source in each transmission line and the equivalent current source, calculated from equation (4), are shown in Table 5. These current sources are then summed (equation 5) to give the nodal current sources [J e ] shown in Table 6.
The nodal current sources are then combined with the inverted matrix (Equation 22) to give the nodal voltages shown in Table   7. For nodes 1-11 these give the bus voltages shown in Table V of Horton et al. (2012). For nodes 12-18, combining the nodal voltage and the substation grounding resistance gives the GIC flow to ground for each substation shown in Table VII of Horton 295 et al., (2012).
The nodal voltages substituted into equations (2) and (17) give the GIC values for the transmission lines and transformer windings shown in Tables VI and VIII of Horton et al. (2012).

Discussion 300
The above example calculation shows that the LPm method provides GIC values that exactly match those for the benchmark model provided by Horton et al. (2012). The values in Tables 2 to 7 can also be used to check intermediate steps in the software used for the LPm calculations.
Any software developed to model GIC should be able to exactly match the values provided by the Horton et al. (2012) paper.
The results presented in that paper are not an average of modelling results nor an approximation to the correct values, but are 305 the identical values obtained using four different software implementations. However, initial calculations involving the four different software implementations provided similar but slightly different results. Further investigation showed that the origin of the differences was in the way that distances between substations were being calculated in the different implementations.
Some versions used formulas based on a spherical earth and some used formulas taking account of its non-spherical shape. It was then decided to standardise on substation latitudes and longitudes based on the WGS84 ellipsoid model of the Earth which 310 is used by the global navigation satellite system (GNSS) for geolocation. After this, all the calculations gave exactly the same results. To get the source voltage values presented in Table 5 (and hence match the GIC results for the benchmark model) thus requires using the formulas presented in the Appendix of Horton et al. (2012) for calculating distances between substations.
Many people have used the LP method for calculating GIC in the highest voltage level of their power networks. With the increasing requirement to calculate GIC in multiple voltage levels of a power network it is hoped that the new LPm method 315 described above will provide an easy way for converting existing LP software. The conversion is a simple process. and do the matrix inversion. This directly gives the nodal voltages which are required to calculate the GIC in the transmission lines and transformer windings. There is no need for any "virtual" nodes or connections. Also, it is more efficient as ([Y e ] + [Y n ]) is symmetric positive definite, so can be solved using Cholesky decomposition, which is a special case of lower-upper 320 (LU) decomposition (Press et al., 2007) with U = L T (i.e., the upper triangular factor is the transpose of the lower triangular factor).
Power networks, on average, have three transmission lines and one or two transformer windings connected to a bus, so a typical row in the admittance matrix has only five or six non-zero elements, independent of the overall network size. Thus for larger networks, where node numbers can be in the thousands, the admittance matrix will have over 99 % of its values equal to zero. 325 Cholesky factorization takes advantage of this fact by making use of sparse matrix methods (Stott and Alsaç, 1987;Press et al, 2007), thus additionally reducing memory usage and computation time. To examine how this affects the GIC modelling, we performed calculations for two networks using both the LP and the LPm methods. The networks modelled were 1) the benchmark network of Horton et al. (2012), which has 18 nodes, and 2) the nation-wide Spanish Power Grid operated by Red Eléctrica de España (REE), which has 1388 nodes. GIC in the 400 kV part of the REE system was considered by Torta et al. 330 (2014); for this study we include both the 400 and 220 kV levels of the REE network (see Torta et al., 2021 for reference).
Tests we did showed that the LP and LPm methods both produce matrices that are sparse, so there is potential for sparse matrix techniques to be applicable. Table 8 shows the calculation times and memory usage for GIC calculations using the LP and LPm methods. These show that memory usage was drastically reduced when using sparse matrix techniques, with the reduction being more significant with the larger REE network. The time for the matrix inversion is significantly affected, as expected, 335 by the size of the matrix involved. For the Horton network (18x18 matrix) the change to sparse techniques actually increased the inversion time. However, the sparse techniques applied to the REE network (1388x1388 matrix) produced an approximately order of magnitude reduction in inversion time. This reduction in inversion time was greatest for the LPm method which was nearly an order of magnitude faster than the LP method.
The column 'Inversion time' in Table 8 reflects the time required to compute [ ] from equations (25), (27) and (28), thus 340 including the Cholesky decomposition and the forward and back substitutions in the LPm method; note, in consequence, that it is not strictly an inversion, though we will refer to it as such. Also note that, when referred to LP, this column reflects the  Table 8 are obtained from GIC modelling using programs in Matlab. Specific inversion times and memory usage will vary with the programming language used, but it is expected that the general results presented here will apply regardless of the programming language used. Guidance is provided for converting software from the LP method to the LPm method and an example calculation using the benchmark model of Horton et al. (2012) is presented to provide a set of values for testing GIC calculation software.
Calculations of GIC using the LPm method involve a matrix that is symmetric positive definite. This enables a solution to be 360 obtained by Cholesky decomposition, (a specific case of LU decomposition) which is numerically more accurate than computing the matrix inversion itself. The factorization of Cholesky decomposition can always be implemented using sparse matrix techniques, speeding up the calculations for large networks.     Figure 3. 485 Table 3. Earthing Admittance matrix [Y e ] for the benchmark model shown in Figure 3.  Figure 3. 490