Average Ionospheric Electric Field Morphologies During Geomagnetic Storm Phases

We utilize principal component analysis to identify and quantify the primary electric potential morphologies during geomagnetic storms. Ordering data from the Super Dual Auroral Radar Network (SuperDARN) by geomagnetic storm phase, we are able to discern changes that occur in association with the development of the storm phases. Along with information on the size of the patterns, the first six eigenvectors provide over ∼80% of the variability in the morphology, providing us with a robust analysis tool to quantify the main changes in the patterns. Studying the first six eigenvectors and their eigenvalues shows that the primary changes in the morphologies with respect to storm phase are the convection potential enhancing and the dayside throat rotating from pointing toward the early afternoon sector to being more sunward aligned during the main phase of the storm. We find that the ionospheric electric potential increases through the main phase and then decreases once the storm phase begins. The dayside convection throat points toward the afternoon sector before the main phase and then as the potential increases throughout the main phase, the dayside throat rotates toward magnetic noon. Furthermore, we find that a two‐cell convection pattern is dominant throughout and that the dusk cell is overall stronger than the dawn cell.

WG19 examined the general trends in the SuperDARN data during geomagnetic storms, such as latitudinal expansion of the ionospheric convection maps, data coverage, data availability, and cross-polar cap potential (i.e., convection strength), in relation to solar wind and geomagnetic conditions. The study also compared statistically the responses of these measured parameters during geomagnetic storm phases, to periods of disturbed geomagnetic activity, irrespective of storm phase, as well as high solar wind driving when no storms occurred. One of the primary results of this paper was that the storm phases, as well as the ionospheric responses measured by SuperDARN are closely tied to the solar wind driving of the system, which matches previous results (e.g., Gillies et al., 2011;Loewe & Prölss, 1997): During the main phase of a geomagnetic storm, higher solar wind driving due to southward interplanetary magnetic field (negative B Z ) enhances the current systems connecting the ionosphere with the magnetosphere. We thus see a higher cross polar cap potential, as well as an enhanced Sym-H index, matching our understanding of how the system works (e.g., Milan et al., 2017). WG19 showed that throughout a geomagnetic storm, there is some asymmetry in the two-cell convection pattern measured by SuperDARN, with the dusk cell being much stronger than the dawn cell, as well as changes throughout the storms in the location where the fastest flows are measured in the ionosphere: This is primarily on the dayside, though in the initial and recovery phase the fastest flows are primarily measured in the noon to early morning sectors whereas during the main phase of a storm, this is shifted toward the afternoon sectors. WG19 also found that the return flow boundary (the latitudinal location where antisunward flows neighbor the sunward flows) and the Heppner-Maynard boundary (HMB) (Heppner & Maynard, 1987) (the boundary where the high-latitude ionospheric convection pattern terminates) move throughout the storm phases, as does the latitudinal distance between them.
Other previous studies using SuperDARN data from geomagnetic storm periods have looked at the number of scattering echoes and line-of-sight velocities in relation to sudden storm commencements and sudden commencements (e.g., Gillies et al., 2012;Kane & Makarevich, 2010), but without a detailed quantitative analysis of ionospheric convection morphologies. A further statistical study by Gabrielse et al. (2019) compared the mesoscale flows measured by SuperDARN during the main phases and recovery phases, as well as coronal mass ejection (CME) and high-speed stream (HSS) storms. Whilst WG19 did not split the data into the exact same categories, the results broadly agree with these previous studies. Here, we only focus on the geomagnetic storm phases to learn about the average ionospheric behavior. Whilst WG19 answers some basic questions on the morphology and latitudinal extent of ionospheric convection during the phases of a geomagnetic storm, we will examine the morphologies of geomagnetic storms in more detail here. Here, we will study these data further to answer the following question: How do ionospheric convection morphologies change throughout the storm phases?
We answer this question by utilizing an objective method for dimensionality reduction (principal component analysis; e.g., Joliffe, 2002), which will tell us what the primary morphologies in the data are with respect to storm phase.

Data
There are two primary data sets used in this study: the geomagnetic storm list and the SuperDARN data, which we describe in this section.

Geomagnetic Storms
The geomagnetic storm list is published by WG19 and can be found in their supplementary material. It is formed by applying an automatic identification algorithm to the Sym-H index, which reflects enhancements in the global ring current (Iyemori, 1990). The algorithm identifies the initial, main, and recovery phases of geomagnetic storms, similar to Hutchinson et al. (2011), which allows us to draw conclusions WALACH ET AL.

10.1029/2020JA028512
2 of 16 about the phenomena associated with the progression of storms. In brief, the initial phase of a geomagnetic storm is classified by a positive excursion in the Sym-H index, associated with an increase in the Ferraro-Chapman currents along the magnetopause, followed by a decrease to below -80 nT during the main phase, where the ring current enhances. The minimum in Sym-H coincides with the end of the main phase, which is followed by a gradual increase to normal values, known as the recovery phase. For further detail, the reader is referred to WG19.
With the expansion of the SuperDARN network to mid-latitudes, we are able to study the dynamics of the high-to-mid-latitude ionospheric convection with unprecedented coverage (Nishitani et al., 2019). One of the findings by WG19 was that the high-latitude convection maps, which can be produced with SuperDARN data, can expand to 40° of geomagnetic latitude during disturbed times, which was not accounted for in previous versions of the SuperDARN Radar Software Toolkit (RST versions <4.2 (SuperDARN Data Analysis Working Group et al., 2018)), which had a cutoff of 50° magnetic latitude. The finding of this expansion matches magnetometer and spacecraft measurements from previous studies (e.g., Kikuchi et al., 2008;Wilson et al., 2001).
The SuperDARN data used here were therefore processed using the RST, which is specifically designed to accommodate SuperDARN observations down to 40° of magnetic latitude. Typically, to make SuperDARN convection maps several steps of processing have to be followed: (1) Using RST, an autocorrelation function is fitted to the raw radar data. This produces fitacf files, which store the line-of-sight velocity data. (2) The data is then gridded onto an equal area latitude-longitude grid (see Equation 1 from Ruohoniemi and Baker [1998]) and combined into 2-min cadence records. (3) Data from different radars are combined and the spherical harmonic fitting algorithm is performed which fits an electrostatic potential in terms of spherical harmonic functions to the data (Ruohoniemi & Baker, 1998;Ruohoniemi & Greenwald, 1996). When this fitting is performed, typically, a background model, parameterized by solar wind conditions is used, to infill information in the case of data gaps (e.g., Thomas & Shepherd, 2018). Alongside this, a HMB (Heppner & Maynard, 1987), the low-latitude boundary of the convection pattern where the flows approach zero, can either be specified or be chosen using the data. This is to constrain the convection pattern when the spherical harmonic fit is applied (Shepherd & Ruohoniemi, 2000). For typical 2-min convection maps, it is appropriate to use the data to find a threshold of three radar velocity measurements of greater than 100 ms −1 for the HMB (Imber et al., 2013).
For the purpose of this study, we make 2-min cadence superposed epoch convection maps, where data from the different storms are combined. This differs slightly from the usual steps outlined above and is explained further in the following section.
We utilize the same storm list and the same gridded SuperDARN data, spanning from 2010 to 2016, as published in WG19. We have 54 storms with a median duration for each storm phase of 19.5 h for the initial phase, 9.1 h for the main phase, and 55.8 h for the recovery phase.

Method
To study the characteristic ionospheric convection morphologies of the storms in detail, we make a superposed epoch analysis. Similarly to Hutchinson et al. (2011) and Wharton et al. (2020), we make a superposed epoch analysis of the storms which treats each storm phase independently. The duration of each storm phase is scaled to the median duration, such that the initial, main, and recovery phases are 19.5-, 9.1-, and 55.8-h long, respectively. We can thus compare average characteristics across storms.

10.1029/2020JA028512
3 of 16 We apply our method to the SuperDARN data to make average storm convection maps, which are parameterized by storm phase and median duration: We use the gridded data from the previous study (WG19), and write new convection maps for each storm phase, which are thus time-normalized and comprise the data from all storms. To make the convection maps, we write files with all the data and run the map-fitting procedure using RST v4.2 (SuperDARN Data Analysis Working Group et al., 2018) and an eighth order spherical harmonic expansion (Ruohoniemi & Greenwald, 1996). This differs slightly from the usual method described earlier: To make the storm maps, no statistical background model was used, as the data coverage is very good when combining data from 7 years of geomagnetic storms. As data coverage at lower latitudes can be sparse, especially during the initial phase, the automatic HMB algorithm can select unrealistic boundaries. We avoid this by forcing the HMB to match the lower quartile of the distribution of HMBs from the individual maps per timestep per phase (this is shown in Figure 8 in WG19 and the second panel from the top in Figure 4 in this study). To minimize unphysical artifacts dominating the dayside potential, we add padding below the HMB on the dayside by adding artificial datapoints with line-of-sight velocities which are equal to zero. We also set all line-of-sight velocities to zero for any backscatter points on the dayside which lie below the HMB. Before fitting the spherical harmonic expansion, we also merge the line-of-sight data, using the MERGE technique (Cerisier & Senior, 1994;André et al., 1999). This resolves all measurements at a given grid point into one vector. It is worth noting that despite the padding and merging of vectors, the fitted electrostatic potentials are not forced to be zero below the HMB (due to the fitting process using a spherical harmonic expansion) and as such, the convection cells do sometimes extend across the HMB.

Intermediary Maps
Examples of these average convection maps are given in Figure 1, which shows a map from the beginning of each storm phase. All other maps are included in the form of animations as supplementary material or can be downloaded as convection map files from Lancaster University's research archive (PURE) (Walach, 2021).
From Figure 1, we see that the convection patterns are different at the beginning of each storm phase: As expected, at the beginning of the initial phase, the convection pattern is relatively small and the ionospheric convection velocities are low, whereas at the beginning of the main phase, the familiar two-cell convection pattern (e.g., Ruohoniemi & Greenwald, 1996) is enhanced and expanded, with fast return flows seen on the dusk-side. From examining these convection maps (see also supporting information), we see that the twocell pattern stays strong and expanded throughout the main phase. Figure 1 and the supporting information shows that this is further enhanced at the beginning of the recovery phase. We see from the supporting information that the fast flows and expanded pattern stays prevalent long into the recovery phase, but start to decrease after the main phase ends.

Principal Component Analysis
Studying these average maps is useful to observe obvious changes in the convection, such as deviations from the two-cell convection regime, expansions, and contractions, or patches of fast flows. To quantify changes in the convection morphologies further, we now utilize principal component analysis on the data. This is a well-known technique for pattern recognition and is also known under different names, such as empirical orthogonal functions, and has been used successfully for geophysical data sets (see Baker et al., 2003, Cousins et al., 2013, 2015Shore et al., 2018;Shi et al., 2020, Kim et al., 2012, andreferences therein). An alternative method is to use the spherical harmonics to examine changes (e.g., Grocott et al., 2012), but in this case, the components are predetermined, which limits their interpretability. In PCA, the components are defined by the data which allows us to find the main constituents which make up the patterns. Overall, this allows us to quantify the main components of the patterns and see how they change over time.
The underlying principle is that the data set can be decomposed into a series of basis functions that reveal underlying correlations within the data. In our case, the data set is made of the electrostatic potential maps, Φ t (where t = 0, …, m), such that m = 2,532 (the median storm duration at a time resolution of 2 min) and each Φ t has n-elements, where n is given by the number of latitude by longitude grid points (2° resolution).

10.1029/2020JA028512 4 of 16
All the observations can be expressed as one m × n matrix (Φ). The covariance matrix Σ is then given by The data Φ t can be expressed (or reconstructed) in terms of eigenvectors, X i , of the covariance matrix Σ and their components, α i , such that This means components at a given time, α i , are given by Applying this method to the convection maps allows us to quantify and detect morphological changes automatically, as well as determine the primary components which make up the ionospheric electric field. To do this, we first scale all the ionospheric convection maps, such that they are the same size. This is necessary for the principal component analysis to work. Using different pattern sizes would involve padding areas with no data with zeros and result with no correlation between the majority of gridpoints, and thus the principal component analysis method would not work. Whilst changing the size of the pattern will make the expansions and contractions invisible for the Principal Component Analysis, this information is kept, so it can be studied in conjunction with the components later. We discuss this again later in the paper and also address the expansions and contractions in WG19. We take the electrostatic potential from each map and resize the potential pattern by scaling by the HMB (Heppner & Maynard, 1987) at midnight to 40° of magnetic latitude. We map the potential to a 2° latitude by 2° longitude grid which allows us to describe each pattern by a one-dimensional 4,500 line matrix (n = 4,500). We then calculate the mean for all storm epochs at each spatial point in the electric potential grid and subtract this from each individual map. On the remaining data set, we perform the eigendecomposition using the Householder method of eigendecomposition (Press et al., 2007). Using only data from geomagnetic storm times for the principle component analysis means that the only bias is in our event selection, which was done using the automatic algorithm from WG19 on the Sym-H index. It is worth noting that whilst selecting by geomagnetic storm times only means we can analyze the storm-time morphologies specifically, we also impose a selection bias: although we include some quieter times during the recovery phase of the storms, this selection bias results in our mean and eigenvector patterns looking different from analyses done in previous studies (e.g., Cousins et al. [2013] used an interval which had very little geomagnetic activity and Milan [2015] used all of the available AMPERE data) and we comment on this further in section 5 (Discussion).

Results
By examining the eigenvalues, we can determine the importance of each of the eigenvectors (i.e., the component patterns that are added or subtracted together to make the convection maps). Figure 2 shows the cumulative explained variance, expressed in percentages. We see immediately that the curve converges fast: The orange dotted and dashed lines show the i-values closest to >80% and >95% cutoff values, respectively. Whilst we have 4,500 eigenvalues and vectors, we see from Figure 2 that we do not need all these values to express the majority of the variability in the electric potential patterns. In fact, the variance converges fast enough that the first six eigenvectors explain over 80% of the variance (this is shown by the green lines).
In the following parts of the manuscript we will thus focus our attention on the first six eigenvectors and components and examine these further.
By adding or subtracting factors of X i (where i = 1, …, 4,500) we are able to thus reconstruct the initial maps. These factors as a function of time are given by the components, α i . i i i f , which represents the approximate CPCP each component holds and we can thus analyze this with respect to time through the storm phase. We now examine these terms (i = 0…6) in more detail. Figure 3 shows the primary electrostatic potential pattern components: the panel in the top left corner shows the mean pattern which was subtracted from all maps before applying the principal component analysis. The other panels show the first six eigenvectors (i.e., the most dominant pattern components). The pattern components  * 1, ,6 X are normalized by their CPCP, such that the color scale approximately represent a range of 1. We will refer to this same normalization factor, f i , again later, as it will aid the interpretation of Figure 4. Each panel shows the eigenvector as a map in the same coordinate system as Figure 1, whereby the magnetic pole is in the center, noon is toward the top of the page, and dusk toward the left. The concentric dashed circles outline equal latitudes at 10° separation. As expected, the mean shows that a clear two-cell electric potential is dominant, with an enhancement in the dusk cell. What is less expected is that we also see an anti-clockwise rotation of the pattern about the pole. We see that * 1 X is able to provide an increase or decrease in the two-cell convection potential by adding or subtracting the asymmetry from the mean pattern due to the similar rotation about the pole in the convection throat.
* 2 X provides morphological asymmetry by being an almost uniformly negative potential, so adding or subtracting this would strengthen one cell and weaken the other, or vice-versa. X provide a motion toward earlier or later local times of the throat and other asymmetries, such as a variation to potential in the center of the pattern.
The top panel of Figure 4 shows a superposed epoch analysis of the interplanetary magnetic field components, B IMF , resolved into the GSM (Geocentric Solar Magnetospheric) coordinates with X in light green, Y in turquoise, and Z in dark blue. The second panel from the top shows the HMB (in black) which the maps were scaled by, as well as the number of backscatter points per average SuperDARN map (in rose). This is followed by the median Sym-H and the CPCP (in yellow). Then we show the first six components of the eigenvectors, all as a function of storm phase-adjusted time, which are shown in gray. The black lines show the low pass filtered curve, using a 60-min centered kernel window to show the large scale changes more clearly. The first vertical dashed blue line marks the end of the initial phase and thus the beginning of the main phase and the second dashed blue line shows the end of the main phase and the beginning of the recovery phase.
We observe that the B Z component is clearly enhanced, especially during the main phase of the storm and that the number of backscatter points per SuperDARN map is high (this can also be seen from the animations MS01-MS03 in the supporting information). for this timestep (see also Equations 1 and 2). The benefit of scaling α i by f i (i.e., the true range of X i ), is that the scaled components  * i approximately represent the CPCP of each pattern and thus aids interpretation.
We see immediately that much of the variability in the components is dominated by what appears to be noise, which we will investigate more quantitatively in the next section. Focusing on the black curves we see a few clear changes in  * i with respect to the geomagnetic storm phases:  * 1 shows a clear change which mirrors the HMB and Sym-H closely. At the start of the main phase, this value decreases abruptly, then stays negative and then starts to increase gradually throughout the recovery phase.  * 2 also decreases as we approach the end of the main phase but then increases quickly into the first part of the recovery phase. During the recovery phase,  * 2 then fluctuates about zero from about 10 normalized hours onwards but remains primarily positive. This is distinctly different to  * 1 which continues to increase throughout the recovery phase.  * 3 is primarily negative throughout the initial phase, then increases to a positive value through the main phase and remains primarily positive throughout the recovery phase.  * 4 to  * 6 remain very small and show no clear deviations from zero with respect to the storm phases.
To analyze these changes further with respect to IMF B Y and B Z and Sym-H, we perform a cross-correlation analysis between each of these parameters and the components. To highlight the variations over larger timescales, we use the smoothed components from Figure 4. The best correlation coefficient, |r|, of each of these and their respective lag times, t, are given in Table 1. We also show p for each correlation pair, which is defined as the significance of the correlation. This is defined by Press et al. (2007) as where erfc is the complementary error function and N is the number of datapoints, which is, as defined earlier, m, the number of 2-min maps. This value expresses the probability that in the null hypothesis of two values being uncorrelated, |r| should be larger than its observed value. A small value of p (i.e., p = 0) thus indicates that the correlation is significant. Table 1 shows that p is practically zero for all the correlations. This means these correlations are statistically significant. We see that the first component in particular is highly correlated with both Sym-H and B Z , with a time lag, t = 0. This means that changes in this component are correlated with changes in Sym-H (i.e., the storm phases) and B Z (i.e., solar wind driving). As i increases, |r| tends to decrease. The correlational pairs with B Y are in general lower than the correlations with B Z , which means the time variability we see in the components tend to correlate better with B Z than B Y . The notable exceptions here are  * 3 , and  * 5 , which are the only components where the correlations with B Y are marginally higher than the correlations with B Z .
The time lags are more difficult to interpret but indicate several patterns: The majority of the convection pattern (i.e.,  * 1 , which holds almost 50% of the variance) shows its best correlation at t = 0, which means this component's contribution is mostly related to Sym-H as this is how the storm phases are defined. We further note, that for any pairs where |r| is very low (<0.3), t tends to be >1 h, which we interpret to not be meaningful and thus do not comment further on these.

Discussion
In Figure 4 we show that, on average, the HMB expands to <50° magnetic latitude approaching the main phase and stays expanded, well into the recovery phase when considering the lower quartile of the distribution shown in WG19. It is possible that in reality, this expansion moves to lower latitudes than 40° for individual storms, but our observations are limited by the geographical location of the SuperDARN radars and our choice of the HMB. This expansion is coincident with the IMF B Z component becoming more southward, leading to a higher dayside reconnection rate and thus more rapid opening of magnetic flux WALACH ET AL.   (Cowley & Lockwood, 1992;Milan et al., 2012;Siscoe & Huang, 1985;Walach et al., 2017). This means an expansion of the open-closed field line boundary occurs, which happens in tandem with the expansion of the convection pattern observed here (see also WG19). The high-latitude ionospheric electric field and thus convection pattern is an important mechanism for plasma transport, and thus its expansion will mean the circulation of plasma at lower latitudes than was previously circulated by the high-latitude convection pattern. Zou et al. (2013) also showed that the convection pattern expanding during geomagnetic storms plays an important role in the generation and propagation of SEDs seen on the dayside at mid-latitudes: Zou et al. (2013) found that there are two parts to SEDs, with the equatorward expansion of the convection pattern being the primary driver for the SED formation.
We find that the first six eigenvalues hold >80% of the variability in the scaled ionospheric electric potential during storms (see Figure 2). As the potential patterns which are analyzed using the principal component analysis are scaled by the HMB, this variability does not include the expansion or contraction of the pattern, which happens in addition to the morphological changes analyzed here. The first and second eigenvectors (see * 1 X and * 2 X in Figure 3) represent a dual-cell convection pattern, associated with the Dungey-cycle (e.g., Dungey, 1961Dungey, , 1963Milan, 2015;Walach et al., 2017); when  * 1 and  * 2 are negative * 1 X and * 2 X are subtracted from the mean, producing a more enhanced dual-cell convection pattern. We see from Figure 4 that this is the case throughout the main phase of the storm, subsiding in the recovery phase and peaking toward the end of the main phase, when solar wind driving is highest. This matches the findings of WG19, which showed that this is also when the cross-polar cap potential is highest. We see from Figure 4 that the CPCP addition from the first component changes from ∼20 kV in the initial phase to ∼−40 during the main phase, which is a step change of 60 kV and slightly higher than the 40 kV step change in CPCP that was seen in WG19. This highlights that whilst this component drives a lot of the storm phase change-related variability, more components need to be added to get an accurate representation of the CPCP. The second component also adds to the potential, in particular during the main phase, where its contribution reaches ∼20 kV. The first component primarily enhances or decreases the duskside of the potentials, whereas the second component primarily enhances or decreases the dawnside potential cell. During the main phase of the storm, when they are both negative, the convection pattern is enhanced and the two cells both increase. A few hours into the recovery phase, however, when  * 1 is still negative and  * 2 is positive (both are at ∼20 kV magnitude), the electric potential increases on the dusk side but decreases on the dawn side, which means the dusk cell is noticeably larger than the dawn cell. We see from Figure 4 that the following components contain slightly lower magnitudes of the potential and decrease with each component.
The third, fourth, fifth, and sixth components only add up to ∼10 kV to the convection pattern at their peak, which is minimal in the context of a CPCP between 50 and 120 kV. It is confirmed by Table 1 that what looks like noise in Figure 4 in some of the higher-order components ( * 4 and  * 5 ), is indeed very weakly correlated with Sym-H, which means these changes are not related to the storm phases. Whilst  * 6 shows a higher correlation (|r| = 0.427), it adds however less to the total CPCP and is thus less important. We see that the correlation between  * 1 and Sym-H is, on the other hand, very high (|r| = 0.706) and significant (p = 0.00), which means this component is clearly correlated with the storm phases. This component is also highly correlated with B Z , which is no surprise, given the high levels of solar wind driving seen during geomagnetic storms.

r and p Values Between Sym-H; B Y ; B Z and Each Component Shown in 4 (Black Smoothed Lines)
The third eigenvector ( * 3 X ) resembles the classic dual cell convection pattern but with a 90° rotation about the pole toward dawn. This component is, therefore, able to add asymmetry to the dual cell pattern in an unconventional way: its addition can move the dayside throat to earlier local times. The fourth and fifth eigenvectors ( * 4 X and * 5 X in Figure 3) represent asymmetric dawn-dusk changes to the patterns, which appear to mainly rotate the convection throat on the dayside, though can rotate the nightside convection throat as well. The sixth eigenvector ( * 6 X ) is very symmetrical and closely resembles the second-order and degree spherical harmonic pattern (e.g., see Figure 2 from Grocott et al. [2012]).
We see from Figures 3 and 4 that the main changes with respect to storm phase are primarily related to the dawn and dusk cells enhancing or decreasing. Figure 4 shows that the third component is primarily negative during the initial phase. Then, going into the main phase of the storm, the third component increases steadily until a change in polarity is seen in this component, right before the end of the main phase. This will not only change the cross-polar cap potential, increasing it during main and recovery phases and decreasing it during the initial phase, but it will also change the location of the dayside throat. It indicates that the convection throat on the dayside reaches across the midnight-noon meridian toward dawn and becomes more noon aligned as the main phase progresses but then jumps back to be more dusk-aligned before the end of the recovery phase. For the rest of the storm time, we see this component varying slightly between positive and negative values, but primarily staying positive, meaning that the dayside throat has a tendency to be noon-aligned.
This may appear to be a result of solar wind driving and a change in the IMF B Y component, which can move the dayside convection throat (e.g., Cowley & Lockwood, 1992;Thomas & Shepherd, 2018). This would be further evidenced as  * 4 shows a mild correlation (0.426) with the IMF B Y component, but this component adds a minor amount of electric potential and  * 3 is much less correlated with B Y (0.228) than Sym-H (0.440). We see however from the top panel in Figure 4 that the average IMF B Y component is near zero for these storms. In fact, 37% of the time the IMF B Y component is positive for these storms, 38% of the time the IMFB Y component is negative, and it is zero the rest of the time. We see that it is the IMF B Z component, which is enhanced during the main phase of the storm. That the average storm does thus not have a strong dusk-dawn component modulating the dayside flows (i.e., neither positive B Y , nor negative B Y are consistently dominant) is also shown in Figure 2 (panel j) in WG19, which shows that during the main phase of the storm, the IMF is overwhelmingly southward for all storms considered here. Usually, when SuperDARN maps are created, base-models, which are in part parameterized by the solar wind are used (e.g., Thomas & Shepherd, 2018) such that datagaps are overcome. In this study, however, no solar wind inputs were used at all as the data coverage is very good when combining data from 7 years of geomagnetic storms. We conclude that some of this rotation in the dayside throat may be due to an IMF B Y component, but we speculate that there are other mechanisms at play due to the inconsistency in the directionality of the B Y component.
We theorize that some of the control in the dayside throat moving toward later local times could be due to a number of factors (or combination thereof): higher solar wind driving and the dayside reconnection rate increasing, or due to feedback through other means (e.g., thermospheric winds  and/or SEDs modulating the location of the throat (Zou et al., 2013(Zou et al., , 2014 and/or the plasmaspheric plume impacting the magnetopause reconnection rate post-noon). Further evidence for the plasmaspheric plume being responsible for this moving of the dayside convection throat is available from comparing our results to those of Wharton et al. (2020): In their study, Wharton et al. (2020) looked at the eigenfrequencies in-ground magnetometer variations on the dayside during the same storm phases as ours. They found that at L-shells <4, the eigenfrequencies in magnetometer measurements increase during the main phase of geomagnetic storms, which is due to the decrease in the plasma mass density caused by plasmaspheric erosion. This approximately corresponds to a geomagnetic latitude of 60° or less (see Table 1 in Wharton et al. [2020]), which corresponds to the dayside throat location we see during the main phase of the storm. Wharton et al. (2020) find that at L > 4 (which maps to higher latitudes and thus inside the convection pattern on the dayside), the eigenfrequencies decrease by ∼50% during the main phase, due to a weaker magnetic field and an enhanced plasma mass density. This may be further evidence of the plasmaspheric plume. Overall, however, to find a conclusive answer for the moving of the dayside throat, further studies are needed. WALACH ET AL. Morphological changes on the nightside are more difficult to analyze and less likely to yield great insight due to the time-averaging that we have done: We know (see Table S1 in WG19) that the minimum and maximum durations of each storm phase can vary vastly (e.g., the recovery phase can be anything from ∼6 to ∼163 h). By combining the data, such that the average convection maps match the median storm phases, we time-shift the data. Whilst the majority of storms are of similar length, it provides a good framework for studying the average storm-time responses; however, other time-dependent phenomena, such as substorms, are averaged out. It is well known that substorms frequently occur during geomagnetic storms and are important for the energization of the ring current (e.g., Daglis, 2006;Sandhu et al., 2019), but Grocott et al. (2009) showed that substorms primarily produce a response in the high-latitude ionospheric convection pattern on the nightside and that ordering by onset location is important when trying to gain insight from the average convection pattern. It thus follows that although substorms commonly occur during geomagnetic storms, we do not see their signatures. We, therefore, cannot say if there is any substorm ordering by storm phase or time throughout the storm phases as no clear substorm signatures are seen in the average maps. Gillies et al. (2011) studied line-of-sight SuperDARN velocity measurements during geomagnetic storms and found that an increase in IMF B Z is accompanied by a speed increase measured with SuperDARN in the noon sector (9-15 MLT) and midnight Sector (21-3 MLT) during the main phase. Gillies et al. (2011) also found a reduction in the measured plasma drift early in the main phase for intense storms, and speculated this either to be due to a reduction in the plasma drift speed or a change in the direction of the drift relative to the SuperDARN radar beam. In this study, we have shown (see Figure 4) that the addition to the convection potential increases during this time (due to the first, second, and third components), which means that the convection potential increases and thus ionospheric convection velocities are likely to be also increasing. This is supported by our previous analysis (WG19) which showed that the cross-polar cap potential increases during this time, and thus the convection should also increase. This provides further evidence that the decrease in the plasma drifts seen by Gillies et al. (2011) during the main phase is due to the change in the direction of the flows relative to the SuperDARN radar beam (i.e., the second of their two theories).
Whilst SuperDARN measures convective flows, from which we can infer the ionospheric electric field, it is useful to compare the patterns to measurements of the Field-Aligned Current system (FACs). We expect the convection electric field to be related to the FACs in strength, but also in location . The FACs close in the ionosphere, such that the location of their footprints is expected to match the HMB (R2 currents) and the flow reversal boundary (R1). We thus expect the strength of the FACs to increase when the electric field increases and vice versa. Similarly, when the convection pattern expands, the location of the FACs' footprints is expected to expand also. Cousins et al. (2015) and Shi et al. (2020) used empirical orthogonal function analysis to describe the modes of the field-aligned currents. Shi et al. (2020) split the data according to different solar wind drivers, including HSS and transient flows related to CMEs, both of which can be drivers of geomagnetic storms. Their patterns reflect the prevalence of the dual cell electrostatic pattern that we also see, but due to different data binning, their modes are different, making a direct comparison difficult. Overall, Shi et al. (2020) found that Sym-H is highly correlated with the modes in the transient flow category, indicating that strong geomagnetic storm activity dominates this category, which gives a strong dual cell convection pattern, as well as expansions and contractions. Both their HSS and transient categories show a mode that gives a strong asymmetry on the dayside (and would result in a similar movement of the dayside throat that we see), which are highly correlated with Sym-H activity, but also the IMF B Y and B X components, and AE and solar wind temperature. Whilst the data presented by Cousins et al. (2015) did not contain any considerable geomagnetic storm activity, their results generally agree with the results from Shi et al. (2020). What does stand out when comparing results, however, is that their first mode shows, similar to Shi et al. (2020), a strengthening of the pattern, which is highly correlated with AE and the IMF B Z component. This is followed by a mode describing the expansions and contractions, which is correlated with B Y , AE, and Sym-H. The third mode from Cousins et al. (2015) describes the cusp shaping, which is also correlated with B Y , AE, and tilt, but not Sym-H. It is worth noting that Cousins et al. (2015) only showed the first few modes, and their chosen time period contains little geomagnetic activity. Cousins et al. (2013), on the other hand, used the EOF analysis to study SuperDARN data. They analyzed 20 months of plasma drift data to study electric field variability WALACH ET AL.

10.1029/2020JA028512
12 of 16 and found that the first component accounted for ∼50% of the observed total squared electric field (which is as a proxy for the electrostatic energy per unit volume) and is primarily responsible for variations on long timescales (∼1 h). It is worth noting that their components look different to ours as they used a different data set (i.e., their K p median was 1, so they used a non-storm time data set) for input but in general find the two-cell convection pattern to be dominant as well.
Comparison between our data, Shi et al. (2020), Cousins et al. (2013), and Cousins et al. (2015) shows that using different data brings out different modes with different properties: the primary EOF in Cousins et al. (2015) strengthens the convection pattern, whereas the second component has a shaping function, followed by expanding and rotating modes. They further find that their top correlation for the first component is at 0.44 for the AE index, which is considerably lower than our top correlation (0.706) coefficient between Sym-H and the first component. The dayside throat in the patterns (mean and components) shown by Cousins et al. (2013) show no movement: their mean is perfectly aligned with noon, which we attribute to the fact that their input data is on average from both positive and negative B Y with no storm effects. Conversely, the mean pattern from Milan (2015), where they applied the principal component analysis to a much larger data set of the Birkeland currents inferred by AMPERE, showed the throat aligned with 11 and 23 MLT. This is comparable to the average conditions, also when studying SuperDARN data (e.g., Thomas & Shepherd, 2018), and indicates that the mean and the components are sensitive to the choice of input data.
As part of this study, we have provided a first analysis of how the dayside throat responds to geomagnetic storms (i.e., internal magnetospheric dynamics), versus IMF B Y conditions (i.e., external magnetospheric dynamics) and studied the timescales of dayside throat changes with respect to geomagnetic storms. To understand this fully requires further study. If the dayside throat is rotated due to the plasmaspheric plume mechanism, we would expect to see the same movement in the throat (away from dusk) in the southern hemisphere, but we would expect to see it moving in the opposite sense in the southern hemisphere for any IMF B Y related effect. We have provided a first-order analysis of this and discussed potential mechanisms here, but to find a more definitive answer, southern hemisphere data will be investigated in more detail in a future study.

Summary
We have utilized SuperDARN line-of-sight ionospheric plasma measurements to study ionospheric electric potential morphologies during geomagnetic storm time and specifically geomagnetic storm phases. We applied a principal component analysis to average ionospheric convection maps to examine the primary morphological features for the first time and using eigenvalue decomposition, we see how dominant patterns change over time (i.e., through the storm phases). The main dynamics in the morphologies that we have uncovered are happening to the ionospheric electric potential pattern on a large scale: the electric potential pattern expands and contracts; the potentials increase and decrease in strength; and the dayside convection throat rotates. We speculate that all these changes are due to the IMF B Z component of the solar wind increasing during the main phase of the storm.