© N.C. Rogers et al., Published by EDP Sciences 2020

1 Introduction

1.1 GICs and magneto-ionospheric currents

Large fluctuating magnetic fields arising from electrical currents in the ionosphere and magnetosphere can cause geomagnetically induced currents (GICs) in any ground-based infrastructure that contains long metal conductors. The damaging effects of GICs have been reported in relation to high voltage (HV) electricity power networks (Boteler et al., 1998; Pirjola et al., 2000; Erinmez et al., 2002; Molinski, 2002; Thomson et al., 2010; Boteler & Pirjola, 2017), trans-oceanic cables (Root, 1979; Lanzerotti et al., 1995), railway signalling systems (Wik et al., 2009; Eroshenko et al., 2010; Qian et al., 2016), railway electrification systems (Liu et al., 2016), and pipelines (Boteler, 2000; Pirjola et al., 2000; Pulkkinen et al., 2001). In many cases, these systems need to be engineered with resilience to the potentially catastrophic GIC damage that occurs only rarely, which can be aided by an assessment of the likelihood of such events over a period of 100 years or more.

The magnitude of GICs in an electricity network is dependent on the geoelectric field, E, induced by a nearby ionospheric or magnetospheric current, and may be calculated using a grid model of electrical impedances and the application of Ohm’s and Kirchhoff’s laws, and Thévenin’s theorem (e.g., Boteler & Pirjola, 2014, 2017). Routine measurements of the geoelectric field are, however, not globally extensive and do not extend over the decades required for accurate climatological prediction of 1/100-year return levels. Instead, we have taken as a proxy the rate of change of the horizontal geomagnetic field $\mathrm{d}{B}_H/\mathrm{d}t$, which, when combined with a model of the local ground conductivity, may be used to predict E (Cagniard, 1953). Digitised measurements of the geomagnetic field with 1-min resolution are readily available and extend over several decades.

To simplify the calculation of E, it is often assumed that both geomagnetic and geoelectric field perturbations result from a downward propagating plane wave and that the Earth conductivity is horizontally layered (see Pirjola, 2002 for a review). Under these assumptions the horizontal component of E follows the relation (Cagniard, 1953):

$$\left(\begin{array}{c}{E}_x\\ E_y\end{array}\right)=\left[\begin{array}{cc}{Z}_{\mathrm{xx}}\left(\omega \right)& Z_{\mathrm{yx}}\left(\omega \right)\\ Z_{\mathrm{xy}}\left(\omega \right)& Z_{\mathrm{yy}}\left(\omega \right)\end{array}\right]\left(\begin{array}{c}\mathrm{d}{B}_x/\mathrm{d}t\\ \mathrm{d}{B}_y/\mathrm{d}t \end{array}\right)$$where the Z_ij(ω) are elements of a frequency-dependent impedance matrix which varies with location (x, y) if the local geology is non-uniform.

An illustrative example of estimating GIC magnitudes from $\left|\mathrm{d}{B}_H/\mathrm{d}t\right|$ statistics was presented by Beggan et al. (2013). These authors applied the “thin-sheets” model of Vasseur & Weidelt (1977) to determine E at the Earth’s surface based on the (frequency-dependent) conductivity in the Earth crust in the United Kingdom. The magnetic field fluctuation was assumed sinusoidal with a fixed period, and the root-mean-square rate-of-change dB_H/dt (required as an input to the model) was taken from 100- and 200-year return value estimates from a study of European magnetometer records (Thomson et al., 2011). The extreme geoelectric field estimates so produced were then input into a model of the UK high-voltage electricity grid to determine the likely magnitude of GICs. This calculation required an assumption of the principal direction for the field fluctuation (e.g., North–South or East–West) to drive various hypothetical scenarios.

In this paper we extend the statistics of the high percentiles and projected extreme values of $|\mathrm{d}{B}_H/\mathrm{d}t|$ to a full global extent and provide a more detailed study of the direction of these fluctuations. We shall further show how the occurrence of extreme $|\mathrm{d}{B}_H/\mathrm{d}t|$ (and hence GIC) at a given latitude may have limited extent in local time and their prevalence may be confined to certain seasons – both factors that could be used to mitigate GIC risk. Fluctuation magnitudes of several hundred nT/min resulting from auroral current intensifications can cause HV electricity network outages, particularly if sustained over a long duration (Erinmez et al., 2002; Knipp, 2011, pp. 638–642) although other power grid impacts have been observed at low- and mid-latitudes even for fluctuations less than 100 nT/min (Kappenman, 2006 and references therein; Gaunt & Coetzee, 2007; Trivedi et al., 2007).

The location and timing of large fluctuations of the magnetic field is, of course, dependent on the climatology of extreme ionospheric and magnetospheric currents. The most intense GICs have been associated with intensifications of the Auroral Electrojets (AEJ), which connect to the magnetosphere via the field-aligned Birkeland currents (Milan et al., 2017). These are intensified during geomagnetically active periods such as substorms (e.g., Pulkkinen et al., 2003). The onset of auroral substorms occurs preferentially in the 20–24 magnetic local time (MLT) sector (Liou et al., 2001; Wang et al., 2005) and is followed by an expansion phase of about 25–40 min (Pothier et al., 2015) during which there are rapid intensifications of the westward electrojet currents flowing along auroral arcs at the equatorward edge of the auroral region, which bulges and expands poleward, eastward, and westward (Kivelson & Russell, 2005, pp. 421–429; Gjerloev et al., 2007; Gjerloev & Hoffman, 2014). The westward current enhancements are recorded as southward deflections of the magnetic field, as observed by ground-based magnetometers. To a good approximation, the direction of any fluctuation, dB_H, measured at the ground is 90° anticlockwise from the ionospheric current direction (as viewed from above the current sheet) assuming planar electromagnetic wave propagation and minimal horizontal gradients in ground conductivity (Viljanen et al., 2001; Belakhovsky et al., 2019).

An illustrative example of substorm geomagnetic activity is presented in Figure 1a which presents the geomagnetic North and East components of the magnetic field (the red and green lines, respectively) at the auroral magnetometer BRW (Utqiagvik, Alaska, 70° N, 109° W in Corrected Geomagnetic (CGM) coordinates; Baker & Wing 1989) and corresponding values for $|\mathrm{d}{B}_H/\mathrm{d}t|$ (blue line), where exceedances of the 99.97-percentile (P_99.97) are highlighted in yellow. The strong southeast-directed fluctuation dB_H at 07:07 – 07:10 UT could well have resulted from a strong southwest-directed electrojet current enhancement along a brightening auroral arc during substorm expansion.

Fig. 1 Three examples of large field fluctuations $|\mathrm{d}{B}_H/\mathrm{d}t|$. (a) An auroral substorm at BRW, Alaska, (b) an SSC at GUA, Guam, 6° N, 144° W CGM), and (c) Pc5 pulsations (0945–1015 UT) at NAQ, Greenland. The 99.97th percentile of $|\mathrm{d}{B}_H/\mathrm{d}t|$ is also shown.

Strong “Region 0” (R0) currents may also flow along magnetic field lines into the polar cusp region, driven by magnetic tension forces on field lines that have recently reconnected with the interplanetary magnetic field (IMF) at the dayside magnetopause (Milan et al., 2017). The cusp lies poleward of the Birkeland currents, typically around 77°–78° CGM latitude and at 12 ± 1.5 h MLT (Campbell, 2003, p. 139; Hunsucker & Hargreaves, 2003, p. 237), although under northward IMF conditions it covers a wider latitude region extending to higher latitudes (Newell et al., 1989). Large, transient (3–10-min) fluctuations of the vertical and horizontal geomagnetic field components have been observed for several decades at or just equatorward of the cusp region (Sibeck, 1993 and references therein). Many possible causes have been postulated, including “flux transfer events” (Russell & Elphic, 1978), enhanced dayside reconnection following a northward turning of the IMF (Pitout et al., 2001), surface wave instabilities such as Kelvin–Helmholtz (K–H) waves (Nykyri & Dimmock, 2016; Masson & Nykyri, 2018), and solar wind plasma injections penetrating the magnetopause (Menietti & Burch, 1988).

In many cases, there is a direct association between the high-latitude transients near the cusp region and geomagnetic storm sudden commencements (SSC) and sudden impulses (SI) (Lanzerotti et al., 1991; Sibeck, 1993; Sibeck & Korotova, 1996). Figure 1b is the magnetogram from a low-latitude magnetometer (Guam) and includes a sharp increase in B_N at 00:47 UT on 18 April 2001 characteristic of an SSC. Sudden commencements (SC) (a general term for both SSC and SI) occur in response to step changes in solar wind dynamic pressure that cause a sudden increase in the Chapman-Ferraro current flowing eastward along the dayside magnetopause boundary (Chapman & Ferraro, 1931; Milan et al., 2017). They typically have a magnitude of several tens of nT/min at low latitudes, but may occasionally exceed 120 nT/min (Fiori et al., 2014; Carter et al., 2015) and have been associated with GIC disturbances in HV electricity networks (Kappenman, 2003, 2004; Marshall et al., 2012; Fiori et al., 2014; Zhang et al., 2015; Belakhovsky et al., 2017, 2019).

Other large fluctuations in B_H arise from intense ultra-low frequency (ULF) geomagnetic pulsations, the largest of which are classed as “Pc5 pulsations” with periods of 2.5–10 min (e.g., Campbell, 2003, p. 168). Under typical mid-latitude conditions and moderate geomagnetic activity, Pc5 waves have amplitudes of approximately 70 nT (Campbell, 2003, p. 171). An example of such waves is evident at 09:45–10:15 UT in the magnetogram for NAQ (Narssarssuaq, Greenland, 66.2° N, 42.5° E CGM) presented in Figure 1c. The most powerful Pc5 pulsations are observed in the period between dawn and noon, and can arise from an Alfvén wave K–H instability, particularly during periods of high velocity solar wind (>500 km s⁻¹) (Engebretson et al., 1998; Vennerstrøm, 1999; Pahud et al., 2009). Alternatively, Pc5 waves may be triggered by shocks in the solar wind associated with Sudden Commencements (Zong et al., 2009; Zhang et al., 2010; Hao et al., 2019).

Belakhovsky et al. (2019) also recently considered the importance of irregular Pi3 ULF waves and Travelling Convection Vortices (TCV) in driving large $|\mathrm{d}{B}_H/\mathrm{d}t|$ events. Pi3 waves are described as “quasi-periodic sequences of magnetic impulses” with timescales between 10 and 20 min, and TCVs have been observed as a twin vortex of Hall currents in the polar cleft region. TCVs are generated by a pair of upward and downward field-aligned currents that can be triggered by a sudden change in the solar wind (Friis-Christensen et al., 1988; Vorobjev et al., 1999; Engebretson et al., 2013).

1.2 Modelling the probability of extreme field fluctuations

The probability distribution of non-extreme values of $|\mathrm{d}{B}_H/\mathrm{d}t|$ at any given location is well approximated by a lognormal distribution (Love et al., 2016; Manpreet, 2018). This distribution may not be suitable for extremes since its purpose is to produce a good fit to the main body of the distribution rather than the tails, and it may not extrapolate well beyond the observed range of the data. Because extreme events are, by definition, rare, there is little evidence from the data with which to validate the appropriateness of the fit of a lognormal to the tail.

To model the tail distribution, extreme value theory (EVT) has been applied (Coles, 2001; Reiss & Thomas, 2007) and a Generalised Pareto (GP) distribution fitted to observations of $|\mathrm{d}{B}_H/\mathrm{d}t|$ above a high threshold. By modelling how the site-specific parameters of the GP distributions vary globally, we shall present model-based predictions of the “return value” of $|\mathrm{d}{B}_H/\mathrm{d}t|$ at any location on the Earth for “return periods” ranging from 5 to 500 years. We have also examined the probability of $|\mathrm{d}{B}_H/\mathrm{d}t|$ threshold exceedances as a function of month, MLT, and the direction of the fluctuation, dB_H. This additional information helps to associate extreme magnetic field fluctuations with their causative extreme ionospheric and magnetospheric electrical currents. Modelling patterns of occurrence probability in this way also improves the modelling of risk when engineering GIC resilient ground infrastructure, noting that electricity usage or signalling error requirements may vary with time-of-day and season (e.g., Molinski, 2002). Individual linear elements of pipelines or electricity cables are also sensitive to fluctuations resolved along only one directional axis.

In previous studies, EVT has been used to examine the likelihood of extreme values of geomagnetic indices including D_st (Silbergleit, 1996; Tsubouchi & Omura, 2007), the Auroral Electrojet indices (AU, AL, and AE) (Nakamura et al., 2015), A_p (Koons, 2001), and the aa and AA* indices (Siscoe, 1976; Silbergleit, 1999), although the location and timing of large $|\mathrm{d}{B}_H/\mathrm{d}t|$ events are not well predicted by geomagnetic index statistics (e.g., Kozyreva et al., 2018). In related space weather fields, EVT has characterised the probability of extreme solar flare X-ray flux (Elvidge & Angling, 2018; Tsiftsi & De la Luz, 2018) and extreme high-energy (>2 MeV) radiation-belt “killer” electron fluxes (Koons, 2001; O’Brien et al., 2007; Meredith et al., 2015).

EVT and other statistical approaches have been applied to the prediction of extreme $|\mathrm{d}{B}_H/\mathrm{d}t|$ in several previous publications (Thomson et al., 2011; Love et al., 2016; Nikitina et al., 2016; Wintoft et al., 2016) and individual case studies have examined extreme $|\mathrm{d}{B}_H/\mathrm{d}t|$ characteristics (e.g., Viljanen et al., 2001; Pulkkinen et al., 2012; Ngwira et al., 2013, Kozyreva et al., 2018). In Section 3 of this paper, we present a new assessment based on a significantly larger dataset than presented in these earlier studies – 1.9 billion field measurements from 125 magnetometers worldwide. This is followed by an analysis of the latitudinal, seasonal, MLT, and directional dependences of the extreme events, presented as occurrence probability distributions for $|\mathrm{d}{B}_H/\mathrm{d}t|$ exceeding the 99.97th percentile (P_99.97). In Section 4 we shall (i) relate our new findings to the earlier studies cited above, (ii) further assess the importance of SC in triggering impulses at high latitudes, and (iii) show how the occurrence probability profiles change with IMF orientation. The paper concludes with a summary of these findings.

2 Data and methods

2.1 Magnetograms

The magnetic field measurements for this study were obtained through SuperMAG, an interface to magnetic measurements made by a global network of geophysical institutes; SuperMAG being a project led by Johns Hopkins University (Gjerloev, 2009). Gjerloev (2012) described the processing of the magnetometer measurements in the SuperMAG data set which may be summarised as follows: The field measurements (1-min averages) were rotated into a local NEZ geomagnetic coordinate system (N = north, E = east, Z = down) using a time-dependent local declination angle obtained using a 17-day sliding window smoothing filter. For each magnetometer, the N, E, and Z components were then individually baselined as follows: (i) Diurnal variations of the field (mainly due to the S_q current system: Matsushita, 1968; Yamazaki & Maute, 2017) were removed by subtracting a diurnal trend, (ii) the secular variation of the Earth’s main field was removed by subtracting the linear trend over each year (with some contribution from neighbouring years), and (iii) the residual scalar offset was subtracted. This method of processing preserves the short-term field fluctuations associated with substorms, as well as those due to geomagnetic storms. Long-period (>1 month) variations in the declination angle are, however, removed by the rotation into local NEZ coordinates.

The magnetometer records in the SuperMAG archive had already been cleaned and manually inspected to remove most sudden changes in the baseline (offsets), spikes, and gradual slopes (Gjerloev, 2012). However, for this study – since the statistics of extreme and rare events may be significantly influenced by a small number of imperfections – all the data in weeks containing $|\mathrm{d}{B}_H/\mathrm{d}t|$ peaks above the 99.97th percentile (P_99.97) were further visually inspected and obvious artefacts replaced by data gaps. Examples of common artefacts requiring correction were: (i) large, isolated, 1–2 min-duration spikes in signal levels (B_N and/or B_E), (ii) large step-changes in signal level, often occurring at the 00:00 UT boundary and sometimes followed by a correction to previous levels after an interval of perhaps 30–60 min, (iii) an obvious saturation in signal level, which occurred for B_N ≈ −1000 nT at several auroral magnetometers on several occasions.

In relation to GICs, the North and East components of the induced geoelectric field are approximately proportional to (Wintoft et al., 2016):

$$\left(E_{\mathrm{N}},E_{\mathrm{E}}\right)=k\left(\frac{\mathrm{d}{B}_{\mathrm{E}}}{\mathrm{d}t}, \hspace{1em}-\frac{\mathrm{d}{B}_{\mathrm{N}}}{\mathrm{d}t}\right)$$(1)where B_N and B_E are the Northward and Eastward components of the geomagnetic field and k is a constant that depends on the local ground conductivity. Equation (1) is valid under the assumptions that the magnetic disturbance propagates as a plane wave and the ground conductivity has no horizontal gradients. We therefore define $|\mathrm{d}{B}_H/\mathrm{d}t|$ as:

$$\left|\frac{\mathrm{d}{B}_H}{\mathrm{d}t}\right|=\frac{1}{\Delta t}\sqrt{{\left(\Delta B_{\mathrm{N}}\right)}^2+{\left(\Delta B_{\mathrm{E}}\right)}^2}$$(2)where $\Delta B_{\mathrm{N}}=B_{\mathrm{N}}\left(t+\Delta t\right)-B_{\mathrm{N}}\left(t\right)$ and $\Delta B_{\mathrm{E}}=B_{\mathrm{E}}\left(t+\Delta t\right)-B_{\mathrm{E}}\left(t\right)$. This definition (which was also adopted by Wintoft et al., 2015, 2016; Falayi et al., 2017; Ngwira et al., 2018; Kozyreva et al., 2018 and others) ensures that statistics of the induced E-field magnitude, $\left|E\right|=\sqrt{E_{\mathrm{N}}^2+E_{\mathrm{E}}^2}$ required for GIC modelling will be directly proportional to the distribution of $|\mathrm{d}{B}_H/\mathrm{d}t|$ given by equation (2). The compass direction of dB_H/dt was also recorded as:

$$D=\arg \left(\Delta B_{\mathrm{N}}+i\Delta B_{\mathrm{E}}\right).$$(3)

There are 130 magnetometer sites in the SuperMAG database for which at least 20 years of data were available. These covered the period from 1 January 1969 to 31 December 2016 (albeit with some data gaps). Five stations at low-latitude locations in the Atlantic longitude sector (HUA, BNG, KOU, MBO, and TAM, shown unshaded in Fig. 2) were excluded from analysis since their locations in the corrected geomagnetic (CGM) coordinate system are undefined (see Fig. 7 of Laundal & Richmond, 2017). Average geomagnetic CGM coordinates were calculated using the software described by Shepherd (2014) and are listed together with site codes, names, and geodetic coordinates in Supplementary Material (folder S1). Figure 2 presents a map of these locations with colours representing P_99.97. Immediately apparent from Figure 2 is the maximum in P_99.97 in the auroral zones (approximately 55°–75° CGM latitude) in both hemispheres, although the density of sites is much greater in the northern hemisphere.

Fig. 2 Magnetometer sites in the SuperMAG collaboration with ≥20 years data, coloured by P99.97 at each site. Contours show CGM latitudes for IGRF model epoch 2000.

2.2 Fitting the GP tail distribution

Taking a sequence of independent and identically distributed (IID) random variables, $\mathit{x}=\left[x_1,x_2,\dots , x_n\right]$, and a sufficiently high threshold, u, it may be shown that the cumulative distribution function of threshold exceedances y = x−u, conditional on x > u, takes the approximate form of a Generalised Pareto (GP) distribution (Coles, 2001, p. 75) given by:

$$H\left(\mathit{y}\right)=1-{\left(1+\frac{\xi \mathit{y}}{\sigma }\right)}^{-\frac{1}{\xi }}$$defined on {y: y > 0 and (1 + $\xi$y/σ) > 0}, where σ > 0 is the scale parameter, and $\xi$ the shape parameter of the distribution. Note that in the limiting case $\xi$ → 0, $H$(y) → 1−exp(−y/σ), that is, the distribution of y is exponential. Where $\xi$ < 0, the distribution is bounded to an upper limit of u−σ/$\xi$, whilst if $\xi$ > 0 the distribution is heavy-tailed, with no upper limit. The occurrence probability of extreme x therefore increases with $\xi$ and/or σ.

Because the above model is based on an asymptotic result that only holds exactly as u → y₊, where:

$$y_{+}=\min \left\{\mathit{y}: H\left(\mathit{y}\right)=1\right\}$$is the upper end-point of the distribution of y, the threshold u must be set high enough that the approximation for finite u can be assumed to be valid, i.e. so the distribution fits the observed threshold exceedances well. Conversely, u should not be so high that the sample size for modelling is too small. A 99.97-percentile threshold was adopted in this study following the assessment of Thomson et al. (2011) based on 28 European magnetometers, and the case study of Nikitina et al. (2016) who suggested thresholds in the region 99.93–99.97% for the site VIC (Victoria, Canada). Based on the assumption that threshold exceedances are a set of independent observations drawn from a GP distribution with common parameters (σ, $\xi$), these parameters can be estimated by maximisation of an appropriate likelihood function (Pawitan, 2001).

Occurrences of peaks over the threshold ($|\mathrm{d}{B}_H/\mathrm{d}t|$ > u) often exhibit short-range temporal dependence, with tight clustering of peaks (e.g., during a geomagnetic storm). This violates the assumption of an IID random sequence required in fitting the GP distribution as described above. The easiest way to deal with this is to ensure temporal independence by extracting independent clusters of extreme events and modelling only the maxima of these clusters. Consequently, the $|\mathrm{d}{B}_H/\mathrm{d}t|$ data were declustered using a “run-length below threshold” method such that any exceedances of P_99.97 within 12 h of the previous exceedance were considered to be part of the same cluster, and only the largest $|\mathrm{d}{B}_H/\mathrm{d}t|$ in each cluster was recorded. Thomson et al. (2011) discussed the merits of alternative declustering thresholds and run lengths, but concluded that the 99.97th percentile and 12-h run-length were a reasonable compromise.

A measure of the short-range temporal dependence in the peaks $\left(\left|\frac{\mathrm{d}{B}_H}{\mathrm{d}t}\right|>u\right)$, called the “extremal index”, was empirically determined as (Coles, 2001, p. 103):

$$\widehat{\theta }=n_c/n_u$$where n_c is the total number of clusters, and n_u is the total number of peaks above the high threshold, u. The index $\widehat{\theta }$ may be interpreted loosely as the reciprocal of the mean number of peaks in each cluster.

The maximum likelihood estimates (MLE) of GP parameters $\xi$ and σ at each site were determined, together with 95% confidence intervals (CI) based on the asymmetric likelihood profiles. Following (Coles, 2001, p. 78ff) diagnostic tests were then conducted to ensure that if the threshold were to be increased to u* (for all u* > u) then both the shape, $\xi$, and the “modified scale”, σ* = σ−$\xi$u*, would remain constant, and the mean excess of the tail distribution, $E\left(\mathit{x}-u^{\mathrm{*}} \right| \mathit{x}>u^{\mathrm{*}})$, would scale in proportion to u*, as required. The empirical and model tail distribution functions were then visually compared to ensure that the model was a good fit to the measurements at the new “manually revised” threshold.

The m-observation return level, x_m (nT/min) was calculated from MLE GP parameters as (Coles, 2001, p. 103):

$$x_m=u+\frac{\sigma }{\xi }\left[{\left(m\widehat{\theta }\widehat{{\zeta }_u}\right)}^{\xi }-1\right]$$(4)where $\widehat{{\zeta }_u}$ is the probability of threshold exceedance, determined empirically as:

$$\widehat{{\zeta }_u}=\frac{n_u}{n}$$and n is the total number of measurements. As an example, for a 100-year return period and 1-min measurements setting m = 100 × 365.25 × 24 × 60 = 52,596,000 in equation (4) yields the 100-year return level.

For each cluster peak ($|\mathrm{d}{B}_H/\mathrm{d}t|$ > u) we recorded the magnetic local time (MLT), defined as (Laundal & Richmond, 2017):

$$\mathrm{MLT}=\frac{\left(\varphi - {\varphi }_{\mathrm{cd}, \mathrm{ss}}\right)}{15\mathrm{°}}+12$$where ϕ_cd,ss is the centred dipole longitude of the sub-solar point, in degrees, (determined using the formulas in Appendix C of Laundal & Richmond, 2017) and ϕ is the geomagnetic longitude of a magnetometer in CGM coordinates. ϕ_cd,ss was determined as:

$${\varphi }_{\mathrm{cd}, \mathrm{ss}}={\mathrm{\Phi }}_{\mathrm{geo},\mathrm{ss}}{- \mathrm{\Phi }}_{\mathrm{N}}$$where Φ_geo,ss is the geographic longitude of the sub-solar point and Φ_N is the centred dipole longitude of the North Pole, determined as:

$${\mathrm{\Phi }}_{\mathrm{N}}=\arg \left({ -g}_1^1-i{h}_1^1\right)$$where $g_1^1,h_1^1$ are Gauss coefficients of the International Geomagnetic Reference Field (IGRF) (Thébault et al., 2015), linearly interpolated in time. The CGM longitude, ϕ, was determined for each site at 1-year intervals using the software described by Shepherd (2014) and then linearly interpolated onto a 1-min-resolution time scale.

IMF components B_y (duskward) and B_z (northward) were also recorded in geocentric solar magnetic (GSM) coordinates (Laundal & Richmond, 2017) for each cluster peak. The IMF measurements were interpolated from 1-h average values provided in the NASA Goddard Space Flight Center’s “OMNI” database (https://omniweb.gsfc.nasa.gov). These measurements from spacecraft at the L₁ Lagrange point have been time-shifted to account for the propagation delay to the magnetospheric bow shock nose using the method of Weimer et al. (2003).

3 Results

3.1 Latitudinal profiles of extreme $|\mathrm{d}{B}_H/\mathrm{d}t|$

Figure 3a presents P_99.97 for the 1-min $|\mathrm{d}{B}_H/\mathrm{d}t|$ as a function of the mean absolute CGM latitude, $\left|\lambda \right|$, for each magnetometer site (the mean is presented since the latitude can vary a few degrees over the period of measurements). The blue crosses in Figure 3a represent P_99.97, whilst the red circles indicate the manually revised thresholds, following the diagnostic tests described above. The outlier at λ = 70° N is for the site at Utqiagvik, (formerly Barrow), Alaska (BRW), which together with nearby sites in Alaska presents anomalously high values for $|\mathrm{d}{B}_H/\mathrm{d}t|$ (Fig. 2). Viljanen et al. (2001) speculated that local anomalies in the strength and direction of the fluctuations (as observed at four out of 25 magnetometers in Fennoscandia) could result from regions of anomalously high conductivity in the Earth’s crust since telluric currents induced in the crust contribute to the measurement of $\mathrm{d}{B}_H/\mathrm{d}t$ at the magnetometer. Kozyreva et al. (2018, p. 10) also noted that coastal stations BRW and AND produced anomalously high $|\mathrm{d}{B}_H/\mathrm{d}t|$ during the intense March 17, 2015 geomagnetic storm. Measurements from these magnetometer sites have, however, been retained in the analysis to avoid any subjective bias in the selection of data.

Fig. 3 GP parameters for 125 magnetometers versus absolute CGM latitude. (a) 99.97-percentiles of $|\mathrm{d}{B}_H/\mathrm{d}t|$ and revised thresholds, u, (b) MLEs of scale, σ (error bars indicate 95% CIs), (c) MLEs of shape, $\xi$, (d) extremal index, $\widehat{\theta }$, (e) number of measurements, n, number exceeding P99.97, nu, number of clusters, nc, and mean nc per year. Red curves indicate fitted smoothed splines (5).

The MLEs for the fitted scale and shape parameters are presented in Figure 3b and c respectively, with error bars representing the 95% confidence intervals. Northern hemisphere (NH) stations are presented in black and southern hemisphere (SH) stations in blue, showing that there are no significant hemispherical differences in the CGM latitude profiles. There is a substantial difference in the functional form between the scale, $\sigma \left(\left|\lambda \right|\right)$, and the threshold, $u\left(\left|\lambda \right|\right)$: There is a positive skew in $\sigma \left(\left|\lambda \right|\right)$ that is not apparent in $u\left(\left|\lambda \right|\right)$, which is more evenly distributed about 67°. Similarly, the functional form of the shape parameter profile,$\xi \left(\left|\lambda \right|\right)$, (Fig. 3c) exhibits a sharp peak at 53° CGM latitude, a broad minimum in the auroral zone (centred about 70°) and a clear increase in value towards the geomagnetic poles. The extremal index, $\widehat{\theta }$, presented in Figure 3d indicates strong latitudinal variation in the level of short-range temporal dependence (or level of clustering). Figure 3e presents the number of measurements for each site (approximately 10–25 million), the number for which $|\mathrm{d}{B}_H/\mathrm{d}t|$ exceeded P_99.97, and the number of clusters for each site (from 146 to 1084), and the mean number of clusters per year of data (from 5 to 56).

The MLE parameters presented in Figure 3 were substituted into (4) to determine the MLE values of return levels presented for each site in Figure 4. The sharp rise in return level at latitudes near $\left|\lambda \right|$ = 53°, typically considered to be sub-auroral under quiet to moderate geomagnetic conditions, arises from large values of both the σ and $\xi$ parameters at this latitude.

Fig. 4 Return levels of $|\mathrm{d}{B}_H/\mathrm{d}t|$ (MLEs, with 95% CIs shown as error bars), for return periods of (a) 100 years, (b) 200 years, and (c) 500 years with fitted smoothed splines ($s\left(\left|\mathrm{\lambda }\right|\right)$ in (5)). Panel (d) shows spline fits for a range of return periods.

The curves in Figures 3 and 4 represent smoothing spline functions, s, found by minimising:

$$p\sum _i{\left(z_i-s\left(\left|{\lambda }_i\right|\right)\right)}^2+\left(1-p\right)\int {\left(\frac{{\mathrm{d}}^{2}s}{\mathrm{d}{\lambda }^2}\right)}^2\mathrm{d}\left|\lambda \right|$$(5)where z_i are the ordinates (e.g., the MLE GP parameter at each site, i), and p is a smoothing parameter. The choice of p = 0.01 provided a good balance between goodness of fit (the first term in (5)) and function smoothness (the second term of (5)). The numerical minimisation algorithm used was based on the MATLAB™ function csaps. A MATLAB™ script to regenerate all of the spline curves of Figures 3 and 4 is provided in Supplementary Material (folder S2).

3.2 Probability of occurrence versus latitude and MLT

When modelling return levels, the probability of threshold exceedance $\widehat{{\zeta }_u}$ in (4) may be modelled as a function of observed covariates, i.e. $\widehat{{\zeta }_u}=\widehat{{\zeta }_u}({\alpha }_1, {\alpha }_2,\dots ),$ where α_i are site-specific, possibly time-varying, covariates of the data (e.g., MLT, month, direction, etc.). This additional information allows a prediction of GIC risk to be refined for systems operating in a restricted range of local times, times of year, etc., or for individual elements of pipeline or cable network with a fixed directional orientation. The latitudinal, seasonal, and diurnal patterns in the occurrences of $|\mathrm{d}{B}_H/\mathrm{d}t|$ exceeding P_99.97 may be used as an approximation to $1-\widehat{{\zeta }_u}({\alpha }_1, {\alpha }_2,\dots )$ for u ≈ P_99.97, and these patterns are also useful in identifying the general ionospheric/magnetospheric current systems driving the extreme geomagnetic fluctuations.

Figure 5a presents the number of cluster peaks for all magnetometers, binned by $\left|\lambda \right|$ and MLT, i.e.,

$$N\left(\left|\lambda \right|,\mathrm{MLT}\right)=\sum _s{n}_{c,s}(\left|\lambda \right|,\mathrm{MLT})$$where $n_{c,s}(\left|\lambda \right|,\mathrm{MLT})$ is the number of cluster peaks in each ($\left|\lambda \right|$, MLT) bin for site s (latitude bins containing no magnetometers are shown in grey). The summation of cluster peaks over multiple sites has the effect of smoothing the distribution and reducing quantisation noise where $n_{c,s}(\left|\lambda \right|,\mathrm{MLT})$ may be small. Since the number of magnetometer measurements varies with latitude, Figure 5b represents the same data after normalisation by the number of cluster peaks recorded in each bin. This is a close approximation to the 2-d distribution of occurrence likelihood for cluster peaks exceeding P_99.97 and it was calculated as:

$$P\left(\left|\lambda \right|,\mathrm{MLT}\right)=\frac{N\left(\left|\lambda \right|,\mathrm{MLT}\right)}{\sum _S{n}_{c,s}\left(\left|\lambda \right|\right)}\times \frac{\Delta \lambda }{90\mathrm{°}}$$where $n_{c,s}\left(\left|\lambda \right|\right)$ is the total number of cluster peaks for a site s that lies within the latitude bin of width Δλ, and we have assumed the observations are uniformly distributed in MLT. Note that, where Δλ is chosen to prevent data gaps,

$$\sum _{\left|\lambda \right|}\sum _{\mathrm{MLT}}P\left(\left|\lambda \right|, \mathrm{MLT}\right)=1.$$

Fig. 5 Total occurrences of cluster peaks binned by geomagnetic latitude $\left|\lambda \right|$ and MLT. (a) Number of clusters, (b) probability distribution of cluster peaks, (c) fitted spherical harmonic expansion (Psph in (6)).

Figure 5b shows that at auroral latitudes (55°–80° CGM) the cluster peaks occur most often at 20–24 MLT, which is the time sector associated with the greatest number of substorm onsets. A cursory inspection of magnetometer records associated with these peaks indicated that many were indeed associated with substorm activity (similar to that in Fig. 1a).

The second most prominent region of high occurrence is at 10–14 MLT above 77° CGM latitude, i.e. the region poleward of the dayside cusp. It should be emphasised that since Figure 5 relates to percentiles P_99.97 that are themselves a function of $\left|\lambda \right|$ (as presented in Fig. 3a) the magnitudes of the $|\mathrm{d}{B}_H/\mathrm{d}t|$ peaks in this high-latitude region are much lower than in the auroral and sub-auroral zones. Nonetheless, the clustering of events about 12 MLT is very clear.

A third region of raised occurrence likelihood is observed in Figure 5b with a well-defined increase in the MLT locus of cluster peaks from 3 to 10 MLT as CGM latitude increases from 60° to 75°. This pattern aligns with patterns of Pc5 wave occurrence reported in the literature: For example, Baker et al. (2003) observed an increase in Pc5 wave power from 7 to 9 MLT as λ increased from 65° N to 74° N, Pahud et al. (2009) showed Pc5 power upper quartiles increasing from 4 to 8 MLT between 61° N and 69° N CGM under conditions of high solar wind speed, and the average wave power distributions in Figure 5 of Vennerstrøm (1999) show a clear pattern of wave power increasing from approximately 67°–75° N invariant latitude over the 4–12 MLT range with a secondary maximum around local midnight (22–2 MLT) attributed to substorms. The effect was also noted in $|\mathrm{d}{B}_H/\mathrm{d}t|$ measurements by Viljanen et al. (2001) for three sites in Scandinavia. Various theories for the asymmetry in Pc5 wave activity about local noon were discussed by Pahud et al. (2009, and references therein) and include a dawn-dusk asymmetry in the wave excitation processes in the magnetosphere, differences in the harmonic modes and polarisations of the waves, and differences in the level of screening by the conductive ionospheric layer.

At latitudes below 50° CGM, there is a distribution of cluster peaks between 6 and 18 MLT forming a V-shaped distribution centred about 11 MLT, with a reduced spread of local times at the lowest latitudes. This is broadly consistent with enhancement in SC amplitudes between 8 and 16 MLT (maximising at 11 MLT) near the equator reported by Shinbori et al. (2009) and Russell et al. (1994) and others. These authors also reported a broad enhancement in SC amplitude in the midnight sector attributed to strong cross-magnetotail and field-aligned currents, but no enhancement of midnight sector cluster peak frequency is evident in Figure 5 indicating that most$\left|\mathrm{d}{B}_H/\mathrm{d}t\right|$ cluster peaks in this local time sector have magnitudes below the 99.97th percentile. Further investigations are required to identify the physical processes associated with the V-shaped occurrence distribution, which has maxima at approximately 8 and 15 MLT for absolute CGM latitudes between 20° and 40°.

The empirical data values from Figure 5b could be applied directly in a model of $\mathbf{1}-\widehat{{\mathit{\zeta }}_{\mathit{u}}}\left(\left|\mathit{\lambda }\right|, \mathbf{MLT}\right)$ with interpolation and extrapolation in regions of missing data. As an alternative, a 2-d spherical harmonic series expansion has been fitted to the data producing the smooth model surface shown in Figure 5c. The spherical harmonic expansion (6) was fitted by regression to $\mathit{P}\left(\left|\mathit{\lambda }\right|,\mathbf{MLT}\right)$ and normalised to ensure a non-negative model probability distribution that integrated to unity:

$$P_{\mathrm{sph}}\left(\lambda , \mathrm{MLT}\right)=\sum _{l=0}^L\sum _{m=-\mathrm{l}}^l{A}_{\mathrm{lm}}{P}_l^m\left(\sin \lambda \right)e^{\mathrm{im}\pi \mathrm{MLT}/12}$$(6)where $P_l^m$ are Schmidt semi-normalized associated Legendre polynomials, and MLT is in hours. A minimum order of L = 19 was required to adequately represent the structure in distribution. The 400 complex A_lm coefficients are provided as Supplementary Material (folder S3), together with example MATLAB™ code to reproduce $P_{\mathrm{sph}}\left(\lambda , \mathrm{MLT}\right)$.

3.3 Seasonal, MLT, and directional distributions

The patterns of cluster peaks with month and MLT are presented in Figure 6. Each panel presents, using a subset of magnetometer sites in a particular band of CGM latitudes, the distribution:

$$P\left(\mathrm{month}, \mathrm{MLT}\right)=\frac{\sum _s{n}_{c,s}(\mathrm{month},\mathrm{MLT})}{\sum _s{n}_{c,s}}$$where n_c,s (month, MLT) is the number of cluster peaks in each (month, MLT) bin for site s, and n_c,s is the total number of clusters for site s. Since the number of field measurements in each (month, MLT) bin is approximately equal, the value P (month, MLT) approximates the 2-d probability distribution of cluster peak occurrences. Figure 7 presents, in a similar format to Figure 6, the distribution:

$$P\left(D, \mathrm{MLT}\right)=\frac{\sum _s{n}_{c,s}\left(D,\mathrm{MLT}\right)}{\sum _s{n}_{c,s}}$$where D is the compass direction of the large geomagnetic fluctuation (3). Note in each case that:

$$\sum _{\mathrm{month}}\sum _{\mathrm{MLT}}P\left(\mathrm{month}, \mathrm{MLT}\right)=\sum _D\sum _{\mathrm{MLT}}P\left(D, \mathrm{MLT}\right)=1.$$

Fig. 6 Occurrence probability of $|\mathrm{d}{B}_H/\mathrm{d}t|$ > P99.97 (after declustering) as a function of month and MLT. Panels (a)–(c) are for NH stations, whilst (d)–(f) represent SH stations. (a) and (d) represent low-latitude stations (<40° CGM), (b) and (e) represent auroral latitudes 55°–75° CGM) and (c) and (f) represent latitudes above 77° CGM.

Fig. 7 Occurrence probabilities in the same format as Figure 6 but shown as a function of dBH direction, D, and MLT. Bin sizes are 1 h × 10°.

Panels a–c in Figures 6 and 7 present the occurrence distribution for NH magnetometer sites (λ > 0), whilst panels d–f are for SH sites (λ < 0). Panels a and d in both figures were produced from a subset of low latitude sites ($\left|\lambda \right|$ < 40°) and display a slight increase in occurrence likelihood on the day side, between 7 and 17 MLT, with little seasonal variation. The associated direction distribution in Figure 7a and d shows a strong “preferred direction” for these large events, principally towards the north (D ≈ 0°) at most times of day, with a secondary maximum for southward fluctuations (D ≈ 180°). A large proportion of the events are associated with SC, which is characterised by a large northward-directed dB_H impulse. However, in the 6–10 MLT region, the directional distribution P(D, MLT) is more isotropic. Repeating the analysis with no declustering produces very similar distributions for P(D, MLT) (Figure S1 in Supplementary Material (folder S4)), albeit with less quantisation noise. This indicates that the choice of declustering parameters does not affect P(D, MLT) substantially.

Panels b and e of Figure 6 represent a subset of stations in a band of auroral latitudes (55° < |λ| < 75°) and present a marked increase in the rate of occurrence in the pre-midnight hours (20–24 MLT) maximising near the equinoxes, with a secondary maximum in the winter months centred approximately 1 h later at 22–24 MLT. This is consistent with the 1-h difference in substorm onset MLTs between summer and winter reported by Liou et al. (2001). Geomagnetic activity is generally increased near the equinoxes when the geomagnetic field is more favourably oriented for reconnection with the IMF (Russell & McPherron, 1973; Zhao & Zong, 2012) and this effect has been observed in several previous studies of geomagnetic fluctuations (Boteler et al., 1998; Viljanen et al., 2001; Beamish et al., 2002). The same panels in Figure 7 show that the principal direction of the dB_H in this MLT sector is southward (D ≈ 180°), which indicates that strong westward auroral electrojet currents are the dominant driver of such fluctuations.

The period 3–11 MLT in Figure 6b and e is also associated with raised occurrence probability and follows the known patterns of enhanced ULF wave activity described in Section 3.2. In the NH (Fig. 6b) the region of higher occurrence likelihood is centred at approximately 4 MLT near the June solstice but increases to 9 MLT towards the December solstice. The seasonal pattern in the southern hemisphere is less clear due to the smaller number of observations in this region. The corresponding panels of Figure 7 show that the direction of these changes in the NH is broadly aligned to the axis ESE – WNW (100°–280°) from 3–7 MLT, but move closer to the N–S axis for events near the December solstice at MLT 7–11. Again the trends in the Southern hemisphere are less clear due in part to the smaller number of magnetometer stations.

Panels c and f in Figure 6 were produced from the subset of magnetometers with $\left|\mathit{\lambda }\right|$ > 77° (i.e., at or poleward of the cusp-region). Near the winter solstices, the maximum occurrence likelihood is centred near 23–24 MLT and may be related to substorm activity in this sector. The corresponding panels of Figure 7 show that the associated direction is predominantly southward for NH stations (indicative of substorm-related westward currents), although curiously this is not the case for SH stations. In the summer months, there is a strong clustering of events in the few hours about 11–12 MLT, which may be attributed to R0 currents and ULF waves and transients in this region. Previous studies have shown that the magnitude of the R0 current has a strong seasonal variation, peaking at the summer solstice (Wang et al., 2008; Milan et al., 2017). Under conditions of northward IMF at the summer solstice near noon, the large dipole tilt creates an “overdraped” magnetospheric tail lobe with a magnetic field topology favourable to field-line merging with the IMF (Crooker, 1992; Watanabe et al., 2005). Such reconnections may drive strong impulsive field-aligned currents into the region that is poleward of the dayside cusp and may explain the prevalence of extreme $|\mathbf{d}{\mathit{B}}_{\mathit{H}}/\mathbf{d}\mathit{t}|$ under these conditions.

4 Discussion

4.1 Predictions of extreme $|\mathrm{d}{B}_H/\mathrm{d}t|$

The Thomson et al. (2011) study of 1-min $|\mathrm{d}{B}_H/\mathrm{d}t|$ on European magnetometers (40°–74° N) predicted that, for a few stations at 55°–59° N geomagnetic (dipole) latitude, 100-year return levels would increase to approximately 4000 nT/min, with 200-year return levels increasing to 6000 nT/min. The highest return levels were predicted for sites BFE and ESK, which are more closely separated in CGM latitude than in dipole latitude, laying between 52° N and 53° N CGM. This was broadly supported by the Wintoft et al. (2016) analysis of 14 European stations (48°–67° N geomagnetic). The return levels (RL) of the present study (Fig. 4), 100-year (200-year) RL (MLE) up to 2885 nT/min (4553 nT/min) at site BFE (λ = 52.0° N), are broadly consistent with these earlier studies, but extend to a full range of latitudes over both hemispheres and more clearly distinguish a sharp skewed peak at $\left|\lambda \right|$ = 53° in both hemispheres. The Love et al. (2016) study of 34 magnetometers worldwide also predicted a sharply peaked distribution of 100-year d|B_H|/dt RLs, maximising at a higher geomagnetic latitude of 60° between 1000 and 2000 nT/min (see their Figure 4d), whilst the Pulkkinen et al. (2012) study of data from two intense geomagnetic storms (March 1989 and October 2003) postulated that there is a sharp “latitude threshold” around 50°–55° geomagnetic latitude, below which the magnitude of |dB_H/dt| dropped by an order of magnitude, which is also observed in Figure 4. A subsequent analysis of twelve large geomagnetic storms by Ngwira et al. (2013) confirmed that this latitude threshold was associated with the limited equatorward expansion of the Auroral Electrojet (AEJ). The most extreme ionospheric currents are likely to be associated with substorm expansions along auroral arcs near the edge of the auroral bulge and this interpretation was recently supported by Kozyreva et al. (2018) who examined a number of substorms associated with the March 2015 geomagnetic storm, observing that the maximum 1-min fluctuations occurred near the edges of the auroral electrojets, the location for which were inferred from the magnitude of the change in B_N. The bulk movement of ionospheric current systems may weaken assumptions of stationarity in the tail distributions of $|\mathrm{d}{B}_H/\mathrm{d}t|$. This may be of particular importance at CGM latitudes just below 53° which may be affected by extreme equatorward displacements of AEJ currents under intense, and as yet unobserved, geomagnetic storm conditions. Such changes in the physical environment could be addressed in future studies through physical modelling.

The MLT distributions of occurrence reported in this study confirm the earlier findings of Viljanen et al. (2001), who presented occurrence distributions of $|\mathrm{d}{B}_H/\mathrm{d}t|$ > 1 nT/s at a number of northern European stations as a function of MLT. Similar MLT and latitude distributions have been observed for Canadian magnetometers (see Fig. 3 of Boteler et al., 1998), and Beamish et al. (2002) noted a local time distribution for UK magnetometers (λ = 53°–62° N) with greater B_H variance in the hours around local midnight, reflecting the increased substorm activity in this local time sector.

Viljanen et al. (2001) observed distinct peaks in the occurrence distributions around local midnight (attributed to the expansion phase of substorms near the Harang discontinuity) and for auroral or higher latitude stations also in the morning sector (attributed to Pc5 pulsations and omega bands), with peaks occurring around 6 MLT at auroral latitudes, or later >10 MLT at the higher-latitude Svalbard location. They also noted the strong southward direction, D, of field fluctuations dB_H for large events in the auroral zone midnight sector, compared to a more East–West alignment in the morning sector. Figure 7b confirms this general pattern and provides further detail as to the MLT distribution of these changes.

In this paper, we chose to fit the extreme $|\mathrm{d}{B}_H/\mathrm{d}t|$ data by a Generalised Pareto distribution. This generalises three classes of distribution, Gumbel, Fréchet, and Weibull, which are specific cases of GP distribution with shape parameter $\xi$ = 0, $\xi$ > 0, and $\xi$ < 0, respectively (Coles, 2001). The fitting is more flexible if it is not constrained to using a particular class of extreme value distribution, since

1. Whilst this data set hints at one class of distribution being most appropriate, this is not to say that any separate but similar data set analysed in a subsequent study would give the same result. Comparison of, say, a Fréchet fit to this data with a Gumbel fit to the second data set would then be harder.

2. Selecting a particular class of distribution effectively fixes the shape parameter a priori; thus standard error and CI estimates on all return level, return period and quantile estimates would not take into account the uncertainty in this choice. Retaining the more flexible GP model, in which the shape is explicitly estimated, does.

Wintoft et al. (2016) observed a decreasing GP distribution shape parameter, $\xi$, up to the highest magnetic latitude site that they studied (Tromsø (λ = 67.0° N)) – with values tending towards zero or slightly negative ($\xi$ = −0.26) at Tromsø. However, our Figure 3c extends to even higher latitudes and indicates a clear local minimum in $\xi$ at 70° with increasing shape parameter, $\xi$ at higher latitudes towards the pole, which was not previously observed. This is largely responsible for the increase in predicted RLs toward the pole for return periods of 100 years or more (Fig. 4). There was some evidence for this trend in the “hourly range” geomagnetic activity predictions presented by Nikitina et al. (2016) from Canadian (54°–87° N) magnetometer records, though it was not clearly evident in the study of Love et al. (2016), which fitted fixed-shape (log-normal) tail distributions. If the data were log-normally distributed in the tail region, then a fit to the GP distribution would yield a shape parameter $\xi$ = 0, but the analyses of Wintoft et al. (2016), Thomson et al. (2011), Weigel & Baker (2003) and our own findings indicate that a heavier-tailed Fréchet distribution ($\xi$ > 0) provides a better model fit, except perhaps for a small number of stations at auroral latitudes.

4.2 Relation to sudden commencement timings

To illustrate the importance of SC in relation to extreme $|\mathrm{d}{B}_H/\mathrm{d}t|$ events, Figure 8 presents the likelihood that cluster peaks of $|\mathrm{d}{B}_H/\mathrm{d}t|>P_{99.97}$ occurred within (a) 5 min, (b) 5 h, and (c) 5 days of a preceding sudden commencement. The SC times were obtained from IAGA bulletins (Ebre Observatory, available online: http://www.obsebre.es/en/rapid) and filtered to retain only those events for which a majority of the five reporting magnetometers attributed a confidence factor of 2 or more (i.e., “unmistakeably a Sudden Commencement”).

Fig. 8 Fraction of cluster peaks occurring within (a) 5 min, (b) 5 h and (c) 5 days of a preceding SC.

For CGM latitudes below 55°, SC occurrences within 5 min account for up to 50% of the cluster peaks, with a greater proportion on the night side. The proportion of events preceded by SC reduces in regions of high occurrence probability (cf. Fig. 5b) and only a small proportion (<5%) of the cluster peaks at latitudes and MLT sectors associated with auroral substorms onsets, the dayside cusp, and Pc5 wave activity, are immediately preceded by SC. However, the extreme $|\mathrm{d}{B}_H/\mathrm{d}t|$ events occurring at CGM latitudes above 80° in the few hours around 6 and 18 MLT are highly likely to have been immediately preceded by a SC (with 30–50% of these events occurring within 5 min).

It is interesting to observe differences between the pre- and post-midnight distributions at 20° < $\left|\lambda \right|$ < 55° (Fig. 8a and b), where the post-midnight distribution has a greater probability of association with SC and is sharply bounded at approximately 6 MLT. For the dayside (6–18 MLT) region at 47° < $\left|\lambda \right|$ < 67°, the strong association with SC is observed within 5 days (Fig. 8c) but not within 5 h of the SC (Fig. 8b), suggesting a possible connection to the later phases of geomagnetic storms.

4.3 Dependence on IMF orientation

A clearer understanding of the magnetospheric drivers of very large $|\mathrm{d}{B}_H/\mathrm{d}t|$ is gained by examining the distribution of cluster peaks ($|\mathrm{d}{B}_H/\mathrm{d}t|>P_{99.97}$) after filtering for the IMF field orientation. Figure 9 presents the cluster peak occurrence probability distribution as shown in Figure 5b but after the dataset has been partitioned into conditions of (a) northward IMF (B_z > 2 nT), (b) low |$B_z$| < 2 nT, and (c) southward IMF (B_z < −2 nT). At low- to mid-latitudes ($\left|\lambda \right|$ < 50°) the occurrence distributions are relatively independent of B_z, as indeed the Chapman-Ferraro currents associated with sudden commencements are little influenced by the IMF B_z orientation. In the auroral region (55° < $\left|\lambda \right|$ < 75°), however, the peaks in the 20–24 MLT sector dominate during negative B_z – the condition necessary for substorms driven by magnetic field line reconnection in the magnetotail. The occurrences at 3–6 MLT near $\left|\lambda \right|$ = 64° are also enhanced under conditions of negative B_z. Positive B_z conditions lead to more frequent occurrences near the dayside cusp around noon (10–14 MLT) and in the 7–10 MLT sector at $\left|\lambda \right|$ = 60–80°, which are regions associated with strong Pc5 wave activity.

Fig. 9 Probability distribution of cluster peaks as for Figure 5b but with data filtered by IMF Bz: (a) Bz > 2 nT, (b) −2 nT ≤ Bz < 2 nT, (c) Bz < −2 nT.

In Figure 10 we examine data from magnetometers in the auroral region ($\left|\lambda \right|$ = 55–75° N) and further partition the cluster peaks by both B_y and B_z. In all MLT sectors, cluster peaks for negative B_z (panels g–i) occur most frequently at the vernal equinox for negative B_y (panel g) and at the autumnal equinox for positive B_y (panel i). The converse is true under positive B_z conditions (panels a–c) with cluster peaks occurring most frequently at the autumnal (vernal) equinox under negative (positive) B_y conditions. These patterns indicate a strong relation to the Russell-McPherron effect, which ascribes seasonal and diurnal changes in magnetic field line reconnection and geomagnetic activity to changes in the orientation of the Earth’s dipole axis with respect to the IMF (Russell & McPherron, 1973).

Fig. 10 Occurrence probability of $|\mathrm{d}{B}_H/\mathrm{d}t|$ > P99.97 (after declustering) as a function of month and MLT, for NH stations between 55° N and 75° N CGM, filtered by IMF orientation: Panels (a)–(c) Bz ≥ 2 nT, (d)–(f) −2 nT ≤ Bz < 2 nT, (g)–(i) Bz ≥ 2 nT, Panels (a), (d), (g) By < −2 nT; panels (b), (e), (h) −2 nT ≤ By < 2 nT; panels (c), (f), (i) By ≥ 2 nT.

Previous studies have shown that the MLT of the cusp region is controlled by the IMF B_y component. In the northern hemisphere, the cusp is shifted pre-noon (post-noon) for negative (positive) B_y, with the opposite behaviour in the southern hemisphere. This shift is more pronounced under negative B_z conditions when the cusp lies further equatorward (Newell et al., 1989). To examine the effect this might have on the distribution of cluster peaks ($|\mathrm{d}{B}_H/\mathrm{d}t|>P_{99.97}$) an analysis of B_y and B_z dependences is presented in Figure 11 for all magnetometers in the NH cusp and polar region, $\left|\lambda \right|$ > 77° N. This shows that under all conditions of B_z, the locus of cluster peaks ($|\mathrm{d}{B}_H/\mathrm{d}t|>P_{99.97}$) shifts from about 10:30 MLT for negative B_y conditions, to about 12:30 MLT for positive B_y.

Fig. 11 As Figure 10, but for NH stations above 77° N CGM.

The polarity of the R0 field-aligned currents entering the cusp region have been previously shown to depend on the IMF B_y component, (being upward (downward) in the northern hemisphere (NH) for B_y > 0 (B_y < 0), and opposite in the southern hemisphere (SH)) (Milan et al., 2017, and references therein). However, we found that the direction, D, of the peak fluctuations for $\left|\lambda \right|$ > 77° N near local noon (Figure S2 in Supplementary Material (folder S4)) was relatively uniform (isotropic) and did not change significantly with B_y or B_z. This indicates that the magnetic fluctuations associated with field-aligned currents are too weak to significantly affect the highest percentiles of $|\mathrm{d}{B}_H/\mathrm{d}t|$, due in part to the shielding effect of the conductive ionospheric layer.

5 Summary and conclusions

In this study we have used Extreme Value Theory (EVT) to fit Generalised Pareto (GP) tail distributions to a large global data set of horizontal magnetic field fluctuations, dB_H/dt from 125 magnetometers. The variation of the fitted GP parameters and derived 5- to 500-year return levels (RL) have been modelled as functions of the corrected geomagnetic (CGM) latitude, λ. A combination of high scale (σ) and shape ($\xi$) parameters in the region of $\left|\lambda \right|$ = 53° (in both hemispheres) creates a sharp maximum in RL at this latitude. Values of $\xi$ increasing above $\left|\lambda \right|$ = 75° lead to a trend towards increasing RLs for locations in the polar cap for return periods greater than 100 years.

The rate of occurrence of $|\mathrm{d}{B}_H/\mathrm{d}t|$ measurements exceeding the 99.97th percentile (P_99.97) has also been modelled as a function of CGM latitude, month, MLT, and the direction of dB_H. The occurrence of large $|\mathrm{d}{B}_H/\mathrm{d}t|$ has been shown to be strongly dependent on CGM latitude and MLT sector, such that a GIC observed in one ground-based system can be expected to differ greatly from that in others well-separated in latitude or longitude (or local time). By fitting analytical surfaces (such as a polynomial or spherical harmonic expansion) to the 2-d occurrence probability profile, e.g. $P\left(\right|\mathrm{d}{B}_H/\mathrm{d}t|>P_{99.97};\lambda , \mathrm{MLT})$ a fully analytical model of return levels may be generated, parameterised by these covariates.

The occurrence probability profile has distinct maxima in the pre-midnight auroral zone ($\left|\lambda \right|$ = 55–75°), with predominantly southward dB_H, that may be attributed to substorm-related westward ionospheric currents. The probability maximises near the equinoxes (with a secondary winter maximum) principally under conditions of negative IMF B_z.

The auroral zone 03–11 MLT sector is characterised by an increase in event occurrences that matches known patterns of ULF Pc5 wave activity. Events occur preferentially at earlier local times and at lower latitudes under negative B_z, with a more isotropic directional distribution centred broadly about the 100°–280° bearing axis. In the northern hemisphere (at least), the earlier MLT occurrences are principally observed near the June solstice, with later MLT occurrences towards December. Both substorm-related and Pc5-related peaks (at all local times) exhibit a seasonal dependence on IMF orientation consistent with the Russell-McPherron effect.

In the region poleward of the polar cusp $\left|\lambda \right|$ > 75° in the summer months there is a raised probability of occurrence at 10–14 MLT, peaking slightly earlier (later) for negative (positive) IMF B_y, with greater occurrence probability under positive B_z conditions. The directional distribution of these changes is nearly isotropic.

At low latitudes ($\left|\lambda \right|$ < 40°), occurrence probability $P\left(\right|\mathrm{d}{B}_H/\mathrm{d}t|>P_{99.97};\lambda , \mathrm{MLT})$ forms a V-shaped distribution in latitude versus MLT, closely distributed about 11 MLT at the equator but extended towards dawn and dusk towards mid-latitudes. Large magnetic field fluctuations in this region are principally northward directed, and may be attributed to sudden commencements (SC). There is a raised probability of events being observed within 5-min of a SC in latitude zones and MLT sectors for which there are relatively few occurrences (night-side, low latitude, or near dawn or dusk at polar-latitudes ($\left|\lambda \right|$ > 80°).

Results from this study may be applied in the evaluation and mitigation of GIC risks. Models of maximum likelihood estimates for $|\mathrm{d}{B}_H/\mathrm{d}t|$ versus geomagnetic latitude (Sect. 3.1) are provided for return periods of up to 500 years. Combining these predictions with an Earth conductivity model would yield predictions of the extreme geoelectric fields induced at the Earth’s surface. When these are input to engineering models, such as HV electricity network models, the return levels for GICs may be estimated. Models of the occurrence probability of $P\left(\right|\mathrm{d}{B}_H/\mathrm{d}t|>P_{99.97})$ with month and MLT may be used to refine RL estimates of $|\mathrm{d}{B}_H/\mathrm{d}t|$for operations limited to certain seasons or times of day. Furthermore, information from the pattern of occurrence probability versus direction (Fig. 7) could be used to refine GIC risk estimates for individual linear components of ground infrastructure, such as sections of pipeline, cables, or elements of an HV electricity network.

It should be noted that the statistics presented in this paper reflect only the minute-to-minute changes in magnetic field. Further studies (not reported here) have examined the directional statistics of magnetic field fluctuations on longer time scales up to 60 min, which are also important for GIC modelling but are governed by different temporal characteristics of the ionospheric and magnetospheric current systems. We shall report these findings in a future paper. Sub-minute geomagnetic field fluctuations are also effective in inducing surface electric fields from shallow depths in the Earth and these will be examined in future studies using a much more limited set of magnetograms recorded at 1- or 10-s cadence.

Supplementary materials

Supplementary Material (folder S1): Geodetic and mean Corrected Geomagnetic (CG) coordinates for the 125 selected magnetometer sites, in comma-separated-variable (CSV) format.

Supplementary Material (folder S2): Spline models fitted to the GPD parameters, 99.97th percentile and return levels for |dB_H/dt|.

Supplementary Material (folder S3): CSV files containing A_lm spherical harmonic function coefficients (6), and associated plotting software.

Supplementary Material (folder S4): Additional figures (Fig. S1 and S2).

Acknowledgments

This work was funded by the UK Natural Environment Research Council (NERC) under grant number NE/P016715/1. The authors are grateful to members of the Space Weather Impacts on Ground-based Systems (SWIGS) consortium led by Alan Thomson (British Geological Survey), with special thanks to Mervyn Freeman (British Antarctic Survey), Tim Yeoman (University of Leicester, UK), Malcolm Dunlop (RAL Space, UK), and Phil Livermore (University of Leeds, UK) for helpful discussions. Sunspot numbers were provided by WDC-SILSO, Royal Observatory of Belgium, Brussels and are available online: http://sidc.oma.be/silso/DATA/SN_ms_tot_V2.0.txt. Interplanetary magnetic field data were provided by the NASA Goddard Space Flight Center Space Physics Data Facility (available from https://omniweb.gsfc.nasa.gov). The ground magnetometer data and the substorm list were provided by SuperMAG (available from http://supermag.jhuapl.edu) and we gratefully acknowledge contributions from the SuperMAG collaborators: Intermagnet; USGS, Jeffrey J. Love; CARISMA, PI Ian Mann; CANMOS; The S–RAMP Database, PI K. Yumoto and Dr. K. Shiokawa; The SPIDR database; AARI, PI Oleg Troshichev; The MACCS program, PI M. Engebretson, Geomagnetism Unit of the Geological Survey of Canada; GIMA; MEASURE, UCLA IGPP and Florida Institute of Technology; SAMBA, PI Eftyhia Zesta; 210 Chain, PI K. Yumoto; SAMNET, PI Farideh Honary; The IMAGE magnetometer network, PI L. Juusola; AUTUMN, PI Martin Connors; DTU Space, PI Anna Willer; South Pole and McMurdo Magnetometer, PI’s Louis J. Lanzarotti and Alan T. Weatherwax; ICESTAR; RAPIDMAG; British Antarctic Survey; McMac, PI Dr Peter Chi; BGS, PI Dr Susan Macmillan; Pushkov Institute of Terrestrial Magnetism, Ionosphere and Radio Wave Propagation (IZMIRAN); GFZ, PI Dr Juergen Matzka; MFGI, PI B. Heilig; IGFPAS, PI J. Reda; University of L’Aquila, PI M. Vellante; BCMT, V. Lesur and A. Chambodut; Data obtained in cooperation with Geoscience Australia, PI Marina Costelloe; AALPIP, co-PIs Bob Clauer and Michael Hartinger; SuperMAG, PI Jesper W. Gjerloev; Sodankylä Geophysical Observatory, PI Tero Raita; Polar Geophysical Institute, Alexander Yahnin and Yarolav Sakharov; Geological Survey of Sweden, Gerhard Schwartz; Swedish Institute of Space Physics, Mastoshi Yamauchi; UiT the Arctic University of Norway, Magnar G. Johnsen; Finnish Meteorological Institute, PI Kirsti Kauristie. The editor thanks Sean Elvidge and an anonymous reviewer for their assistance in evaluating this paper.

References

All Figures

	Fig. 1 Three examples of large field fluctuations $|\mathrm{d}{B}_H/\mathrm{d}t|$. (a) An auroral substorm at BRW, Alaska, (b) an SSC at GUA, Guam, 6° N, 144° W CGM), and (c) Pc5 pulsations (0945–1015 UT) at NAQ, Greenland. The 99.97th percentile of $|\mathrm{d}{B}_H/\mathrm{d}t|$ is also shown.
In the text

	Fig. 2 Magnetometer sites in the SuperMAG collaboration with ≥20 years data, coloured by P99.97 at each site. Contours show CGM latitudes for IGRF model epoch 2000.
In the text

Fig. 3 GP parameters for 125 magnetometers versus absolute CGM latitude. (a) 99.97-percentiles of $|\mathrm{d}{B}_H/\mathrm{d}t|$ and revised thresholds, u, (b) MLEs of scale, σ (error bars indicate 95% CIs), (c) MLEs of shape, $\xi$, (d) extremal index, $\widehat{\theta }$, (e) number of measurements, n, number exceeding P99.97, nu, number of clusters, nc, and mean nc per year. Red curves indicate fitted smoothed splines (5).

In the text

	Fig. 4 Return levels of $|\mathrm{d}{B}_H/\mathrm{d}t|$ (MLEs, with 95% CIs shown as error bars), for return periods of (a) 100 years, (b) 200 years, and (c) 500 years with fitted smoothed splines ($s\left(\left|\mathrm{\lambda }\right|\right)$ in (5)). Panel (d) shows spline fits for a range of return periods.
In the text

	Fig. 5 Total occurrences of cluster peaks binned by geomagnetic latitude $\left|\lambda \right|$ and MLT. (a) Number of clusters, (b) probability distribution of cluster peaks, (c) fitted spherical harmonic expansion (Psph in (6)).
In the text

	Fig. 6 Occurrence probability of $|\mathrm{d}{B}_H/\mathrm{d}t|$ > P99.97 (after declustering) as a function of month and MLT. Panels (a)–(c) are for NH stations, whilst (d)–(f) represent SH stations. (a) and (d) represent low-latitude stations (<40° CGM), (b) and (e) represent auroral latitudes 55°–75° CGM) and (c) and (f) represent latitudes above 77° CGM.
In the text

	Fig. 7 Occurrence probabilities in the same format as Figure 6 but shown as a function of dBH direction, D, and MLT. Bin sizes are 1 h × 10°.
In the text

	Fig. 8 Fraction of cluster peaks occurring within (a) 5 min, (b) 5 h and (c) 5 days of a preceding SC.
In the text

	Fig. 9 Probability distribution of cluster peaks as for Figure 5b but with data filtered by IMF Bz: (a) Bz > 2 nT, (b) −2 nT ≤ Bz < 2 nT, (c) Bz < −2 nT.
In the text

	Fig. 10 Occurrence probability of $|\mathrm{d}{B}_H/\mathrm{d}t|$ > P99.97 (after declustering) as a function of month and MLT, for NH stations between 55° N and 75° N CGM, filtered by IMF orientation: Panels (a)–(c) Bz ≥ 2 nT, (d)–(f) −2 nT ≤ Bz < 2 nT, (g)–(i) Bz ≥ 2 nT, Panels (a), (d), (g) By < −2 nT; panels (b), (e), (h) −2 nT ≤ By < 2 nT; panels (c), (f), (i) By ≥ 2 nT.
In the text

	Fig. 11 As Figure 10, but for NH stations above 77° N CGM.
In the text