© Owned by the authors, Published by EDP Sciences 2012

1 Introduction

The EURISGIC project (European Risk from Geomagnetically Induced Currents) will produce the first European-wide real-time 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; 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 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.

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):

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.

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$$E_x\left(\omega \right)=\frac{Z\left(\omega \right)}{{\mu }_0}{B}_y\left(\omega \right), E_y\left(\omega \right)=-\frac{Z\left(\omega \right)}{{\mu }_0}{B}_x\left(\omega \right),$$(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:$${\mathbf{I}}^e=(U+{\mathbf{Y}}^n{\mathbf{Z}}^e{)}^{-1} {\mathbf{J}}^e,$$(2)where U is a unit matrix and Yⁿ 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ⁿ 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ⁿ is defined by$$Y_{\mathrm{ij}}^n=-\frac{1}{R_{\mathrm{ij}}^n} \left(i\ne j\right),$$(3) $$Y_{\mathrm{ii}}^n=\sum _{k\ne i}\mathrm{} \frac{1}{R_{\mathrm{ik}}^n},$$where denotes the conductor line resistance between nodes i and j. The elements of the N × 1 column matrix J^e are defined by:$$J_i^e=\sum _{j\ne i}\mathrm{}{J}_{\mathrm{ji}}^n,$$(4) $$J_{\mathrm{ji}}^n=\frac{V_{\mathrm{ji}}}{R_{\mathrm{ji}}^n},$$(5) $$V_{\mathrm{ji}}={\int }_{s_{\mathrm{ji}}}\mathrm{}\mathbf{E}\cdot \mathbf{ds}.$$(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$$I_{\mathrm{ij}}^n{R}_{\mathrm{ij}}^n=V_{\mathrm{ij}}+\sum _{k,l}\mathrm{}{J}_l^e(Z_{\mathrm{ik}}^e-Z_{\mathrm{jk}}^e)\left[\right(U+{\mathbf{Y}}^n{\mathbf{Z}}^e{)}^{-1}{]}_{\mathrm{kl}}.$$(7)

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 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 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 Ω.

• 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 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.

Fig. 1 Prototype model of high-voltage 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 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 μ₀ everywhere and a constant permittivity ε = 5ε₀, where ε₀ 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 1-D 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 (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.

2.2.3 Interpolation of the magnetic field

We used 1-min 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 quiet-time baselines by subtracting the average of a 2-day sequence starting from the noon and lasting 48 h. The actual event under study is the full day in the middle of the 48-h 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” quiet-time 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 quick-look 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:$$\mathrm{AE}\left(t\right)=\max \left(\right\{X\left(t\right)\left\}\right)-\min \left(\right\{X\left(t\right)\left\}\right)$$(8)where {X(t)} refers to the value of the field at all sites at time t.

2. Largest horizontal difference field:$$\mathrm{d}H\left(t\right)=\max \left(\sqrt{\left\{\mathrm{d}X\right(t{)}^2+\mathrm{d}Y\left(t{)}^2\right\}}\right)$$(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 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.

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$$\mathrm{GIC}=a\cdot E_x+b\cdot E_y,$$(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 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.

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 three-layer 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 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 Ω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 two-layer 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 North-West 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). One-minute averages are shown. A two-layer 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 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 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).

• One-minute magnetic field data of 2003 from about 40 European magnetometer sites.

• Twenty-two different ground conductivity blocks.

• Ionospheric grid of 283 SECS poles.

• 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.

3 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.

Acknowledgments

The research leading to these results has received funding from the European Community’s Seventh Framework Programme (FP7/2007-2013) 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 long-term 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/2007-2013) 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:$$\mathbf{E}=E_{\theta }{\mathbf{e}}_{\theta }+E_{\varphi }{\mathbf{e}}_{\varphi }.$$(A.1)

Consider a transmission line between nodes 1 and 2. The voltage V₁₂ between these points is the line integral of the electric field:$$V_{12}={\int }_1^2\mathrm{}\mathbf{E}\cdot \mathbf{dr}.$$(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$$V_{12}=R{\int }_{{\varphi }_1}^{{\varphi }_2}\mathrm{}\left(E_{\theta }\left(\theta \left(\varphi \right),\varphi \right)\frac{\mathrm{d}\theta \left(\varphi \right)}{\mathrm{d}\varphi }+E_{\varphi }\left(\theta \left(\varphi \right),\varphi \right)\sin \theta \left(\varphi \right)\right)\mathrm{d}\varphi ,$$(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:$$\tan \theta =-\frac{C}{A\cos \varphi +B\sin \varphi },$$(A.4)where the constants A, B, C depend on the (θ, ϕ) coordinates of the points 1 and 2:$$\begin{array}{c}A=\sin {\theta }_1\cos {\theta }_2\sin {\varphi }_1-\cos {\theta }_1\sin {\theta }_2\sin {\varphi }_2,\\ B=-\sin {\theta }_1\cos {\theta }_2\cos {\varphi }_1+\cos {\theta }_1\sin {\theta }_2\cos {\varphi }_2 \\ C=\sin {\theta }_1\sin {\theta }_2\sin \left({\varphi }_2-{\varphi }_1\right).\end{array},$$(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₀e_θ (E₀ > 0, e_θ 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₁₂ = −E₀R(θ₂ − θ₁). The apparently simple uniform eastward field E_e = E₀e_ϕ has a non-zero vertical component of the curl, so the voltage is path-dependent. Instead, an eastward field E_e = E₀/sin θe_ϕ is associated with an path-independent voltage V₁₂ = −E₀R(ϕ₂ − ϕ₁).

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₁, r₂) at the Earth’s surface with polar coordinates (R, θ₁, ϕ₁) and (R, θ₂, ϕ₂) (Fig. A.1). We denote their difference vector by$$\mathbf{s}={\mathbf{r}}_2-{\mathbf{r}}_1.$$(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$$\mathbf{u}=\frac{\mathbf{s}\times {\mathbf{r}}_1}{|\mathbf{s}\times {\mathbf{r}}_1|}$$(A.7)is perpendicular both to r₁ and s, and tangential to the surface at r₁ = (R, θ₁, ϕ₁). It follows that the unit vector$${\mathbf{v}}_1=\frac{{\mathbf{r}}_1\times \mathbf{u}}{R}$$(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$$E_{\mathrm{par}1}={\mathbf{E}}_{\mathbf{1}}\cdot {\mathbf{v}}_1=\frac{R{\mathbf{E}}_1\cdot {\mathbf{r}}_2}{|{\mathbf{r}}_2\times {\mathbf{r}}_1|}.$$(A.9)

Similarly, we obtain at the point 2$$E_{\mathrm{par}2}={\mathbf{E}}_{\mathbf{2}}\cdot {\mathbf{v}}_2=- \frac{R{\mathbf{E}}_2\cdot {\mathbf{r}}_1}{|{\mathbf{r}}_2\times {\mathbf{r}}_1|}.$$(A.10)

The angle γ between the radius vectors r₁, r₂ is defined by$$\cos \gamma =\cos {\theta }_1\cos {\theta }_2+\sin {\theta }_1\sin {\theta }_2\cos \left({\varphi }_2-{\varphi }_1\right).$$(A.11)

In all practical situations, 0 ≤ γ ≤ π, so sin γ ≥ 0, and |r₂ × r₁| = R²sin γ. Then the average electric field E_par = (E_par1 + E_par2)/2 along the great circle arc between points 1 and 2 is:$$E_{\mathrm{par}}=\left[E_{\theta 1}\right(\cos {\theta }_1\sin {\theta }_2\cos ({\varphi }_2-{\varphi }_1)-\sin {\theta }_1\cos {\theta }_2)-E_{\theta 2}(\cos {\theta }_2\sin {\theta }_1\cos ({\varphi }_2-{\varphi }_1)-\sin {\theta }_2\cos {\theta }_1)+(E_{\varphi 1}\sin {\theta }_2+E_{\varphi 2}\sin {\theta }_1\left)\sin \right({\varphi }_2-{\varphi }_1\left)\right]/\left(2\sin \gamma \right).$$(A.12)

The voltage is$$V_{12}=\mathrm{R\gamma }{E}_{\mathrm{par}}.$$(A.13)

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):$${\mathbf{I}}^{\mathbf{e}}=(\mathbf{U}+{\mathbf{Y}}^{\mathbf{n}}{\mathbf{Z}}^{\mathbf{e}}{)}^{-1}{\mathbf{J}}^{\mathbf{e}}=\mathbf{M}{\mathbf{J}}^{\mathbf{e}},$$(A.14)where the matrix M = (U+YⁿZ^e)⁻¹ 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$$J_i^e=\sum _{j\ne i}\mathrm{}{J}_{\mathrm{ji}}^n=\sum _{j\ne i}\mathrm{}\frac{V_{\mathrm{ji}}}{R_{\mathrm{ji}}^n}.$$(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 ⁿ in spherical geometry:$$J_{\mathrm{ji}}^n=c_{\mathrm{ji}}\left(a_{\mathrm{ij}}{E}_{\theta j}-a_{\mathrm{ji}}{E}_{\theta i}+b_{\mathrm{ij}}{E}_{\varphi j}-b_{\mathrm{ji}}{E}_{\varphi i}\right),$$(A.16)where$$\begin{array}{c}{a}_{\mathrm{ij}}=\cos {\theta }_j\sin {\theta }_i\cos \left({\varphi }_i-{\varphi }_j\right)-\sin {\theta }_j\cos {\theta }_i,\\ b_{\mathrm{ij}}=\sin {\theta }_i\sin \left({\varphi }_i-{\varphi }_j\right),\\ c_{\mathrm{ij}}=\frac{-R{\gamma }_{\mathrm{ij}}}{2\sin {\gamma }_{\mathrm{ij}}}\frac{1}{R_{\mathrm{ij}}^n}.\end{array}$$(A.17)

Note that c_ij = c_ji, but a_ij ≠ ±a_ji and b_ij ≠ ±b_ji. Coefficients c_ij are non-zero 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$$V_{\mathrm{ij}}=(E_{\mathrm{xi}}+E_{\mathrm{xj}})\cdot (x_j-x_i)/2+(E_{\mathrm{yi}}+E_{\mathrm{yj}})\cdot (y_j-y_i)/2.$$(A.18)

Substituting Eqs. (A.18) and (A.15) shows that the earthing current vector is$${\mathbf{I}}^e=\mathbf{M}\left({\mathbf{A}}_{\mathbf{x}}+{\mathbf{A}}_{\mathbf{y}}\right),$$(A.19)where M is defined in Eq. (A.14), and the kth element of vector A_x is$${\mathbf{A}}_{\mathit{x}k}=\sum _{l=1,l\ne k}^N\mathrm{}\frac{(E_{\mathrm{xk}}+E_{\mathrm{xl}})(x_k-x_l)}{2R_{\mathrm{lk}}^n}$$(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₁)/R_lkⁿ 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.

References

All Tables

All Figures

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

	Fig. 2 Map of conductivity blocks in Europe. Each number denotes a different 1-D 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 three-layer 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 two-layer 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). One-minute averages are shown. A two-layer 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