Climatological Statistics of Extreme Geomagnetic Fluctuations with Periods from 1 s to 60 min

Using a global database of 125 magnetometers covering several decades we present occurrence statistics for fluctuations of the horizontal geomagnetic field ( 𝑑𝐵 ℎ /𝑑𝑡 ) exceeding the 99.97 th percentile ( P 99.97 ) for both ramp changes ( R n ) and the root-mean-square ( S n ) of fluctuations over periods, 𝜏 , from 1 to 60 min and describe their variation with geomagnetic latitude and magnetic local time (MLT). Rates of exceedance are explained by reference to the magneto-ionospheric processes dominant in different latitude and MLT sectors, including ULF waves, interplanetary shocks, auroral substorm currents, and travelling convection vortices. By fitting Generalised Pareto tail distributions above P 99.97 we predict return levels (RLs) for R n and S n over return periods of between 5 and 500 years. P 99.97 and RLs increase monotonically with frequency ( 1/𝜏 ) (with a few exceptions at auroral latitudes) and this is well modelled by quadratic functions whose coefficients vary smoothly with latitude. For UK magnetometers providing 1-s cadence measurements, the analysis is extended to cover periods from 1 to 60 seconds and empirical Magnetotelluric Transfer functions are used to predict percentiles and return levels of the geoelectric field over a wide frequency range ( 2 × 10 −4 to 4 × 10 −2 Hz) assuming a sinusoidal field fluctuation. These results help identify the principal causes of field fluctuations leading to extreme geomagnetically induced currents (GIC) in ground infrastructure over a range of timescales and they inform the choice of frequency dependence to use with 𝑑𝐵 ℎ /𝑑𝑡 as a GIC proxy.


Introduction
Large electrical currents are occasionally induced in ground-based infrastructure as a result of rare and intense currents in the ionosphere or magnetosphere. These Geomagnetically Induced Currents (GIC) have been identified as a substantial hazard to national infrastructure (Cannon et al., 2013;Hapgood et al., 2021) since they may cause catastrophic failure in high-voltage electricity supply networks (Gaunt, 2016;Oyedokun & Cilliers, 2018;Thomson et al., 2010), damage long-cable communication systems  and cause railway signalling errors (Boteler, 2021;Eroshenko et al., 2010;Wik et al., 2009). The cumulative effect of GICs above a certain threshold may also cause corrosion in oil and gas pipelines (Boteler, 2000;Pulkkinen et al., 2001). The science of GICs and their effects is reviewed in Knipp, 2011, Chapter 13, andBuzulukova, 2017, Chapter 8. Modelling the risk of extreme GICs requires a statistical characterisation of the geoelectric field, , induced by electrical currents in the ionosphere and magnetosphere. This information may, for example, be combined with a model of electrical impedances in a high-voltage (HV) electricity network ) to determine the 'return level' (RL) of GIC expected in a 'return period' of 100 years or more. Direct measurements of are often subject to contamination from anthropogenic electromagnetic interference and require an experienced expert to remove noise and biases (Kelbert et al., 2017). They are also not global in extent, and do not cover the decades required for accurate prediction over long return periods. For climatological studies it is therefore expedient to instead use an archive of measurements of the rate of change of the horizontal component of the geomagnetic field, ℎ ⁄ , measured at ground level. Using Faraday's law of induction (Faraday, 1832) and magneto-telluric (MT) theory (Cagniard, 1953;Chave & Jones, 2012) these may be combined with a model of the local ground conductivity to determine climatological statistics for . Alternatively, may be derived using collocated measurements of ground impedance at a magnetometer site.
The calculation of E requires knowledge of both the temporal spectrum of geomagnetic oscillations and the frequency dependence of the surface impedance. Databases of impedance tensors are increasingly available for public use (e.g. Kelbert et al. 2011;Kelbert et al. 2018) and can cover a wide frequency range corresponding to periods from milliseconds to hours. The most effective source of geoelectric fields producing damaging GIC in power transmission lines lie in 1-1000 s period oscillations (Kappenman, 2004;Barnes et al., 1991) and electricity companies have identified that fluctuations on timescales from tens of seconds to over an hour have led to vulnerability of high-voltage (HV) electricity networks to GIC (e.g., NERC, 2017; Girgis & Vedante, 2012). A well-reported example is the geomagnetic storm of 13 March 1989 in which the 21 GW Hydro-Québec power supply failed for nine hours following horizontal geomagnetic field fluctuations | ℎ / | of approximately 500 nT/minute (p.640, Knipp, 2011).
The frequency of the induced E field fluctuations and consequent GICs is much less than the frequency of highvoltage electricity networks (50 or 60 Hz) and so is often modelled as a quasi-direct current. Currents of more than a few amperes sustained over periods similar to the thermal time constants of the components of a highvoltage transformertypically 30-45 minutesmay cause irreversible damage resulting in power failures (p.8, IEEE, 2015;Girgis & Vedante, 2012;Erinmez et al., 2002;Molinski, 2002;NERC, 2017). GICs generated by field fluctuations with periods longer than 1 hour have amplitudes too small to be of concern, whilst sub-1-s

Accepted Article
This article is protected by copyright. All rights reserved.
fluctuations are heavily damped by inductances in electric power systems .
Understanding the climatology of extreme | ℎ / | over periods from 1 s to 1 hour should, therefore, help to quantify the GIC risk to electrical power systems.
Large-scale statistical surveys often exploit measurements at 1-min resolution, in large part enabled by the successful SuperMAG project (Gjerloev, 2011), thus many have examined only the 1-minute changes in ℎ , (denoted R1), with this metric being adopted as a proxy for GICs (e.g. Viljanen et al., 2001Viljanen et al., , 2015Thomson et al., 2011). However, probability distributions of | ℎ / | are observed to depend strongly on the time resolution (or sample averaging period) of the B field measurements, with lower amplitudes at longer sampling intervals due to the effect of smoothing. In recent years, an increasing number of magnetometer operators have offered users measurements at 1-s cadence and so the question arises as to which temporal resolution to apply when using | ℎ / | as a proxy indicator for GIC. Modelling by  showed that smoothing the B-field components from their native resolution of 1 s up to 60 s reduced the amplitude of | ℎ / | by 80% whilst the computed peak E-field amplitudes were reduced by only 20%, the inference being that a 60-s (but no more) sample interval is acceptable as a proxy to use for E-field (and hence GIC) calculations. Other studies have noted that rather than taking R1 as a proxy for GIC, a better performing indicator was obtained by taking an average ℎ / over 20-minutes (Tõth et al., 2014) or 30 minutes , whilst others have used the hourly range or standard deviation (Beamish et al. 2002;Nikitina et al., 2016;Danskin & Lotz, 2015) or 3-hourly range indices as a proxy (Trichtchenko & Boteler, 2004).
In several cases, the magnitude of ℎ relative to its quiet-day value (often denoted ) has provided a better proxy for GIC than ℎ / (Pulkkinen et al. , 2010;Tõth et al. , 2014;Watari et al. , 2009). Pirjola (2010) showed how this is more likely to arise in regions for which there is an upper, highly conductive layer overlying a deeper layer of low conductivity. Heyns et al. (2020) presented examples of GIC amplitudes and phases matching closely to the 20-min period fluctuations of the field ( ) which were poorly represented by high-cadence ℎ / indicators, whilst ℎ / was a better indicator of the rapid field variation that occurred during Sudden Commencements, which often initiate geomagnetic storms. Heyns et al. (2020) explained that this is because the B field (or ) has low-frequency components that are deweighted when taking the time derivativefor example, if ℎ = 0 exp( ) then frequency components of | ℎ / | are weighted by the factor 1/ . Consequently, 1-s resolution ℎ / measurements (R1/60) would be even less effective as a proxy for GIC (compared to R1) for GIC caused by field fluctuations of a much longer period. Power networks can respond strongly to B-field fluctuations over tens of minutes, indicative of finite reactive impedances in the network components, and assumptions that the geomagnetic driving is d.c. in nature may be insufficient to replicate the observed GIC (Heyns et al., 2020;Jankee et al., 2020).
The study of extreme geomagnetic fluctuations over a range of periods yields much information about the causes and impacts of GIC as well as the drivers of these fluctuations. The ionospheric and magnetospheric processes contributing to ℎ / over a 1-min period will differ greatly from those at 60 min and will depend on the latitude, magnetic local time (MLT), season, and other factors. The principal drivers of short transients (timescales of minutes) may be categorised into the following phenomena:

Accepted Article
This article is protected by copyright. All rights reserved.

Sudden Commencements (SC):
Interplanetary shocks arriving in the solar wind, which generate a sudden eastward (dusk-to-dawn) Chapman-Ferraro current at the dayside magnetopause, are observed as Sudden Commencements (SC) in magnetograms (Fiori et al., 2014;Kappenman, 2003;Smith et al., 2019). The characteristic rapid magnetic field variation may be short-lived, lasting several minutes or up to an hour (Knipp, 2011, p.496), and are associated with ℎ / of up to 30 nT/min at low geomagnetic latitudes (< 40°) or up to a maximum of 270 nT/min in the auroral zone (approximately 65° geomagnetic latitude) (Fiori et al., 2014).

Auroral substorm onsets:
A substorm is the sudden brightening and expansion of auroral arcs resulting from bursts of energetic electron precipitation from the magnetotail (Akasofu, 2017;Ieda et al., 2018).
This enhances the ionisation and electrical conductivity of the ionospheric E region allowing strong Hall currents to flow, most often in a westward direction which manifest in magnetograms as a rapid decline in the north component of the geomagnetic field, . Substorm onsets have been categorised by Newell & Gjerloev (2011) from the SML geomagnetic index (which measures the lower envelope of ) as a reduction of at least 45 nT over 3 minutes followed by a mean level at least 100 nT below the initial value during the half-hour following onset.

Day-time Magnetic Impulse Events (MIE):
Pairs of up-and down-field aligned currents generated by a pulse in dynamic pressure at the dayside magnetopause couple into the ionosphere as Travelling Convection Vortices (TCV) at latitudes in the vicinity of the dayside cusp/cleft (approximately 77-78° magnetic) (Zesta et al., 2002;Kataoka, 2003;Engebretson et al., 2013;Friis-Christensen et al., 1988;Lanzerotti et al., 1991). Magnetometers in this region observe the ionospheric Hall current loops (a pair of vortices) as isolated magnetic impulse events (MIE) in magnetograms, lasting typically 5-15 minutes with amplitudes of typically 50-200 nT or up to a maximum of 400 nT (Kataoka et al., 2003;Lanzerotti et al., 1991). Several mechanisms have been postulated to explain the generation of TCVs near the dayside magnetopause, including bursts of magnetic field line reconnection (flux transfer events), solar wind pressure pulses, plasma injections into the low-latitude boundary layer, Kelvin-Helmholtz instabilities, and perturbations of the ion foreshock upstream of the Earth's bow shock (see references in Kataoka et al., 2003 andEngebretson et al., 2013). In general, TCVs are defined so as to exclude sudden commencement perturbations associated with a large interplanetary shock (e.g. Pilipenko et al., 2019).
A significant number of GIC events occur under geomagnetic storm conditions at auroral and mid-latitudes due to sustained ULF pulsations in the Pc5 band (2.5-10 min period field oscillations (Baker et al., 2003;McPherron,
B-field fluctuations over tens of minutes may also arise from the expansion and recovery phases of substorms in the auroral zone Pothier et al., 2015): The substorm expansion phase typically lasts 25-40 minutes (Pothier et al., 2015) followed by a more gradual recovery phase. Changes over an hour or more may arise from slow changes and movements of an electrojet over a magnetometer station or from gradual changes of the magnetospheric inner ring current intensity during the main and recovery phases of a geomagnetic storm.
At very high latitudes (poleward of the dayside cusp) and under conditions of northward interplanetary magnetic field (IMF) and large dipole tilt (e.g. at summer noon), magnetic fluctuations may be associated with the merging of 'overdraped' tail-lobe field lines with the IMF (Crooker, 1992;Watanabe et al. 2005). Rogers et al. (2020) postulated that field-line reconnections may drive impulsive 'Region-0' field-aligned currents (Wang et al. 2008;Milan et al. 2017) into this region that could manifest as large | ℎ / | fluctuations at the surface.
In this paper we have extended a global climatological statistical model of extreme 1-minute fluctuations, R1, (Rogers et al., 2020) to include the magnitude and frequency of occurrence of extreme | ℎ / | over sampling periods between 1 and 60 minutes, both as ramp changes (applying a moving average of the geomagnetic field measurements) and as a root-mean-square (RMS) of the R1 values over n-minute periods that we denote for n = 1-60 (defined explicitly in Section 2). The latter is a measure of the sustained power in extreme geomagnetic field fluctuations, which is important in modelling the risk to transformer components due to heating, for example.
Our study complements that of Love et al. (2016a) who provided an analysis of extreme | ℎ / | over 1-and 10-minute periods, (R1 and R10), and the RMS of R1 over 10 minutes (S10). Wintoft (2005) and Wintoft et al. (2005) also chose to study S10 as a predictor of the RMS GIC amplitude. Part of our study will focus on three UK magnetometer sites, and as such complements the work of Beamish et al. (2002) who examined the hourly standard deviation of 1-min B-field north and east components (independently), a measure similar to the S60 calculated in this paperand the works of Beggan et al. (2013) and Beggan (2015), who estimated extreme Efield and GICs for the UK national grid at 100-and 200-year return periods using UK ground conductivity models for 2 and 10-min period fluctuations of the inducing B-field, with amplitudes inferred from predicted extremes of R1 presented by Thomson et al., (2011).
In section 2 we describe the processing of magnetometer measurements data set and the determination of extreme values for | ℎ / | as both ramp changes and RMS fluctuations. Section 3 presents the latitude and MLT distributions of large percentiles and projected extreme values for a range of sampling frequencies and develops

Accepted Article
This article is protected by copyright. All rights reserved.
a global model to characterise the dependences on sampling frequency. The frequency range is extended up to 1 Hz sampling for three UK sites, and for these locations empirical MT transfer functions (or surface impedance matrices) are used to predict high percentiles and extreme values of the geoelectric field.

Measurements
Magnetic field measurements (magnetograms) were obtained from 125 magnetometers in the global SuperMAG collaboration (Gjerloev, 2011) at sites for which at least 20 years of data was available, with an average of 28 years' data per site. Table 1 provides the locations of these magnetometer sites in geodetic and corrected geomagnetic (CG) coordinates (Laundal & Richmond, 2017;Shepherd, 2014). Due to the secular variation of the Earth's main field, CG coordinates are given as averages over all years in which magnetometer data was available at each site. In this paper we consider only the north and east components of the magnetic field ( and , respectively) in local magnetic coordinates (Gjerloev, 2012) neglecting the downward vertical field component, , which contributes little to GICs in surface-based infrastructure. The magnetograms provided by SuperMAG had already been cleaned and manually inspected to remove most artificial sudden changes in the baseline (offsets), spikes, and gradual slopes (Gjerloev, 2012). Nonetheless, as a further check, all data in weeks containing R1 peaks above the 99.97 th percentile ( 99.97 ) were visually inspected and obvious artefacts (such as large spikes, step changes, and instrument saturation effects) were replaced by data gaps, as described in (Rogers et al., 2020).
At each magnetometer, the 'ramp' change in the horizontal component of over n-minute intervals was defined as = { ( ): = 1,2,3, … , } (1) where is the number of field measurements, and Δ = 1 minute was the cadence of the measurements. For computational efficiency, the n-minute backward difference values, , were calculated using n-point movingaverage filters on the 1-minute first differences of and . Intervals containing missing data were excluded from the analysis. The statistics of (1-min field fluctuations) were modelled in (Rogers et al., 2020). The definition in (2) ensures that statistics of the induced E-field magnitude, | | = √ 2 + 2 will be approximately proportional to (with exact proportionality for an idealised half-space model of surface conductivitysee Annex A). The expression for 1 is the same as that adopted by The root-mean-square of over n-minute periods was defined as with

Accepted Article
This article is protected by copyright. All rights reserved.
and this was implemented in software using a convolution filter. Since we are only interested in extreme values, a high threshold for and was set at the 99.97 th percentile level, P99.97. The application of extreme value statistics (Coles, 2001) requires an assumption that exceedances of this threshold are temporally independent rather than clustered together. Therefore, the threshold exceedances were declustered to ensure a minimum 12 hours between clusters and only the peak value in each cluster was recorded. The magnetic local times (MLT) (Laundal & Richmond, 2017) associated with each peak were also calculated as described by Rogers et al. (2020).
Declustered exceedances ( > 99.97 ) were then fitted to a Generalised Pareto (GP) 'tail' distribution and the fitted GP profile was used to predict return levels (RL) expected over return periods (RP) of up to 500 years (see (Rogers et al., 2020) and (Coles, 2001) for mathematical details). The analysis of extreme field fluctuations at 28 European magnetometer sites by Thomson et al. (2011) showed that the choice of a P99.97 threshold and 12-hour declustering provides relatively stable GP coefficients whilst ensuring temporal independence of the extreme events. For consistency of approach we have therefore adopted these thresholds for our analysis of magnetometer data worldwide.
A further set of magnetometer measurements at 1-s cadence were obtained for three sites in the UK operated by the British Geological Survey, namely, HAD (Hartland, southern England, CG latitude = 47.55°), ESK (Eskdalemuir, southern Scotland, = 52.65°), and LER (Lerwick, Shetland Is, northern Scotland, =57.97°) (see Table 1). Data at 1-s resolution were available from 1 Jan 2001 to 14 Sep 2016 for all three sites, whilst the 1-min SuperMAG data set extended from 1 January 1983 to 31 December 2016 for all three sites. The data were visually inspected for weeks containing 1-s | ℎ / | ( 1/60 ) exceeding the 99.97 th percentile, and obvious artefacts removed in the same manner as for the 1-minute SuperMAG data set described above. When fitting GP distributions to predict return levels for the 1-s datasets, and for all averaging periods, ≡ Δ , we have used the same consistent percentile threshold (99.97 th percentile) and declustering run-length (12-h) as selected in the study of 1-min cadence measurements by Thomson et al. 2011 using a selection of visual diagnostics. Thomson et al. (2011) noted that, for most geomagnetic observatories in their study, the return level was "only weakly dependent on the decluster length". Historically, in extreme value statistical analyses, justification of the threshold selected has been through visual diagnostics combined with any available scientific insight or expert knowledge on the process of interest. The use of visual diagnostics becomes infeasible as the number of data sets (in this case sitefrequency combinations) grows. Consequently, we took a pragmatic approach and defined the same proportion of observations to be the tail sample for each site and for each value of τ. For a given τ, this permits comparison of return levels across sites, and for a given site it permits comparison across all values of τ. As the duration of the 1-s datasets obtained from the UK observatory sites differed from those in the 1-min data set in the SuperMAG archive, we ran additional visual diagnostic checks for thresholds at 99.95, 99.97 and 99.99 percentiles (with and without 12-h declustering applied) as were performed in the analysis by Thomson et al. (2011). These checks confirmed that for the Δ = 1 s dataset, 12-h declustering and the 99.97th percentile threshold remained the most appropriate for return level estimation at all three UK observatories.

Accepted Article
This article is protected by copyright. All rights reserved.
The measurement of the ground magnetic field has a long established tradition in many countries and data quality and standards are set to a high level, e.g. through INTERMAGNET (Thomson & Flower, 2021;Love & Chulliat, 2013). In contrast, long-term observations of the ground electric field are relatively rare (Beggan et al., 2021 and references therein) and more influenced by man-made electromagnetic noise due to a low signal-noise ratio.
Available data sets are scarce and often discontinuous. In the UK, the ground electric field has been monitored at the three geomagnetic observatories (HAD, ESK, and LER) since 2015 with non-polarizable electrodes along north-south (N-S) and east-west (E-W) oriented baselines. (Some recent examples of these measurements are available online: http://www.geomag.bgs.ac.uk/data_service/space_weather/geoelectric.html.) To obtain estimates of the geoelectric field for times when no data was recorded, the horizontal geoelectric field spectrum, ( ) = ( ) may be estimated from the horizontal magnetic field spectrum ( ) = ( ) via where is the permeability, and ( ) = ( ) is the impedance (with units of Ω), where x and y refer to north and east components, respectively (e.g., Chave & Jones, 2012, Simpson & Bahr, 2005. Fourier transforms may be used to convert between the frequency (f) and time domains. The frequency-dependent term / is called the Magnetotelluric (MT) transfer function (with units of V/km/nT) and is informative of the electrical conductivity structure of the subsurface that is useful in deep geophysical exploration.
/ was estimated from simultaneous measurements of the horizontal components of the ground electric and magnetic field using robust statistical approaches to minimize the influence of noise. For the estimation of / at HAD, ESK and LER, we used six months of electric and magnetic field measurements from 2015 and applied the impedance estimation algorithm of Smirnov (2003). Further details of the procedure are given in (Beggan et al., 2021). Due to the sampling cadence of 1 s and the frequency response of the fluxgate magnetometers at the observatory sites, the impedance estimates cover a period range of 20 to 20,000 s (or 5×10 -2 -5×10 -5 Hz).
3 Latitude, MLT, and Seasonal distribution of large R n and S n on timescales from 1 to 60 min

Accepted Article
This article is protected by copyright. All rights reserved. To gain a better understanding of the physical drivers of these large fluctuations, we first examine their magnetic local time (MLT) dependence. Figure 2 presents the probability of (declustered) peaks of | ℎ / | exceeding 99.97 as a function of | | and MLT. This was calculated by counting the number of peaks in 1-hour bins of MLT and 3.3° bins of | |, where data from multiple magnetometers were aggregated where they lay within the same latitude bin. (Bin sizes were chosen as a compromise between resolution and quantisation noise.) The bin counts were then normalised by the total number of field measurements in each bin. Panels (a-d) present the distributions for ramp changes over 1, 10, 30 and 60 min, respectively, whilst panels (e-h) present the distributions for the RMS magnitudes over 1, 10, 30 and 60 min, respectively. We have used absolute latitude on the vertical axes since the distributions of occurrence probability against (signed , MLT) were, to a close approximation, symmetric about = 0. Note that panels a) and e) are identical, which may be noted from Equation (4) with =1.
When interpreting the distributions in Figure 2 it is important to remember that the threshold 99.97 itself varies with | | (see Figure 1) and as such it is simplest to focus on the MLT distribution in each individual latitude band.
It is also important to note that, due to the method of declustering, peaks occurring within 12-hours of a larger peak are not represented. However, it was observed that if the peaks over threshold were not declustered, then the general shape and form of the probability distribution in Figure 2 remained largely unaltered; for | | > 40°, with no declustering, the occurrence probabilities were slightly reduced in the hours 12-24 MLT and slightly raised in the hours 0-12 MLT, indicating a greater clustering of peaks associated with events occurring pre-noon.

Accepted Article
This article is protected by copyright. All rights reserved.
At the highest latitudes (| | > 80°), poleward of the dayside cusp, there is an occurrence maximum in the few hours about noon MLT, which persists over all timescales (1-60 min). For > 1 min, the maximum is much more sharply peaked for ramp changes than for RMS fluctuations, and as increases towards 60 min the MLT of the maximum occurs slightly later (towards 14 MLT). (Note that the timestamps and MLTs associated with each cluster peak of | ℎ / | refers to the end of the -minute period in question (from Equations (2) and (4)) but this is not sufficient to account for the apparent shift of the maximum towards the post-noon.) Analysis of the R1 distribution by Rogers et al. (2020) showed that these peaks near noon occur predominantly under northward IMF conditions during the summer months (i.e. under conditions of greatest dipole tilt angle), suggesting a possible relation to impulsive field line reconnection between the IMF and an 'overdraped' tail lobe (Wang et al., 2008;Milan et al., 2017;Crooker, 1992;Watanabe et al., 2005). The MLT distribution of occurrence probability at dayside cusp latitudes does not match the distributions of MIEs observed by Lanzerotti et al. (1991) and Kataoka et al. (2003) who reported a relatively flat distribution over 06-18 MLT with a minimum around 11 MLT, although these MIE distributions were not thresholded at a very high percentile. Nonetheless, the MIE amplitude distribution presented in Fig. 5d of (Kataoka et al., 2003) indicates perturbations approaching 400 nT (over ~5-15 min) in the 07-11 MLT period, which is not observed in the MLT profile of 99.97 exceedances of Figure 2a.
Such discrepancies indicate that it is less likely that MIEs (caused by TCVs) provide a significant contribution to the extremes of | ℎ / | in this region.
At low latitudes | | < 40°, for R1 and R10, and Sn for all n, the occurrence probabilities increase on the day-side at 07-16 MLT, although for 20°< | | < 43° the distribution is double-peaked with a dip in occurrence in the few hours around noon, creating a Y-shaped pattern most clearly discernible in the 1-min data (panels (a) or (e)). The distributions for R10 and R30 also have a night-time maximum in the period (19-03 MLT). Rogers et al. (2020) showed (in their Fig. 8) that approximately 25-70% of the R1 peaks at these latitudes occurred at or within 30 minutes of a sudden commencement, as recorded with high confidence in IAGA bulletins (http://www.obsebre.es/en/rapid). However, the lower figure (25%) was associated with the largest occurrence probabilities near noon, suggesting that alternative or delayed driving processes may be contributing to the largest R1 at these times.

Accepted Article
This article is protected by copyright. All rights reserved.

Ramp
RMS a) e)

Accepted Article
This article is protected by copyright. All rights reserved. of R1 peaks exceeding P99.97 were associated with the expansion or recovery phase of a substorm.
A secondary peak of occurrence is observed in the dawn-noon sector. Some of these peaks below 70°N may be associated with MPEs since they are consistent with the 02-06 MLT distribution observed by Engebretson et al. (2021) for the station at 65°N geomagnetic, as noted above. However, this is also a region in which Pc5 pulsations are the dominant wave activity (e.g. Engebretson et al., 1998;Pulkkinen & Kataoka, 2006). The R1 occurrence probabilities maximise at around 03 MLT at | |=60°, increasing to 12 MLT at | | =80°, and similar patterns have been reported in the distribution of Pc5 wave power (compare, for example, Fig. 5 of Vennerstrøm (1999), Fig. 2 and 4b of Baker et al., 2003, or Fig. 1 of Weigel et al., 2002). The rate of occurrence for longer-period ramp changes, R10, R30 and R60, is suppressed in the latitude band | | =70-77°, although this may be an effect of declustering where the peaks occur within 12 hours of larger amplitude fluctuations in the pre-midnight sector.
In contrast to the distribution of ramp changes, the occurrence patterns of large RMS fluctuations (Figure 2e-h) show that as the period, increases, the probability of occurrence ( > P99.97) in the auroral zone increases strongly in the dawn sector (03-07 MLT). A cursory inspection of magnetograms for the largest peaks of indicated that many are indeed associated with ULF wave activity lasting tens of minutes (see, for example, Fig.   1c of (Rogers et al., 2020). To examine this further, an analysis of the probability of occurrence vs (month, MLT) is presented in Figure 3 for the 26 sites at latitudes = 60°-70°N. This figure shows that in the pre-midnight hours the frequency of occurrence is greatest near the equinoxes, when the geomagnetic field is more favourably oriented for reconnection with the IMF (Russell & McPherron, 1973;Zhao & Zong, 2012). However, for RMS fluctuations ( Figure 3e-h), as increases from 1 min to 60 min, the greatest frequency of occurrence occurs on the dawn side (03-09 MLT). We also note, for both 1 and 1 distributions, a change in the locus of peak occurrence from 04-05 MLT near the summer solstice to 07-08 MLT near the winter solstice, which may be associated with changes in the position of the dawn terminator at these latitudes and the seasonal changes in the geometry of the geomagnetic field relative to the IMF. For >= 10 min, however, the frequency of occurrence in the winter months (December and January) is reduced relative to that for =1 min, in both and , and this also limits the time zones of occurrence in the late morning.

Accepted Article
This article is protected by copyright. All rights reserved.

Accepted Article
This article is protected by copyright. All rights reserved. The best-fit quadratic coefficient, 1 linear coefficient, 2 , and constant term, 3 , are presented in Figure 5

Accepted Article
This article is protected by copyright. All rights reserved.  The goodness of the quadratic fits at each magnetometer site are presented in Figure 6 for 99.97 th percentiles of (left panels, a and b) and (right panels, c and d). Panels (a) and (c) present the coefficients of determination, 2 . To better illustrate values of 2 close to 1, the vertical axis scaling in panels a and c is "inverse logarithmic" such that a set of values, = 1 − 10 − would be uniformly spaced for uniformly spaced . Panels (b) and (

Accepted Article
This article is protected by copyright. All rights reserved.

Modelling return levels
Generalised Pareto (GP) distribution functions were fitted to exceedances of above a 99.97 threshold (after 12-h run-length declustering above the same threshold) independently for each magnetometer site. 100-year return levels of were then determined from the GP distribution at a probability level equivalent to a 1-in-100 years of observations. A numerical method was used to determine a maximum likelihood estimate (MLE) for the return level with 95% confidence intervals determined from the (asymmetric) log-likelihood profile, as described in (Gilleland & Katz, 2016). This procedure was repeated for all 125 magnetometer sites and the results are plotted against | | in Figure 7. Panels a, b, c and d, present 100-year return levels of for n = 1, 10, 30, and 60 (minutes), respectively; points represent MLEs (coloured blue for southern hemisphere sites, black for northern hemisphere) with error bars indicating the 95% CI. The red curve in each panel is a smoothing-spline interpolation to the MLE values. In Figure 8 the interpolating spline curves are presented for return periods from 5 to 500 years.
The 100-year return levels for 1 (Figure 7a) (1-minute ramp changes) are distinctly elevated for sites around | | ≅ 52-54° and reference to Figure 8 indicates that the latitude of this maximum decreases with increasing return period. This indicates that the extreme 1 events (declustered threshold exceedances) that occur less frequently (i.e. with longer return periods) have greater amplitude and occur at lower absolute latitudes. This pattern of behaviour could indicate that largest and rarest auroral current fluctuations occur during substorm expansions associated with brightening auroral arcs at the equatorward edge of a greatly expanded auroral oval (i.e. following a large substorm growth phase). Over 10-60 minute timescales (Figure 7b to d) the peak near 53° is still present but less pronounced, and Figure 8b-d shows that it has similar or lower magnitude than the broad peak around | | ≅ 67° that was observed in the P99.97 profiles (Figure 1a).
The same procedure of fitting GP distribution functions was used to determine extreme values for the RMS variation over n-minute periods, . Figure 9 presents the 100-year return levels and Figure 10 presents the smoothed-spline fits for 5-500 year return periods for the , again for periods of = 1, 10, 30, and 60 minutes, The shape of these distributions are very similar to those of the fluctuations although the reduction in level with increasing is much less pronounced.
For both and metrics (Figure 8 and Figure 10, respectively) there is an increase in RLs towards the equator, potentially associated with activity in the equatorial electrojet current systems, and for return periods greater than 100 years there is a predicted increase in RL as latitude | | increases above 74°, for = 1 and 10 min.

Accepted Article
This article is protected by copyright. All rights reserved.

Accepted Article
This article is protected by copyright. All rights reserved.

Accepted Article
This article is protected by copyright. All rights reserved.

Accepted Article
This article is protected by copyright. All rights reserved. geomagnetic (Cliver & Dietrich, 2013). The measurements at low-latitude sites contain no observations of such extreme conditions and therefore may indicate misleadingly low return levels. Similar situations (in which more severe events tend to be spatially more localised) are frequently encountered in environmental and geospatial data sets and advanced methods for analysing such 'spatial extremes' are reviewed by Hüser & Wadsworth (2020).
We now present models of the MLEs of 100-year return levels of and as functions of sampling frequency, , and absolute CG latitude, | |, following the same procedure as for the P99.97 levels developed in Section 4.1. Figure 11 presents 100-year RLs for a) and b) , in the same format as Figure 4. The coefficients of the polynomials (6) fitted to the return levels are presented in Figure 12. The 95% CI of the fitted coefficients (error bars in Figure 12) are larger than for the P99.97 model, but the profiles remain approximately symmetric about =

Accepted Article
This article is protected by copyright. All rights reserved.
0°. It is interesting to note for the ramp changes, , there is a pronounced change from positive to negative curvature as | | increases, which can be seen in the profiles of Figure 11 and the change in quadratic coefficient, p1, in Figure 12. For both and , the gradients (or the linear coefficients, 2 ) are significantly higher at lower latitude | |. Figure 13 provides goodness-of-fit metrics for the polynomials (6) fitted to MLE of RL100, presented in the same format as Figure 6. Not unexpectedly, RL100 shows greater variation from the polynomial model than 99.97 (cf. Figure 6) but in the vast majority of cases the RMS errors are still less than 15% and have a coefficient of determination greater than 0.9.

Accepted Article
This article is protected by copyright. All rights reserved.

Predictions of extreme geoelectric fields in the UK
We shall now focus on the statistics for three UK magnetometer sites, HAD, ESK and LER. Figure 14 presents, for each site, the 99.97 th percentiles of a) , and b) , and 100-year return levels for c) , and d) . For the ramp changes, (Figure 14a and c) the frequency scale is extended up to 1 Hz using the 1-s cadence dataset. Whilst the length of the datasets differ for 1-s and 1-min data, the discontinuities in the P99.97 curves (Figure 14a) at = 1 min are negligible, although a larger discontinuity arises from the RL estimates (Figure 14c). Statistics for (Figure 14b and d) could not be extended to 1 Hz since they are defined from 1-min cadence measurements (Equation (4)), but they are presented here for = 1-60 min to illustrate that whilst the 99.97 th percentile varies little with sample frequency (panel b), their 100-year RLs (panel d) have a much more significant frequency dependence, albeit with large 95% confidence intervals (illustrated by the shaded regions).

Accepted Article
This article is protected by copyright. All rights reserved. To derive estimates of the 99.97 and 100-year return levels of the geoelectric field from statistics of the geomagnetic field, we make use of the MT transfer functions, ( )/ measured at each UK site, as described in Section 2. In Figure 15 (panels a, c, e) we present, for each site, the 'apparent resistivity' associated with each of the four components of the observed / , defined as where ( = , , , ) are the components of the impedance matrix, . Figure 15 panels b, d and f show the phases of . Apparent resistivity is the resistivity of an electrically homogeneous and isotropic half-space of permeability = 0 (the permeability of free space) that would be consistent with the measured E and B fields. Cagniard (1953) and Pirjola (1982) showed that, using a simple half-space model of the surface, an electromagnetic wave polarised in the N-S plane, at the surface (z = 0), the magnetic field would induce a geoelectric field, E, at the surface ( = 0) of

Accepted Article
This article is protected by copyright. All rights reserved.
in the east direction, where is the conductivity of the ground. Equation (9) is known as the "basic equation of magnetotellurics" and is valid under the assumptions that the permittivity ≪ /2 and the conductivity of the air above the surface is negligible. As discussed by Wait (1962), the plane wave approximation (8) may be used provided there is negligible change in the incident wave field amplitude over a lateral scale equal to the 'skin depth' of the ground. Considering an additional orthogonal component of the magnetic field , we may write (9) more generally as Thus, for a uniform half-space model earth, would be invariant with frequency, , and the components of Z would have a constant 45° phase for all . The measurements in Figure 15 show that the apparent resistivity of the ground differs greatly between sites as is expected from the very different geological settings that give rise to the electrical response.

Accepted Article
This article is protected by copyright. All rights reserved. For each site, the off-diagonal components ( ) and ( ) are not of equal magnitude, which indicates that the MT transfer function introduces 'directional anisotropy' (i.e. from Equation (10), | | ≠ | | when | | = | |).
The diagonal terms ( ) and ( ) are non-zero (notably for Lerwick), suggesting some deviation from the simple half-space model (i.e. measurements imply a fully three-dimensional distribution of electrical resistivity).
Noting from (8) that and similarly for , we may estimate the geoelectric field from the rate-of-change of the magnetic field: where is determined from the empirical MT transfer function ( / ) using the approximation ≅ 0 . In the ideal case of homogenous ground conductivity, equations (13) and (11) indicate that the spectrum of | | is proportional to −0.5 times the spectrum of | / | (i.e. it is low-pass filtered).

Accepted Article
This article is protected by copyright. All rights reserved.
To estimate the amplitude of the geoelectric field expected to result from the 99.97 th percentile of Rn for a range of frequencies, we modelled the waveform as a vertically propagated sinusoid 0 sin(2 ) with amplitude 0 = 99.97 /2 and frequency = 1/(2 ). This required a linear interpolation of the 99.97 (Figure 14a) to frequencies recorded in the MT transfer function at each site ( Figure 15). The resulting estimates of the magnitude | | = √ 2 + 2 from Equation (13) are presented in Figure 16a where circles represent B-field fluctuations confined to the N-S plane (| | = 0 ; = 0) and asterisks represent B-field fluctuations in the E-W plane ( = 0; | | = 0 ). The two polarisations yield E-fields that differ in magnitude by a factor of up to 2 because the MT transfer function is not directionally isotropic and the ground impedance depends on all three coordinates (x, y, z). At each UK site, exceedances of 99.97 , after declustering, occurred on average every 0.1 to 0.35 years over the range = 1 to 60 min, and so should be considered as large, but not extreme values. Figure 16b presents the E-field magnitude for sinusoids with peak-to-peak amplitude (2 0 ) equal to the 100-year return levels, (from To put these values in context, an E field of 1-2 V/km over large distances can, depending on the grid topology, produce GIC that saturates the steel core of a high-voltage transformer, which may lead to heating and potential failure of core components and the introduction of harmonics in the power system (Barnes et al., 1991). Winter et al. (2017) estimated the E field at UK latitudes associated with the 1859 stormthe largest geomagnetic storm on record (Carrington, 1859;Cliver & Dietrich, 2013) to be approximately 9 V/km, and it is estimated that the nine-hour Hydro-Québec electricity blackout of March 1989 resulted from E fields of about 10 V/km (Barnes et al., 1991).
The predicted frequency dependences for the 99.97 th percentile of | | (denoted 99.97 ) take a very different form to those for the 100-year return levels (denoted 100 ): -field amplitudes at the 99.97 th percentile (occurring several times a year) are greatest for sinusoid periods of approximately 20 min, whilst 1/100 year events have greatest amplitude for periods between 30 s and 2 min. The observation that 100-year RL predictions vary greatly with sinusoid frequency has important implications when comparing and contrasting statistical studies evaluating extremes of |E| which may have been based on different sinusoid frequencies.

Accepted Article
This article is protected by copyright. All rights reserved.
Model estimates of the E field based on single-frequency components of the geomagnetic fluctuation have been reported by several authors (Beggan et al., 2013;Beggan, 2015;Bedrosian and Love, 2015;Love et al., 2016b). Love et al. (2016b) examined the amplitude of 4-min period sinusoids fitted to geomagnetic measurements (over sliding 10-min windows) and estimated extreme E-field amplitudes using empirical MT transfer functions at sites in the contiguous USA ( ≅ 40-60°N). Only at the northern limit, in the northern mid-west states, did they find 100 exceeding 3 V/km, which is similar to the 3-5 V/km predicted in Figure 16 for LER ( = 58°N) for a 4min sinusoid period. However, direct comparisons between sites cannot be made without considering differences in the surface impedance and its gradients. Bedrosian and Love (2015) illustrated this point by simulating the E fields generated by sinusoids with 10-, 100-, and 1000-s periods using MT transfer functions from the EarthScope MT array in the Midwest USA and showed that a constant-amplitude 0 = 500 nT, 100-s period B field would induce |E| of 2.7 V/km, averaged across all sites, but with values ranging from 0.15 to 16.8 V/km depending on site. Similarly, Pulkkinen et al. (2012), by extrapolating a log-normal distribution of 10-s field data from 23 European sites (55°-75°N geomagnetic), predicted 100 ranging from 5V/km with a high-conductivity ground model, to 20 V/km for poor-conductivity ground. Beggan et al. (2013) and Beggan (2015) also modelled the extreme E-field in the UK based on a conductivity model and B-fields modelled as sinusoids with periods, , of 2, 10 and 30 minutes and amplitudes based on the 30-, 100-, and 200-year return levels of 1-min ℎ / predicted by Thomson et al. (2011). The 2-min 100 prediction of Beggan et al. (2013) shown in their Fig. 6 (middle column) shows not only the high level of localisation of the E field intensity, ranging from around 2 to 7 V/km, but also the importance of the direction of the inducing B-field (whether N-S or E-W aligned) for some locations.
We intend to report further on the importance of directionality in extreme ℎ / statistics in a forthcoming publication.
There are, of course, limitations to 'narrowband' models of geomagnetic events since, in practice, fluctuations will be broadband in nature and the frequency spectrum of any individual geomagnetic event will be unique. We have noted that many of the extreme events (exceeding 99.97 ) identified in our dataset occur simultaneously (within hours of each other) over a wide range of timescales (or frequencies), but our results should not be used to infer a frequency spectrum of B or E fields for any given extreme geomagnetic event. For this information the reader may refer to several studies of extreme values that have taken the approach of analysing the E field produced during rare and intense geomagnetic storm periods and in some cases scaling up their effect to simulate 100-year return levels (e.g. Ngwira et al., 2013;Pulkkinen et al. 2012;Lotz & Danskin, 2017).

Conclusion
The importance of ULF waves in driving extreme geoelectric fields and GICs has received a great deal of interest in recent years (Hartinger et al., 2020;Belakhovsky et al. 2019;Heynes et al. 2020;Pulkkinen & Kataoka, 2006) and there is a need for better understanding of the frequency dependence of the B and E field fluctuations driving GICs (e.g. Pulkkinen et al., 2017). Most previous statistical climatological studies of extreme values for E and ℎ / have been based on sampling at just one or two frequencies. In this paper, however, we have presented statistics of large ( 99.97 ) and extreme (e.g. 1/100-year) values for | ℎ / | on a wide range of timescales, , from 1 to 60 min. At latitudes above the dayside cusp ( > 80°), for example, we find that occurrences of | ℎ / | ramp changes above 99.97 become tightly clustered in the few hours about local noon, and the effect is greatest

Accepted Article
This article is protected by copyright. All rights reserved.
for longer timescales ( ≥ 30 min). We have contrasted the statistics of ramp changes with those of the RMS of 1-min fluctuations over the same range of timescales and find, in particular, that in the auroral zone, for > 10 min the MLT of greatest occurrence of large RMS variation is from dawn to noon, indicative of strong ULF wave activity in this local time sector. The frequency (1/ ) dependences (for both ramp changes and RMS variations) are found to be not a simple power law, but are well modelled by quadratic functions whose three coefficients vary predictably with geomagnetic latitude.
For three UK locations we extended the data set to 1 Hz sampling frequency and, using a plane wave approximation and measured MT transfer functions, we derived the frequency dependence of the 99.97 th percentile and 100-year return levels of the geoelectric field, E at those sites. For events occurring several times a year (at the 99.97 th percentile) the induced E fields were greatest for fluctuations of 20-min period, whilst the 1-in-100year return levels were greatest for 0.5-2 min period fluctuations.
These statistics may be useful when inferring the likely extremes of | ℎ / | or E over a wide frequency range based on studies that used a single sampling cadence. The distributions of extreme occurrence rates with latitude, local time and season may also improve our understanding of the main ionospheric and magnetospheric drivers of GICs.

Annex A
It is here demonstrated that the horizontal geoelectric field magnitude |E| is proportional to as defined in Equation (2) when using a half-space model for the impedance matrix. From Equation (10), the idealised halfspace impedance is given by and for a sinusoidal horizontal magnetic field, ≡ ( ), with frequency (Hz), the horizontal geoelectric field is given by (Equation (13)) By comparing components it is observed that | | ∝ | / | (where ∝ denotes proportionality) and | | ∝ | / |, and so | | = √ 2 + 2 ∝ , where is defined using the expression in Equation (2). Approximate

Accepted Article
This article is protected by copyright. All rights reserved.
proportionality may be observed when the impedance Z differs only slightly from the half-space model of impedance in (A.1).

Accepted Article
This article is protected by copyright. All rights reserved. Sites in bold provided 1-s resolution data for this study.

Data Availability Statement
The 1-minute cadence magnetometer data used in this paper are available from https://supermag.jhuapl.edu and described in (Gjerloev, 2012). 1-second cadence UK magnetometer data are available from the British Geological Survey: http://www.geomag.bgs.ac.uk/data_service/data/home.html. The electric field data used to calculate the