Determining ULF Wave Contributions to Geomagnetically Induced Currents: The Important Role of Sampling Rate

Past studies found that large‐amplitude geomagnetically induced current (GIC) related to magnetospheric Ultra Low Frequency (ULF) waves tend to be associated with periods >120 s at magnetic latitudes >60°, with comparatively (a) smaller GIC amplitudes at lower latitudes and shorter wave periods and (b) fewer reports of waves associated with GIC at lower latitudes. ULF wave periods generally decrease with decreasing latitude; thus, we examine whether these trends might be due, in part, to the undersampling of ULF wave fields in commonly available measurements with 60 s sampling intervals. We use geomagnetic field (B), geoelectric field (E), and GIC measurements with 0.5–10 s sampling intervals during the 29–31 October 2003 geomagnetic storm to show that waves with periods <∼120 s were present during times with the largest amplitude E and GIC variations. These waves contributed to roughly half the maximum E and GIC values, including during times with the maximum GIC values reported over a 14‐year monitoring interval in New Zealand. The undersampling of wave periods <120 s in 60 s measurements can preclude identification of the cause of the GIC during some time intervals. These results indicate (a) ULF waves with periods ≤120 s are an important contributor to large amplitude GIC variations, (b) the use of 0.1–1.0 Hz sampling rates reveals their contributions to B, E, and GIC, and (c) these waves' contributions are likely strongest at magnetic latitudes <60° where ULF waves often have periods <120 s.

Several studies have explored the impact of B sampling rate on estimates of extreme E values and GIC modeling, with mixed results. Some studies found that sampling intervals of 60 s can adequately capture B/E/GIC variations (e.g., A. , while others (e.g., Grawe & Makela, 2021;Grawe et al., 2018;Trichtchenko, 2021) found that 60 s measurements are not adequate due to the presence of variations with frequencies (f) higher than the Nyquist frequency (f Nyquist60 = 0.0083 Hz for measurements sampled every 60 s). Rogers et al. (2021) found for three locations in the United Kingdom that B fluctuations with roughly 1,200 s period (f ≪ f Nyquist60 ) produced the most intense E when considering events occurring several times per year, whereas the 1-in-100 years return levels were greatest for 30-120 s period fluctuations (f ≥ f Nyquist60 ). Trichtchenko (2021) found that a 0.1 Hz sampling rate is generally needed to understand power system responses and posed the question, "What should be the sampling rate (or Nyquist frequency) in order to provide an adequate representation of the fast geomagnetic variations (such as storm sudden commencement (SSC), substorms, pulsations, and rapid variations during the main phase of the storm) for use in the calculations of the extreme geoelectric field values and in the GIC modeling?" Motivated by this question, the present study will focus on the recommended sampling rate for one particular source of B listed by Trichtchenko (2021), pulsations, considering several factors: expected pulsation frequencies, the amplitude of B/E related to pulsations at different frequencies, and the ability of pulsations to cause potentially hazardous GIC. Pulsations and other sources of fast geomagnetic variations (e.g., SSCs, substorms) are not mutually exclusive and can occur at the same time. We return to this point in later sections, particularly Section 3.
Pulsations are also known as Ultra Low Frequency (ULF) waves (McPherron, 2005;Orr, 1973), with ULF being a frequency designation for magnetospheric plasma waves ranging from 0.001 to 5 Hz (Jacobs et al., 1964). Hereafter, we only use the term "ULF waves." At the lower end of the ULF band, waves are usually well-described by an ideal magnetohydrodynamic (MHD) approximation (Southwood & Hughes, 1983). The bounded, inhomogeneous medium of the Earth's magnetosphere supports a variety of MHD modes analogous to standing waves: 1. Standing Alfvén waves: These are analogous to waves on a string where the natural frequency of the wave depends on the length of the magnetic field line and wave speed (Sugiura & Wilson, 1964). The frequencies of these waves vary with magnetic latitude, generally increasing in frequency as magnetic latitude decreases. Standing Alfvén waves can be observed at magnetic latitudes corresponding to closed magnetic field lines and in all geomagnetic conditions. The fundamental frequency varies depending on several factors including geomagnetic latitude, geomagnetic activity, and local time (Takahashi & Anderson, 1992;Wild et al., 2005). 2. Trapped/Partially trapped magnetosonic waves: Magnetosonic waves can become partly radially trapped between various boundaries in the magnetosphere, with their frequencies related to the wave speeds and the size of the region in which they are trapped. These wave modes include cavity modes, waveguide modes, and virtual resonances (Lee & Takahashi, 2006;Wright & Mann, 2006). They can be observed at a wide range of magnetic latitudes during all geomagnetic conditions (Keiling & Takahashi, 2011). The wave period does not change with latitude, and the same period can be observed at a wide range of latitudes and local times (e.g., M. D. Hartinger et al., 2017;Pilipenko et al., 2010;Shi et al., 2017). The plasmasphere region has frequently been linked to wave trapping, with theory, modeling, and observation confirming the typical fundamental plasmaspheric cavity mode/virtual resonance frequency is in the ∼0.007-0.022 Hz range (M. D. Hartinger et al., 2017;Keiling & Takahashi, 2011;Shi et al., 2017;Takahashi et al., 2018); these plasmaspheric waves are most often observed at magnetic latitudes ∼<60°.
These two broad standing wave categories can occur simultaneously. For example, Takahashi et al. (2018) showed observations and numerical simulations of a cavity mode driven by a solar wind pressure pulse that, in turn, feeds energy into standing Alfvén waves via a process known as field line resonance (Glassmeier et al., 1999;Tamao, 1965). Figure 1 shows typical frequencies for these different wave modes based on past observations, including during the 29-31 October 2003 storm (Kale et al., 2009;Marin et al., 2014;Takasaki et al., 2006) and similar conditions such as the aftermath of interplanetary shocks/SSC (Takahashi et al., 2018). As seen in Figure 1, both plasmaspheric virtual resonance/cavity modes and standing Alfvén waves are undersampled in studies relying on measurements sampled every 60 s at magnetic latitudes <∼60°, as their typical frequencies are >∼f Nyquist60 .
Past studies have linked ULF waves to B (see above), E (e.g., M. D. Hartinger et al., 2020;Heyns et al., 2021;, and GIC variations (e.g., Belakhovsky et al., 2019;Viljanen et al., 2006;Yagova et al., 2021). These studies have generally found that the largest ULF wave B and GIC amplitudes are observed at magnetic latitudes >∼60° for waves with f < f Nyquist60 (Pilipenko, 2021). Several studies  (Figure 3 of that study, YOR-CRK and NUR-OUJ pairs-comparable to a separate study at a different local time by Chi et al. (2005)) and Takasaki et al. (2006) (Figure 2 of that study); the bars are for the range of observed values during the 3-day interval. The green line is for an observed global waveguide mode on 31 October 2003 (Marin et al., 2014); the blue line is for an observed and simulated plasmaspheric cavity mode during an interplanetary shock event on 15 August 2015 (Takahashi et al., 2018); and the cyan shaded region is for a typical range of plasmaspheric cavity mode/virtual resonances based on review papers (Keiling & Takahashi, 2011) and wave transit time calculations (Takahashi et al., 2003, Figure 11). Dashed lines indicate the Nyquist frequencies for data sampled at 10 and 60 s, while black circles on the bottom of the plot indicate the magnetic latitudes where most observations were collected for this study (see Table 1; the KAK station is at 28.99° and not shown here). have been conducted at magnetic latitudes >∼60° using measurements from the Kola peninsula in Russia (e.g., Apatenkov et al., 2020;Belakhovsky et al., 2019;Yagova et al., 2021), where waves with f < f Nyquist60 drive the largest amplitude GIC, in some cases as large as ∼100 A. Additional studies have examined GIC in gas pipelines at magnetic latitudes near 60°, with the largest amplitude GIC often originating from waves with f < f Nyquist60 (Pulkkinen & Kataoka, 2006). Based on these past observations, in particular the tendency for larger wave amplitudes and larger GIC for f < f Nyquist60 , a few studies have recently suggested that waves with periods of 120-600 s drive the most intense GIC (Pilipenko, 2021;Yagova et al., 2021).
Reasons to further explore whether this tendency for smaller B and GIC for waves with f > f Nyquist60 generalizes to magnetic latitudes <60° include: 1. There are limited B measurements with sampling intervals less than 60 s, particularly during historical geomagnetic storms or long monitoring intervals (Trichtchenko, 2021). Thus, extreme ULF wave B amplitudes for wave f > f Nyquist60 are not as well constrained as for f < f Nyquist60 . 2. As shown in Figure 1, several categories of ULF waves are undersampled when using 60 s measurements at magnetic latitudes <60°. Since their frequency content is concentrated in a narrow frequency band, they will be removed from time series that have been low-pass filtered and resampled to obtain 60 s measurements.
In contrast, waves with f < f Nyquist60 and other phenomena with a broadband frequency spectrum are fully or partially retained in 60 s data. 3. Case studies, such as the 24 March 1991 SSC, show the tendency for larger amplitude B and E variations at f < f Nyquist60 does not always apply to magnetic latitudes <60°. During this event, waves with f > ∼f Nyquist60 were observed globally with B amplitudes of ∼100-200 nT (Yumoto et al., 1994;Araki et al., 1997, Figure 5a). Concerning GIC,  presented measurements following an SSC of damped, sinusoidal GIC at multiple locations in the United States power grid (Western Virginia, New York) with initial amplitudes ranging from 66 to 130 A (Figures 2 and 3 in that study) and a wave f > ∼f Nyquist60 . According to , this SSC, "produced some of the largest GICs ever measured in the United States at midlatitude locations" (page 2 of ).
Motivated by these facts and past observations, we examine additional B, E, and GIC measurements at magnetic latitudes <60° and with sampling intervals less than 60 s during the 29-31 October 2003 geomagnetic storm (Section 2). We examine these measurements as well as high-pass filtered measurements (f > ∼f Nyquist60 ), low-pass filtered measurements (f < f Nyquist60 ), and measurements with 60 s sampling intervals, comparing maximum and minimum values in the time rate of change of B, E, and GIC. We further determine whether ULF waves are present during periods with maximum/minimum E and GIC values using standard visual identification criteria (Jacobs et al., 1964). We show with examples that waves with f > ∼f Nyquist60 drive or contribute significantly to B, E, and GIC variations with amplitudes comparable to the maximum values detected during extended monitoring intervals in Japan and New Zealand. In Section 3, we use these results to justify a recommended sampling rate of 0.1-1.0 Hz for ULF waves. Finally, we summarize our results in Section 4.

ULF B, E, and GIC During the 29-31 October 2003 Storm
The 29-31 October 2003 geomagnetic storm, also known as the Halloween storm, was driven by two coronal mass ejections (CME). The first CME arrived at Earth on 29 October 2003 just after 06:10 UT, and the other arrived on 30 October at 16:00 UT; each resulted in an SSC (Balch et al., 2004). During the storm, there were numerous reports of GIC in power grids (Heyns et al., 2021;Kappenman, 2005;Ngwira et al., 2008;A. Pulkkinen et al., 2005;Rodger et al., 2017), gas pipelines , and oil pipelines (Hejda & Bochníček, 2005 Potapov et al. (2006), Pilipenko et al. (2010), and Marin et al. (2014) examined wave activity during several phases of the storm, finding several time intervals with ∼0.003 Hz wave amplitudes of ∼100-300 nT at a wide range of magnetic latitudes and longitudes. Wave amplitudes were typically found to peak in the magnetic latitude range from 58° to 65°, and the waves were associated with both standing Alfvén waves and global waveguide modes (Pilipenko et al., 2010). Love and Gannon (2010) observed similar period waves at low latitude stations during several intervals during the storm.
In this section, we present B (2.2), E (2.3), and GIC (2.4) measurements associated with ULF waves at magnetic latitudes <60° for several time intervals during the 29-31 October 2003 storm, placing them in context with extreme values obtained during extended monitoring intervals in Japan and New Zealand.

Data Sets and Methodology
We use B, E, and GIC measurements from several regions around the world with measurements available at sampling rates ≥0.1 Hz. Station codes, coordinates, measurement type, and sampling rate are shown in Table 1 (see Figure S1 in Supporting Information S1 for a map). The geographic coordinates provided are in degrees north and degrees east. The geomagnetic (quasi-dipole) latitudes provided were obtained using the International Geomagnetic Reference Field (IGRF) model and calculator available from the British Geological Survey (Emmert et al., 2010) assuming a fractional year of 2003.8.
For the United States locations, we use B measurements from the United States Geological Survey (USGS) Geomagnetism Program sampled at 1.0 Hz (Love & Finn, 2011;Rigler & USGS Geomagnetism Program, 2023) at the Fredericksburg (FRD) and Shumagin (SHU) sites as well as B and E measurements from the University of California, Berkeley, Parkfield (PKD) monitoring site sampled at 1.0 Hz (Kappler et al., 2010).
Outside of the United States, we analyze B (Minamoto, 2013) and E (Fujii et al., 2015) measurements from the Kakioka Magnetic Observatory (KAK) station in Japan. These data are sampled at 1.0 Hz during the period of interest. Fujii et al. (2015) noted that unique ground conductivity structure near KAK tends to amplify the eastwest component of E significantly relative to the north-south component at KAK, and relative to the east-west component at other nearby stations. Thus, we will not emphasize the amplitude of E variations at KAK relative to other locations, but rather the amplitudes observed at KAK during the 29-31 October 2003 storm relative to other intervals.
We also analyze measurements from two locations in Europe near southern Sweden, where GIC and power system disruptions were reported (e.g., A. Pulkkinen et al., 2005). In Sweden, we investigate B measurements from the IMAGE network (Tanskanen, 2009); specifically, B sampled at 0.1 Hz from the Uppsala (UPS) observatory in Sweden operated by the Geological Survey of Sweden (SGU). In Germany, we use B and E measurements from the German Research Centre for Geosciences (GFZ)-Potsdam Niemegk (NGK) observatory sampled at 2 (B) and 1.0 Hz (E).
We further investigate B measurements from the South Island of New Zealand Eyrewell (EYR) observatory due to its long record of measurements sampled at >0.2 Hz and the availability of nearby GIC measurements. This site is operated by GNS Science. The highest sampling rate available from EYR has changed over time  These procedures do not completely remove slowly varying sensor responses related to daily variations in air temperature and other factors unrelated to magnetosphere-ionosphere current systems (Fujii et al., 2015) which may be present in some of the measurements we analyze, but none of these factors affect our analysis or conclusions in this study; future work could remove these trends to study a wider frequency range (e.g., daily variations related to the solar quiet [Sq] current system). To show the relative contributions of B, E, and GIC variations with frequencies above and below f Nyquist60 , we apply digital high-pass and low-pass filters with 0.0083 Hz cutoff frequencies to these data in some figures. In particular, we use a forward-backward eighth order Butterworth filter constructed using the publicly available Python SciPy signal module ("butter" and "sosfiltfilt"). These filters are constructed to have a flat response in the pass band at frequencies near f Nyquist60 , thus providing more reliable estimates for maximum/minimum values during time intervals when variations have frequencies close to f Nyquist60 . The sum of the high-pass and low-pass filtered measurements also yields the original time series.
We use 60 s measurements to show how commonly available data that have frequently been used to study the 29-31 October 2003 storm and other historical storms differ from the 0.5-10 s samples obtained for this study. The 60 s data from FRD, EYR, KAK, UPS, and NGK are obtained from the INTERMAGNET database (Love & Chulliat, 2013;St. Louis, 2020); data from each site are processed using one of several different filter types and resampling methods to obtain the INTERMAGNET samples. For example, FRD and UPS sites use a Gaussian filter, while the KAK and NGK sites use a boxcar average. In some cases, the samples are centered at second 00 of the minute (e.g., KAK, UPS), whereas in other cases they are centered at 30 s (e.g., NGK). SHU and PKD are not available from INTERMAGNET for these dates, thus 60 s data are obtained for these sites via a 60 s boxcar average similar to that used for the KAK site and by SuperMAG (Gjerloev, 2012). The fact that different sites listed in Table 1 use different filters and resampling methods to obtain 60 s time series does not change our results, as all of the 60 s data sets are affected by the undersampling issues we discuss throughout this study. However, these different methods can result in different amplitudes for variations with frequencies near f Nyquist60 , and we briefly discuss these effects in later sections.

Geomagnetic Field Measurements
In this section, we review B measurements from several locations (see Table 1 and Figure S1 in Supporting Information S1) at magnetic latitudes <60°, including two time intervals during the 29-31 October 2003 storm that correspond to periods with ULF waves and with the largest ISL M6 GIC values during 14-year monitoring intervals. Figure 2 is for the horizontal magnitude of B, or H. In all panels, the black curves are for 0.5-10 s measurements (see Table 1 for sampling rate at each location); light gray curves are for 60 s measurements obtained from either INTERMAGNET (FRD, EYR, KAK, UPS, and NGK) or by applying a centered boxcar average (PKD, SHU) as described in Section 2.1; blue curves are for high-pass filtered measurements (f > f Nyquist60 ); and red curves are for low-pass filtered measurements (f < f Nyquist60 ). To emphasize amplitude differences between the low-and high-pass filtered time series at each site rather than differences in amplitude between different sites, we use different y-axis ranges for each panel; see Figure S2 in Supporting Information S1 for data presented with the same y-axis range. Figures 2a and 2b show a step-like increase (SSC) starting just after 06:11 UT in H (black curves) that is mostly captured in low-pass filtered measurements at FRD and PKD (red curves). There are smaller amplitude H variations with f > f Nyquist60 (blue curves) that are not present in the 60 s samples (compare black and gray lines). Figure 2c is for a station (SHU) at a higher magnetic latitude and closer to local midnight.  Table 1; panel (a) is for Fredericksburg, (b) for Parkfield, (c) for Shumagin, (d) for Eyrewell, (e) for Kakioka, (f) for Uppsala, and (g) for Niemegk. In all panels, black curves are for measurements sampled at 0.1-2 Hz (see Table 1 for specific sampling rate at each site), blue curves are for high-pass filtered measurements (f > 0.0083 Hz), and red curves are for low-pass filtered measurements (f < 0.0083 Hz). Light gray curves are for measurements sampled at 0.016 Hz obtained from either INTERMAGNET (all panels except (b, c)) or from a centered 60 s boxcar average (panels (b, c)).
Here, wave activity is more evident in the original 1 s samples and high-pass filtered measurements; several wave cycles can be seen in the black and blue curves. Figures 2d and 2e for the EYR station in New Zealand and KAK station in Japan are similar to SHU, including a step-like increase along with H variations with f > ∼f Nyquist60 (blue line) that are not captured in 60 s data from INTERMAGNET. Finally, Figures 2f and 2g are for the UPS (Sweden) and NGK (Germany) stations; the step-like increase is not visible, and H primarily has variations with f > f Nyquist60 . Note that the UPS and NGK stations bracket the locations where disruptions occurred in the Swedish power grid at 06:11 and 06:12 UT on 29 October 2003 (A. Pulkkinen et al., 2005), and the UPS site is at a similar magnetic latitude and longitude to a Finnish gas pipeline where GIC variations were observed with similar periodicities near this time , black curves in their Figure 6).
In Figure 3, we place these H measurements in context with Halloween storm measurements from later time periods at EYR, examining the entire storm as well as shorter intervals corresponding to periods with intense GIC at the nearby ISL site (Section 2.4). Figures 3a and 3b show EYR H measurements over the 3 day period from 29-31 October 2003. Figure 3a shows measurements sampled at 0.2 Hz (black line) and the same data low-pass filtered (red line, f < f Nyquist60 ), while Figure 3b Figure 3d. ULF wave activity is seen in both panels. Note that maximum and minimum values differ across panels because different baseline values were removed (corresponding to the mean value for each time range) to make it easier to compare variations; the range of variation is the same in each panel. These shorter time intervals correspond to the storm periods with largest GIC amplitudes (described in Section 2.4) and largest f > f Nyquist60 wave amplitudes (blue curves in Figure 3b), but not to the periods with the largest overall magnetic disturbance amplitude (Figure 3a). Some of the f > f Nyquist60 variations shown in Figure 3 are caused by disturbances with a broadband frequency spectrum rather than ULF waves. More generally, broadband disturbances such as the step-like increase in H caused by the SSC can occur at the same or nearly the same time as ULF waves and make significant contributions to H, E, and GIC, making it difficult to quantify the relative contributions from ULF waves and other disturbances with frequency content above f Nyquist60 . We return to these points in Section 3. Figure 4 is for a widely used proxy for |E|, , hereafter referred to as H′; this proxy for EYR |E| has been shown to correlate well with nearby ISL GIC measurements . Figure 4a is for H′ measurements from EYR sampled at 5 s (black) and subsequently low-pass filtered (red, f < f Nyquist60 ) during the 29-31 October 2003 storm interval. Figure 4b is for high-pass filtered 5 s H′ measurements (blue, f > f Nyquist60 ). Comparing  Figures 4a and 4b, much of the frequency content in H′ is contained in f > f Nyquist60 , particularly during periods with the largest H′. This is seen, for example, by noting that the maximum and minimum values of the black line in Figure 4a are often much larger than the corresponding values in the red line (f < f Nyquist60 ). More direct, quantitative evidence that much of the frequency content in H′ is contained in f > f Nyquist60 can be seen in power spectra and integrated power comparisons for f < f Nyquist60 and f > f Nyquist60 in Figures S3 and S4 in Supporting Information S1. For example, Figure S4b in Supporting Information S1 shows that the power contained in frequencies above f Nyquist60 is greater than or equal to power from frequencies below the Nyquist frequency during much of the storm; the ratio of higher frequency power to lower frequency power is usually greater than 1. In each panel, 5 s (black), low-pass filtered (red), and high-pass filtered (blue) H′ observations are shown; a gray line corresponding to H′ calculated with INTERMAGNET 60 s measurements is also shown for reference. In Figure 4c, H′ variations with f > f Nyquist60 make up the largest contribution to the peak value of 583 nT/min at 06:12 UT (note the correspondence between the blue and black lines at that time) and also have roughly equal contributions with variations having f < f Nyquist60 during a sustained period with ±100 nT/min amplitude ULF wave activity with f ∼ f Nyquist60 starting well after the SSC at 06:40 UT (compare red and blue lines from 06:40 to 07:00 UT). Figure 4d shows similar results for the second interval. Here, variations with a mix of frequencies are present (consistent with power spectra results shown in Figures S3 and S4 in Supporting Information S1), with sinusoidal H′ variations with f > f Nyquist60 evident from 02:00 to 02:30 UT producing maximum/ minimum values in H′ of roughly ±500 nT/min. The largest contributions to H′ during these two intervals are not present in the INTERMAGNET 0.016 Hz (1 min) data. Table 2 summarizes maximum and minimum values for H′ across the entire 3-day interval and for the two shorter intervals shown in Figures 3 and 4; recall that the two shorter intervals correspond to times with the largest ISL M6 GIC amplitudes. Across the entire storm, the maximum and minimum values of H′ are 583 and −521 nT/ min for the 5 s sampling interval, 236 and −198 nT/min for low-pass filtered 5 s data, and 166 and −171 nT/min for INTERMAGNET 60 s data. Substantial differences are also found in the two shorter intervals. For example, during Interval 2 on 31 October 2003 from 02:00 to 03:00 UT, H′ has a range of 993 nT/min in the original 5 s measurements and 173 nT/min in the 60 s INTERMAGNET measurements. Further results for H′ are listed in Table S1 in Supporting Information S1, including a comparison of values when using different resampling methods. Table 2 of Rodger et al. (2017) showed that the maximum H′ recorded in the 5 s EYR data during geomagnetically disturbed periods between 2001 and 2015 was 583 nT/min; this occurred on 29 October 2003 at 06:11 UT, the same period shown in Figure 2. Thus, the H′ measurements examined in Figures 2-4 have values equal or comparable to the maximum H′ recorded at EYR during geomagnetically disturbed periods across that 14-year interval.

Geoelectric Field Measurements
We next review E measurements at a few of the locations discussed in Table 1 with available data. Figures 5a and 5b are for the x and y component (x is for geographic north and y is for geographic east) of E at the PKD station in the United States. As before, the black curves are for the original E measurement sampled at 1.0 Hz, and red and blue curves are for low-pass and high-pass filtered 1.0 Hz measurements (cutoff frequency at f Nyquist60 ). Figures 5c and 5d are for E measurements at the NGK station in Germany; the primary contributor to NGK E variations is from f > f Nyquist60 (variations in the black curves mostly match the blue curves). Figures 5e and 5f are for KAK E measurements in Japan. At both KAK and NGK, the wave amplitudes in the high-pass filtered measurements often exceed the wave amplitudes in the low-pass filtered measurements (compare red and blue curves in the bottom four panels). Panels (e, f) are in the same format as Panels (a, b), but for NGK and KAK, respectively. Figure 6 is for the geographic y component of E at KAK for the entire three day-interval, presented in a similar format as Figures 3 and 4 for EYR magnetic field measurements (the x component of E at KAK is show in Figure  S5 in Supporting Information S1). Figure 6a is for 1 s (black) and subsequently low-pass filtered (red, f < f Nyquist60 ) measurements, while Figure 6b is for high-pass filtered (blue, f > f Nyquist60 ) measurements. The low-and highpass filtered time series have comparable amplitudes during many portions of the storm (compare red curve in Figure 6a to blue curve in Figure 6b), particularly during the periods with the largest E y amplitudes. This is also consistent with dynamic power spectra and integrated power shown in Figure S4 in Supporting Information S1; for example, Figure S4d in Supporting Information S1 shows that while integrated power at f > f Nyquist60 is usually smaller than integrated power at f < f Nyquist60 during the entire 3-day interval, the ratio of the two quantities is usually close to 1 at times when the power (E y amplitude) is largest. Figures 6c and 6d for the shorter time intervals on 29 October 2003 from 06:00 to 07:00 UT and 31 October 2003 from 02:00 to 03:00 UT also show this. Figure 6c shows that the higher frequencies continue to make significant contributions to the overall E y variations well after the SSC (note the comparable amplitudes of red and blue curves). Figure 6d shows that in the later period, lower frequencies have larger amplitudes, but the higher frequencies have amplitudes within a factor of two or three. The maximum value of E y of 0.644 V/km at KAK shown in Figure 5 (black line, see also Table S2 in Supporting Information S1 for |E| which is 0.658 V/km) is the largest value measured at that location during the 3-day interval from 29-31 October 2003 when considering measurements sampled at 1.0 Hz. Unfortunately, there are no long-term statistics available for the 1 s E measurements at KAK to benchmark this value, but we can compare with long-term statistics based on longer sampling intervals: 1. Fujii et al. (2015) considered E measurements with 3,600 s (1 hr) sampling interval during an 11-year monitoring interval from 1 January 2000 to 28 February 2011 (see their Figure 2), finding maximum E variation amplitudes more than a factor of two lower than 0.644 V/km. 2. Zhang and Ebihara (2022) (2022)); this value was exceeded in the 29 October 2003 storm when considering 1 s data ( Figure 5, Table S2 in Supporting Information S1 of this study, 0.658 V/km).
As shown in these examples and Table 2 and Table S2 in Supporting Information S1, the 1 s data can produce roughly a factor of two or larger peak values of E when compared to 3,600 s (1 hr) and 60 s (1 min) data. The smaller peak values in 3,600 and 60 s data are due in part to the undersampling or removal of ULF waves with f > f Nyquist60 expected at KAK's magnetic latitude (Figure 1). The SSC observed at KAK on 29 October 2003 was large but not the largest ever recorded (e.g., Araki et al., 1997); thus it is somewhat surprising that the 100-year peak value of SSC-related E predicted by Zhang and Ebihara (2022) (when using 60 s data) was realized for the 29 October 2003 SSC (when using 1 s data). This result is at least partly related to ULF waves that are concurrent with SSC (Araki et al., 1997;Saito & Matsushita, 1967) being (a) undersampled in 60 s measurements and (b) averaged out during superposed epoch analysis due to wave phase/frequency varying from event to event.

Geomagnetically Induced Current Measurements
We next examine ISL M6 GIC measurements from New Zealand (located very close to EYR [ Table 1]). Figure 7a shows H′ sampled at 0.2 Hz at EYR and is provided for reference to compare against GIC measurements. Figure 7b is for ISL M6 GIC measured throughout the 29-31 October 2003 storm. The largest amplitude GIC tends to occur at the same time as large values of H′ at EYR, consistent with the finding from Rodger et al. (2017) that these two measurements are well-correlated. Figure 7c shows GIC measurements during the 29 October 2003 SSC. As in Figure 6, the black line is for the 4 s measurements; the red for low-pass filtered measurements (f < f Nyquist60 ); the gray dotted line for low-pass filtered measurements resampled to 60 s centered at second 00 of the minute; and the blue for high-pass filtered measurements (f > f Nyquist60 ). The largest GIC magnitude of 34.1 A during a 14-year monitoring interval ( Table 1 in Rodger et al. (2017)) corresponds to the period with the spike and subsequent variations in H′ seen in Figure 7a. Starting from that initial increase, GIC variations with f ∼ f Nyquist60 continue with decaying intensity until approximately 06:40 UT, when they begin increasing in intensity again with amplitudes of ∼10 A. This trend of decaying then increasing variation amplitude is also seen in H′ (Figure 4c, see also Figure S6 in Supporting Information S1). For much of the interval, frequencies above and below f Nyquist60 have comparable GIC amplitudes, as evidenced by the comparable amplitudes of the red and blue curves.  A, Figure 7c), the range of GIC values is similar in both events, including the rapid change from −17.7 to +23.2 A (40.9 A range, Figure 7d). Similar to H′ and E y results in Figure 4 and Figures S3 and S4 in Supporting Information S1, there are multiple frequencies present in the GIC measurements shown in Figure 7d (variations are present in both the red and blue curves); those with f > ∼f Nyquist60 make significant contributions throughout the interval, including a large negative GIC that is not captured in the low-pass filtered measurements. The GIC correlate best with the 5 s H′ in Figures 4c and 4d (black curves) rather than the 60 s H′ from INTERMAGNET (gray curves) or the low-pass filtered 5 s H′ measurement. The bottom two rows of Table 2 further compare maximum and minimum GIC values before and after low-pass filtering during this interval. In the original measurements, the maximum and minimum values are 23.2 and −17.7 A, whereas for f < f Nyquist60 they are 18.3 and −8.1 A. The range of GIC variation is a factor of 1.55 larger in the original measurements, (40.9 A compared to 26.4 A), again indicating that the GIC has significant frequency content above f Nyquist60 . Rodger et al. (2017) reported statistical results pertaining to ISL M6 GIC measurements from 2001 to 2015. Their Table 2 shows that the maximum GIC at this location occurred during the 29 October 2003 SSC at 06:11 UT (34.1 A, interval shown in Figure 4c of this study). All entries in their Table 2 correspond to periods when the EYR 60 s H′ exceeds 40 nT/min, and the range of GIC intensities for the other 30 events in their list is 4.6-33.1 A. Thus, to the extent that any source of ISL M6 GIC during the 14-year monitoring interval could be regarded as significant (i.e., produce large amplitude GIC), the waves with frequencies >∼f Nyquist60 are significant. They can generate GIC variation amplitudes of roughly 10-20 A (blue line in Figure 7d), comparable to the maximum value of 34.1 A over the 14-year monitoring interval. They also likely contribute to the peak value of 34.1 A itself, as they may occur at the same time as the magnetopause current intensification that gives rise to the step-like increase in B during the SSC (interval 1 in Figure 3, see further discussion in Section 3). These waves can also be regarded as potentially hazardous, as these GIC intensities are comparable to the ISL M6 measurement of 33.1 A when a transformer failure occurred farther south in Dunedin/Halfway Bush, New Zealand, on 6 November 2001 (Marshall et al., 2012). No GIC measurements were available at the location of the transformer that failed, but it was roughly estimated to experience a 100 A magnitude GIC . When examining time series over a 27-hr period during the 29-31 October 2003 storm, Rodger et al. (2017) noted that the B measurements sampled at 0.2 Hz, "have a more similar time variation to that seen in the GIC measurement…suggesting that the higher time resolution magnetic field measurements are better at capturing the driving of the GIC." As we have shown in Figure 7, this is due at least in part to the removal of ULF wave activity in 60 s measurements. Although 60 s B (H′) measurements didn't match the time variation of the ISL M6 GIC measurements during many periods of peak GIC for the 29-31 October 2003 storm, there were other periods of the storm with elevated GIC where they matched well. In general, the suitability of H′ (which tends to overweight higher frequencies, Heyns et al., 2021) as a proxy for E and the correlation of H' with GIC will vary by location and event depending on local ground conductivity, the frequency content of H in a given event, and the power system of interest.

Discussion
In this section, we return to the question from Trichtchenko (2021), "What should be the sampling rate (or Nyquist frequency) in order to provide an adequate representation of…pulsations…for use in the calculations of the extreme geoelectric field values and in the GIC modeling?" Recall from Section 1 that results from previous studies suggested that the largest amplitude ULF-wave driven B and GIC would occur in the 120-600 s (2-10 min) period range, thus a 60 s sampling interval would be appropriate for extreme value analysis and GIC modeling related to ULF waves. These results were based in large part on B and GIC measurements from magnetic latitudes >∼60°.
As this study shows, the situation is different at magnetic latitudes <60° for two reasons: (a) past theory, modeling, and observational work show that most ULF wave modes have f > ∼f Nyquist60 (Section 1, Figure 1) and (b) the observations shown in Section 2 indicate waves with f > ∼f Nyquist60 are associated with H′, E, and GIC amplitudes comparable to the largest values recorded during extended monitoring intervals in Japan and New Zealand. Thus, we argue that the appropriate range for sampling ULF waves at magnetic latitudes below 60° is 0.1-1.0 Hz (sampling interval 1-10 s), with the particular sampling rate depending on the magnetic latitude range of interest due to the magnetic latitude dependence of standing Alfvén wave frequency (Section 1). This frequency range applies to all storm phases since these waves can occur at any time and are not limited to, for example, SSC's or storm main phase. Considering the example of the 29-31 October 2003 storm, the results shown in Section 2 indicate a sampling rate of 0.1 Hz would have worked well at most locations, even the low latitude KAK station, during periods when the largest E and GIC were observed; however, there were some intervals at EYR with significant power at 0.1 Hz, suggesting that a sampling rate of at least 0.2 Hz was needed (Figure 7 and Figures S3 and S4 in Supporting Information S1). During the 24 March 1991 geomagnetic storm, SSC high frequency waves that drove intense H′ and E likely would have required a 1.0 Hz sampling rate, as noted by Araki et al. (1997) and .
The H variations shown in Figures 2-4 are associated with a variety of magnetosphere-ionosphere current systems and waves. The transient, step-like increase in H seen at multiple locations at 06:11 UT on 29 October 2003 was caused by the intensification of magnetopause currents and a magnetosonic wave, both triggered in response to a solar wind dynamic pressure increase; together they generate the step-like ground magnetic disturbance seen at low latitudes (Araki, 1994). There are also several types of ULF waves present in these figures. As seen in Figure 1, the frequencies of the ULF wave activity seen in the blue curves in Figures 2-4 are consistent with standing Alfvén waves and plasmaspheric cavity modes/virtual resonances, whereas the wave activity seen in the red curves may also include global cavity/waveguide modes reported in past studies (e.g., Marin et al., 2014). These current systems and waves can occur simultaneously. For example, the magnetosonic wave triggered by the solar wind pressure increase can itself trigger (i.e., form the first wave cycle of) a standing magnetosonic wave if sufficient wave energy is reflected in the magnetospheric cavity (e.g., Takahashi et al., 2018;Yu & Ridley, 2011), and standing magnetosonic waves and standing Alfvén waves can occur simultaneously and with multiple harmonics as discussed in Section 1. ULF waves are commonly observed in conjunction with SSC, similar to the 29 October 2003 SSC (Araki et al., 1997;Saito & Matsushita, 1967;Wedeken et al., 1986).
When selecting an appropriate sampling rate for ULF waves for GIC studies, there are other important factors to consider besides the magnetic latitude: 1. The conducting structure of the Earth: As shown in Bedrosian and Love (2015), even a spatially uniform, sinusoidal B results in a non-uniform E with regional variations in polarization and intensity that depend on the frequency of B. Different regions may respond more (or less) to higher frequency geomagnetic variations (Grawe et al., 2018). 2. The properties of the power system of interest: Recent modeling and observational work suggests that sinusoidal E oscillations may couple to GIC in different ways when compared to more impulsive or irregular variations with longer periods. The coupling is frequency dependent, and the coupling at different frequencies depends on the properties of the power system/network such as the network orientation and reactive power response of transformers in the network (Heyns et al., , 2021Oyedokun et al., 2020). When considering the impact of GIC on a power system, it is also clearly important to allow for factors that control GIC-produced half-cycle saturation, leading to the production of even order harmonic distortion (e.g., Clilverd et al., 2018Clilverd et al., , 2020Rodger et al., 2020). 3. The choice of filter: To reduce aliasing of higher frequency signals, many 60 s measurements provided via INTER-MAGNET use an anti-aliasing filter that attenuates signals close to f Nyquist60 , leading to differences in maximum values of H′ when a digital filter with a relatively flat frequency response near f Nyquist60 is used ( Table 2, Table  S1 in Supporting Information S1, compare INTERMAGNET values to resampled low-pass filtered values). As shown in Figures 6 and 7 (compare red and gray curves in Figures 6c, 6d, 7c, and 7d) and Table S2 in Supporting Information S1, even when a filter with a relatively flat response is used, resampling the filtered data to 60 s can lead to differences in maximum/minimum values depending on the sampling convention (centered at second 00 of minute, second 30 of minute, etc.) due to the short duration that the maximum/minimum values occur. In general, the filter properties should be carefully considered in conjunction with the sampling rate when estimating extreme values of B, H′, E, and GIC related to ULF waves with f ∼ f Nyquist (Figure 1).
As discussed in Section 1 and shown in Sections 2.2-2.4, ULF waves with f > f Nyquist60s are not captured in measurements with sampling intervals of 60 s, whereas lower frequency variations and impulses/step-like variations with a broadband frequency spectrum are fully or at least partially captured in these data. When these higher frequency waves are present, event-specific conclusions concerning the magnetosphere-ionosphere current systems that drive B, E, and GIC may differ significantly when examining measurements collected with different sampling intervals. Figure 8 illustrates these points further. Three different situations are shown where the magnetopause current intensification related to a solar wind pressure pulse dominates H (Figures 8a and 8d) (Figures 2 and 3). The time series shown in Figures 8d-8f indicate that the contributions of the waves are removed when using 60 s data; in contrast, variations related to the step-like increase in magnetic field caused by the intensification of magnetopause currents are mostly preserved as they have a broadband frequency spectrum. If 60 s B measurements are subsequently used for calculations of E and predictions for GIC, the contributions of the waves with f > f Nyquist60 to E and GIC would not be included.
These results have implications for choosing appropriate models of the underlying magnetosphere-ionosphere phenomena that drive the E and hence GIC. For example, while several studies have used self-consistent global MHD simulation codes to study B related to magnetopause current intensifications and transient high-latitude field-aligned current systems (e.g., Fujita et al., 2003;Ozturk et al., 2017;Yu & Ridley, 2011), these same simulation codes are not always well-suited to examining ULF waves. There are several reasons for this, including over-damping due to the grid resolution being insufficient (e.g., Claudepierre et al., 2009;M. Hartinger et al., 2014M. Hartinger et al., , 2015, lack of appropriate boundary conditions such as plasmaspheric density (e.g., Claudepierre et al., 2016), and time steps that are too long to capture wave activity and transients with f > f Nyquist60 (Shi, Lin, et al., 2022), all of which significantly impact the ability to model the frequency and amplitude of ULF waves that occur at magnetic latitudes below 60°. When examining 60 s time series similar to the results in Figures 2-8, one may conclude that global MHD simulations with coarse spatial grid resolution and 60 s time steps would be sufficient to model the B during, for example, the SSC, whereas examination of 0.1-1.0 Hz measurements would indicate that higher grid resolutions, shorter time steps, and constraints on plasmaspheric mass density are needed to capture the wave activity.
As discussed in Section 1, estimates for extreme B, E, and GIC increase when higher sampling rates are used, though the amount of the increase depends on location. We cannot say how large a contribution ULF waves with f > ∼f Nyquist60 make to previously reported differences, though the results presented in this study suggest this is an important question that should be tested with further measurements sampled at 0.1-1.0 Hz collected during major geomagnetic storms, supplemented by numerical simulations of extreme events. For example, as noted in Section 2.2, large amplitude waves with frequencies >∼f Nyquist60 were observed following SSC in southern Sweden in the Halloween storm; other SSC events with very similar wave activity have occurred and should also be investigated to determine when and where these waves represent a hazard (Wedeken et al., 1986). Rosenqvist et al. (2022) recently simulated an extreme SSC event in this region, finding that waves with similar or higher frequencies led to modeled GIC of roughly ±50-100 A (Figure 8 in that study). Similar simulations are needed for other extreme driving conditions and other regions to better estimate wave amplitudes during extreme events, given the limited availability of 0.1-1.0 Hz data during past time intervals with major geomagnetic storms.

Summary
In this study, we examined measurements during the 29-31 October 2003 storm and past ULF wave theory, modeling, and observational work to show (a) the appropriate sampling rate for capturing B, E, and GIC variations related to ULF waves is 0.1-1.0 Hz, with the specific value in that range dependent on magnetic latitude and (b) waves with frequencies >∼f Nyquist60 can drive significant and potentially hazardous GIC at magnetic latitudes below 60°: 1. Theory, modeling, and observation work all indicate that many ULF wave modes (plasmaspheric fast mode resonances, standing Alfvén waves) occurring at magnetic latitudes below 60° have f > f Nyquist60 (Section 1, Figure 1). 2. In this study, ULF wave variations in B, E, and GIC measurements with f > ∼f Nyquist60 are present during several portions of the 29-31 October 2003 geomagnetic storm. These waves are undersampled in 60 s time series, leading to underestimates in maximum/minimum values of H′, E, and GIC variations by roughly a factor of two or more when using 60 s data (Table 2). For example, during a period of ULF wave activity on 31 October 2003 02:00-03:00 UT, the range of EYR H′ was 993 and 173 nT/min in 5 and 60 s (INTERMA-GNET) measurements, respectively. 3. ULF waves drive or contribute significantly to ISL M6 GIC variations in New Zealand with amplitudes comparable to the maximum value of 34.1 A occurring during a 14-year monitoring interval and to values related to a transformer failure. These waves also likely contribute to the peak value of 34.1 A itself, which occurred during the 29 October 2003 SSC. Thus, ULF waves with frequencies >∼f Nyquist60 are associated with significant and potentially hazardous GIC.
These results differ from past reports of B and GIC from magnetic latitudes >60°, where 2-10 min wave periods (f Nyquist60 > f > 0.0017 Hz) were associated with the largest amplitude ULF wave B and GIC. At magnetic latitudes <60°, waves with periods <∼2 min can drive or contribute significantly to some of the largest E and GIC, for example, the largest GIC reported during a 14-year monitoring interval in New Zealand. A significant part of the world's population and industry fall in the geographic region where these shorter period waves occur.
The use of 60 s data misses contributions from ULF waves that occur at magnetic latitudes <60° (Figure 1) since these waves have their frequency content concentrated above f Nyquist60 . In contrast, B and E variations at lower frequencies or with a broadband frequency spectrum are fully or partially retained in 60 s data. The absence of these waves in 60 s data can prevent the identification of these waves as a cause of extreme GIC and thus affect the choice of appropriate methods for modeling the GIC. It is likely these waves affect extreme value estimates in at least some locations (e.g., discussion in Section 2.3).
Our results indicate that measurements sampled at 0.1-1.0 Hz and numerical simulations are needed to determine extreme values of B and E that may be associated with these waves. In the future, B, E, and GIC measurements consistently recorded at 0.1-1.0 Hz are needed at more locations to determine when, where, and how often these waves may represent a hazard to power systems at different magnetic latitudes, with different power system configurations, and with different local ground conductivities. Such data can be used, for example, to tailor global MHD simulation configurations used for space weather forecasts (T. Pulkkinen et al., 2021) to include contributions from ULF waves.

Data Availability Statement
The geomagnetic (quasi-dipole) latitudes provided were obtained using the IGRF model and calculator available from the British Geological Survey (http://www.geomag.bgs.ac.uk/data_service/models_compass/coord_calc. html). The filtering software were obtained from the publicly available Python SciPy signal module (https:// docs.scipy.org/doc/scipy/reference/signal.html#module-scipy.signal). The 60 s ground-based magnetic field data from FRD, EYR, KAK, UPS, and NGK are available from the INTERMAGNET repository (https://www. intermagnet.org In addition, we are unable to directly provide the New Zealand LEM DC data or the derived GIC observations. Requests for access to the measurements need to be made to Transpower New Zealand. At this time the contact point is Michael Dalzell (michael.dalzell@transpower.co.nz). We are very grateful for the substantial data access they have provided, noting this can be a challenge in the space weather field (Hapgood et al., 2016).