Evaluating the relationship between strong geomagnetic storms and electric grid failures in Poland using the geoelectric field as a GIC proxy

We study intense geomagnetic storms (Dst < 100nT) during the first half of the solar cycle 24. This type of storm appeared only a few times, mostly associated with southwardly directed heliospheric magnetic field Bz. Using various methods such as self-organizing maps, statistical and superposed epoch analysis, we show that during and right after intense geomagnetic storms, there is growth in the number of transmission line failures. We also examine the temporal changes in the number of failures during 20102014 and find that the growing linear tendency of electrical grid failure occurrence is possibly connected with solar activity. We compare these results with the geoelectric field calculated for the region of Poland using a 1-D layered conductivity Earth model.


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(Tozzi et al., , 2019Piersanti 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 Topical Issue -Geomagnetic Storms and Substorms: a Geomagnetically Induced Current perspective 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 partcycle semisaturation (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.

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

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 A. Gil et al.: J. Space Weather Space Clim. 2021, 11, 30 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.

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 0 : m 1 = m 2 , where m 1 and m 2 express the means of the number of DEF failures of the first and second group, respectively, whereas H 1 : m 1 6 ¼ m 2 is the alternative hypothesis. From the t-test distribution table, one takes a critical value t a,n for a given level of significance a and n degrees of freedom. The critical set for the t-test is (À1, Àt a,n ] [ [t a,n , 1). If t < |t a,n | we do not reject hypothesis H 0 , otherwise we accept H 1 . The value of p, given in the t-test expresses the probability of error involved in accepting hypothesis H 1 about the existence of a difference. A. Gil et al.: J. Space Weather Space Clim. 2021, 11, 30 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).

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 & 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 r 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;: 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 l 0 = 4p Â 10 À7 H/m . The initial value in equation (1) is the case when the layer n = N is a uniform half-space and K N ¼ , while the final value K 1 is the transfer function relating E and B at the Earth's surface (Trichtchenko & Boteler, 2002;. 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/r): 6 km, 5 Xm; 105 km, 1000 Xm; 300 km, 100 Xm; above a half space of 10 Xm. This Earth-layered model for the European regions has already been used in a number of studies (e.g., Viljanen et al., 2013Viljanen et al., , 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, Finally, we took the inverse Fourier transform to obtain the geoelectric field E(t) in the time domain (Boteler, 2012;. 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;. 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.

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.

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: where w ij is the weight j of neuron i, a(t) is the learning rate, h ci (t) is the neighborhood function (Caldas et al., 2017) that A. Gil et al.: J. Space Weather Space Clim. 2021, 11, 30 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).

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.
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 0 hypothesis. In our calculations, we put a = 0.05, for a < p we accept H 0 hypothesis, otherwise, when a > p we reject H 0 and accept H 1 which is the alternative hypothesis, moreover, when a < p 0.1 we can also reject H 0 and accept H 1 . 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). 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).

Computed geoelectric field
We are aware that the induced E field can vary widely even during a given event 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. A. Gil et al.: J. Space Weather Space Clim. 2021, 11, 30

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 Figure 3 displays case 1 circumstances when B Z < À10 nT for more than 3 h. Figure 3, panel (a) presents that on the zeroday 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 .31 p = 0.04 p = 0.23 * In the case of the assumption of normality of our distributions, we formulate hypothesis H 0 : the distribution of data is normal for mean 24 h+ as for mean 24 hÀ. Using the Shapiro-Wilk test (for a small sample) we observe that the assumption on normality is not fulfilled for mean 24 hÀ, because p = 0.04, < 0.05. Thus, we reject H 0 and accept H 1 , which says that the distribution of data is not normal. Analogically, we check assumptions for the Levene test. (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. 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.

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] 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., 2019Gil et al., , 2020aŠ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  A. Gil et al.: J. Space Weather Space Clim. 2021, 11, 30 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.
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.  A. Gil et al.: J. Space Weather Space Clim. 2021, 11, 30 However, the greater distance from the neurons related to electrical grid failures expresses their lower impact.

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 networktherefore, more investigations are needed to determine the effects of GIC in technical infrastructure.

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.