Continental scale modelling of geomagnetically induced currents

The EURISGIC project (European Risk from Geomagnetically Induced Currents) aims at deriving statistics of geomagnetically induced currents (GIC) in the European high-voltage power grids. Such a continent-wide system of more than 1500 substations and transmission lines requires updates of the previous modelling, which has dealt with national grids in fairly small geographic areas. We present here how GIC modelling can be conveniently performed on a spherical surface with minor changes in the previous technique. We derive the exact formulation to calculate geovoltages on the surface of a sphere and show its practical approximation in a fast vectorised form. Using the model of the old Finnish power grid and a much larger prototype model of European high-voltage power grids, we validate the new technique by comparing it to the old one. We also compare model results to measured data in the following cases: geoelectric field at the Nagycenk observatory, Hungary; GIC at a Russian transformer; GIC along the Finnish natural gas pipeline. In all cases, the new method works reasonably well.


Introduction
The EURISGIC project (European Risk from Geomagnetically Induced Currents) will produce the first European-wide realtime prototype forecast service of geomagnetically induced currents (GIC) in power grids, based on in-situ solar wind observations and comprehensive simulations of the Earth's magnetosphere Viljanen (2011).By utilising geomagnetic recordings, we will also derive a map of the statistical occurrence of GIC throughout Europe.To reach these goals, we have extended previous methods to be applicable to continental scales.
There are a large number of GIC studies performed in Europe, for example: Finland: Lehtinen & Pirjola (1985) and Pirjola & Lehtinen (1985) have developed the basic methods of GIC modelling in power grids, and Pulkkinen et al. (2001b) generalised the approach by Boteler (1997) for pipelines.Several studies have been made for the Finnish high-voltage power grid (Pirjola 1989;Elovaara et al. 1992;Ma ¨kinen 1993) and the natural gas pipeline (Pulkkinen et al. 2001a;Viljanen et al. 2006).Sweden: Wik et al. (2008) performed a statistical study of GIC in South Sweden, and Wik et al. (2009) and Pulkkinen et al. (2005) considered the most extreme GIC events observed in Sweden in July 1982and October 2003. United Kingdom: Thomson et al. (2005) modelled GIC during the extreme storm of October 2003.Pulkkinen et al. (2003b) performed a comparative study of a large event with simultaneous observations from UK and Finland.
In this paper, we will first summarise the general methodology and then describe in more details the specific new features: the model of European high-voltage power grids, the model of the ground conductivity and calculation of GIC in a large region in spherical geometry.We verify the results by comparisons to previous studies on smaller geographic regions and also show some examples of comparisons between models and measurements.

Basic method
Calculation of GIC in a given conductor system is convenient to divide into two steps (Pirjola 2002): 1. Determine the geoelectric field associated with geomagnetic variations.This is a purely geophysical problem, which is independent of the technological system.
2. Determine GIC due to the given geoelectric field in a conductor system whose topology and resistances are known.
In this paper, we will only consider power grids, except for one example of a pipeline event.For the modelling technique of GIC in pipelines, see Pulkkinen et al. (2001b) and Trichtchenko & Boteler (2001), for example.

Geoelectric field
We apply the method of spherical elementary current systems (SECS) to calculating the geoelectric field as outlined by Viljanen et al. (2004).The first task is to determine ionospheric equivalent currents in terms of elementary systems.This approach was introduced by Amm (1997) and Amm & Viljanen (1999) and extensively demonstrated by Pulkkinen et al. (2003a).Tests for sparse magnetometer networks by McLay & Beggan (2010) further confirmed the usefulness of the method.For GIC calculations, we do not explicitly need ionospheric currents, but only the geomagnetic variation field interpolated in the area of the power grid.This interpolation is easy to perform with the SECS method.
GIC depends on the electric field integrated along power lines.Because integration is a spatially smoothing operation, we do not need to know the electric field at a kilometre scale, but a resolution of about 100 km is acceptable.We can thus reasonably assume locally a layered ground conductivity structure.Then the horizontal electric field (E x , E y ) is related to the magnetic field (B x , B y ) as where Z is the plane wave surface impedance, x is the angular frequency and (x, y) refer to the north and east components.This approach, called the local plane wave method, is validated by Viljanen et al. (2004).

GIC in a power grid
A power grid is a discretely grounded system with N nodes (earthing points).Lehtinen & Pirjola (1985) present the following formula for the N • 1 column matrix I e consisting of the earthing currents at the nodes: where U is a unit matrix and Y n and Z e are the network admittance matrix and earthing impedance matrix, respectively.GIC vary slowly compared to the 50 Hz AC, so a DC treatment is sufficient, and Y n and Z e are real.The earthing impedance matrix relates the earthing currents with voltages between the nodes and a remote earth.If the nodes are far from each other then Z e is diagonal with the elements equal to the earthing resistances (Z e ii ¼ R e i ).The network admittance matrix Y n is defined by where R n ij denotes the conductor line resistance between nodes i and j.The elements of the N • 1 column matrix J e are defined by: The conductor line from node j to i is denoted by s ji , and E is the geoelectric field.Thus V ji is the voltage affecting the line s ji .
The transmission line current from node i to j can be solved from the equation We follow the convention that the three-phase lines are treated as a single circuit element.For the special cases of autotransformers and full-wound transformers, see Ma ¨kinen (1993) or Pirjola (2005).(2008a, 2008b) indicate that a reasonable understanding of GIC can be obtained even with incomplete data.
The following assumptions and simplifications are used: -Not all substations or lines are included.
-The locations of the stations are given with a varying accuracy of up to 10-20 km (0.1-0.2 degrees in latitude and longitude).-All transmission lines are straight and follow the great circle between the end points.If this were not the case then the line could be constructed of shorter straight sections.
In reality, the lines are slightly longer than the great circle distance.This means that their resistances are larger than assumed here, so we obtain somewhat too high GIC estimates.-All stations are far enough from each other to ignore the off-diagonal elements of the earthing impedance matrix Z e .-We do not have information about the number of transformers per substation.We simply assume that all diagonal elements of Z e are equal to 0.5 X.
-The value used for the line resistance per unit length is 0.008 X/km for 400 kV lines (Elovaara et al. 1992;Elovaara 2007;Pirjola 2008a) and 0.022 X/km for 220 kV lines (typical in Finland).For 300-330 kV lines (Baltic countries and Norway), we assume an interpolated value of 0.015 X/km and for 275 kV lines (UK) 0.0185 X/km.These values take into account the reduction of three phases into a single conductor.
Figure 1 shows the prototype model used in this study based on various sources such as ENTSO-E maps, individual maps of smaller areas and previous models used in local studies.The model consists of five separate main parts: (1) South and Central Europe; (2) Nordic countries (Finland, Sweden, Norway, East Denmark); (3) United Kingdom; (4) Baltic countries (Estonia, Latvia, Lithuania); (5) North-West Russia.Although there are high-voltage DC links between these areas, they do not have a galvanic connection for the AC current.Consequently, GIC in one part cannot directly flow into another part.

Ground conductivity model
The ground conductivity model EURHOM (EUropean RHO Model) is based on extensive literature study and inversion of published magnetotelluric sounding curves into one-dimensional layered structures.Details of the construction of the model will be a topic of a separate paper, so we only outline the main results here.We note that some modifications of the conductivity map are foreseen whenever new data from geomagnetic and magnetotelluric soundings and from GIC recordings become available.
Figure 2 shows the locations of rectangular blocks comprising the conductivity model.The different numbers refer to different layered structures defined by layer thicknesses and resistivities (inverse of conductivities), and these parameters are explicitly listed on http://real.mtak.hu/2957/.We assume the vacuum permeability l 0 everywhere and a constant permittivity e = 5e 0 , where e 0 is the vacuum permittivity.The exact value of the permittivity has no practical role due to vanishingly small displacement currents compared to Ohmic currents.
Using the resistivities and layer thicknesses, we calculated the surface impedance Z in Eq. ( 1) for each different block of Figure 2 using a fixed length of an event of 2048 timesteps (1-min values).This allows for applying the fast Fourier transform in transformations between the time and frequency domains.The surface impedance is determined as if the block were an infinite half-space with a flat surface.This is clearly acceptable taking into account the sizes of blocks in Figure 2.
1. Local AE (auroral electrojet) index is defined by the difference of the upper and lower envelope curves of the northward (X) component of the magnetic field: where {X(t)} refers to the value of the field at all sites at time t.
2. Largest horizontal difference field: where dX and dY are the differences between two successive values of the northward and eastward components.Magnetometer data were used to determine ionospheric equivalent currents with the SECS method.More precisely, we calculated the amplitudes of the divergence-free elementary systems (Eq.( 6) of Amm & Viljanen 1999).The poles of the elementary systems were located as follows: a denser grid above North Europe and a sparser one elsewhere as shown in Figure 4. We also tested if the configuration of SECS poles affects calculated GIC values.The conclusion is that the SECS grid should cover the region of magnetometers with some outward extension.
Given the poles of SECS and their amplitudes, we can calculate the magnetic field at any location at the earth's surface (Eqs.( 9) and (10) of Amm & Viljanen 1999).We need the horizontal field only at the nodes of the power grid.

GIC in spherical geometry
The general expression for GIC (Eq.( 2)) is valid both in planar and spherical geometries.It is more natural to use spherical geometry, since it does not require any planar projection that might distort distances between power grid nodes.Geometry appears explicitly when the voltage is calculated from Eq. ( 6).As a simplification mentioned in Sect.2.2.1, we assume that the transmission line between two nodes follows the   corresponding great circle arc.It is then possible to evaluate the integral in Eq. ( 6), but the exact result is not practical.The distance between the end points of a transmission line is at most a few hundred kilometres (cf.Fig. 1).Assuming that the electric field does not vary much between these points, we can approximate the field by a constant value along the line.Explicit results and a matrix formulation are derived in Appendix A.

Differences between spherical and planar methods
If the electric field E ¼ E x e x þ E y e y is spatially uniform then GIC at a given site is where (a, b) depend only on the power grid parameters (Pirjola & Lehtinen 1985).We compared the method of the planar geometry (Sect.A.1.4) with the spherical method by giving a constant value (1 V/km north or east) for the electric field at all nodes of the power grid.An independent older code was used to derive directly the coefficients a and b yielding the same results.
Concerning spherical geometry, we can similarly assume a constant electric field.However, there is a principal difference to the planar case when considering an eastward field, since such a vector field parallel to latitudinal lines is not curl-free (Sect.A.1.1).In a small region this is not a problem in practice.As the basic model, we considered the Finnish 400 kV power grid in 1978-1979 (see Appendix A).Some differences in GIC calculated in the two geometries appear when a spatially uniform electric field is assumed.In general, both geometries yield very similar results (Table 1).
In another test with the Finnish grid in 1978-1979, we used the ground conductivity map of Figure 2. The geomagnetic event considered occurred on 30 October 2003, and we use 1-min data of all available IMAGE magnetometer stations.In the planar geometry, we centred the stereographic projection at (61.5 N, 24.9 E).The difference between GIC in spherical and planar methods is small as shown in Figure 5 for the Huutokoski station.
Next, we considered the South and Central European part of the grid, in which there are about 1100 nodes and 1600 transmission lines.We again used the conductivity map of Figure 2 and geomagnetic data of 29 October 2003 from all available European sites (Fig. 3).In the planar geometry, we centred the stereographic projection at 46.1 N, 10.0 E. Comparison of the results between spherical and planar methods is shown in Figure 6 for the sum of absolute GIC values at all nodes.The relative difference between the two geometries is again   small.In the following examples, we use the spherical method only.

Measured and modelled electric fields at Nagycenk
Recordings of the geoelectric field are available from the Nagycenk Observatory (NCK, 47.6 N, 16.7 E) in Hungary since 1957 and in digital form since 1996 (Kis et al. 2007).
The ground model of this region consists of three layers (thicknesses 3 km and 57 km above an infinitely deep layer), and their resistivities are 5, 1000 and 10 Xm (block 3 in Fig. 2).
To demonstrate the robustness of the SECS method, we excluded NCK and the nearby Tihany (THY, 46.9 N, 17.5 E) and Hurbanovo (HRB, 47.9 N, 18.2 E) observatories in the input magnetic field data.We interpolated the magnetic field at the location of NCK and calculated the electric field there using the local surface impedance.The result is nearly perfect (Fig. 7).
We also performed a more radical test by excluding all magnetometer stations south of the latitude 55 degree.Even then we obtained a reasonable result (figure not shown).The opposite test with including only stations south of the latitude 55 degree gave an excellent result.It is evident that magnetic field variations are spatially quite smooth in South and Central Europe, so a sparse magnetometer network can still provide a sufficient accuracy in interpolation and even in extrapolation (cf.McLay & Beggan 2010).

GIC in the North-West Russian power grid
In 2011, five GIC recording sites were (re)installed in the 330 kV power grid in North-West Russia.An initial DC description of the grid has also been prepared with more precise parameters than in the European prototype model.Figure 8 shows an example of a large event at the substation Vykhodnoy (VKH) (68.3 N, 33.1 E).In this case, we used all available magnetometer data from 10 Finnish sites, one Estonian site, one Norwegian site and Lovozero (67.97 N, 35.08 E), which is the closest observatory to VKH.Then we interpolated the magnetic field at VKH and calculated GIC from the electric field applying Eq. ( 10) (this simplification is a good approximation as shown by Viljanen et al. 2004).We tested two different conductivity models.In the first case, we applied a two-layer model, where the upper layer is 30 km with a resistivity of 5000 Xm and the lower infinite layer has a resistivity of 200 Xm.The second case, used in the figure, has also two layers with parameters 300 km, and 3000 and 3 Xm.
The modelled GIC corresponds qualitatively quite well to the measured one.However, there are some obvious deficiencies: power grid parameters may be inaccurate, the ground model needs adjustments and some noise seems to appear in GIC (especially the spike around 02:30 UT looks erroneous).A more systematic analysis of the North-West Russian grid and GIC data will be a topic of a separate paper after a larger set of recordings will be available.

GIC in the Finnish natural gas pipeline
In the example shown in Figure 9, we used data from nine Finnish, one Norwegian and one Estonian IMAGE magnetometer station, but excluded the Nurmija ¨rvi observatory, which is close to the pipeline GIC recording site at Ma ¨ntsa ¨la ¨(60.6 N, 25.2 E).We calculated GIC using Eq. ( 10) with parameters by Pulkkinen et al. (2001a).We also performed the same calculation with data from Nurmija ¨rvi included, and obtained a nearly perfect fit especially just after 17 UT (figure not shown).
In this case, we used the following ground model (block 25 in Fig. 2): the upper layer is 30 km with a resistivity of 5000 Xm and the lower infinitely deep layer has a resistivity of 20 Xm.Although this is different from the model applied previously by Viljanen et al. (2006), it gives very accurate results.It is evident that a large variety of apparently very different ground models can yield a satisfactory result (e.g., Parker 1983).This is a topic to be studied in more detail later when assessing the sensitivity of GIC statistics to the ground model.

Solar cycle test: Finnish grid in 1978-1979
To test the new codes with a long data series, we calculated GIC in the Finnish power grid in 1978-1979 using 1-min IMAGE magnetometer data of the years 1996-2006 (4018 days).We applied the ground model of Figure 2. The long-term GIC activity is illustrated in Figure 10.The quantity plotted is based on the sum of absolute GIC through all nodes of the power grid.Its daily maxima have been used to calculate the monthly indicator.GIC activity has no prominent correlation with the sunspot number.It can also be at comparable levels both during sunspot minimum and maximum, and large GIC values may occur at any time.The most active GIC month was October 2003 during the descending phase of the solar cycle 23.The most active day of the cycle was 30 October, when the Swedish city Malmo ¨experienced a blackout due to GIC (Pulkkinen et al. 2005).

Computational performance
We tested the efficiency of the new codes with the following parameters: -Power grid of 1107 nodes and 1557 transmission lines (South and Central Europe).-One-minute magnetic field data of 2003 from about 40 European magnetometer sites.-Twenty-two different ground conductivity blocks.
-GIC were calculated at all nodes as 1-min values for the whole year 2003.
GNU Octave was used (the code runs in MatLab too).Computation time in a modern laptop was about 30 min.Consequently, modelling a full sunspot cycle takes only a few hours which makes it feasible to perform extensive test runs with different model parameters.

Conclusions and future work
The main parts of the updated GIC calculation are 1.General input: block model of the ground conductivity, continental power grid model, magnetic field data from tens of sites.2. Calculation of the ionospheric equivalent current system from measured geomagnetic variations using the method of spherical elementary current systems.3. Interpolation of the magnetic field to the nodes of the power grid by utilising the SECS method.4. Calculation of the geoelectric field at the nodes using the local plane wave method.5. Calculation of voltages between nodes on a spherical surface.6. Calculation of GIC.
Comparisons with the results by the old code as well as comparisons to measurements validate the success of the updated method.The computational performance allows for using geomagnetic data of a solar cycle.
The next extensive application of the new code is to produce statistics of the occurrence of GIC in the European high-voltage power grid.Parallel to this investigation, we will establish a prototype GIC forecast server for the European power grid.It will use magnetic field values obtained by magnetohydrodynamic simulations starting from real-time solar wind observations.

Appendix A: Calculation of GIC in spherical geometry
We present first the general expression of the voltage due to the geoelectric field on a spherical surface.A practical approximation is also derived, and the corresponding result in planar geometry is given for reference.
A.1.Voltage along a great circle arc Assume that the horizontal geoelectric field E is known at the earth's surface and expressed in spherical components: ðA:1Þ Consider a transmission line between nodes 1 and 2. The voltage V 12 between these points is the line integral of the electric field: Generally, it is not possible to compute this integral analytically.
If we assume that the transmission line follows the path expressed by h = h(/) then the voltage takes the form ðA:3Þ where R is the Earth's radius.Assuming now the special case of a great circle arc, the dependence of the polar angle h on the azimuthal angle / can be expressed implicitly as follows: where the constants A, B, C depend on the (h, /) coordinates of the points 1 and 2: ; ðA:5Þ The integral in Eq. (A.3) is difficult to evaluate for a general electric field.The only trivial example is a uniform northward electric field E n = ÀE 0 e h (E 0 > 0, e h points to the south).In this case the line integral in Eq. (A.2) is path-independent, and the voltage between two points depends only on their latitudinal difference: V 12 = ÀE 0 R(h 2 À h 1 ).The apparently simple uniform eastward field E e = E 0 e / has a non-zero vertical component of the curl, so the voltage is path-dependent.

A.1.2. Approximate method
A realistic assumption is that the distance between nodes 1 and 2 is at most a few hundred kilometres.We assume that the elec-tric field does not vary much in such distance scales, so we can approximate the field by a constant value along the transmission line from 1 to 2. A natural choice for this constant is the mean value of the field components parallel to the great circle at the end points 1 and 2. The voltage is then equal to this mean field multiplied by the arc length.Consider two points (r 1 , r 2 ) at the Earth's surface with polar coordinates (R, h 1 , / 1 ) and (R, h 2 , / 2 ) (Fig. A.1).We denote their difference vector by s ¼ r 2 À r 1 : ðA:6Þ The task is to determine the component of E pointing to the direction given by s, but being tangential to the Earth's surface.The unit vector u ¼ s Â r 1 js Â r 1 j ðA:7Þ is perpendicular both to r 1 and s, and tangential to the surface at r 1 = (R, h 1 , / 1 ).It follows that the unit vector is tangential to the surface and points to the direction of the great circle arc.The component of the electric field parallel to the great circle arc is then Similarly, we obtain at the point 2 The angle c between the radius vectors r 1 , r 2 is defined by cos c ¼ cos h 1 cos h 2 þ sin h 1 sin h 2 cos / 2 À / 1 ð Þ : ðA:11Þ In all practical situations, 0 c p, so sin c ! 0, and |r 2 • r 1 | = R 2 sin c.Then the average electric field E par = (E par1 + E par2 )/2 along the great circle arc between points 1 and 2 is:  Table A.1.Parameters of the Finnish power grid in 1978-1979(Pirjola & Lehtinen 1985).Upper part: nodes, lower part: transmission lines.
Locations of the stations are approximate.The earthing and line resistances R e are given in ohms.Pirjola & Lehtinen (1985) and Pirjola (2009) assumed a zero R e at the Swedish border stations Letsi and Messaure.We use the value of 0.50 X to roughly approximate the effect of the rest of the Swedish power grid.Note that there are presently many more stations in Finland, and their parameters differ from the ones given here.

Fig. 1 .
Fig. 1.Prototype model of high-voltage power grids in Europe.

Fig. 4 .
Fig. 4. Locations of the poles of spherical elementary current systems used for calculating GIC in Europe.

Fig. 5 .
Fig. 5. Upper panel: calculated GIC at the Huutokoski 400 kV station using data of 30 October 2003, assuming the Finnish grid model in 1978-1979, the conductivity map of Figure 2, and using spherical geometry.Lower panel: difference in GIC by spherical and planar methods.Note the different scales on the vertical axes.A positive GIC flows into the earth.

Fig. 6 .
Fig. 6.Upper panel: sum of the calculated GIC at all stations on 29 October 2003, using the South and Central European part of the grid model of Figure 1, the conductivity map of Figure 2, and the spherical method.Lower panel: difference of the GIC sum by spherical and planar methods.

Fig. 7 .
Fig. 7. Measured (black) and modelled (blue) electric fields at the location of Nagycenk (NCK) when all available European magnetometer data (except NCK and its nearest neighbours Tihany and Hurbanovo) are used as input to interpolate the magnetic field at NCK. Upper panel: northward component, lower panel: eastward component.A three-layer conductivity model was applied (block 3 in Fig. 2).

Fig. 8 .
Fig. 8. Measured (black) and modelled (blue) GIC at the Russian substation Vykhodnoy on 15 March 2012.A two-layer conductivity model is used as explained in the text.

Fig. 10 .
Fig. 10.Bars: monthly means of the daily maximum GIC sum using the Finnish grid model of 1978-1979.Red curve: monthly sunspot number.

Fig. 9 .
Fig. 9. GIC along the Finnish natural gas pipeline at the Ma ¨ntsa ¨la compressor station on 15 March 2012.Black curve: measured value, blue: modelled (without data from the nearby Nurmija ¨rvi magnetic observatory).One-minute averages are shown.A two-layer conductivity model is used (block 25 in Fig. 2).
Fig. A.1.Vectors related to two nodes on a spherical surface.

Fig. A. 2 .
Fig. A.2. Schematic map of the Finnish 400 kV power grid in 1978-1979 with transmission lines approximated by straight lines.

Table 1 .
GIC (amperes) due to a spatially constant northward or eastward electric field (1 V/km) in the Finnish high-voltage power grid in 1978-1979.''Plane'' refers to the planar method and ''sphere'' to the spherical method.
Viljanen et al.: Continental scale modelling of geomagnetically induced currents