Wavelet and Network Analysis of Magnetic Field Variation and Geomagnetically Induced Currents During Large Storms

During geomagnetic storms rapid magnetic variations cause large, sharp enhancements of the magnetic and geoelectric field at mid‐latitudes. These present a potential hazard to grounded technology such as high voltage transformers, pipelines and railway systems. Spatiotemporal quantification can provide insight into the magnitude and configuration of their potential hazard. We use the Haar wavelet transform to localize in time and frequency the storm‐time response in both European ground‐based magnetometer measurements and modeled geomagnetically induced currents (GICs) from the high voltage grid of Great Britain (GB). Wavelet cross‐correlation of the GIC in the grounded nodes is then used to build a time‐varying network of GIC flow around the GB grid during storms including the 2003 Halloween storm. We find a highly intermittent (few tens of minutes duration) long‐range coherent response that can span the entire physical grid at most intense times. The spatial pattern of response seen in the GIC flow network does not simply follow that of the amplitude of the rate of change of B field. Coherent response is excited across spatially extended clusters comprised of a subset of nodes that are highly connected to each other, with a tendency for east‐west linkages following that of the physical grid, simultaneous with the overhead presence of the auroral electrojet and the inducing component of the magnetic field. This can quantify the spatial and temporal location of increased hazard in specific regions during large storms by including effects of both the geophysical and engineering configuration of the high voltage grid.

geoelectric field during geomagnetic storms which exhibit a range of different peak values at both local and regional scales (Dimmock et al., 2019).
It is well-established the main cause of damage in grounded technological infrastructure such as the high voltage (HV) power grid or gas pipeline network is from geomagnetically induced currents (GICs) generated during large storms (e.g., Albertson et al., 1981;Kappenman, 2003). The magnitude of GICs are dependent on the rate of change of the magnetic field in combination with the local ground conductivity structure and the resistance and topological parameters of the grid (e.g., Pulkkinen et al., 2001;Thomson et al., 2005). As modern society becomes increasingly dependent on continuous electrical and energy supply, the potential economic cost of prolonged widespread loss of power is estimated at tens of billions of dollars per day (Eastwood et al., 2017;Oughton et al., 2018). This strongly motivates the desire to better understand the various components of the hazard, in particular the extent and exposure of regions to rapid rates of magnetic field change and the connection between these events and the impact on grounded technology.
Quantifying the hazard for a given location is non-trivial due to the interaction between time-varying magnetic fields and regional geological structures. For example, Ngwira et al. (2015) found that large differences in the magnetic and geoelectric field peaks can arise across different latitudes during a geomagnetic storm while Beggan (2015) demonstrated the instantaneous location of the auroral electrojet is the leading factor in driving modeled GIC in Great Britain (GB) high voltage grid. Analysis of recent large storms at high latitudes such as the September 2017 event by Dimmock et al. (2019) shows that, even with comparatively low Dst values of around 150 nT, unexpectedly large variations in the rate of change of the magnetic field were recorded in Scandinavia. Variometer records in and around the North Sea showed similar large rates of change were experienced further south between  50 N and  60 N (Beggan & Marple, 2018;Simpson & Bahr, 2020b).
Quantifying the extremes of the magnetic field change during the era of digital magnetic records (from the 1980s onward) gives an insight into the potential magnitude of variation both in Europe (Thomson et al., 2011) and globally (Rogers et al., 2020) though only in the context of single observatory records rather than across multiple observation sites. From these studies, the geomagnetic storms of March 1989 and October 2003 rank as two of the largest within the past 40 years, though many larger instrumentally recorded storms have exceeded these, such as the Carrington storm of September 1859 (e.g., Siscoe et al., 2006) or May 1921 (Hapgood, 2019). Although there is much debate regarding the true magnitude of the Carrington storm, for example, whether it represents a truly extreme outlier (or "dragon-king") event or is more mundanely in the category of a large "super" storm (Chapman, Horne, & Watkins, 2020;Siscoe et al., 2006), an event of a similar magnitude today would likely cause serious issues for the safe operation of modern technology. For instance, Baker et al. (2013) document a "near-miss" in July 2012 from a Coronal Mass Ejection sampled by Stereo-A, directed away from Earth, which produced an IMF of almost 100 nT, larger than any other satellite measurement of the past few decades.
The likelihood of extreme events can be assessed from the long-running aa-index with Chapman, McIntosh, et al. (2020) noting that severe storms can occur at any time although their likelihood is modulated by the solar cycle (around 1%-3% of all severe storms have occurred around minimum). Thus, there remains a large level of uncertainty around the hazard posed by space weather in terms of intensity, frequency and the pathways in which it could cause social and economic damage. In particular the effect of simultaneous enhancement across a wide region is potentially more difficult to manage than local enhancements; Boteler (2019) and Rajput et al. (2021) offer insights into issues within power grids during large storms.
The link between the extent of the auroral oval and the amplitude of the fastest rates of change of the magnetic field is of increasing interest, particularly as these parameters are important in determining the magnitude and extent of the potential hazard to grounded technology. An investigation of the rate of change of magnetic field by Freeman et al. (2019) at GB observatories showed that the substorm expansion phase produces the most extreme rates of change. Though the expansion phase of the auroral oval is usually a relatively short interval within the overall duration of a substorm, these intervals contain a high percentage of incidences of the largest rates of change of the magnetic field. Indeed, in any given storm, several of these expansion phase intervals may occur creating an extended period of risk. Studying the manner in which ORR ET AL. 10.1029/2021SW002772 2 of 30 auroral oval expansion leads to large rates of change of the magnetic field is useful in characterizing the general properties of the field rate of change, as well as offering insight into the behavior of large storms.
In this study we first use wavelet analysis to investigate the location and intensity of ground magnetic perturbations during large geomagnetic storms. Wavelet transforms have previously been used in magnetospheric physics to analyze the sunspot cycle (Ochadlick et al., 1993), geostationary satellite data of substorm breakup (Holter et al., 1995), current disruption events (Lui & Najmi, 1997), magnetospheric-ionospheric coupling (Smith et al., 2020) and Pi 1/2 pulsations (Milling et al., 2008;Nosé et al., 1998). We will use a Haar wavelet decomposition that provides a time-frequency localization of fast changes in the signal. This is applied to SuperMAG ground-based magnetometer observations over northern Europe; however the method is applicable globally. Several intense geomagnetic storms are analyzed in order to determine the location, and in particular the latitude, of the fastest rate of change of ground magnetic field. We next apply these ground magnetic perturbations as input to a model for the GB high voltage grid. The model incorporates details of the physical grid structure and ground conductivity and as output models the time-varying GIC at HV grid transformers. We will quantify this spatiotemporal GIC response using network analysis. Network analysis based on Fourier cross-correlation have been applied to SuperMAG data previously to characterize substorm dynamics (Dods et al., 2015;Orr et al., 2019Orr et al., , 2021 and the response to northward and southward turnings of the interplanetary magnetic field (Dods et al., 2017). To perform this network analysis on the GIC flow across the GB power grid we use the Haar wavelet decomposition and cross-correlation. We obtain the wavelet cross-correlation between all pairs of transformers in the HV grid model. This generates a cross-correlation matrix between all pairs of observations, with matrix elements that vary with time. Thresholding this matrix gives a time dependent, or dynamical network which then characterizes the extent of spatially coherent GIC response. This characterizes how different, distant parts of a highly connected power grid can become coupled during intervals of intense geomagnetic activity.
The next section details the ground-based magnetometer observations and modeled GIC used in the study. Section 3 presents the results of the wavelet decomposition of multiple ground-based magnetometer observations for the 2003 Halloween storm across Europe (the analysis is repeated for the March 1989 and September 2017 storms in the appendix). Section 4 catalogs the output of the network analysis of GIC for a synthetic storm and for the October 2003 storm. We place the results in context in Section 5. Technical details of the wavelet decomposition and network analysis are given in the appendix.

Data Sets
High quality variometer and magnetometer data are freely available from a number of sources including IN-TERMAGNET, the World Data Centers for Geomagnetism and the SuperMAG database. SuperMAG collates magnetometer observations and processes them with a common baseline removal technique, standardized coordinate system and sampling to a 1-min cadence (Gjerloev, 2012). We use the SuperMAG database to analyze the intensity of the rate of change of magnetic field at the Earth's surface during the Halloween storms of 2003, using data from all available European magnetometers (the analysis is repeated for two other examples, the March 1989 and September 2017 storms, see Appendix A3). Figure 1 (left) shows the spatial configuration of the 50 magnetometer stations used in the analysis of the Halloween storms.
For the GIC analysis, we generated time series of GIC flow in the GB high voltage grid based on the methodology of Kelly et al. (2017). GIC data (in Amperes) were calculated using a high voltage grid model for the year 2017, containing 398 transformers each with a fixed grounding resistance of 0.5 Ohm at the earthing point (Figure 1, right). The transformers are connected by approximately 1,200 power lines ranging from 400 m to 189 km in length, with a mean of 29 km. The resistance of each power line is determined from information provided in the National Grid UK Ten Year Statement (National Grid ESO, 2020).
The GICs are driven by geoelectric fields computed using minute mean magnetic data extrapolated over the region using the Spherical Elementary Currents System (SECS) method (Amm, 1997) in combination with a thin-sheet model which incorporates a detailed 2D conductivity model with 1D variation at depth (Beamish & White, 2012;Beggan, 2015). The SECS method uses magnetic data from a number of magnetometers around the North Sea region to extrapolate the field over the entire region onto a 10 × 10 km grid.
ORR ET AL.

10.1029/2021SW002772
3 of 30 Depending on the storm, weak spatial smoothing is applied to the magnetic field extrapolation to reduce the effects of the limited number of stations available (McLay & Beggan, 2010).
The differences of minute mean values are used to compute the regional geoelectric field across the British Isles using a fixed frequency of 120 s. The use of a single frequency does mean energy associated with other frequencies can be missed (as described in Kelly et al. [2017]) but for the purpose of demonstrating our method it is sufficient for the GIC values computed to be representative rather than exact. The choice of frequency is somewhat arbitrary but with 1 min cadence for the magnetic field data, the shortest periodic signal we could capture is 2 min. Beggan et al. (2021) shows periods of 2 min tend to generate larger geoelectric fields which correspond well to observatory measurements. However, it is possible to include a wider spectral representation of the magnetic field using the method of Divett et al. (2020), for example, and on-going improvements to the UK conductivity model should allow a better approximation in future. The output model consists of geoelectric field values on a 10 km grid cell map covering the area shown in Figure 1 (right). Examples of the conductivity model for the GB with thin-sheet code can be found from Beamish (2013) and Ádám et al. (2012). The geoelectric field values from the thin-sheet have been compared against measurements at three GB observatories and found to correlate consistently highly though generally tend to underestimate the true size of the field during storms (Beggan et al., 2021). The outputs of the high voltage grid model have, at present, been validated in eastern Scotland through comparison with indirect measurements of line GIC using the Differential Magnetometer Method (Hübert et al., 2020). Again, on-going work continues to validate the remainder of the network.
We produce a series of snapshot maps of the geoelectric field across Britain each minute. We concentrate our analysis on the October 2003 ("Halloween") storm as it was both intense and well sampled by a large number of variometers at the time, we have also computed GIC models for the March 1989 and September 2017 storms, but their analysis did not reveal significant additional information.
ORR ET AL.

Mapping Magnetic Field Rate of Change Using Wavelets
As the hazard posed to grounded technology is directly linked to rapid magnetic field change, we first demonstrate that wavelets can effectively map the spatial variation of the field's rate of change. We focus on October 29-31, 2003 as this represents the largest storm in the past 20 years for which a large number of magnetometer stations are available. We have also examined on March 13 and 14, 1989 storm, which was larger in magnitude but has fewer stations available, and on September 7-8, 2017 storm, which was smaller but has more stations available (shown in Appendix A3).  the peak of the storm in the European sector. The wavelet decomposition of the magnetic field is shown in panels (b, d, and f). These plot the value of C j k ,  where  1,2,3 j is the wavelet scale. These panels capture the latitudinal distribution of enhanced fast change in ground magnetic field. The wavelet analysis clearly resolves the sharp changes associated with the storm. Before 06:00 UT there is very little activity below 60° magnetic latitude, but shortly after midnight sharp changes in magnetic field (i.e., high C j k ,  ) are measured at increasingly lower latitudes. A peak across all scale lengths is observed at ∼21:00 on 30th.
Panels (c, e, and g) overlay the C j k ,  from each magnetometer (gray, left ordinate) and the average C j k ,  of all magnetometers (orange, right ordinate). Note the most intense rates of change are in the  1 j , 2-min scale length and can be highly spatially localized, the  1 j peak value for a single station (gray) can, for the sharpest changes, be several times that of the average (orange). The figure also shows that the rate of change remains high across a wide latitude and spectral band for a prolonged period of time. The largest rates of change (wavelet power) are at the sudden storm commencement and onset, but large rate of change also occurs in the middle and the end of the first phase of the storm during the 29/30th and is also embedded in the second main phase later on during October 30. Having demonstrated that our wavelet analysis provides an effective procedure for time-frequency localization of the rate of change of magnetic field at each station, we will now apply it to the analysis of GIC flow in the GB high voltage grid in conjunction with a network cross-correlation analysis. Note, we will use the term "network" exclusively to describe spatially coherent GIC response as connections in the network; it does not refer to the physical connections of the high voltage grid itself.

Network Analysis of Geomagnetically Induced Currents
The flow of GICs in a power grid is a function of the spatial and temporal distribution of the geoelectric field in conjunction with the topological layout of the overlying high voltage grid. The geoelectric field is influenced by the magnetic field and the regional geology which can give rise to large variations over tens of kilometers (Lucas et al., 2020;Simpson & Bahr, 2020a). The response of the high voltage grid is then dependent upon the geoelectric field and the location and number of earthing points, the length and orientation of the power lines and the resistance parameters of the wires and transformers.
In the UK, the highest voltage power lines operate at 400 kV with a line resistance of around 0.006 /km. The 275 and 132 kV lines have slightly larger resistances at 0.01 /km and 0.03 /km, respectively. The grounding resistance of the individual transformers is also important to know as lower resistance earthing points experience higher GIC flow and vice versa. Modeling of GICs is based on combining these features together in a matrix formulation to produce an instantaneous snapshot of the flow in and out of the ground across the grid; we use the formulation of Lehtinen and Pirjola (1985) and have verified our modeling code against the benchmark of Horton et al. (2012).
An important feature for understanding the hazard to transformers is to examine how the GIC flow in each node co-varies with its neighbors and within its local region. Large GICs (>25 A) can begin to affect the behavior of older types of transformers (e.g., Divett et al., 2018) though it is difficult to generalize given the complex nature of the power grid. The GB grid is considered to be resilient to small to medium space weather effects given its strongly connected nature (Cannon, 2013;Oughton et al., 2018). However, as the connection lines between transformers are often complicated it may not be apparent how GIC flows across regions or, in the case of large and widespread geoelectric fields, between regions. An unknown factor is how GIC travel over longer physical distances; for example, whether locations in the south of England become more exposed if transformers in central Scotland experience large GIC. We investigate this using a network analysis based on wavelet cross-correlation of the signals from spatially distributed GIC grounded nodes, to identify strong connections and pathways between nodes.

Application to a Synthetic Storm
To understand how network analysis of GIC modeled from a geomagnetic storm should be interpreted, we first model the GIC response of the grid to a synthetic geomagnetic storm. A synthetic auroral electrojet ORR ET AL.

10.1029/2021SW002772
6 of 30 intensification is represented with a simple cosine half-wave shaped spatial variation in the north-south direction with a total width fixed of 400 km. The strength of the magnetic field at the center of the auroral electrojet varies sinusoidally between 1,000 nT with a period of 61 time-steps where each time-step nominally represents a minute. The electrojet location dwells successively at three locations over the UK: northern Scotland, central England and then the southern coast of England, moving instantaneously between each of these locations after each 61 time-step cycle. At each time-step the geoelectric field is then computed using the thin-sheet model. A GIC value (in Amperes) for each node in the high voltage grid is then calculated from the geoelectric field at each instant. We then perform our wavelet analysis of the GIC at the 398 nodes for this test scenario. Figure 3 shows the analysis of the modeled GIC. The electrojet is initially over northern Scotland such that at time index 1, the model magnetic field is at zero, rising to 1,000 nT by time index 16, crossing through zero at time index 31 and returning to zero at time index 61. This is then repeated above central England from time index 62-122 and over south England from time index 123-183. Panel (a) plots the modeled GIC data for each grounded HV grid node ordered by latitude (non-uniformly spaced). The mean and maximum magnitudes of the synthetic GIC are 1.4 and 26.5 A, respectively. To emphasize the high magnitude GIC we have used a natural log (ln) color scale in panel (a), where the mean corresponds to 0.3 A and the max to 3.3 A. The wavelet coefficient estimates of the rate of change of the field at scales  1,2,3 j are plotted in panels (b, e, and h); these panels are in the same format as panels (b, d, and f) of Panels (c,f,i) overlay the C j k ,  value at each HV grid node (gray, left ordinate) and the average C j k ,  over all grid nodes (light red, right ordinate) for the  1,2,3 j wavelet coefficients respectively. These panels are in the same format as panels (c, e, and g) of Figure 2. We again see a significant spatial variation across the grid nodes, the rate of change at individual grid nodes (gray) can significantly exceed that of the average (red), the largest by a factor of 2-3. In panel (c) the stepwise change in electrojet latitude is picked up by the highest time resolution,  1 j wavelet coefficient where we see a dip in intensity.
The final panels in Figure 3, panels (d, g, and j), combine the latitudinal dependence of the modeled C j k ,  at each grid node with an indication of the spatial coherence between the nodes. Spatial coherence is quantified by calculating the time-varying wavelet cross-correlation, using a 30 "minute" running window, pairwise between all pairs of HV grid nodes using the modeled C j k , Since the wavelet transform provides a decomposition of the signal that is, time-localized, the cross-correlation output is on the same time resolution as the wavelet scale. The full network construction method is described in Appendix A2. Panels (d, g, and j) plot as colored circles the C j k ,  for the subset of most active nodes that satisfy two criteria.
First, for each node, the , j k C  | | must exceed the 80th quantile of all values that occurred at that grid node throughout the simulated time interval. Second, the wavelet cross-correlation between the node and at least one other must exceed a threshold (here | 0.995|). In order to highlight sharp changes with time, timings refer to the cross-correlation window leading edge. The cross-correlation network is at the time resolution of the wavelet, 2 j (2, 4, and 8 time units for panels (d, g, and j), respectively).
These threshold criteria identify the grid nodes, and the connections between them, that form the cross-correlation network. Spatial coherence between the GIC at grid nodes is then captured by the local network degree at each node, that is, the number of other nodes in the network that connect to each node. For example, if a node has 1 connection it will have a degree of 1. The size of the circles in panels (d, g, and j) indicate the degree of each node in the network, and the color plots the natural log of C  Nodes which have particularly high levels of connectivity, that is, a degree that exceeds the 90th quantile of the degree distribution for the entire simulated storm (here a degree of 114, 116, and 115 for panels (d, g, and j), respectively) are plotted with black edged circles. Figure 3 shows that both the individual wavelet coefficients (panel b, e, and h), and the cross-correlation as captured by the network degree (panel d, g, and j) respond strongly to high rates of change of the signal and its latitudinal dependence. Enhanced cross-correlation of the transformers in space coincides with the location of the synthetic electrojet. The network degree also reveals underlying latitudinal variation in the modeled HV grid response to the GIC, the greatest level of spatial correlation occurs when the electrojet is sited across the center of the island (panels d, g, and j). The wavelet coefficients, which capture the power in the fast changing part of the signal, track the latitudinal dependence. Figure 3 indicates the latitudes at which a threshold for spatial coherence (threshold on local network degree) is exceeded. We will now look in more detail at how spatial coherence varies geographically, by mapping it for the times indicated by the vertical colored lines.
The connectivity of the network is mapped spatially in Figure 4. This figure shows snapshots of the network calculated from the most time-localized  1 j wavelet cross-correlation of the GIC at the 398 grounded nodes modeled in the test scenario. The plots are selected for time indices 32, 92, and 152, indicated by colored vertical lines on Figure 3, where the modeled GIC is most rapidly changing. Grid nodes are included in the network if the wavelet cross-correlation with at least one other node exceeds a threshold and we plot the network obtained for three increasingly high thresholds from left to right. The size of the circles again indicates the degree of each node in the network, black circles indicate the most connected nodes (that have a local degree exceeding the 90th quantile for the entire simulated storm) and the color plots the rate of change natural logarithm C k 1,  . Note the connecting lines are not physical power lines (as is the case in Figure 1) but indicate strong correlation between the nodes occurs within the preceding fixed 30 min window in time.
From this sequence of plots we can see the spatial pattern of network connections vary as the synthetic auroral electrojet moves southward over time. In the upper panels when the modeled electrojet is over Scotland, nodes within the central region of Scotland are strongly connected to each other, but there is no strong connectivity elsewhere. When the center of the electrojet has moved to central England the network fully connects both the Scotland and central England regions to each other, down to around  52.5 latitude. Finally, when the electrojet is sited on the southern coast the strong connections only extend as far north as the central part of England.
The modeled electrojet is fully spatially coherent and Figure 4 shows that the response of the HV grid model is also spatially coherent. A given cross-correlation threshold then effectively determines a spatial region within which all the nodes are fully connected to each other as they respond to a signal which is in phase across all the nodes. Any localized phase delays in the modeled HV grid are therefore short compared to the time variation of the modeled driving electrojet. Within each of these networks, the nodes all have the same degree, the value of the degree simply reflects the spatial extent (number of nodes encompassed by the network) of the HV grid response. Each time window chosen here includes half a period of the modeled cosine-shaped wave, therefore if the first network criteria ( , j k C  | | exceeds 80th quantile) of a node is met, the remaining nodes all have a cross-correlation of 1.
As seen in the previous figure, the rate of change of the GIC response of the HV grid is largest when the electrojet is over central England. This corresponds to a network which extends across Scotland and the middle of England and as a consequence, all nodes with local degree above the 90th quantile of the degree distribution taken over the entire modeled storm duration (black circles) are found in this network.
These details of the cross-correlation response depend on the details of the HV grid and ground conductivity model. This provides a benchmark with which we can now compare the effect of a physically observed event. Figures 3 and 4 plot the response to an electrojet which varies smoothly in latitude, has no longitudinal variation, and, once centered on a given region, varies smoothly in time. With this in mind, we apply the method to GIC computed from the October 2003 storm.
ORR ET AL.

Application to Halloween 2003 Storm
We now apply the same methodology to the Halloween 2003 storm, which reached a Dst value of 383 nT at the peak of the activity around 21:00 UT on October 30. The north and east components of the magnetic field over the British Isles were modeled using Spherical Elementary Current Systems (SECS) to extrapolate measurements from nine observatories across a 1,680 by 1,200 km grid to match the input required for the thin-sheet modeling code (McLay & Beggan, 2010). The output geoelectric field maps for every minute were used to compute GICs during the storm. We point out that modeled GIC during the geomagnetic storms in question are indicative of where large GIC flow and the regional variation across the power grid rather than being the precise value expected to be experienced in any given node. The wavelet analysis clearly identifies the sharp changes associated with the storm on the time scales of minutes. Prior to 06:00 UT there is little activity, but after 06:00 rapid and large changes in magnetic field (high C k 1,  ) are measured. Shortly after 06:00 nodes across a range of latitudes have a high degree (i.e., are highly correlated). Due to the inherent bandpass filtering of the geoelectric data for the thin-sheet modeling code used here, the very rapid variation during the sudden storm commencement are damped, though it is clear that a large pulse of correlation across all latitudes occurs at 06:00 UT on 29th. As the storm progresses, the nodes with the highest degree become centered around a latitude of   56 N. For much of the storm, intervals of strong correlation between nodes are established in the Midland Valley of Scotland (  56 N) shown in the closed black circles. Much of the network is correlated to some degree during intermittent intervals over the main phases of the storm, and importantly, high spatial correlation (high network degree) does not simply correspond to peak wavelet power (peak rate of change); this can also be seen in Figures A3 and A5 which show the same analysis for the 1989 and 2017 storms.
The geographical map of the network at different snapshots in time is shown in Figure 6, which follows the format of Figure 4 (the same information is plotted in Figures A4 and A6 for the 1989 and 2017 storms). The times correspond to the vertical colored lines in Figure 5. The panels show the network calculated with increasingly high correlation threshold from left to right. Again, the size of the circles indicate the degree of each node in the network and the color plots the natural logarithm of  connections. Both the degree and wavelet coefficient, C j k ,  , of the HV grid response can now be seen to vary geographically across the correlation network, reflecting the response to the physical, rather than modeled, electrojet. In the analysis that follows we will see that the spatial pattern of cross-correlation is highly variable. We begin by showing snapshots that illustrate the different network topologies that can arise.
The uppermost panels (at 06:13 UT) show the initial onset of the storm and we can see that the entire network is highly, but not fully, connected. Almost all of the network nodes have a degree which is above the 90th quantile for the storm and there is a prevalence of long-range connections extending across the entire HV grid. There is strong correlation between parts of the HV grid that are geographically distant. This dominance of long-range connections at the onset of the storm is also seen in the 1989 and 2017 events.
ORR ET AL.  In the second set of panels (at 14:25 UT) we can see that the character of the network has changed significantly. The network is dominated by connections spanning the midlands and extending south, there is a second, much more weakly connected, network in Scotland, which at high correlation thresholds is centered mainly on large populations in Scotland. Long-range connections (Scotland-midlands-south) are only found at lower thresholds for cross-correlation so that unlike at (06:13 UT), the strongest correlation is geographically more localized. The intensity of the  1 j wavelet coefficient is highest for the dominant (midlands-south spanning) network which contains many nodes with local degree above the 90th quantile. Comparing these snapshots it is clear that once a threshold is exceeded, long-range connections do not necessarily follow from high levels of wavelet coefficient, that is, localized fast rate of change of GIC.
This distinction between localized wavelet intensity and long-range correlation is also apparent from the remaining two snapshots. The most intense point of the storm as seen in Dst is around 21:23 UT which is shown in the next set of panels. We can see that the wavelet coefficient intensity has indeed reached its highest value across the network. However the network is not as globally connected as seen earlier in the storm. Long-range connections are seen at lower thresholds; at the strongest connections are east-west and at particular latitude bands. These east-west connections include the most highly connected nodes seen at earlier times. The bottom set of panels (at 00:23) also shows these east-west connections, but notably, the regions of strongest correlation (highest degree) are now at higher latitude. The prevalence of these east-west long-range connections in the network are also seen in the 1989 and 2017 events. They correspond to the structure of the 400 kV HV grid which contains long E-W connections at roughly constant latitudes where population centers lie: (a) passing through central Scotland, (b) extending from Anglesey, through N Wales, Midlands to eastern England, and (c) extending from South Wales to Thames Estuary.

Mapping Network Connections Across a Storm
The rapid change and morphology of the magnetic field and corresponding geoelectric fields are reflected in the highly dynamic nature of the network connectivity. We can parameterize the time-varying network topology by counting the number of connections within and between distinct geographical regions. The island of Britain was divided into three regions: north, middle, and south. Nodes in the power grid north of  55 latitude were grouped together in the "north" category while those south of  53 latitude are in the "south" category. "Middle" nodes are considered to lie between  53 to  Figure 7 plots the number of connections (edges) in the network as a function of time for the duration of the storm (again, for the  1 j wavelet cross-correlation). The time resolution for the network is that of the original data (1 min, details in Appendix A2). The total number of possible connections within and between each spatial region is different for the distinct spatial regions; the modeled HV grid has a relatively high density of spatial sampling in the North compared to other regions. We therefore plot both the relative and absolute connectivity on the left and right axis respectively. The relative connectivity, that is, the number of connections between and within each region, normalized to the total possible number, is plotted on the left ordinate, whereas the absolute connectivity, that is, the raw number of connections, is plotted on the right ordinate. A relative connectivity of 1 in a given region would be a 100% connected network in which every node is connected to all others and the grid nodes respond fully coherently to a fully spatially coherent driving electrojet. The left ordinate then has the same scale for the different regions, whereas the right ordinate varies by about a factor of five. The network is calculated using the same method as described above where (a) for each node, the , j k C  | | must exceed the 80th quantile of all values that occurred at that grid node throughout the simulated time interval and (b) the magnitude of the wavelet cross-correlation between the node and at least one other must exceed a threshold. The number of connections obtained for different network cross-correlation thresholds of 0.7, 0.8, and 0.9 are overplotted, (colors indicate the threshold level). The thresholds used to construct the network are lower than those above ( 0.995), resolving a broader range of network structure and consequently coherent activity across the HV grid.
ORR ET AL.

10.1029/2021SW002772
14 of 30 Higher thresholds of connectivity systematically yield smaller numbers of connections. Comparing the number of edges at different thresholds gives an indication of how variable the level of coherence is across the network at any given time. The level of connectivity is highly intermittent in time, and a given network state (as shown in the snapshots of Figure 6) generally persists for 10 min or less. The variability in coherence across the network is also highly time dependent. Nevertheless, there are overall trends within and between these regions throughout the storm.
The upper three panels show connections within each of the three regions. The "North" panel has the highest absolute connectivity (right scale) but a relatively low relative connectivity, it is less than 20% connected. The "South" region is the most highly connected, peaking at more than 50% connectivity, with absolute connectivity less than 3,000. The "Middle" region has lowest relative connectivity, it is less than 10% connected. Similar trends can be seen for the synthetic storm shown in Figure A7, suggesting the primary driver will most likely be the location of the electrojet (as noted in Beggan [2015]). The level of network connectivity tracks geomagnetic activity and matches well to the overall variation of wavelet power in Figure 2. Recall for a node to be included in the network it must have wavelet power, C j k ,  , of more than its 80th quantile over the event. This means at the peak in the South more than 50% of nodes have a rate of change in the top 20% of the event, as well as a high correlation with other nodes.
The interconnections between regions informs how GIC risk is potentially transferred across the high voltage power system; though we emphasize that connections within the network imply coherent response but not necessarily direct flow between nodes. This is plotted in the lower three panels in Figure 7, where the y-axis is approximately half that of the top three panels. During the most intense activity, the north-middle ORR ET AL.  and north-south regions are about 10% connected indicating a lower risk that GIC from a geomagnetic storm located in the northern part of the island would affect the southern regions. The relative connectivity between the middle and south regions rises to 20% however, indicating a greater coherence in response across these regions as a percentage of the total number of nodes. The total number of connections differ in each zone but relate to the topology of the grid in each region (e.g., from population density and generator location). Figure 8 explores the topological properties of the networks at the same times (06:13, 14:25, 21:23, and 00:23 UT) as the network maps shown in Figure 6. The variation of network connectivity with threshold is explored in the lower panels of Figure 8 which plots the relative (left ordinate) and absolute (right ordinate) connectivity as a function of wavelet cross-correlation threshold for the times of the snapshots of the time series shown in the previous figures (Figures 5 and 6).
The panels share the same ordinate axis scaled for relative connectivity. The snapshots consistently show that within each region, the number of edges does not vary strongly with cross-correlation threshold so that these sub-networks/spatial regions are roughly uniformly coherent. Between the regions we also see weak variation with threshold at the onset of the storm when we see widespread long-range connections across GB in Figure 6, again suggesting the same (high) level of coherence across a broad spatial region. In later snapshots, the number of edges between the regions decreases as we increase the cross-correlation threshold. This is consistent with a few long-range connections which are highly coherent.
The upper panels of Figure 8 plot the degree distributions. In detail, these are highly dynamic, but overall they confirm that the networks are not fully connected. These degree distributions for connections within each region are highly peaked around a single characteristic degree; they do not simply follow that of the physical HV grid which has been shown (Rosas-Casals et al., 2007) to have a much broader degree distribution which is approximately exponential, with a characteristic degree (typical number of connections per node) of less than 5. The physical HV grid is also "small world," that is, dominated by local, nearest-neighbor connections (Rosas-Casals et al., 2007). Instead we can see from the degree distributions for connections between regions in Figure 8, and in the network maps (Figures 6, A4, and A6), a significant fraction of connections are not between neighboring physical nodes in the HV grid, instead they represent long-range spatial coherence between a subset of nodes. Highly peaked degree distributions are consistent with a cluster, or subset of nodes that are almost all connected to each other so that (almost all) m connected nodes in the cluster have the same number  1 m connections.

Discussion
The goal of this study is to provide insight into the hazard posed by large rates of change of magnetic field during storms and the manner in which GIC flows through a connected power grid, particularly a complex and highly developed one. Wavelet analysis is used to provide a time-frequency localization of both the rapidly changing magnetic field and the GIC response; showing how the wavelet power is distributed spectrally and spatially from the magnetic field variation and the resulting GIC. As described in Appendix A1, the Haar wavelet provides a complete orthogonal decomposition of the time series so that each wavelet coefficient captures a fraction of the total power in the signal (see Equations A4 and A5). As first described for Figure 2, for signals with   1 t minute resolution, each wavelet scale j resolves the rate of change in the signal at time resolution    2 j j t t minutes. The spectral energy density is concentrated in the shortest time periods (responding to a rate of change of magnetic field on the fastest resolvable 2 min timescale, captured by the  1 j Haar wavelet) which is often linked to the rapid equatorward motion of the aurora. Measurements show the combined impact of the Earth and inductive power network is to act as a bandpass filter for which GICs respond best to short period variation (for example, see Figure 3 of Hübert et al. [2020]). Using the wavelet decomposition in combination with network analysis provides a new tool to understand both the connections between magnetic field and GIC, and how the GIC response can become coherent across the HV grid.
In simple models of power grids, co-variation of GIC in nodes can be easily observed, but more complex and realistic models for HV grid networks do not easily allow these deductions to be made. The GB high voltage grid is an amalgam of several individual networks developed over decades. They were formally brought ORR ET AL.

10.1029/2021SW002772
16 of 30 under central control in the mid-2000s. The grid thus offers a large amount of redundancy and alternative pathways for GIC to flow making it difficult to interpret the interaction of topology and the geophysical drivers. We first considered the response of a realistic model for the HV grid to a simplified synthetic GIC model. The modeled electrojet was fully spatially coherent and varying uniformly and relatively slowly in time. To quantify this response, we calculated the wavelet cross-correlation between all pairs of nodes on the HV grid. In our analysis, wavelet cross-correlation above a threshold between a pair of nodes corresponds to these nodes being connected in a network. Our wavelet-based cross-correlation network then quantifies spatial coherence between nodes on the HV grid induced by the electrojet, rather than the properties of the physical HV grid. When driven by a fully spatially coherent model electrojet, the network is fully connected. The spatial extent of the network during much of the synthetic storm shows that there is a general tendency for connections to establish in discrete, fully connected clusters or "islands" (Figure 4) located below the electrojet, within which each network node is connected to every other. This is in contrast to the properties of the physical HV grid which are small world and exponentially distributed, so that connections to each node are dominated by those to their nearest physical neighbors.
We then replaced the simple model electrojet driver with that sampled from observations of real events. This more realistic driving differs from the model electrojet in two important ways (a) it is not necessarily spatially coherent and (b) it can be highly time-variable. We note that due to the limitations of the thin-sheet model used that the magnitude of the geoelectric field (and hence GIC) may not match the true values experienced during the modeled storms though the regional variation is realistic. In the October 2003 storm, the network connections and maps in Figures 6 and 7 are not fully connected and show a variety of spatial structure. At storm onset, and intermittently throughout the storm, there is a prevalence of long-range connections spanning the entire physical grid, with up to 5% of nodes becoming connected between north, middle and south regions. The initial part of the storm is associated with a sudden storm commencement (SSC) which is a depression of the magnetopause boundary (due to the increased speed and density of the CME at it arrives). This creates a rapid and global increase in current flow on the sun-side of the planet which leads to a strong increase in the horizontal magnitude of the external field. This is the reason for the initial coherence across the region. After the SSC, the magnetosphere loads up with energy which then drives substorms and the auroral electrojet at high latitudes which are locally intense thus creating smaller networks of connections.
Quiescent periods are interspersed with short (<10 min) intervals where there are connections forming distinct sub-networks, with a tendency for east-west linkages which follow long-line features in the physical structure of the HV grid. These have strongly peaked degree distributions, consistent with the coherent response of a subset of network nodes which can be spatially separated but all become highly connected to each other, embedded in a much larger system with weak or no connections. The GB grid layout typically reflects the population centers of Britain and their connection to power generation which often run along a latitude axis. The predisposition for east-west connections may also relate to the induction of an east-west geoelectric field by the north-south component of the magnetic field, particularly if there is a large rate of change due to motion of the auroral oval equatorwards, though there are other storm-time related phenomena that also generate sizable GICs (e.g., Belakhovsky et al., 2019).
As can be seen in Figure 7 the network is actually very dynamic during the periods of high geomagnetic variation. The time dynamics of the network response does not simply follow the overall power level in the driving rate of change of magnetic field (that is, the wavelet coefficient amplitude), and we do not find a one-to-one correspondence between instantaneous power and number of connections in the network. The connections can partition into distinct and spatially extended "islands" or clusters during the more intense phases of the storms but with gross overall connections occurring during strong pulses of auroral motion. Hence, this confirms that GIC hazard in any particular location is not directly linked to trans-island flow but remains beholden to the location and spatial extent of the auroral electrojet in conjunction with the ORR ET AL.  Figure 7; corresponding to the times: 06:13, 14:25, 21:23, and 00:23 UT. The histograms in the top panels show the degree distribution at each time within and between latitude bands. The top two panels of histograms are for a threshold of 70% and the following two for a threshold of 90%. Plotted below are the normalized number of edges in a latitude band (y-axis) where the network is calculated with increasingly high correlation threshold abscissa. underlying conductivity of the region. This correlation analysis does suggest that strategic monitoring of GIC at selected well connected nodes can provide additional insight into the expected variations at other local transformers.

Conclusions
We used wavelets to examine, in a time and frequency localized manner, the spatial and temporal extent of rapid magnetic field changes over the European sector during large geomagnetic storms based on minute mean data from SuperMAG. The wavelet decomposition maps the time series into a time-localized distribution of power in specific frequency bands, resolving the short time intervals when there are large rates of change of the magnetic field during periods of equatorward motion of the auroral oval.
We apply the same wavelet decomposition to modeled GICs in the GB power grid for examples of large storms. We can then combine this with network analysis based on wavelet cross-correlation to understand the flow of GIC in the high voltage grid. To construct a network response we calculated the wavelet cross-correlation between all pairs of nodes on the HV grid. In our analysis, wavelet cross-correlation above a threshold between a pair of nodes corresponds to these nodes being connected in a network. Our wavelet-based cross-correlation network then quantifies spatial coherence between nodes on the HV grid induced by the electrojet, rather than the properties of the physical HV grid.
We applied this methodology to characterize the HV grid response to a synthetic storm. The synthetic storm has a fully spatially coherent electrojet which drives a fully spatially coherent response below the synthetic electrojet, that is, a fully connected sub-network. We then study the HV grid response to several real events, most notably the October 2003 storm. We find that: 1. Coherent networks form with long-range connections that can span the entire physical grid at most intense times. 2. These networks are highly intermittent, each configuration is of a few tens of minutes duration, interspersed with quiescent intervals. 3. The networks are not fully connected which is distinct from the response seen for a fully spatially correlated synthetic electrojet. 4. The networks are not a simple response to the amplitude of the rate of change of magnetic field as resolved by the amplitude of the wavelet coefficients. There is no one-to-one relationship between number of connections in the network and overall amplitude of rate of change of magnetic field. 5. At storm onset, and intermittently throughout the storm, networks form where there is a prevalence of long-range connections spanning the entire physical grid, with up to 5% of nodes becoming connected between north, middle and south regions of the GB. 6. Quiescent periods are interspersed with short (<10 min) intervals where there are connections forming distinct sub-networks, with a tendency for east-west linkages which follow long-line east-west features in the physical structure of the HV grid which will also be aligned with the maximum driving due to a rate of change of north-south magnetic field by motion of the auroral oval. 7. These networks have strongly peaked degree distributions, consistent with the coherent response of a subset of network nodes forming a spatially extended "island" or cluster in which the nodes all become connected to each other, embedded in a much larger system with weak or no connections. 8. The highly peaked network degree distributions are distinct from that of the physical HV network which is exponential and "small world," with many connections between nodes which are in close physical proximity to each other.
The network connections between nodes indicate a time-localized coherent response to the electrojet and as such can be interpreted as correlations between transformers and a proxy for GIC hazard in locations where no direct monitoring is made. The GIC response networks that we have determined here have significantly different properties to that of the physical HV grid. This is important since it implies that previous studies that focus on stability of the physical grid to the failure of individual network connections may not fully inform the assessment of space weather risk. From this study we can infer where further monitoring for GIC either directly or indirectly should be provided in order to improve real-time assessments or forecasts of where issues may arise during large geomagnetic storms.
Here, all the , j k C are obtained over the same finite time window  N t (N even). As we move from one scale j, to the next,  1 j , the time resolution of is a single value which for the Haar wavelet is just the N point average of the signal. Equation A2 is a particular choice of normalizing prefactor for the wavelet transform. This prefactor is not unique, it defines that required for the inverse transform so that the transform pair satisfies Parseval's theorem. Practically, Equation A2 is not evaluated directly, instead, a recursive "filter bank" approach is used to construct the , j k C . We will use the Haar wavelet which is complete and orthogonal. The sum of the average energy contained in each wavelet coefficient is just the average energy contained in the original signal so that using Equation A4 we have: Finally, we then have the rate of change of magnetic field estimated at each wavelet scale which is simply and as here, the signal x is magnetic field (nT ) so that Equation A5 has physical units 1 nTs . Note that the total energy of the signal, and the total energy contained in each wavelet scale, will increase with N, and the fraction will decrease correspondingly with N. The Haar wavelet is a step function, with mother wavelet: and under operation Equation A1 effectively "convolves" the original signal with a ramp of width  2 j t.
ORR ET AL.
10.1029/2021SW002772 20 of 30 An additional wavelet coefficient we will use for our network construction (described below in Appendix A2) is the maximal overlap discrete wavelet transform (MODWT) (Percival & Mofjeld, 1997;Percival & Walden, 2000). The MODWT is a non-decimated modification of the DWT which is defined for the full sample size N and can be used for multi-resolution analysis of geophysical processes (Percival & Guttorp, 1994;Percival & Mofjeld, 1997). The N wavelet coefficients (MODWT) can be thought of as circularly filtering the original time series ( ) x t using a rescaled mother wavelet,     1 ( ) 2 j t , with no subsampling (Percival & Mofjeld, 1997). MODWT wavelet coefficients are defined as for   1, , t N,   m t l modulo N and     (2 1)( 1) 1 j j L L when L is the width of the wavelet (Percival & Walden, 2000). If the sample length N is an integer multiple of the scale length 2 j , DWT wavelet coefficients, , j k C Equation A2 can be recovered by rescaling a subsample of every 2 j th MODWT coefficients,  , j t W (Percival & Guttorp, 1994).