© A. Gil et al., Published by EDP Sciences 2021

1 Introduction

Powerful solar phenomena are the drivers of large-scale time variations of currents flowing in the Earth’s ionosphere and magnetosphere (e.g., Pulkkinen et al., 2012; Carter et al., 2015). Magnetic variations during geomagnetic storms tend to be larger at higher latitudes. However, very powerful space weather events have caused intensive geomagnetic disturbances also at the middle and low latitudes. As an example, we can mention the damages to some transformers, e.g., in the South African region as a result of the Halloween storm (Gaunt & Coetzee, 2007) or in New Zealand during the November 2001 storm (Marshall et al., 2012). Lotz & Danskin (2017) showed that for the middle-latitude stations, the Dst geomagnetic index can be used as a proxy for the cumulative induced E field.

Recently, the functioning of electrical power networks under the influence of space weather effects was extensively studied also in Europe: in the Czech Republic (Výbošt’oková & Švanda, 2019; Švanda et al., 2020), in Italy (Tozzi et al., 2018, 2019; Piersanti et al., 2020), in Greece (Zois, 2013), in Spain (Torta et al., 2012) and in Austria (Bailey et al., 2018). Besides European studies, Schrijver et al. (2014) found an increase in the number of insurance claims in United States around time of space weather events.

Švanda et al. (2020) and Výbošt’oková & Švanda (2019) showed that the 5–10% increase of failures in the Czech power grid recorded in the 5-day period following the enhanced geomagnetic activity is possibly linked to geomagnetically induced current (GIC) effects. Švanda et al. (2020) presented that failures in power lines were growing immediately (within 1 day). Failures regarding transformers were generally delayed by 2–3 days, still allowing the operation of power systems without registering problems until the cumulative damage reached a sufficient level (more details in Kappenman, 2013; Wang et al., 2015).

Tozzi et al. (2018) reconstructed a latitudinal profile of maximum GIC indices for several geomagnetic storms for different locations. They showed that the latitudinal variability of GIC intensity is almost the same for all considered storms. Tozzi et al. (2019) studied the possible impact of space weather events on the Italian territory by computing the GIC index for the 2015 St. Patrick’s Day storm. Detrimental GICs could flow through the power network, especially at the highest Italian latitudes characterized by a low conductivity lithosphere. Zois (2013) found the correlation between solar activity and the annual number of transformer failures in Greece for 1989–2010.

The main impact of GIC is system instability. It appears due to high levels of VAR (volt-amp reactive) and significant current harmonics as a result of the transformer core part – cycle semi – saturation (Girgis & Vedante, 2012). GIC as a function of time can be coupled with transformer temperature models to follow the winding hot spot temperatures throughout the storm. It is reasonable to expect that most of the transformer failures during an extreme geomagnetic storm would come from the population that has exceeded a critical temperature threshold (Oughton et al., 2019). The hotspot temperature profile simulation for specific transformers impacted by GIC based on the fitted analytical function was proposed by Marti et al. (2013).

Another critical issue is the loss of transformer life due to overheating during a storm. Even a transformer whose winding hotspot remained below the critical temperature threshold may sustain insulation damage, which reduces the transformer’s lifetime. This additional loss-of-life in years is essential as the reason for aging electrical equipment. Threats of aging electrical equipment and components reduce their efficiency, being potentially dangerous. Overheating in electrical enclosures can shorten the life of the equipment or worse, leading to costly downtime (Výbošt’oková & Švanda, 2019).

In our previous paper (Gil et al., 2020a), we computed the geoelectric field for Poland using the method by Viljanen & Pirjola (1989). Basing on the data of electrical network failures in Southern Poland in 2010 and the first seven months of 2014 and using various geomagnetic storm indicators (e.g., sudden storm commencement, Kp ≥ 5 and fast halo coronal mass ejecta appearance) we showed a grow in the number of failures around that time. In this paper, we study the influence of intense geomagnetic storms on the functioning of electrical power grids in Poland for the period 2010–2014 (the first half of the solar cycle 24). Information about failures in power networks is sensitive data. Therefore, we have access to a limited time span of data. We analyse here the relation of transmission line (network) failures with heliospheric and geomagnetic data. What is more, the temporal evolution of geoelectric field estimates are determined using a 1-D layered conductivity model. It has to be underlined that the first computation of geoelectric field for Poland based on a 1-D layered conductivity model was performed by Viljanen et al. (2014), but for a different period, namely years 1996–2008 (solar cycle 23). Our results reveal growth in the number of transmission network failures during and right after considered classes of geomagnetic storms. This allows us to speculate about the sensitivity of electrical power grids in Poland to solar activity conditions and to underline the importance of further studies on this topic.

2 Data

2.1 Intense geomagnetic storms in 2010–2014

We considered here data characterising the overall situation near and on the Earth in 2010–2014. Special attention was paid to data of Bz component of the heliospheric magnetic field (HMF) and Dst-index, with 1-hour resolution between January 2010 and July 2014. The Bz component of the HMF used in the study is in the Geocentric Solar Magnetospheric System (GSM), with X-axis pointed from the Earth to the Sun, the Y-axis being perpendicular to the Earth magnetic dipole (the XZ-plane contains the dipole axis), and the positive Z-axis is at the same meaning as the North magnetic pole.

Following Gonzalez & Tsurutani (1987), we extracted from the whole set of data all intense geomagnetic storms, i.e., the cases when Bz < −10 nT and Dst < −100 nT for more than 3 h. During January 2010–July 2014 were registered 16 events with Bz < −10 nT for longer than 3 h (it is designated as case 1), 8 events with Dst < −100 nT for longer than 3 h (case 2), and 5 events fulfilling both conditions simultaneously, i.e., Bz < −10 nT and Dst < −100 nT during more than 3 h (case 3). All events are listed in Table 1. When an event’s duration was intermittent, we denote its duration as a sum (+) of hours. Some of these storms were discussed in the context of geomagnetic events with intensified GIC risk (Bailey & Leonhardt, 2016). Although the Dst-index is not the best indicator of GIC hazard (Cid et al., 2014), we use it here combined with the HMF Bz component for the severity of geomagnetic storm assessment (Gonzalez et al., 1994).

In this work frame, we also considered two components of the geomagnetic field B_X, B_Y measured in Belsk, the Polish INTERMAGNET observatory. In summary, 5-years of measurements at 1-minute intervals have been considered. Additionally, the first derivative of the geomagnetic field horizontal () component B_X, dB_X/dt has been considered as another indicator of strong magnetic storms (e.g., Adebesin et al., 2016; Love et al., 2016; Ngwira et al., 2018).

Moreover, we analyse heliospheric parameters, a strength of HMF B [nT] and its B_Z [nT] component, solar wind speed SWs [km/s], and proton density SWd [N/cm³], as well geomagnetic indices: ap-index [nT], AE-index [nT], local K-index from Belsk. The maximal values of the above-mentioned parameters during the analysed geomagnetic storms are listed in Table 1.

2.2 Transmission line failures

In 2019, electricity production in Poland was falling, with the most significant decrease seen in production from lignite and hard coal. The share of coal in electricity production in 2019 was 73.6%, 4.8 percentage points less than in 2018. In 2019, the import of electricity to Poland almost doubled, amounting to 10.6 TWh (e.g., Macuk, 2020). There are three largest producers in the Polish energy market (which are parts of the groups: PGE Polska Grupa Energetyczna S.A., TAURON Polska Energia S.A., ENEA S.A.) having in total almost 2/3 of the installed capacity and being responsible for almost 70% of domestic electricity production. TAURON Capital Group (TCG) is the largest Polish electricity distributor, which data we are using here (Annual Report, 2019).

We utilize around half of the million transmission line/network operator logs about electrical grid failures (EGFs), starting from January 2010 to July 2014. EGFs causes were aggregated into six general groups denoted by A–F (details in Gil et al., 2019). Comparable categorization was used by Zois (2013). The EGFs occurrence rate in years 2010–2013 and during the first seven months of 2014 is presented in Table 2. Table 2 shows the percentage of EGFs in each group from A to F in individual years, and collectively for the whole period of data availability, i.e., in 2010–2014.

3 Methodology

3.1 Statistical relation

Among statistical methods that can be used to compare two sets of data, there are, for example, statistical tests that verify the hypothesis of the equality of means or variances of these sets. Here, we applied the t-test (see also the t-Student test) (Shao, 2003), as the most common and easy method to evaluate the differences in means between two groups of data. The first group of our data determines the number of DEF failures that occurred in time after such geomagnetic events from 2010 to 2014, in which the maximal value of time derivative of B_X is not less than 25 nT/min, i.e. max(dB_X/dt) ≥ 25 nT/min. The second group determines the number of DEF failures which occurred in times before max(dB_X/dt) ≥ 25 nT/min (so in the second group we observe the failures of the type DEF that max(dB_X/dt) < 25 nT/min). Let us assume that there are no significant differences between their values in the two groups of data. Therefore, H₀: m₁ = m₂, where m₁ and m₂ express the means of the number of DEF failures of the first and second group, respectively, whereas H₁: m₁ ≠ m₂ is the alternative hypothesis. From the t-test distribution table, one takes a critical value t_α,n for a given level of significance α and n degrees of freedom. The critical set for the t-test is (−∞, −t_α,n] ∪ [t_α,n, ∞). If t < |t_α,n| we do not reject hypothesis H₀, otherwise we accept H₁. The value of p, given in the t-test expresses the probability of error involved in accepting hypothesis H₁ about the existence of a difference.

Theoretically, the t-test can be applied even if the samples contain only a few observations, as long as the two assumptions are satisfied, i.e., the variables have normal distributions within each group, and the variance of scores in both groups is not substantially different. The first assumption on the “normality” of the data can be estimated by looking at the data distribution histograms or performing a normality test. Here we have applied the Shapiro - Wilk test (for small statistical samples). The second one can be verified by using the F-test (see also the test of F – Fischer) or by using the Levene test, as well as the Brown-Forsythe test (which is a modification of the Levene test) (Shao, 2003).

3.2 Calculation of the geoelectric field

In the next step of the analysis to determine geoelectric field {E_X, E_Y} a layered conductivity Earth model has been applied (see, e.g., Boteler et al., 2019; Boteler & Pirjola, 2019, and references therein). This 1-D model considers the change of conductivity with depth in the Earth, representing it by N horizontal layers with specified conductivities σ_n and thicknesses l_n (n = 1, ..., N). The recursive relation expresses the layered case of transfer function K (in frequency domain f) (Weaver, 1994; Boteler & Pirjola, 2019):

(1)where K_n denotes the ratio of E to B at the top surface of layer n, while K_n+1 at the top surface of the underlying layer n + 1, , and μ₀ = 4π × 10⁻⁷ H/m (Boteler et al., 2019). The initial value in equation (1) is the case when the layer n = N is a uniform half-space and , while the final value K₁ is the transfer function relating E and B at the Earth’s surface (Trichtchenko & Boteler, 2002; Boteler et al., 2019).

To perform analysis of data of geomagnetic horizontal field from Belsk station (see Sect. 2.1) we used a Earth model number 45 of Ádam et al. (2012), dedicated for Poland and consisted of 4 layers with, from the top down, thicknesses and resistivities (1/σ): 6 km, 5 Ωm; 105 km, 1000 Ωm; 300 km, 100 Ωm; above a half space of 10 Ωm. This Earth-layered model for the European regions has already been used in a number of studies (e.g., Viljanen et al., 2013, 2014). Differences concerning various methods were emphasized by (e.g., Viljanen et al., 2012; Kelbert & Lucas, 2020).

In the next step, by taking the Fourier transform we decomposed the 1-minute geomagnetic field {B_X, B_Y} data into its frequency components {B_X(f), B_Y(f)}. Then, each of these components has been multiplied by the corresponding transfer function values giving the geoelectric field frequency components, namely, E_X(f) = K(f)B_Y(f) and E_Y(f) = −K(f)B_X(f). Finally, we took the inverse Fourier transform to obtain the geoelectric field E(t) in the time domain (Boteler, 2012; Boteler & Pirjola, 2019). It is worth stressing that each step of the presented here procedure has been implemented in Matlab software and checked using the synthetic magnetic field variation considered in (Boteler, 2012; Boteler & Pirjola, 2019). What is more, to perform the annual data analysis, the numerical calculations were carried out in the 3-day moving window. However, only results for the mid-day have been finally saved and used in further studies.

3.3 Superposed epoch analysis

Superposed epoch analysis is a method manifesting links between studied time series (Chree, 1913). This method application allowed to prove the existence of a 27-day variation in geomagnetic storm appearance (Chree, 1913; Chree & Stagg, 1927), further showing that it is only the case of weaker geomagnetic storms and the intense storms occur sporadically (Greaves & Newton, 1929). Denton et al. (2005) analysing relations between the geomagnetic storm phase with temporal fluctuations of plasma found at geosynchronous orbit, presented that for the plasma sheet density, one of the key-factor is the solar cycle. Superposed epoch analysis was used in the Forbush decreases of galactic cosmic rays studies (Singh & Badruddin, 2006; Okike & Umahi, 2019). Liemohn et al. (2008) explored an occurrence time or magnetic storms strength as their features. Gil et al. (2019) showed an increase in the superposed averaged number of EGF appearing around one day after the fast halo coronal mass ejection (CME) occurrence, on the day of sudden storm commencement, as well as around zero-day or the day after when the daily maximal Kp geomagnetic index value was greater or equal to 5.

Here we defined a pivot time based on the Table 1 as the day when there appeared the B_Z < −10 nT (case 1), Dst < −100 nT (case 2) and both, B_Z < −10 nT and Dst < −100 nT (case 3) during more than 3 h. Next, we extracted subsets of daily sums of electrical network failures (of a particular type) during the whole week around the pivot events, i.e., 3 days before and after each event. Subsequently, we superposed all extracted subsets of failures, synchronizing all zero-days and normalising daily values to the average during the whole week being equal to 0. We use this normalisation to highlight the percentage increase in the number of failures.

3.4 Self-organizing maps

Self-organizing maps (SOM) are a type of artificial neural network trained with unsupervised learning. The goal is to display training data’s input space as a low-dimensional, usually one or two-dimensional discretized map (Kohonen, 1990). SOM networks have only two layers: input and output layers. The output layer is represented as a map showing the topological proximity of mapped data points in a given space as a function of the number of neurons (Somervuo & Kohonen, 1999). Similar inputs are grouped in close neurons. This allows finding links between input data. Each input layer neuron is connected to each output layer neuron. The connections between these neurons are assigned a certain weight. Unlike supervised learning, where the winning neuron “takes it all”, the winning neuron’s neighborhood is also updated in self-organizing maps (Hu et al., 2019). The Kohonen rule is used for updating, as follows:

(2)where w_ij is the weight j of neuron i, α(t) is the learning rate, h_ci(t) is the neighborhood function (Caldas et al., 2017) that determines the neighborhood of the winner. The neighborhood is estimated using a distance function between neighbors, which may be the Euclidean function, the Gauss function, the Manhattan distance function, etc.

SOMs are mainly used for data clustering and classification problems because a trained SOM aggregates similar data elements into clusters (Lampinen & Oja, 1992; Caldas et al., 2017). In our case, the SOM has been used to group selected heliospheric and geomagnetic parameters (see Sect. 2.1) with transmission line failures (see Table 2). We focused on failures originating from the aging of the infrastructure elements (group D), linked to electronic devices (group E), and having unknown reasons (group F).

4 Analysis results

4.1 Statistical analysis results

Firstly, we checked if there is an overall trend in the number of electrical grid failures, being (to some extent) of solar origin. Thus, we have concentrated on failures originating from the aging of the infrastructure elements (group D), linked to electronic devices (group E), and having unknown reasons (group F). Figure 1 presents the annual percentage rate of EGFs from the three groups D–F (a sum of these three groups we denote as DEF failures). Error bars are calculated as RMSE from monthly values for each year (Table 2). One can see that the percentage fraction of EGFs from groups DEF has a rising trend similar to the solar activity cycle in 2010–2014, represented by yearly changes of sunspot number (visualised by orange squares). DEF’s percentage rate concerning all EGFs changes from (23 ± 2)% near the solar minimum in 2010 up to (30 ± 5)% near solar maximum in 2014.

Fig. 1 Linear regression of the annual percentage rate of EGFs from groups DEF in January 2010–July 2014 (circles, left axis) and yearly changes of sunspot number (orange squares, right axis).

Next, to compare dB_X/dt and failures of transmission line elements caused by selected reasons (groups D, E and F from Table 2) we applied the t-test (see Sect. 3.1). In our analysis we considered 3 h intervals of dB_X/dt data and failures caused by aging, electronic devices breakdowns, or having an unidentified reasons. We have taken into account such the geomagnetic events from 2010 to 2014 that max(dB_X/dt) ≥ 25 nT/min (e.g., Adebesin et al., 2016) and by the first group of our data, we have set these values of DEF failures, which observed in times directly after storm events. The second group of the data consisted of the aging, electronic breakdowns, and unidentified type of failures that occurred immediately before the geomagnetic storm, i.e., before the dB_X/dt max value was more or equal to 25 nT/min. In this way, we have got the sets of statistical parameters such as sums, means, medians, standard deviations, coefficient of variations, skewnesses, kurtosis, etc., in the intervals of time 24, 48, and 72 h. We have denoted them, for example, as mean 24 h− for a mean of the number of DEF failures that occurred during 24 h before max(dB_X/dt) ≥ 25 nT/min and mean 24 h+ for the mean of the number of DEF failures occurred during 24 h after the whole time interval with max(dB_X/dt) ≥ 25 nT/min. Similarly, 48 h− and 48 h+, 72 h− and 72 h+ work in the same manner and standard deviations are included alongside the mean. Using Statistica software, we compared the sets of DEF failures in both intervals of time 24 h− and 24 h+, 48 and 72 h. The summary of the obtained results is given in Tables 3 and Table 4.

In Tables 3 and 4, the values of: the Shapiro - Wilk test, W, the Levene test, F and t-test were obtained for our data by using corresponding formulas. Their values are not as significant as the p-values which are responsible for accepting or rejecting the H₀ hypothesis. In our calculations, we put α = 0.05, for α < p we accept H₀ hypothesis, otherwise, when α > p we reject H₀ and accept H₁ which is the alternative hypothesis, moreover, when α < p ≤0.1 we can also reject H₀ and accept H₁. From Tables 3 and 4, we can see that both assumptions for the t-test are fulfilled if we have 48 h or 72 h intervals. The comparison of means of DEF failures in 2010–2014 years for 48 h has shown that mean 48 h+ > mean 48 h− (see, mean 48 h+ = 9.17 and mean 48 h− = 5.91). It was confirmed by the t-test for means of DEF failures for 48 h (t = 2.21 with p = 0.034). Similarly, we considered the means of DEF failures denoted for 72 h, and we have obtained mean 72 h+ > mean 72 h− (see, mean 72 h+ = 7.23 and mean 72 h− = 4.31). The t-test gave t = 2.13 with p = 0.04. Analogically, for standard deviation we obtained SD 48 h+ > SD 48 h−, as well mean 72 h+ > mean 72 h− (see t = 1.88 with p = 0.07 and t = 1.92 with p = 0.06, respectively).

4.2 Computed geoelectric field

Figure 2 shows the geoelectric field components E_X and E_Y calculated based on 5 years measurements of the geomagnetic field and by applying 1-D layered conductivity Earth model (see Sect. 3.2). One can see E_X and E_Y in mV/km calculated for particular years of 2010–2014. Similar scales on Y-axis allow us to quickly compare the state of the geoelectric field between considered years. In particular, we see that year 2012 (fifth and sixth panels of Figure 2) was one of the most variable periods. This is in agreement with information from Table 1, where the geomagnetic storms have been listed mostly for the year 2012. The more detailed analysis of Figure 2 reveals the increase of E_X and E_Y values in days 69, 114–115, 197–198, 275 of the year 2012. Selected dates correspond to the most intensive geomagnetic storms (Dst < −100 nT and HMF B_Z < −10 nT).

Fig. 2 Geoelectric field components EX and EY for years 2010–2014 calculated using a 1-D layered-Earth model (Ádam et al., 2012).

We are aware that the induced E field can vary widely even during a given event (Ngwira et al., 2015; Pulkkinen et al., 2015). Therefore, a quantitative and more detailed analysis of computed geoelectric field E for selected geomagnetic storms has been performed in the next section.

4.3 Superposed epoch analysis results

Here we defined a pivot time based on the Table 1 as the days when there appeared the HMF southward component B_Z < −10 nT (case 1), Dst < −100 nT (case 2) and both, B_Z < −10 nT and Dst < −100 nT (case 3) during more than 3 h. Next, we extracted subsets of daily data of electrical network failures (of a particular type) during a whole week around the pivot events, i.e., 3 days before and after each event. Subsequently, we superposed all extracted subsets of failures synchronizing all zero-days and normalising values to the average during the whole week being equal to 0. Figures 3–5, panels (a) show that around a pivot event, huge growth in the EGF in all three groups is observed. To visualise the overall situation during the considered geomagnetic storms in Figures 3–5, panels (b) we display the temporal behaviour of superposed daily mean of Dst-index (averaged over the number of events, details in Table 1), panel (c) presents superposed daily mean of B_Z HMF component changes (averaged over the number of events) and panel (d) depicts superposed daily mean of computed geoelectric field [mV/km] (details in the Sect. 3.2).

Fig. 3 % above normalised sum of daily electrical grid failures (a), averaged over the number of events (b) Dst-index [nT], (c) BZ [nT] component of HMF and (d) computed geoelectric field E [mV/km], during the geomagnetic storms with BZ < −10 nT (case 1), for years 2010–2014.

Figure 3 displays case 1 circumstances when B_Z < −10 nT for more than 3 h. Figure 3, panel (a) presents that on the zero-day EGF caused by the aging of infrastructure elements (group D) increased to 52% above the average value during the whole week and rose to 44% above this average on the second day after the pivot time. The EGF growth connected to the unreliability of electronic devices (group E) on the 0 day was 263% above the average, on the day after reaching 434% and 193% above the average on the second day. EGF having unknown reasons (group F) grew from the mean value by 148% on the 0 day.

Figure 4 presents the conditions for the case 2, when Dst < −100 nT remained for more than 3 h. Figure 4, panel (a) presents that EGF caused by the aging of infrastructure elements (group D) increased to 247% above the average value during the whole week on the second day after the pivot time. The EGF growth connected to the unreliability of electronic devices (group E) was 86% on the day after the 0 day and still growing on the second day after the pivot time to 235% above the average. EGF having unknown reasons (group F) grew from the mean value by 104% on the 0 day and still being 50% above the average on the day after.

Fig. 4 % above normalised sum of daily electrical grid failures (a), averaged over the number of events (b) Dst-index [nT], (c) BZ [nT] component of HMF and (d) computed geoelectric field E [mV/km], during the geomagnetic storms with Dst < −100 nT (case 2), for years 2010–2014.

Figure 5 depicts the situation for the case 3, when B_Z < −10 nT and Dst < −100 nT lasted more than 3 h. Figure 5, panel (a) presents that EGF caused by the aging of infrastructure elements (group D) increased to 152% above the average value during the whole week on the second day after the pivot time. The EGF growth connected to the unreliability of electronic devices (group E) also appeared on the second day after the pivot time, reaching 520% above the average. EGF having unknown reasons (group F) grew from the mean value by 111% on the 0 day.

Fig. 5 % above normalised sum of daily electrical grid failures (a), averaged over the number of events (b) Dst-index[nT], (c) BZ [nT] component of HMF and (d) computed geoelectric field E [mV/km], during the intense geomagnetic storms with Dst < −100 nT and BZ < −10 nT (case 3), for years 2010–2014.

4.4 Self-organizing maps analysis results

SOM were applied to investigate the impact on the number of electrical grid failures of the selected heliospheric and geomagnetic data. Twelve parameters were selected for the analysis. Four of them are the number of electrical grid failures: caused by the aging of the infrastructure elements (group D), electronic devices (group E), having unknown reasons (group F), and the sum of groups D-F. The remaining 8 are the following data: strength of HMF B [nT] and its B_Z [nT] component, solar wind speed SWs [km/s] and proton density SWd [N/cm³], as well geomagnetic indices: Dst-index [nT], ap-index [nT], AE-index [nT], K-index.

We considered here the intense geomagnetic storms between January 2010 and July 2014 (case 3 described in the Sect. 2.1, details in the Table 1). Because of various delays between geomagnetic disturbances and transmission line failures appearance described in the literature (e.g., Zois, 2013; Gil et al., 2019, 2020a,b; Švanda et al., 2020), we considered here data sets where EGFs are: zero, one, two, and three days delayed in relation to the heliospheric and geomagnetic parameters. We compared our results with computations for an undisturbed period at the beginning of 2010 (control data set).

The above-described data sets with a 3-hour resolution were selected as an input. Various sizes of the map, the different number of iterations, and the neighborhood distance function were tested. The optimal set of computations parameters were a 4 × 4 neuron grid, 100 iterations, and the link distance function.

Our analysis showed that the most substantial relationship between the number of failures and considered heliospheric and geomagnetic data was observed for the data set without delay. Figure 6 presents an example of our analysis results for the intense geomagnetic storm on 15 July 2012. The SOM neighbor weight distances describe the weight assigned to the adjacent neurons’ distances, labeled as blue hexagons. The red lines show which particular neurons are connected and the colors in the areas containing the red lines show the weight of the connection between neurons. The lighter the color, the shorter the distance between adjacent neurons. Figure 6 also displays numerical representation of the adjacent neurons’ connection weights. A weight close to 1 means a shorter distance between neurons and more robust grouping into clusters in the studied space. The connection weights between neurons close to 0 mean that they are located far apart in the studied space. The parameters related to these neurons are not correlated with each other.

Fig. 6 SOM neighbor weight distances with weight values of connections between neighboring neurons for data without delay for the intense geomagnetic storm, 15 July 2012. The blue hexagons represent neurons, and the red lines show which particular neurons are connected. Colors from black to yellow display the weight values of the connection between neighboring neurons.

In Figure 6, the parameters related to the individual neurons are plotted. They represent the distribution of data in the studied space. The Dst geomagnetic index and B_Z component of HMF related to the first neuron in the lower-left corner have been grouped into one cluster located at some distance in the studied space from the other parameters constituting the second cluster. This means that the Dst and B_Z parameters are strongly connected, but not so with the others. In the lower-right corner are grouped parameters related to the number of electrical grid failures caused by aging of the infrastructure elements, group D (eighth neuron counting from the lower-left corner by rows), failures related to electronic devices, group E (seventh neuron), failures having unknown reasons, group F (third neuron) and the sum of groups DEF. The proton density is located near neurons related to the number of electrical grid failures at a distance of 0.82. Additionally, the geomagnetic ap-index, HMF strength B, and AE-index are located in the upper-right corner, while the solar wind speed and K-index are located in the upper-left corner. The light color of the connections between the neurons related to the parameters above proves that the SWd, B, and ap-index parameters significantly impact the number of electrical grid failures, especially on groups D and E. The AE-index, K-index, and SWs parameters also impact the occurrence of failures due to the grouping of these parameters in the same cluster. However, the greater distance from the neurons related to electrical grid failures expresses their lower impact.

5 Discussion

The timeline of the yearly EGFs from groups DEF, which might have the solar origin, in 2010–2014 presents the rising linear trend connected with the more frequent intense geomagnetic storms in the rising phase of solar activity cycle 24.

Superposed epoch analysis results show that on the day of and just after the intense geomagnetic storms, there is observed a significant growth in the number of electrical grid failures in South Poland. The growth is visible in the groups of failures that might be connected with space weather effects, i.e., EGFs triggered by the aging of infrastructure elements, associated with electronic devices’ unreliability and having unknown causes. What is more, group F’s failures (of unknown reasons) seem to occur almost immediately (days + 0), while the others have an apparent delay of 1–2 days.

Moreover, SOM analysis shows that proton density SWd had the most decisive impact on the number of transmission line failures. Moreover, the strength of HMF B and ap-index had a strong influence on the number of failures, mainly those from groups D and E. The most significant impact of the remaining heliospheric and geomagnetic parameters was observed for data without delay.

The GIC related disturbances have an additional effect in the power network related to the response of protection devices such as differential relays. The GIC flows in the electrical grid are not only connected to the transformers overheating, but the spiked-narrow pulse of extra energy can also produce higher harmonics in a transformer. As a result, transient states and resonance in the power grid may occur, differential relays may enable protection, and the grid stability may be compromised (Girgis & Vedante, 2012; Power System Relaying and Control Committee, 2019). Harmonics from GIC may affect the operation of power generators (Meliopoulos et al., 2018) and propagate into other generators and transformers. The power grid response to the transients and surges is unpredictable. There is a possibility of uncontrolled propagation of GIC origin transients effects into other segments of the power network – therefore, more investigations are needed to determine the effects of GIC in technical infrastructure.

6 Conclusions

In this work, we performed a systematic study of the influence of the most intense geomagnetic storms on the functioning of electrical power grids in Poland in 2010–2014, which is around a half of the solar cycle.

The main findings of this work can be summarized as follows:

• yearly average number of transmission line failures in South Poland between January 2010 (near the solar minimum) and July 2014 (around the solar maximum) show a rising trend, which can be the indication of solar cycle phase dependency;

• during and right after the intense geomagnetic storms, there is observed a massive increase in the number of electrical grid failures in South Poland, in the groups of failures which might be connected with space weather effects, namely transmission line failures caused by the aging of infrastructure elements, connected to the unreliability of electronic devices and having unknown reasons (groups DEF);

• t-test was used to study the difference between means and standard deviations of two groups of DEF failures (which might be related with space weather): the first for periods when max(dB_X/dt) ≥ 25 nT/min and the second for max(dB_X/dt) < 25 nT/min, in 2010–2014, with a delay of two (48 h) and three days (72 h), showed that there exists the growth of the number of aging, electronic breakdowns and unidentified type of failures for the group with max(dB_X/dt) ≥ 25 nT/min;

• increases of geoelectric field, computed for Poland region by using 1-D layered Earth model, correspond to the intense geomagnetic storms occurrence;

• proton density SWd seems to have the most significant impact on the number of power grid failures in the case without delay concerning the intense geomagnetic storms;

• other analyzed heliospheric and geomagnetic data (ap, AE, K indices, as well as solar wind speed and HMF strength) also show a varying degree of impact on the number of power grid failures.

Summarizing, we see that during and right after intense geomagnetic storms, an increase in the number of transmission line/network failures, which might be of solar origin, occurs.

We plan to continue the further detailed analysis of selected intense geomagnetic storms. In particular, further studies based on a wider time span are still needed and we are convinced that they can shed more light on the considered problem. Moreover, we plan to apply the determined geoelectric field for the systematic studies of GIC. There is a need for future engineering analysis for quantitative failure studies to understand their possible direct relationship to geomagnetically induced currents. Another exciting aspect of analysis will be the improvement of the method of geoelectric field determination. We are aware that 1-D layered Earth model conductivity used in this study changes vertically but not horizontally. Two-dimensional and three-dimensional Earth conductivity structure introduces some features not seen with 1-D models (Marti et al., 2014; Boteler & Pirjola, 2017). Therefore, the forthcoming study will take into account the verification of other Earth modeling assumptions.

Acknowledgments

Heliospheric, solar wind parameters and geomagnetic indices are from http://omniweb.gsfc.nasa.gov, geomagnetic field components from Belsk Observatory, a part of INTERMAGNET, http://rtbel.igf.edu.pl/. We acknowledge the Polish National Science Centre’s financial support, grant number 2016/22/E/HS5/00406. The editor thanks Ari Viljanen and two anonymous reviewers for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1 Linear regression of the annual percentage rate of EGFs from groups DEF in January 2010–July 2014 (circles, left axis) and yearly changes of sunspot number (orange squares, right axis).
In the text

	Fig. 2 Geoelectric field components EX and EY for years 2010–2014 calculated using a 1-D layered-Earth model (Ádam et al., 2012).
In the text

	Fig. 3 % above normalised sum of daily electrical grid failures (a), averaged over the number of events (b) Dst-index [nT], (c) BZ [nT] component of HMF and (d) computed geoelectric field E [mV/km], during the geomagnetic storms with BZ < −10 nT (case 1), for years 2010–2014.
In the text

	Fig. 4 % above normalised sum of daily electrical grid failures (a), averaged over the number of events (b) Dst-index [nT], (c) BZ [nT] component of HMF and (d) computed geoelectric field E [mV/km], during the geomagnetic storms with Dst < −100 nT (case 2), for years 2010–2014.
In the text

	Fig. 5 % above normalised sum of daily electrical grid failures (a), averaged over the number of events (b) Dst-index[nT], (c) BZ [nT] component of HMF and (d) computed geoelectric field E [mV/km], during the intense geomagnetic storms with Dst < −100 nT and BZ < −10 nT (case 3), for years 2010–2014.
In the text

	Fig. 6 SOM neighbor weight distances with weight values of connections between neighboring neurons for data without delay for the intense geomagnetic storm, 15 July 2012. The blue hexagons represent neurons, and the red lines show which particular neurons are connected. Colors from black to yellow display the weight values of the connection between neighboring neurons.
In the text