Issue 
J. Space Weather Space Clim.
Volume 2, 2012
EUFP7 funded space weather projects



Article Number  A17  
Number of page(s)  11  
DOI  https://doi.org/10.1051/swsc/2012017  
Published online  28 September 2012 
Continental scale modelling of geomagnetically induced currents
^{1}
Finnish Meteorological Institute, Helsinki, Finland
^{2}
Natural Resources Canada, Ottawa, Canada
^{3}
NeuroSpace, Malmö, Sweden
^{4}
Research Centre for Astronomy and Earth Sciences of the Hungarian Academy of Sciences, Sopron, Hungary
^{5}
Polar Geophysical Institute, Apatity, Russia
^{*} corresponding author: email: ari.viljanen@fmi.fi
Received:
23
May
2012
Accepted:
17
September
2012
The EURISGIC project (European Risk from Geomagnetically Induced Currents) aims at deriving statistics of geomagnetically induced currents (GIC) in the European highvoltage power grids. Such a continentwide 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 highvoltage 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.
Key words: 7904: geomagnetically induced currents (GIC)
© Owned by the authors, Published by EDP Sciences 2012
This is an Open Access article distributed under the terms of creative Commons Attribution 3.0 Unported License
1 Introduction
The EURISGIC project (European Risk from Geomagnetically Induced Currents) will produce the first Europeanwide realtime prototype forecast service of geomagnetically induced currents (GIC) in power grids, based on insitu 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 highvoltage power grid (Pirjola 1989; Elovaara et al. 1992; Mä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 1982 and 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.

Elsewhere in the world, the largest interest in GIC has been in North America as demonstrated, for example, by Albertson & Thorson (1974), Lanzerotti (1979), Bolduc (2002), Trichtchenko & Boteler (2004), Kappenman (2005) and Pulkkinen et al. (2010). In 2000s, and especially motivated by the October 2003 event, GIC has reached attention also at lower geomagnetic latitudes: see Bernhardi et al. (2008) (South Africa), Trivedi et al. (2007) (Brazil), Liu et al. (2009) (China), Watari et al. (2009) (Japan) and Marshall et al. (2011) (Australia).
In this paper, we will first summarise the general methodology and then describe in more details the specific new features: the model of European highvoltage 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.
2 Calculation of GIC
2.1 Basic method
Calculation of GIC in a given conductor system is convenient to divide into two steps (Pirjola 2002):

Determine the geoelectric field associated with geomagnetic variations. This is a purely geophysical problem, which is independent of the technological system.

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.
2.1.1 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(1)where Z is the plane wave surface impedance, ω 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).
2.1.2 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:(2)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 (). The network admittance matrix Y^{n} is defined by(3) where denotes the conductor line resistance between nodes i and j. The elements of the N × 1 column matrix J^{e} are defined by:(4) (5) (6)
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(7)
We follow the convention that the threephase lines are treated as a single circuit element. For the special cases of autotransformers and fullwound transformers, see Mäkinen (1993) or Pirjola (2005).
2.2 Details of continental scale modelling
2.2.1 Power grid model
The European power grids at the voltage level of 220 kV and above contain a few thousand substations and transmission lines with frequent temporary and permanent changes. Detailed information of the grid is not publicly available, so simplifying assumptions are necessary. One purpose of the EURISGIC project is to demonstrate that GIC modelling in a large area is feasible. To reach this goal, we consider a “prototype” model, which covers Europe and has a realistically large number of nodes and lines, and qualitatively mimics the true situation. If more accurate information becomes available then it is straightforward to update this template model. The examples by Pirjola (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 offdiagonal 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 Ω.

The value used for the line resistance per unit length is 0.008 Ω/km for 400 kV lines (Elovaara et al. 1992; Elovaara 2007; Pirjola 2008a) and 0.022 Ω/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 Ω/km and for 275 kV lines (UK) 0.0185 Ω/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 ENTSOE 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) NorthWest Russia. Although there are highvoltage 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.
Fig. 1 Prototype model of highvoltage power grids in Europe. 
2.2.2 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 onedimensional 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 μ_{0} everywhere and a constant permittivity ε = 5ε_{0}, where ε_{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.
Fig. 2 Map of conductivity blocks in Europe. Each number denotes a different 1D ground model, whose parameters (layer thicknesses and resistivities) are listed on http://real.mtak.hu/2957/. Blocks outside of the presently modelled power grid are marked by 0. 
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 (1min 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 halfspace with a flat surface. This is clearly acceptable taking into account the sizes of blocks in Figure 2.
2.2.3 Interpolation of the magnetic field
We used 1min averaged European magnetometer data in 1996–2008 from stations shown in Figure 3. The number of available sites increases from about 35 in 1996 to about 50 in 2008. We determined quiettime baselines by subtracting the average of a 2day sequence starting from the noon and lasting 48 h. The actual event under study is the full day in the middle of the 48h sequence. This is the simplest objective method, which is appropriate when a large number of stations are considered. Because we apply the Fourier transform to the interpolated field (cf. Eq. (1)), a possible small shift from the “true” quiettime level is insignificant.
Fig. 3 Magnetometer stations used in this study. Blue dots: IMAGE magnetometer network (http://space.fmi.fi/image/), black dots: sites available at the World Data Centre for Geomagnetism (Edinburgh) (http://www.wdc.bgs.ac.uk/) (also belonging to INTERMAGNET). 
The data have a high quality, so large error spikes are rare. We calculated two quicklook indicators to assess the quality:
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:(8)where {X(t)} refers to the value of the field at all sites at time t.
2. Largest horizontal difference field:(9) 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 divergencefree 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.
Fig. 4 Locations of the poles of spherical elementary current systems used for calculating GIC in Europe. 
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.
2.2.4 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.
2.3 Test results
2.3.1 Differences between spherical and planar methods
If the electric field is spatially uniform then GIC at a given site is(10)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 curlfree (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).
GIC (amperes) due to a spatially constant northward or eastward electric field (1 V/km) in the Finnish highvoltage power grid in 1978–1979. “Plane” refers to the planar method and “sphere” to the spherical method.
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 1min 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.
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. 
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.
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. 
2.3.2 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 Ωm (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).
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 threelayer conductivity model was applied (block 3 in Fig. 2). 
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).
2.3.3 GIC in the NorthWest Russian power grid
In 2011, five GIC recording sites were (re)installed in the 330 kV power grid in NorthWest 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 twolayer model, where the upper layer is 30 km with a resistivity of 5000 Ωm and the lower infinite layer has a resistivity of 200 Ωm. The second case, used in the figure, has also two layers with parameters 300 km, and 3000 and 3 Ωm.
Fig. 8 Measured (black) and modelled (blue) GIC at the Russian substation Vykhodnoy on 15 March 2012. A twolayer conductivity model is used as explained in the text. 
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 NorthWest Russian grid and GIC data will be a topic of a separate paper after a larger set of recordings will be available.
2.3.4 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 Nurmijärvi observatory, which is close to the pipeline GIC recording site at Mäntsälä (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 Nurmijärvi included, and obtained a nearly perfect fit especially just after 17 UT (figure not shown).
Fig. 9 GIC along the Finnish natural gas pipeline at the Mäntsälä compressor station on 15 March 2012. Black curve: measured value, blue: modelled (without data from the nearby Nurmijärvi magnetic observatory). Oneminute averages are shown. A twolayer conductivity model is used (block 25 in Fig. 2). 
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 Ωm and the lower infinitely deep layer has a resistivity of 20 Ωm. 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.
2.3.5 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 1min IMAGE magnetometer data of the years 1996–2006 (4018 days). We applied the ground model of Figure 2. The longterm 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 Malmö experienced a blackout due to GIC (Pulkkinen et al. 2005).
Fig. 10 Bars: monthly means of the daily maximum GIC sum using the Finnish grid model of 1978–1979. Red curve: monthly sunspot number. 
2.3.6 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).

Oneminute magnetic field data of 2003 from about 40 European magnetometer sites.

Twentytwo different ground conductivity blocks.

Ionospheric grid of 283 SECS poles.

GIC were calculated at all nodes as 1min 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.
3 Conclusions and future work
The main parts of the updated GIC calculation are

General input: block model of the ground conductivity, continental power grid model, magnetic field data from tens of sites.

Calculation of the ionospheric equivalent current system from measured geomagnetic variations using the method of spherical elementary current systems.

Interpolation of the magnetic field to the nodes of the power grid by utilising the SECS method.

Calculation of the geoelectric field at the nodes using the local plane wave method.

Calculation of voltages between nodes on a spherical surface.

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 highvoltage 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 realtime solar wind observations.
Acknowledgments
The research leading to these results has received funding from the European Community’s Seventh Framework Programme (FP7/20072013) under Grant Agreement No. 260330. We thank the large number of institutes providing magnetometer data to the World Data Centre for Geomagnetism (Edinburgh) and INTERMAGNET, and to the IMAGE magnetometer network. Lasse Häkkinen (Finnish Meteorological Institute) prepared tools to combine these data sets. We are also grateful to Fingrid Oyj and Gasum Oy for a longterm interest and support to GIC studies in Finland. Software related to the calculation of equivalent currents and interpolation of the magnetic field was improved in collaboration with the ECLAT project (European Cluster Assimilation Technology), which has received funding from the European Community’s Seventh Framework Programme (FP7/20072013) under Grant Agreement No. 263325. Special thanks go to Drs. Liisa Juusola and Heikki Vanhamäki (Finnish Meteorological Institute) for developing the code.
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
A.1.1 General expression
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:(A.2)
Generally, it is not possible to compute this integral analytically. If we assume that the transmission line follows the path expressed by θ = θ(ϕ) 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 θ on the azimuthal angle ϕ can be expressed implicitly as follows:(A.4)where the constants A, B, C depend on the (θ, ϕ) 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_{θ} (E_{0} > 0, e_{θ} points to the south). In this case the line integral in Eq. (A.2) is pathindependent, and the voltage between two points depends only on their latitudinal difference: V_{12} = −E_{0}R(θ_{2} − θ_{1}). The apparently simple uniform eastward field E_{e} = E_{0}e_{ϕ} has a nonzero vertical component of the curl, so the voltage is pathdependent. Instead, an eastward field E_{e} = E_{0}/sin θe_{ϕ} is associated with an pathindependent voltage V_{12} = −E_{0}R(ϕ_{2} − ϕ_{1}).
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 electric 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, θ_{1}, ϕ_{1}) and (R, θ_{2}, ϕ_{2}) (Fig. A.1). We denote their difference vector by(A.6)
Fig. A.1 Vectors related to two nodes on a spherical surface. 
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(A.7)is perpendicular both to r_{1} and s, and tangential to the surface at r_{1} = (R, θ_{1}, ϕ_{1}). It follows that the unit vector(A.8)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(A.9)
Similarly, we obtain at the point 2(A.10)
The angle γ between the radius vectors r_{1}, r_{2} is defined by(A.11)
In all practical situations, 0 ≤ γ ≤ π, so sin γ ≥ 0, and r_{2} × r_{1} = R^{2}sin γ. Then the average electric field E_{par} = (E_{par1} + E_{par2})/2 along the great circle arc between points 1 and 2 is:(A.12)
A.1.3 Vectorised form
The geometrical factors containing trigonometric functions need to be computed only once for a given power grid. Then the final formulas for calculating GIC reduce to efficient matrix operations, for example, in MatLab or Octave.
The earthing currents at N nodes of a power grid follow from Eq. (2):(A.14)where the matrix M = (U+Y^{n}Z^{e})^{−1} depends on the electric parameters of the power grid. The elements of the vector J^{e} are defined according to Eqs. (4) and (5) by(A.15)
Until this point, we have not made any assumptions about geometry.
Combining ^{Eqs. (A.12)}, ^{(A.13)} and ^{(A.15)}, we obtain the following expression for J_{ji} ^{n} in spherical geometry:(A.16)where(A.17)
Note that c_{ij} = c_{ji}, but a_{ij} ≠ ±a_{ji} and b_{ij} ≠ ±b_{ji}. Coefficients c_{ij} are nonzero only if there is a transmission line between the nodes i and j. The corresponding matrices are thus sparse in practice.
A.1.4 Vectorised form in planar geometry
The approximation derived above for spherical geometry can be made for planar geometry too. We again assume that the electric field along a transmission line between nodes i and j equals the average of the fields at the nodes. It follows that the voltage between i and j is(A.18)
Substituting Eqs. (A.18) and (A.15) shows that the earthing current vector is(A.19)where M is defined in Eq. (A.14), and the kth element of vector A_{x} is(A.20)and a similar formula is obvious for A_{y}. In the practical computation, we can make use of the antisymmetry of (E_{xk} + E_{xl})(x_{k} – x_{1})/R_{lk}^{n} with respect to the indices k, l.
A.2 Finnish test grid
The old Finnish 400 kV grid in 1978–1979 serves as a realistic test model (Fig. A.2, Table A.1). Parameter values are based on Pirjola & Lehtinen (1985) and Pirjola (2009) with some modifications as explained in the caption of Table A.1. We also note that the present Finnish grid has many more stations and lines, and the total resistances are generally larger than in the old grid due to neutral point reactors in many transformers. There are also some series capacitors which block the flow of GIC.
Fig. A.2 Schematic map of the Finnish 400 kV power grid in 1978–1979 with transmission lines approximated by straight lines. 
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 Ω 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.
References
 Albertson, V.D., and J.M. Thorson Jr., Power system disturbances during a K8 geomagnetic storm: August 4, 1972, IEEE Trans. Power Apparatus Syst., PAS93, 1025–1030, 1974. [CrossRef] [Google Scholar]
 Amm, O., Ionospheric elementary current systems in spherical coordinates and their application, J. Geomag. Geoelectr., 49, 947–955, 1997. [CrossRef] [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 and Space, 51, 431–440, 1999. [Google Scholar]
 Bernhardi, E.H., P.J. Cilliers, and C.T. Gaunt, Improvement in the modelling of geomagnetically induced currents in southern Africa, South African J. Sci., 104, 265–272, 2008. [Google Scholar]
 Bolduc, L., GIC observations and studies in the HydroQuébec power system, J. Atmos. Sol.Terr. Phys., 64, 1793–1802, 2002. [CrossRef] [Google Scholar]
 Boteler, D., Distributed source transmission line theory for electromagnetic induction studies, Supplement of the Proceedings of the 12th International Zurich Symposium and Technical Exhibition on Electromagnetic Compatibility, 401–408, 1997. [Google Scholar]
 Elovaara, J., P. Lindblad, A. Viljanen, T. Mäkinen, R. Pirjola, S. Larsson, and B. Kielén, Geomagnetically induced currents in the Nordic power system and their effects on equipment, control, protection and operation, Paper 36301 of the Proceedings of the CIGRE 1992 Session, 30 August–5 September, 1992, Paris, 10 pp, 1992. [Google Scholar]
 Elovaara, J., Finnish experiences with grid effects of GICs. In: Space Weather – Research Towards Applications in Europe, ed. J. Lilensten, Springer, Dordrecht, The Netherlands, pp. 311–326, 2007. [CrossRef] [Google Scholar]
 Kappenman, J.G., An overview of the impulsive geomagnetic field disturbances and power grid impacts associated with the violent SunEarth connection events of 29–31 October 2003 and a comparative evaluation with other contemporary storms, Space Weather, 3, S08C01, 2005. [CrossRef] [Google Scholar]
 Kis, A., A. Koppan, I. Lemperger, T. Prodan, J. Szendröi, J. Verö, and V. Wesztergom, Longterm variation of the geoelectric activity index T, Publs. Inst. Geophys. Pol. Acad. Sci., C99(398), 353–360, 2007. [Google Scholar]
 Lanzerotti, L.J., Geomagnetic influences on manmade systems, J. Atmos. Terr. Phys., 41, 787–796, 1979. [CrossRef] [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]
 Liu, C.M., L.G. Liu, and R. Pirjola, Geomagnetically induced currents in the highvoltage power grid in China, IEEE Trans. Power Deliv., 24, 2368–2374, 2009. [CrossRef] [Google Scholar]
 Marshall, R., E. Smith, M. Francis, C. Waters, and M. Sciffer, A Preliminary Risk Assessment of the Australian region Power Network to Space Weather, Space Weather, 9, S10004, 2011. [CrossRef] [Google Scholar]
 McLay, S.A., and C.D. Beggan, Interpolation of externallycaused magnetic fields over large sparse arrays using spherical elementary current systems, Ann. Geophys., 28, 1795–1805, 2010. [CrossRef] [Google Scholar]
 Mäkinen, T. Geomagnetically induced currents in the Finnish power transmission system, Finnish Meteorological Institute Geophysical Publications, No. 32, 101 p, 1993. [Google Scholar]
 Parker, R.L., The magnetotelluric inverse problem, Geophys. Surv., 6, 5–25, 1983. [CrossRef] [Google Scholar]
 Pirjola, R., Geomagnetically induced currents in the Finnish 400 kV power transmission line, Phys. Earth Planet. Inter., 53, 214–220, 1989. [CrossRef] [Google Scholar]
 Pirjola, R., Review of the calculation of surface electric and magnetic fields and of geomagnetically induced currents in groundbased technological systems, Surv. Geophys., 23, 71–90, 2002. [CrossRef] [Google Scholar]
 Pirjola, R., Effects of space weather on highlatitude ground systems, Adv. Space Res., 36, 2231–2240, 2005. [CrossRef] [Google Scholar]
 Pirjola, R., Effects of interactions between stations on the calculation of geomagnetically induced currents in an electric power transmission system, Earth Planets Space, 60, 743–751, 2008a. [Google Scholar]
 Pirjola, R., Study of effects of changes of earthing resistances on geomagnetically induced currents in an electric power transmission system, Radio Sci., 43, RS1004, 2008b. [CrossRef] [Google Scholar]
 Pirjola, R., Properties of matrices included in the calculation of geomagnetically induced currents (GICs) in power systems and introduction of a test model for GIC computation algorithms, Earth Planets Space, 61, 263–272, 2009. [Google Scholar]
 Pirjola, R., and M. Lehtinen, Currents produced in the Finnish 400 kV power transmission grid and in the Finnish natural gas pipeline by geomagneticallyinduced electric fields, Ann. Geophys., 3, 485–491, 1985. [Google Scholar]
 Pulkkinen, A., A. Viljanen, K. Pajunpää, and R. Pirjola, Recordings and occurrence of geomagnetically induced currents in the Finnish natural gas pipeline network, J. Appl. Geophys., 48, 219–231, 2001a. [CrossRef] [Google Scholar]
 Pulkkinen, A., R. Pirjola, D. Boteler, A. Viljanen, and I. Yegorov, Modelling of space weather effects on pipelines, J. Appl. Geophys., 48, 233–256, 2001b. [CrossRef] [Google Scholar]
 Pulkkinen, A., O. Amm, A. Viljanen, and BEAR Working Group, Ionospheric equivalent current distributions determined with the method of spherical elementary current systems, J. Geophys. Res., 108 (A2), 1053, 2003a. [CrossRef] [Google Scholar]
 Pulkkinen, A., A. Thomson, E. Clarke, and A. McKay, April 2000 geomagnetic storm: ionospheric drivers of large geomagnetically induced currents, Ann. Geophys., 21, 709–717, 2003b. [CrossRef] [Google Scholar]
 Pulkkinen, A., S. Lindahl, A. Viljanen, and R. Pirjola, October 29–31, 2003 geomagnetic storm: geomagnetically induced currents and their relation to problems in the Swedish highvoltage power transmission system, Space Weather, 3, S08C03, 2005. [CrossRef] [Google Scholar]
 Pulkkinen, A., M. Hesse, S. Habib, L. Van der Zel, B. Damsky, F. Policelli, D. Fugate, W. Jacobs, and E. Creamer, Solar shield: forecasting and mitigating space weather effects on highvoltage power transmission systems, Nat. Hazards, 53, 333–345, 2010. [CrossRef] [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, 2005. [CrossRef] [Google Scholar]
 Trichtchenko, L., and D.H. Boteler, Specification of geomagnetically induced electric fields and currents in pipelines, J. Geophys. Res., 106, 21039–21048, 2001. [CrossRef] [Google Scholar]
 Trichtchenko, L., and D.H. Boteler, Modeling Geomagnetically Induced Currents Using Geomagnetic Indices and Data, IEEE Trans. Plasma Sci., 32, 1459–1467, 2004. [CrossRef] [Google Scholar]
 Trivedi, N.B., I. Vitorello, W. Kabata, S.L.G. Dutra, A.L. Padilha, M.S. Bologna, M.B. Padua, A.P. Soares, G.S. Luz, R. Pirjola, and A. Viljanen, Geomagnetically induced currents in an electric power transmission system at low latitudes in Brazil: a case study, Space Weather, 5, S04004, 2007. [CrossRef] [Google Scholar]
 Viljanen, A., European project to improve models of geomagnetically induced currents, Space Weather, 9, S07007, 2011. [CrossRef] [Google Scholar]
 Viljanen, A., A. Pulkkinen, O. Amm, R. Pirjola, 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. [CrossRef] [Google Scholar]
 Viljanen, A., A. Pulkkinen, R. Pirjola, K. Pajunpää, P. Posio, and A. Koistinen, Recordings of geomagnetically induced currents and a nowcasting service of the Finnish natural gas pipeline system, Space Weather, 4, S10004, 2006. [CrossRef] [Google Scholar]
 Watari, S., M. Kunitake, K. Kitamura, T. Hori, T. Kikuchi, K. Shiokawa, N. Nishitani, R. Kataoka, Y. Kamide, T. Aso, Y. Watanabe, and Y. Tsuneta, Measurements of geomagnetically induced current in a power grid in Hokkaido, Japan, Space Weather, 7, S03002, 2009. [CrossRef] [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, 2008. [CrossRef] [Google Scholar]
 Wik, M., R. Pirjola, H. Lundstedt, A. Viljanen, P. Wintoft, and A. Pulkkinen, Space weather events in July 1982 and October 2003 and the effects of geomagnetically induced currents on Swedish technical systems, Ann. Geophys., 27, 1775–1787, 2009. [CrossRef] [Google Scholar]
All Tables
GIC (amperes) due to a spatially constant northward or eastward electric field (1 V/km) in the Finnish highvoltage power grid in 1978–1979. “Plane” refers to the planar method and “sphere” to the spherical method.
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 Ω 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.
All Figures
Fig. 1 Prototype model of highvoltage power grids in Europe. 

In the text 
Fig. 2 Map of conductivity blocks in Europe. Each number denotes a different 1D ground model, whose parameters (layer thicknesses and resistivities) are listed on http://real.mtak.hu/2957/. Blocks outside of the presently modelled power grid are marked by 0. 

In the text 
Fig. 3 Magnetometer stations used in this study. Blue dots: IMAGE magnetometer network (http://space.fmi.fi/image/), black dots: sites available at the World Data Centre for Geomagnetism (Edinburgh) (http://www.wdc.bgs.ac.uk/) (also belonging to INTERMAGNET). 

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

In the text 
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. 

In the text 
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. 

In the text 
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 threelayer conductivity model was applied (block 3 in Fig. 2). 

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

In the text 
Fig. 9 GIC along the Finnish natural gas pipeline at the Mäntsälä compressor station on 15 March 2012. Black curve: measured value, blue: modelled (without data from the nearby Nurmijärvi magnetic observatory). Oneminute averages are shown. A twolayer conductivity model is used (block 25 in Fig. 2). 

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

In the text 
Fig. A.1 Vectors related to two nodes on a spherical surface. 

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

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.