Analysis of geomagnetically induced currents at a low-latitude region over the solar cycles 23 and 24 : comparison between measurements and calculations

Geomagnetically Induced Currents (GIC) are a space weather effect, which affects ground-based technological structures at all latitudes on the Earth’s surface. GIC occurrence and amplitudes have been monitored in power grids located at high and middle latitudes since 1970s and 1980s, respectively. This monitoring provides information about the GIC intensity and the frequency of occurrence during geomagnetic storms. In this paper, we investigate GIC occurrence in a power network at low latitudes (in the central Brazilian region) during the solar cycles 23 and 24. Calculated and measured GIC data are compared for the most intense geomagnetic storms (i.e. 50 < Dst < 50 nT) of the solar cycle 24. The results obtained from this comparison show a good agreement. The success of the model employed for the calculation of GIC leads to the possibility of determining GIC for events during the solar cycle 23 as well. Calculated GIC in one transformer reached ca. 30 A during the ‘‘Halloween storm’’ in 2003 whilst most frequent intensities lie below 10 A. The normalized inverse cumulative frequency for GIC data was calculated for the solar cycle 23 in order to perform a statistical analysis. It was found that a q-exponential Tsallis distribution fits the calculated GIC frequency distribution for more than 99% of the data. This analysis provides an overview of the long-term GIC monitoring at low latitudes and suggests new insight into critical phenomena involved in the GIC generation.


Introduction
Geomagnetically Induced Currents (GIC) are ground effects due to space weather phenomena.The GIC have widely been monitored and studied in ground-based technological structures located at high latitudes, after the Carrington event that occurred on 1-2 September 1859 (Tsurutani et al. 2003).Another significant GIC event occurred in March 1989, when the entire Canadian Hydro-Québec power system was disrupted for several hours during a geomagnetic storm (Bolduc 2002).Despite the occurrence of intense and extreme space weather events at low latitudes during the last decades, GIC effects at mid-to-low latitudes have remained less explored.
Measurements of GIC magnitudes during different geomagnetic storm phases show that values obtained in power networks located at low and middle latitudes can reach the same levels as those observed at high latitudes.For example, Trivedi et al. (2007) present GIC recordings around 15 A in a Brazilian power network during 2004, which were also recently reproduced by calculations (Barbosa et al. 2015).These GIC events lasted for around 5 hours on 7-9 November 2004.Watari et al. (2009) observe about 15 GIC events in Japan, from 2007 to 2009.Marshall et al. (2012) report about a failure in a transformer located in New Zealand during the sudden impulse on 6 November 2001.Liu et al. (2009) present GIC values of 47.2 A and 75.5 A in China associated with the geomagnetic storms on 7-9 November 2004, respectively.Moreover, Gaunt & Coetzee (2007) and Kappenman (2005) verify a series of transformer damages ascribed to high intensity GIC in South Africa during the ''Halloween storm '' on 29-31 October 2003.Consequently, meaningful GIC magnitudes can be observed at both high and low-latitude regions.Ground-based technological infrastructures, including especially power networks, are exposed to GIC effects regardless of their location.Additionally, in countries with extensive territories (compared to typical continental dimensions), the power network lengths are becoming larger, which contributes to enhancing GIC values (Zheng et al. 2014).Finally, it is worth stressing that the vulnerability of a power system due to the GIC depends on many technical details including the transformer types used.Thus, certain GIC magnitudes could be harmless in one power network whilst catastrophic in another.
In this work, we present an analysis of the GIC data measured in the neutral of the 500 kV transformer at the Itumbiara substation in central Brazil during the solar cycle 24 (Fig. 1).The GIC data recorded during geomagnetic storms (À150 < Dst < À50 nT) are compared to GIC values obtained from calculations.We use the comparison between these values to validate the calculation method, including the power network model.Once the validation succeeds, the GIC are calculated for the same power grid considering geomagnetic storms (Dst < À100) during the solar cycle 23.Finally, a statistical analysis is performed by fitting a Tsallis distribution (Tsallis 1988) to the normalized inverse cumulative frequency (NICF) for GIC data.Implications of such distribution are discussed.

Measurements of GIC
The GIC data have been obtained in the neutral of the Itumbiara 500 kV transformer.The whole power network located in the central territory of Brazil is shown in Figure 1.There are no series capacitors in the transmission lines, and the ground is highly resistive in that region.The topology of the network was different in 2003 (during the solar cycle 23) from that in 2014 (solar cycle 24), since new transmission lines were installed in 2006 (Fig. 1).The details of the network before 2006 are reported by Trivedi et al. (2007) and Barbosa et al. (2015).
In order to provide GIC data with a high signal-to-noise ratio (SNR), a sensor for measuring GIC based on the Hall effect was mounted at the neutral of the Itumbiara transformer.The sensor package was specifically designed for a correct mechanical positioning around the neutral lead; the Hall sensors were in pairs allowing the sensors to be 1 mm apart from the lead.The sampling rate of the GIC data was six samples per minute, and the measurements were continuous from 2009 until 2013, except for certain periods of maintenance.Temperature variations were monitored daily.
GIC data were obtained from 32 geomagnetic storms (Table 1) identified by using the Dst definition of geomagnetic storms (Gonzalez et al. 1994).In this paper, only the three storms that produced GIC amplitudes higher than 4 A are analysed.They are highlighted in bold face in Table 1.

Calculation of GIC
Calculation methods of the GIC in a power network are based on Kirchhoff's laws.Two main approaches can be used: the mesh impedance matrix method and the nodal admittance matrix method (Pirjola 2007;Boteler 2014;Boteler & Pirjola 2014).The former is constructed representing the network as a mesh of loops and considering the voltage sources in series with the resistances in the transmission lines.In the second method, the network is regarded as a set of nodes grounded by earthing admittances.In this method, the equivalent current sources are assumed to be located in parallel with the admittances of the transmission lines connecting the nodes.The Lehtinen-Pirjola (LP) method, developed independently of the load-flow computation techniques, by Lehtinen & Pirjola (1985), is a variant of the nodal admittance method.Since the mesh impedance matrix method leads to long computation times, the nodal admittance matrix method and the equivalent LP method are preferred in practical GIC computations.In this paper, we choose the LP method due to its wide application, especially in the geophysical community.
The use of the LP method presumes that the geoelectric field affecting the power network is known.Therefore, the GIC computation contains the following two steps (e.g.Pirjola 2002).First, we determine the geoelectric field associated with the geomagnetic variation.Second, we apply the LP method with the geoelectric field as input.

Calculation of the geoelectric field
The time variations of the geomagnetic field and so the geoelectric field (GEOF) at the Earth's surface are primarily affected by currents in the magnetosphere and the ionosphere and secondarily by currents and charges induced in the ground.The electromagnetic fields are determined by solving Maxwell's equations and applying appropriate boundary conditions.Different techniques have been proposed to determine the geoelectric fields, e.g. the complex image method (CIM) (Boteler & Pirjola 1998;Pirjola & Viljanen 1998) and the plane wave method in conjunction with geomagnetic field measurements (Cagniard 1953).The results obtained using the second technique have shown to be in accordance with measurements at different latitudes.Therefore, in the present analysis, the GEOF calculation is based on the plane wave method.For that, we consider a small region of the Earth's surface, which defines the xy plane of a coordinate system, with the x, y and z axes being northward, eastward and downward, respectively.The geomagnetic disturbance produced by ionospheric-magnetospheric currents is represented by an electromagnetic plane wave, which propagates vertically downwards.The Earth's conductivity structure is laterally uniform (i.e.layered).The use of the plane wave method thus strictly presumes that there are no lateral changes in the Earth's conductivity or in the geoelectric and geomagnetic fields over the area considered (e.g. the area of a power network).In practice, however, this is often not the case.Then the region is divided into blocks, with different ground conductivity structures, and the plane wave technique is applied separately in each block (Viljanen et al. 2004).This approach is called the ''local plane wave method'' (or a ''piecewise layered-Earth model'' (Marti et al. 2014).
In the plane wave method, the time series corresponding to the geomagnetic field data are Fourier transformed to the frequency domain, and then the geoelectric northward (E x ) and eastward (E y ) components are obtained from the relation where Z(x) is the surface impedance, which is a function of the angular frequency and the Earth's conductivity structure.B y and B x are the measured geomagnetic field eastward and northward components, respectively.
In this work, we use a ground conductivity model obtained from magnetotelluric surveys, as described by Bologna et al. (2001).It consists of four layers, with the deepest one being infinitely thick.Such a model represents the local geology in the area of the 500 kV power network shown in Figure 1.
The magnetic observatory closest to the Itumbiara substation is located at Vassouras, Rio de Janeiro State (VSS, 22.4°S, 43.6°W, see Fig. 1).The distance between VSS and Itumbiara is around 740 km.The large distance between the area considered and the VSS observatory can introduce errors in the GEOF estimation.We evaluate the inaccuracy of the GEOF calculation by comparing the GEOF results obtained from VSS geomagnetic field recordings to those obtained using interpolated geomagnetic field data at the closest point to Itumbiara (Fig. 2).For the interpolation, we employ the ''Spherical Elementary Current Systems (SECS)'' method (Amm 1997;Amm & Viljanen 1999;Caraballo et al. 2013) using data from four geomagnetic stations belonging to the EMBRACE/INPE network (Denardini et al. 2013).The operating times and locations of all EMBRACE/INPE stations used in the interpolation are shown in Table 2.
We find a root mean square (rms) value of 0.0465 V/km for the absolute values of GEOF calculated via VSS data and via SECS-interpolated data at Itumbiara on 8 October 2013.In addition, for other events during the solar cycle 24, the rms values between GEOF calculated by the two methods are lower than 0.050 V/km.These small rms values in GEOF data indicate that, despite the great distance between VSS and the power network area, VSS geomagnetic field data can be used to calculate GEOF at the power network at all times, The geovoltage between two points i and j at the Earth's surface is calculated by integrating GEOF along a trajectory from i to j, which, in the case of a power network, is the transmission line between i and j: The power grid is modelled as a network of N earthed substations (i.e.nodes) interconnected by transmission lines with voltage sources in series with the line resistances.In the LP method, the GIC flowing into the ground at the nodes are expressed as a N • 1 column matrix I e (Lehtinen & Pirjola 1985): where Id is the N • N identity matrix and Y n is the N • N network admittance matrix with their elements given by: where R ij is the resistance of the transmission line between the nodes i and j (being infinite if i and j are not connected by a transmission line), and Z e is the N • N earthing impedance matrix with the diagonal elements being the earthing resistances (r i ) of the nodes: If the earthing points are far apart (as in the case of different substations), the off-diagonal elements of Z e are set equal to zero.The elements of the N • 1 column matrix J e are defined by:  J e thus includes the information about the geovoltages V ij .Since I e = J e in the case of perfect earthings of the nodes (i.e.Z e = 0), the elements of J e can be interpreted as ''perfect earthing currents''.

Analysis method of the GIC data and
Tsallis statistics A combination of the method to determine the GEOF discussed in Section 2.3 with the LP method described in Section 2.4 is used to calculate the GIC level in the power grid shown in Figure 1.First, we use some of the most intense geomagnetic storms during the solar cycle 24 to validate this modelling through the comparison between calculated and measured GIC data at the Itumbiara substation.Second, we use the model to estimate the GIC values reached during major storms of the solar cycle 23.The geomagnetic storms registered during the solar cycle 23 are more intense than those during the solar cycle 24.
A comparison between GIC intensities obtained during disturbed periods of the solar cycles 23 and 24 provides an estimate of the impact of great geomagnetic storms on the Itumbiara power network.The GIC intensities calculated for the solar cycle 23 provide a large number of data.Since they encompass the complete solar cycle, we can determine the frequency distribution of GIC intensities during the solar cycle 23.GIC are calculated for three days during geomagnetic storms (i.e. for the day of the Dst minimum and for the day after and the day before).
Statistical distributions have previously been proposed to model the frequency distribution of the occurrence of the geoelectric field, e.g. a power law (Langlois et al. 1996;Watari et al. 2009;Myllys et al. 2014) and a log-normal distribution (Pulkkinen et al. 2008).In these previous works, the statistical distributions proposed succeed in fitting only a subset of the geoelectric field dataset, as the tail of the distribution (e.g.Myllys et al. 2014) or the middle values (e.g.Watari et al. 2009).All the proposed distributions fail when they are applied to the complete dataset.
For the purpose of fitting the complete dataset, we apply a generalization of the power law distribution, namely a Tsallis distribution, to the GIC calculated for the Itumbiara substation.The criteria to choose a Tsallis distribution function to represent the GIC frequency distribution are based on two physical aspects.Firstly, to calculate the induced GEOF (and, thus GIC) in the time domain, a time-retarded geomagnetic field can be used.This considers a memory effect in the system.Besides, the ground effects that come from space weather events are also a result of long-range interaction with magnetospheric currents.Secondly, the GEOF produced by induction is driven by an out-of-equilibrium system (i.e.magnetospheric and ionospheric currents perturbed during geomagnetic storms) (Balasis et al. 2009).Typical systems following Tsallis distributions are described as being originated by complex and out-of-equilibrium processes.They show long-range spatial connections between their elements including a memory effect.Besides, power law distributions are special cases of Tsallis distributions in the limit of extreme distribution values (Tsallis 2009).
In order to investigate the goodness of the fitting in detail, it has been proposed that the cumulative frequency of the calculated GIC at Itumbiara, during the solar cycle 23, follows a Tsallis distribution.The Tsallis distribution is a generalization of the Boltzmann-Gibbs (B-G) distribution developed by Tsallis (1988), to represent physical systems out of equilibrium.The generalized Tsallis entropy S q is defined by: where k is a suitable scaling constant, the p i quantities are the probabilities and W is the total number of states.The Tsallis statistics is non-extensive (i.e. the entropy is not proportional to the number of states in a system), and the entropic parameter q in equation ( 7) can be understood as a measure of the non-extensivity.For the limit q = 1, S q tends towards the B-G entropy.The Tsallis frequency distribution p(x) can be obtained by maximizing the Tsallis entropy given by equation ( 7): where B is a scale parameter and C is a normalization constant.The parameter q has a physical meaning only when it reaches the interval 1 q 3. The q parameter can be interpreted as a measurement of the correlation level in the physical system that produces the GIC effect.When q ! 1, the system can be purely uncorrelated whilst as q departs from unity, the system becomes more correlated.If the variable is constrained to the positive domain in the maximum entropy procedure, the q-exponential distribution (Eq.( 8)) is derived.The main features of systems that follow Tsallis distributions are spatial long-range interaction, temporal memory effect, absence of characteristic size and multifractal nature.

Comparison between measured and calculated GIC in the
Itumbiara 500 kV transformer GIC measurements are performed in the neutral of the Itumbiara 500 kV transformer during the solar cycle 24, as described in Section 2.1.We choose three geomagnetic storm events with high GIC intensities for the purpose of comparing the measured and calculated GIC.These results are shown in Figure 3, panels (a) 16 June 2012, (b) 2 October 2013 and (c) 8 October 2013, where the blue and red curves depict the calculated and measured GIC, respectively.We note that the calculated and measured GIC peaks are in good agreement.
Regarding the event on 2 October 2013 (Fig. 3b), the curves are also similar when the GIC values are smaller, so a good agreement between calculations and measurements holds true for most of the GIC dataset.The small discrepancies at lower GIC amplitudes can be related to local ionospheric instabilities, which occur in the vicinity of equatorial latitudes (see Abdu 2005), where Itumbiara and Vassouras are located.Thus, comparing calculated and measured GIC data indicates that the C. Barbosa et al.: Analysis of geomagnetically induced currents at a low-latitude calculation procedure is reliable and so it can be applied to the determination of GIC at Itumbiara during periods when measurements are not available.

Frequency distribution of GIC at the Itumbiara 500 kV transformer for the solar cycle 23
The GIC values at the Itumbiara 500 kV transformer for large magnetic storms (Dst < À100 nT, 90 events in total) that occurred during the solar cycle 23 (Pandey & Dubey 2009) are calculated for the day when the Dst reached its minimum value, and for the previous day and the day after.
Figure 4 shows the minimum Dst value used as a measure of the storm's intensity and the daily maximum values of calculated GIC amplitude at Itumbiara during the solar cycle 23.The Dst peak and maximum GIC amplitudes need not coincide in time.We observe that for the majority of the storms the Dst peak was À100 nT < Dst < À200 nT and the GIC values were between 5 and 15 A. However, the GIC measured at Itumbiara substation reached about 30 A during the ''Halloween storm'' in 2003 and the corresponding Dst peak was À401 nT.
Figure 4 also provides the GIC level reached during geomagnetic storms with different intensities.However, it is worth mentioning that the GIC level reached during a geomagnetic storm also has contributions from other magnetospheric sources.Those sources may include magnetopause currents (related to sudden impulses), pulsations (due to magnetospheric plasma instabilities), changes in ionospheric solar quiet currents (observed during ionospheric storms and substorms) and a disturbed dynamo.Some of these phenomena are not properly described by the Dst index.In the data presented in Figure 4 we can find an example of these contributions.There are two storms represented by Dst < À400 nT and each of them is related to a different GIC value, namely 30 A and  15 A. This substantial difference may be due to the mechanism of electromagnetic induction.Since GIC are a result of regional ionospheric current system and the Dst index is a global measure of the geomagnetic field disturbance, the local ground magnetic field perturbation will be different for storms entering the main phase at different local times.
With the calculation of GIC for a long period (including about 90 events), it is possible to propose a probability distribution function to fit the GIC occurrence frequency.At least two relevant pieces of information can be obtained with the knowledge of the GIC frequency distribution.First, the obtained function can be used to make a rough estimate of the GIC level for extreme geomagnetic storms and determine its occurrence frequency (Thomson et al. 2011).In the second place, the knowledge of the probability density function of a physical quantity can give us insight into the nature of the processes that generate it.Pulkkinen et al. (2008) calculated the geoelectric field at sites in northern Europe for the period from January 1993 to December 2006.They fitted the GEOF frequency distribution to a lognormal distribution.However, by using the Shapiro-Wilk test, they found that the lognormal distribution failed in representing the dataset.Myllys et al. (2014) calculated the GEOF at sites in Norway between 1994 and 2011 and realized that for extreme GEOF values (i.e. the tail of the distribution) a power law distribution fit succeeds.Nevertheless, for small and intermediate values of GEOF a discrepancy between the fit and the data is observed.As the power law distribution is a particular case of Tsallis distributions, we tested if a class of Tsallis distributions succeeds in fitting the entire range.
Figure 5 shows the frequency distribution of GIC calculated for each geomagnetic storm during the solar cycle 23.The GIC values are in the 1 A to 30 A range.The limits of the GIC magnitudes applied to the statistical investigation are defined by the characteristics of the analysed sample.The lower limit represents minimum GIC values observed when the geomagnetic field was active.On the other hand, the upper limit for the GIC statistical study was chosen as 10 A. GIC data comprise 270,658 values; from this total only 219 exceed 10 A, corresponding to less than 0.1% of the total GIC dataset.This implies that ~10 À3 above the interval diminish the normalized inverse cumulative frequency (NICF).This value corresponds to the GIC intensity for which the frequency data is diminished by three orders of magnitude with respect to the complete dataset.The NICF distribution is defined as the number of GIC values higher than a given GIC bin.The data are binned in intervals of 0.1 A. Since around 99.5% of the calculated GIC amplitudes are lower than 10 A, we consider this value as the limit of the statistically representative dataset for the NICF.For these conditions, we verify that the best fit to the NICF is the q-exponential Tsallis distribution given by equation ( 8) with parameters q = 1.25, B = 1.77A À1 and C = 1.18, since this minimizes the v 2 quantity, and gives a correlation coefficient of R 2 = 0.999.These results suggest that GIC are physical phenomena compatible with extensive statistical mechanics.In other words, the physical mechanisms responsible for generating GIC might belong to a nonequilibrium system.As a consequence, small variations on the geomagnetic field such as those observed during magnetic storms (i.e.small perturbations in the solar-terrestrial system) can drive extreme GIC events, suggesting that GIC phenomena could be a part of the non-extensive statistical mechanics domain.

Conclusions
In this work, we present the analysis of the GIC measurements performed during 2009-2013 (the solar cycle 24) in the neutral of a high-voltage transformer at the Itumbiara substation located at low latitudes in central Brazil.In this 5-year period, three geomagnetic storms caused GIC amplitudes higher than 4 A. For those events, we calculate GIC by using the Lehtinen-Pirjola method.The geoelectric field is determined by two methods.In the first, the geomagnetic field is considered uniform and equal to the recorded field data at the Vassouras magnetic observatory located at a distance of ~700 km from Itumbiara.In the second, the geomagnetic field is interpolated at Itumbiara substation using the ''Spherical Elementary Current Systems'' method.The results show that the calculated and measured GIC values are in good agreement.Thus, the calculation method is accurate enough to be applied at lowlatitude regions in Brazil using geomagnetic field data from Vassouras.
Consequently, we can obtain the GIC data at Itumbiara (and in the whole high-voltage power network) from calculations for the solar cycle 23 when the GIC recordings are not operational.We note that the geomagnetic storms observed during the solar cycle 23 were more intense than those during the solar cycle 24.The calculations show that GIC may have reached values as high as 30 A during the most intense event, i.e. during the ''Halloween storm '' on 29-31 October 2003.Based on the calculated GIC data, a statistical distribution is fitted to the NICF of GIC occurrence at Itumbiara.It was found that a q-exponential class of the Tsallis statistical function fits more than 99% of the data in the range from 1 A to 10 A with a correlation coefficient of 0.999.This result has implications on the understanding of the physical mechanisms resulting in GIC generation.These implications suggest GIC are caused by magnetospheric current systems in a nonequilibrium state having a spatial long-range correlation and a memory effect.In this context, a non-equilibrium state means that a small perturbation in the system could drive a significant GIC effect.Spatial long-range correlation indicates the strong coupling including elements that compose the Sun-Earth system and their effects at the Earth's surface.Finally, memory effect means that the GIC value at each instant is correlated to its past values and depends on the past values of geomagnetic variations.The Tsallis distribution that fits the Itumbiara GIC data needs to be extended to other datasets from different locations.In this way, we will be able to confirm if our observations correspond to the global understanding of the physical processes involved in GIC phenomena regardless of the location and the power network configuration.

Fig. 1 .
Fig. 1.The 500 kV power network investigated in this paper and located in central Brazil.The transmission lines shown in blue were installed in 2006 and those in red earlier.The locations of the VSS magnetic observatory and of the nearest EMBRACE/INPE magnetic stations (SJC, CXP) used to interpolate geomagnetic data are presented by green marks.

Fig. 2 .
Fig. 2. Geoelectric field modulus calculated using the VSS data (red curve) and by using the SECS-interpolated data (black curve) at Itumbiara on 8 October 2013 (upper panel) and their difference (lower panel).The rms value equals 0.0465 V km À1 .

Fig. 3 .
Fig. 3. Calculated (blue) and measured (red) GIC in the neutral of the Itumbiara 500 kV transformer on (a) 16 June 2012, (b) 2 October 2013 and (c) 8 October 2013.The right panel shows a zoom at the intervals where peak of GIC occur.

Fig. 4 .
Fig. 4. Daily maximum amplitude of calculated GIC in the neutral of the Itumbiara 500 kV transformer versus the intensity of the magnetic storm given by the Dst index during the solar cycle 23.

Fig. 5 .
Fig. 5. Normalized inverse cumulative frequency (NICF) distribution of calculated GIC in the neutral of the Itumbiara 500 kV transformer for strong geomagnetic storms during the solar cycle 23.The black dots show the data.The red curve represents the q-exponential Tsallis distribution fit given by Eq. (8) with q = 1.25, B = 1.77A À1 and C = 1.18.

Table 1 .
Solar cycle 24 geomagnetic storms during which GIC was measured in the neutral of the Itumbiara 500 kV transformer (see Figure1).The three storms that produced GIC amplitudes above 4 A are shown in bold face.

Table 2 .
Location of the EMBRACE/INPE magnetic stations and their periods of operation during the solar cycle 24.The negative sign denotes the southward and westward latitude and longitude, respectively.Additionally, the Itumbiara substation coordinates where the GIC are recorded, are shown.