Geomagnetically Induced Currents: Frequency spectra and threats to voltage stability

Geomagnetically Induced Currents (GICs) are induced after a complex interaction between the sun and earth’s magnetic field. GICs may cripple power systems leading to voltage collapse and a drop in power quality. While the commonly-used dc model for GICs is valid for steady-state analyses of power systems and transformer responses, GICs must be represented with their continuously varying magnitudes and frequency components to understand the dynamic power system’s response to solar storms. In this paper, we used measured GIC data from the 2003 Halloween storm to compare several signal processing techniques. The Fast Fourier Transform (FFT) with a fixed window size, Short-Time Fourier Transform (STFT) with a variable window size, wavelet transform, and Particle-Swarm Optimization (PSO) methods were used to identify the GIC frequency characteristics. To understand the effects of these GIC characteristics on voltage stability, a multi-machine network was modelled on a digital real-time simulator (OPAL-RT). The results reveal the limitations of fixed window FFT in capturing the dominant high-frequency GIC bandwidth arising during peak solar activity. Wavelet and STFT are recommended for GIC data analysis as they provide time-frequency localization. High-frequency GICs during sudden storm commencements and pulsation activity are more of a concern for voltage stability than lower frequency components. Novel results from this study suggest that limiting contingency analysis to only high magnitudes of dc/GIC does not cater for the underlining GIC dynamics. In addition to magnitude thresholds, contingency analysis should incorporate GIC characteristics such as frequency which affect voltage dynamics and thus power system stability.


I. INTRODUCTION
Geomagnetically Induced Currents (GICs) are initiated by increased solar activity. When active regions of the sun (often associated with sunspots) become unstable and magnetic reconnection occurs, large amounts of solar plasma can be ejected from the surface of the sun. These impulsive solar storms or coronal mass ejections (CMEs) can drive high magnitude GICs whicch may cause damage and disruptions to grounded electric power systems. The solar wind carries these ejected charged particles throughout the heliosphere [1], interacting with the Earth's magnetic field along way. The coupling between the solar wind and the Earth's magnetosphere-ionosphere system gives rise to sudden and rapid variations in the geomagnetic field (B-field) on the Earth's surface. This phenomenon is called a geomagnetic disturbance (GMD) [1], [2].
The ground geomagnetic field response associated with a GMD is highly dependent on the coupling between the near-Earth current systems and the solar wind. The most effective GMDs are typically a result of the interplanetary magnetic field opposing the Earth's magnetic field, allowing solar wind plasma to couple directly to near-Earth current systems. In these cases, the GMD is associated with a geomagnetic storm. The global morphology of a geomagnetic storm can be described as an initial intensification of the geomagnetic field, known as the storm sudden commencement (SSC). Following an SSC is the main phase of the geomagnetic storm. During the main phase, the geomagnetic field is highly disturbed and the main component generally weakens, driven by strengthening near-Earth current systems. The period following the peak of a geomagnetic storm, as the system returns to equilibrium and geomagnetic activity returns to quiet-time levels, is known as the recovery phase of a geomagnetic storm. According to Maxwell-Faraday's law of electromagnetism (∇ × = − / ), the varying geomagnetic field, , induces an electric field (E-field), , on the ground surface [1], [2]. The E-field then drives GICs into power systems through the neutrals of grounded wye-connected transformers. Transformers are designed to operate in the linear region of their B-H characteristics, drawing very small magnetizing currents from the network. GICs shift the operating point of transformers into the saturation region leading to an increase in the magnetizing current, generation of even and odd harmonics, transformer heating and gassing [3], [4], [5], [6], [7]. Core and winding losses in the transformer increase due to harmonic current flow. Noise levels may increase due to the asymmetric magnetization of the core during GIC flow [8], [9]. Harmonics including those of non-conventional sequences may also cause rotor heating in synchronous generators [10], [11]. The flow of harmonic currents increases system losses as these harmonic currents experience higher impedances along three-phase transmission line conductors due to the 'skin effect'.
The presence of even and triplen harmonics due to saturating transformers under GMD conditions causes the misoperation of relays and improper operation of breakers controlling SVCs and capacitor banks. As a result, a threat is posed to the voltage stability of power systems [12], [13], [14], [15], [16]. The cumulative effects of transformer saturation and harmonic current flow result in reduced efficiency in the transfer of useful (real) power to the loads.
Significant progress has been made over the years in understanding the geophysical processes leading to GICs. In power system analysis, GICs have traditionally been modelled as dc currents. While this well-established dc approximation can be valid for certain studies (such as effects of GIC on transformer heating, static voltage stability analysis using PV and QV curveswhich all come with the unrealistic assumptions of balanced conditions in the phases at a steady state), it cannot be used to study the actual voltage dynamics that arise due to the underlying GIC dynamics. Therefore, there is a need to either use time-domain measured GIC data or more representative frequencydependent GIC models for dynamic voltage stability analysis. In the absence of measured GIC data, lowfrequency ac GIC models may be used for GIC modelling instead of the dc approximation [17], [18], [19], [20]. In such cases, the GIC magnitudes and frequencies during different phases of the storm need to be identified to accurately replicate their effects on the power system.
In this paper, we use measured GIC data from the 2003 Halloween storm to compare four methods to extract the frequency components of GICs revealing the limitations and advantages of the different approaches. These methods are: Fast Fourier Transform (FFT), Short-Time Fourier Transform (STFT), the wavelet transform and Particle Swarm Optimisation (PSO). We explore the following questions to establish a preliminary framework that considers the frequency dependence in the GIC waveform and its implications on voltage stability: 1. Which techniques are suitable to identify the frequency components of GIC? 2. What are the typical GIC frequencies and their associated rates of voltage decline? 3. Which GIC driving regime is more of a concern for voltage stability of power systems? The novel results from this study, therefore, serve as a preliminary understanding of voltage stability threats to power networks during GMDs.

II. LITERATURE REVIEW
In this section, we review the approaches used in previous research to identify the frequency components of GICs and discuss papers that investigated the voltage stability of power systems under GIC conditions.

A. FREQUENCY COMPONENTS OF GICS
In the frequency domain, the B-field and E-field are related via the frequency-dependent surface impedance of the Earth [1], [2], [21], [22]. This surface impedance tensor can be a 1D or 3D model made up of many geophysical variables characterised by continuous frequency spectra. Since the Efield drives the GIC, the latter is also characterised by a spectrum of frequencies. The geophysical processes may drive large GICs, especially during storm sudden commencements (SSCs), geomagnetic pulsations and auroral substorms [23]. This suggests that modelling GICs as dc does not cater for the temporal variations in realistic GIC characteristics (magnitude and frequency).
Trichtchenko [24] highlights the need for high sampling rates in GIC and geomagnetic field data measurements to capture the variability in GIC. A 1 s data cadence is recommended to capture the peaks and short term events of 10-15 s duration. Watari et al. recommend a higher sampling to pick up shorter-term dynamics [25]. Harmonic analysis of the transformer neutral current revealed that there are differences between the GIC and harmonic variation [24]. Delays between 25-30 s were observed in the harmonic response. Such delays are introduced due to the non-linear transformer inductance [26].
Pulkkinen and Kataoka [27] provided a time-frequency analysis of GICs induced during superstorms using the Stransform -an extended STFT. One distinct feature of all Stransform plots showed a "flame-like or spotty" appearance originating from the noise-like character of the GIC signal. They show that the spectral content of GICs depends on the storm phase. The main phase of the storm is characterised by wide-band fluctuations with frequencies between 0.1 mHz to 4.2 mHz [27]. Geomagnetic pulsation activity characterized by narrow-band fluctuations may be observed during the recovery phase of the storm [27]. Even in the presence of noise within the GIC data, the S-transform can bring out the storm phase-dependency in characterising GICs. This is important to understand the threats to voltage stability of power systems under different storm phases and to plan contingency measures for such events.
Mandrikova et al. [28] investigated the periodicities in magnetic field data using the wavelet transform. Analysis of the data over different time intervals confirmed that frequency components of magnetic field variations during solar storms are irregularly distributed in time. This finding suggests that GICs are non-stationary currents whose frequency components do not exist indefinitely but instead vary with time. FFT which is based on stationary signals is, therefore, an approximation of the frequency components of GICs for a defined interval. Based on this, the authors show the effectiveness of the wavelet transform in identifying periods of increased geomagnetic activity, durations of such events, and the temporal variations in GIC characteristics.
Zaourar et al. [29] used a wavelet-based multiscale analysis to obtain quantitative information on the short-term dynamics of B-fields. B-field spectral contents from three storms including the 2003 Halloween storm were analyzed for their variation in a scaling exponent. This scaling exponent indicates the intensity and properties of geomagnetic disturbances especially the temporal evolution of storm activity. Results showed that the B-field spectra are continuous and broadband. Therefore, the E-field spectra and GIC spectra exhibit similar non-periodic trends. As a result, FFT analysis of GIC data using a fixed window for the full duration of the storm may provide misleading results because FFT assumes periodic signals of infinite duration. In the absence of high frequency GIC data measurements, FFT might also fail to capture the high frequency components associated with peak GICs.
Falayi et al. [30] used statistical methods to calculate the correlation between measured GIC and ⁄ . To extract the wavelet energy coefficients and periodicities from the measured GIC, wavelet analysis was used with the Morlet wavelet as the mother wavelet. It was observed that high energy wavelet coefficients of ⁄ correlated well with GIC. Similar to Falayi et al., another study applies Continuous Wavelet Transform (CWT) and Discrete Wavelet Transform (DWT) to GIC and ⁄ showing very similar spectrograms [31]. The wavelet coefficient of ⁄ and GIC are very similar indicating strong correlation between them. These results depend the directionality of the network with respect to the induced Efield and the surface impedance of Earth modelled either in 1D or 3D. No continuous periodicities in the GIC spectra were observed.
The rate of change of the magnetic field, / , is proportional to ( ) in the frequency domain, where represents frequency. As a result, dB/dt has a flatter power spectrum compared to the B-field [21]. The GIC frequency spectrum lies between the two, i.e. ⁄ has a high frequency bias relative to GICs and the B-field has a low frequency bias relative to GICs. Since impulsive peaks are associated with high frequencies, ⁄ may track such peak events better than the B-field. While ⁄ may provide a preliminary understanding of the expected GIC magnitudes, Heyns et al. [22] show that ⁄ may not be an appropriate GIC proxy for geomagnetic pulsation driving. Oscillations of the geomagnetic field within the ultra-lowfrequency band (1 mHz -1 Hz) may lead to active power oscillations in power networks leading to a need for appropriate control measures to maintain the grid's stability. Such active power fluctuations will be investigated in this paper.
To investigate the periodicities in GIC, Adhikari et al. [32] used wavelet analysis showing the magnitude and period (and, therefore, frequency components) of GIC. The results of the CWT revealed that GIC demonstrates a highly variable nature with time. Distinct periodicities exist when the horizontal B-field component is highly perturbed. When the magnetosphere is under quiet conditions, the GIC can be represented by smooth functions, and consequently, the wavelet coefficients show very small amplitudes. On the other hand, when GIC is highly perturbed, the wavelet coefficients are significantly large. These coefficients can identify the sudden variations in GIC which are attributed to high-frequency components.
Oyedokun et al. [21] investigated the dominant frequencies of GICs showing that 80 % of the energy within the GIC sits below 50 mHz. The authors also highlight the importance of high frequency GIC measurements for accurate replication of peak GIC values used for contingency analysis.
Most of the studies discussed previously only look into the frequency components of the B-field, E-field and the GIC. However, there are very few studies that relate the frequency components and their effects on the voltage stability of power systems. In an adjacent study [20], we use GIC data sets collected from different storms and locations to investigate the dynamic response of power systems and the impact of generator excitation control in reducing the voltage dip caused by GIC. The results show the importance of timefrequency localization in characterizing GICs for a realistic GIC-impact analysis, highlight the need for dynamic voltage stability studies and recommend adding the GMD pheonomena to the "N-1" contingency planning. The work is further extended in this paper by comparing frequencyextraction techniques and the threats to voltage stability in the presence of GICs.

B. VOLTAGE STABILITY
The IEEE PES-TR77 technical report defines voltage stability as "the ability of a power system to maintain steadystate voltages at all buses in the system after being subject to a disturbance. It depends on the ability of combined generation and transmission systems to provide the power requested by loads. This ability is constrained by the maximum power transfer to a specific set of buses and linked to the voltage drop that occurs when active and/or reactive power flows through the inductive reactances of the transmission network" [33]. In the context of GICs, harmonics flowing in the transmission network cause increased losses in the delivery system. As a result, the loadability of the network is reduced and voltage stability is more probable. The loadability of the network is affected worse in very long transmission lines having bigger series inductances and a lower power transfer capability. In such cases, series or shunt compensation is usually employed to support the system voltage.
Voltage stability can be classified according to the time scale or the size of the disturbance. Methods used to analyse the voltage stability of a network include static analysis and dynamic analysis.

1) STATIC VOLTAGE STABILITY ANALYSIS
Static analysis techniques often use P-V and Q-V curves whereby the loading is increased at a bus and the real power, reactive power and bus voltage are measured. PV curve defines two terms namely the operational margin and maximum margin. The operational margin refers to the point at which the PV curve crosses 0.9 pu voltage which represents a 10% voltage drop. The maximum margin lies at the nose point on the PV curve where a further increase in active power results in a decrease in bus voltage. With the QV curve approach, the proximity to voltage collapse can be determined. The nose point of QV curves represents the reactive power margin. When the QV curve reaches a stationary point, the stability limit is reached. The V-Q sensitivity at a bus is often used to assess the voltage stability of a system. The elements of the Jacobian matrix give the sensitivity of power with bus voltage changes. The main problem with static voltage stability analysis is that it does not provide detailed information on the mechanism of the voltage instability/collapse.

2) DYNAMIC VOLTAGE STABILITY ANALYSIS
A power system is a highly non-linear system made up of slow-and fast-acting devices. The dynamics associated with these devices need to be rigorously modelled to improve power system stability analysis. Dynamic voltage stability is concerned with the voltage instability including collapse mechanism following a disturbance, and the interaction between the different power system components leading to the continuous decline in bus voltages [33], [34], [35], [36]. The disturbance can be dynamic and is typical of geomagnetic field variations. The time constants associated with different power system equipment are also important to understand the time response of the network [12].
Time-domain simulations are used for dynamic voltage stability analysis. At each simulation time step, the power system differential-algebraic equations are solved. This method of analysis is the closest replication of the actual dynamics of voltage instability [37]. The complexity of timedomain simulations lies in the computational efforts of the CPU and the time taken for the simulations. In this regard, dynamic approaches tend to be less favourable when compared to fast-acting static methods, although dynamic approaches give a better insight into the voltage collapse phenomenon.
Mathematically, dynamic voltage stability can be assessed by using the overall power system equations expressed as a set of first-order differential equations. The algebraic equations are represented by (1) and (2).
is a vector with states of the system. is a vector of bus voltages. is a vector of current injections at the buses. is the network admittance matrix.
Equations (1) and (2) can be solved in the time domain using numerical integration methods and power flow analysis.
Over the years, GICs have been modelled as dc. To evaluate the impact of GICs on the grid, power system models are reduced to their dc equivalent with a representation of the network resistances only. As such, steady-state analyses (static voltage stability analyses) of power systems under GIC conditions have been carried out in many studies [38], [39], [40].
PV curve analysis reveals that the maximum power transfer is reduced when GICs flow in power systems. QV curve analysis identifies the amount of reactive power support required to maintain bus voltages at acceptable levels [41]. However, such analysis only shows the pre-GIC and post-GIC steady state operating condition of the power system. The mechanism of voltage collapse during GIC conditions cannot be studied using static analysis. There is a need for more detailed studies looking at transient stability analysis and dynamic voltage stability studies.
Using a time-invariant GMD, Zhang et al. [42] showed how GICs may reduce the power system's transient stability margin. The use of a constant E-field in this study was validated by the fact that GICs have frequencies much less than 1 Hz and do not vary much during transient time frames. However, control systems with time constants in the order of 5 s to 10 s such as Static Var Compensators (SVCs) and voltage controllers may become unstable if influenced by short-term GIC dynamics as depicted from the Hydro-Québec blackout during which the power grid collapsed in 92 seconds following multiple SVC tripping [43], [44]. Therefore, it is deemed necessary to investigate the shortterm voltage stability of power networks under GIC conditions using dynamic time-domain simulations.
Overbye et al. [12] investigated the effects of different GIC rise times and ramp rates on the voltage stability of a 5bus network. The results showed that the decline in bus voltage depends on the GIC dynamics. Voltage stability is influenced by load behavior especially when induction motor loads and non-unity power factor loads are concerned, tap changer action and control systems' response during short time periods in the order of seconds.
A limitation of previous voltage stability studies is that the induced E-field arising from the GMD, is modelled using a dc voltage source in series with the transmission lines. Modelling the GIC in this way neglects the frequency dependence of power system coefficients [45]. This frequency dependence is necessary because the superposition of different GIC frequency components gives rise to short-term variations in the GIC waveform that can be detrimental to the power system. During sudden commencements (fast GIC transients [46] -often associated with peak GICs), the time constants associated with network equipment and control systems may reach unstable states leading to instability.

III. THEORY DEVELOPMENT
In this section, the different methods used to identify the frequency components of GICs will be discussed. The advantages and limitations of each technique will be identified.

A. FAST FOURIER TRANSFORM (FFT)
Fourier analysis can be used to extract the magnitude and frequency components of GIC by decomposing the GIC into a sum of sinusoids. These sinusoids are assumed to have infinite duration. The main issue with the Fast Fourier Transform (FFT) is that it uses a fixed window size which loses time localisation at high frequencies. At low frequencies, frequency information is lost if the window size is too small. Fourier analysis applies to periodic and stationary signals; GICs being the opposite. GICs require a variable window size to capture the time-frequency localisation, especially during high-frequency events such as sudden storm commencements. The data sampling rate is also a major factor in accurately capturing the highfrequency components during peak GIC events using FFT and the network response [22], [21], [24].

B. SHORT-TIME FOURIER TRANSFORM (STFT)
In 1946, Dennis Gabor introduced the STFT which addresses some of the limitations of FFT. Rather than computing the Fourier transform using the same window size over the full duration of the GIC profile, the latter is segmented into narrow time windows assumed stationary within the length of the window. The GIC signal is multiplied by a window function and the Fourier transform is computed for each successive window. To pick up the frequency evolution with time, windows are overlapped and a spectrogram is generated. In this way, the frequency components appearing during specific time intervals (defined by the window) may be identified. However, the window size and window function must be carefully chosen to avoid spectral leakage giving inconsistent results.

C. WAVELET TRANSFORM
A wavelet is a short-duration waveform with an average value of zero. With wavelet analysis, the GIC is decomposed into a shifted and scaled version of the mother wavelet to extract the spectral content and its evolution with time. An irregular shifted and scaled wavelet can pick up sharp changes in GIC including the high frequency components. Unlike a sine wave of infinite duration used in Fourier transform, the wavelet decays rapidly and oscillates in time to pick up sharp changes or peaks in GIC. Moreover, wavelets offer flexibility in terms of the periodicity and shape of the underlying function (mother wavelet) to capture the time evolution of the different frequency components.
The choice of the mother wavelet is fundamental for wavelet transforms and should be based on the similarity of the wavelet to the signal being analysed [47]. For studying geomagnetic data and GICs, the "Morlet" wavelet has commonly been used due to its Gaussian-modulated shape that can capture high frequencies within GIC data [29], [30], [47].
In this study, we will focus on continuous wavelet transform (CWT). CWT tiles the time-frequency-magnitude plane with variable-sized windows, similar to STFT. For low frequency phenomena, the CWT window widens in time whereas for high frequency phenomena, the CWT window narrows.
Continuous Wavelet Transform (CWT) involves a convolution of GIC(t) with the mother wavelet which in this paper, is the Morlet wavelet.
defined the shape of the wavelet. ( ) is calculated by successively applying the change in time scale and shift such that the mother wavelet is well localised in time-frequency oscillating so that the integral = 0. Small is used for high frequency components and large is used for low frequency components.
The output of CWT are coefficients which are functions of scale and frequency with time. The scale is discretized using the number of "scales per octave".

D. PARTICLE SWARM OPTIMISATION (PSO)
The Particle Swarm Optimization (PSO) computational algorithm is a swarm intelligence method developed by Eberhart and Kennedy to solve optimization problems [48]. PSO simulates the behaviour of swarms as a simplified social system. In a similar manner to fish schools, bird flocks, and herds of animals, which use 'information sharing' approaches to adapt to their environments, find food, and avoid predators, the social behaviour of these organisms can be modelled as an optimization process.
PSO system is a multidimensional search in which each particle finds the optimal position with time. The particles adjust their position based on their own experience, their best position attained in the past ( ), and that of their neighbours, the particle with the best fitness value ( ), considering their current position and velocity.
In the PSO optimization algorithm, vectors are used to represent the individual particles whose length represents the degree of freedom of the optimization problem. PSO algorithm is done in the following steps: • First, the population of NS particles is initially distributed at random positions, 0 , and at random velocities, 0 .
• An objective function (OF), which has been defined based on the optimization problem, is evaluated for each particle of the swarm.
• According to (5), (6) and (7), the position and velocity of each particle is updated: +1 : Velocity of i th particle at (k+1) th iteration : Velocity of i th particle at k th iteration 1 , 2 : Constants having values between 0 and 2.5 1 , 2 : Randomly numbers between 0 and 1 : Inertia weight of the particle : Inertia weight at the beginning of iterations : Inertia weight at the end of iterations : Maximum number of iterations : Current iteration number : the best position of the i th particle obtained based upon its own experience : Global best position of the particle in the population

•
The new positions are checked to be within the minimum and maximum boundary conditions • The OF for each particles with updated positions is calculated and new for each particle and new global is identified.

•
Repeat 3, 4, and 5 steps for times, or until an acceptable minimum amount of OF is obtained.
• The final vector is the final solution to the problem.
PSO can be applied to identify the magnitude and frequency components of GIC. The objective is to search for the optimal set of sinusoids, that if summed up, produces the original GIC waveform.
In our approach to identify the magnitude and frequency components of GIC, we defined a mathematical time-based sinusoidal function for GIC as follows: ( ) = 1 cos(2 1 + 1 ) + 2 cos(2 2 + 2 ) + ⋯ + cos(2 + ) The following parameters need to be identified: Based on the measured GIC data and the defined function with its parameters, we defined an Objective Function (OF) to calculate the deference between the measured data and the optimized function: Ultimately, we applied the PSO optimization method to minimize the OF by finding the best values for the parameters. Based on the magnitudes of the sinusoidal terms, the algorithm finds the optimal number of sinusoids to limit the presence of low magnitude and low power terms that might represent noise. This makes it practical to simulate the optimized function in the laboratory for GIC experiments. The following have been confided in implementing the PSO algorithm: Population size, Ns : 100 Maximum number of iterations, : 10000 Learning constants, 1 , 2 : 2 Inertia weight at the beginning of iterations, : 0.9 Inertia weight at the end of iterations, ∶ 0.4

A. GIC DATA USED AND SELECTION OF TIME INTERVALS OF INTEREST
The 2003 Halloween storm was the largest of storm of solar cycle 23 and consisted of multiple CMEs and storm phases.
Overheating and damage of transformer windings in the South African network were reported in [7]. Therefore, GIC data with 2 s cadence measured in the neutral of a 400 kV transformer at Grassridge Substation in South Africa during the 2003 Halloween storm, was used for this study. Using the full representative GIC profile, four different areas (1-4) were chosen as shown in Table I

B. SIGNAL PROCESSING AND ANALYSIS IN MATLAB
Firstly, the GIC power spectrum for the full storm duration was computed using a fixed-window FFT. Secondly, by applying a fixed window size of 1000 s to capture frequencies as low as 1 mHz, a windowed FFT was carried out for each time interval or sub-areas A to P. Then, a running windowed-FFT was carried out by manually specifying frequency components of interest. After realizing the limitations of FFT in time localisation, a spectrogram of the GIC data was generated using STFT followed by a wavelet transform of the GIC data using the Morlet wavelet as the mother wavelet. Finally, Particle Swarm Optimisation (PSO) technique was used to extract the magnitudes and frequencies of the sinusoids that make up the trend in the GIC. The results from the four different analysis approaches were compared to highlight the advantages and disadvantages of each technique.

C. POWER SYSTEM MODEL AND SIMULATION
We used a similar multi-machine power system model and simulation approach to [20].The time-domain GIC dataset representing Areas A to P were imported to MATLAB Simulink and converted into a current source. The latter is driven by the GIC data points. Note that the GIC data points were scaled up by a factor 25 to increase the GIC magnitude to a level which is within the sensitivity of MATLAB transformers to saturate. There is a limitation in the MATLAB transformer models as identified in [18]. Firstly, real three-phase multi-limb transfomers have inter-phase magnetic coupling which is sensitive to dc bias. All transformers modelled in MATLAB, (three-phase five-limb, three-phase three-limb and three-phase banks) are made of three single phase units. Even though they are electrically coupled, they do not fully represent the unbalance and asymmetrical saturation that occur in reality. Other EMTP software have alternative transformer models which use the Unified Magentic Equivalent Circuit (UMEC) model but even those slightly better representations of a real transformer are still an approximation of the saturation response [49].
In this paper, we focus on the frequency components of GICs in power systems and realistic saturation of the available (limited) transformer models is achieved with the scaling factor of 25. This scaling factor has been carefully chosen to cater for high magnitude GICs of up to 225 A being the effective GIC threshold identified by FERC and the TPL-007-2 report from NERC for transformer thermal limit assessment [50], [51]. Using such a scaling factor also caters for worse-case scenarios in a power system subjected to estimated peak GICs [52], [53]. It is to be noted that there is currently no standard GIC threshold for voltage stability concerns to power systems since all networks are have different power reserves, topologies, equipment and geophysical locations.
The power system model used for the simulation is shown in Figure 6. Parameters used to model the generator, transformer, transmission line and loads are given in the Appendix. The power system consists of two areas linked by an 80 km, 400 kV line. The GIC was injected in the neutral of the Generator Step-Up transformers (GSU1 and GSU4). The load bus voltage and active power profile were captured for each simulation.
Each area consisted of two synchronous generators modelled with AC1A excitation system and turbine/governor controls.
Default MATLAB parameters were used for the steam turbine and governor system. The generator terminal voltage was initialized to 1.0 pu. The drop in the frequencydependent inductance of a transformer subjected to GICs depends on their non-linear magnetization curves. Hence, the magnetization characteristic of each transformer was specified using a piecewise linear relationship with eight different ⁄ slopes representing the non-linear This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Due to the computational complexity of the system, OPAL-RT, a digital real-time simulator (based on MATLAB/Simulink) was used to run the simulations. The simulator considers the frequency dependence of the network model, inductive delays and can, therefore, reflect realistic power systems' responses without any overruns (or errors).

V. RESULTS AND DISCUSSION
In this section, we show the limitations of fixed-window Fourier-based methods. We highlight the advantages of STFT and wavelet transform in characterizing GICs. Moreover, the role of PSO in creating a synthetic GIC for laboratory and simulation work is shown. The signal processing techniques allow us to provide clarify on questions 1 and 2 posed in section I. Most importantly, we present novel results that reveal the influence of GIC characteristics on voltage response and stability with the aim to answer question 3 from section I. Figure 7 shows the power spectrum of the GIC data for the entire duration of the storm. From Figure 7, the FFT results show that as the GIC frequency increases, the energy content in the GIC decreases. The power spectrum shows that the highest energy GICs tend to have frequencies less than 25 mHz confirming the findings from [21]. Higher frequency GICs have relatively lower energy content and therefore it is sensible to assume that the higher frequency components will not disturb power systems significantly due to their low magnitudes. However, the power spectrum uses Fourier transform with a fixed window size to compute the square of the GIC magnitudes. Using a fixed window size for the full duration of the storm does not reveal all frequencies in the GIC. For example, if a window size of 1000 s is chosen, then high frequency components such as 100 mHz components will not be identified. To capture the 100 mHz component, the window size must be at most 10 s.

FIGURE 8. GIC profile and corresponding scalogram during time interval A (Area 1 -Sub Area A) representing the sudden storm commencement
High frequency components between 10 mHz and 116 mHz with relatively high magnitudes were identified. While FFT underestimates the GIC magnitude during the SSC, wavelet analysis reveals that high frequency GICs can also have relatively high magnitudes.
If a windowed FFT is applied during time interval A starting at t = 200 s for a fixed window size of 1000 s with a fundamental frequency specified at 1 mHz, the Fourier spectrum reveals dominant frequencies below 50 mHz as shown in Figure 9. High frequency components such as the dominant 80 mHz frequencies identified by the wavelet transform, are not captured by FFT since the latter is limited by the data sampling rate and the errors introduced by FFT. It is possible to generate smaller time intervals or FFT windows (of let say 20 s) to extract frequency components of GIC greater than 50 mHz. However, the low frequency components of GIC (<1 mHz) require longer FFT windows greater than 1000 s. Varying the FFT window manually reduces the accuracy of the results in terms of identifying both low and high frequency components simultaneously including their evolution with time.
It is also possible to use a running window Fourier transform using user specified frequency components that one wants to examine. This is akin to a STFT. The window length for identifying each frequency component is taken as 1 cycle of that frequency (From t = 0 s to t = T s). After each cycle, the running FFT re-calculates the magnitude of each frequency component (From T s to 2T s and so on). Using a set of 7 frequencies in the low, moderate and high frequency range as an example, a comparison of the measured GIC and This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ sum of FFT-calculated sinusoids was made as shown in Figure  10.

FIGURE 10. Reconstruction of GIC waveform using sum of 7 sinusoids generated from FFT (top plot) and errors introduced by FFT (bottom plot).
The error margin is significant between the measured GIC and sum of sinusoids calculated from FFT. This error margin may be reduced by using more sinusoids. However, the error is non-negligible during transients and peak events such as the SSC. Although not shown in this paper, analysis of all the different areas reveals that during low frequency events (slow variations in GIC) the FFT results are in good agreement with the wavelet transform. For example, analysis of Area 2 shows that during the main phase of the storm, FFT does not pose a limitation to capture the frequency components of GICs due to the absence of transient peaks. However, FFT fails to pick up the high frequency components during high frequency events (rapid transients). FFT results lead to an underestimation of peak GICs expected to flow in power networks. Without an proper knowledge of the peak GIC values, contingency analysis may not reveal worse-case scenarios.
To find the time localization of the frequency spectra, STFT and wavelet transform may be applied as shown in Figure 11. To compute the STFT, a Hann window was used with a window size was set to 100 s with an overlap between windows of 50 %. The Morlet mother wavelet was used in the case of the wavelet analysis with specified sampling rate of 0.5 Hz (based on GIC sampling frequency).
There is a good agreement between the STFT results and wavelet transform results. Both methods are able to identify dominant high frequency components arising during peak events such as the SSC. Therefore, STFT and wavelets offer a better visualisation of the time-frequency localisation in GIC data compared to fixed window FFT. However, the window size and window function chosen determines the accuracy of the results in STFT.
Geomagnetic storms are sudden and instantaneous providing almost little lead time for recognizing the imminent threats and corrective actions by system operators [43], [54], [55]. Forecasts by the National Oceanic and Atmospheric Administration (NOAA) only provide short 15-30 minutes lead time forecast to system operators [56]. If a GMD is imminent, then system operators may need to alter the normal operation of the power system to enhance voltage stability. Identifying the temporal changes in GIC characteristics is helpful for system operators and engineers in forecasting trends in GIC, evaluating the threats to voltage stability of power systems, planning contingency measures and determining appropriate control actions required from system operators to minimise the risks of voltage collapse. Such control actions may include activation of compensators, increasing excitation of synchronous generators to provide more spinning reserve, altering Load Tap Changer (LTC) time delay deadbands and modifying relay pick up settings to prevent induction motors from stalling [12]. Other mitigation strategies include pre-determining the critical nodes in the power system using simulations and installing GIC blocking devices at these nodes to limit the impacts on the power system [13]. Such contingency planning requires precise knowledge of the temporal variations in GIC characteristics which therefore highlights the importance of identifying appropriate signal processing techniques in GIC data analysis [56].  Figure 12 shows four examples of comparison between the measured GIC and recreated GIC profile from the PSOderived sinusoids. The PSO algorithm determines the set of sinusoids that can be summed up to reproduce the original GIC waveform. If the percentage difference between the measured GIC and PSO-derived GIC is high, then, increasing n, the number of sinusoids, improves the accuracy of the results. This is shown in Figure 13.

B. ROLE OF PSO IN GENERATING A SYNTHETIC GIC WAVEFORM
Even though PSO is computationally intensive, it is a useful tool for practical tests involving multi-frequency GIC injection, for example, for testing the response of transformers and power systems to 'real' GICs. Fluctuations in a PSOderived real GIC in a laboratory test setup may aid in understanding the response of power systems to GIC. Moreoever, utilities may also use PSO-derived sinusoids to inject representative GIC models into simulation software for contingency analysis such as the "N-1" criterion. Such contingency plans and situational awareness are vital to avoid events such as the 1989 blackout during which the power network went from normal conditions to a situation where seven contingencies (N-7) were sustained within 57 seconds before collapsing 35 seconds after [43].

C. EFFECT OF GIC DYNAMICS ON VOLTAGE STABILITY
From a dynamic voltage stability perspective, we can confirm that the GIC characteristics determine the voltage response. Both the magnitude and frequency of the GIC are important. The magnitude determines the % voltage drop whereas the GIC frequencies determine the GIC dynamics such as rate of voltage decline and therefore the resulting shape of the voltage waveform. The superposition of different frequency components including high magnitude, high frequency components may disturb power system dynamics as the GIC changes magnitude and polarity. Since the GIC drives the power system's response, voltage profiles depend on the GIC characteristics during any time interval. Hence, incorporating these GIC characteristics is vital for a realistic voltage stability assessment. Figure 14 shows the GIC profile and load bus voltage during time interval A which represents the sudden storm commencement. The scalogram showing the dominant frequency components can be observed on Figure 8.

1) STORM SUDDEN COMMENCEMENT (SSC)
When the storm commences, the GIC peaks in the positive half cycle at 130.3 A (scaled GIC value). This causes an increase in the magnetising current of GSU 1 causing it to operate in the saturation region. As GSU 1 saturates, it generates even and odd harmonics that flow into the transmission line. Hence, the active power losses along the line are increased leading to a reduction in power transfer to the load as shown in Figure 15. The voltage at the load end depends on the angle between the two buses that make the transmission corridor. The load angles (between the generator buses and load bus) are shown in Figure 15. The reduction in power transfer can be explained by the drop in the angle between the buses and decline in bus voltage. In a more realistic power system model inclusive of tap changers, the OLTCs will try to restore the load power by increasing their taps. This adds additional stress on the supply system and further erodes system integrity [12].

MW load during time interval A (top plot) and load angles between generator buses and load A (bottom plot)
GIC dynamics are usually considered slow enough for the changes to be sensed by control systems or voltage regulators. However, high frequency components (short period) appearing during SSC cause rapid transients in the GIC waveform and voltage response. For example, during the SSC, the GIC can switch polarity from positive to This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022.3182237 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ negative with high peaks within less 10 -15 s. Control systems such as exciters, protection relays may easily track such changes. However, the cumulative response of different control systems and total time to issue a control action is where the problem arises. Moreover, in the case of failure of primary or secondary protection, backup protection is employed but with a time delay. Backup protection usually trips an entire substation which could lead to further voltage drop due to the lack of enough power to sustain bus voltages [43].
Sudden changes in the voltage during SSC change the operating point of control systems and may drift the states of the system towards instability. If no appropriate control action is taken to bring the state of the system back to normal, then voltage collapse may be imminent. If the voltage drops below the 10 % limit for extended durations, under-voltage relays may trip network equipment such as var compensators, during moments when their presence is highly required. The speed of voltage collapse can vary from a few seconds to tens of minutes. System dynamics during voltage collapse time frames are heavily dependent on the dynamic power reserves, timing of control devices such as LTCs and load response to voltage changes. System operators must be able to recognize voltage instability conditions and act promptly to prevent a collapse. Protection engineers need to ensure reliability of mitigation techniques by proper design, coordination and control of protective relays during conditions of harmonic pollution caused by saturating transformers from GIC [57].

2) MAIN PHASE OF STORM
The GIC profile and its corresponding scalogram for time interval B (within Area 1) are shown in Figure 16. Time interval B represents a region just after the SSC towards the beginning of the main phase of the storm. The highest magnitude GICs (up to 380 A in scaled version) were observed during time interval B. The high GIC causes increasing levels of harmonics leading to a reduction in power transfer. Therefore, the voltage drop observed at the load bus exceeds the allowed 10 % voltage drop limit. This is identified as an undervoltage condition leading to tripping of protective relays. As relays trip, network equipment such as generators and compensators may be disconnected leading to a progressive decline in bus voltages.
There is no sharp transient in the GIC waveform but rather slow fluctuations and changes in polarity. These can be interpreted from the scalogram. The scalogram shows that the dominant GIC frequencies lie between 0.92 mHz and 5.5 mHz. A relatively small transient in the GIC causes the appearance of frequencies between 6.28 mHz and 14.8 mHz. This small transient was driven solely by the geomagnetic field dynamics as the elements of the power system are unchanged. Due to the dominance of low frequency components, the rate of the decline in bus voltage is slower compared to a situation with dominant high frequency components. Figure 17 shows the GIC profile, corresponding scalogram and load voltage during time interval C. The voltage profile shows sharp voltage dips during the highlighted time intervals. The low frequency components between 0.30 mHz and 4.01 mHz cause a slow varying voltage decline. The higher frequency range of 8.38 mHz to 25.0 mHz causes a sharp and rapid voltage decline causing a shift in the state of the power system variables such as excitation controller states, tap changer control states amongst others. If control systems have time constants larger than the time it takes for the voltage to decline, then voltage instability issues may arise. The GIC profile, its corresponding scalogram and load bus voltage for time interval E (within Area 2) is shown in Figure  18. The slow variation in the GIC causes a relatively slow voltage decline. The voltage follows the GIC profile. When the GIC peaks, the voltage drop is the highest, as expected. The maximum drop in the voltage is 1 % which is not a problem for voltage stability as the voltage is within standard This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/ACCESS.2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ limits (10 % for undervoltage). No sharp changes or transients were observed in the GIC profile during time interval E. The GIC slowly increases in magnitude to a maximum and decreases slowly. The scalogram reflects the gradual slow change in the GIC as no dominant high frequency components were identified.
Protective relays are designed to operate based on conventional threats such as faults which last for a short duration compared to GIC events. When low frequency components dominate the GIC profile, the rise times in currents take much longer than the slowest backup time delay for protection schemes to operate. Hence, backup protection may be activated especially when provoked by the even and odd harmonics [43]. Backup protection is usually configured to trip an entire substation or portion of a network further aggravating the instability. Operational decisions during such events need to be clear and pre-determined from contingency planning by practicing engineers.  Figure 19 shows the GIC profile, scalogram and timedomain FFT results for Area K which is a subset of Area 3, representing the main phase of the storm.
It is observed that the main phase of the storm can also have rapid, sudden fluctuations in the GIC pattern, but not as sharp as the SSC as depicted from Table II showing the maximum rate of voltage drop during different storm phases. Therefore, there is less probability of voltage instability during the main phase of the storm compared to the SSC if protection schemes operate correctly. Figure 20 shows the GIC profile, scalogram and load bus voltage during time interval L. During this period, geomagnetic pulsations were observed that can drive significant GICs causing rapid and sudden GIC dynamics. During such events, voltage stability becomes a threat, like the SSC. The sharp transients observed in the GIC profile causes the appearance of relatively higher frequency components compared to Areas E-J. However, if one compares Area K with the SSC, the sharp transients observed during time interval K are not as "peaky" as the SSC. Therefore, during the main phase of the storm, we do not expect to see high frequency components as observed during SSC.  Figure 21 shows the GIC profile, scalogram and load bus voltage during time interval N which represents the recovery period of the 2003 Halloween storm. The second trough is more effective in the voltage response due to cumulative driving and the voltage drops to a lower value than the first trough. The frequency components around the first and second trough are within the same band. However, the magnitude of the GIC components has increased around the second trough. Therefore, the voltage drops to a lower value due to the higher GIC magnitude.   Figure 23 shows the effect of GIC sampling frequency on the load bus voltage during three representative intervals within the Halloween storm, each with different frequency components dominating: high frequency event (for example the SSC), low frequency event (such as the main phase) and pulsation activity. These novel results show that, during highly disturbed intervals where high frequency components dominate (such as the SSC), 1 minute GIC sampling data may not fully capture the short-term GIC dynamics and hence voltage fluctuations arising from the temporal changes in GIC characteristics. Instead, 2 s GIC data provides more useful information on the transient peaks in GIC and drops in voltageOtherwise, low GIC sampling rates may underestimate the voltage drops during SSCs.

D. EFFECT OF GIC SAMPLE RATES ON VOLTAGE
Intervals dominated by low frequency GIC components are associated with slow changes in voltage. In such storm intervals, low sampling rates for GIC suffices in reflecting the effects on the voltage.
Intervals associcated with geomagnetic pulsations (narrow-band fluctuations) are characterized by similar voltage oscillations . During such instances, GIC sampling rates lower than the pulsation activity lead to an underestimation of the voltage drop in the network as observed from Figure 23.
Conventional GIC modelling for voltage stability analyses uses a dc model to represent the GIC. This dc approximation appears to stem from the low GIC sampling rates due to hardware limitations in the past. Threats to dynamic voltage stability of power systems during SSCs may be underestimated when using the conventional dc approximation or low frequency GIC sampling as proved for the first time in this paper.

E. SUMMARY OF RESULTS
In this section, we provide answers to some pertinent questions raised in Section I: 1. Which technique is most suitable to identify the frequency components of GIC? In this paper, we focused on the temporal evolution of GIC frequencies and their effects on voltage stability. Four signal processing techniques were used to identify the frequency components of GICs. Applying FFT using a fixed window 2 s data 1 min data 2 s data 1 min data 2 s data 1 min data This article has been accepted for publication in IEEE Access. This is the author's version which has not been fully edited and content may change prior to final publication. The four signal processing techniques provide the basis to explore mitigation and risk management strategies based on the temporal evolution of GIC characteristics including magnitude and frequency. Implications relate to the training of system operators to heighten their situation awareness under GMD situations [55]. Moreover, understanding the time-frequency localization of GICs provides much better perspective to the power utility in revising equipment design, protection schemes including time delays between primary, secondary and backup protection, and contingency measures based on the "N-1" criterion. Table II summarises the range of dominant frequencies, and rates of voltage decline during the three main storm phases investigated in this paper. The SSC is characterized by high frequency components that may have relatively high magnitudes. The high frequencies lead to rapid fluctuations in the voltage which can change polarity within small time intervals of 10 s. Hence, the rate of voltage decline during the SSC is much larger compared to the main phase and recovery period. Pulsation activity during the recovery phase may also cause faster rate of change of voltage compared to the main phase in the data analysed.

What are the typical GIC frequencies and rates of voltage decline during different storm phases?
3. What instant of the storm is more of a concern for power systems? Both high frequency and low frequency components of GICs are a problem to power systems. High frequency components cause rapid, sudden transients in the GIC which reflect on bus voltages. If the GIC suddenly increases in magnitude up to an extent that conventional var-support is not available, then undervoltage relays may trip network equipment as the voltage goes below permissible limits. Moreover, control systems with large time constants may reach unstable states leading to incorrect or unnecessary control actions. Therefore, forecasting GICs more accurately will allow system operators to manage risk more effectively [56]. During sudden storm manifestations, operators are usually left with no time to realise the imminent threats to voltage stability of the power system. NOAA only provides 15-30 minutes lead time forecast to network operators for operational control actions [56]. Signal processing techniques discussed in this paper may be used to forecast the trend in GIC allowing enough time for system operators to take appropriate corrective actions especially during the SSC.
Low frequency components do not pose a problem to the voltage stability directly. The low frequency components (longer period) cause gradual slow changes in GIC which causes a slow voltage decline. If corrective control actions can be taken within such time frames, then voltage stability will not be at threat. The relatively high magnitude, low frequency components are more of a problem to transformers since the low frequency GICs persist for long and cause transformer heating, insulation degradation and cumulative damage as reported in [7] If the transformer goes out of service either due to a transformer protection relay trip or physical damage, then the low frequency components indirectly contribute to voltage instability. Therefore, the GIC/voltage transients (caused by high frequency GIC components) such as those observed during SSC and cumulative driving from geomagnetic pulsations, are more of a concern for voltage stability whereas the main phase of the storm is a concern for power transformer thermal assessment [44] .

VI. CONCLUSION
In this paper, we compared four different techniques used to identify the frequency components of GICs. The impacts of both low frequency and high frequency components of GICs on voltage stability were analysed using time-domain simulations. We looked at different driving regimes within the different storm phases of an extreme storm, the Halloween storm. To do this, we probed different spectral characteristics throughout the storm. The results reveal that wavelet analysis and STFT brings about a more comprehensive view of the time-evolution of GIC frequency components including their varying amplitudes. GICs with both high magnitudes and frequencies are picked up by wavelets, which are otherwise hidden by FFT due to its limitations of fixed window size, data cadence and assumption of stationary signals. Moderate to severe voltage instability are positively correlated to high frequencies such as during SSC and geomagnetic pulsations. While it may be difficult to generalize the results, the latter advances understanding of voltage stability threats in the presence of GICs. Signal processing techniques proposed in this paper may be used to forecast the trend in GIC. Dynamic time-domain simulations provide information about the critical buses and expected voltage drops. Using such information for contingency planning and implementing appropriate control measures provides a more representative assessment of the dangers from GMDs. Future work investigates harmonic analysis, protection and delays in network responses with real GIC and saturable transformer models.